/*
 * Program to take a model and produce another model, describing the
 *  same solid but with every triangle edge shared by exactly two
 *  triangles, with existing triangles split as necessary to accomplish
 *  this.
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include "3d.h"
#include "model.h"

typedef struct edge EDGE;
typedef struct tv TV;

struct tv {
  TRIANGLE *t;
  int v[2];
  } ;

struct edge {
  EDGE *link;
  int e[2];
  TV tv;
  } ;

static MODEL *m;
static EDGE *edges;

static PT3 xyz_to_pt3(XYZ) __attribute__((__const__));
static PT3 xyz_to_pt3(XYZ xyz)
{
 return((PT3){.x=xyz.x,.y=xyz.y,.z=xyz.z});
}

static void errf(const char *s)
{
 fprintf(stderr,"error in model, offset %lu:\n%s\n",ftell(stdin),s);
}

static void read_model(void)
{
 m = load_model(stdin,&errf);
 if (! m) exit(1);
}

static void dump_model(void)
{
 int i;
 int n;

 n = n_locations(m);
 for (i=0;i<n;i++)
  { printf("location %d: (%g,%g,%g)\n",i,get_location(m,i)->l.x,get_location(m,i)->l.y,get_location(m,i)->l.z);
  }
 n = n_points(m);
 for (i=0;i<n;i++)
  { printf("point %d: location %d (%g,%g,%g)\n", i,
	get_point(m,i)->loc,
	get_location(m,get_point(m,i)->loc)->l.x,
	get_location(m,get_point(m,i)->loc)->l.y,
	get_location(m,get_point(m,i)->loc)->l.z );
  }
 n = n_triangles(m);
 for (i=0;i<n;i++)
  { printf("triangle %d: %d-%d-%d (%g,%g,%g) - (%g,%g,%g) - (%g,%g,%g)\n", i,
	get_triangle(m,i)->corner[0],
	get_triangle(m,i)->corner[1],
	get_triangle(m,i)->corner[2],
	get_location(m,get_point(m,get_triangle(m,i)->corner[0])->loc)->l.x,
	get_location(m,get_point(m,get_triangle(m,i)->corner[0])->loc)->l.y,
	get_location(m,get_point(m,get_triangle(m,i)->corner[0])->loc)->l.z,
	get_location(m,get_point(m,get_triangle(m,i)->corner[1])->loc)->l.x,
	get_location(m,get_point(m,get_triangle(m,i)->corner[1])->loc)->l.y,
	get_location(m,get_point(m,get_triangle(m,i)->corner[1])->loc)->l.z,
	get_location(m,get_point(m,get_triangle(m,i)->corner[2])->loc)->l.x,
	get_location(m,get_point(m,get_triangle(m,i)->corner[2])->loc)->l.y,
	get_location(m,get_point(m,get_triangle(m,i)->corner[2])->loc)->l.z );
  }
}

static void queue_edge(int p1, int p2, TV tv)
{
 EDGE *e;

 e = malloc(sizeof(EDGE));
 e->e[0] = p1;
 e->e[1] = p2;
 e->tv = tv;
 e->link = edges;
 edges = e;
}

static void queue_edges(OBJECT *o)
{
 typedef struct e E;
 struct e {
   int p[2];
   int n;
   TV tv;
   } ;
 E *es;
 int tx;
 TRIANGLE *t;
 int n;
 int j;
 XYZ l;
 int fatal;

 static void save_edge(TRIANGLE *t, int v1, int v2)
  { int i;
    int e1;
    int e2;
    e1 = get_point(m,t->corner[v1])->loc;
    e2 = get_point(m,t->corner[v2])->loc;
    printf("edge (%g,%g,%g) - (%g,%g,%g): ",
	get_location(m,e1)->l.x,
	get_location(m,e1)->l.y,
	get_location(m,e1)->l.z,
	get_location(m,e2)->l.x,
	get_location(m,e2)->l.y,
	get_location(m,e2)->l.z );
    for (i=n-1;i>=0;i--)
     { if ( ((es[i].p[0] == e1) && (es[i].p[1] == e2)) ||
	    ((es[i].p[0] == e2) && (es[i].p[1] == e1)) )
	{ es[i].n ++;
	  printf("found %d, n now %d\n",i,es[i].n);
	  return;
	}
     }
    i = n++;
    printf("new %d\n",i);
    es[i].p[0] = e1;
    es[i].p[1] = e2;
    es[i].n = 1;
    es[i].tv.t = t;
    es[i].tv.v[0] = v1;
    es[i].tv.v[1] = v2;
  }

 fatal = 0;
 n = n_triangles(m) * 3;
 es = malloc(n*sizeof(*es));
 n = 0;
 for (tx=o->ntriangles-1;tx>=0;tx--)
  { t = get_triangle(m,o->triangles[tx]);
    save_edge(t,0,1);
    save_edge(t,1,2);
    save_edge(t,2,0);
  }
 for (j=n-1;j>=0;j--)
  { switch (es[j].n)
     { case 1:
	  printf("queueing edge %d (%g,%g,%g) - %d (%g,%g,%g): %d\n",
		es[j].p[0],
		get_location(m,es[j].p[0])->l.x,
		get_location(m,es[j].p[0])->l.y,
		get_location(m,es[j].p[0])->l.z,
		es[j].p[1],
		get_location(m,es[j].p[1])->l.x,
		get_location(m,es[j].p[1])->l.y,
		get_location(m,es[j].p[1])->l.z,
		es[j].n );
	  printf("  t %d (%g,%g,%g) - (%g,%g,%g) - (%g,%g,%g) [%d, %d]\n",
		es[j].tv.t->inx,
		get_location(m,get_point(m,es[j].tv.t->corner[0])->loc)->l.x,
		get_location(m,get_point(m,es[j].tv.t->corner[0])->loc)->l.y,
		get_location(m,get_point(m,es[j].tv.t->corner[0])->loc)->l.z,
		get_location(m,get_point(m,es[j].tv.t->corner[1])->loc)->l.x,
		get_location(m,get_point(m,es[j].tv.t->corner[1])->loc)->l.y,
		get_location(m,get_point(m,es[j].tv.t->corner[1])->loc)->l.z,
		get_location(m,get_point(m,es[j].tv.t->corner[2])->loc)->l.x,
		get_location(m,get_point(m,es[j].tv.t->corner[2])->loc)->l.y,
		get_location(m,get_point(m,es[j].tv.t->corner[2])->loc)->l.z,
		es[j].tv.v[0], es[j].tv.v[1] );
	  queue_edge(es[j].p[0],es[j].p[1],es[j].tv);
	  printf("  queued edge %p\n",(void *)edges);
	  break;
       case 2:
	  break;
       default:
	  fprintf(stderr,"edge %d-%d (",es[j].p[0],es[j].p[1]);
	  l = get_location(m,es[j].p[0])->l;
	  fprintf(stderr,"(%g,%g,%g)-",l.x,l.y,l.z);
	  l = get_location(m,es[j].p[1])->l;
	  fprintf(stderr,"(%g,%g,%g)): %d\n",l.x,l.y,l.z,es[j].n);
	  fatal = 1;
	  break;
     }
  }
 free(es);
 if (fatal) exit(1);
}

static int collinear_test(PT3 p0, PT3 v, int px, double *dp)
{
 PT3 p;
 double d;
 PT3 perpendicular;

 p = sub3(xyz_to_pt3(get_location(m,get_point(m,px)->loc)->l),p0);
 d = dot3(v,p);
 perpendicular = sub3(p,smul3(d,v));
 if (norm23(perpendicular) < 1e-16)
  { *dp = d;
    return(1);
  }
 return(0);
}

/*
 * Given an unshared edge, we look for other unshared edges having
 *  points collinear with it and split triangles as necessary to ensure
 *  edges are shared properly.
 */
static void process_edges(OBJECT *o)
{
 typedef struct ei EI;
 struct ei {
   EI *link;
   TV tv;
   double l;
   double h;
   } ;
 EDGE *e;
 EDGE *e2;
 EDGE **ep;
 double d1;
 double d2;
 EI *work;
 PT3 p;
 PT3 v;
 PT3 ep1;
 PT3 ep2;
 EI *ei;
 EI *ei2;
 EI **eip;
 EI *eipref[2];
 EI **eiprefp[2];

 void eisave(EDGE *e, double l, double h)
  { EI *ei;
    ei = malloc(sizeof(EI));
    ei->tv = e->tv;
    if (l <= h)
     { ei->l = l;
       ei->h = h;
     }
    else
     { ei->tv.v[0] = e->tv.v[1];
       ei->tv.v[1] = e->tv.v[0];
       ei->l = h;
       ei->h = l;
     }
    ei->link = work;
    work = ei;
    printf("eisave %p %g %g -> %p\n",(void *)e,l,h,(void *)ei);
  }

 void free_ei(EI *ei)
  { free(ei);
  }

 void dump_ei(FILE *to, EI *ei)
  { fprintf(to,"%p: t %p (%d), (%g,%g,%g) - (%g,%g,%g) - (%g,%g,%g)\n",
	(void *)ei, (void *)ei->tv.t, ei->tv.t->inx,
	get_location(m,get_point(m,ei->tv.t->corner[0])->loc)->l.x,
	get_location(m,get_point(m,ei->tv.t->corner[0])->loc)->l.y,
	get_location(m,get_point(m,ei->tv.t->corner[0])->loc)->l.z,
	get_location(m,get_point(m,ei->tv.t->corner[1])->loc)->l.x,
	get_location(m,get_point(m,ei->tv.t->corner[1])->loc)->l.y,
	get_location(m,get_point(m,ei->tv.t->corner[1])->loc)->l.z,
	get_location(m,get_point(m,ei->tv.t->corner[2])->loc)->l.x,
	get_location(m,get_point(m,ei->tv.t->corner[2])->loc)->l.y,
	get_location(m,get_point(m,ei->tv.t->corner[2])->loc)->l.z );
    fprintf(to,"  v1 %d v2 %d l %g h %g\n",ei->tv.v[0],ei->tv.v[1],ei->l,ei->h);
  }

 EI *ei_insert(EI *list, EI *el)
  { EI **lp;
    EI *t;
    lp = &list;
    while ((t = *lp))
     { if (t->l >= el->l) break;
       lp = &t->link;
     }
    el->link = t;
    *lp = el;
    return(list);
  }

 EI *eisort(EI *list)
  { EI *a;
    EI *b;
    EI *t;
    EI **lp;
    if (!list || !list->link) return(list);
    a = 0;
    b = 0;
    while (list)
     { t = list;
       list = t->link;
       t->link = a;
       a = b;
       b = t;
     }
    a = eisort(a);
    b = eisort(b);
    lp = &list;
    while (1)
     { if (a && (!b || (a->l < b->l)))
	{ t = a;
	  a = t->link;
	}
       else if (b)
	{ t = b;
	  b = t->link;
	}
       else
	{ break;
	}
       *lp = t;
       lp = &t->link;
     }
    *lp = 0;
    return(list);
  }

 /*
  * Split ei.  The low part is split off and becomes a new triangle;
  *  the existing ei is updated to describe the remainder.  The split
  *  occurs a one end of ei2; the low end if end2=0, high if end2=1.
  *  ei2 is not modified.
  */
 void split_triangle(EI *ei, EI *ei2, int end2)
  { int t1rep;
    int t2rep;
    int spt;
    int ntx;
    TRIANGLE *nt;
    printf("split_triangle:\n");
    printf("end2 %d\n",end2);
    dump_ei(stdout,ei);
    t1rep = ei->tv.t->corner[ei->tv.v[1]];
    t2rep = ei->tv.t->corner[ei->tv.v[0]];
    printf("new replace [%d]=%d (%g,%g,%g)\n", ei->tv.v[1], t1rep,
	get_location(m,t1rep)->l.x,
	get_location(m,t1rep)->l.y,
	get_location(m,t1rep)->l.z );
    printf("old replace [%d]=%d (%g,%g,%g)\n", ei->tv.v[0], t2rep,
	get_location(m,t2rep)->l.x,
	get_location(m,t2rep)->l.y,
	get_location(m,t2rep)->l.z );
    dump_ei(stdout,ei2);
    spt = ei2->tv.t->corner[ei2->tv.v[end2]];
    printf("with [%d]=%d (%g,%g,%g)\n", ei2->tv.v[end2], spt,
	get_location(m,spt)->l.x,
	get_location(m,spt)->l.y,
	get_location(m,spt)->l.z );
    ntx = new_triangle(m);
    nt = get_triangle(m,ntx);
    *nt = *ei->tv.t;
    nt->inx = ntx;
    nt->corner[ei->tv.v[1]] = spt;
    ei->tv.t->corner[ei->tv.v[0]] = spt;
    o->triangles = realloc(o->triangles,(o->ntriangles+1)*sizeof(*o->triangles));
    o->triangles[o->ntriangles++] = ntx;
    dump_model();
    exit(0);
  }

 work = 0;
 while (edges)
  { e = edges;
    edges = e->link;
    printf("considering edge %p\n",(void *)e);
    p = xyz_to_pt3(get_location(m,get_point(m,e->e[0])->loc)->l);
    ep1 = sub3(xyz_to_pt3(get_location(m,get_point(m,e->e[1])->loc)->l),p);
    d1 = norm3(ep1);
    v = smul3(1/d1,ep1);
    printf("p=(%g,%g,%g) v=(%g,%g,%g)\n",p.x,p.y,p.z,v.x,v.y,v.z);
    eisave(e,0,d1);
    ep = &edges;
    while ((e2 = *ep))
     { if (collinear_test(p,v,e2->e[0],&d1) && collinear_test(p,v,e2->e[1],&d2))
	{ *ep = e2->link;
	  eisave(e2,d1,d2);
	}
       else
	{ ep = &e2->link;
	}
     }
    work = eisort(work);
    printf("work list:\n");
    for (ei=work;ei;ei=ei->link) printf("  %p\n",(void *)ei);
    while (1)
     { printf("loop, work %p\n",(void *)work);
       if (! work) break;
       ei = work;
       ei2 = ei->link;
       printf("ei %p ei2 %p\n",(void *)ei,(void *)ei2);
       if (!ei2 || (ei2->l >= ei->h))
	{ // No second triangle, or it doesn't overlap the first
	  // Just complain and drop first triangle
	  fprintf(stderr,"Unmatched triangle edge: ");
	  ep1 = add3(p,smul3(ei->l,v));
	  fprintf(stderr,"(%g,%g,%g) - ",ep1.x,ep1.y,ep1.z);
	  ep2 = add3(p,smul3(ei->h,v));
	  fprintf(stderr,"(%g,%g,%g)\n",ep2.x,ep2.y,ep2.z);
	  free_ei(ei);
	  work = ei2;
	  continue;
	}
       if (fabs(ei2->l-ei->l) > 1e-8)
	{ // Second triangle's start point doesn't match the first's
	  // Split off unmatched part of first triangle.
	  fprintf(stderr,"Unmatched triangle edge: ");
	  ep1 = add3(p,smul3(ei->l,v));
	  fprintf(stderr,"(%g,%g,%g) - ",ep1.x,ep1.y,ep1.z);
	  ep2 = add3(p,smul3(ei2->l,v));
	  fprintf(stderr,"(%g,%g,%g)\n",ep2.x,ep2.y,ep2.z);
	  dump_ei(stdout,ei);
	  dump_ei(stdout,ei2);
	  split_triangle(ei,ei2,0);
	  continue;
	}
       // First two triangles have identical start points.
       // Find the two with the lowest end points.
       if (ei->h < ei2->h)
	{ eipref[0] = ei;
	  eiprefp[0] = &work;
	  eipref[1] = ei2;
	  eiprefp[1] = &ei->link;
	}
       else
	{ eipref[1] = ei;
	  eiprefp[1] = &work;
	  eipref[0] = ei2;
	  eiprefp[0] = &ei->link;
	}
       eip = &ei2->link;
       while ((ei2 = *eip))
	{ if (fabs(ei2->l-ei->l) > 1e-8) break;
	  if (ei2->h < eipref[0]->h)
	   { eipref[1] = eipref[0];
	     eiprefp[1] = eiprefp[0];
	     eipref[0] = ei2;
	     eiprefp[0] = eip;
	   }
	  else if (ei2->h < eipref[1]->h)
	   { eipref[1] = ei2;
	     eiprefp[1] = eip;
	   }
	}
       printf("eipref %p %p\n",(void *)eipref[0],(void *)eipref[1]);
       // Remove those two triangles from the list.
       *eiprefp[0] = eipref[0]->link;
       if (eiprefp[1] == &eipref[0]->link) eiprefp[1] = eiprefp[0];
       *eiprefp[1] = eipref[1]->link;
       // Are the end points identical?
       if (fabs(eipref[0]->h-eipref[1]->h) < 1e-8)
	{ // Yes.  Just drop both triangles and carry on.
	  printf("h match, dropping and done\n");
	  free_ei(eipref[0]);
	  free_ei(eipref[1]);
	}
       else
	{ // No.  Split the one with the higher endpoint, link the
	  //  remainder back in, drop the other one, and continue.
	  printf("[0] %g..%g [1] %g..%g\n",eipref[0]->l,eipref[0]->h,eipref[1]->l,eipref[1]->h);
	  printf("p1=%d (%g,%g,%g) p2=%d (%g,%g,%g), t %p (%d) v1 %d v2 %d\n",
		e->e[0],
		get_location(m,get_point(m,e->e[0])->loc)->l.x,
		get_location(m,get_point(m,e->e[0])->loc)->l.y,
		get_location(m,get_point(m,e->e[0])->loc)->l.z,
		e->e[1],
		get_location(m,get_point(m,e->e[1])->loc)->l.x,
		get_location(m,get_point(m,e->e[1])->loc)->l.y,
		get_location(m,get_point(m,e->e[1])->loc)->l.z,
		(void *)e->tv.t, e->tv.t->inx, e->tv.v[0], e->tv.v[1] );
	  dump_ei(stdout,eipref[0]);
	  dump_ei(stdout,eipref[1]);
	  split_triangle(eipref[1],eipref[0],1);
	  work = ei_insert(work,eipref[1]);
	  free_ei(eipref[0]);
	}
     }
  }
}

static void sub_triangulate(void)
{
 int ox;
 OBJECT *o;

 for (ox=n_objects(m)-1;ox>=0;ox--)
  { o = get_object(m,ox);
    queue_edges(o);
    process_edges(o);
  }
}

int main(void);
int main(void)
{
 read_model();
 dump_model();
 sub_triangulate();
// dump_model();
 return(0);
}