home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / DISCRETE.PAK / COMPUTAT.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  62.5 KB  |  1,637 lines
  1.  
  2. (*:Mathematica Version: 2.0 *)
  3.  
  4. (*:Context: DiscreteMath`ComputationalGeometry` *)
  5.  
  6. (*:Title: ComputationalGeometry *)
  7.  
  8. (*:Author:
  9.     E.C.Martin (Wolfram Research) October 1990, April 1991. 
  10. *)
  11.  
  12. (*:Legal:
  13.   Copyright (c) 1990, 1991 Wolfram Research, Inc.
  14. *)
  15.  
  16. (*:Summary:
  17.   Selected functions from computational geometry relevent to nearest
  18.   neighbor point location. *)
  19.  
  20. (*:Keywords: Delaunay triangulation, Voronoi diagram, nearest neighbor,
  21.         convex hull, computational geometry
  22. *)
  23.  
  24. (*:Package Version: 1.0 *)
  25.  
  26. (*:Requirements: No special system requirements. *)
  27.  
  28. (*:Warnings:  none. *)
  29.  
  30. (*:Limitations: Computational efficiency is maximized for non-colinear points.
  31. *)
  32.     
  33.  
  34. (*Sources:
  35.    F. R. Preparata & M. I. Shamos, Computational Geometry: An Introduction,
  36.     Springer-Verlag, 1985.
  37.    D. T. Lee & B. J. Schachter, "Two Algorithms for constructing a Delaunay
  38.     triangulation," Int. J. of Computer & Information Sciences, Vol. 9,
  39.     No. 3, pp. 219-242, 1980.
  40. *)
  41.  
  42. (*:Discussion:
  43.  
  44.    The Delaunay triangulation is represented as a vertex adjacency list
  45.     {{v1,{v11,v12,...}},...,{vn,{vn1,vn2,...}} where vi is the label of 
  46.     the ith vertex of the triangulation and indicates the position of the
  47.     point in the original coordinate list (this may not be position i if 
  48.     there are duplicates in the original list).  The list {vi1,vi2,...}
  49.     are the vertices adjacent to vertex vi, listed in counterclockwise order.
  50.     If vi is on the convex hull, then vi1 is the next counterclockwise point
  51.     on the convex hull.
  52.  
  53.    The Voronoi diagram is represented as (1) a list of Voronoi polygon vertex
  54.     coordinates and (2) a vertex adjacency list: {{q1,...,qm}, {{v1,{u11,u12,
  55.     ...}},...,{vn,{un1,un2,...}}}.  The Voronoi polygons at the periphery of
  56.     the diagram are unbounded and this means that the qj may be a true Voronoi
  57.     vertex {rj,sj} or a "quasi" vertex Ray[{rj,sj},{r1j,s1j}].  A closed 
  58.     Voronoi polygon will be defined by a list of true vertices, while an open
  59.     Voronoi polygon will be defined by a list containing two or more true 
  60.     vertices and precisely two Rays located adjacently in the list.
  61.     Ray[{rj,sj},{r1j,s1j}] indicates a ray having origin {rj,sj} (a true
  62.     vertex) and containing point {r1j,s1j}.
  63.     In the vertex adjacency list vi is the label of the ith vertex of the
  64.     dual Delaunay triangulation which is also the defining  point of the ith
  65.     Voronoi polygon.  The list {ui1,ui2,...} contains the Voronoi vertices
  66.     that compose the ith polygon, listed in counterclockwise order.
  67.  
  68. *)
  69.  
  70. BeginPackage["DiscreteMath`ComputationalGeometry`"]
  71.  
  72. DelaunayTriangulation::usage =
  73. "DelaunayTriangulation[{{x1,y1},{x2,y2},...,{xn,yn}}] yields the (planar)
  74. Delaunay triangulation of the points.  The triangulation is represented as a
  75. vertex adjacency list, one entry for each unique point in the original
  76. coordinate list indicating the adjacent vertices in counterclockwise order."
  77. (* A duplicate point is indexed according to the location of the first instance
  78.     of the point in the original input list. *)
  79. (* O(n log(n)) time *)
  80.  
  81. VoronoiDiagram::usage =
  82. "VoronoiDiagram[{{x1,y1},{x2,y2},...,{xn,yn}}] yields the (planar) Voronoi
  83. diagram of the points.  The diagram is represented as two lists: (1) a Voronoi
  84. vertex coordinate list, and (2) a vertex adjacency list, one entry for each
  85. unique point in the original coordinate list indicating the associated Voronoi
  86. polygon vertices in counterclockwise order.  VoronoiDiagram[{{x1,y1},{x2,y2},
  87. ...,{xn,yn}},val] takes val to be the vertex adjacency list of the dual 
  88. Delaunay triangulation.  VoronoiDiagram[{{x1,y1},{x2,y2},...,{xn,yn}},val,hull]
  89. takes hull to be the convex hull of the unique points." 
  90. (* Supplying more arguments to VoronoiDiagram reduces computation time. *)
  91. (* A duplicate point is indexed according to the location of the first instance
  92.     of the point in the original input list. *)
  93. (* O(n log(n)) time *)
  94.  
  95. ConvexHull::usage =
  96. "ConvexHull[{{x1,y1},{x2,y2},...,{xn,yn}}] yields the (planar) convex hull
  97. of the n points, represented as a list of point indices arranged in counter
  98. clockwise order."
  99. (* A duplicate point is indexed according to the location of the first instance
  100.     of the point in the original input list. *)
  101. (* O(n log(n)) time *)
  102.  
  103. NearestNeighbor::usage =
  104. "NearestNeighbor[{a,b},{{x1,y1},{x2,y2},...,{xn,yn}}] yields the label of the
  105. nearest neighbor of {a,b} out of the points {{x1,y1},{x2,y2},...,{xn,yn}}.
  106. NearestNeighbor[{a,b},{q1,q2,...,qm},val] yields the label of the nearest
  107. neighbor of {a,b} using the Voronoi vertices qj and the Voronoi vertex
  108. adjacency list val.  NearestNeighbor[{{a1,b1},...,{ap,bp}},{{x1,x2},...,
  109. {xn,yn}}] or NearestNeighbor[{{a1,b1},...,{ap,bp}},{q1,q2,...,qm},val] yields
  110. the nearest neighbor labels for the list {{a1,b1},...,{ap,bp}}."
  111. (* O(n log(n)) time *)
  112.  
  113. TriangularSurfacePlot::usage =
  114. "TriangularSurfacePlot[{{x1,y1,z1},...,{xn,yn,zn}}] plots the zi according 
  115. to the planar Delaunay triangulation established by the {xi,yi}.  
  116. TriangularSurfacePlot[{{x1,y1,z1},...,{xn,yn,zn}},val] plots the zi according
  117. to the planar triangulation of the {xi,yi} stipulated by the vertex adjacency
  118. list val."
  119.  
  120. PlanarGraphPlot::usage =
  121. "PlanarGraphPlot[{{x1,y1},{x2,y2},...,{xn,yn}}] plots the planar graph 
  122. corresponding to the Delaunay triangulation established by the {xi,yi}.
  123. PlanarGraphPlot[{{x1,y1},{x2,y2},...,{xn,yn}},list] plots the planar graph of 
  124. the {xi,yi} as stipulated by list; list may have the form {{l1,{l11,l12,...},
  125. ...,{ln,{ln1,ln2,...}}} (vertex adjacency list) or {l1,...,ln} (representing a 
  126. single circuit) where li is the position of the ith unique point in the input
  127. list."
  128.  
  129. DiagramPlot::usage =
  130. "DiagramPlot[{{x1,y1},{x2,y2},...,{xn,yn}}] plots the {xi,yi} and the polygons
  131. corresponding to the Voronoi diagram established by the {xi,yi}.  
  132. DiagramPlot[{{x1,y1},{x2,y2},...,{xn,yn}},{q1,q2,...,qm},val] plots the polygon
  133. 'centers' {xi,yi} and polygon vertices qj as stipulated by the vertex adjacency
  134. list val." 
  135.  
  136. DelaunayTriangulationQ::usage =
  137. "DelaunayTriangulationQ[{{x1,y1},{x2,y2},...,{xn,yn}},val] returns True if the
  138. triangulation of the {xi,yi} represented by the vertex adjacency list val is a
  139. Delaunay triangulation, and False otherwise.  DelaunayTriangulationQ[{{x1,y1},
  140. {x2,y2},...,{xn,yn}},val,hull] takes hull to be the convex hull of the unique
  141. points.  The val must be such that a point on the hull lists first the hull
  142. neighbor that follows the point on a counterclockwise traversal of the hull."
  143.  
  144. LabelPoints::usage =
  145. "LabelPoints is an option of the plotting functions within
  146. ComputationalGeometry.m.  LabelPoints->True means that the points will be
  147. labeled according to their position in the input list.  Repeated instances of
  148. a point will be plotted once and labeled according to the position of the
  149. first instance in the list."
  150.  
  151. Ray::usage =
  152. "Ray[{x1,y1},{x2,y2}] is a means of representing the infinite rays found in
  153. diagrams containing open polygons.  Here {x1,y1} indicates the head of the ray
  154. (a polygon vertex), and {x2,y2} indicates a point along the ray."
  155.  
  156. AllPoints::usage =
  157. "Allpoints is an option of ConvexHull indicating whether all distinct points
  158. on the hull or only the minimum set of points needed to define the hull are
  159. returned."
  160.  
  161. Hull::usage =
  162. "Hull is an option of DelaunayTriangulation indicating whether the convex hull
  163. is to be returned, in addition to the vertex adjacency list val describing the
  164. triangulation.  Hull->False implies that val is returned, while Hull->True
  165. implies that {val,hull} is returned."
  166.  
  167. TrimPoints::usage =
  168. "TrimPoints is an option of DiagramPlot indicating the order of the diagram
  169. oulier vertex that lies on the PlotRange limit.  The default TrimPoints -> 0
  170. implies that all diagram vertices lie within PlotRange.  TrimPoints -> n
  171. implies that PlotRange is such that the nth largest outlier vertex lies on the
  172. PlotRange limit."
  173. (* Regardless of the value of TrimPoints, PlotRange is never reduced below
  174. what is needed to plot the original set of points (i.e., polygon "centers"). *)
  175.  
  176.  
  177. Begin["Private`"]
  178.  
  179.  
  180.  
  181.  
  182. (*================================ ConvexHull ===============================*)
  183. (* The Graham scan method cannot be generalized beyond two dimensions because 
  184.    it is based on polar angle. *)
  185.  
  186. Options[ConvexHull] = {AllPoints->True} 
  187.  
  188. (* Don't use Graham scan for colinear points.  Graham scan returns the hull of
  189.     colinear points in a non-sequential order due to the polar angle sort. *)
  190. ConvexHull[set:{{_,_}..},opts___Rule]:=
  191.   Module[{distinct = Distinct[N[set]], sorted, label,
  192.         allpoints = AllPoints /. {opts} /. Options[ConvexHull]},
  193.     (* Make sure numerical coordinates are distinct and labeled
  194.        according to position in original list.  Then sort points
  195.        according to either x or y coordinate. *)
  196.         sorted = If[ distinct[[1,2,1]] != distinct[[2,2,1]],
  197.            (* points have different x coordinates *)
  198.                Sort[distinct,(#1[[2,1]] > #2[[2,1]])& ],
  199.            (* points have different y coordinates *)
  200.                Sort[distinct,(#1[[2,2]] < #2[[2,2]])& ]
  201.                  ];
  202.         label = Map[#[[1]]&,sorted];
  203.     If[!TrueQ[allpoints],
  204.       (* Return the minimum number of vertices needed to define the hull. *)
  205.       {First[label],Last[label]},
  206.       (* Return all hull vertices. *)
  207.       label
  208.         ]
  209.   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] && ColinearQ[set]
  210.  
  211.  
  212. ConvexHull[set:{{_,_}..},opts___Rule]:=
  213.     Module[{orig=N[set],sorted,label,
  214.         allpoints = AllPoints /. {opts} /. Options[ConvexHull]},
  215.     (* Make sure numerical coordinates are distinct and labeled
  216.        according to position in original list.  Then sort points
  217.        according to polar angle about the centroid *)
  218.     sorted = PolarSort[Distinct[orig]];
  219.     (* Save original labels in label and compute hull assuming points are
  220.         labeled by position in 'sorted'. *)
  221.         {label,sorted} = Transpose[sorted];
  222.     (* Compute convex hull & restore old labels indicating position
  223.         in original list. *)
  224.     Map[label[[#]]&,convexhull[sorted,allpoints]]
  225.     ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]]
  226.  
  227.  
  228. (* It is assumed that 'sorted' is a list of distinct points sorted according
  229.     to polar angle. *)
  230. convexhull[sorted:{{_,_}..},allpoints_]:=
  231.    Module[{lsort = Length[sorted], dlcl, maxx, miny, maxxlist,
  232.         start, v, successor, succsucc, out},
  233.     (* Generate doubly-linked circular list where each element is of the
  234.        form {prev,next}. *)
  235.         dlcl = Map[{pred[#,lsort],succ[#,lsort]}&,Range[lsort]];
  236.     (* Start Graham scan at rightmost point having the smallest ordinate;
  237.         this point is guarranteed on hull. *)
  238.     maxx = Apply[Max, Map[#[[1]]&,sorted]];
  239.     maxxlist = Select[sorted, (#[[1]]==maxx)&]; 
  240.         miny = Apply[Min, Map[#[[2]]&,maxxlist]];
  241.     start = Position[sorted,{maxx,miny}][[1,1]];
  242.     (* Graham scan *)
  243.     v=start; out={start};
  244.     While[dlcl[[v,2]]!=start,
  245.        successor=dlcl[[v,2]];
  246.        succsucc=dlcl[[successor,2]];
  247.        soa = SignOfArea[sorted[[v]],sorted[[successor]],
  248.             sorted[[succsucc]]];
  249.        If[ soa==1 || (soa==0 && TrueQ[allpoints]),  
  250.           (* If 'left turn' or 'no turn, but all points wanted',
  251.              advance scan, constructing output as you scan. *)
  252.           v=successor;
  253.           AppendTo[out,v],
  254.           (* Otherwise, delete successor of v by updating successor
  255.             pointer of v and updating predecessor pointer of succsucc *)
  256.           dlcl[[v,2]]=succsucc;
  257.           dlcl[[succsucc,1]]=v;
  258.           (* Also update 'out' and backtrack to predecessor of v
  259.              if (v != start). *)
  260.           If[v!=start,
  261.             out = Drop[out,-1];
  262.             v=dlcl[[v,1]] ]
  263.        ] (* end If *)
  264.     ]; (* end While *)
  265.     out
  266.    ]
  267.  
  268.  
  269. (*=============================== SignOfArea =================================*)
  270.  
  271. (* left turn if positive, right turn if negative *)
  272. SignOfArea[{x1_,y1_},{x2_,y2_},{x3_,y3_}]:=
  273.     Sign[
  274.      Chop[x1(y2-y3) - x2(y1-y3) + x3(y1-y2)]
  275.     ]        /; Apply[And,Map[NumberQ,N[ {x1,y1,x2,y2,x3,y3} ]]]
  276.           
  277.  
  278. (*=============================== ColinearQ =================================*)
  279.  
  280. ColinearQ[set:{{_,_}..}] :=
  281.   Module[{colinearq, length = Length[set], nset = N[set]},
  282.       colinearq = Scan[Module[{next = succ[#,length]},
  283.                 If[ Apply[SignOfArea,
  284.                   nset[[ {#,next,succ[next,length]} ]] ] != 0,
  285.                   Return[False]
  286.                                 ]
  287.                        ]&,
  288.                Range[length]
  289.              ];
  290.         If[ SameQ[colinearq,False],
  291.         False,
  292.         True
  293.         ]
  294.   ] /; Apply[And,Map[NumberQ,Flatten[N[set]]]]
  295.  
  296.  
  297.  
  298. (*=============================== Distinct =================================*)
  299.  
  300. (* returns an object of the form {{l1,{x1,y1}},{l2,{x2,y2}},...,{ln,{xn,yn}}} *)
  301. Distinct[orig:{{_,_}..}]:=
  302.    Module[{union,distinct={}},
  303.     (* {0,0}, {0,0.}, and {0.,0.} must be specifically mapped to {0.,0.}
  304.        before Union recognises duplicates *)
  305.     union=Union[Map[If[#=={0,0},{0.,0.},#]&,orig]];
  306.     (* Find positions of duplicates & generate unsorted list of unique
  307.     coordinates labeled according to their position in the original list. *)
  308.     Scan[
  309.       (If[#=={0.,0.},
  310.          (* {0,0}, {0,0.}, {0.,0} & {0.,0.} represents a special
  311.          case which must be checked for *)
  312.          position=Flatten[Join[
  313.         Position[orig,{0,0}],Position[orig,{0,0.}],
  314.             Position[orig,{0.,0}],Position[orig,{0.,0.}]]],
  315.          position=Flatten[Position[orig,#]]
  316.        ];
  317.        AppendTo[distinct,{position[[1]],orig[[position[[1]]]]}]
  318.       )&,
  319.       union
  320.     ];
  321.     distinct
  322.    ]
  323.  
  324.  
  325. (*================================ PolarSort ===============================*)
  326.  
  327. PolarSort[l:{{_,{_,_}}..}]:=
  328.    Module[{n=Length[l],origin,p1,p2,in,sorted},
  329.     (* The centroid of the points is interior to the convex hull. *)
  330.     origin=Apply[Plus,Map[#[[2]]&,l]]/n;
  331.     (* 1st component of elements of 'in' is label,
  332.        2nd component of elements of 'in' is centered coordinate *)
  333.     in=Map[{#[[1]],#[[2]]-origin}&,l];
  334.     (* 3rd component of elements of 'in' is polar angle *)
  335.     in = Map[Append[#,PolarAngle[#[[2]]]]&,
  336.          in];
  337.     sorted = Sort[in,
  338.                   Function[{p1,p2},
  339.                         p1[[3]]<p2[[3]] ||
  340.                        (p1[[3]]==p2[[3]] && 
  341.                        (p1[[2,1]]^2 + p1[[2,2]]^2 <  
  342.                         p2[[2,1]]^2 + p2[[2,2]]^2))
  343.               ] (* end Function *)
  344.          ]; (* end Sort *)
  345.         Map[{#[[1]],#[[2]]+origin}&,sorted]
  346.     ] /;        Apply[And, Map[NumberQ,Flatten[N[l]] ]]
  347.  
  348.  
  349. (*============================== PolarAngle ===============================*)
  350.  
  351. PolarAngle[{x_,y_}]:=
  352.     If[ y>=0,
  353.         N[Arg[x + I y]],
  354.                 N[ 2 Pi + Arg[x + I y]]
  355.         ] /; NumberQ[N[x]] && NumberQ[N[y]]
  356.  
  357.  
  358. (*============================== CommonTangents ===========================*)
  359.  
  360. (* CommonTangents[set,hullL,hullR,rmL,lmR] *)
  361. (* Returns {{l1,r1},{l2,r2}}.  l1 indicates the position in hullL of the
  362.    leftmost point of the lower common tangent (of hullL & hullR).  r1 indicates
  363.    the position in hullR of the rightmost point of the lower common tangent.
  364.    l2 indicates the position in hullL of the leftmost point of the upper common
  365.    tangent.  r2 indicates the position in hullR of the rightmost point of the
  366.    upper common tangent. *)
  367. (* Below 'rightmostLeft' is a pointer to the rightmost point in hullLeft,
  368.    while 'leftmostRight' is a pointer to the leftmost point in hullRight. *)
  369.  
  370.  
  371. CommonTangents[set:{{_,_}..}, hullLeft:{_Integer..}, hullRight:{_Integer..},
  372.            rightmostLeft_Integer, leftmostRight_Integer] :=
  373.   Module[{x, y, predx, succy, succx, predy, lowertan={}, uppertan={}, 
  374.       lengthLeft = Length[hullLeft], lengthRight = Length[hullRight],
  375.       soaR, soaL},
  376.     
  377.     (*  LOWER TAN *)
  378.     {x,y} = {rightmostLeft,leftmostRight};
  379.     {predx,succy} = {pred[x,lengthLeft],succ[y,lengthRight]};
  380.     While[ SameQ[lowertan,{}],
  381.       soaR = SignOfArea[ set[[hullLeft[[x]]]], set[[hullRight[[y]]]],
  382.             set[[hullRight[[succy]]]] ];
  383.       If[ soaR <= 0,
  384.       If[ soaR == 0, Return[comtan[set,hullLeft,hullRight]] ];
  385.           (* update y (by moving ccw on hullRight to succy) *)
  386.           {succy,y} = {succ[succy,lengthRight],succy}, 
  387.       soaL = SignOfArea[ set[[hullLeft[[predx]]]], set[[hullLeft[[x]]]],
  388.                 set[[hullRight[[y]]]] ];
  389.           If[ soaL <= 0,
  390.         If[ soaL == 0, Return[comtan[set,hullLeft,hullRight]] ];
  391.                 (* update x (by moving cw on hullLeft to predx) *)
  392.                 {predx,x} = {pred[predx,lengthLeft],predx},
  393.         (* otherwise {x,y} points to the lower common tangent *)
  394.         lowertan = {x,y}
  395.           ]
  396.       ]
  397.     ];
  398.  
  399.     (* UPPER TAN *)
  400.     {x,y} = {rightmostLeft,leftmostRight};
  401.     {succx,predy} = {succ[x,lengthLeft],pred[y,lengthRight]};
  402.     While[ SameQ[uppertan,{}], 
  403.       soaR = SignOfArea[ set[[hullLeft[[x]]]], set[[hullRight[[y]]]],
  404.              set[[hullRight[[predy]]]] ];
  405.       If[ soaR >= 0,
  406.       If[ soaR == 0, Return[comtan[set,hullLeft,hullRight]] ];
  407.           (* update y (by moving cw on hullRight to predy) *)
  408.           {predy,y} = {pred[predy,lengthRight],predy},
  409.       soaL = SignOfArea[ set[[hullLeft[[succx]]]], set[[hullLeft[[x]]]],
  410.                  set[[hullRight[[y]]]] ];
  411.           If[ soaL >= 0,
  412.         If[ soaL == 0, Return[comtan[set,hullLeft,hullRight]] ];
  413.                 (* update x (by moving ccw on hullLeft to succx) *)
  414.                 {succx,x} = {succ[succx,lengthLeft],succx},
  415.         (* otherwise {x,y} points to the upper common tangent *)
  416.         uppertan = {x,y}
  417.           ]
  418.       ]
  419.     ];
  420.  
  421.     {lowertan,uppertan}
  422.   ] 
  423.  
  424.  
  425.   (*** comtan is inefficient, but it will not be called unless points are  
  426.          colinear.  This inefficiency is not necessary. *)
  427. comtan[set:{{_,_}..}, hullLeft:{_Integer..}, hullRight:{_Integer..}] :=
  428.    Module[{hull = ConvexHull[ Join[ set[[hullLeft]], set[[hullRight]] ] ],
  429.        leng1 = Length[hullLeft], leng2 = Length[hullRight],
  430.        leng, uppertan, lowertan},
  431.       leng = Length[hull];
  432.       While[ !(hull[[1]] < leng1+1 && hull[[leng]] > leng1),
  433.         hull = RotateLeft[hull]];
  434.       (* Now the first points in hull belong to hullLeft. *)
  435.       uppertan = { hull[[1]], hull[[leng]] - leng1 }; 
  436.       While[ !(hull[[1]] > leng1 && hull[[leng]] < leng1+1),
  437.         hull = RotateLeft[hull]];
  438.       (* Now the first points in hull belong to hullRight. *)
  439.       lowertan = { hull[[leng]], hull[[1]] - leng1 };
  440.       {lowertan,uppertan}
  441.    ]
  442.     
  443.     
  444.  
  445.  
  446.  
  447. (*======================== DelaunayTriangulation ===========================*)
  448.  
  449. Options[DelaunayTriangulation] = {Hull->False}
  450.  
  451. DelaunayTriangulation::colin = "Points `` are colinear."
  452.  
  453. DelaunayTriangulation[set:{{_,_}..},opts___Rule]:=
  454.    Module[{orig=N[set], distinct = Distinct[N[set]], sorted, label, unlabeled,
  455.        delaunay, delval, hull, leftmost, rightmost,
  456.        returnhull = Hull /. {opts} /. Options[DelaunayTriangulation]},
  457.      (
  458.      (* Make sure numerical coordinates are distinct and labeled
  459.         according to position in original list.  Also, order 
  460.         according to x coordinate before recursive triangulation. *)
  461.      sorted = Sort[distinct,
  462.              (#1[[2,1]]<#2[[2,1]] ||
  463.               (#1[[2,1]]==#2[[2,1]] && #1[[2,2]]>#2[[2,2]]))&
  464.                   ];
  465.          (* Save labels, but use points sans labels in calculating
  466.           triangulation. *)
  467.          {label,unlabeled} = Transpose[sorted];
  468.      (* Recursive triangulation. *)
  469.      delaunay = Del[unlabeled];
  470.      If[ SameQ[delaunay,$Failed],
  471.         Return[$Failed],
  472.             {delval,hull,leftmost,rightmost} = delaunay
  473.      ];
  474.      (* Add vertex label field to delval and relabel delval with positions
  475.         of points in original set. *)
  476.          delval = Map[Module[{v=#[[1]],adjlist=#[[2]]},
  477.                         {label[[v]],
  478.                          Map[label[[#]]&,adjlist]}
  479.                 ]&,
  480.                       Transpose[{Range[Length[delval]],delval}]
  481.                   ];
  482.      (* Sort vertex adjacency list according to vertex label. *)
  483.      delval = Sort[delval,(#1[[1]] < #2[[1]])&];
  484.          If[returnhull,
  485.         (* Relabel convex hull with positions of points in original set. *)
  486.         hull = Map[label[[#]]&,hull];
  487.         {delval,hull},
  488.         delval
  489.          ]
  490.      ) /; Length[distinct] >= 3
  491.    ] /;    Apply[And,Map[NumberQ,N[Flatten[set]]]] 
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498. (************* 2 point triangulation *************)
  499. (***** returns {delval,convexhull,leftmost,rightmost} *****)
  500. Del[s:{{x1_,y1_},{x2_,y2_}}] :=
  501.   Module[{hull},
  502.      (* The point that has the smallest ordinate out of all rightmost points
  503.     is listed first. *)
  504.      hull = If[ x1 > x2 || (x1 == x2 && y1 < y2),
  505.         {1,2},
  506.         {2,1}];
  507.      {{{2},{1}},hull,2,1}
  508.   ]
  509.  
  510.  
  511.  
  512. (************* 3 point triangulation (non-colinear) *************)
  513. (***** returns {delval,convexhull,leftmost,rightmost} *****)
  514. Del[s:{{_,_},{_,_},{_,_}}] :=
  515.   Module[{sorted,hull,maxx,miny,maxxlist,val,leftmost},
  516.      sorted = PolarSort[Transpose[{Range[3],s}]];
  517.      {hull,sorted} = Transpose[sorted];
  518.      maxx = Apply[Max, Map[#[[1]]&,sorted]];
  519.      maxxlist = Select[sorted, (#[[1]]==maxx)&];
  520.      miny = Apply[Min, Map[#[[2]]&,maxxlist]];
  521.      (* Rotate hull so that rightmost point having the smallest ordinate
  522.         is listed first. *)
  523.      While[ s[[ First[hull] ]] != {maxx,miny} ,
  524.         hull = RotateLeft[hull]
  525.      ];
  526.      val = Map[Module[{p = Position[hull,#][[1,1]]},
  527.             Drop[ RotateLeft[hull,p], -1]
  528.                   ]&,
  529.               Range[3]];
  530.      leftmost = If[ s[[hull[[2]],1]] < s[[hull[[3]],1]] ||
  531.                 (s[[hull[[2]],1]] == s[[hull[[3]],1]] &&
  532.              s[[hull[[2]],2]] > s[[hull[[3]],2]]),
  533.             2,
  534.             3];
  535.      {val,hull,leftmost,1}
  536.   ]
  537.  
  538.  
  539.  
  540. (**************** >=6 point triangulation *****************)
  541. (***** returns {delval,convexhull,leftmost,rightmost} *****)
  542. Del[s:{{_,_}..}]:=
  543.    Module[{leng0 = Length[s], leng1, leng2, val, hull, lenghull, order,
  544.     leftptr, rightptr, xcoord, minx, maxx, lowertan, uppertan,
  545.     s1, val1, hull1, l1, r1, lenghull1, num1, part1, lpart1,
  546.     s2, val2, hull2, l2, r2, lenghull2, num2, part2, lpart2,
  547.     LTl, LTr, UTl, UTr, colinearq1 = False, colinearq2 = False},
  548.        If[ ColinearQ[s],
  549.         Message[DelaunayTriangulation::colin,s];
  550.         Return[$Failed]
  551.        ];
  552.        leng1 = Ceiling[leng0/2];
  553.      (** s[[i]] is mapped into s1[[i]] for i=1,...,leng1 *)
  554.        s1 = Take[s,leng1]; 
  555.        leng2 = leng0 - leng1;
  556.      (** s[[i]] is mapped into s2[[i-leng1]] for i=leng1+1,...,leng0 *)
  557.        s2 = Drop[s,leng1];
  558.      (** Delaunay triangulation of the subsets s1 & s2. *)
  559.      (** If either subset is colinear, need information about overall hull,
  560.         in order to compute a suitable triangulation for the colinear
  561.         subset. *)
  562.        colinearq1 = leng1 > 2 && ColinearQ[s1];
  563.        colinearq2 = leng2 > 2 && ColinearQ[s2];
  564.      (** Computing ConvexHull is time consuming; presumably there will be few
  565.         cases where (colinearq1 || colinearq2) == True. *)
  566.        If[ colinearq1 || colinearq2,
  567.     (************ "Time-consuming" procedure for colinear points. ******)
  568.         hull = ConvexHull[s];
  569.         lenghull = Length[hull];
  570.         While[ !(hull[[1]] < leng1+1 && hull[[lenghull]] > leng1),
  571.             hull = RotateLeft[hull]
  572.         ];
  573.           (** Now points in s1 listed first in hull. *)
  574.         lpart1 = Length[Select[hull,(# < leng1+1)&]];
  575.         uppertan = {hull[[1]],hull[[lenghull]]};
  576.         While[ !(hull[[1]] > leng1 && hull[[lenghull]] < leng1+1),
  577.             hull = RotateLeft[hull]
  578.         ];
  579.           (** Now points in s2 listed first in hull. *)
  580.         lpart2 = Length[Select[hull,(# > leng1)&]];
  581.         lowertan = {hull[[lenghull]],hull[[1]]};
  582.           (** Calculate ptrs to leftmost & rightmost points in hull. *)
  583.         xcoord = Map[#[[1]]&, s[[hull]] ];
  584.         {minx,maxx} = {Min[xcoord],Max[xcoord]};
  585.         leftptr = Position[xcoord,minx][[1,1]];
  586.         rightptr = Position[xcoord,maxx][[1,1]];
  587.           (** Calculate subset adjacency lists such that merge will
  588.           work properly as it adds edges from lowertan to uppertan. *)
  589.         val1 = If[ colinearq1,
  590.               If[ lpart1 == 1 ||
  591.                 (lpart1 > 1 && lowertan[[1]] != 1),
  592.                   Join[{{2}},
  593.                       Map[{#+1,#-1}&,Range[2,leng1-1]],
  594.                        {{leng1-1}}],
  595.                   Join[{{2}},
  596.                       Map[{#-1,#+1}&,Range[2,leng1-1]],
  597.                        {{leng1-1}}]
  598.               ],
  599.               Del[s1][[1]]
  600.                ];
  601.                 val2 = If[ colinearq2,
  602.               If[ lpart2 == 1 ||
  603.                 (lpart2 > 1 && lowertan[[2]] == leng1+1),
  604.                   Join[{{leng1+2}},
  605.                        Map[{#+1,#-1}&,
  606.                      Range[leng1+2,leng1+leng2-1]],
  607.                        {{leng1+leng2-1}}],
  608.                   Join[{{leng1+2}},
  609.                        Map[{#-1,#+1}&,
  610.                      Range[leng1+2,leng1+leng2-1]],
  611.                        {{leng1+leng2-1}}]
  612.               ],
  613.               Del[s2][[1]] + leng1
  614.               ],
  615.         (************ "Normal" procedure for non-colinear points. ******)
  616.                {val1,hull1,l1,r1} = Del[s1];
  617.                {val2,hull2,l2,r2} = Del[s2];
  618.                (** Offset labels in val2 & hull2 by number of points in s1. *)
  619.                val2+=leng1; hull2+=leng1;  
  620.                (** Compute ptrs to lower & upper tangents of the two hulls. *)
  621.                lenghull1 = Length[hull1];
  622.         lenghull2 = Length[hull2];
  623.                {{LTl,LTr},{UTl,UTr}} = CommonTangents[s,hull1,hull2,r1,l2];
  624.                lowertan = {hull1[[LTl]],hull2[[LTr]]};
  625.                uppertan = {hull1[[UTl]],hull2[[UTr]]};
  626.           (** num1 is # of points in hull1 also in hull *)
  627.         num1 = Mod[LTl-UTl,lenghull1]+1;
  628.           (** part1 is the part of hull1 also in hull *)
  629.         part1 = Take[RotateLeft[hull1,UTl-1],num1];
  630.         leftptr = Mod[l1-UTl,lenghull1]+1;
  631.           (** num2 is # of points in hull2 also in hull *)
  632.         num2 = Mod[UTr-LTr,lenghull2]+1;
  633.           (** part2 is the part of hull2 also in hull *)
  634.         part2 = Take[RotateLeft[hull2,LTr-1],num2];
  635.         rightptr = Mod[r2-LTr,lenghull2]+1+num1;
  636.         hull = Join[part1,part2]
  637.        ];
  638.      (** Merge vertex adjacency lists val1 & val2 into val.  *)
  639.        val = MergeDel[s,
  640.               Join[val1,val2],
  641.               lowertan,         (* lowertan points *)
  642.               uppertan          (* uppertan points *)
  643.              ];
  644.      (** Return new vertex adjacency list, convex hull, and pointers to 
  645.      leftmost & rightmost points in hull. *)
  646.        {val,hull,leftptr,rightptr}
  647.    ]
  648.  
  649.  
  650.  
  651. MergeDel[s:{{_,_}..},val:{{_Integer..}..},
  652.      {LTl_Integer,LTr_Integer},{UTl_Integer,UTr_Integer}] :=
  653.    Module[{edge = {LTl,LTr}, uppertan = {UTl,UTr}, l0 = LTl, r0 = LTr,
  654.            merge = val, l0r0r1LeftTurn, r0l0l1RightTurn,
  655.        l1, l2, r1, r2, l1val, r1val,  
  656.        r0ptr, l0ptr, r1ptr, l1ptr, r2ptr, l2ptr, lengr0, lengl0},
  657.        
  658.        (* Lengths of r0's val & l0's val AFTER first edge is added. *)
  659.        lengr0 = Length[merge[[r0]]] + 1;
  660.        lengl0 = Length[merge[[l0]]] + 1;
  661.        (* Insert lower tangent edge {l0,r0}. *)
  662.        r0ptr = 1;                   (* r0ptr will point to r0 in l0's val *)
  663.        l0ptr = lengr0;              (* l0ptr will point to l0 in r0's val *)
  664.        (* insert r0 in l0's val at position r0ptr;
  665.       insert l0 in r0's val at position l0ptr *)
  666.        {merge[[l0]],merge[[r0]]} = insert[merge[[{l0,r0}]],l0,r0,l0ptr,r0ptr];
  667.      
  668.      (* Merge the vertex adjacency lists by adding & deleting edges in val;
  669.     work from lower tangent edge to upper tangent edge. *)
  670.      While[ edge != uppertan,
  671.        (* Initially {l0,r0,r1} form a left turn & {r0,l0,l1} form a right. *)
  672.        l0r0r1LeftTurn = r0l0l1RightTurn = True;
  673.  
  674.        (********************* Derive points r1 and r2. ********************)
  675.        r1ptr = pred[l0ptr,lengr0];     (* r1ptr will point to r1 in r0's val *)
  676.        r1 = merge[[r0,r1ptr]];
  677.        If[ SignOfArea[s[[l0]],s[[r0]],s[[r1]]] == 1,  (* left turn *)
  678.       r2ptr = pred[r1ptr,lengr0]; (* r2ptr will point to r2 in r0's val *)
  679.       r2 = merge[[r0,r2ptr]];
  680.       (* Delete {r0,r1} only if edge {r1,r2} still exists. *)
  681.       While[ MemberQ[merge[[r1]],r2] &&
  682.          CircleMemberQ[ s[[l0]], s[[{r0,r1,r2}]] ],
  683.             (* Delete edge {r0,r1}. *)
  684.         (* r1ptr is position of r1 in r0's val. *)
  685.         {merge[[r0]],merge[[r1]]} = delete[merge[[{r0,r1}]],
  686.                     Position[merge[[r1]],r0][[1,1]],r1ptr];
  687.         If[ l0ptr > r1ptr, l0ptr--];    (* update l0ptr *)
  688.         If[ r2ptr > r1ptr, r2ptr--];    (* update r2ptr *)
  689.         lengr0--;                (* update lengr0 *)
  690.         r1ptr = r2ptr; r2ptr = pred[r1ptr,lengr0]; 
  691.         r1 = merge[[r0,r1ptr]]; r2 = merge[[r0,r2ptr]]
  692.           ],
  693.       l0r0r1LeftTurn = False
  694.        ];
  695.  
  696.        (********************* Derive points l1 & l2. **********************)
  697.        l1ptr = succ[r0ptr,lengl0];     (* l1ptr will point to l1 in l0's val *)
  698.        l1 = merge[[l0,l1ptr]];
  699.        If[ SignOfArea[s[[r0]],s[[l0]],s[[l1]]] == -1,  (* right turn *)
  700.       l2ptr = succ[l1ptr,lengl0]; (* l2ptr will point to l2 in l0's val *)
  701.       l2 = merge[[l0,l2ptr]];
  702.       (* Delete {l0,l1} only if edge {l1,l2} still exists. *)
  703.       While[ MemberQ[merge[[l1]],l2] &&
  704.          CircleMemberQ[ s[[r0]], s[[{l0,l1,l2}]] ],
  705.         (* Delete edge {l0,l1}. *)
  706.         {merge[[l0]],merge[[l1]]} = delete[merge[[{l0,l1}]],
  707.                     Position[merge[[l1]],l0][[1,1]],l1ptr];
  708.         If[ r0ptr > l1ptr, r0ptr--];    (* update r0ptr *)
  709.         If[ l2ptr > l1ptr, l2ptr--];    (* update l2ptr *)
  710.         lengl0--;                (* update lengl0 *)
  711.         l1ptr = l2ptr; l2ptr = succ[l1ptr,lengl0];
  712.         l1 = merge[[l0,l1ptr]]; l2 = merge[[l0,l2ptr]]
  713.           ],
  714.       r0l0l1RightTurn = False
  715.        ];
  716.  
  717.  
  718.  
  719.        (*********** Derive new l0 or new r0, & update l0ptr & r0ptr. **********)
  720.        If[ !l0r0r1LeftTurn,
  721.       (* Connect r0 to l1. *)
  722.       l1val = merge[[l1]]; lengr0++; lengl0 = Length[l1val]+1;
  723.       {l0,r0ptr} = {merge[[l0,l1ptr]],
  724.                 succ[Position[l1val,l0][[1,1]],lengl0]},
  725.       If[ !r0l0l1RightTurn,
  726.          (* Connect l0 to r1. *)
  727.          r1val = merge[[r1]]; lengl0++; lengr0 = Length[r1val]+1;
  728.          {r0,l0ptr,r0ptr} = {merge[[r0,r1ptr]],
  729.                  Position[r1val,r0][[1,1]],
  730.                  succ[r0ptr,lengl0]}, 
  731.          (* Neither {l0,r0,r1} nor {r0,l0,l1} are colinear. *)
  732.          If[ CircleMemberQ[ s[[l1]], s[[{l0,r0,r1}]] ],
  733.         (* Connect r0 to l1. *) 
  734.         l1val = merge[[l1]]; lengr0++; lengl0 = Length[l1val]+1;
  735.         {l0,r0ptr} = {merge[[l0,l1ptr]],
  736.                   succ[Position[l1val,l0][[1,1]],lengl0]},
  737.         (* Connect l0 to r1. *)
  738.         r1val = merge[[r1]]; lengl0++; lengr0 = Length[r1val]+1;
  739.         {r0,l0ptr,r0ptr} = {merge[[r0,r1ptr]],
  740.                     Position[r1val,r0][[1,1]],
  741.                     succ[r0ptr,lengl0]},
  742.              ]
  743.           ]
  744.        ];
  745.  
  746.  
  747.  
  748.        (***************** Update edge & insert into merge. ******************)
  749.        edge = {l0,r0};
  750.        (* Insert edge {l0,r0} by inserting r0 in l0's val at position r0ptr
  751.       and inserting l0 in r0's val at position l0ptr. *)
  752.        {merge[[l0]],merge[[r0]]} =
  753.       insert[ merge[[{l0,r0}]],l0,r0,l0ptr,r0ptr ]
  754.      ];
  755.      merge
  756.    ]
  757.  
  758.  
  759.         
  760.  
  761. (*================= Auxiliery Triangulation Functions ======================*)
  762.  
  763. succ[i_Integer,mod_Integer] := Mod[i,mod]+1
  764.  
  765. succ[{{i_Integer}},mod_Integer] := Mod[i,mod]+1
  766.  
  767. pred[i_Integer,mod_Integer] := Mod[i-2,mod]+1
  768.  
  769. pred[{{i_Integer}},mod_Integer] := Mod[i-2,mod]+1
  770.  
  771.  
  772.    
  773. (**** Insert v2 in position v1's adjacency list at position v2ptr. 
  774.       Insert v1 in v2's adjacency list at position v1ptr. *)
  775. (**** Returns object of the form {val1,val2}. *)
  776. insert[{val1:{_Integer..},val2:{_Integer..}},
  777.     v1_Integer, v2_Integer, v1ptr_Integer, v2ptr_Integer] :=
  778.     {Insert[val1, v2, v2ptr ],
  779.      Insert[val2, v1, v1ptr ]
  780.     } 
  781.  
  782.  
  783. (**** Delete v2 (pointed to by v2ptr) from v1's adjacency list val1.
  784.       Delete v1 (pointed to by v1ptr) from v2's adjacency list val2. *)   
  785. delete[{val1:{_Integer..},val2:{_Integer..}},
  786.     v1ptr_Integer, v2ptr_Integer] :=
  787.    {Delete[val1,v2ptr],
  788.     Delete[val2,v1ptr]
  789.    } 
  790.    
  791.  
  792. (**** The default is to not include membership in circle boundary. *)
  793.  
  794. CircleMemberQ[q:{_,_},circle:{{_,py1_},{_,py2_},{_,py3_}},
  795.         boundary_:False] :=
  796.     CircleMemberQ[{q},circle,boundary] /; !(py1 == py2 == py3)
  797.  
  798. CircleMemberQ[q:{{_,_}..},{p1:{_,py1_},p2:{_,py2_},p3:{_,py3_}},
  799.         boundary_:False] :=
  800.    Module[{x, y, u1, u2, u3, m2, m3, bx2, by2, bx3, by3, r2},
  801.     (* Insure that there is no "/0" problem. *)
  802.     {u1,u2,u3} = If[SameQ[py1,py2],
  803.                     {p3,p2,p1},
  804.                     If[SameQ[py1,py3],
  805.                        {p2,p1,p3},
  806.                        {p1,p2,p3}]];
  807.     m2 = -Apply[Divide,u1-u2]; {bx2,by2} = (u1+u2)/2;
  808.     m3 = -Apply[Divide,u1-u3]; {bx3,by3} = (u1+u3)/2;
  809.     (* Center is intersection of bisectors of two pairs of points. *)
  810.     x = (by2-by3+m3 bx3-m2 bx2)/(m3-m2);
  811.     y = ((m3 m2)(bx3-bx2) + m3 by2 - m2 by3)/(m3-m2);
  812.     r2 = (u1[[1]] - x)^2 + (u1[[2]] - y)^2;
  813.     (* Result is True if one or more points listed in q belong to the circle. *)
  814.     If[ boundary,
  815.     Apply[Or,Map[
  816.         ( Chop[(#[[1]]-x)^2+(#[[2]]-y)^2 - r2] <= 0 )&,
  817.          q]],
  818.         Apply[Or,Map[
  819.               ( Chop[(#[[1]]-x)^2+(#[[2]]-y)^2 - r2] < 0 )&,
  820.                q]]
  821.     ]
  822.   ] /; !(py1 == py2 == py3)
  823.  
  824.  
  825.  
  826. (*======================= DelaunayTriangulationQ ===========================*)
  827.  
  828. (*** It is assumed that the val is such that a vertex on the hull lists first
  829.      the neighbor that follows it on a counterclockwise traversal of the convex
  830.      hull. *)
  831.  
  832. DelaunayTriangulationQ::inval = "Triangle `` is not a valid Delaunay triangle."
  833.  
  834.  
  835. (*** Only point set and Delaunay vertex adjacency list available. ***)
  836. DelaunayTriangulationQ[set:{{_,_}..},val:{{_Integer,{_Integer..}}..}] :=
  837.   Module[{orig =  N[set], maxx, maxxlist, miny, start, hull},
  838.     (* "convexhull" is found by traversing the vertex adjacency list under
  839.        the assumption that this is faster than calculating it from scratch. *)
  840.     (* Find one point on hull, let's say rightmost point having the smallest
  841.         ordinate. *)
  842.      maxx = Apply[Max, Map[#[[1]]&,orig]];
  843.      maxxlist = Select[orig, (#[[1]]==maxx)&];
  844.      miny = Apply[Min, Map[#[[2]]&,maxxlist]];
  845.      start = Scan[(If[orig[[#[[1]]]]=={maxx,miny},
  846.              Return[#[[1]]]])&,
  847.                   val];
  848.      hull = TriangulationToHull[val,start];
  849.      If[ SameQ[hull,$Failed],
  850.         False,
  851.              DelaunayTriangulationQ[set,val,hull]
  852.      ]
  853.    ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&  (* numerical points *)
  854.     (Max[Flatten[val]] <= Length[set]) &&  (* val refers to valid pts *)
  855.     Length[val] >= 3 &&        (* val must have at least 1 triangle *)
  856.     !ColinearQ[set]            (* set must not be colinear *)
  857.  
  858.  
  859. (*** Only point set and vertex adjacency list (val) available. ***)
  860. DelaunayTriangulationQ[set:{{_,_}..},val:{{_Integer,{_Integer..}}..},
  861.                    hull:{_Integer..}]:=
  862.    Module[{orig =  N[set], circle, vertices = Map[#[[1]]&,val], delaunayq},
  863.      (* Compute list of triples defining Delaunay triangles. *)
  864.      circle = Map[
  865.          (Module[{vert1=#[[1]],vertlist=#[[2]],
  866.               vertNum=Length[#[[2]]],circle},
  867.        (* Create list of triples representing triangles adjacent to vert1. *)
  868.                      circle = Map[{vert1, vertlist[[#]],
  869.                        vertlist[[Mod[#,vertNum]+1]]}&,
  870.                            Range[vertNum]];
  871.                      If[MemberQ[hull,vert1],
  872.                 circle = Drop[circle,-1]];
  873.                      Map[Sort,circle]
  874.                  ])&,
  875.          val];
  876.     (* Eliminate duplicates. *)
  877.     circle = Union[Flatten[circle,1]];
  878.     delaunayq = If[ Apply[Or,Map[ (Apply[SignOfArea,orig[[#]]] == 0)&, circle]],
  879.             (* One or more of the triangles are colinear! *)
  880.             False,    
  881.             (* Check that the circumcircles do not include
  882.                other vertices *)
  883.                     Scan[(Module[{tri = #, compare, x},
  884.                     compare = Complement[
  885.                      Apply[Union,
  886.                       Map[Module[{x = #},
  887.                            Select[val,
  888.                              MatchQ[#,{x,_List}]&,
  889.                              1][[1,2]]
  890.                           ]&,
  891.                           tri]],
  892.                      tri];    
  893.                          If[ CircleMemberQ[orig[[compare]],orig[[tri]]], 
  894.                    Message[DelaunayTriangulationQ::inval,tri];
  895.                    Return[False]]
  896.                        ])&,
  897.                    circle]
  898.     ];
  899.     If[ SameQ[delaunayq,Null],
  900.     True,
  901.     False
  902.     ]
  903.    ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&    (* numerical points *)
  904.     (Max[Flatten[val]] <= Length[set]) &&  (* val refers to valid pts *)
  905.     (Max[hull] <= Length[set]) &&      (* hull refers to valid pts *)
  906.     Length[val] >= 3 &&        (* val includes at least 1 triangle *)
  907.     !ColinearQ[set]            (* set must not be colinear *)
  908.  
  909.  
  910.  
  911. (*====================== TriangulationToHull ==============================*)
  912.  
  913. (* It is assumed that the val is such that for a vertex on the hull, the first
  914.    neighbor listed is that one lying next on a counterclockwise traversal of
  915.    the hull. *)
  916. TriangulationToHull[val:{{_Integer,{_Integer..}}..},start_Integer]:=
  917.   Module[{label1, relabel1, start1, val1, hull, next, hullLength},
  918.     (* Relabel vertex labels to make tracing hull easier.  This wouldn't
  919.        be necessary if all the points in the original set were unique. *)
  920.      label1 = Map[#[[1]]&,val];
  921.      relabel1 = MapIndexed[(#1->#2[[1]])&,label1];
  922.      start1 = start /. relabel1;
  923.      val1 = Map[(#[[2]] /. relabel1)&,val];
  924.      hull = {start};
  925.      hullLength = 1;
  926.      next = val1[[start,1]];
  927.     (* Trace points on hull using relabeled adjacency list. *)  
  928.      While[next != start,
  929.        (* Guard against infinite loops. *)
  930.        If[ hullLength > Length[val],
  931.         Return[$Failed]];
  932.        AppendTo[hull,next];
  933.        hullLength++;
  934.        next = val1[[next,1]]
  935.      ];
  936.     (* Label convexhull with original labels. *)
  937.      label1[[hull]]
  938.   ]
  939.  
  940.  
  941. (*=========================== VoronoiDiagram ================================*)
  942.  
  943. (*** Only point set available. ***)
  944. VoronoiDiagram[set:{{_,_}..}]:=
  945.    Module[{delaunay=DelaunayTriangulation[set,Hull->True]},
  946.      If[ SameQ[delaunay,$Failed],
  947.     $Failed,
  948.         VoronoiDiagram[set,delaunay[[1]],delaunay[[2]]]
  949.      ]        /; FreeQ[delaunay,DelaunayTriangulation]
  950.    ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]]  (* numerical points *)
  951.  
  952.  
  953. (*** Both point set and Delaunay vertex adjacency list available. ***)
  954. VoronoiDiagram[set:{{_,_}..},delval:{{_Integer,{_Integer..}}..}]:=
  955.    Module[{hull = Module[{orig = N[set], maxx, maxxlist, miny, start},
  956.                       (* "convexhull" is found by traversing the Delaunay
  957.              vertex adjacency list under the assumption that
  958.                          this is faster than calculating it from scratch. *)
  959.                       (* Find one point on hull, let's say rightmost point
  960.              having the smallest ordinate. *)
  961.                       maxx = Apply[Max, Map[#[[1]]&,orig]];
  962.                       maxxlist = Select[orig, (#[[1]]==maxx)&];
  963.                       miny = Apply[Min, Map[#[[2]]&,maxxlist]];
  964.                       start = Scan[(If[orig[[#[[1]]]]=={maxx,miny},
  965.                          Return[#[[1]]]])&,
  966.                                    delval];
  967.                       TriangulationToHull[delval,start]
  968.                  ]},
  969.      If[ SameQ[hull,$Failed] || !FreeQ[hull,TriangulationToHull] ,
  970.         $Failed,
  971.              VoronoiDiagram[set,delval,hull]
  972.      ]
  973.    ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&  (* numerical points *)
  974.     (Max[Flatten[delval]] <= Length[set]) &&  (* val refers to valid pts *)
  975.     Length[delval] >= 3        (* val must have at least 1 triangle *)
  976.  
  977.  
  978. (*** Point set, Delaunay vertex adjacency list, and convexhull available. ***)
  979. VoronoiDiagram[set:{{_,_}..},delval:{{_Integer,{_Integer..}}..},
  980.            hull:{_Integer..}]:=
  981.    Module[{orig =  N[set], lhull = Length[hull], ltrue, relabel,
  982.        vorval, vorvertices, truevertices, quasivertices,
  983.        truecoordinates, quasicoordinates},
  984.     (* Compute Voronoi vertex adjacency list. *)
  985.      vorval = Map[
  986.          (Module[{delvert1=#[[1]],delvertlist=#[[2]],
  987.               delvertNum=Length[#[[2]]],vorpoly},
  988.        (*  Each vertex of a Voronoi polygon corresponds to a Delaunay vertex
  989.        triple.  Exceptions: a "quasi" Voronoi vertex corresponds to a
  990.        Delaunay vertex pair; a "true"  Voronoi vertex can correspond to
  991.        two Delaunay vertex triples. *)
  992.        (* Create list of triples representing triangles adjacent to delvert1.
  993.       These determine Voronoi vertices. *)
  994.                      vorpoly = Map[{delvert1, delvertlist[[#]],
  995.                         delvertlist[[Mod[#,delvertNum]+1]]}&,
  996.                            Range[delvertNum]];
  997.        (* If delvert1 is on hull, add information needed for determining
  998.       infinite rays of associated Voronoi polygon. *)    
  999.                      vorpoly = If[MemberQ[hull,delvert1],
  1000.                   Module[{truevert=Map[Sort,Drop[vorpoly,-1]]},
  1001.                     Join[truevert,
  1002.                          { {Sort[{vorpoly[[delvertNum,1]],
  1003.                               vorpoly[[delvertNum,2]]}],
  1004.                         truevert[[delvertNum-1]]},
  1005.                        {Sort[{vorpoly[[delvertNum,1]],
  1006.                          vorpoly[[delvertNum,3]]}],
  1007.                         truevert[[1]]}
  1008.                      } ]
  1009.                   ],
  1010.                                  Map[Sort,vorpoly]
  1011.                    ]
  1012.                  ])&,
  1013.          delval];
  1014.  
  1015.     vorvertices = Union[Flatten[vorval,1]];
  1016.     quasivertices = Take[vorvertices,lhull];
  1017.     truevertices = Drop[vorvertices,lhull];
  1018.     coordinates = Map[Apply[circlecenter,orig[[#]]]&,truevertices];
  1019.  
  1020.     (* True Voronoi vertices. *)
  1021.     truecoordinates = Union[coordinates];
  1022.     ltrue = Length[truecoordinates];
  1023.     (* lables for true Voronoi vertices *)
  1024.     relabel = Flatten[
  1025.         MapIndexed[Module[{pos = Flatten[Position[coordinates,#1],1],
  1026.                  label = #2[[1]]},
  1027.                  Map[(truevertices[[#]] -> label)&,pos]
  1028.                            ]&,
  1029.                truecoordinates],
  1030.               1];
  1031.     (* lables for quasi Voronoi vertices *)
  1032.     relabel = Join[relabel,
  1033.            MapIndexed[(#1 -> #2[[1]]+ltrue)&,
  1034.                   quasivertices
  1035.                    ]
  1036.               ];
  1037.     (* Label Voronoi vertices in Voronoi val uniquely, and eliminate
  1038.     duplicates. *)
  1039.     vorval = Map[Module[{val = # /. relabel, unique, current},
  1040.            current = val[[1]];
  1041.            unique = {current};
  1042.                Scan[(If[ # != current,
  1043.                  current = #;
  1044.                  AppendTo[unique,#]
  1045.                          ])&,
  1046.             Drop[val,1]];
  1047.                    unique
  1048.                  ]&,
  1049.          vorval];
  1050.     (* Add a field indicating Delaunay vertex label to each Voronoi val. *)
  1051.     vorval = MapIndexed[{delval[[#2[[1]],1]],#1}&,vorval];
  1052.  
  1053.     (* Quasi Voronoi vertices. *)
  1054.     (* Re-establish the ccw direction of the hull edges that define quasi
  1055.     Voronoi vertices.  (Direction was lost when unique vertices were
  1056.     found using "Union".)  Then determine ray "tail" from ray "head"
  1057.     and the associated Delaunay hull edge. *)
  1058.     quasicoordinates = Map[Module[{edge = #[[1]],
  1059.                    truevertex = #[[2]] /. relabel,
  1060.                    head, ptr1},
  1061.                ptr1 = Position[hull,edge[[1]]][[1,1]];
  1062.                If[edge[[2]] != hull[[ succ[ptr1,lhull] ]],
  1063.                     edge = edge[[{2,1}]] ];
  1064.                            head = truecoordinates[[ truevertex ]];
  1065.                Ray[head, raytail[ head, orig[[edge]] ]]
  1066.                         ]&,
  1067.             quasivertices
  1068.             ];
  1069.  
  1070.     {Join[truecoordinates,quasicoordinates],vorval}   
  1071.    ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&    (* numerical points *)
  1072.     (Max[Flatten[delval]] <= Length[set]) &&  (* val refers to valid pts *)
  1073.     (Max[hull] <= Length[set]) &&      (* hull refers to valid pts *)
  1074.     Length[delval] >= 3        (* val includes at least 1 triangle *)
  1075.  
  1076.  
  1077. (*** Compute the center of the circumcircle of triangle {p1,p2,p3}. *)
  1078. circlecenter[p1:{_,py1_}, p2:{_,py2_}, p3:{_,py3_}] :=
  1079.    Module[{u1, u2, u3, m2, m3, bx2, by2, bx3, by3},
  1080.      (* Insure that there is no "/0" problem. *)
  1081.      {u1,u2,u3} = If[SameQ[py1,py2],
  1082.             {p3,p2,p1},
  1083.             If[SameQ[py1,py3],
  1084.                 {p2,p1,p3},
  1085.                 {p1,p2,p3}]];
  1086.      m2 = -Apply[Divide,u1-u2]; {bx2,by2} = (u1+u2)/2;
  1087.      m3 = -Apply[Divide,u1-u3]; {bx3,by3} = (u1+u3)/2;
  1088.      { (by3-by2) + m2 bx2 - m3 bx3,
  1089.        m2 m3 (bx2-bx3) + m2 by3 - m3 by2} / (m2 - m3)
  1090.    ] /; !(py1 == py2 == py3)
  1091.  
  1092.  
  1093. (*** Compute the tail of the ray having head "head" and bisecting segment
  1094.      {e1,e2}.  The tail is defined by a point a distance dist(e1,e2) from
  1095.      the ray head in a direction such that {e1,head,tail} is a right turn
  1096.      OR a distance dist(e1,e2) from the edge midpoint {xm,ym} in a direction
  1097.      such that {e1,{xm,ym},tail} is a right turn.  The distance is computed
  1098.      from the head or from {xm,ym} depending on which will look better in
  1099.      DiagramPlot. *)
  1100. raytail[head:{xh_,yh_},{e1:{_,_},e2:{_,_}}] :=
  1101.    Module[{m, out, xm, ym, distfromhead, soa, 
  1102.        d0 = Sqrt[ Apply[Plus,(e1-e2)^2] ]},
  1103.       {xm,ym} = (e1+e2)/2; (* edge midpoint *)
  1104.       (* Determine whether d0 will be measured from head or {xm,ym}; 
  1105.      this is an aesthetic consideration. *)
  1106.       If[ SignOfArea[e1,{xm,ym},head] == -1,
  1107.         distfromhead = True,
  1108.         distfromhead = False;
  1109.         d0 = d0 + Sqrt[ Apply[Plus,(head - {xm,ym})^2] ]
  1110.       ];
  1111.       out = If[ e1[[2]] === e2[[2]],
  1112.             (* ray is vertical *)
  1113.             {{xh,yh-d0},{xh,yh+d0}},
  1114.             (* ray has slope m *)
  1115.         m = -Apply[Divide,e1-e2];
  1116.         Map[{#,m(#-xm)+ym}&,
  1117.                 Module[{b = (m(ym-m xm-yh)-xh)/(1+m^2),sqrt},
  1118.               sqrt = Sqrt[b^2 - (xh^2 + (ym-m xm-yh)^2 - d0^2)/(1+m^2)];
  1119.               -b + {sqrt,-sqrt}
  1120.                     ]
  1121.         ] 
  1122.       ];
  1123.       (* Choose the correct direction of the ray.
  1124.      {e1,head,tail} or {e1,{xm,ym},tail} should form a right turn. *)
  1125.       soa = If[ distfromhead,
  1126.         SignOfArea[e1,head,out[[1]]],
  1127.         SignOfArea[e1,{xm,ym},out[[1]]]
  1128.             ];
  1129.       If[ soa == -1,
  1130.          out[[1]], (* right turn *)
  1131.          out[[2]]] (* left or no turn *)
  1132.    ]
  1133.  
  1134.  
  1135.  
  1136.  
  1137. VoronoiDiagram::notimp =
  1138.   "VoronoiDiagram is not yet implemented for dimension ``."
  1139.  
  1140. VoronoiDiagram[set_List]:=
  1141.    Module[{out=Module[{d=Length[set[[1]]]},
  1142.           Message[VoronoiDiagram::notimp,d];
  1143.           $Failed]},
  1144.     out /; !SameQ[out,$Failed]
  1145.    ] /; Apply[And,Map[NumberQ,N[set]]] ||
  1146.     (Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
  1147.      Apply[Equal,Map[Length[#]&,set]] &&
  1148.      Apply[And,Map[SameQ[Head[#],List]&,set]])
  1149.  
  1150.  
  1151.  
  1152. (*=========================== NearestNeighbor ==============================*)
  1153.  
  1154. NearestNeighbor[input:{{_,_}..},vorpts_List,
  1155.         vorval:{{_Integer,{_Integer..}}..}]:=
  1156.   Module[{ninput = N[input], closed = {}, open = {}},
  1157.     (* If VoronoiDiagram produced vorval, the open polygons are guarranteed to
  1158.        be listed after the closed polygons.  The following general approach
  1159.        to identifying open and closed polygons is useful if some other method
  1160.        generated the diagram of convex polygons represented by vorval. *)
  1161.     Scan[(If[FreeQ[vorpts[[#[[2]]]],Ray],
  1162.          AppendTo[closed,#],
  1163.          AppendTo[open,#]
  1164.       ])&,
  1165.       vorval];
  1166.     (* If VoronoiDiagram produced vorval, each open polygon is guarranteed to
  1167.        list the two Ray "vertices" last.  The following scheme for putting
  1168.        open polygons in this format is useful if some other method generated
  1169.        the diagram of convex polygons represented by vorval.  A valid open
  1170.        polygon will list the two rays adjacently. *)
  1171.     open = Map[Module[{vertexlist = #[[2]], length=Length[#[[2]]]},
  1172.            While[!(SameQ[Head[vorpts[[#[[2,length]]]]],Ray] &&
  1173.                SameQ[Head[vorpts[[#[[2,length-1]]]]],Ray]),
  1174.                   vertexlist = RotateLeft[vertexlist]];
  1175.                    {#[[1]],vertexlist}
  1176.                ]&,
  1177.            open];
  1178.     (* First compute an internal point for each closed convex polygon. *)
  1179.     closed = Map[Append[#,
  1180.             Apply[Plus, vorpts[[ Take[#[[2]],3] ]] ]/3
  1181.                  ]&,
  1182.                  closed];
  1183.     (* For each input point, determine the label of the nearest neighbor. *)
  1184.     Map[Module[{query = #, class},
  1185.       (* Scan closed polygons for location of query. *)
  1186.       class = Scan[Module[{internal = #[[3]]},
  1187.              If[PolygonMemberQ[query-internal,
  1188.                        Map[(#-internal)&,
  1189.                        vorpts[[ #[[2]] ]] ]],
  1190.                 Return[#[[1]]]
  1191.                          ]
  1192.                        ]&,
  1193.                        closed];
  1194.           (* If query not located yet, scan the open polygons. *)
  1195.           If[SameQ[class,Null],
  1196.         class = Scan[Module[{ray = vorpts[[ Take[#[[2]],-2] ]]},
  1197.                            If[
  1198.                   (* Include in the open polygon those points on
  1199.                  ray forming the boundary first encountered
  1200.                  when one traverses the open polygons in a
  1201.                  counterclockwise direction. *)
  1202.                   (* left turn or no turn *)
  1203.                   SignOfArea[query,ray[[1,1]],ray[[1,2]]] >= 0 &&
  1204.                   (* right turn *)
  1205.                   SignOfArea[query,ray[[2,1]],ray[[2,2]]] == -1,
  1206.                  Return[#[[1]]]
  1207.                            ]
  1208.                          ]&,
  1209.              open]
  1210.           ];
  1211.       If[SameQ[class,Null],$Failed,class]
  1212.     ]&, 
  1213.     ninput]
  1214.   ] /; Apply[And,Map[NumberQ,N[Flatten[input]]]] &&
  1215.     ValidDiagramQ[vorpts,vorval]
  1216.  
  1217.  
  1218. NearestNeighbor[{a_,b_},vorpts_List,vorval:{{_Integer,{_Integer..}}..}]:=
  1219.   NearestNeighbor[{{a,b}},vorpts,vorval][[1]] /;
  1220.     (NumberQ[N[a]] && NumberQ[N[b]] && ValidDiagramQ[vorpts,vorval])
  1221.  
  1222. NearestNeighbor[input:{{_,_}..},pts:{{_,_}..}]:=
  1223.   Module[{voronoi,vorvertices,vorval},
  1224.     voronoi = VoronoiDiagram[pts];
  1225.     If[ SameQ[voronoi,$Failed],
  1226.     Return[$Failed]];
  1227.     {vorvertices,vorval} = voronoi;
  1228.     NearestNeighbor[input,vorvertices,vorval]
  1229.   ] /;  Apply[And,Map[NumberQ,N[Flatten[input]]]] &&
  1230.     Apply[And,Map[NumberQ,N[Flatten[pts]]]]
  1231.  
  1232. NearestNeighbor[{a_,b_},pts:{{_,_}..}]:=
  1233.   NearestNeighbor[{{a,b}},pts][[1]] /;
  1234.     NumberQ[N[a]] && NumberQ[N[b]] && Apply[And,Map[NumberQ,N[Flatten[pts]]]]
  1235.  
  1236.  
  1237.  
  1238. (*============================= PolygonMemberQ ===============================*)
  1239.  
  1240. (* Binary search of wedges that comprise a closed convex polygon containing
  1241.    the origin.  Returns True if the polygon contains the query. Polygon must
  1242.    contain origin because this makes use of PolarAngle. *)
  1243. PolygonMemberQ[query:{_,_},polygon:{{_,_}..}] :=
  1244.     Module[{median,a1,am,aq,p = Append[polygon,First[polygon]]},
  1245.       While[ Length[p]>2,
  1246.         median = Floor[Length[p]/2]+1;
  1247.     {a1,am,aq} = Map[PolarAngle,{p[[1]],p[[median]],query}];
  1248.     If[ (a1 <= aq < am) || ((a1 > am) && !(aq < a1 && aq >= am)),
  1249.        (* query lies in wedge between p[[1]] & p[[median]] *)
  1250.        p = Take[p,median],
  1251.        (* query lies in wedge between p[[median]] & p[[1]] *)
  1252.        p = Drop[p,median-1]
  1253.         ]
  1254.       ]; (* here p is of the form {{x1,y1},{x2,y2}} *)
  1255.       (* Include points on the boundary... this means that lower-labled
  1256.      polygons will be favored over higher-labled polygons when it
  1257.      comes to points on the boundary. *)
  1258.       (* test for left turn or no turn *)
  1259.       SignOfArea[query,p[[1]],p[[2]]] >= 0 
  1260.     ]
  1261.  
  1262.  
  1263. (*============================= ValidDiagramQ ===============================*)
  1264.  
  1265. ValidDiagramQ[vert_List,val:{{_Integer,{_Integer..}}..}] :=
  1266.     (* vertices must have head List or Ray *)
  1267.     Apply[And,Map[(SameQ[Head[#],List] || SameQ[Head[#],Ray])&,vert]] &&
  1268.     (* vertices must have length 2 *)
  1269.     Apply[And,Map[(Length[#]==2)&,vert]] &&
  1270.     (* each polygon must have either 0 rays (closed) or 2 rays (open) *)
  1271.      Apply[And,Map[Module[{raynum=Length[Position[ vert[[ #[[2]] ]], Ray]]},
  1272.              raynum == 0 || raynum == 2
  1273.            ]&,
  1274.            val]] &&
  1275.         (* adjacency list should not refer to vertices not in vertex list *)
  1276.     (Length[vert] == Length[Union[Flatten[Map[#[[2]]&,val]]]])
  1277.  
  1278.  
  1279.  
  1280. (*=================== Auxiliery Plotting Functions ===========================*)
  1281.  
  1282. (* In computational geometry, it's important to plot things without any
  1283.     scaling or stretching. *)
  1284. absolutePlotRange[set:{{_,_}..}] :=
  1285.   Module[{xcoord, ycoord, minx, maxx, miny, maxy,
  1286.       xhalfrange, yhalfrange, xmid, ymid},
  1287.       {xcoord,ycoord} = Transpose[set]; 
  1288.       {minx,maxx} = {Min[xcoord],Max[xcoord]};
  1289.       {miny,maxy} = {Min[ycoord],Max[ycoord]};
  1290.       {xhalfrange,yhalfrange} = {maxx-minx,maxy-miny}/2;
  1291.       If[ xhalfrange > yhalfrange,
  1292.      ymid = miny+yhalfrange;
  1293.      {{minx,maxx},{ymid-xhalfrange,ymid+xhalfrange}},
  1294.      xmid = minx+xhalfrange;
  1295.      {{xmid-yhalfrange,xmid+yhalfrange},{miny,maxy}}
  1296.       ]
  1297.  ]
  1298.  
  1299.  
  1300. (*============================= DiagramPlot ================================*)
  1301.  
  1302. Options[DiagramPlot] =  Join[{LabelPoints->True, TrimPoints->0},
  1303.                  Options[Graphics]]
  1304.  
  1305. DiagramPlot::trim =
  1306. "Warning: TrimPoints -> `` is not a nonnegative integer.  TrimPoints -> 0 used
  1307. instead."
  1308.             
  1309.  
  1310.  
  1311. (* Polygon diagram is already computed and polygon center coordinates are
  1312.    available. *)
  1313. DiagramPlot[set:{{_,_}..},vorvert_List,val:{{_Integer,{_Integer..}}..},
  1314.     opts___] :=
  1315.   Module[{labelpoints, trimpoints, orig = N[set], distinct},
  1316.     {labelpoints,trimpoints} = {LabelPoints,TrimPoints} /. {opts} /.
  1317.                     Options[DiagramPlot];
  1318.     If[ !(IntegerQ[trimpoints] && !Negative[trimpoints]),
  1319.         Message[DiagramPlot::trim,trimpoints];
  1320.         trimpoints = 0
  1321.     ];
  1322.     distinct = Map[{#[[1]],orig[[#[[1]]]]}&,val];
  1323.     diagramplot[distinct,vorvert,val,labelpoints,trimpoints,opts] 
  1324.   ] /;  Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
  1325.     ValidDiagramQ[vorvert,val] &&
  1326.     (*  adjacency list should not refer to points not in set *)
  1327.        (Max[Map[#[[1]]&,val]] <= Length[set])
  1328.  
  1329.  
  1330. (* Compute Voronoi diagram before plotting. *)
  1331. DiagramPlot[set:{{_,_}..},opts___Rule] :=
  1332.   Module[{voronoi = VoronoiDiagram[set], distinct, vorvert, val},
  1333.     If[ SameQ[voronoi,$Failed], Return[$Failed]];
  1334.     {vorvert,val} = voronoi;
  1335.     DiagramPlot[set,vorvert,val,opts]
  1336.   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] 
  1337.  
  1338.  
  1339. (* diagramplot *)
  1340. diagramplot[distinct:{{_,{_,_}}..},vorvert_List,
  1341.     val:{{_Integer,{_Integer..}}..},
  1342.     labelpoints_Symbol,trimpoints_Integer,opts___] :=
  1343.   Module[{centerlist,vertexlist={PointSize[.012]},
  1344.       linelist={Thickness[.003]}, leng = Length[distinct],
  1345.       trim = trimpoints, original, absolutepoints, centroid,
  1346.       max2, dist2, pos, xmin, xmax, ymin, ymax, delx, dely},
  1347.     (* diagram polygon "centers" *)
  1348.     centerlist = If[labelpoints,
  1349.                     Map[Apply[Text,#]&,distinct],
  1350.                      Prepend[Map[Point,Map[#[[2]]&,distinct]],
  1351.               PointSize[.012]]
  1352.          ];
  1353.     (* diagram vertices and infinite rays *)
  1354.     Scan[(If[FreeQ[#,Ray],
  1355.         AppendTo[vertexlist,Point[#]],
  1356.         AppendTo[linelist,Line[Apply[List,#]]],
  1357.       ])&,
  1358.      vorvert];
  1359.     (* closed diagram polygons *)
  1360.     Scan[(Module[{list = vorvert[[#[[2]]]], closedpoly},
  1361.          closedpoly = FreeQ[list,Ray];
  1362.          list = Select[list,(!SameQ[Head[#],Ray])&];
  1363.          If[ Length[list] > 1, (* val contains at least one segment *)
  1364.         If[ closedpoly,
  1365.           (* add in the segment that closes the poly *)
  1366.               AppendTo[linelist,Line[Join[list,{list[[1]]}]]],
  1367.           (* just add in explicit segments *)
  1368.           AppendTo[linelist,Line[list]]
  1369.                 ]
  1370.              ]
  1371.           ])&,
  1372.      val];
  1373.     (* Compute points that will determine PlotRange. *)
  1374.     original = Map[#[[2]]&,distinct];
  1375.     If[ trim == 0,
  1376.     (* include all diagram vertices AND ray tails *)
  1377.     absolutepoints = Join[ original,
  1378.                    Map[(If[ SameQ[Head[#],Ray],
  1379.                     #[[2]],
  1380.                     #])&,
  1381.                        vorvert]
  1382.                              ],
  1383.         (* exclude ray tails and (trim-1) of the furthest outlying diagram
  1384.        vertices *)
  1385.         absolutepoints = Join[ original,
  1386.                    Select[vorvert, !SameQ[Head[#],Ray]&] ];
  1387.         If[ trim > 1,
  1388.             centroid = Apply[Plus,original]/Length[original];
  1389.         max2 = Max[Map[Apply[Plus,(#-centroid)^2]&,
  1390.                        original
  1391.                           ]];
  1392.             dist2 = Map[Apply[Plus,(#-centroid)^2]&,
  1393.                     Drop[absolutepoints,leng]
  1394.                         ];
  1395.             pos = Position[dist2,Max[dist2]][[1,1]];
  1396.                 (* Delete diagram vertices from 'absolutepoints' as long as 
  1397.            they lie further from the centroid then does the furthest
  1398.            polygon "center". *)
  1399.                  While[ trim!=1 && dist2[[pos]] > max2,
  1400.             dist2 = Delete[dist2,pos];
  1401.             absolutepoints = Delete[absolutepoints, pos+leng];
  1402.             pos = Position[dist2,Max[dist2]][[1,1]];
  1403.             trim--
  1404.                  ]
  1405.     ]
  1406.     ];
  1407.     {{xmin,xmax},{ymin,ymax}} = absolutePlotRange[absolutepoints];
  1408.     {delx,dely} = .05/1.9 {xmax-xmin,ymax-ymin};
  1409.     optslist = Select[{opts},!(SameQ[#[[1]],LabelPoints] ||
  1410.                    SameQ[#[[1]],TrimPoints])&];
  1411.     Show[Graphics[
  1412.         {centerlist,
  1413.          vertexlist,
  1414.          linelist},
  1415.         PlotRange -> {{xmin-delx,xmax+delx},{ymin-dely,ymax+dely}},
  1416.     AspectRatio->1
  1417.     ],optslist]
  1418.   ]
  1419.  
  1420.  
  1421.  
  1422. (*=========================== PlanarGraphPlot ===============================*)
  1423.  
  1424. Options[PlanarGraphPlot] = Join[{LabelPoints->True},
  1425.                 Options[Graphics]]
  1426.  
  1427.  
  1428. (* This is the form for plotting convex hull. *)
  1429. PlanarGraphPlot[set:{{_,_}..}, circuit:{_Integer..}, opts___]:=
  1430.    Module[{orig = N[set], pointlist, linelist, distinct,
  1431.        xmax, xmin, ymax, ymin, delx, dely, optslist,
  1432.        labelpoints = LabelPoints /. {opts} /. Options[PlanarGraphPlot]},
  1433.       If[labelpoints,
  1434.      distinct = Distinct[orig];
  1435.      pointlist = Map[Text[#,orig[[#]]]&,Transpose[distinct][[1]]],
  1436.      pointlist = Prepend[Map[Point,orig],PointSize[.012]] ];
  1437.       linelist = {Thickness[.003],
  1438.           Line[Append[ orig[[circuit]],orig[[circuit[[1]]]] ]]};
  1439.       {{xmin,xmax},{ymin,ymax}} = absolutePlotRange[orig];
  1440.       {delx,dely} = .05/1.9 {xmax-xmin,ymax-ymin};
  1441.       optslist = Select[{opts},!SameQ[#[[1]],LabelPoints]&];
  1442.       Show[Graphics[
  1443.       {pointlist,
  1444.        linelist},
  1445.       PlotRange -> {{xmin-delx,xmax+delx},{ymin-dely,ymax+dely}},
  1446.       AspectRatio->1
  1447.       ],optslist]
  1448.    ]    /;  (Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
  1449.          (Max[circuit] <= Length[set]))
  1450.  
  1451.  
  1452. (* Vertex adjacency list is already computed. *)
  1453. PlanarGraphPlot[set:{{_,_}..},val:{{_Integer,{_Integer..}}..},opts___] :=
  1454.   Module[{orig = N[set], pointlist, linelist, distinct, pairs, 
  1455.       xmax, xmin, ymax, ymin, delx, dely, optslist,
  1456.       labelpoints = LabelPoints /. {opts} /. Options[PlanarGraphPlot]},
  1457.       If[labelpoints,
  1458.      distinct = Distinct[orig];
  1459.      pointlist = Map[Text[#,orig[[#]]]&,Transpose[distinct][[1]]],
  1460.      pointlist = Prepend[Map[Point,orig],PointSize[.012]] ];
  1461.       pairs = Map[(Outer[List,{#[[1]]},#[[2]]][[1]])&,
  1462.               val];
  1463.       pairs = Union[Map[Sort,Flatten[pairs,1]]];
  1464.       linelist = Prepend[Map[Line[orig[[#]]]&,pairs],Thickness[.003]];
  1465.       {{xmin,xmax},{ymin,ymax}} = absolutePlotRange[orig];
  1466.       {delx,dely} = .05/1.9 {xmax-xmin,ymax-ymin};
  1467.       optslist = Select[{opts},!SameQ[#[[1]],LabelPoints]&];
  1468.       Show[Graphics[
  1469.       {pointlist,
  1470.        linelist},
  1471.       PlotRange -> {{xmin-delx,xmax+delx},{ymin-dely,ymax+dely}},
  1472.       AspectRatio->1
  1473.       ],optslist]
  1474.   ]   /;   Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
  1475.             (Max[Flatten[val]] <= Length[set])
  1476.  
  1477.  
  1478. (* Compute Delaunay triangulation before plotting. *)
  1479. PlanarGraphPlot[set:{{_,_}..},opts___] :=
  1480.   Module[{delaunay = DelaunayTriangulation[set]},
  1481.     If[ SameQ[delaunay,$Failed],
  1482.     $Failed,
  1483.         PlanarGraphPlot[set,delaunay,opts]
  1484.     ]
  1485.   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]]
  1486.  
  1487.  
  1488.  
  1489.  
  1490.  
  1491.  
  1492. (*======================= TriangularSurfacePlot =============================*)
  1493.  
  1494. Options[TriangularSurfacePlot] = Options[Graphics3D]
  1495.  
  1496.  
  1497. (* Vertex adjacency list is already computed. *)
  1498. TriangularSurfacePlot[set:{{_,_,_}..},val:{{_Integer,{_Integer..}}..},
  1499.     opts___] :=
  1500.   Module[{orig = N[set], distinct, label, hull, polygonlist},
  1501.     (* If two or more 3d points have the same {x,y} coordinates, but different
  1502.        z coordinates, only the first instance of the {x,y} coordinates in the
  1503.        original set is considered.  Other values for z are ignored. *)
  1504.      distinct = Map[{#[[1]],orig[[ #[[1]] ]]}&,val];
  1505.     (* In the case of Delaunay triangulation, convexhull may be extracted by
  1506.        traversing the vertex adjacency list val.  However, since val may
  1507.        represent some other triangulation, the convex hull is computed from
  1508.        scratch. *)
  1509.      label = Map[#[[1]]&,val];
  1510.      hull = Map[ label[[#]]&,
  1511.                convexhull[Map[Drop[#[[2]],-1]&,distinct],True]
  1512.          ];
  1513.     (* find all triangles surrounding interior points (i.e. not on hull) *)
  1514.     polygonlist = Map[Module[{vertex=#[[1]],adjvert=#[[2]],
  1515.                   vertNum=Length[#[[2]]]},
  1516.             Map[{vertex,adjvert[[#]],
  1517.                  adjvert[[ Mod[#,vertNum]+1 ]]}&,
  1518.                  Range[vertNum]]
  1519.                       ]&,
  1520.               Select[val,!MemberQ[hull,#[[1]]]&]
  1521.           ];
  1522.     (* Map Polygon onto coordinates of unique triangles  *)
  1523.     polygonlist = Map[Polygon[Module[{triangle = #},
  1524.                     Map[distinct[[#,2]]&,
  1525.                         triangle]
  1526.                               ]]&,
  1527.               Union[Map[Sort,Flatten[polygonlist,1]]]
  1528.                   ];
  1529.     Show[Graphics3D[polygonlist],opts]
  1530.   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
  1531.     (Max[Flatten[val]] <= Length[set])
  1532.  
  1533. (* Compute Delaunay triangulation before plotting. *)
  1534. TriangularSurfacePlot[set:{{_,_,_}..},opts___] :=
  1535.   Module[{delaunay = DelaunayTriangulation[Map[Drop[#,-1]&,set]]},
  1536.     If[ SameQ[delaunay,$Failed],
  1537.     $Failed,
  1538.         TriangularSurfacePlot[set,delaunay,opts]
  1539.     ]
  1540.   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]]
  1541.  
  1542.  
  1543.  
  1544.  
  1545.  
  1546. End[]
  1547.  
  1548. EndPackage[] 
  1549.  
  1550.  
  1551. (* :Examples:
  1552.  
  1553. input = {{4.4,14},{6.7,15.25},{6.9,12.8},{2.1,11.1},{9.5,14.9},{13.2,11.9},
  1554.     {10.3,12.3},{6.8,9.5},{3.3,7.7},{.6,5.1},{5.3,2.4},{8.45,4.7},
  1555.     {11.5,9.6},{13.8,7.3},{12.9,3.1},{11,1.1}};
  1556. query = Map[(#+{.001,.001})&,input]
  1557. input3D = Map[{#[[1]],#[[2]],Sqrt[64-(#[[1]]-8)^2-(#[[2]]-8)^2]}&,input]
  1558.  
  1559. *****************************************************************************
  1560.   convex hull of "input":
  1561.         ConvexHull[input]
  1562.  
  1563.   plot convex hull of "input":
  1564.         PlanarGraphPlot[input,ConvexHull[input]]
  1565.  
  1566. *****************************************************************************
  1567.   Delaunay triangulation of "input":
  1568.         val = DelaunayTriangulation[input]
  1569.  
  1570.   plot Delaunay triangulation of "input" where the Delaunay vertex adjacency
  1571.     list is not available:
  1572.         PlanarGraphPlot[input]
  1573.  
  1574.   plot triangulation of "input" where the triangulation vertex adjacency list
  1575.     is "val" ("val" need not represent a Delaunay triangulation):
  1576.         PlanarGraphPlot[input,val]
  1577.  
  1578.   plot triangular surface where the vertex adjacency list representing the
  1579.     Delaunay triangulation of the projected input is not available:
  1580.         TriangularSurfacePlot[input3D]
  1581.  
  1582.   plot triangular surface where the vertex adjacency list representing the
  1583.     triangulation of the projected input is "val" ("val" need not
  1584.     represent a Delaunay triangulation):
  1585.         val = DelaunayTriangulation[Map[Drop[#,-1]&,input3D]];
  1586.         TriangularSurfacePlot[input3D,val]
  1587.  
  1588. *****************************************************************************
  1589.   Voronoi diagram of "input":
  1590.         {diagvert,diagval} = VoronoiDiagram[input]
  1591.  
  1592.   Voronoi diagram of "input" where the vertex adjacency list representing the
  1593.     dual Delaunay triangulation is "val":
  1594.         VoronoiDiagram[input,val]
  1595.  
  1596.   Voronoi diagram of "input" where the vertex adjacency list representing the
  1597.     dual Delaunay triangulation is "val" and the convex hull of "input"
  1598.     is "hull":
  1599.         {val,hull} = DelaunayTriangulation[input,Hull->True];
  1600.         VoronoiDiagram[input,val,hull]
  1601.  
  1602.   plot Voronoi diagram of "input" where the Voronoi vertices and vertex
  1603.     adjacency list are not available:
  1604.         DiagramPlot[input]
  1605.  
  1606.   plot diagram of "input" where the diagram vertices are "diagvert" and the
  1607.     diagram vertex adjacency list is "diagval" ("diagvert" and "diagval"
  1608.     need not represent a Voronoi diagram):
  1609.         DiagramPlot[input,diagvert,diagval]
  1610.  
  1611.   nearest neighbors of "query" where the Voronoi diagram vertices and vertex
  1612.     adjacency list associated with "input" are not available:
  1613.         NearestNeighbor[query,input]
  1614.  
  1615.   nearest neighbors of "query" where the Voronoi diagram of "input" is
  1616.     represented by the vertices "diagvert" and the vertex adjacency list
  1617.     "diagval":
  1618.         NearestNeighbor[query,diagvert,diagval]
  1619.  
  1620. *****************************************************************************
  1621.   Delaunay triangulation query:
  1622.     DelaunayTriangulationQ[input,val]
  1623.  
  1624.     notdelaunayval = Join[{{1,{4,2}},{2,{1,4,3,5}},{3,{2,4,8,7,5}},
  1625.                 {4,{10,9,3,2,1}}}, Drop[val,4]];
  1626.     DelaunayTriangulationQ[input,notdelaunayval]
  1627.  
  1628.   Plotting a non-Delaunay triangulation:
  1629.     PlanarGraphPlot[input,notdelaunayval]
  1630.  
  1631. *)
  1632.  
  1633.  
  1634.  
  1635.  
  1636.  
  1637.