/* triangle.c - Don Yang (uguu.org)
   Triangle / vertex management routines.

   07/15/01
                                                                           */

#include<assert.h>
#include<stdio.h>
#include<stdlib.h>
#include"triangle.h"


/* Global edge/vertex numbering indexes */
int Edge1[3] = {1, 0, 0};
int Edge2[3] = {2, 2, 1};

/* Local functions */
static int AddReference(Triangle *t);
static int Circumscribe(Vertex *v[3], double *cx, double *cy, double *cr2);
static void DeleteReference(Triangle *t);

/************************************************************* AddTriangle */
int AddTriangle(Triangle **root, Triangle **n, Vertex *v[3])
{
   Triangle *t, **node;
   double cx, cy, cr2;

   assert(root);

   /* Create new triangle */
   if( (t = (Triangle*)malloc(sizeof(Triangle))) == NULL )
      return 1;

   /* Create circumcircle */
   if( Circumscribe(v, &cx, &cy, &cr2) )
   {
      /* No triangle is added, but operation is still successful */
      free(t);
      if( n )
         *n = NULL;
      return 0;
   }

   /* Add triangle to BST */
   t->v[0] = v[0];   t->cx = cx;
   t->v[1] = v[1];   t->cy = cy;
   t->v[2] = v[2];   t->cr2 = cr2;
   t->left = t->right = t->n[0] = t->n[1] = t->n[2] = NULL;
   t->key = v[0]->x + v[1]->x + v[2]->x;
   t->id = t->flags = 0;

   /* Add reference to triangle */
   if( AddReference(t) )
      return 1;

   for(node = root; *node;)
   {
      node = (t->key > (*node)->key)
         ?  &((*node)->right)
         :  &((*node)->left);
   }
   *node = t;

   if( n )
      *n = t;
   return 0;

} /* AddTriangle() */

/*************************************************************** AddVertex */
int AddVertex(Vertex **list, Vertex **v, double x, double y)
{
   Vertex *node;

   assert(list);

   if( (node = (Vertex*)malloc(sizeof(Vertex))) == NULL )
      return 1;

   node->x = x;
   node->y = y;
   node->next = *list;
   node->tlist = NULL;
   *list = node;

   if( v )
      *v = node;
   return 0;

} /* AddVertex() */

/******************************************************************* Apply */
void Apply(Triangle *root, void (*op)(Triangle*))
{
   while( root )
   {
      Apply(root->left, op);
      op(root);
      root = root->right;
   }

} /* Apply() */

/********************************************************** AssignVertexID */
void AssignVertexID(Vertex *list)
{
   Vertex *node;
   int id;

   id = 0;
   for(node = list; node; node = node->next)
      node->id = node->tlist ? id++ : -1;

} /* AssignVertexID() */

/********************************************************** CountTriangles */
int CountTriangles(Triangle *root)
{
   int count;

   count = 0;
   while( root )
   {
      count += CountTriangles(root->left) + 1;
      root = root->right;
   }
   return count;

} /* CountTriangles() */

/********************************************************** DeleteTriangle */
void DeleteTriangle(Triangle **root, Triangle *d)
{
   Triangle **node, *parent, *newroot;

   assert(root);
   assert(d);

   /* Find node to be deleted */
   for(node = root; *node != d;)
   {
      node = (d->key > (*node)->key)
         ?  &((*node)->right)
         :  &((*node)->left);
   }

   /* Move subtree */
   if( d->left == NULL )
   {
      /*   X          O
            \    ->  ...
             O
            ...          */
      *node = d->right;
   }
   else if( d->right == NULL )
   {
      /*   X        O
          /    ->  ...
         O
        ...             */
      *node = d->left;
   }
   else
   {
      parent = NULL;
      for(newroot = d->left; newroot->right; newroot = newroot->right)
         parent = newroot;
      /*      parent != NULL                parent == NULL

             X               N             X             N
            / \             /             / \           / \
           O   ...         O             N   ...  ->  ...  ...
          / \             / \           /
       ...  ...    ->  ...  ...       ...
              \               \
               N               S
              /               ...
             S
            ...                                                */
      *node = newroot;
      newroot->right = d->right;
      if( parent )
      {
         parent->right = newroot->left;
         newroot->left = d->left;
      }
   }

   /* Remove references to this triangle */
   DeleteReference(d);

   free(d);

} /* DeleteTriangle() */

/****************************************************** DeleteTriangleList */
void DeleteTriangleList(Triangle **root)
{
   assert(root);

   if( *root == NULL )
      return;
   DeleteTriangleList(&((*root)->left));
   DeleteTriangleList(&((*root)->right));

   /* Remove references */
   DeleteReference(*root);

   /* Delete triangle */
   free(*root);
   *root = NULL;

} /* DeleteTriangleList() */

/************************************************************ DeleteVertex */
void DeleteVertex(Vertex **list, Vertex *d)
{
   Vertex **node, *del;

   assert(list);

   for(node = list; *node; node = &((*node)->next))
   {
      if( *node == d )
         break;
   }
   if( *node )
   {
      del = *node;

      /* There should be no triangles that are using this vertex */
      assert(del->tlist == NULL);

      *node = del->next;
      free(del);
   }

} /* DeleteVertex() */

/******************************************************** DeleteVertexList */
void DeleteVertexList(Vertex **list)
{
   Vertex *node, *next;

   assert(list);

   for(node = *list; node; node = next)
   {
      /* All triangles should be deleted before deleting the vertices,
         no triangle should be referencing this vertex.                 */
      assert(node->tlist == NULL);

      next = node->next;
      free(node);
   }
   *list = NULL;

} /* DeleteVertexList() */

/**************************************************************** Enclosed */
int Enclosed(Vertex *v[3], Vertex *point)
{
   double d, dx1, dx2, dy1, dy2, s, t;

   /* (x, y) = (x0, y0) + s * (x1 - x0, y1 - y0) + t * (x2 - x0, y2 - y0)

      Solve for s and t, point is in triangle if s>=0 && t>=0 && s+t<=1
      (point can be located from (x0, y0) using two edges as basis vectors,
      the parameters s=0, t=0, and s+t=1 specifies the three edges).       */
   d = (dx1 = v[1]->x - v[0]->x) * (dy2 = v[2]->y - v[0]->y)
     - (dy1 = v[1]->y - v[0]->y) * (dx2 = v[2]->x - v[0]->x);
   if( d == 0 )
      return 0;
   s = ((point->x - v[0]->x) * dy2 - (point->y - v[0]->y) * dx2) / d;
   t = ((point->y - v[0]->y) * dx1 - (point->x - v[0]->x) * dy1) / d;

   return (s >= 0) && (t >= 0) && (s + t <= 1);

} /* Enclosed() */


/************************************************************ AddReference */
static int AddReference(Triangle *t)
{
   TList *ref;
   int i;

   assert(t);

   for(i = 0; i < 3; i++)
   {
      assert(t->v[i]);

      if( (ref = (TList*)malloc(sizeof(TList))) == NULL )
         return 1;
      ref->t = t;
      ref->next = t->v[i]->tlist;
      t->v[i]->tlist = ref;
   }

   return 0;

} /* AddReference() */

/************************************************************ Circumscribe */
static int Circumscribe(Vertex *v[3], double *cx, double *cy, double *cr2)
{
   double dx1, dy1, dx2, dy2, d, a, b, x2, y2;

   /* Reject degenerate triangles */
   dx1 = v[1]->x - v[2]->x;
   dy1 = v[1]->y - v[2]->y;
   if( (dx1 * dx1 + dy1 * dy1) < MIN_VERTEX_DISTANCE_2 )
      return 1;
   dx1 = v[1]->x - v[0]->x;
   dy1 = v[1]->y - v[0]->y;
   if( (dx1 * dx1 + dy1 * dy1) < MIN_VERTEX_DISTANCE_2 )
      return 1;
   dx2 = v[2]->x - v[0]->x;
   dy2 = v[2]->y - v[0]->y;
   if( (dx2 * dx2 + dy2 * dy2) < MIN_VERTEX_DISTANCE_2 )
      return 1;

   /* Find circumcenter (point equidistance to all three vertices)
         (xc - x0)^2 + (yc - y0)^2 = (xc - x1)^2 + (yc - y1)^2

         xc^2 - 2*x0*xc + x0^2 + yc^2 - 2*y0*yc + y0^2
      =  xc^2 - 2*x1*xc + x1^2 + yc^2 - 2*y1*yc + y1^2

         -2*x0*xc + x0^2 + -2*y0*yc + y0^2
      =  -2*x1*xc + x1^2 + -2*y1*yc + y1^2

      2 * (x1 - x0) * xc + 2 * (y1 - y0) * yc = x0^2 - x1^2 + y0^2 - y1^2
      2 * (x2 - x0) * xc + 2 * (y2 - y0) * yc = x0^2 - x2^2 + y0^2 - y2^2

      dx1 = x1 - x0        dx2 = x2 - x0     a = x0^2 - x1^2 + y0^2 - y1^2
      dy1 = y1 - y0        dy2 = y2 - y0     b = x0^2 - x2^2 + y0^2 - y2^2

      d = 2*dx1 * 2*dy2 - 2*dx2 * 2*dy1 = 4 * (dx1 * dy2 - dx2 * dy1)
      xc = (a * 2*dy2 - b * 2*dy1) / d
      yc = (2*dx1 * b - 2*dx2 * a) / d
                                                                           */
   if( (d = 2 * (dx2 * dy1 - dx1 * dy2)) == 0 )
      return 1;
   a = (x2 = v[0]->x * v[0]->x) - v[1]->x * v[1]->x
     + (y2 = v[0]->y * v[0]->y) - v[1]->y * v[1]->y;
   b = x2 - v[2]->x * v[2]->x + y2 - v[2]->y * v[2]->y;
   *cx = (a * dy2 - b * dy1) / d;
   *cy = (dx1 * b - dx2 * a) / d;

   /* Calculate radius */
   dx2 = v[0]->x - *cx;
   dy2 = v[0]->y - *cy;
   *cr2 = dx2 * dx2 + dy2 * dy2;
   return 0;

} /* Circumscribe() */

/********************************************************* DeleteReference */
static void DeleteReference(Triangle *t)
{
   TList **ref, *del;
   int i;

   assert(t);

   for(i = 0; i < 3; i++)
   {
      assert(t->v[i]);

      for(ref = &(t->v[i]->tlist); *ref; ref = &((*ref)->next))
      {
         if( (*ref)->t == t )
            break;
      }
      if( *ref )
      {
         del = *ref;
         *ref = del->next;
         free(del);
      }
   }

} /* DeleteReference() */