/*
 * 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"

// Define this to get internal traces on stderr.
//#define VERBOSE

typedef struct edge EDGE;
typedef struct tv TV;
typedef struct px PX;
typedef struct lx LX;
typedef struct tcx TCX;

struct tcx {
  int tc;
  } ;

struct px {
  int p;
  } ;

struct lx {
  int l;
  } ;

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

struct edge {
  EDGE *link;
  LX 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 POINT *gpoint(PX x)
{
 return(get_point(m,x.p));
}

static LX pointloc(POINT *p)
{
 return((LX){.l=p->loc});
}

static LOCATION *glocation(LX x)
{
 return(get_location(m,x.l));
}

static PX tcorner(TRIANGLE *t, TCX tc)
{
 return((PX){.p=t->corner[tc.tc]});
}

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);
}

#ifdef VERBOSE
static void dump_model(void)
{
 int i;
 int n;

 n = n_locations(m);
 for (i=0;i<n;i++)
  { fprintf(stderr,"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++)
  { fprintf(stderr,"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++)
  { fprintf(stderr,"triangle %d: %d-%d-%d (%g,%g,%g) - (%g,%g,%g) - (%g,%g,%g)\n", i,
	tcorner(get_triangle(m,i),(TCX){0}).p,
	tcorner(get_triangle(m,i),(TCX){1}).p,
	tcorner(get_triangle(m,i),(TCX){2}).p,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){0}).p)->loc)->l.x,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){0}).p)->loc)->l.y,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){0}).p)->loc)->l.z,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){1}).p)->loc)->l.x,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){1}).p)->loc)->l.y,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){1}).p)->loc)->l.z,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){2}).p)->loc)->l.x,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){2}).p)->loc)->l.y,
	get_location(m,get_point(m,tcorner(get_triangle(m,i),(TCX){2}).p)->loc)->l.z );
  }
}
#endif

static void queue_edge(LX p1, LX 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 {
   LX l[2];
   int n;
   TV tv;
   } ;
 E *es;
 int tx;
 TRIANGLE *t;
 int n;
 int j;
 int fatal;

 static void save_edge(TRIANGLE *t, TCX v1, TCX v2)
  { int i;
    LX e1;
    LX e2;
    e1 = pointloc(gpoint(tcorner(t,v1)));
    e2 = pointloc(gpoint(tcorner(t,v2)));
#ifdef VERBOSE
    fprintf(stderr,"edge (%g,%g,%g) - (%g,%g,%g): ",
	glocation(e1)->l.x,
	glocation(e1)->l.y,
	glocation(e1)->l.z,
	glocation(e2)->l.x,
	glocation(e2)->l.y,
	glocation(e2)->l.z );
#endif
    for (i=n-1;i>=0;i--)
     { if ( ((es[i].l[0].l == e1.l) && (es[i].l[1].l == e2.l)) ||
	    ((es[i].l[0].l == e2.l) && (es[i].l[1].l == e1.l)) )
	{ es[i].n ++;
#ifdef VERBOSE
	  fprintf(stderr,"found %d, n now %d\n",i,es[i].n);
#endif
	  return;
	}
     }
    i = n++;
#ifdef VERBOSE
    fprintf(stderr,"new %d\n",i);
#endif
    es[i].l[0] = e1;
    es[i].l[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,(TCX){0},(TCX){1});
    save_edge(t,(TCX){1},(TCX){2});
    save_edge(t,(TCX){2},(TCX){0});
  }
 for (j=n-1;j>=0;j--)
  { switch (es[j].n)
     { case 1:
#ifdef VERBOSE
	  fprintf(stderr,"queueing edge %d (%g,%g,%g) - %d (%g,%g,%g): %d\n",
		es[j].l[0].l,
		glocation(es[j].l[0])->l.x,
		glocation(es[j].l[0])->l.y,
		glocation(es[j].l[0])->l.z,
		es[j].l[1].l,
		glocation(es[j].l[1])->l.x,
		glocation(es[j].l[1])->l.y,
		glocation(es[j].l[1])->l.z,
		es[j].n );
	  fprintf(stderr,"  t %d (%g,%g,%g) - (%g,%g,%g) - (%g,%g,%g) [%d, %d]\n",
		es[j].tv.t->inx,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){0}))))->l.x,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){0}))))->l.y,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){0}))))->l.z,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){1}))))->l.x,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){1}))))->l.y,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){1}))))->l.z,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){2}))))->l.x,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){2}))))->l.y,
		glocation(pointloc(gpoint(tcorner(es[j].tv.t,(TCX){2}))))->l.z,
		es[j].tv.v[0].tc, es[j].tv.v[1].tc );
#endif
	  queue_edge(es[j].l[0],es[j].l[1],es[j].tv);
#ifdef VERBOSE
	  fprintf(stderr,"  queued edge %p\n",(void *)edges);
#endif
	  break;
       case 2:
	  break;
       default:
	  fprintf(stderr,"edge %d-%d (%g,%g,%g)-(%g,%g,%g): %d\n",
		es[j].l[0].l, es[j].l[1].l,
		glocation(es[j].l[0])->l.x,
		glocation(es[j].l[0])->l.y,
		glocation(es[j].l[0])->l.z,
		glocation(es[j].l[1])->l.x,
		glocation(es[j].l[1])->l.y,
		glocation(es[j].l[1])->l.z,
		es[j].n );
	  fatal = 1;
	  break;
     }
  }
 free(es);
 if (fatal) exit(1);
}

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

 p = sub3(xyz_to_pt3(glocation(lx)->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;
#ifdef VERBOSE
    fprintf(stderr,"eisave %p %g %g -> %p\n",(void *)e,l,h,(void *)ei);
#endif
  }

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

#ifdef VERBOSE
 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,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){0}))))->l.x,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){0}))))->l.y,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){0}))))->l.z,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){1}))))->l.x,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){1}))))->l.y,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){1}))))->l.z,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){2}))))->l.x,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){2}))))->l.y,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){2}))))->l.z );
    fprintf(to,"  v1 %d v2 %d l %g h %g\n",ei->tv.v[0].tc,ei->tv.v[1].tc,ei->l,ei->h);
  }
#endif

 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)
  { PX t1rep;
    PX t2rep;
    PX spt;
    int ntx;
    TRIANGLE *nt;
#ifdef VERBOSE
    fprintf(stderr,"split_triangle:\n");
    fprintf(stderr,"end2 %d\n",end2);
    dump_ei(stderr,ei);
#endif
    t1rep = tcorner(ei->tv.t,ei->tv.v[1]);
    t2rep = tcorner(ei->tv.t,ei->tv.v[0]);
#ifdef VERBOSE
    fprintf(stderr,"new replace [%d]=%d (%g,%g,%g)\n", ei->tv.v[1].tc, t1rep.p,
	glocation(pointloc(gpoint(t1rep)))->l.x,
	glocation(pointloc(gpoint(t1rep)))->l.y,
	glocation(pointloc(gpoint(t1rep)))->l.z );
    fprintf(stderr,"old replace [%d]=%d (%g,%g,%g)\n", ei->tv.v[0].tc, t2rep.p,
	glocation(pointloc(gpoint(t2rep)))->l.x,
	glocation(pointloc(gpoint(t2rep)))->l.y,
	glocation(pointloc(gpoint(t2rep)))->l.z );
    dump_ei(stderr,ei2);
#endif
    spt = tcorner(ei2->tv.t,ei2->tv.v[end2]);
#ifdef VERBOSE
    fprintf(stderr,"with [%d]=%d (%g,%g,%g)\n", ei2->tv.v[end2].tc, spt.p,
	glocation(pointloc(gpoint(spt)))->l.x,
	glocation(pointloc(gpoint(spt)))->l.y,
	glocation(pointloc(gpoint(spt)))->l.z );
#endif
    ntx = new_triangle(m);
    nt = get_triangle(m,ntx);
    *nt = *ei->tv.t;
    nt->inx = ntx;
    nt->corner[ei->tv.v[1].tc] = spt.p;
    ei->tv.t->corner[ei->tv.v[0].tc] = spt.p;
    ei->l = end2 ? ei2->h : ei2->l;
    o->triangles = realloc(o->triangles,(o->ntriangles+1)*sizeof(*o->triangles));
    o->triangles[o->ntriangles++] = ntx;
#ifdef VERBOSE
    fprintf(stderr,"new triangle %d = (%g,%g,%g) - (%g,%g,%g) - (%g,%g,%g)\n",
	ntx,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){0}))))->l.x,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){0}))))->l.y,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){0}))))->l.z,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){1}))))->l.x,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){1}))))->l.y,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){1}))))->l.z,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){2}))))->l.x,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){2}))))->l.y,
	glocation(pointloc(gpoint(tcorner(nt,(TCX){2}))))->l.z );
    fprintf(stderr,"old triangle %d now = (%g,%g,%g) - (%g,%g,%g) - (%g,%g,%g)\n",
	ei->tv.t->inx,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){0}))))->l.x,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){0}))))->l.y,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){0}))))->l.z,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){1}))))->l.x,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){1}))))->l.y,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){1}))))->l.z,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){2}))))->l.x,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){2}))))->l.y,
	glocation(pointloc(gpoint(tcorner(ei->tv.t,(TCX){2}))))->l.z );
#endif
  }

 work = 0;
 while (edges)
  { e = edges;
    edges = e->link;
#ifdef VERBOSE
    fprintf(stderr,"considering edge %p\n",(void *)e);
#endif
    p = xyz_to_pt3(glocation(e->e[0])->l);
    ep1 = sub3(xyz_to_pt3(glocation(e->e[1])->l),p);
    d1 = norm3(ep1);
    v = smul3(1/d1,ep1);
#ifdef VERBOSE
    fprintf(stderr,"p=(%g,%g,%g) v=(%g,%g,%g)\n",p.x,p.y,p.z,v.x,v.y,v.z);
#endif
    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);
#ifdef VERBOSE
    fprintf(stderr,"work list:\n");
    for (ei=work;ei;ei=ei->link) fprintf(stderr,"  %p\n",(void *)ei);
#endif
    while (1)
     {
#ifdef VERBOSE
       fprintf(stderr,"loop, work %p\n",(void *)work);
#endif
       if (! work) break;
       ei = work;
       ei2 = ei->link;
#ifdef VERBOSE
       fprintf(stderr,"ei %p ei2 %p\n",(void *)ei,(void *)ei2);
#endif
       if (!ei2 || (ei2->l >= ei->h))
	{ // No second triangle, or it doesn't overlap the first
	  // Just complain and drop first triangle
	  ep1 = add3(p,smul3(ei->l,v));
	  ep2 = add3(p,smul3(ei->h,v));
	  fprintf(stderr,"unmatched triangle edge: (%g,%g,%g) - (%g,%g,%g)\n",
		ep1.x, ep1.y, ep1.z, 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.
	  ep1 = add3(p,smul3(ei->l,v));
	  ep2 = add3(p,smul3(ei2->l,v));
	  fprintf(stderr,"unmatched triangle edge part: (%g,%g,%g) - (%g,%g,%g)\n",
		ep1.x, ep1.y, ep1.z, ep2.x, ep2.y, ep2.z );
#ifdef VERBOSE
	  dump_ei(stderr,ei);
	  dump_ei(stderr,ei2);
#endif
	  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;
	   }
	}
#ifdef VERBOSE
       fprintf(stderr,"eipref %p %p\n",(void *)eipref[0],(void *)eipref[1]);
#endif
       // 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.
#ifdef VERBOSE
	  fprintf(stderr,"h match, dropping and done\n");
#endif
	  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.
#ifdef VERBOSE
	  fprintf(stderr,"[0] %g..%g [1] %g..%g\n",eipref[0]->l,eipref[0]->h,eipref[1]->l,eipref[1]->h);
	  fprintf(stderr,"p1=%d (%g,%g,%g) p2=%d (%g,%g,%g), t %p (%d) v1 %d v2 %d\n",
		e->e[0].l,
		glocation(e->e[0])->l.x,
		glocation(e->e[0])->l.y,
		glocation(e->e[0])->l.z,
		e->e[1].l,
		glocation(e->e[1])->l.x,
		glocation(e->e[1])->l.y,
		glocation(e->e[1])->l.z,
		(void *)e->tv.t, e->tv.t->inx, e->tv.v[0].tc, e->tv.v[1].tc );
	  dump_ei(stderr,eipref[0]);
	  dump_ei(stderr,eipref[1]);
#endif
	  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();
#ifdef VERBOSE
 dump_model();
#endif
 sub_triangulate();
#ifdef VERBOSE
 dump_model();
#endif
 save_model(m,stdout);
 free_model(m);
 return(0);
}