(* Delaunay triangulation on a sphere. *)

INTERFACE SPTriang;

(* See 

   "Primitives for the Manipulation of General Subdivisions 
   and the Computation of Voronoi Diagrams"

   L. Guibas, J. Stolfi, ACM TOG, April 1985

*)

IMPORT SPQuad, HLR3, LR3x3, Wr, Rd;

TYPE 
  Arc = SPQuad.Arc;
  Coords = ARRAY [0..2] OF LONGREAL;
  SiteNumber = CARDINAL;
  Site = RECORD num: SiteNumber; c: Coords END;
  
  TriangleNumber = CARDINAL;
  Triangle = RECORD num: TriangleNumber END;
  
TYPE
  T = RECORD
      out: REF ARRAY OF Arc;
        (*
          Element "out[i]" is one arc out of site number "i". *)
      arc: REF ARRAY OF Arc;
        (*
          Elements "arc[2*i]" and "arc[2*i+1]" are
          the same edge in both orientations. *)
      side: REF ARRAY OF Arc;
        (*
          Element "side[i]" is one arc whose left face is triangle "i". *)
      b2c: REF ARRAY OF LR3x3.T;
      c2b: REF ARRAY OF LR3x3.T;
        (*
          "MapCol(b2c[i],a)" computes the Cartesian coordinates
          of a point given its barycentric coordinates "a" 
          relative to the corners of triangle "i".
          The corners are, in order, "Dest(Onext(e))",
          "Org(e)", "Dest(e)", where "e = side[i]".
          These fields are set by "ComputeTriangleMatrices". *)
    END;
  
PROCEDURE BuildTetrahedron (READONLY site: ARRAY [0..3] OF REF Site): Arc; 
  (*
    Builds a tetrahedron with vertices "site[0],..site[3]".
    Returns one arc of the tetrahedron. *)

PROCEDURE InsertSite (READONLY s: REF Site; e: Arc): Arc;
  (*
    Adds a new site "s" to the convex hull "e".
    Returns an arc of the new hull. *)

PROCEDURE BuildInc(READONLY site: ARRAY OF REF Site): Arc;
  (*
    Builds the convex hull of the given sites incrementally, by
    creating a tetrahedron with the first four sites and inserting the
    others one by one.  Returns an arc of the hull. *)

PROCEDURE Locate (READONLY p: Coords; e: Arc): Arc;
  (* 
    Returns an edge "f" of the polyhedron "e" such that the 
    face "Left(f)" is visible from the point "p" (which is 
    assumed to be outside the polyhedron). *)

PROCEDURE SamePoint(READONLY p, q: Coords): BOOLEAN;
  (*
    True iff "p = q". *)

PROCEDURE RemoveSite (b: Arc);
  (* 
    Removes Dest(b) and all its incident edges. *)

PROCEDURE Retriangulate(
    e: Arc; 
    deg: INTEGER;
    VAR tr: ARRAY OF Arc;
    VAR nTr: CARDINAL;
  );
  (* 
    Retriangulates the face Left(e), resulting from the elimination of
    a vertex, which must have "deg" sides.
    Also returns in "tr[0..nTr]" one side arc out of each new
    triangle. *)

PROCEDURE Check(READONLY arc: ARRAY OF Arc);
  (*
    Checks consistency of Delaunay triangulation; aborts if 
    inconsistent. *)
  
PROCEDURE PositiveSide(READONLY a, b, c, d: Coords): BOOLEAN;
  (* 
    TRUE if "d" is on the positive side of plane Join(a,b,c) *)

PROCEDURE Degree(e: Arc): CARDINAL;
  (*
    Number of arcs out of Org(e). *)

(* Output procedures: *)

PROCEDURE NumberVertices(READONLY aa: ARRAY OF Arc): REF ARRAY OF Arc;
  (*
    Renumbers the origin vertices of the arcs in "aa[..]" 
    consecutively from 0.  Returns a vector where element "i"
    is one arc out of vertex number "i". *)

PROCEDURE CollectVertices(e: Arc; VAR out: ARRAY OF Arc);
  (*
    Enumerates all vertices of the hull "e", and returns in 
    "out[i]" one arc out site number "i" (or NullArc,
    if site "i" does not appear in the hull). *)
    
PROCEDURE NumberTriangles(READONLY arc: ARRAY OF Arc): REF ARRAY OF Arc;
  (*
    Creates one Triangle record for each face of the
    triangulation reachable from "arc[...]", and sets the 
    Ldata fields of all edges.  The triangles
    will be numbered consecutively starting from 0.
    Returns a vector with one Arc out of each triangle.  *)

PROCEDURE CollectTriangles(e: Arc; VAR side: ARRAY OF Arc);
  (*
    Stores in "side[i]" one Arc whose Ldata is triangle 
    number "i", if triangle "i" is reachable;
    or "NullArc" if triangle "i" is unreachable.
    Assumes NumberTriangles has been called. *)

PROCEDURE ComputeTriangleMatrices(VAR tri: T);
  (*
    Computes the matrices "b2c[i]" and "c2b[i]" for all triangles
    in "tri". *)
  
PROCEDURE BuildTables(e: Arc): T;
  (*
    Calls NumberArcs(e), NumberVertices(e), NumberTriangles(e),
    and ComputeTriangleMatrices(e), and packs the resulting vectors
    into a "T" record. *)

PROCEDURE PerspDraw(
    READONLY arc: ARRAY OF Arc; (* Result of NumberArcs *)
    name: TEXT; 
    obs, upp: HLR3.Point;
  );
  (*
    Writes a file "name.ps" containing a perspective view of 
    the triangulation, as seen from "obs" and looking towards 
    the center of the sphere.  The point "upp"
    defines the vertical plot axis. *)
  
PROCEDURE OutputX3D(
    READONLY side: ARRAY OF Arc; (* Result of NumberTriangles *)
    name: TEXT;
  );
  (*
    Writes a file "name.poly" containing a 3D representation of 
    the triangulation, suitable for viewing with the "convert-poly"
    and "x3d" tools. *)
    
PROCEDURE Write(wr: Wr.T; READONLY t: T);
  (*
    Writes a description of "t" into "wr", in a format that 
    can be read back with "Read" below. *)
    
PROCEDURE Read(rd: Rd.T): T;
  (*
    Reads a triangulation "t" from the reader "rd".
    Assumes the format used by "Write" above. *)

PROCEDURE Print(READONLY t: T; name: TEXT);
  (*
    Writes a file "name.txt" with a readable description of the
    triangulation. *)
  
PROCEDURE PrintCoords(
    f: Wr.T; 
    READONLY p: Coords; 
    lp: TEXT := "( ";
    rp: TEXT := " )";
    cm: TEXT := ", ";
  );
  (*
    Prints the given site coordinates to "f". *)
  
(* Quad-edge data pointers: *)

PROCEDURE Org (e: Arc): REF Site;
PROCEDURE Dest (e: Arc): REF Site;

PROCEDURE Left(e: Arc): REF Triangle;
PROCEDURE Right(e: Arc): REF Triangle;

END SPTriang.