home *** CD-ROM | disk | FTP | other *** search
/ Resource Library: Graphics / graphics-16000.iso / general / raytrace / renaisnc / delaunay.lha / Delaunay / DCEL.h < prev    next >
C/C++ Source or Header  |  1992-01-04  |  12KB  |  406 lines
  1. /*
  2.  File: DCEL.h
  3.  Author: K.R. Sloan
  4.  Last Modified: 30 August 1988
  5.  Purpose: Doubly-Connected Edge List (improved),
  6.           along with a few utilities for triangulating planar
  7.           point sets
  8.  
  9.  */
  10. #define ArraySize 200000
  11. #include <stdio.h>
  12. #define true (-1)
  13. #define false (0)
  14.  
  15.  /*
  16.     These forward declarations are needed to make the compiler
  17.     happy.
  18.  */
  19. typedef struct dcelPoint Point;
  20. typedef struct dcelPSE   PointSetElement;
  21. typedef PointSetElement  *PointSetPointer;
  22. typedef struct dcelEE    EdgeEnd;
  23. typedef struct dcelE     Edge;
  24. typedef Edge             *EdgePointer;
  25. typedef struct dcelT     Triangle;
  26. typedef Triangle         *TrianglePointer;
  27.  
  28.   /*
  29.    Points are stored in an array, so that an application can pun
  30.    on the array index and keep auxiliary information.  We only care about
  31.    the x,y coordinates, and the Edges emanating from this vertex
  32.   */ 
  33. struct dcelPoint
  34.         {
  35.          REAL x, y;
  36.          Edge *Edges;
  37.         };  
  38.  
  39.   /*
  40.    On the other hand, we occasionally need a "set of points".  For that,
  41.    we use PointSetElements, arranged in a circular, doubly-linked list.
  42.   */
  43. struct dcelPSE
  44.         {
  45.          int P;
  46.          PointSetElement *Previous, *Next;
  47.         };
  48.  
  49.   /*
  50.    Edges have two ends.  Each end has a vertex, links for a circular,
  51.    doubly-linked list, and an associated Triangle.
  52.    In addition, we will need sets of Edges - yet another circular, 
  53.    doubly-linked list.
  54.   */
  55. struct dcelEE
  56.         {
  57.          int V;       /* Vertex at this End */
  58.          Triangle *T; /* Triangle on the LEFT, seen from this End */ 
  59.          Edge *Previous, *Next;
  60.         };
  61. struct dcelE
  62.         {
  63.          EdgeEnd Head;
  64.          EdgeEnd Tail;
  65.          Edge *Previous, *Next;
  66.         };
  67.  
  68.   /*
  69.    Triangles have three vertices, three edges, and a figure-of-merit.
  70.    Like everything else, they can be collected on circular, doubly-connected
  71.    linked lists.
  72.  
  73.    For the truly ambitious: generalize the Triangle to Facet.  Give up
  74.    the efficiency of immediate access to the vertices and edges, in favor
  75.    of arbitrary k-gons (yup, more circular, doubly connected linked lists.
  76.  
  77.    Me, I like triangles... 
  78.   */
  79. struct dcelT
  80.         {
  81.          REAL Merit;             /* See Rho, and Shrink */
  82.          int P1, P2, P3;           /* turning to the RIGHT */
  83.          Edge *E1, *E2, *E3;       /* Ei = <Pi,P(i+1)mod3> */
  84.          Triangle *Previous, *Next;
  85.         }; 
  86.  
  87.   /*
  88.    Thus, by convention:
  89.  
  90.          LOOKING FROM THE TAIL               LOOKING FROM THE HEAD
  91.        ------------------------              ---------------------
  92.                Head.V                               Tail.V
  93.                  ^                                    |
  94.                  |                                    |
  95.        Tail.T    |   Head.T                  Head.T   |   Tail.T
  96.                  |                                    |
  97.                  |                                    \/ 
  98.                Tail.V                               Head.V  
  99.   */
  100.  
  101.   /*
  102.     Allocation routines
  103.    */
  104. /*
  105. extern PointSetPointer NewPoint();
  106. extern Edge *NewEdge();
  107. extern Triangle *NewTriangle();
  108. */
  109.  
  110.   /* for Animation */
  111.  extern int ANIMATING;
  112.  
  113.   /* for error messages */
  114.  extern int DCELVERBOSE;
  115.   /* data */
  116.  extern Point Points[ArraySize];
  117.   /*
  118.    Triangulation information -  a Hull, a set of Triangles, a set of
  119.    "Suspect" edges, and a "ProbeRoot" (a vertex used as a hint to the
  120.    Search procedure.
  121.  
  122.    Only one region here.  You want more?  Looks easy, go do it.
  123.  
  124.   */
  125.  extern PointSetPointer Hull;
  126.  extern TrianglePointer Triangles;
  127.  extern EdgePointer Suspects;
  128.  extern int ProbeRoot;
  129.  extern REAL AnimateColor, AnimateErase;
  130.  
  131.  extern void Animate();
  132.  
  133.  /*
  134.   Add and Delete Suspect handle a queue of edges.
  135.  
  136.   AddSuspect allows the same edge to be added multiple times. 
  137.   The first addition puts the edge at the end of the queue.
  138.   Subsequent additions do nothing.
  139.  
  140.   RemoveSuspect may be called even when the edge is not in the queue.
  141.  */ 
  142.  extern void AddSuspect();
  143.  extern void RemoveSuspect();
  144.  
  145.  /*
  146.   FindEdge tries to find an existing edge (P1,P2) in the edge list for P1.
  147.  */
  148.  extern Edge *FindEdge ();
  149.  
  150.  /*
  151.   Does the path A -> B -> C turn to the left?
  152.  */
  153.  extern int LeftTurn();
  154.  
  155.  /*
  156.   This little gem is due to Rob Fowler.  Given two points A, and B
  157.   we calculate a number in [0.0, 8.0) which is a monotonic function of
  158.   the direction from A to B. 
  159.  
  160.   (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0) correspond to
  161.   (  0,  45,  90, 135, 180, 225, 270, 315, 360) degrees, measured
  162.   counter-clockwise from the positive x axis.
  163.  */
  164.  extern REAL FowlerAngle();
  165.  
  166.  /*
  167.   FindNeighbors locates the edges which would be the left
  168.   and right neighbors of the edge O-P, in the list of edges
  169.   emanating from O.
  170.  
  171.   ASSUMPTION: there must be at least 1 Edge emanating from O; check
  172.               before you call...
  173.  */
  174.  extern void FindNeighbors ();
  175.  
  176.  /*
  177.   AddEdge inserts the edge P1-P2, and the information that this edge
  178.   has the Triangle T on its left.
  179.  */
  180.  extern Edge *AddEdge ();
  181.  
  182.  /*
  183.   RemoveEdge removes the Triangle T as support for Edge E.
  184.   If there is no triangle on the other side of E, then E is destroyed.
  185.  */
  186.  extern void RemoveEdge();
  187.  
  188.  /*
  189.   AddTriangle creates the triangle A,B,C - with figure of merit R.
  190.   (See Rho and Shrink for the definition and use of R)
  191.  */
  192.  extern void AddTriangle ();
  193.  
  194.  /*
  195.   RemoveTriangle destroys T
  196.  */
  197.  extern void RemoveTriangle();
  198.  
  199.  /*
  200.    InsideTriangle decides if a point P is Inside, or OnEdge
  201.    of the triangle defined by A, B, C.
  202.  
  203.    If P is on an edge, V1 and V2 identify that edge
  204.  */
  205.  extern void InsideTriangle ();
  206.  
  207.  
  208.  /*
  209.   Given (the GLOBAL variable) ProbeRoot - a vertex in a 
  210.   CONVEX, TRIANGULATED region
  211.   of the plane, Search attempts to find the Triangle containing P.
  212.   This is accomplished by traversing the triangulation from ProbeRoot
  213.   towards P until an enclosing triangle is found, or P is shown to be
  214.   outside the convex hull. 
  215.  
  216.   If an enclosing triangle is found, Search also determines if P is
  217.   strictly Inside, or OnEdge.  If OnEdge, then V1 and V2 identify
  218.   the edge.  
  219.  
  220.   LeftNeighbor and RightNeighbor are two edges emanating from a vertex
  221.   in the triangulated region which bracket the point P.  If P is Inside
  222.   or OnEdge, then these are edges of the enclosing triangle.  If no 
  223.   enclosing triangle is found, then these are edges from the convex hull
  224.   of the triangulated region.  In the latter case, V1 identifies the
  225.   shared vertex.
  226.  
  227.   Thus, if P is known to be inside the convex hull, then Search will
  228.   always return the enclosing triangle.
  229.  
  230.   If P is OUTSIDE the convex hull, then V1 is visible from P, and the
  231.   Neighboring edges are the hull edges emanating from V1.
  232.  */
  233.  extern Triangle *Search ( );
  234.  
  235.  /*
  236.   Rho calculates 4r^2, where r is the radius of the circle passing 
  237.   through A, B, C
  238.  
  239.   The calculation of Rho depends on the observations that:
  240.  
  241.                             2
  242.                 (   abc    )     2 2 2             2
  243.    Rho = 4*r*r= (-------   )  = a b c / (a cross c)
  244.                 (a cross c )
  245.  
  246.              2   2   2
  247.    and that a , b , c , and (a cross c) are all easy to calculate.
  248.  
  249.     2
  250.    a  is the square of the length of the side opposite vertex A
  251.  
  252.     2
  253.    b  is the square of the length of the side opposite vertex B
  254.  
  255.     2
  256.    c  is the square of the length of the side opposite vertex C
  257.    
  258.  */
  259.  extern REAL Rho();
  260.  
  261.  /*
  262.    Shrink is the absolute HOT SPOT - so we've optimized as much as
  263.    is seemly.  Given a pair of triangles ABC and BAD, we want to know if
  264.    the pair ADC, BCD would be better.  If it is, then we make the switch.
  265.    
  266.    To determine which local configuration is best, we calculate
  267.  
  268.     Rabc = the square of the radius of the circle through ABC
  269.     Rbad = the square of the radius of teh circle through BAD
  270.  
  271.     Radc = the square of the radius of the circle through ADC
  272.     Rbcd = the square of the radius of the circle through BCD
  273.     
  274.     Rabc and Rbad are associated with the current configuration;
  275.     Radc and Rbcd are associated with the proposed configuration.
  276.  
  277.     The calculation of the  R's depends on the observations that:
  278.  
  279.                             2
  280.                (   abc    )
  281.         4*r*r= (-------   )
  282.                (a cross c )
  283.  
  284.     so, we get:
  285.                2 2 2              2
  286.      Rabc =   a c e  / (a cross c)
  287.  
  288.                2 2 2              2
  289.      Racd =   a d f  / (a cross d)
  290.  
  291.                2 2 2              2 
  292.      Rbad =   b d e  / (b cross d)
  293.  
  294.                2 2 2              2
  295.      Rbcd =   b c f  / (b cross c)
  296.  
  297.    where
  298.         a is the edge connecting A and C
  299.         b is the edge connecting B and D
  300.         c is the edge connecting C and B
  301.         d is the edge connecting D and A  
  302.         e is the existing diagonal  connecting A and B
  303.         f is the candidate diagonal connecting C and D
  304.  
  305.   Of the four circles, the smallest one is guaranteed to EXCLUDE the fourth 
  306.   point - that is the configuration we want.  So, we first pick the minimum
  307.   from each pair ([Rabc, Rabd] and [Racd, Rbcd]) and then compare the minima.
  308.   if min(Rabc,Rabd) > min(Racd, Rbcd) then SWAP else DO NOT SWAP.
  309.  
  310.   Finally, this criterion is not valid when the quadrilateral ADBC is not
  311.   convex. So, we first check the angles DAC and CBD: 
  312.   If (a cross d) and (b cross c) are both positive
  313.   then all is well.  If either (a cross d) of (b cross c) are negative,
  314.   then all of the above is moot, and no swap is possible.  
  315.   Here, it is important that:
  316.         
  317.          a is the DIRECTED edge from A to C
  318.          b is the DIRECTED edge from B to D
  319.          c is the DIRECTED edge from C to B
  320.          d is the DIRECTED edge from D to A
  321.  
  322.   Rabc and Rbad have already been calculated.
  323.   Radc and Rbcd need to be calculated - we do it ourself, for speed.
  324.  
  325.  */
  326.  extern void Shrink ();
  327.  
  328.  /*
  329.    Optimize is called to clean up a triangulation.
  330.  
  331.    On entry: the queue pointed at by Suspect contains "suspect" edges.  
  332.  
  333.    Optimize considers each suspect edge, in turn.  For each such edge,
  334.    the local quadrilateral is examined.  If the suspect edge is OK,
  335.    then it is simply removed from the suspect queue.  If not, then the
  336.    two triangles which share this edge are removed, and two other triangles
  337.    are inserted.  The four outer edges of these two triangles then become 
  338.    suspect, and are added to the end of the queue.
  339.  
  340.    Parameters?  WE DON'T GOT TO SHOW NO STINKING PARAMETERS!
  341.  */
  342.  extern void Optimize();
  343.  
  344.  /*
  345.    InsertPoint adds a new point.
  346.  
  347.    The current version only understands how to do that when P is
  348.    inside an already triangulated convex region.
  349.  
  350.    In case of error, return false.
  351.  
  352.    Coming soon, a version that will do on-line insertions
  353.  */ 
  354.  extern int InsertPoint ();
  355.  
  356.  /*
  357.   ConvexTriangles will triangulate a CONVEX polygon.
  358.  */
  359.  extern void ConvexTriangles();
  360.  
  361.  /*
  362.   Is A in PointSet P ?
  363.  */
  364.  extern int InPointSet();
  365.  
  366.  /*
  367.   Is A Inside Convex Polygon H
  368.   */
  369.  extern int InSideHull();
  370.  
  371.  /*
  372.     given three points (x1,y1,x2,y3,x3,y3) calculate the center
  373.     of the circle thru those points.  
  374.  
  375.      WARNING: if the points are colinear, we print a message, and 
  376.               return something reasonble.  Best not to call this routine
  377.               with colinear points, eh? 
  378.  
  379.   */
  380.  extern void VoronoiPt();
  381.  
  382. /* ------------------------------------------------------------------------ */
  383.     /*  Allocation routines  */
  384.  
  385. #define POINT_GLOB_SIZE 1024 
  386. extern PointSetPointer PointFreeList;
  387. #define NewPoint()   \
  388. (PointSetPointer) New( sizeof(PointSetElement),POINT_GLOB_SIZE,&PointFreeList )
  389. #define DisposePoint(p)  \
  390.   Dispose( p, &PointFreeList );
  391.  
  392.  
  393. #define EDGE_GLOB_SIZE 1024
  394. extern EdgePointer EdgeFreeList;
  395. #define NewEdge() \
  396.   (EdgePointer) New( sizeof(Edge), EDGE_GLOB_SIZE, &EdgeFreeList )
  397. #define DisposeEdge(e) \
  398.   Dispose( e, &EdgeFreeList );
  399.  
  400. #define TRIANGLE_GLOB_SIZE 1024
  401. extern TrianglePointer TriangleFreeList;
  402. #define NewTriangle() \
  403.   (TrianglePointer) New( sizeof(Triangle),TRIANGLE_GLOB_SIZE,&TriangleFreeList)
  404. #define DisposeTriangle(t) \
  405.   Dispose( t, &TriangleFreeList );
  406.