
(*:Mathematica Version: 2.0 *)

(*:Context: DiscreteMath`ComputationalGeometry` *)

(*:Title: ComputationalGeometry *)

(*:Author:
	E.C.Martin (Wolfram Research) October 1990, April 1991. 
*)

(*:Legal:
  Copyright (c) 1990, 1991 Wolfram Research, Inc.
*)

(*:Summary:
  Selected functions from computational geometry relevent to nearest
  neighbor point location. *)

(*:Keywords: Delaunay triangulation, Voronoi diagram, nearest neighbor,
		convex hull, computational geometry
*)

(*:Package Version: 1.0 *)

(*:Requirements: No special system requirements. *)

(*:Warnings:  none. *)

(*:Limitations: Computational efficiency is maximized for non-colinear points.
*)
	

(*Sources:
   F. R. Preparata & M. I. Shamos, Computational Geometry: An Introduction,
	Springer-Verlag, 1985.
   D. T. Lee & B. J. Schachter, "Two Algorithms for constructing a Delaunay
	triangulation," Int. J. of Computer & Information Sciences, Vol. 9,
	No. 3, pp. 219-242, 1980.
*)

(*:Discussion:

   The Delaunay triangulation is represented as a vertex adjacency list
    {{v1,{v11,v12,...}},...,{vn,{vn1,vn2,...}} where vi is the label of 
    the ith vertex of the triangulation and indicates the position of the
    point in the original coordinate list (this may not be position i if 
    there are duplicates in the original list).  The list {vi1,vi2,...}
    are the vertices adjacent to vertex vi, listed in counterclockwise order.
    If vi is on the convex hull, then vi1 is the next counterclockwise point
    on the convex hull.

   The Voronoi diagram is represented as (1) a list of Voronoi polygon vertex
    coordinates and (2) a vertex adjacency list: {{q1,...,qm}, {{v1,{u11,u12,
    ...}},...,{vn,{un1,un2,...}}}.  The Voronoi polygons at the periphery of
    the diagram are unbounded and this means that the qj may be a true Voronoi
    vertex {rj,sj} or a "quasi" vertex Ray[{rj,sj},{r1j,s1j}].  A closed 
    Voronoi polygon will be defined by a list of true vertices, while an open
    Voronoi polygon will be defined by a list containing two or more true 
    vertices and precisely two Rays located adjacently in the list.
    Ray[{rj,sj},{r1j,s1j}] indicates a ray having origin {rj,sj} (a true
    vertex) and containing point {r1j,s1j}.
    In the vertex adjacency list vi is the label of the ith vertex of the
    dual Delaunay triangulation which is also the defining  point of the ith
    Voronoi polygon.  The list {ui1,ui2,...} contains the Voronoi vertices
    that compose the ith polygon, listed in counterclockwise order.

*)

BeginPackage["DiscreteMath`ComputationalGeometry`"]

DelaunayTriangulation::usage =
"DelaunayTriangulation[{{x1,y1},{x2,y2},...,{xn,yn}}] yields the (planar)
Delaunay triangulation of the points.  The triangulation is represented as a
vertex adjacency list, one entry for each unique point in the original
coordinate list indicating the adjacent vertices in counterclockwise order."
(* A duplicate point is indexed according to the location of the first instance
	of the point in the original input list. *)
(* O(n log(n)) time *)

VoronoiDiagram::usage =
"VoronoiDiagram[{{x1,y1},{x2,y2},...,{xn,yn}}] yields the (planar) Voronoi
diagram of the points.  The diagram is represented as two lists: (1) a Voronoi
vertex coordinate list, and (2) a vertex adjacency list, one entry for each
unique point in the original coordinate list indicating the associated Voronoi
polygon vertices in counterclockwise order.  VoronoiDiagram[{{x1,y1},{x2,y2},
...,{xn,yn}},val] takes val to be the vertex adjacency list of the dual 
Delaunay triangulation.  VoronoiDiagram[{{x1,y1},{x2,y2},...,{xn,yn}},val,hull]
takes hull to be the convex hull of the unique points." 
(* Supplying more arguments to VoronoiDiagram reduces computation time. *)
(* A duplicate point is indexed according to the location of the first instance
	of the point in the original input list. *)
(* O(n log(n)) time *)

ConvexHull::usage =
"ConvexHull[{{x1,y1},{x2,y2},...,{xn,yn}}] yields the (planar) convex hull
of the n points, represented as a list of point indices arranged in counter
clockwise order."
(* A duplicate point is indexed according to the location of the first instance
	of the point in the original input list. *)
(* O(n log(n)) time *)

NearestNeighbor::usage =
"NearestNeighbor[{a,b},{{x1,y1},{x2,y2},...,{xn,yn}}] yields the label of the
nearest neighbor of {a,b} out of the points {{x1,y1},{x2,y2},...,{xn,yn}}.
NearestNeighbor[{a,b},{q1,q2,...,qm},val] yields the label of the nearest
neighbor of {a,b} using the Voronoi vertices qj and the Voronoi vertex
adjacency list val.  NearestNeighbor[{{a1,b1},...,{ap,bp}},{{x1,x2},...,
{xn,yn}}] or NearestNeighbor[{{a1,b1},...,{ap,bp}},{q1,q2,...,qm},val] yields
the nearest neighbor labels for the list {{a1,b1},...,{ap,bp}}."
(* O(n log(n)) time *)

TriangularSurfacePlot::usage =
"TriangularSurfacePlot[{{x1,y1,z1},...,{xn,yn,zn}}] plots the zi according 
to the planar Delaunay triangulation established by the {xi,yi}.  
TriangularSurfacePlot[{{x1,y1,z1},...,{xn,yn,zn}},val] plots the zi according
to the planar triangulation of the {xi,yi} stipulated by the vertex adjacency
list val."

PlanarGraphPlot::usage =
"PlanarGraphPlot[{{x1,y1},{x2,y2},...,{xn,yn}}] plots the planar graph 
corresponding to the Delaunay triangulation established by the {xi,yi}.
PlanarGraphPlot[{{x1,y1},{x2,y2},...,{xn,yn}},list] plots the planar graph of 
the {xi,yi} as stipulated by list; list may have the form {{l1,{l11,l12,...},
...,{ln,{ln1,ln2,...}}} (vertex adjacency list) or {l1,...,ln} (representing a 
single circuit) where li is the position of the ith unique point in the input
list."

DiagramPlot::usage =
"DiagramPlot[{{x1,y1},{x2,y2},...,{xn,yn}}] plots the {xi,yi} and the polygons
corresponding to the Voronoi diagram established by the {xi,yi}.  
DiagramPlot[{{x1,y1},{x2,y2},...,{xn,yn}},{q1,q2,...,qm},val] plots the polygon
'centers' {xi,yi} and polygon vertices qj as stipulated by the vertex adjacency
list val." 

DelaunayTriangulationQ::usage =
"DelaunayTriangulationQ[{{x1,y1},{x2,y2},...,{xn,yn}},val] returns True if the
triangulation of the {xi,yi} represented by the vertex adjacency list val is a
Delaunay triangulation, and False otherwise.  DelaunayTriangulationQ[{{x1,y1},
{x2,y2},...,{xn,yn}},val,hull] takes hull to be the convex hull of the unique
points.  The val must be such that a point on the hull lists first the hull
neighbor that follows the point on a counterclockwise traversal of the hull."

LabelPoints::usage =
"LabelPoints is an option of the plotting functions within
ComputationalGeometry.m.  LabelPoints->True means that the points will be
labeled according to their position in the input list.  Repeated instances of
a point will be plotted once and labeled according to the position of the
first instance in the list."

Ray::usage =
"Ray[{x1,y1},{x2,y2}] is a means of representing the infinite rays found in
diagrams containing open polygons.  Here {x1,y1} indicates the head of the ray
(a polygon vertex), and {x2,y2} indicates a point along the ray."

AllPoints::usage =
"Allpoints is an option of ConvexHull indicating whether all distinct points
on the hull or only the minimum set of points needed to define the hull are
returned."

Hull::usage =
"Hull is an option of DelaunayTriangulation indicating whether the convex hull
is to be returned, in addition to the vertex adjacency list val describing the
triangulation.  Hull->False implies that val is returned, while Hull->True
implies that {val,hull} is returned."

TrimPoints::usage =
"TrimPoints is an option of DiagramPlot indicating the order of the diagram
oulier vertex that lies on the PlotRange limit.  The default TrimPoints -> 0
implies that all diagram vertices lie within PlotRange.  TrimPoints -> n
implies that PlotRange is such that the nth largest outlier vertex lies on the
PlotRange limit."
(* Regardless of the value of TrimPoints, PlotRange is never reduced below
what is needed to plot the original set of points (i.e., polygon "centers"). *)


Begin["Private`"]




(*================================ ConvexHull ===============================*)
(* The Graham scan method cannot be generalized beyond two dimensions because 
   it is based on polar angle. *)

Options[ConvexHull] = {AllPoints->True} 

(* Don't use Graham scan for colinear points.  Graham scan returns the hull of
    colinear points in a non-sequential order due to the polar angle sort. *)
ConvexHull[set:{{_,_}..},opts___Rule]:=
  Module[{distinct = Distinct[N[set]], sorted, label,
		allpoints = AllPoints /. {opts} /. Options[ConvexHull]},
	(* Make sure numerical coordinates are distinct and labeled
	   according to position in original list.  Then sort points
	   according to either x or y coordinate. *)
        sorted = If[ distinct[[1,2,1]] != distinct[[2,2,1]],
		   (* points have different x coordinates *)
	           Sort[distinct,(#1[[2,1]] > #2[[2,1]])& ],
		   (* points have different y coordinates *)
	           Sort[distinct,(#1[[2,2]] < #2[[2,2]])& ]
                 ];
        label = Map[#[[1]]&,sorted];
	If[!TrueQ[allpoints],
	  (* Return the minimum number of vertices needed to define the hull. *)
	  {First[label],Last[label]},
	  (* Return all hull vertices. *)
	  label
        ]
  ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] && ColinearQ[set]


ConvexHull[set:{{_,_}..},opts___Rule]:=
	Module[{orig=N[set],sorted,label,
		allpoints = AllPoints /. {opts} /. Options[ConvexHull]},
	(* Make sure numerical coordinates are distinct and labeled
	   according to position in original list.  Then sort points
	   according to polar angle about the centroid *)
	sorted = PolarSort[Distinct[orig]];
	(* Save original labels in label and compute hull assuming points are
		labeled by position in 'sorted'. *)
        {label,sorted} = Transpose[sorted];
	(* Compute convex hull & restore old labels indicating position
		in original list. *)
	Map[label[[#]]&,convexhull[sorted,allpoints]]
	] /; Apply[And,Map[NumberQ,N[Flatten[set]]]]


(* It is assumed that 'sorted' is a list of distinct points sorted according
	to polar angle. *)
convexhull[sorted:{{_,_}..},allpoints_]:=
   Module[{lsort = Length[sorted], dlcl, maxx, miny, maxxlist,
		start, v, successor, succsucc, out},
	(* Generate doubly-linked circular list where each element is of the
	   form {prev,next}. *)
        dlcl = Map[{pred[#,lsort],succ[#,lsort]}&,Range[lsort]];
	(* Start Graham scan at rightmost point having the smallest ordinate;
		this point is guarranteed on hull. *)
	maxx = Apply[Max, Map[#[[1]]&,sorted]];
	maxxlist = Select[sorted, (#[[1]]==maxx)&]; 
        miny = Apply[Min, Map[#[[2]]&,maxxlist]];
	start = Position[sorted,{maxx,miny}][[1,1]];
	(* Graham scan *)
	v=start; out={start};
	While[dlcl[[v,2]]!=start,
	   successor=dlcl[[v,2]];
	   succsucc=dlcl[[successor,2]];
	   soa = SignOfArea[sorted[[v]],sorted[[successor]],
			sorted[[succsucc]]];
	   If[ soa==1 || (soa==0 && TrueQ[allpoints]),  
		  (* If 'left turn' or 'no turn, but all points wanted',
		     advance scan, constructing output as you scan. *)
		  v=successor;
		  AppendTo[out,v],
		  (* Otherwise, delete successor of v by updating successor
		    pointer of v and updating predecessor pointer of succsucc *)
		  dlcl[[v,2]]=succsucc;
		  dlcl[[succsucc,1]]=v;
		  (* Also update 'out' and backtrack to predecessor of v
		     if (v != start). *)
		  If[v!=start,
			out = Drop[out,-1];
			v=dlcl[[v,1]] ]
	   ] (* end If *)
	]; (* end While *)
	out
   ]
 

(*=============================== SignOfArea =================================*)

(* left turn if positive, right turn if negative *)
SignOfArea[{x1_,y1_},{x2_,y2_},{x3_,y3_}]:=
    Sign[
     Chop[x1(y2-y3) - x2(y1-y3) + x3(y1-y2)]
    ]		/; Apply[And,Map[NumberQ,N[ {x1,y1,x2,y2,x3,y3} ]]]
	      

(*=============================== ColinearQ =================================*)

ColinearQ[set:{{_,_}..}] :=
  Module[{colinearq, length = Length[set], nset = N[set]},
      colinearq = Scan[Module[{next = succ[#,length]},
				If[ Apply[SignOfArea,
				  nset[[ {#,next,succ[next,length]} ]] ] != 0,
				  Return[False]
                                ]
                       ]&,
		       Range[length]
             ];
        If[ SameQ[colinearq,False],
		False,
		True
        ]
  ] /; Apply[And,Map[NumberQ,Flatten[N[set]]]]



(*=============================== Distinct =================================*)

(* returns an object of the form {{l1,{x1,y1}},{l2,{x2,y2}},...,{ln,{xn,yn}}} *)
Distinct[orig:{{_,_}..}]:=
   Module[{union,distinct={}},
	(* {0,0}, {0,0.}, and {0.,0.} must be specifically mapped to {0.,0.}
	   before Union recognises duplicates *)
	union=Union[Map[If[#=={0,0},{0.,0.},#]&,orig]];
	(* Find positions of duplicates & generate unsorted list of unique
	coordinates labeled according to their position in the original list. *)
	Scan[
	  (If[#=={0.,0.},
	     (* {0,0}, {0,0.}, {0.,0} & {0.,0.} represents a special
	     case which must be checked for *)
	     position=Flatten[Join[
		Position[orig,{0,0}],Position[orig,{0,0.}],
	        Position[orig,{0.,0}],Position[orig,{0.,0.}]]],
	     position=Flatten[Position[orig,#]]
	   ];
	   AppendTo[distinct,{position[[1]],orig[[position[[1]]]]}]
	  )&,
	  union
	];
	distinct
   ]


(*================================ PolarSort ===============================*)

PolarSort[l:{{_,{_,_}}..}]:=
   Module[{n=Length[l],origin,p1,p2,in,sorted},
	(* The centroid of the points is interior to the convex hull. *)
	origin=Apply[Plus,Map[#[[2]]&,l]]/n;
	(* 1st component of elements of 'in' is label,
	   2nd component of elements of 'in' is centered coordinate *)
	in=Map[{#[[1]],#[[2]]-origin}&,l];
	(* 3rd component of elements of 'in' is polar angle *)
	in = Map[Append[#,PolarAngle[#[[2]]]]&,
		 in];
	sorted = Sort[in,
	              Function[{p1,p2},
		                p1[[3]]<p2[[3]] ||
		               (p1[[3]]==p2[[3]] && 
		               (p1[[2,1]]^2 + p1[[2,2]]^2 <  
		                p2[[2,1]]^2 + p2[[2,2]]^2))
		      ] (* end Function *)
		 ]; (* end Sort *)
        Map[{#[[1]],#[[2]]+origin}&,sorted]
	] /;		Apply[And, Map[NumberQ,Flatten[N[l]] ]]


(*============================== PolarAngle ===============================*)

PolarAngle[{x_,y_}]:=
	If[ y>=0,
		N[Arg[x + I y]],
                N[ 2 Pi + Arg[x + I y]]
        ] /; NumberQ[N[x]] && NumberQ[N[y]]


(*============================== CommonTangents ===========================*)

(* CommonTangents[set,hullL,hullR,rmL,lmR] *)
(* Returns {{l1,r1},{l2,r2}}.  l1 indicates the position in hullL of the
   leftmost point of the lower common tangent (of hullL & hullR).  r1 indicates
   the position in hullR of the rightmost point of the lower common tangent.
   l2 indicates the position in hullL of the leftmost point of the upper common
   tangent.  r2 indicates the position in hullR of the rightmost point of the
   upper common tangent. *)
(* Below 'rightmostLeft' is a pointer to the rightmost point in hullLeft,
   while 'leftmostRight' is a pointer to the leftmost point in hullRight. *)


CommonTangents[set:{{_,_}..}, hullLeft:{_Integer..}, hullRight:{_Integer..},
	       rightmostLeft_Integer, leftmostRight_Integer] :=
  Module[{x, y, predx, succy, succx, predy, lowertan={}, uppertan={}, 
	  lengthLeft = Length[hullLeft], lengthRight = Length[hullRight],
	  soaR, soaL},
    
    (*  LOWER TAN *)
    {x,y} = {rightmostLeft,leftmostRight};
    {predx,succy} = {pred[x,lengthLeft],succ[y,lengthRight]};
    While[ SameQ[lowertan,{}],
      soaR = SignOfArea[ set[[hullLeft[[x]]]], set[[hullRight[[y]]]],
			set[[hullRight[[succy]]]] ];
      If[ soaR <= 0,
	  If[ soaR == 0, Return[comtan[set,hullLeft,hullRight]] ];
          (* update y (by moving ccw on hullRight to succy) *)
          {succy,y} = {succ[succy,lengthRight],succy}, 
	  soaL = SignOfArea[ set[[hullLeft[[predx]]]], set[[hullLeft[[x]]]],
			    set[[hullRight[[y]]]] ];
          If[ soaL <= 0,
		If[ soaL == 0, Return[comtan[set,hullLeft,hullRight]] ];
                (* update x (by moving cw on hullLeft to predx) *)
                {predx,x} = {pred[predx,lengthLeft],predx},
		(* otherwise {x,y} points to the lower common tangent *)
		lowertan = {x,y}
          ]
      ]
    ];

    (* UPPER TAN *)
    {x,y} = {rightmostLeft,leftmostRight};
    {succx,predy} = {succ[x,lengthLeft],pred[y,lengthRight]};
    While[ SameQ[uppertan,{}], 
      soaR = SignOfArea[ set[[hullLeft[[x]]]], set[[hullRight[[y]]]],
			 set[[hullRight[[predy]]]] ];
      If[ soaR >= 0,
	  If[ soaR == 0, Return[comtan[set,hullLeft,hullRight]] ];
          (* update y (by moving cw on hullRight to predy) *)
          {predy,y} = {pred[predy,lengthRight],predy},
	  soaL = SignOfArea[ set[[hullLeft[[succx]]]], set[[hullLeft[[x]]]],
			     set[[hullRight[[y]]]] ];
          If[ soaL >= 0,
		If[ soaL == 0, Return[comtan[set,hullLeft,hullRight]] ];
                (* update x (by moving ccw on hullLeft to succx) *)
                {succx,x} = {succ[succx,lengthLeft],succx},
		(* otherwise {x,y} points to the upper common tangent *)
		uppertan = {x,y}
          ]
      ]
    ];

    {lowertan,uppertan}
  ] 


  (*** comtan is inefficient, but it will not be called unless points are  
     	colinear.  This inefficiency is not necessary. *)
comtan[set:{{_,_}..}, hullLeft:{_Integer..}, hullRight:{_Integer..}] :=
   Module[{hull = ConvexHull[ Join[ set[[hullLeft]], set[[hullRight]] ] ],
	   leng1 = Length[hullLeft], leng2 = Length[hullRight],
	   leng, uppertan, lowertan},
      leng = Length[hull];
      While[ !(hull[[1]] < leng1+1 && hull[[leng]] > leng1),
		hull = RotateLeft[hull]];
      (* Now the first points in hull belong to hullLeft. *)
      uppertan = { hull[[1]], hull[[leng]] - leng1 }; 
      While[ !(hull[[1]] > leng1 && hull[[leng]] < leng1+1),
		hull = RotateLeft[hull]];
      (* Now the first points in hull belong to hullRight. *)
      lowertan = { hull[[leng]], hull[[1]] - leng1 };
      {lowertan,uppertan}
   ]
	
	



(*======================== DelaunayTriangulation ===========================*)

Options[DelaunayTriangulation] = {Hull->False}

DelaunayTriangulation::colin = "Points `` are colinear."

DelaunayTriangulation[set:{{_,_}..},opts___Rule]:=
   Module[{orig=N[set], distinct = Distinct[N[set]], sorted, label, unlabeled,
	   delaunay, delval, hull, leftmost, rightmost,
	   returnhull = Hull /. {opts} /. Options[DelaunayTriangulation]},
     (
	 (* Make sure numerical coordinates are distinct and labeled
	    according to position in original list.  Also, order 
	    according to x coordinate before recursive triangulation. *)
	 sorted = Sort[distinct,
			 (#1[[2,1]]<#2[[2,1]] ||
			  (#1[[2,1]]==#2[[2,1]] && #1[[2,2]]>#2[[2,2]]))&
                  ];
         (* Save labels, but use points sans labels in calculating
	      triangulation. *)
         {label,unlabeled} = Transpose[sorted];
	 (* Recursive triangulation. *)
	 delaunay = Del[unlabeled];
	 If[ SameQ[delaunay,$Failed],
		Return[$Failed],
	        {delval,hull,leftmost,rightmost} = delaunay
	 ];
	 (* Add vertex label field to delval and relabel delval with positions
	    of points in original set. *)
         delval = Map[Module[{v=#[[1]],adjlist=#[[2]]},
	                	{label[[v]],
		                 Map[label[[#]]&,adjlist]}
			    ]&,
                      Transpose[{Range[Length[delval]],delval}]
                  ];
	 (* Sort vertex adjacency list according to vertex label. *)
	 delval = Sort[delval,(#1[[1]] < #2[[1]])&];
         If[returnhull,
	    (* Relabel convex hull with positions of points in original set. *)
	    hull = Map[label[[#]]&,hull];
	    {delval,hull},
	    delval
         ]
     ) /; Length[distinct] >= 3
   ] /;	Apply[And,Map[NumberQ,N[Flatten[set]]]] 






(************* 2 point triangulation *************)
(***** returns {delval,convexhull,leftmost,rightmost} *****)
Del[s:{{x1_,y1_},{x2_,y2_}}] :=
  Module[{hull},
     (* The point that has the smallest ordinate out of all rightmost points
	is listed first. *)
     hull = If[ x1 > x2 || (x1 == x2 && y1 < y2),
		{1,2},
		{2,1}];
     {{{2},{1}},hull,2,1}
  ]



(************* 3 point triangulation (non-colinear) *************)
(***** returns {delval,convexhull,leftmost,rightmost} *****)
Del[s:{{_,_},{_,_},{_,_}}] :=
  Module[{sorted,hull,maxx,miny,maxxlist,val,leftmost},
     sorted = PolarSort[Transpose[{Range[3],s}]];
     {hull,sorted} = Transpose[sorted];
     maxx = Apply[Max, Map[#[[1]]&,sorted]];
     maxxlist = Select[sorted, (#[[1]]==maxx)&];
     miny = Apply[Min, Map[#[[2]]&,maxxlist]];
     (* Rotate hull so that rightmost point having the smallest ordinate
		is listed first. *)
     While[ s[[ First[hull] ]] != {maxx,miny} ,
		hull = RotateLeft[hull]
     ];
     val = Map[Module[{p = Position[hull,#][[1,1]]},
		    Drop[ RotateLeft[hull,p], -1]
                  ]&,
	          Range[3]];
     leftmost = If[ s[[hull[[2]],1]] < s[[hull[[3]],1]] ||
		    	(s[[hull[[2]],1]] == s[[hull[[3]],1]] &&
			 s[[hull[[2]],2]] > s[[hull[[3]],2]]),
			2,
			3];
     {val,hull,leftmost,1}
  ]



(**************** >=6 point triangulation *****************)
(***** returns {delval,convexhull,leftmost,rightmost} *****)
Del[s:{{_,_}..}]:=
   Module[{leng0 = Length[s], leng1, leng2, val, hull, lenghull, order,
	leftptr, rightptr, xcoord, minx, maxx, lowertan, uppertan,
	s1, val1, hull1, l1, r1, lenghull1, num1, part1, lpart1,
	s2, val2, hull2, l2, r2, lenghull2, num2, part2, lpart2,
	LTl, LTr, UTl, UTr, colinearq1 = False, colinearq2 = False},
       If[ ColinearQ[s],
		Message[DelaunayTriangulation::colin,s];
		Return[$Failed]
       ];
       leng1 = Ceiling[leng0/2];
     (** s[[i]] is mapped into s1[[i]] for i=1,...,leng1 *)
       s1 = Take[s,leng1]; 
       leng2 = leng0 - leng1;
     (** s[[i]] is mapped into s2[[i-leng1]] for i=leng1+1,...,leng0 *)
       s2 = Drop[s,leng1];
     (** Delaunay triangulation of the subsets s1 & s2. *)
     (** If either subset is colinear, need information about overall hull,
		in order to compute a suitable triangulation for the colinear
		subset. *)
       colinearq1 = leng1 > 2 && ColinearQ[s1];
       colinearq2 = leng2 > 2 && ColinearQ[s2];
     (** Computing ConvexHull is time consuming; presumably there will be few
		cases where (colinearq1 || colinearq2) == True. *)
       If[ colinearq1 || colinearq2,
	(************ "Time-consuming" procedure for colinear points. ******)
		hull = ConvexHull[s];
		lenghull = Length[hull];
		While[ !(hull[[1]] < leng1+1 && hull[[lenghull]] > leng1),
			hull = RotateLeft[hull]
		];
	      (** Now points in s1 listed first in hull. *)
		lpart1 = Length[Select[hull,(# < leng1+1)&]];
		uppertan = {hull[[1]],hull[[lenghull]]};
		While[ !(hull[[1]] > leng1 && hull[[lenghull]] < leng1+1),
			hull = RotateLeft[hull]
		];
	      (** Now points in s2 listed first in hull. *)
		lpart2 = Length[Select[hull,(# > leng1)&]];
		lowertan = {hull[[lenghull]],hull[[1]]};
	      (** Calculate ptrs to leftmost & rightmost points in hull. *)
		xcoord = Map[#[[1]]&, s[[hull]] ];
		{minx,maxx} = {Min[xcoord],Max[xcoord]};
		leftptr = Position[xcoord,minx][[1,1]];
		rightptr = Position[xcoord,maxx][[1,1]];
	      (** Calculate subset adjacency lists such that merge will
		  work properly as it adds edges from lowertan to uppertan. *)
		val1 = If[ colinearq1,
			  If[ lpart1 == 1 ||
				(lpart1 > 1 && lowertan[[1]] != 1),
				  Join[{{2}},
				      Map[{#+1,#-1}&,Range[2,leng1-1]],
				       {{leng1-1}}],
				  Join[{{2}},
				      Map[{#-1,#+1}&,Range[2,leng1-1]],
				       {{leng1-1}}]
			  ],
			  Del[s1][[1]]
		       ];
                val2 = If[ colinearq2,
			  If[ lpart2 == 1 ||
				(lpart2 > 1 && lowertan[[2]] == leng1+1),
				  Join[{{leng1+2}},
				       Map[{#+1,#-1}&,
					 Range[leng1+2,leng1+leng2-1]],
				       {{leng1+leng2-1}}],
				  Join[{{leng1+2}},
				       Map[{#-1,#+1}&,
					 Range[leng1+2,leng1+leng2-1]],
				       {{leng1+leng2-1}}]
			  ],
			  Del[s2][[1]] + leng1
		      ],
        (************ "Normal" procedure for non-colinear points. ******)
       		{val1,hull1,l1,r1} = Del[s1];
       		{val2,hull2,l2,r2} = Del[s2];
     	      (** Offset labels in val2 & hull2 by number of points in s1. *)
       		val2+=leng1; hull2+=leng1;  
     	      (** Compute ptrs to lower & upper tangents of the two hulls. *)
       		lenghull1 = Length[hull1];
		lenghull2 = Length[hull2];
       		{{LTl,LTr},{UTl,UTr}} = CommonTangents[s,hull1,hull2,r1,l2];
       		lowertan = {hull1[[LTl]],hull2[[LTr]]};
       		uppertan = {hull1[[UTl]],hull2[[UTr]]};
	      (** num1 is # of points in hull1 also in hull *)
		num1 = Mod[LTl-UTl,lenghull1]+1;
	      (** part1 is the part of hull1 also in hull *)
		part1 = Take[RotateLeft[hull1,UTl-1],num1];
		leftptr = Mod[l1-UTl,lenghull1]+1;
	      (** num2 is # of points in hull2 also in hull *)
		num2 = Mod[UTr-LTr,lenghull2]+1;
	      (** part2 is the part of hull2 also in hull *)
		part2 = Take[RotateLeft[hull2,LTr-1],num2];
		rightptr = Mod[r2-LTr,lenghull2]+1+num1;
		hull = Join[part1,part2]
       ];
     (** Merge vertex adjacency lists val1 & val2 into val.  *)
       val = MergeDel[s,
		      Join[val1,val2],
		      lowertan, 		(* lowertan points *)
		      uppertan  		(* uppertan points *)
             ];
     (** Return new vertex adjacency list, convex hull, and pointers to 
	 leftmost & rightmost points in hull. *)
       {val,hull,leftptr,rightptr}
   ]



MergeDel[s:{{_,_}..},val:{{_Integer..}..},
	 {LTl_Integer,LTr_Integer},{UTl_Integer,UTr_Integer}] :=
   Module[{edge = {LTl,LTr}, uppertan = {UTl,UTr}, l0 = LTl, r0 = LTr,
           merge = val, l0r0r1LeftTurn, r0l0l1RightTurn,
	   l1, l2, r1, r2, l1val, r1val,  
	   r0ptr, l0ptr, r1ptr, l1ptr, r2ptr, l2ptr, lengr0, lengl0},
       
       (* Lengths of r0's val & l0's val AFTER first edge is added. *)
       lengr0 = Length[merge[[r0]]] + 1;
       lengl0 = Length[merge[[l0]]] + 1;
       (* Insert lower tangent edge {l0,r0}. *)
       r0ptr = 1;                   (* r0ptr will point to r0 in l0's val *)
       l0ptr = lengr0;              (* l0ptr will point to l0 in r0's val *)
       (* insert r0 in l0's val at position r0ptr;
	  insert l0 in r0's val at position l0ptr *)
       {merge[[l0]],merge[[r0]]} = insert[merge[[{l0,r0}]],l0,r0,l0ptr,r0ptr];
     
     (* Merge the vertex adjacency lists by adding & deleting edges in val;
	work from lower tangent edge to upper tangent edge. *)
     While[ edge != uppertan,
       (* Initially {l0,r0,r1} form a left turn & {r0,l0,l1} form a right. *)
       l0r0r1LeftTurn = r0l0l1RightTurn = True;

       (********************* Derive points r1 and r2. ********************)
       r1ptr = pred[l0ptr,lengr0];     (* r1ptr will point to r1 in r0's val *)
       r1 = merge[[r0,r1ptr]];
       If[ SignOfArea[s[[l0]],s[[r0]],s[[r1]]] == 1,  (* left turn *)
	  r2ptr = pred[r1ptr,lengr0]; (* r2ptr will point to r2 in r0's val *)
	  r2 = merge[[r0,r2ptr]];
	  (* Delete {r0,r1} only if edge {r1,r2} still exists. *)
	  While[ MemberQ[merge[[r1]],r2] &&
		 CircleMemberQ[ s[[l0]], s[[{r0,r1,r2}]] ],
            (* Delete edge {r0,r1}. *)
	    (* r1ptr is position of r1 in r0's val. *)
	    {merge[[r0]],merge[[r1]]} = delete[merge[[{r0,r1}]],
					Position[merge[[r1]],r0][[1,1]],r1ptr];
	    If[ l0ptr > r1ptr, l0ptr--];	(* update l0ptr *)
	    If[ r2ptr > r1ptr, r2ptr--];	(* update r2ptr *)
	    lengr0--;				(* update lengr0 *)
	    r1ptr = r2ptr; r2ptr = pred[r1ptr,lengr0]; 
	    r1 = merge[[r0,r1ptr]]; r2 = merge[[r0,r2ptr]]
          ],
	  l0r0r1LeftTurn = False
       ];

       (********************* Derive points l1 & l2. **********************)
       l1ptr = succ[r0ptr,lengl0]; 	(* l1ptr will point to l1 in l0's val *)
       l1 = merge[[l0,l1ptr]];
       If[ SignOfArea[s[[r0]],s[[l0]],s[[l1]]] == -1,  (* right turn *)
	  l2ptr = succ[l1ptr,lengl0]; (* l2ptr will point to l2 in l0's val *)
	  l2 = merge[[l0,l2ptr]];
	  (* Delete {l0,l1} only if edge {l1,l2} still exists. *)
	  While[ MemberQ[merge[[l1]],l2] &&
		 CircleMemberQ[ s[[r0]], s[[{l0,l1,l2}]] ],
	    (* Delete edge {l0,l1}. *)
	    {merge[[l0]],merge[[l1]]} = delete[merge[[{l0,l1}]],
					Position[merge[[l1]],l0][[1,1]],l1ptr];
	    If[ r0ptr > l1ptr, r0ptr--];	(* update r0ptr *)
	    If[ l2ptr > l1ptr, l2ptr--];	(* update l2ptr *)
	    lengl0--;				(* update lengl0 *)
	    l1ptr = l2ptr; l2ptr = succ[l1ptr,lengl0];
	    l1 = merge[[l0,l1ptr]]; l2 = merge[[l0,l2ptr]]
          ],
	  r0l0l1RightTurn = False
       ];



       (*********** Derive new l0 or new r0, & update l0ptr & r0ptr. **********)
       If[ !l0r0r1LeftTurn,
	  (* Connect r0 to l1. *)
	  l1val = merge[[l1]]; lengr0++; lengl0 = Length[l1val]+1;
	  {l0,r0ptr} = {merge[[l0,l1ptr]],
		        succ[Position[l1val,l0][[1,1]],lengl0]},
	  If[ !r0l0l1RightTurn,
	     (* Connect l0 to r1. *)
	     r1val = merge[[r1]]; lengl0++; lengr0 = Length[r1val]+1;
	     {r0,l0ptr,r0ptr} = {merge[[r0,r1ptr]],
				 Position[r1val,r0][[1,1]],
				 succ[r0ptr,lengl0]}, 
	     (* Neither {l0,r0,r1} nor {r0,l0,l1} are colinear. *)
	     If[ CircleMemberQ[ s[[l1]], s[[{l0,r0,r1}]] ],
		(* Connect r0 to l1. *) 
		l1val = merge[[l1]]; lengr0++; lengl0 = Length[l1val]+1;
		{l0,r0ptr} = {merge[[l0,l1ptr]],
			      succ[Position[l1val,l0][[1,1]],lengl0]},
		(* Connect l0 to r1. *)
		r1val = merge[[r1]]; lengl0++; lengr0 = Length[r1val]+1;
		{r0,l0ptr,r0ptr} = {merge[[r0,r1ptr]],
				    Position[r1val,r0][[1,1]],
				    succ[r0ptr,lengl0]},
             ]
          ]
       ];



       (***************** Update edge & insert into merge. ******************)
       edge = {l0,r0};
       (* Insert edge {l0,r0} by inserting r0 in l0's val at position r0ptr
	  and inserting l0 in r0's val at position l0ptr. *)
       {merge[[l0]],merge[[r0]]} =
	  insert[ merge[[{l0,r0}]],l0,r0,l0ptr,r0ptr ]
     ];
     merge
   ]


        

(*================= Auxiliery Triangulation Functions ======================*)

succ[i_Integer,mod_Integer] := Mod[i,mod]+1

succ[{{i_Integer}},mod_Integer] := Mod[i,mod]+1

pred[i_Integer,mod_Integer] := Mod[i-2,mod]+1

pred[{{i_Integer}},mod_Integer] := Mod[i-2,mod]+1


   
(**** Insert v2 in position v1's adjacency list at position v2ptr. 
      Insert v1 in v2's adjacency list at position v1ptr. *)
(**** Returns object of the form {val1,val2}. *)
insert[{val1:{_Integer..},val2:{_Integer..}},
	v1_Integer, v2_Integer, v1ptr_Integer, v2ptr_Integer] :=
    {Insert[val1, v2, v2ptr ],
     Insert[val2, v1, v1ptr ]
    } 


(**** Delete v2 (pointed to by v2ptr) from v1's adjacency list val1.
      Delete v1 (pointed to by v1ptr) from v2's adjacency list val2. *)   
delete[{val1:{_Integer..},val2:{_Integer..}},
	v1ptr_Integer, v2ptr_Integer] :=
   {Delete[val1,v2ptr],
    Delete[val2,v1ptr]
   } 
   

(**** The default is to not include membership in circle boundary. *)

CircleMemberQ[q:{_,_},circle:{{_,py1_},{_,py2_},{_,py3_}},
		boundary_:False] :=
	CircleMemberQ[{q},circle,boundary] /; !(py1 == py2 == py3)

CircleMemberQ[q:{{_,_}..},{p1:{_,py1_},p2:{_,py2_},p3:{_,py3_}},
		boundary_:False] :=
   Module[{x, y, u1, u2, u3, m2, m3, bx2, by2, bx3, by3, r2},
    (* Insure that there is no "/0" problem. *)
    {u1,u2,u3} = If[SameQ[py1,py2],
                    {p3,p2,p1},
                    If[SameQ[py1,py3],
                       {p2,p1,p3},
                       {p1,p2,p3}]];
    m2 = -Apply[Divide,u1-u2]; {bx2,by2} = (u1+u2)/2;
    m3 = -Apply[Divide,u1-u3]; {bx3,by3} = (u1+u3)/2;
    (* Center is intersection of bisectors of two pairs of points. *)
    x = (by2-by3+m3 bx3-m2 bx2)/(m3-m2);
    y = ((m3 m2)(bx3-bx2) + m3 by2 - m2 by3)/(m3-m2);
    r2 = (u1[[1]] - x)^2 + (u1[[2]] - y)^2;
    (* Result is True if one or more points listed in q belong to the circle. *)
    If[ boundary,
	Apply[Or,Map[
		( Chop[(#[[1]]-x)^2+(#[[2]]-y)^2 - r2] <= 0 )&,
		 q]],
        Apply[Or,Map[
      		( Chop[(#[[1]]-x)^2+(#[[2]]-y)^2 - r2] < 0 )&,
      		 q]]
    ]
  ] /; !(py1 == py2 == py3)



(*======================= DelaunayTriangulationQ ===========================*)

(*** It is assumed that the val is such that a vertex on the hull lists first
     the neighbor that follows it on a counterclockwise traversal of the convex
     hull. *)

DelaunayTriangulationQ::inval = "Triangle `` is not a valid Delaunay triangle."


(*** Only point set and Delaunay vertex adjacency list available. ***)
DelaunayTriangulationQ[set:{{_,_}..},val:{{_Integer,{_Integer..}}..}] :=
  Module[{orig =  N[set], maxx, maxxlist, miny, start, hull},
    (* "convexhull" is found by traversing the vertex adjacency list under
       the assumption that this is faster than calculating it from scratch. *)
    (* Find one point on hull, let's say rightmost point having the smallest
		ordinate. *)
     maxx = Apply[Max, Map[#[[1]]&,orig]];
     maxxlist = Select[orig, (#[[1]]==maxx)&];
     miny = Apply[Min, Map[#[[2]]&,maxxlist]];
     start = Scan[(If[orig[[#[[1]]]]=={maxx,miny},
		     Return[#[[1]]]])&,
                  val];
     hull = TriangulationToHull[val,start];
     If[ SameQ[hull,$Failed],
		False,
     		DelaunayTriangulationQ[set,val,hull]
     ]
   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&  (* numerical points *)
	(Max[Flatten[val]] <= Length[set]) &&  (* val refers to valid pts *)
	Length[val] >= 3 &&		(* val must have at least 1 triangle *)
	!ColinearQ[set]			(* set must not be colinear *)


(*** Only point set and vertex adjacency list (val) available. ***)
DelaunayTriangulationQ[set:{{_,_}..},val:{{_Integer,{_Integer..}}..},
	               hull:{_Integer..}]:=
   Module[{orig =  N[set], circle, vertices = Map[#[[1]]&,val], delaunayq},
     (* Compute list of triples defining Delaunay triangles. *)
     circle = Map[
		 (Module[{vert1=#[[1]],vertlist=#[[2]],
			  vertNum=Length[#[[2]]],circle},
       (* Create list of triples representing triangles adjacent to vert1. *)
                     circle = Map[{vert1, vertlist[[#]],
			           vertlist[[Mod[#,vertNum]+1]]}&,
		                   Range[vertNum]];
                     If[MemberQ[hull,vert1],
			    circle = Drop[circle,-1]];
                     Map[Sort,circle]
                 ])&,
		 val];
    (* Eliminate duplicates. *)
    circle = Union[Flatten[circle,1]];
    delaunayq = If[ Apply[Or,Map[ (Apply[SignOfArea,orig[[#]]] == 0)&, circle]],
		    (* One or more of the triangles are colinear! *)
		    False,	
		    (* Check that the circumcircles do not include
		       other vertices *)
                    Scan[(Module[{tri = #, compare, x},
		  	      compare = Complement[
					 Apply[Union,
					  Map[Module[{x = #},
					       Select[val,
					         MatchQ[#,{x,_List}]&,
					         1][[1,2]]
					      ]&,
					      tri]],
					 tri];	
       			      If[ CircleMemberQ[orig[[compare]],orig[[tri]]], 
				   Message[DelaunayTriangulationQ::inval,tri];
				   Return[False]]
          	         ])&,
	  	         circle]
    ];
    If[ SameQ[delaunayq,Null],
	True,
	False
    ]
   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&	(* numerical points *)
	(Max[Flatten[val]] <= Length[set]) &&  (* val refers to valid pts *)
	(Max[hull] <= Length[set]) &&	  (* hull refers to valid pts *)
	Length[val] >= 3 &&		(* val includes at least 1 triangle *)
	!ColinearQ[set]			(* set must not be colinear *)



(*====================== TriangulationToHull ==============================*)

(* It is assumed that the val is such that for a vertex on the hull, the first
   neighbor listed is that one lying next on a counterclockwise traversal of
   the hull. *)
TriangulationToHull[val:{{_Integer,{_Integer..}}..},start_Integer]:=
  Module[{label1, relabel1, start1, val1, hull, next, hullLength},
    (* Relabel vertex labels to make tracing hull easier.  This wouldn't
       be necessary if all the points in the original set were unique. *)
     label1 = Map[#[[1]]&,val];
     relabel1 = MapIndexed[(#1->#2[[1]])&,label1];
     start1 = start /. relabel1;
     val1 = Map[(#[[2]] /. relabel1)&,val];
     hull = {start};
     hullLength = 1;
     next = val1[[start,1]];
    (* Trace points on hull using relabeled adjacency list. *)  
     While[next != start,
       (* Guard against infinite loops. *)
       If[ hullLength > Length[val],
		Return[$Failed]];
       AppendTo[hull,next];
       hullLength++;
       next = val1[[next,1]]
     ];
    (* Label convexhull with original labels. *)
     label1[[hull]]
  ]


(*=========================== VoronoiDiagram ================================*)

(*** Only point set available. ***)
VoronoiDiagram[set:{{_,_}..}]:=
   Module[{delaunay=DelaunayTriangulation[set,Hull->True]},
     If[ SameQ[delaunay,$Failed],
	$Failed,
        VoronoiDiagram[set,delaunay[[1]],delaunay[[2]]]
     ]		/; FreeQ[delaunay,DelaunayTriangulation]
   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]]  (* numerical points *)


(*** Both point set and Delaunay vertex adjacency list available. ***)
VoronoiDiagram[set:{{_,_}..},delval:{{_Integer,{_Integer..}}..}]:=
   Module[{hull = Module[{orig = N[set], maxx, maxxlist, miny, start},
                      (* "convexhull" is found by traversing the Delaunay
			 vertex adjacency list under the assumption that
                         this is faster than calculating it from scratch. *)
                      (* Find one point on hull, let's say rightmost point
			 having the smallest ordinate. *)
                      maxx = Apply[Max, Map[#[[1]]&,orig]];
                      maxxlist = Select[orig, (#[[1]]==maxx)&];
                      miny = Apply[Min, Map[#[[2]]&,maxxlist]];
                      start = Scan[(If[orig[[#[[1]]]]=={maxx,miny},
		     			Return[#[[1]]]])&,
                                   delval];
                      TriangulationToHull[delval,start]
                 ]},
     If[ SameQ[hull,$Failed] || !FreeQ[hull,TriangulationToHull] ,
		$Failed,
     		VoronoiDiagram[set,delval,hull]
     ]
   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&  (* numerical points *)
	(Max[Flatten[delval]] <= Length[set]) &&  (* val refers to valid pts *)
	Length[delval] >= 3		(* val must have at least 1 triangle *)


(*** Point set, Delaunay vertex adjacency list, and convexhull available. ***)
VoronoiDiagram[set:{{_,_}..},delval:{{_Integer,{_Integer..}}..},
	       hull:{_Integer..}]:=
   Module[{orig =  N[set], lhull = Length[hull], ltrue, relabel,
	   vorval, vorvertices, truevertices, quasivertices,
	   truecoordinates, quasicoordinates},
    (* Compute Voronoi vertex adjacency list. *)
     vorval = Map[
		 (Module[{delvert1=#[[1]],delvertlist=#[[2]],
			  delvertNum=Length[#[[2]]],vorpoly},
       (*  Each vertex of a Voronoi polygon corresponds to a Delaunay vertex
	   triple.  Exceptions: a "quasi" Voronoi vertex corresponds to a
	   Delaunay vertex pair; a "true"  Voronoi vertex can correspond to
	   two Delaunay vertex triples. *)
       (* Create list of triples representing triangles adjacent to delvert1.
	  These determine Voronoi vertices. *)
                     vorpoly = Map[{delvert1, delvertlist[[#]],
			            delvertlist[[Mod[#,delvertNum]+1]]}&,
		                   Range[delvertNum]];
       (* If delvert1 is on hull, add information needed for determining
	  infinite rays of associated Voronoi polygon. *)	
                     vorpoly = If[MemberQ[hull,delvert1],
				  Module[{truevert=Map[Sort,Drop[vorpoly,-1]]},
				    Join[truevert,
				         { {Sort[{vorpoly[[delvertNum,1]],
					          vorpoly[[delvertNum,2]]}],
					    truevert[[delvertNum-1]]},
					   {Sort[{vorpoly[[delvertNum,1]],
					     vorpoly[[delvertNum,3]]}],
					    truevert[[1]]}
					 } ]
				  ],
                                 Map[Sort,vorpoly]
			       ]
                 ])&,
		 delval];

    vorvertices = Union[Flatten[vorval,1]];
    quasivertices = Take[vorvertices,lhull];
    truevertices = Drop[vorvertices,lhull];
    coordinates = Map[Apply[circlecenter,orig[[#]]]&,truevertices];

    (* True Voronoi vertices. *)
    truecoordinates = Union[coordinates];
    ltrue = Length[truecoordinates];
    (* lables for true Voronoi vertices *)
    relabel = Flatten[
		MapIndexed[Module[{pos = Flatten[Position[coordinates,#1],1],
				 label = #2[[1]]},
				 Map[(truevertices[[#]] -> label)&,pos]
                           ]&,
			   truecoordinates],
              1];
    (* lables for quasi Voronoi vertices *)
    relabel = Join[relabel,
		   MapIndexed[(#1 -> #2[[1]]+ltrue)&,
			      quasivertices
                   ]
              ];
    (* Label Voronoi vertices in Voronoi val uniquely, and eliminate
	duplicates. *)
    vorval = Map[Module[{val = # /. relabel, unique, current},
		   current = val[[1]];
		   unique = {current};
	           Scan[(If[ # != current,
			     current = #;
			     AppendTo[unique,#]
                         ])&,
			Drop[val,1]];
                   unique
                 ]&,
		 vorval];
    (* Add a field indicating Delaunay vertex label to each Voronoi val. *)
    vorval = MapIndexed[{delval[[#2[[1]],1]],#1}&,vorval];

    (* Quasi Voronoi vertices. *)
    (* Re-establish the ccw direction of the hull edges that define quasi
	Voronoi vertices.  (Direction was lost when unique vertices were
	found using "Union".)  Then determine ray "tail" from ray "head"
	and the associated Delaunay hull edge. *)
    quasicoordinates = Map[Module[{edge = #[[1]],
				   truevertex = #[[2]] /. relabel,
				   head, ptr1},
			   ptr1 = Position[hull,edge[[1]]][[1,1]];
			   If[edge[[2]] != hull[[ succ[ptr1,lhull] ]],
					edge = edge[[{2,1}]] ];
                           head = truecoordinates[[ truevertex ]];
			   Ray[head, raytail[ head, orig[[edge]] ]]
                        ]&,
			quasivertices
		    ];

    {Join[truecoordinates,quasicoordinates],vorval}   
   ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&	(* numerical points *)
	(Max[Flatten[delval]] <= Length[set]) &&  (* val refers to valid pts *)
	(Max[hull] <= Length[set]) &&	  (* hull refers to valid pts *)
	Length[delval] >= 3		(* val includes at least 1 triangle *)


(*** Compute the center of the circumcircle of triangle {p1,p2,p3}. *)
circlecenter[p1:{_,py1_}, p2:{_,py2_}, p3:{_,py3_}] :=
   Module[{u1, u2, u3, m2, m3, bx2, by2, bx3, by3},
     (* Insure that there is no "/0" problem. *)
     {u1,u2,u3} = If[SameQ[py1,py2],
			{p3,p2,p1},
			If[SameQ[py1,py3],
				{p2,p1,p3},
				{p1,p2,p3}]];
     m2 = -Apply[Divide,u1-u2]; {bx2,by2} = (u1+u2)/2;
     m3 = -Apply[Divide,u1-u3]; {bx3,by3} = (u1+u3)/2;
     { (by3-by2) + m2 bx2 - m3 bx3,
       m2 m3 (bx2-bx3) + m2 by3 - m3 by2} / (m2 - m3)
   ] /; !(py1 == py2 == py3)


(*** Compute the tail of the ray having head "head" and bisecting segment
     {e1,e2}.  The tail is defined by a point a distance dist(e1,e2) from
     the ray head in a direction such that {e1,head,tail} is a right turn
     OR a distance dist(e1,e2) from the edge midpoint {xm,ym} in a direction
     such that {e1,{xm,ym},tail} is a right turn.  The distance is computed
     from the head or from {xm,ym} depending on which will look better in
     DiagramPlot. *)
raytail[head:{xh_,yh_},{e1:{_,_},e2:{_,_}}] :=
   Module[{m, out, xm, ym, distfromhead, soa, 
	   d0 = Sqrt[ Apply[Plus,(e1-e2)^2] ]},
      {xm,ym} = (e1+e2)/2; (* edge midpoint *)
      (* Determine whether d0 will be measured from head or {xm,ym}; 
	 this is an aesthetic consideration. *)
      If[ SignOfArea[e1,{xm,ym},head] == -1,
		distfromhead = True,
		distfromhead = False;
		d0 = d0 + Sqrt[ Apply[Plus,(head - {xm,ym})^2] ]
      ];
      out = If[ e1[[2]] === e2[[2]],
	        (* ray is vertical *)
	        {{xh,yh-d0},{xh,yh+d0}},
	        (* ray has slope m *)
		m = -Apply[Divide,e1-e2];
		Map[{#,m(#-xm)+ym}&,
	            Module[{b = (m(ym-m xm-yh)-xh)/(1+m^2),sqrt},
		      sqrt = Sqrt[b^2 - (xh^2 + (ym-m xm-yh)^2 - d0^2)/(1+m^2)];
		      -b + {sqrt,-sqrt}
                    ]
		] 
      ];
      (* Choose the correct direction of the ray.
	 {e1,head,tail} or {e1,{xm,ym},tail} should form a right turn. *)
      soa = If[ distfromhead,
		SignOfArea[e1,head,out[[1]]],
		SignOfArea[e1,{xm,ym},out[[1]]]
            ];
      If[ soa == -1,
	 	out[[1]], (* right turn *)
	 	out[[2]]] (* left or no turn *)
   ]




VoronoiDiagram::notimp =
  "VoronoiDiagram is not yet implemented for dimension ``."

VoronoiDiagram[set_List]:=
   Module[{out=Module[{d=Length[set[[1]]]},
		  Message[VoronoiDiagram::notimp,d];
		  $Failed]},
	out /; !SameQ[out,$Failed]
   ] /; Apply[And,Map[NumberQ,N[set]]] ||
	(Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
	 Apply[Equal,Map[Length[#]&,set]] &&
	 Apply[And,Map[SameQ[Head[#],List]&,set]])



(*=========================== NearestNeighbor ==============================*)

NearestNeighbor[input:{{_,_}..},vorpts_List,
		vorval:{{_Integer,{_Integer..}}..}]:=
  Module[{ninput = N[input], closed = {}, open = {}},
    (* If VoronoiDiagram produced vorval, the open polygons are guarranteed to
       be listed after the closed polygons.  The following general approach
       to identifying open and closed polygons is useful if some other method
       generated the diagram of convex polygons represented by vorval. *)
    Scan[(If[FreeQ[vorpts[[#[[2]]]],Ray],
	     AppendTo[closed,#],
	     AppendTo[open,#]
	  ])&,
	  vorval];
    (* If VoronoiDiagram produced vorval, each open polygon is guarranteed to
       list the two Ray "vertices" last.  The following scheme for putting
       open polygons in this format is useful if some other method generated
       the diagram of convex polygons represented by vorval.  A valid open
       polygon will list the two rays adjacently. *)
    open = Map[Module[{vertexlist = #[[2]], length=Length[#[[2]]]},
		   While[!(SameQ[Head[vorpts[[#[[2,length]]]]],Ray] &&
			   SameQ[Head[vorpts[[#[[2,length-1]]]]],Ray]),
			      vertexlist = RotateLeft[vertexlist]];
                   {#[[1]],vertexlist}
               ]&,
	       open];
    (* First compute an internal point for each closed convex polygon. *)
    closed = Map[Append[#,
			Apply[Plus, vorpts[[ Take[#[[2]],3] ]] ]/3
                 ]&,
                 closed];
    (* For each input point, determine the label of the nearest neighbor. *)
    Map[Module[{query = #, class},
	  (* Scan closed polygons for location of query. *)
	  class = Scan[Module[{internal = #[[3]]},
			 If[PolygonMemberQ[query-internal,
					   Map[(#-internal)&,
					   vorpts[[ #[[2]] ]] ]],
			    Return[#[[1]]]
                         ]
                       ]&,
                       closed];
          (* If query not located yet, scan the open polygons. *)
          If[SameQ[class,Null],
	    class = Scan[Module[{ray = vorpts[[ Take[#[[2]],-2] ]]},
                           If[
			      (* Include in the open polygon those points on
				 ray forming the boundary first encountered
				 when one traverses the open polygons in a
				 counterclockwise direction. *)
			      (* left turn or no turn *)
			      SignOfArea[query,ray[[1,1]],ray[[1,2]]] >= 0 &&
			      (* right turn *)
			      SignOfArea[query,ray[[2,1]],ray[[2,2]]] == -1,
				 Return[#[[1]]]
                           ]
                         ]&,
			 open]
          ];
	  If[SameQ[class,Null],$Failed,class]
	]&, 
	ninput]
  ] /; Apply[And,Map[NumberQ,N[Flatten[input]]]] &&
	ValidDiagramQ[vorpts,vorval]


NearestNeighbor[{a_,b_},vorpts_List,vorval:{{_Integer,{_Integer..}}..}]:=
  NearestNeighbor[{{a,b}},vorpts,vorval][[1]] /;
    (NumberQ[N[a]] && NumberQ[N[b]] && ValidDiagramQ[vorpts,vorval])

NearestNeighbor[input:{{_,_}..},pts:{{_,_}..}]:=
  Module[{voronoi,vorvertices,vorval},
    voronoi = VoronoiDiagram[pts];
    If[ SameQ[voronoi,$Failed],
	Return[$Failed]];
    {vorvertices,vorval} = voronoi;
    NearestNeighbor[input,vorvertices,vorval]
  ] /;  Apply[And,Map[NumberQ,N[Flatten[input]]]] &&
	Apply[And,Map[NumberQ,N[Flatten[pts]]]]

NearestNeighbor[{a_,b_},pts:{{_,_}..}]:=
  NearestNeighbor[{{a,b}},pts][[1]] /;
    NumberQ[N[a]] && NumberQ[N[b]] && Apply[And,Map[NumberQ,N[Flatten[pts]]]]



(*============================= PolygonMemberQ ===============================*)

(* Binary search of wedges that comprise a closed convex polygon containing
   the origin.  Returns True if the polygon contains the query. Polygon must
   contain origin because this makes use of PolarAngle. *)
PolygonMemberQ[query:{_,_},polygon:{{_,_}..}] :=
    Module[{median,a1,am,aq,p = Append[polygon,First[polygon]]},
      While[ Length[p]>2,
        median = Floor[Length[p]/2]+1;
	{a1,am,aq} = Map[PolarAngle,{p[[1]],p[[median]],query}];
	If[ (a1 <= aq < am) || ((a1 > am) && !(aq < a1 && aq >= am)),
	   (* query lies in wedge between p[[1]] & p[[median]] *)
	   p = Take[p,median],
	   (* query lies in wedge between p[[median]] & p[[1]] *)
	   p = Drop[p,median-1]
        ]
      ]; (* here p is of the form {{x1,y1},{x2,y2}} *)
      (* Include points on the boundary... this means that lower-labled
	 polygons will be favored over higher-labled polygons when it
	 comes to points on the boundary. *)
      (* test for left turn or no turn *)
      SignOfArea[query,p[[1]],p[[2]]] >= 0 
    ]


(*============================= ValidDiagramQ ===============================*)

ValidDiagramQ[vert_List,val:{{_Integer,{_Integer..}}..}] :=
	(* vertices must have head List or Ray *)
	Apply[And,Map[(SameQ[Head[#],List] || SameQ[Head[#],Ray])&,vert]] &&
	(* vertices must have length 2 *)
	Apply[And,Map[(Length[#]==2)&,vert]] &&
	(* each polygon must have either 0 rays (closed) or 2 rays (open) *)
     Apply[And,Map[Module[{raynum=Length[Position[ vert[[ #[[2]] ]], Ray]]},
		     raynum == 0 || raynum == 2
		   ]&,
		   val]] &&
        (* adjacency list should not refer to vertices not in vertex list *)
	(Length[vert] == Length[Union[Flatten[Map[#[[2]]&,val]]]])



(*=================== Auxiliery Plotting Functions ===========================*)

(* In computational geometry, it's important to plot things without any
	scaling or stretching. *)
absolutePlotRange[set:{{_,_}..}] :=
  Module[{xcoord, ycoord, minx, maxx, miny, maxy,
	  xhalfrange, yhalfrange, xmid, ymid},
      {xcoord,ycoord} = Transpose[set]; 
      {minx,maxx} = {Min[xcoord],Max[xcoord]};
      {miny,maxy} = {Min[ycoord],Max[ycoord]};
      {xhalfrange,yhalfrange} = {maxx-minx,maxy-miny}/2;
      If[ xhalfrange > yhalfrange,
	 ymid = miny+yhalfrange;
	 {{minx,maxx},{ymid-xhalfrange,ymid+xhalfrange}},
	 xmid = minx+xhalfrange;
	 {{xmid-yhalfrange,xmid+yhalfrange},{miny,maxy}}
      ]
 ]


(*============================= DiagramPlot ================================*)

Options[DiagramPlot] =  Join[{LabelPoints->True, TrimPoints->0},
			     Options[Graphics]]

DiagramPlot::trim =
"Warning: TrimPoints -> `` is not a nonnegative integer.  TrimPoints -> 0 used
instead."
			


(* Polygon diagram is already computed and polygon center coordinates are
   available. *)
DiagramPlot[set:{{_,_}..},vorvert_List,val:{{_Integer,{_Integer..}}..},
	opts___] :=
  Module[{labelpoints, trimpoints, orig = N[set], distinct},
    {labelpoints,trimpoints} = {LabelPoints,TrimPoints} /. {opts} /.
					Options[DiagramPlot];
    If[ !(IntegerQ[trimpoints] && !Negative[trimpoints]),
		Message[DiagramPlot::trim,trimpoints];
		trimpoints = 0
    ];
    distinct = Map[{#[[1]],orig[[#[[1]]]]}&,val];
    diagramplot[distinct,vorvert,val,labelpoints,trimpoints,opts] 
  ] /;  Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
	ValidDiagramQ[vorvert,val] &&
	(*  adjacency list should not refer to points not in set *)
       (Max[Map[#[[1]]&,val]] <= Length[set])


(* Compute Voronoi diagram before plotting. *)
DiagramPlot[set:{{_,_}..},opts___Rule] :=
  Module[{voronoi = VoronoiDiagram[set], distinct, vorvert, val},
    If[ SameQ[voronoi,$Failed], Return[$Failed]];
    {vorvert,val} = voronoi;
    DiagramPlot[set,vorvert,val,opts]
  ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] 


(* diagramplot *)
diagramplot[distinct:{{_,{_,_}}..},vorvert_List,
	val:{{_Integer,{_Integer..}}..},
	labelpoints_Symbol,trimpoints_Integer,opts___] :=
  Module[{centerlist,vertexlist={PointSize[.012]},
	  linelist={Thickness[.003]}, leng = Length[distinct],
	  trim = trimpoints, original, absolutepoints, centroid,
	  max2, dist2, pos, xmin, xmax, ymin, ymax, delx, dely},
    (* diagram polygon "centers" *)
    centerlist = If[labelpoints,
                    Map[Apply[Text,#]&,distinct],
     	            Prepend[Map[Point,Map[#[[2]]&,distinct]],
		      PointSize[.012]]
		 ];
    (* diagram vertices and infinite rays *)
    Scan[(If[FreeQ[#,Ray],
	    AppendTo[vertexlist,Point[#]],
	    AppendTo[linelist,Line[Apply[List,#]]],
	  ])&,
	 vorvert];
    (* closed diagram polygons *)
    Scan[(Module[{list = vorvert[[#[[2]]]], closedpoly},
	     closedpoly = FreeQ[list,Ray];
	     list = Select[list,(!SameQ[Head[#],Ray])&];
	     If[ Length[list] > 1, (* val contains at least one segment *)
		If[ closedpoly,
		  (* add in the segment that closes the poly *)
	          AppendTo[linelist,Line[Join[list,{list[[1]]}]]],
		  (* just add in explicit segments *)
		  AppendTo[linelist,Line[list]]
                ]
             ]
          ])&,
	 val];
    (* Compute points that will determine PlotRange. *)
    original = Map[#[[2]]&,distinct];
    If[ trim == 0,
	(* include all diagram vertices AND ray tails *)
	absolutepoints = Join[ original,
			       Map[(If[ SameQ[Head[#],Ray],
					#[[2]],
					#])&,
    				   vorvert]
                             ],
        (* exclude ray tails and (trim-1) of the furthest outlying diagram
	   vertices *)
        absolutepoints = Join[ original,
			       Select[vorvert, !SameQ[Head[#],Ray]&] ];
        If[ trim > 1,
        	centroid = Apply[Plus,original]/Length[original];
		max2 = Max[Map[Apply[Plus,(#-centroid)^2]&,
		               original
                          ]];
        	dist2 = Map[Apply[Plus,(#-centroid)^2]&,
			        Drop[absolutepoints,leng]
                        ];
        	pos = Position[dist2,Max[dist2]][[1,1]];
                (* Delete diagram vertices from 'absolutepoints' as long as 
		   they lie further from the centroid then does the furthest
		   polygon "center". *)
                 While[ trim!=1 && dist2[[pos]] > max2,
			dist2 = Delete[dist2,pos];
			absolutepoints = Delete[absolutepoints, pos+leng];
			pos = Position[dist2,Max[dist2]][[1,1]];
			trim--
                 ]
	]
    ];
    {{xmin,xmax},{ymin,ymax}} = absolutePlotRange[absolutepoints];
    {delx,dely} = .05/1.9 {xmax-xmin,ymax-ymin};
    optslist = Select[{opts},!(SameQ[#[[1]],LabelPoints] ||
			       SameQ[#[[1]],TrimPoints])&];
    Show[Graphics[
        {centerlist,
         vertexlist,
         linelist},
        PlotRange -> {{xmin-delx,xmax+delx},{ymin-dely,ymax+dely}},
	AspectRatio->1
    ],optslist]
  ]



(*=========================== PlanarGraphPlot ===============================*)

Options[PlanarGraphPlot] = Join[{LabelPoints->True},
				Options[Graphics]]


(* This is the form for plotting convex hull. *)
PlanarGraphPlot[set:{{_,_}..}, circuit:{_Integer..}, opts___]:=
   Module[{orig = N[set], pointlist, linelist, distinct,
	   xmax, xmin, ymax, ymin, delx, dely, optslist,
	   labelpoints = LabelPoints /. {opts} /. Options[PlanarGraphPlot]},
      If[labelpoints,
	 distinct = Distinct[orig];
	 pointlist = Map[Text[#,orig[[#]]]&,Transpose[distinct][[1]]],
	 pointlist = Prepend[Map[Point,orig],PointSize[.012]] ];
      linelist = {Thickness[.003],
		  Line[Append[ orig[[circuit]],orig[[circuit[[1]]]] ]]};
      {{xmin,xmax},{ymin,ymax}} = absolutePlotRange[orig];
      {delx,dely} = .05/1.9 {xmax-xmin,ymax-ymin};
      optslist = Select[{opts},!SameQ[#[[1]],LabelPoints]&];
      Show[Graphics[
	  {pointlist,
	   linelist},
	  PlotRange -> {{xmin-delx,xmax+delx},{ymin-dely,ymax+dely}},
	  AspectRatio->1
      ],optslist]
   ]	/;  (Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
	     (Max[circuit] <= Length[set]))


(* Vertex adjacency list is already computed. *)
PlanarGraphPlot[set:{{_,_}..},val:{{_Integer,{_Integer..}}..},opts___] :=
  Module[{orig = N[set], pointlist, linelist, distinct, pairs, 
	  xmax, xmin, ymax, ymin, delx, dely, optslist,
	  labelpoints = LabelPoints /. {opts} /. Options[PlanarGraphPlot]},
      If[labelpoints,
	 distinct = Distinct[orig];
	 pointlist = Map[Text[#,orig[[#]]]&,Transpose[distinct][[1]]],
	 pointlist = Prepend[Map[Point,orig],PointSize[.012]] ];
      pairs = Map[(Outer[List,{#[[1]]},#[[2]]][[1]])&,
	          val];
      pairs = Union[Map[Sort,Flatten[pairs,1]]];
      linelist = Prepend[Map[Line[orig[[#]]]&,pairs],Thickness[.003]];
      {{xmin,xmax},{ymin,ymax}} = absolutePlotRange[orig];
      {delx,dely} = .05/1.9 {xmax-xmin,ymax-ymin};
      optslist = Select[{opts},!SameQ[#[[1]],LabelPoints]&];
      Show[Graphics[
	  {pointlist,
	   linelist},
	  PlotRange -> {{xmin-delx,xmax+delx},{ymin-dely,ymax+dely}},
	  AspectRatio->1
      ],optslist]
  ]   /;   Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
			(Max[Flatten[val]] <= Length[set])


(* Compute Delaunay triangulation before plotting. *)
PlanarGraphPlot[set:{{_,_}..},opts___] :=
  Module[{delaunay = DelaunayTriangulation[set]},
    If[ SameQ[delaunay,$Failed],
	$Failed,
        PlanarGraphPlot[set,delaunay,opts]
    ]
  ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]]






(*======================= TriangularSurfacePlot =============================*)

Options[TriangularSurfacePlot] = Options[Graphics3D]


(* Vertex adjacency list is already computed. *)
TriangularSurfacePlot[set:{{_,_,_}..},val:{{_Integer,{_Integer..}}..},
	opts___] :=
  Module[{orig = N[set], distinct, label, hull, polygonlist},
    (* If two or more 3d points have the same {x,y} coordinates, but different
       z coordinates, only the first instance of the {x,y} coordinates in the
       original set is considered.  Other values for z are ignored. *)
     distinct = Map[{#[[1]],orig[[ #[[1]] ]]}&,val];
    (* In the case of Delaunay triangulation, convexhull may be extracted by
       traversing the vertex adjacency list val.  However, since val may
       represent some other triangulation, the convex hull is computed from
       scratch. *)
     label = Map[#[[1]]&,val];
     hull = Map[ label[[#]]&,
		       convexhull[Map[Drop[#[[2]],-1]&,distinct],True]
	     ];
    (* find all triangles surrounding interior points (i.e. not on hull) *)
    polygonlist = Map[Module[{vertex=#[[1]],adjvert=#[[2]],
			      vertNum=Length[#[[2]]]},
			Map[{vertex,adjvert[[#]],
			     adjvert[[ Mod[#,vertNum]+1 ]]}&,
			     Range[vertNum]]
                      ]&,
		      Select[val,!MemberQ[hull,#[[1]]]&]
		  ];
    (* Map Polygon onto coordinates of unique triangles  *)
    polygonlist = Map[Polygon[Module[{triangle = #},
			        Map[distinct[[#,2]]&,
			            triangle]
                              ]]&,
		      Union[Map[Sort,Flatten[polygonlist,1]]]
                  ];
    Show[Graphics3D[polygonlist],opts]
  ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]] &&
	(Max[Flatten[val]] <= Length[set])

(* Compute Delaunay triangulation before plotting. *)
TriangularSurfacePlot[set:{{_,_,_}..},opts___] :=
  Module[{delaunay = DelaunayTriangulation[Map[Drop[#,-1]&,set]]},
    If[ SameQ[delaunay,$Failed],
	$Failed,
        TriangularSurfacePlot[set,delaunay,opts]
    ]
  ] /; Apply[And,Map[NumberQ,N[Flatten[set]]]]





End[]

EndPackage[] 


(* :Examples:

input = {{4.4,14},{6.7,15.25},{6.9,12.8},{2.1,11.1},{9.5,14.9},{13.2,11.9},
	{10.3,12.3},{6.8,9.5},{3.3,7.7},{.6,5.1},{5.3,2.4},{8.45,4.7},
	{11.5,9.6},{13.8,7.3},{12.9,3.1},{11,1.1}};
query = Map[(#+{.001,.001})&,input]
input3D = Map[{#[[1]],#[[2]],Sqrt[64-(#[[1]]-8)^2-(#[[2]]-8)^2]}&,input]

*****************************************************************************
  convex hull of "input":
		ConvexHull[input]

  plot convex hull of "input":
		PlanarGraphPlot[input,ConvexHull[input]]

*****************************************************************************
  Delaunay triangulation of "input":
		val = DelaunayTriangulation[input]

  plot Delaunay triangulation of "input" where the Delaunay vertex adjacency
	list is not available:
		PlanarGraphPlot[input]

  plot triangulation of "input" where the triangulation vertex adjacency list
	is "val" ("val" need not represent a Delaunay triangulation):
		PlanarGraphPlot[input,val]

  plot triangular surface where the vertex adjacency list representing the
	Delaunay triangulation of the projected input is not available:
		TriangularSurfacePlot[input3D]

  plot triangular surface where the vertex adjacency list representing the
	triangulation of the projected input is "val" ("val" need not
	represent a Delaunay triangulation):
		val = DelaunayTriangulation[Map[Drop[#,-1]&,input3D]];
		TriangularSurfacePlot[input3D,val]

*****************************************************************************
  Voronoi diagram of "input":
		{diagvert,diagval} = VoronoiDiagram[input]

  Voronoi diagram of "input" where the vertex adjacency list representing the
	dual Delaunay triangulation is "val":
		VoronoiDiagram[input,val]

  Voronoi diagram of "input" where the vertex adjacency list representing the
	dual Delaunay triangulation is "val" and the convex hull of "input"
	is "hull":
		{val,hull} = DelaunayTriangulation[input,Hull->True];
		VoronoiDiagram[input,val,hull]

  plot Voronoi diagram of "input" where the Voronoi vertices and vertex
	adjacency list are not available:
		DiagramPlot[input]

  plot diagram of "input" where the diagram vertices are "diagvert" and the
	diagram vertex adjacency list is "diagval" ("diagvert" and "diagval"
	need not represent a Voronoi diagram):
		DiagramPlot[input,diagvert,diagval]

  nearest neighbors of "query" where the Voronoi diagram vertices and vertex
	adjacency list associated with "input" are not available:
		NearestNeighbor[query,input]

  nearest neighbors of "query" where the Voronoi diagram of "input" is
	represented by the vertices "diagvert" and the vertex adjacency list
	"diagval":
		NearestNeighbor[query,diagvert,diagval]

*****************************************************************************
  Delaunay triangulation query:
	DelaunayTriangulationQ[input,val]

	notdelaunayval = Join[{{1,{4,2}},{2,{1,4,3,5}},{3,{2,4,8,7,5}},
				{4,{10,9,3,2,1}}}, Drop[val,4]];
	DelaunayTriangulationQ[input,notdelaunayval]

  Plotting a non-Delaunay triangulation:
	PlanarGraphPlot[input,notdelaunayval]

*)





