///////////////////////////////////////////////////////////////////////////////
//
// tin.c
//
///////////////////////////////////////////////////////////////////////////////



#include "tin.h"
#define DEBUG if(0)

// dummy pointer for all triangles with maxE < e
COORD DONE[3] = {0,0,0};

// dummy pointer for pointlist when triangle should not be drawn
QUEUE NO_DRAW;

// dummy value for corners with nodata values
COORD NODATA_CORNER = -32767;





//
// add trinagle to a TIN and return pointer to it
//
TRIANGLE* addTri(COORD* p1, COORD* p2, COORD* p3, 
		 TRIANGLE *t12, TRIANGLE *t13, TRIANGLE *t23) {
  assert(p1 && p2 && p3);
  
  // Validate that the triangle is not collinear
  assert(areaSign(p1,p2,p3));


  // Allocate space for this triangle and give it values
  TRIANGLE *tn;
  tn = (TRIANGLE*)malloc(sizeof(TRIANGLE));
  assert(tn);
  
  // Assign values to the triangle
  tn->maxE = NULL;
  tn->points = NULL;
  tn->p1=p1;
  tn->p2=p2;
  tn->p3=p3;
  tn->p1p2 = t12;
  tn->p1p3 = t13;
  tn->p2p3 = t23;
 
  // point neighbor triangles to this triangle
  pointNeighborTo(tn,p1,p2,t12);
  pointNeighborTo(tn,p1,p3,t13);
  pointNeighborTo(tn,p2,p3,t23);

  return tn;
}


//
// Given a triangle and 2 pts return the neighbor triangle to that edge
//
TRIANGLE* whichTri(TRIANGLE *tn,COORD *pa,COORD *pb){
  assert(tn && pa && pb);
  if( (tn->p1 == pa || tn->p2 == pa) && (tn->p1 == pb || tn->p2 == pb) ){
    assert(tn->p1p2 != tn);
    return tn->p1p2;
  }
  if( (tn->p1 == pa || tn->p3 == pa) && (tn->p1 == pb || tn->p3 == pb) ){
    assert(tn->p1p3 != tn);
    return tn->p1p3;
  }
  else if ( (tn->p2 == pa || tn->p3 == pa) && (tn->p2 == pb || tn->p3 == pb) ){
    assert(tn->p2p3 != tn);
    return tn->p2p3;
  }
  else {
    assert(0);
    exit(1);
    return NULL;
  }
}


//
// Find neighbor (t) to tn and point proper edge to t.
//
void pointNeighborTo(TRIANGLE *t,COORD *pa,COORD *pb,TRIANGLE *tn){
  assert(t);
  if(tn != NULL){
    if( (tn->p1 == pa || tn->p2 == pa) && (tn->p1 == pb || tn->p2 == pb) )
      tn->p1p2 = t;
    else if( (tn->p1 == pa || tn->p3 == pa) && (tn->p1 == pb || tn->p3 == pb) )
      tn->p1p3 = t;
    else if ( (tn->p2 == pa || tn->p3 == pa) && (tn->p2 == pb || tn->p3 == pb))
      tn->p2p3 = t;
    else{
      assert(0);
      exit(1);
    }
  }
}


//
// Print the triangles in a TIN for debugging
//
void printTin(TIN* tin){
  TRIANGLE *curT = tin->t;
  TRIANGLE *prevT = curT;
  
  // Create an edge which 
  EDGE curE;
  curE.type = IN;
  curE.t1 = NULL;
  curE.t2 = tin->t;
  curE.p1 = tin->e.p1;
  curE.p2 = tin->e.p2;
  printf("\n\n PRINTING TIN \n\n");

  do{
    // Print tri if we are on its IN edge
    if(curE.type == IN && prevT != NULL){
      printTriangle(prevT);
    }

    // Go to next edge
    prevT=curT;
    curT = nextEdge(curT,tin->v,&curE);
    assert(curT);
    
  }while((curE.p1 != tin->e.p1 || curE.p2 != tin->e.p2) && 
	 (curE.p1 != tin->e.p2 || curE.p2 != tin->e.p1));
  
}


//
// Print a tin list for debugging
//
void printTinList(TINLIST tl){
  // Skip dummy head
  tl = tl->next;
  while(tl->next != NULL){
    printf("\n\n Print TNODE \n\n");
    printTin(tl->t);
    tl = tl->next;
  }
}


//
// Print the points of a triangle and pointers for debugging
//
void printTriangle(TRIANGLE* t){
  printf("Tri: t:%p Points: %p(%d,%d) %p(%d,%d) %p(%d,%d) Tris: %p %p %p\n",
	 t,t->p1,t->p1[X],t->p1[Y],t->p2,t->p2[X],t->p2[Y],t->p3,t->p3[X],
	 t->p3[Y],t->p1p2,t->p1p3, t->p2p3);
  //printPointList(t);
}

void printTriangleCoords(TRIANGLE* t){
  printf("Triangle t: Points: (%d,%d,%d) (%d,%d,%d) (%d,%d,%d)\n",
	 t->p1[X],t->p1[Y],t->p1[Z], t->p2[X],t->p2[Y],t->p2[Z],
	 t->p3[X], t->p3[Y],t->p3[Z]);
  //printPointList(t);
}


//
// This can be used to traverse the tin. Since this will cross every
// triangle 3 times, it is is important to only use next tri when the
// returned edge is of type IN
//
TRIANGLE *nextEdge(TRIANGLE *t, COORD* v, EDGE *edge){
  assert(t && v && edge && edge->p1 && edge->p2);
  EDGE e12, e13, e23;
  short inCount,outCount;
  int areaOp, areaV;
  EDGE *in, *out, *inback, *outback;
  
  in = out = inback = outback = NULL;
  
  inCount = outCount = 1;
  
  //
  // Define the three edges of this triangle to work with them. If t1
  // is NULL then t1 should point to this triangle since the traversal
  // calls for boundary edges to point to the triangle they came from
  //
  
  if(t->p1p2 == NULL)
    e12.t1 = t;
  else
    e12.t1 = t->p1p2;
  e12.t2 = t;
  e12.p1 = t->p1;
  e12.p2 = t->p2;
  if(t->p1p3 == NULL)
    e13.t1 = t;
  else
    e13.t1 = t->p1p3;
  e13.t2 = t;
  e13.p1 = t->p1;
  e13.p2 = t->p3;
  if(t->p2p3 == NULL)
    e23.t1 = t;
  else
    e23.t1 = t->p2p3;
  e23.t2 = t;
  e23.p1 = t->p2;
  e23.p2 = t->p3;
  
  //
  // This code figures out which edges are IN and which are OUT. It
  // assumes that if areaV is 0 (collinear) then the edge is IN, we
  // will handle this later in the next section of code
  //
  
  // Edge p1,p2
  areaOp = areaSign(t->p1,t->p2,t->p3);
  areaV = areaSign(t->p1,t->p2,v);
  if(areaOp > 0 && areaV <= 0)
    e12.type = IN;
  else if(areaOp < 0 && areaV >= 0 )
    e12.type = IN;
  else
    e12.type = OUT;
  // Edge p1,p3
  areaOp = areaSign(t->p1,t->p3,t->p2);
  areaV = areaSign(t->p1,t->p3,v);
  if(areaOp > 0 && areaV <= 0)
    e13.type = IN;
  else if(areaOp < 0 && areaV >= 0 )
    e13.type = IN;
  else
    e13.type = OUT;
  // Edge p2,p3
  areaOp = areaSign(t->p2,t->p3,t->p1);
  areaV = areaSign(t->p2,t->p3,v);
  if(areaOp > 0 && areaV <= 0)
    e23.type = IN;
  else if(areaOp < 0 && areaV >= 0 )
    e23.type = IN;
  else
    e23.type = OUT;

  //
  // At this point edges are just marked in or out, we need to handle
  // 2 outs or 2 ins. If there are 2 in's mark the left one as the
  // real in, change the count and backward link the variables. There
  // is also the chance that one of two degenerate cases. The first is
  // that v may be one of the point of the trianle which means that we
  // are opperating on the lower left most triangle. This is unique
  // because there will be two collinear edge with v. The other case
  // is where only one of the edges is collinear with v. In the either
  // case there is a perscribed solution to the problem in the code
  // below.
  //

  if(e12.type == IN && e13.type == IN){
    areaV = areaSign(v,t->p1,t->p2);
    areaOp = areaSign(v,t->p1,t->p3);
    // e12 is the real in edge and there are no collinears
    if( areaV > 0  && areaOp < 0){
      inback = &e13;
      in = &e12;
      out = &e23;
      e13.type = INBACK;
      inCount++;
    }
    // e13 is the real in edge and there are no collinears
    else if ( areaV < 0 && areaOp > 0){
      inback = &e12;
      in = &e13;
      out = &e23;
      e12.type = INBACK;
      inCount++;
    }
    // e13 is collinear and p2 is right of e13
    else if ( areaV < 0 && areaOp == 0){
      in = &e12;
      out = &e13;
      outback = &e23;
      e13.type = OUT;
      e23.type = OUTBACK;
      outCount++;
    }
    // e13 is collinear and p2 is left of e13
    else if ( areaV > 0 && areaOp == 0){
      inback = &e13;
      in = &e12;
      out = &e23;
      e13.type = INBACK;
      inCount++;
    }
    // e12 is collinear and p3 is right of e12
    else if ( areaV == 0 && areaOp < 0){
      in = &e13;
      out = &e12;
      outback = &e23;
      e12.type = OUT;
      e23.type = OUTBACK;
      outCount++;
    }
    // e12 is collinear and p3 is left of e12
    else if ( areaV == 0 && areaOp > 0){
      inback = &e12;
      in = &e13;
      out = &e23;
      e12.type = INBACK;
      inCount++;
    }
    // They are both collinear, this must be the lower left corner
    else {
      assert(areaV == 0 && areaOp == 0);
      // Determine which edge is left most
      areaV = areaSign(t->p1,t->p2,t->p3);
      // areaV should not be 0
      assert(areaV); 
      // if p3 is left of p1p2 then e13 is OUT
      if(areaV > 0){
	in = &e12;
	out = &e13;
	outback = &e23;
	e13.type = OUT;
	e23.type = OUTBACK;
	outCount++;
      }
      // if p3 is right of p1p2 then e12 is OUT
      else {
	in = &e13;
	out = &e12;
	outback = &e23;
	e12.type = OUT;
	e23.type = OUTBACK;
	outCount++;
      }
    }
  }

  if(e12.type == IN && e23.type == IN){
    areaV = areaSign(v,t->p2,t->p1);
    areaOp = areaSign(v,t->p2,t->p3);
    // e12 is the real in edge and there are no collinears
    if( areaV > 0  && areaOp < 0){
      inback = &e23;
      in = &e12;
      out = &e13;
      e23.type = INBACK;
      inCount++;
    }
    // e23 is the real in edge and there are no collinears
    else if ( areaV < 0 && areaOp > 0){
      inback = &e12;
      in = &e23;
      out = &e13;
      e12.type = INBACK;
      inCount++;
    }
    // e23 is collinear and p1 is right of e23
    else if ( areaV < 0 && areaOp == 0){
      in = &e12;
      out = &e23;
      outback = &e13;
      e23.type = OUT;
      e13.type = OUTBACK;
      outCount++;
    }
    // e23 is collinear and p1 is left of e23
    else if ( areaV > 0 && areaOp == 0){
      inback = &e23;
      in = &e12;
      out = &e13;
      e23.type = INBACK;
      inCount++;
    }
    // e12 is collinear and p3 is right of e12
    else if ( areaV == 0 && areaOp < 0){
      in = &e23;
      out = &e12;
      outback = &e13;
      e12.type = OUT;
      e13.type = OUTBACK;
      outCount++;
    }
    // e12 is collinear and p3 is left of e12
    else if ( areaV == 0 && areaOp > 0){
      inback = &e12;
      in = &e23;
      out = &e13;
      e12.type = INBACK;
      inCount++;
    }
    // They are both collinear, this must be the lower left corner
    else {
      assert(areaV == 0 && areaOp == 0);
      // Determine which edge is left most
      areaV = areaSign(t->p2,t->p1,t->p3);
      // areaV should not be 0
      assert(areaV); 
      // if p3 is left of p2p1 then e23 is OUT
      if(areaV > 0){
	in = &e12;
	out = &e23;
	outback = &e13;
	e23.type = OUT;
	e13.type = OUTBACK;
	outCount++;
      }
      // if p3 is right of p2p1 then e12 is OUT
      else {
	in = &e23;
	out = &e12;
	outback = &e13;
	e12.type = OUT;
	e13.type = OUTBACK;
	outCount++;
      }
    }
  }

  if(e13.type == IN && e23.type == IN){
    areaV = areaSign(v,t->p3,t->p1);
    areaOp = areaSign(v,t->p3,t->p2);
    // e13 is the real in edge and there are no collinears
    if( areaV > 0  && areaOp < 0){
      inback = &e23;
      in = &e13;
      out = &e12;
      e23.type = INBACK;
      inCount++;
    }
    // e23 is the real in edge and there are no collinears
    else if ( areaV < 0 && areaOp > 0){
      inback = &e13;
      in = &e23;
      out = &e12;
      e13.type = INBACK;
      inCount++;
    }
    // e23 is collinear and p1 is right of e23
    else if ( areaV < 0 && areaOp == 0){
      in = &e13;
      out = &e23;
      outback = &e12;
      e23.type = OUT;
      e12.type = OUTBACK;
      outCount++;
    }
    // e23 is collinear and p1 is left of e23
    else if ( areaV > 0 && areaOp == 0){
      inback = &e23;
      in = &e13;
      out = &e12;
      e23.type = INBACK;
      inCount++;
    }
    // e13 is collinear and p2 is right of e13
    else if ( areaV == 0 && areaOp < 0){
      in = &e23;
      out = &e13;
      outback = &e12;
      e13.type = OUT;
      e12.type = OUTBACK;
      outCount++;
    }
    // e13 is collinear and p2 is left of e13
    else if ( areaV == 0 && areaOp > 0){
      inback = &e13;
      in = &e23;
      out = &e12;
      e13.type = INBACK;
      inCount++;
    }
    // They are both collinear, this must be the lower left corner
    else {
      assert(areaV == 0 && areaOp == 0);
      // Determine which edge is left most
      areaV = areaSign(t->p3,t->p1,t->p2);
      // areaV should not be 0
      assert(areaV); 
      // if p2 is left of p3p1 then e23 is OUT
      if(areaV > 0){
	in = &e13;
	out = &e23;
	outback = &e12;
	e23.type = OUT;
	e12.type = OUTBACK;
	outCount++;
      }
      // if p2 is right of p3p1 then e13 is OUT
      else {
	in = &e23;
	out = &e13;
	outback = &e12;
	e13.type = OUT;
	e12.type = OUTBACK;
	outCount++;
      }
    }
  }


  // If there are 2 out's mark the left one as the real out
  if(e12.type == OUT && e13.type == OUT){
    areaV = areaSign(v,t->p1,t->p2);
    outCount++;
    if( areaV >= 0 ){
      outback = &e13;
      out = &e12;
      in = &e23;
      e13.type = OUTBACK;
    }
    else{
      outback = &e12;
      out = &e13;
      in = &e23;
      e12.type = OUTBACK;
    }
  }
  if(e12.type == OUT && e23.type == OUT){
    areaV = areaSign(v,t->p2,t->p1);
    outCount++;
    if( areaV >= 0 ){
      outback = &e23;
      out = &e12;
      in = &e13;
      e23.type = OUTBACK;
    }
    else{
      outback = &e12;
      out = &e23;
      in = &e13;
      e12.type = OUTBACK;
    }
  }
  if(e13.type == OUT && e23.type == OUT){
    areaV = areaSign(v,t->p3,t->p1);
    outCount++;
    if( areaV >= 0 ){
      outback = &e23;
      out = &e13;
      in = &e12;
      e23.type = OUTBACK;
    }
    else{
      outback = &e13;
      out = &e23;
      in = &e12;
      e13.type = OUTBACK;
    }
  }

  
  // Make sure in/out count is not messed up
  assert((inCount == 1 && outCount == 2) || (inCount == 2 && outCount == 1));

  // Determine the new type of edge
  if((t->p1 == edge->p1 || t->p1 == edge->p2) && 
     (t->p2 == edge->p1 || t->p2 == edge->p2))
    edge->type = e12.type; 
  if((t->p1 == edge->p1 || t->p1 == edge->p2) && 
     (t->p3 == edge->p1 || t->p3 == edge->p2))
    edge->type = e13.type; 
  if((t->p2 == edge->p1 || t->p2 == edge->p2) && 
     (t->p3 == edge->p1 || t->p3 == edge->p2))
    edge->type = e23.type; 

  DEBUG{ 
    printf("e12: %d e13: %d e23: %d edge: %d t: %p\n",e12.type,e13.type,
	   e23.type,edge->type,t);
    fflush(stdout);}
  // Decide which triangle is next
  
  // Only one in edge
  if(inCount == 1 && edge->type == IN){
    // Change e to be the out edge of t but keep t the same so we
    // return the next tri
    edge->p1 = out->p1;
    edge->p2 = out->p2;
    DEBUG{printf("Case1 - 1 IN\n");fflush(stdout);}
    return out->t1;
  }
  else if(edge->type == IN){
    edge->p1 = out->p1;
    edge->p2 = out->p2;
    DEBUG{printf("Case2 - 2 IN\n");fflush(stdout);}
    return out->t1;
  }
  else if(edge->type == INBACK){
    DEBUG{printf("Case3 - INBACK\n");fflush(stdout);}
    return inback->t1;
  }
  else if(outCount == 2 && edge->type == OUT){
    edge->p1 = outback->p1;
    edge->p2 = outback->p2;
    DEBUG{printf("Case4 - 2 OUT\n");fflush(stdout);}
    return outback->t1;
  }
  else if(outCount == 2 && edge->type == OUTBACK){
    edge->p1 = in->p1;
    edge->p2 = in->p2;
    DEBUG{printf("Case5 - 2 OUT: Outback\n");fflush(stdout);}
    return in->t1;
  }
  // if there is only one out edge
  else{
    edge->p1 = in->p1;
    edge->p2 = in->p2;
    DEBUG{printf("Case6 - 1 Out\n");fflush(stdout);}
    return in->t1;
  }
}


//
// Remove triangle from TIN and free memory
//
void removeTri(TRIANGLE* t) {
  assert(t);
  free(t);
}


//
// printPointList - given a triangle,print its point list
//
void printPointList(TRIANGLE* t){
  // loop through each node in the list and free the points
  if(t != NULL){
    if(t->points == NULL)
      printf("\t Triangle: %p has no point list to print\n",t);
    else {
      printf("Point list for: %p \n",t);
      QUEUE h;
      QNODE *n;
      h = t->points;
      assert(h);
      n = Q_first(h);
      
      while(n){
	printf("\t x: %d y: %d z: %d ptr: %p\n",
	       n->e[X],
	       n->e[Y],
	       n->e[Z],
	       n->e);
	n = Q_next(h,n);
      }
    }
  }
}


//
// getTileLength finds the dimesions for a sub grid given a memory allocation
//
int getTileLength (double MEM){

  // At anyone time in memory one may need:
  // - 1 Grid with R points
  // - 1 TIN with R points and 2R Triangles
  // - 1 PQ with 2R elements
  //
  // Max memory for one tile:
  // Let R be number of points in tile and sqrt(R) be length of tile side
  // Let TL be tile length
  // MEM = 2R * (sizeOf(Triangle) + sizeOf(PQelement) + sizeOf(Point))
  //
  // TL = sqrt(MEM/(2*(sizeOf(Triangle) + sizeOf(PQelement) + sizeOf(Point))))

  double TL;
  MEM = MEM * 1048576.0; //convert MB into B
  TL = sqrt(MEM/
	    (2*((double)sizeof(TRIANGLE)+
		(double)sizeof(PQ_elemType)+
		(double)sizeof(POINT)))); 
  
  return TL;
}

///////////////////////////////////////////////////////////////////////////////
// areaSign - computes the signed area of a triangle, positive if the
// third point is off to the left of ab
//
int areaSign(COORD *a, COORD *b, COORD *c){
  
  // Rouding error for area computation. areaX should always be an integer
  double error = 0.5;

  double areaX;
  areaX = (b[X]-a[X]) * (double)(c[Y] - a[Y]) -
          (c[X]-a[X]) * (double)(b[Y] - a[Y]);

  // c is left of ab
  if ( areaX > error)
    return 1;
  // c is right of ab
  else if ( areaX < -error)
    return -1;
  // c is on ab
  else 
    return 0;
}

//
// inTri2d - is point z in triangle abc?
//
int inTri2D(COORD* a, COORD* b, COORD* c, COORD* z){
  
  int area0, area1, area2;
  
  // Assume no collinear triangles
  assert(areaSign(a, b, c) != 0);

  area0 = areaSign( a, b, z);
  area1 = areaSign( b, c, z);
  area2 = areaSign( c, a, z);
  
  if ( ((area0 == 0) && (area1 > 0) && (area2 > 0)) ||
       ((area1 == 0) && (area0 > 0) && (area2 > 0)) ||
       ((area2 == 0) && (area0 > 0) && (area1 > 0)) )
      // on an edge
      return 1;

    if ( ((area0 == 0) && (area1 < 0) && (area2 < 0)) ||
	 ((area1 == 0) && (area0 < 0) && (area2 < 0)) ||
	 ((area2 == 0) && (area0 < 0) && (area1 < 0)) )
      // on an edge
      return 1;

    if ( ((area0 < 0) && (area1 < 0) && (area2 < 0)) ||
	 ((area1 > 0) && (area0 > 0) && (area2 > 0)) )
      // inside tri
      return 1;
    
    if ( (area0 == 0) && (area1 == 0) && (area2 == 0)){
      // a,b,c,z are collinear
      printf("failure");
      assert(0);
      exit(1);
    }
    
    if ( ((area0 == 0) && ( area1 == 0)) ||
	 ((area0 == 0) && ( area2 == 0)) ||
	 ((area2 == 0) && ( area1 == 0)) )
      // z is either a, b, or c this should not happen as we remove points from
      // point lists when they become part of the triangle.
      return 1;

    else 
      return 0;
}


//
// Validate that all points in a triangle's point list are actually in
// that triangle
//
void checkPointList(TRIANGLE* t){
  if(t != NULL){
    if(t->points == NULL)
      printf("triangle %p is DONE\n",t);
    else{
      QNODE *h, *n;
      h = t->points;
      assert(h);
      n = Q_first(h);
      
      int badpoints=0, npoints=0;
      while(n !=  NULL){
	npoints++;
	assert(inTri2D(t->p1, t->p2, t->p3, n->e));
	if (!inTri2D(t->p1, t->p2, t->p3, n->e)) {
	  badpoints++;
	}
	//printf("Point %p in triangle %p\n",n->e,t);
	n = Q_next(h,n);
      }
      printf("triangle %p: npoints=%d, badpoints=%d maxE=%p \n",t, npoints, 
	     badpoints,t->maxE);
    }
  }
}

//
// Find the third point given two known points
//
COORD* findThirdPoint(POINT p1, POINT p2, POINT p3, POINT pa, POINT pb){
  if((pa == p1 && pb == p2) ||
     (pa == p2 && pb == p1))
    return p3;
  else if((pa == p1 && pb == p3) ||
     (pa == p3 && pb == p1))
    return p2;
  else
    return p1;

}


//
// Is a point an end point of a particular triangle
//
int isEndPoint(TRIANGLE* t, POINT p){
  if(t->p1 == p || t->p2 == p || t->p3 == p)
    return 1;
  else
    return 0;
}


//
// Find the edge opposite to a point p
//
EDGE findOpposite(TRIANGLE* t, POINT p){
  EDGE e;
  if(t->p1 == p){
    e.p1 = t->p2;
    e.p2 = t->p3;
    return e;
  }
  if(t->p2 == p){
    e.p1 = t->p1;
    e.p2 = t->p3;
    return e;
  }
  else{
    assert(t->p3 == p);
    e.p1 = t->p1;
    e.p2 = t->p2;
    return e;
  }
}


//
// Is a given edge in a triangle t?
//
int edgeInTriangle(TRIANGLE* t, EDGE e){ 
  if((t->p1 == e.p1 && t->p2 == e.p2) || 
     (t->p1 == e.p2 && t->p2 == e.p1) || // edge12
     (t->p1 == e.p1 && t->p3 == e.p2) || 
     (t->p1 == e.p2 && t->p3 == e.p1) || // edge13
     (t->p2 == e.p1 && t->p3 == e.p2) || 
     (t->p2 == e.p2 && t->p3 == e.p1))   // edge23
    return 1;
  else
    return 0;
}


//
// Do two edges have equal valued points?
//
int edgePointsEqual(EDGE e1, EDGE e2){
  if((e1.p1 == e2.p1 || e1.p1 == e2.p2) &&
     (e1.p2 == e2.p1 || e1.p2 == e2.p2) )
    return 1;
  else
    return 0;
}


//
// Validate that at least t1 and t2 are part of t. This is used to
// debug that the triangles being created inside t are actually inside
// triangle t
//
void triangleCheck(TRIANGLE* t,TRIANGLE* t1,TRIANGLE* t2,TRIANGLE* t3){
  assert(t1 != NULL && t2 != NULL);

  // Check s with 3 internal tris
  if(t3 != NULL){
    // Find center point
    COORD* cp;
    cp = (COORD*) malloc(3*sizeof(COORD));
    if (!isEndPoint(t,t1->p1))
      cp = t1->p1;
    if (!isEndPoint(t,t1->p2))
      cp = t1->p2;
    if (!isEndPoint(t,t1->p3))
      cp = t1->p3;

    // Find opposite edges to Center point
    EDGE e1,e2,e3;
    e1 = findOpposite(t1,cp);
    e2 = findOpposite(t2,cp);
    e3 = findOpposite(t3,cp);

    //for every t1, t2, t3 the edge opposite cp is an edge of t
    // Bad asserts!!!
    assert(edgeInTriangle(t, e1));
    assert(edgeInTriangle(t, e2));
    assert(edgeInTriangle(t, e3));
    // No two edges are the same
    assert(!edgePointsEqual(e1,e2) && !edgePointsEqual(e1,e3) && !edgePointsEqual(e2,e3));

    printf("\nTriangles %p %p %p in %p \n",t1,t2,t3,t);
  }
  else {
    // Find center point (the point on t1 not on s)
    COORD* cp;
    cp = (COORD*) malloc(3*sizeof(COORD));
    if (!isEndPoint(t,t1->p1))
      cp = t1->p1;
    if (!isEndPoint(t,t1->p2))
      cp = t1->p2;
    if (!isEndPoint(t,t1->p3))
      cp = t1->p3;

    // Find opposite edges to Center point
    EDGE e1,e2;
    e1 = findOpposite(t1,cp);
    e2 = findOpposite(t2,cp);

    //for every t1, t2, t3 the edge opposite cp is an edge of t
    assert(edgeInTriangle(t, e1));
    assert(edgeInTriangle(t, e2));
    // No two edges are the same
    assert(!edgePointsEqual(e1,e2));

    printf("\nTriangles %p %p in %p \n",t1,t2,t);
  }
  
}