/* See SPTriang.h. */
/* Last edited on 2023-02-12 07:55:17 by stolfi */

#include <SPTriang.h>
#include <SPQuad.h>
#include <SPBasic.h>
#include <SPIntegral.h>

#include <vec.h>
#include <r3.h>
#include <r3x3.h>
#include <fget.h>
#include <nget.h>
#include <filefmt.h>
#include <affirm.h>
#include <nat.h>
#include <bool.h>

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

#define NullArc SPQuad_NullArc
#define Splice SPQuad_Splice
#define MakeEdge SPQuad_MakeEdge

Arc_vec_t ArcsFromEdges(Arc_vec_t edge);
  /* Given an array {edge.e[..]} with one arc from each undirected 
    primal edge, returns an array {res} with two symmetric arcs for
    each edge.  The arcs for {edge.e[i]} are {res.e[2*i]} and 
    {res.e[2*i+1]}, in no particular order. */

#define Triangulation_FileFormat "1997-12-15"

Arc_vec_t RenumberVertices(Arc_vec_t arc);
  /* Renumbers the origin vertices of the arcs in {arc[..]} 
    consecutively from 0.  Returns a vector where element {i}
    is one arc out of vertex number {i}. */

Sign TriangSegmSide(R3Point *u, R3Point *v,  R3Point *p, R3Point *q, R3Point *r);
  /* TRUE if vectors {p,q,r} are all on the non-positive side
    of the plane through the origin generated by vectors {u} and {v}. */

void Project2(R3Point *u, R3Point *v, R3Point *p, double *au, double *av);
  /* Computes {au} and {av} such that {au*u + av*v} is the orthogonal
    projection of {p} on the plane spanned by {u} and {v}.
    Assumes that {u} and {v} are linearly independent. */

R3Point FourCenter(Arc e);
  /* Computes the barycenter of the four vertices in {Left(e)} and
     {Right(e)}. */

nat_t OrgDegree(Arc e)
  { Arc t = e;
    nat_t n = 0;
    do 
      { n++; t = Onext(t); }
    while (t != e);
    return n;
  }

nat_t LeftDegree(Arc e)
  { Arc t = e;
    nat_t n = 0;
    do 
      { n++; t = Lnext(t); }
    while (t != e);
    return n;
  }

Site *Org (Arc e)
  { return Odata(e); }
  
Site *Dest (Arc e)
  { return Ddata(e); }

Face *Left (Arc e)
  { return Ldata(e); }
  
Face *Right (Arc e)
  { return Rdata(e); }
  
void SPTriang_SwapEdge(Arc e)
  { Arc a = Lnext(Sym(e));
    Arc b = Lnext(e);
    /* Required conditions: */
    affirm((Onext(e) != e) && (Dnext(e) != e), "edge is not swappable");
    affirm(b != e, "Left(e) is a monogon");
    affirm(Lnext(Lnext(b)) == e, "Left(e) is not a triangle");
    affirm(a != e, "Right(e) is a monogon");
    affirm(Lnext(Lnext(a)) == Sym(e), "Right(e) is not a triangle");
    /* 
    assert(Org(e) != Dest(e), {inconsistent swap args});
    assert(Dest(a) != Dest(b), {inconsistent swap args});
    assert(Dest(a) != Org(e), {inconsistent swap args});
    assert(Dest(b) != Org(e), {inconsistent swap args});
    assert(Dest(a) != Dest(e), {inconsistent swap args});
    assert(Dest(b) != Dest(e), {inconsistent swap args});
    */
    /* Disconnect {e} from triangulation: */
    Splice(e, a);
    Splice(Sym(e), b);
    /* Connect it as the other diagonal: */
    a = Lnext(a);
    b = Lnext(b);
    Splice(e, a);
    Splice(Sym(e), b);
    /* Fix {Org} and {Dest}: */
    Odata(e) = Org(a);
    Odata(Sym(e)) = Org(b);
    /* Fix face pointers, if any: */
    Ldata(Lnext(a)) = Right(e);
    Ldata(Lnext(b)) = Left(e);
    /* Some consistency checks: */
    affirm(Lnext(Lnext(Lnext(e))) == e, "Left(e) is not a triangle");
    affirm(Lnext(Lnext(Lnext(Sym(e)))) == Sym(e), "Right(e) is not a triangle");
  }

Arc SPTriang_Connect(Arc a, Arc b)
  { Arc e = MakeEdge();
    Odata(e) = Dest(a);
    Ddata(e) = Org(b);
    Splice(e, Lnext(a));
    Splice(Sym(e), b);
    return e;
  }
  
void SPTriang_Disconnect(Arc a)
  { if (a != Onext(a)) { Splice(a, Oprev(a)); }
    a = Sym(a);
    if (a != Onext(a)) { Splice(a, Oprev(a)); }
  }

void SPTriang_RemoveSite(Arc b)
  { Site *pri = Org(b);
    Arc e;
    do 
      { e = Oprev(b); SPTriang_Disconnect(b); b = Lnext(e); }
    while (Dest(e) != pri);
  }
  
void SPTriang_InsertSiteInFace(Arc e, Site *x)
  { Arc base;
    Site *pri;

    base = MakeEdge();
    pri = Org(e);
    Odata(base) =  pri;
    Ddata(base) =  x; 
    Splice(base, e);
    do 
      { base = SPTriang_Connect(e, Sym(base));
        e = Oprev(base);
      }
    while (Dest(e) != pri);
  }

Arc_vec_t ArcsFromEdges(Arc_vec_t edge)
  { Arc_vec_t arc = SPQuad_Arc_vec_new(2*edge.ne);
    int i;
    for (i = 0; i < edge.ne; i++)
      { arc.e[2*i] = edge.e[i];
        arc.e[2*i+1] = Sym(edge.e[i]);
      }
    return arc;
  }
  
Arc_vec_t RenumberVertices(Arc_vec_t arc)
  { auto void VisitOrg(Arc a);
      /* Visits vertex {Org(a)}. */

    auto bool_t VisitedOrg(Arc a);
      /* TRUE if {Org(a)} has been visited already. */
   
    Arc_vec_t r = SPQuad_Arc_vec_new(1);
    int nVertices = 0;
    int i;
   
    bool_t VisitedOrg(Arc a)
      { Site *v = Org(a);
        return (v != NULL) && (v->num < nVertices) && (Org(r.e[v->num]) == v);
      }

    void VisitOrg(Arc a)
      { Site *v = Org(a);
        Arc b = a;
        v->num = nVertices;
        do 
          { affirm(Org(b) == v, "inconsistent Org");
            b = Onext(b);
          }
        while (b != a);
        SPQuad_Arc_vec_expand(&(r), nVertices);
        r.e[nVertices] = a;
        nVertices++;
      }

    for (i = 0; i < arc.ne; i++)
      { if (! VisitedOrg(arc.e[i])) { VisitOrg(arc.e[i]); } }

    SPQuad_Arc_vec_trim(&(r), nVertices);
    return r;
  }
  
Triangulation *SPTriang_RegularOctahedron(int smpOrder)
  { 
    auto S2Point oct_site(int i);
    S2Point oct_site(int i)
      { int k = i / 2;
        int s = i % 2;
        r3_t p = (r3_t){{ 0.0, 0.0, 0.0 }};
        p.c[k] = (s == 0 ? 1.0 : -1.0); 
        return p;
      }
      
    auto void set_Onext(Arc a, Arc b);
    void set_Onext(Arc a, Arc b)
      { if (Onext(a) != b) { Splice(a, Oprev(b)); } }
      
    SiteRef_vec_t site = SPTriang_MakeSites(6, oct_site);

    Arc x[2][2], y[2][2], z[2][2];
    int i, j;
    for (i = 0; i < 2; i++)
      { for (j = 0; j < 2; j++)
        { x[i][j] = MakeEdge();
          Odata(x[i][j]) = site.e[4+i];
          Ddata(x[i][j]) = site.e[0+j];
          y[i][j] = MakeEdge();
          Odata(y[i][j]) = site.e[0+i];
          Ddata(y[i][j]) = site.e[2+j];
          z[i][j] = MakeEdge();
          Odata(z[i][j]) = site.e[2+i];
          Ddata(z[i][j]) = site.e[4+j];
        }
      }

    set_Onext(x[0][0], Sym(z[0][0])); set_Onext(Sym(x[0][0]), y[0][1]);
    set_Onext(x[0][1], Sym(z[1][0])); set_Onext(Sym(x[0][1]), y[1][0]);
    set_Onext(x[1][0], Sym(z[1][1])); set_Onext(Sym(x[1][0]), y[0][0]);
    set_Onext(x[1][1], Sym(z[0][1])); set_Onext(Sym(x[1][1]), y[1][1]);
                                                                       
    set_Onext(y[0][0], Sym(x[0][0])); set_Onext(Sym(y[0][0]), z[0][1]);
    set_Onext(y[0][1], Sym(x[1][0])); set_Onext(Sym(y[0][1]), z[1][0]);
    set_Onext(y[1][0], Sym(x[1][1])); set_Onext(Sym(y[1][0]), z[0][0]);
    set_Onext(y[1][1], Sym(x[0][1])); set_Onext(Sym(y[1][1]), z[1][1]);
                                                                       
    set_Onext(z[0][0], Sym(y[0][0])); set_Onext(Sym(z[0][0]), x[0][1]);
    set_Onext(z[0][1], Sym(y[1][0])); set_Onext(Sym(z[0][1]), x[1][0]);
    set_Onext(z[1][0], Sym(y[1][1])); set_Onext(Sym(z[1][0]), x[0][0]);
    set_Onext(z[1][1], Sym(y[0][1])); set_Onext(Sym(z[1][1]), x[1][1]);
    
    /* Checking... */
    for (i = 0; i < 2; i++)
      { for (j = 0; j < 2; j++)
        { affirm(OrgDegree(x[i][j]) == 4, "wrong OrgDegree (x)");
          affirm(OrgDegree(y[i][j]) == 4, "wrong OrgDegree (y)");
          affirm(OrgDegree(z[i][j]) == 4, "wrong OrgDegree (z)");
          affirm(LeftDegree(x[i][j]) == 3, "wrong LeftDegree (x)");
          affirm(LeftDegree(y[i][j]) == 3, "wrong LeftDegree (y)");
          affirm(LeftDegree(z[i][j]) == 3, "wrong LeftDegree (z)");
        }
      }

    { Triangulation *tri = SPTriang_FromTopology(x[0][0], smpOrder);
      SPTriang_CheckTopology(tri, FALSE);
      return tri;
    }
  }

SiteRef_vec_t SPTriang_MakeSites
  ( int nSites, 
    S2Point vertex_pos(int i)
  )
  { SiteRef_vec_t site = SiteRef_vec_new(nSites);
    int i;
    for (i = 0; i < nSites; i++)
      { Site *s = (Site *)notnull(malloc(sizeof(Site)), "no mem for Site");
        s->num = i;
        s->pos = vertex_pos(i);
        site.e[i] = s;
      }
    return site;
  }

#define INIT_VISITED_SIZE 1024

Arc_vec_t SPTriang_CreateFaceRecords(Arc_vec_t arc)
  {
    auto bool_t VisitedLeft(SPQuad_Arc a);
      /* TRUE if face {Left(a)} has been visited before. */

    auto void VisitLeft(Arc a);
      /* If the left face of {a} has not been visited yet, 
        create its face record, set all pointers to it,
        and mark that face as visited. */

    Arc_vec_t side = SPQuad_Arc_vec_new(INIT_VISITED_SIZE);
    int nFaces = 0; 
      /* A face has been visited if all arcs around it have the {Left}
        field pointing to the same triangle record {t}, and
        {Left(side[f.num]) == t}. Outside of {VisitLeft}, either all
        arcs around the face satisfy this condition, or none of them
        does. */

    void VisitLeft(Arc a)
      { if (! VisitedLeft(a))
          { Face *t = (Face *)malloc(sizeof(Face));
            Arc b = a;
            int order = 0;
            t->num = nFaces;
            r3x3_zero(&(t->b2c));
            r3x3_zero(&(t->c2b));
            t->sp = (S2Point_vec_t){ 0, NULL };
            t->wp = (double_vec_t){ 0, NULL };
            do 
              { Ldata(b) =  t;
                b = Lnext(b);
                order++;
              }
            while (b != a);
            affirm(order == 3, "non-triangular face");
            SPQuad_Arc_vec_expand(&side, nFaces);
            side.e[nFaces] = a;
            nFaces++;
          }
      }

    bool_t VisitedLeft(Arc a)
      { Face *t = Left(a);
        return (t != NULL) && (t->num < nFaces) && 
          (Left(side.e[t->num]) == t);
      }

    int i;
    int nClosed = 0; 
      /* Faces {Left(side.e[0..nClosed-1])} are the visited faces whose 
        children were already visited.  The children of {Left(e)}, by 
        definition, are all faces that share a vertex or edge with 
        {Left(e)}. */

    /* Put the given faces on the visited (minus repetitions): */
    for (i = 0; i < arc.ne; i++) { VisitLeft(arc.e[i]); }
      
    /* Visit descendants of visited faces in BFS order: */
    while (nClosed < nFaces )
      { Arc s = side.e[nClosed];
        /* Enumerate edges around {Left(s)} including {s}: */
        Arc e = s;
        do
          { /* Enumerate edges around {Dest(e)} except {e}: */
            Arc a = Dnext(e);
            while (a != e) { VisitLeft(a); a = Dnext(a); }
            e = Lnext(e);
          }
        while (e != s);
        nClosed = nClosed + 1;
      }
    /* Return the triangles which were found: */
    SPQuad_Arc_vec_trim(&side, nFaces);
    return side;
  }
  
void SPTriang_ComputeGeometryDataOfFaces(Triangulation *tri, int smpOrder) 
  { Arc_vec_t side = tri->side;
    int i;

    for (i = 0; i < side.ne; i++) 
      { Arc ei = side.e[i];
        SPTriang_ComputeFaceMatrices(ei);
        SPTriang_ComputeFaceSamples(ei, smpOrder);
      }
    tri->smpOrder = smpOrder;
  }
  
void SPTriang_ComputeFaceMatrices(Arc e)
  { Face *f = Left(e);
    r3x3_t *b2c = &(f->b2c);
    r3x3_t *c2b = &(f->c2b);
    
    S2Point u = Dest(Onext(e))->pos;
    S2Point v = Org(e)->pos;
    S2Point w = Dest(e)->pos;

    *b2c = (r3x3_t){{
      {u.c[0], v.c[0], w.c[0]},
      {u.c[1], v.c[1], w.c[1]},
      {u.c[2], v.c[2], w.c[2]}
    }};
    r3x3_inv(b2c, c2b);
  }

void SPTriang_ComputeFaceSamples(Arc e, int smpOrder)
  { Face *f = Left(e);
    S2Point_vec_t *sp = &(f->sp);
    double_vec_t *wp = &(f->wp);
    
    int ns = 0;
    S2Point *u = &(Dest(Onext(e))->pos);
    S2Point *v = &(Org(e)->pos);
    S2Point *w = &(Dest(e)->pos);
    
    SPIntegral_SuperSampleTriangle(u, v, w, smpOrder, sp, wp, &ns);
    S2Point_vec_trim(sp, ns);
    double_vec_trim(wp, ns);
  }

Triangulation *SPTriang_FromTopology(Arc e, int smpOrder)
  { void *trit = malloc(sizeof(Triangulation));
    Triangulation *tri = (Triangulation *)notnull(trit, "no mem for triangulation");
    Arc_vec_t edge = SPQuad_RenumberEdges((Arc_vec_t){1, &e}, 0);
    tri->arc = ArcsFromEdges(edge);
    tri->out = RenumberVertices(tri->arc);
    tri->side = SPTriang_CreateFaceRecords((Arc_vec_t){1, &e});
    SPTriang_ComputeGeometryDataOfFaces(tri, smpOrder);
    return tri;
  }

Arc SPTriang_Locate(R3Point *p, Arc e)
  { R3Point q = FourCenter(e);
    while (1)
      { S2Point *a = &(Org(e)->pos);
        S2Point *b = &(Dest(e)->pos);
        S2Point *c = &(Dest(Onext(e))->pos);
        if (SamePoint(p, a))
          { return e; }
        else if (SamePoint(p, b))
          { return Sym(e); }
        else if (SamePoint(p, c))
          { return Sym(Onext(e)); }
        else if (PositiveSide(&(q), b, a, p))
          { e = Sym(e); }
        else if (! PositiveSide(&(q), c, a, p))
          { e = Onext(e); }
        else if (! PositiveSide(&(q), b, c, p))
          { e = Dprev(e); }
        else 
          { return e; };
      }
  }
  
R3Point FourCenter(Arc e)
  { R3Point c = Org(e)->pos;
    r3_add(&(c), &(Dest(e)->pos), &(c));
    r3_add(&(c), &(Dest(Onext(e))->pos), &(c));
    r3_add(&(c), &(Dest(Oprev(e))->pos), &(c));
    r3_scale(0.25, &(c), &(c));
    return c;
  }

bool_t PositiveSide(R3Point *a, R3Point *b, R3Point *c, R3Point *d)
  { R3Point ba, ca, da;
    double det;
    r3_sub(b, a, &(ba));
    r3_sub(c, a, &(ca));
    r3_sub(d, a, &(da));
    det = r3_det(&(ba), &(ca), &(da));
    return det > 0.0;
  }

double TriDist(S2Point *p, r3x3_t *c2b)
  { int i;
    double d = 0.0;
    r3_t b;
    r3x3_map_col(c2b, p, &b);
    for (i = 0; i < 3; i++)
      { double bi = b.c[i];
        if ((bi < 0.0) && (-bi > d)){ d = -bi; }
      }
    return d;
  }

double SphericalTriangleArea(S2Point *u, S2Point *v, S2Point *w)
  { r3_t uv, vw, wu;
    r3_cross(u, v, &uv);
    if (r3_norm_sqr(&uv) == 0) { return 0.0; } 
    r3_cross(v, w, &vw);
    if (r3_norm_sqr(&vw) == 0) { return 0.0; } 
    r3_cross(w, u, &wu);
    if (r3_norm_sqr(&wu) == 0) { return 0.0; } 
    { double au = acos(- r3_cos(&wu, &uv)); 
      double av = acos(- r3_cos(&uv, &vw)); 
      double aw = acos(- r3_cos(&vw, &wu));
      return au + av + aw - PI;
    }
  }   

double TriangleArea(R3Point *u, R3Point *v, R3Point *w)
  { r3_t uv, uw, n;
    r3_sub(v, u, &uv);
    r3_sub(w, u, &uw);
    r3_cross(&uv, &uw, &n);
    return 0.5 * r3_norm(&n);
  } 
  
void Project2(R3Point *u, R3Point *v, R3Point *p, double *au, double *av)
  { double duu = r3_dot(u, u);
    double duv = r3_dot(u, v);
    double dvv = r3_dot(v, v);
    double dup = r3_dot(u, p);
    double dvp = r3_dot(v, p);
    double det = duu*dvv - duv*duv;
    (*au) = (dup*dvv - dvp*duv)/det;
    (*av) = (duu*dvp - duv*dup)/det;
  }
  
void SPTriang_ConvexCone(int *N, r3_t *v)
  { int NV = *N; /* Number of input vectors. */
    int NH = 0;  /* Number of extremal vectors on convex cone. */
    int i;
    for (i = 0; i < NV; i++)
      { /* Add {v[i]} to the cone: */
        r3_t vi = v[i];  /* Copy it, we may need the space: */
        if (r3_L_inf_norm(&vi) == 0.0)
          { /* Null vector, ignore. */ }
        else if (NH == 0)
          { v[NH] = vi; NH++; }
        else if (NH == 1)
          { double c = r3_cos(&vi, &(v[0]));
            if (c == 1.0)  
              { /* Ignore, {vi} is collinear with {v[0]}. */ }
            else if (c == -1.0)
              { /* {vi} is antipodal to {v[0]}: */ (*N) = -1; return; }
            else
              { v[NH] = vi; NH++; }
          }
        else if (NH == 2)
          { /* Check {vi} against the 2-dim cone {v[0],v[1]}: */
            double d = r3_det(&vi, &(v[0]), &(v[1]));
            if (d > 0.0)
              { v[NH] = vi; NH++; }
            else if (d < 0.0)  
              { v[2] = v[1]; v[1] = vi; NH++; }
            else 
              { double a0, a1;
                Project2(&(v[0]), &(v[1]), &vi, &a0, &a1);
                if ((a0 <= 0.0) && (a1 <= 0.0)) 
                  { /* {vi} is antipode to the cone: */ (*N) = -1; return; }
                else if (a0 < 0.0) 
                  { v[1] = vi; }
                else if (a1 < 0.0) 
                  { v[0] = vi; }
                else
                  { /* {vi} lies in the cone of {v[0],v[1]}, ignore. */ }
              }
          }
        else
          { /* Look for first & last cone face {v[j],v[j+1]} that excludes {vi}: */
            int r, last = -1, frst = -1, nneg = 0; 
            double dprev = -1.0, dini;
            for (r = 0; r <= NH; r++)
              { int j = (r % NH);
                int k = (r + 1) % NH;
                double dijk = (r == NH ? dini : r3_det(&vi, &(v[j]), &(v[k])));
                if ((dprev > 0) && (dijk <= 0.0)) { last = j; }
                if ((dprev <= 0) && (dijk > 0.0)) { frst = k; }
                if ((last != -1) && (frst != -1)) { break; }
                if (r == 0) { dini = dijk; }
                if ((r < NH) && (dijk < 0)) { nneg++; }
                dprev = dijk;
              }
            if (nneg == NH) { /* {vi} antipodal to cone. */ (*N) = -1; return; }
            if (nneg == 0) { /* {vi} already is inside cone: */ break; }
            /* Ensure one free slot for {vi} between {v[last]} and {v[frst]}: */
            { int gap = ((NH + frst) - (last+1)) % NH;
              if (gap == 0)
                { if (frst > last)
                    { for (r = NH; r > frst; r--) { v[r] = v[r-1]; }
                      frst++;
                    }
                }
              else if (gap >= 2)
                { if (frst > last)
                    { for (r = frst; r < NH ; r++) { v[r-gap+1] = v[r]; }
                      frst = frst-gap+1; 
                    }
                  else
                    { for (r = frst; r <= last ; r++) { v[r-frst] = v[r]; }
                      frst = 0; last = last - frst;
                    }
                }
              NH = NH - gap + 1; 
            }
            /* Now insert {vi} in its place: */
            v[last+1] = vi;
          }
      }
    (*N) = NH;
  }

bool_t TrianglesOverlap
  ( S2Point *u, S2Point *v, S2Point *w, 
    S2Point *p, S2Point *q, S2Point *r
  )
  { /* Two (convex) spherical triangles are disjoint if and only if
      the corresponding trihedra are disjoint. Two convex trihedra T1,
      T2 are disjoint iff there is a plane that separates them. A
      plane separates a set T1 from a set T2 iff it bounds T1 and the
      antipode of T2. There is a plane bounding two sets iff their
      convex enclosing cone is defined. A convex cone encoses a convex
      trihedron iff it encloses its defining vectors. Thus, we try to
      build the convex enclosing cone of {u,v,w,-p,-q,-r} --- if that
      succeeds, then the triangles don't overlap. */
      
    r3_t t[6];
    int N = 6;
    t[0] = *u; t[1] = *v; t[2] = *w;
    r3_neg(p, &(t[3])); r3_neg(q, &(t[4])); r3_neg(r, &(t[5])); 
    SPTriang_ConvexCone(&N, &(t[0]));
    return (N == 0); 
  }

bool_t SPTriang_EdgeIsDegnerate(Arc e)
{ if ((Onext(e) == e) || (Dnext(e) == e))
      { /* {Left(e)} and {Right(e)} are the same: */ return TRUE; }
    else
      { Site *p = Org(e); 
        Site *q = Dest(e);
        Site *u = Dest(Onext(e));
        Site *v = Dest(Oprev(e));
        affirm(Lnext(e) != e, "Left(e) is monogon");
        affirm(Lnext(Sym(e)) != Sym(e), "Right(e) is monogon");
        affirm(Lnext(Lnext(Lnext(e))) == e, "Left(e) not triangle");
        affirm(Lnext(Lnext(Lnext(Sym(e)))) == Sym(e), "Right(e) not triangle");
        return 
          (u == v) || (u == p) || (u == q) || 
          (v == q) || (v == p) || (p == q);
      }
  }

Sign SPTriang_FaceOrientation(Arc e)
  { S2Point *p = &(Org(e)->pos); 
    S2Point *q = &(Dest(e)->pos);
    S2Point *u = &(Dest(Onext(e))->pos);
    double du = r3_det(p, q, u);
    return (du > 0.0 ? +1 : (du < 0.0 ? -1 : 0));
  }

bool_t SPTriang_EdgeIsOutwards(Arc e)
  { return SPTriang_FaceOrientation(e) > 0 && SPTriang_FaceOrientation(Sym(e)) > 0; }

bool_t SPTriang_EdgeIsSilhouette(Arc e)
  { return SPTriang_FaceOrientation(e)*SPTriang_FaceOrientation(Sym(e)) <= 0; }

void SPTriang_CheckTopology(Triangulation *tri, bool_t verbose)
  { 
    int an, vn, fn;
    int NA = tri->arc.ne;
    int NV = tri->out.ne;
    int NF = tri->side.ne;
    
    /* Check if all arcs are correctly stored: */
    for (an = 0; an < NA; an += 2)
      { Arc e = tri->arc.e[an];
        Arc f = tri->arc.e[an + 1];
        affirm(Edge(e)->num == an/2, "inconsistent edge nums");
        affirm(Sym(e) == f, "inconsistent arc numbers");
      }

    /* Check if all {Org} pointers are correctly set: */
    if (verbose) { fprintf(stderr, "vertex degrees: "); }
    for (vn = 0; vn < NV; vn++)
      { Arc e = tri->out.e[vn];
        Arc a = e;
        int order = 0;
        do 
          { affirm(Org(a)->num == vn, "inconsistent vertex nums");
            a = Onext(a);
            order++;
          } 
        while (a != e);
        if (verbose) { fprintf(stderr, " %d", order); }
      }
    if (verbose) { fprintf(stderr, "\n"); }

    /* Check if all {Left} pointers are correctly set: */
    for (fn = 0; fn < NF; fn++)
      { Arc e = tri->side.e[fn];
        Arc a = e;
        int order = 0;
        do 
          { affirm(Left(a)->num == fn, "inconsistent face nums");
            order++;
            a = Lnext(a);
          } 
        while (a != e);
        affirm(order == 3, "non-triangular faces");
      }
  }

bool_t Adjacent(Arc e, Arc f)
  { Site *u = Org(e);
    Arc b = f;
    do 
      { if (Dest(b) == u) { return TRUE; }
        b = Onext(b);
      } 
    while (b != f);
    return FALSE;
  }

/* TRIANGULATION I/O */

void SPTriang_Print(Triangulation *t, char *name)
  { auto void PrintTriangEdge(Arc e);
    auto void PrintArc(Arc e);

    FILE *f = open_write(txtcat(name, ".txt"), TRUE);
    int i;

    void PrintArc(Arc e)
      { fprintf(f, "%d:%d", Edge(e)->num, QuadBits(e)); }

    void PrintTriangEdge(Arc e)
      { Site *op = Org(e); 
        Site *dp = Dest(e); 
        PrintArc(e); fputc(' ', f);
        PrintArc(Onext(e)); fputc(' ', f);
        PrintArc(Oprev(e)); fputc(' ', f); fputc(' ', f);
        fprintf(f, "%4d - %4d ", op->num, dp->num);
        fprintf(f, "   ");
        { Face *tr = (Face *)Ldata(e); 
          if (tr != NULL)
            { fprintf(f, "%4d", tr->num); }
          else
            { fprintf(f, "NULL"); }
        }
        fprintf(f, " - ");
        { Face *tr = (Face *)Rdata(e);
          if (tr != NULL)
            { fprintf(f, "%4d", tr->num); }
          else
            { fprintf(f, "NULL"); }
          fprintf(f, "\n");
        }
      }

    fprintf(f, "# %d Vertices:\n", t->out.ne);
    for (i = 0; i < t->out.ne; i++)
      { if (t->out.e[i] != NullArc)
          { Site *v = Org(t->out.e[i]); 
            fprintf(f, "%4d  ", v->num); 
            r3_gen_print(f, &(v->pos), "12.9f", "( ", ", ", " )");
            fprintf(f, "\n");
          }
      }

    fprintf(f, "# %d Edges:\n", t->arc.ne / 2);
    for (i = 0; i < t->arc.ne; i += 2)
      { if (t->arc.e[i] != NullArc)
        { PrintTriangEdge(t->arc.e[i]); }
      }

    /* Should print triangles, too */
    fclose(f);
  }

void SPTriang_Write(FILE *wr, Triangulation *t)
  { 
    auto void WritePoint(S2Point *p);
    auto void WriteArc(Arc a);
    auto void WriteOrg(Arc a);
    auto void WriteLeft(Arc a);
    auto void WriteEdgeData(Arc e);

    int NV = t->out.ne;
    int NE = t->arc.ne / 2;
    int NF = t->side.ne;

    void WritePoint(S2Point *p)
      { int i;
        for (i = 0; i <= 2; i++)
          { if (i != 0){ fputc(' ', wr); }
            fprintf(wr, "%.16g", p->c[i]);
          }
      }

    void WriteArc(Arc e)
      { fprintf(wr, "%d:%d", Edge(e)->num, QuadBits(e)); }

    void WriteOrg(Arc a)
      { fputc('v', wr); fprintf(wr, "%d", Org(a)->num); }

    void WriteLeft(Arc a)
      { fputc('f', wr); fprintf(wr, "%d", Left(a)->num); }

    void WriteEdgeData(Arc e)
      { WriteArc(e); fputc(' ', wr);
        fputc(' ', wr);

        WriteArc(Onext(e)); fputc(' ', wr);
        WriteArc(Onext(Sym(e))); fputc(' ', wr);
        fputc(' ', wr);

        WriteOrg(e); fputc(' ', wr);
        WriteOrg(Sym(e)); fputc(' ', wr);
        WriteLeft(e); fputc(' ', wr);
        WriteLeft(Sym(e)); 
        fputc('\n', wr);
      }

    int i;

    filefmt_write_header(wr, "triangulation", Triangulation_FileFormat);

    fprintf(wr, "vertices = %d\n", NV);
    fprintf(wr, "edges = %d\n", NE);
    fprintf(wr, "faces = %d\n", NF);

    fprintf(wr, "vertices:\n");
    for (i = 0; i < NV; i++)
      { Arc a = t->out.e[i];
        Site *v; 
        affirm(t->out.e[i] != NullArc, "vertex with no arcs");
        v = Org(a);
        fprintf(wr, "%d ", i);
        WriteArc(a); fputc(' ', wr);
        WritePoint(&(v->pos)); fputc('\n', wr);
      }

    fprintf(wr, "edges:\n");
    for (i = 0; i < t->arc.ne; i += 2)
      { affirm(t->arc.e[i] != NullArc, "edge with no arcs");
        fprintf(wr, "%d ", i / 2);
        WriteEdgeData(t->arc.e[i]);
      }

    fprintf(wr, "faces:\n");
    for (i = 0; i < t->side.ne; i++)
      { affirm(t->side.e[i] != NullArc, "face with no arcs");
        fprintf(wr, "%d ", i);
        WriteArc(t->side.e[i]);
        fputc('\n', wr);
      }

    filefmt_write_footer(wr, "triangulation");
    fflush(wr);
  }

Triangulation *SPTriang_Read(FILE *rd, int smpOrder)
  { int NV, NE, NF;
    Triangulation *t = (Triangulation*)notnull(malloc(sizeof(Triangulation)), "out of mem");
    filefmt_read_header(rd, "triangulation", Triangulation_FileFormat);
    NV = nget_int32(rd, "vertices"); fget_eol(rd);
    affirm(NV >= 0, "bad NV");
    NE = nget_int32(rd, "edges"); fget_eol(rd);
    affirm(NE >= 0, "bad NE");
    NF = nget_int32(rd, "faces"); fget_eol(rd);
    affirm(NF >= 0, "bad NF");
    t->out = SPQuad_Arc_vec_new(NV);
    t->arc = SPQuad_Arc_vec_new(2*NE);
    t->side = SPQuad_Arc_vec_new(NF);
    {
      auto Arc GetArc(void);
      auto Site *GetVertex(void);
      auto Face *GetFace(void);
      auto void ReadPoint(r3_t *c);
     
      SiteRef_vec_t site = SiteRef_vec_new(NV);
      Arc_vec_t edge = SPQuad_Arc_vec_new(NE);
      FaceRef_vec_t face = FaceRef_vec_new(NF);
     
      Arc GetArc(void)
        { int qnum, rnum; 
          Arc a;
          fget_skip_spaces(rd);
          qnum = fget_int32(rd);
          affirm((qnum >= 0) && (qnum < NE), "bad quad number");
          fget_match(rd, ":");
          rnum = fget_int32(rd);
          affirm((rnum >= 0) & (rnum <= 3), "inconsistent rot number");
          a = edge.e[qnum];
          while (QuadBits(a) != rnum) { a = Rot(a); }
          return a;
        }

      Site *GetVertex(void)
        { int vnum;
          fget_skip_spaces(rd);
          fget_match(rd, "v");
          vnum = fget_int32(rd);
          affirm((vnum >= 0) && (vnum < NV), "bad vertex number");
          return site.e[vnum];
        }

      Face *GetFace(void)
        { int fnum;
          fget_skip_spaces(rd);
          fget_match(rd, "f");
          fnum = fget_int32(rd);
          affirm((fnum >= 0) && (fnum < NF), "bad vertex number");
          return face.e[fnum];
        }

      void ReadPoint(r3_t *c)
        { int i;
          for (i = 0; i <= 2; i++)
            { fget_skip_spaces(rd); 
              c->c[i] = fget_double(rd);
            }
        }

      /* Create all vertices, edges, and faces: */
      int i;
      for (i = 0; i < NV; i++)
        { Site *v = (Site*)notnull(malloc(sizeof(Site)), "out of mem");
          site.e[i] = v; 
          v->num = i;
        }
      for (i = 0; i < NE; i++)
        { Arc a = MakeEdge();
          Edge(a)->num = i;
          edge.e[i] = a;
        }
      for (i = 0; i < NF; i++)
        { Face *f = (Face *)notnull(malloc(sizeof(Face)), "out of mem");
          face.e[i] = f; 
          f->num = i; 
        }

      /* Read vertex coordinates: */
      fget_match(rd, "vertices:"); fget_eol(rd);
      for (i = 0; i < NV; i++)
        { Site *v = site.e[i];
          int j = fget_int32(rd);
          affirm(j == i, "inconsistent vertex number");
          fget_skip_spaces(rd);
          t->out.e[i] = GetArc();
          fget_skip_spaces(rd);
          ReadPoint(&(v->pos)); fget_eol(rd);
        }

      /* Read edges pointers: */
      fget_match(rd, "edges:"); fget_eol(rd);
      for (i = 0; i < NE; i++)
        { int j = fget_int32(rd);
          affirm(j == i, "inconsistent edge number");
          fget_skip_spaces(rd);
          { Arc a = GetArc();
            Arc b = GetArc(); 
            Arc c = GetArc();
            Site *o = GetVertex();
            Site *d = GetVertex();
            Face *l = GetFace();
            Face *r = GetFace();
            t->arc.e[2*i] = a; 
            t->arc.e[2*i+1] = Sym(a);
            Splice(a, Oprev(b));
            if (b != Sym(a)){ Splice(Sym(a), Oprev(c)); }
            Odata(a) =  o;
            Ddata(a) =  d;
            Ldata(a) =  l;
            Rdata(a) =  r;
          }
          fget_eol(rd);
        }

      /* Read face data: */
      fget_match(rd, "faces:"); fget_eol(rd);
      for (i = 0; i < NF; i++)
        { int j = fget_int32(rd);
          affirm(j == i, "inconsistent face number");
          fget_skip_spaces(rd);
          t->side.e[i] = GetArc();
          fget_eol(rd);
        }
    }
    filefmt_read_footer(rd, "triangulation");
    SPTriang_ComputeGeometryDataOfFaces(t, smpOrder);
    return t;
  }

static char **nameCache = NULL;
static Triangulation **triCache = NULL;
static int szCache = 0; /* Allocated size of {nameCache}, {triCache}. */
static int nCached = 0; /* Occupied entries in {nameCache}, {triCache}. */

Triangulation *SPTriang_ReadCached(char *fileName, int smpOrder)
  { int i;
    Triangulation *tri;
    /* Look up in cache */
    for (i = 0; i < nCached; i++)
      { if (strcmp(fileName, nameCache[i]) == 0) 
          { tri = triCache[i];
            affirm(tri->smpOrder == smpOrder, "inconsistent smpOrders");
            return tri;
          }
      }
    /* Not found, read from file: */
    tri = SPTriang_Read(open_read(fileName, TRUE), smpOrder);
    /* Ensure there is space in cache: */
    if (nCached >= szCache)
      { /* Expand cache: */
        int sztmp = 2*szCache + 10;
        char **tmpc = (char **)notnull(malloc(sztmp*sizeof(char*)), "no mem for tmpc");         Triangulation **tmpt = (Triangulation **)notnull(malloc(sztmp*sizeof(Triangulation*)), "no mem for tmpt");
        for (i = 0; i < nCached; i++) 
          { tmpc[i] = nameCache[i]; tmpt[i] = triCache[i]; }
        if (nameCache != NULL) { free(nameCache); }
        if (triCache != NULL) { free(triCache); }
        nameCache = tmpc;
        triCache = tmpt;
        szCache = sztmp;
      }
    /* Store in cache: */
    nameCache[nCached] = fileName;
    triCache[nCached] = tri;
    nCached++;
    return tri;
  }

/* Arrays of Site*: */

vec_typeimpl(SiteRef_vec_t,SiteRef_vec,Site *);

/* Arrays of Face*: */

vec_typeimpl(FaceRef_vec_t,FaceRef_vec,Face *);