#include <vector>
using namespace std;

#include <d3clipping.h>

typedef point3<double> pt3;
typedef point3<double> const pt3c;

void d3clipping::reprocessStack
(
  uintc nstacksz,
  uintc s,
  uintc neib0,
  uintc neib1
)
{
  processSimplex(s);

  if (neib0!=0)
    processSimplex(neib0);
  if (neib1!=0)
    processSimplex(neib1);

  uintc nstack1 = vi.size();
  for (uint i=nstacksz; i<nstack1; ++i)
    processSimplex(i);

}


void d3clipping::addintersectionpoint( pt3c & A, pt3c & B )
{
  assert(bv.size()==pt.size());
  // Any point added is assumed to be black.

  static pt3 p;
  H.intersection(p,A,B);

  pt.push_back(p);
  bv.push_back(true);
}


d3clipping::d3clipping
(
  d3tess & tess_, 
  partitionspace< pt3, double > const & H_
)
  : tess(tess_), vi(tess.vi), pt(tess.pt), H(H_)
{
}


void d3clipping::processB1W2
( 
  uintc s, 
  uintc b0, 
  uintc w0, 
  uintc w1
)
{
  assert(s>0);
  uintc nstack=vi.size();
  assert(s<nstack);
  assert(b0<pt.size());
  assert(w0<pt.size());
  assert(w1<pt.size());
  assert(b0>0);
  assert(w0>0);
  assert(w1>0);

  pt3c & B0(pt[b0]);
  pt3c & W0(pt[w0]);
  pt3c & W1(pt[w1]);

  if ( H.isOnBoundary(B0,zero) )
  {
    tess.simplexdelete(s);
    return;
  }

  // Neighboring simplexes are reprocessed.

  uintc neib0 = vi[s].nifrom(w0);
  uintc neib1 = vi[s].nifrom(w1);

  // The simplex is cut in two.

  uintc z0 = pt.size();
  uintc z1 = z0 + 1;
  addintersectionpoint(B0,W0);
  addintersectionpoint(B0,W1);

  tess.splitwithline(s,w1,z0,w0,z1);

  reprocessStack(nstack,s,neib0,neib1);
}


void d3clipping::processW1B2
( 
  uintc s, 
  uintc w0, 
  uintc b0, 
  uintc b1 
)
{
  assert(s>0);
  uintc nstack = vi.size();
  assert(s<nstack);
  assert(w0<pt.size());
  assert(b0<pt.size());
  assert(b1<pt.size());
  assert(w0>0);
  assert(b0>0);
  assert(b1>0);

  pt3c & W0(pt[w0]);
  pt3c & B0(pt[b0]);
  pt3c & B1(pt[b1]);

  if ( H.isOnBoundary(B0,zero) )
  {
    if ( H.isOnBoundary(B1,zero) )
    {
      tess.simplexdelete(s);
      return;
    }

    uintc neib0 = vi[s].nifrom(b0);

    addintersectionpoint(B1,W0);

    tess.split1D(s,b0,pt.size()-1);

    reprocessStack(nstack,s,neib0);

    return;
  }

  if ( H.isOnBoundary(B1,zero) )
  {
    uintc neib1 = vi[s].nifrom(b1);

    addintersectionpoint(B0,W0);

    tess.split1D(s,b1,pt.size()-1);

    reprocessStack(nstack,s,neib1);

    return;
  }

  // Line split

  // Neighboring simplexes are reprocessed.

  uintc neib0 = vi[s].nifrom(b0);
  uintc neib1 = vi[s].nifrom(b1);

  // The simplex is cut in two. 

  uintc z0 = pt.size();
  addintersectionpoint(B0,W0);
  addintersectionpoint(B1,W0);
  uintc z1 = z0 + 1;

  tess.splitwithline(s,b0,z1,b1,z0);

  reprocessStack(nstack,s,neib0,neib1);
}

void d3clipping::processSimplex(uintc s)
{
  assert(s>0);
  assert(s<vi.size());

  // Initilly classify the point without considering
  //   if the point is on the boundary.

  d3tri const * t = & vi[s];
  if (t->isnull())
    return;

  uintc p0 = t->pi[0];
  uintc p1 = t->pi[1];
  uintc p2 = t->pi[2];
    
  uint k=0;
  if (bv[p0])
    k += 1;
  if (bv[p1])
    k += 2;
  if (bv[p2])
    k += 4;

  switch (k)
  {
    // All black points, accept the triangle.
    case 7:  
      break; 

    // All white points, delete the triangle.
    case 0:  
      tess.simplexdelete(s);
      break;

    // One black two white
    case 1: 
      processB1W2(s,p0,p1,p2);
      break;

    case 2: 
      processB1W2(s,p1,p0,p2);
      break;

    case 4: 
      processB1W2(s,p2,p0,p1);
      break;

    // Two black on white
    case 3:
      processW1B2(s,p2,p0,p1);
      break;

    case 5:
      processW1B2(s,p1,p0,p2);
      break;

    case 6:
      processW1B2(s,p0,p1,p2);
      break;
  }
}


//  Note:  I could merge processSimplex(uintc s) into
//    eval() however I will wait for a while to see how
//    the algorithm evolves.  The merge I believe would
//    produce faster code as memory does not have to be
//    thrashed around.
void d3clipping::eval()
{
  uint ptsz = pt.size();
  bv.resize(ptsz);
  // Ignore the first element as point 0 is defined as no point.
  bv[0]=true;
  // Because the partition H is subtracting from the mesh its
  // balls are white.  So the H.isInside function has to be
  // inverted.  If the point is on H's boundary there is a possibiliby
  // that it could be classified wrongly.  In this case I believe
  // the algorithm is robust enough to handle it by re interpreting any
  // 2 black 1 white cases to 3 black again because such cases are 
  // re processed anyway.
  for (uint i=1; i<ptsz; ++i)
    bv[i] = ! H.isInside( pt[i] );

  for (uint i=1; i<vi.size(); ++i )
    processSimplex(i);
}