(***********************************************************************

                    Mathematica-Compatible Notebook

This notebook can be used on any computer system with Mathematica 4.0,
MathReader 4.0, or any compatible application. The data for the notebook 
starts with the line containing stars above.

To get the notebook into a Mathematica-compatible application, do one of 
the following:

* Save the data starting with the line of stars above into a file
  with a name ending in .nb, then open the file inside the application;

* Copy the data starting with the line of stars above to the
  clipboard, then use the Paste menu command inside the application.

Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode.  Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).

NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing the 
word CacheID, otherwise Mathematica-compatible applications may try to 
use invalid cache data.

For more information on notebooks and Mathematica-compatible 
applications, contact Wolfram Research:
  web: http://www.wolfram.com
  email: info@wolfram.com
  phone: +1-217-398-0700 (U.S.)

Notebook reader applications are available free of charge from 
Wolfram Research.
***********************************************************************)

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[     42372,       1232]*)
(*NotebookOutlinePosition[     43334,       1262]*)
(*  CellTagsIndexPosition[     43290,       1258]*)
(*WindowFrame->Normal*)



Notebook[{

Cell[CellGroupData[{
Cell["Higher order Delauany triangulations.", "Section",
  PageWidth->WindowWidth],

Cell[CellGroupData[{

Cell["k-th Voronoi code.", "Subsection"],

Cell["Initialize  the programs fromn the  kvoronoi package.", "Text"],

Cell["\<\
<<DiscreteMath`Combinatorica`
Off[General::spell1];
<<Graphics`Colors`
<<DiscreteMath`ComputationalGeometry`
<<VertexEnum.m
<<FaceLattice.m

ca[cl_] := Module[
\t\t(* Determinants for circle center *)
        {m,i},
        m = Table[Append[cl[[i]],1],{i,1,3}];
        Det[m]
]
cd[cl_] := Module[
        {m,i},
         m = 
      Table[{Plus @@{cl[[i]][[1]]^2, cl[[i]][[2]]^2}, cl[[i]][[2]], 1},
                {i, 1, 3}];
        -Det[m]
]
ce[cl_] := Module[
        {m,i},
         m = 
      Table[{Plus @@{cl[[i]][[1]]^2, cl[[i]][[2]]^2}, cl[[i]][[1]], 1},
                {i, 1, 3}];
        Det[m]
]
cf[cl_] := Module[
        {m,i},
         m = Table[{Plus @@{cl[[i]][[1]]^2, cl[[i]][[2]]^2}, cl[[i]][[1]],
                 cl[[i]][[2]]}, {i, 1, 3}];
        -Det[m]
]


circlecenter[xyz_List] := Module[
\t\t\t(* computes the centre of the circle, defined by the points <xyz> *)
\t\t{a,d,e},
\t\ta = ca[xyz];
\t\tIf[ a==0,
\t\t\tPrint[\"Dividing by zero\"];
\t\t];
\t\td = cd[xyz];
\t\te = ce[xyz];
\t\t{-d/(2a),-e/(2a)}
\t]

distance[l_List] := Module[
\t\t(* square of euclidean distance *)
\t\t{},
\t(l[[1,1]]-l[[2,1]])^2 +(l[[1,2]]-l[[2,2]])^2
\t\t\t]

allkvoronoi[points_List] := Module[
\t\t(* gets all voronoi circles and the points they contain *)
\t\t{l, rr, lc, circlelist, thiscircle, i, center, radius, otherpoints, k,h, \
j},
    
\t\tl = Length[points];
\t\trr = Range[l];
\t\t(* lc =  #of different circles *)
\t\tlc = Binomial[l,3];
\t\tcirclelist = {};
\t\tthiscircle = {1,2,3};
\t\t(* Use the undocumented function `NextKSubset' *)
\t\tFor[ i=1 , i <= lc,
\t\t\tcenter = circlecenter[ points[[thiscircle]] ];
\t\t\tradius = distance[{points[[thiscircle]][[1]], center}];
\t\t\totherpoints = Complement[ rr, thiscircle];
\t\t\tk = {};
\t\t\tFor[ j=1, j <= (l-3),
\t\t\t\th = otherpoints[[j]];
\t\t\t\tIf[ distance[{points[[h]],center}] < radius,
\t\t\t\t\t\tk=Append[k,h];
\t\t\t\t\t];
\t\t\t\tj++;
\t\t\t];
\t\t\tcirclelist = Append[circlelist, {i, thiscircle, N[center], k}];
\t\t\tthiscircle = NextKSubset[rr,thiscircle];
\t\t\ti++;
\t];
\tcirclelist
]\t

aroundvpoint[circle_,k_] := Module[
\t\t(* gets the edges and faces incident with a <
        k>-th order voronoi vertex (<circle>) *)
\t\t{nr,order, points, edges, nrinpoints, inpoints, regions},
\t\tnr = circle[[1]];
\t\torder = {1,3,2};
\t\torder2={3,1,2};
\t\tpoints = circle[[2]];
\t\tedges= KSubsets[ points,2];
\t\tinpoints = circle[[4]];
\t\tnrinpoints = Length[circle[[4]]];
\t\tIf[ nrinpoints== k-2,
\t\t\tregions = Table[Sort[Join[edges[[order2[[i]]]], inpoints]],{i,3}];
\t\t];
\t\tIf[ nrinpoints == k-1,
\t\t\tregions = Table[Sort[Append[inpoints, points[[order[[i]] ]] ] ],{i,3}];
\t\t];
\t\t{nr, edges, regions}
\t]

edgesandregions[struc_, k_] := Module[
\t\t(* Get all circles that contain k-2or k-1 points *)
\t\t{nr, kstruc,lci, i, ci},
\t\tnr = Length[struc];
\t\tkstruc = {};
\t  For[ i=1, i <= nr,
\t\t\tci = struc[[i]];
\t\t\tlci = Length[ci[[4]]];
\t\t\tIf[ lci== k-2  || lci == k-1,
\t\t\t\tkstruc = Append[kstruc, aroundvpoint[ci, k]];
\t\t\t];
\t\t\ti++;
  \t];
  \tkstruc
]

buildedges[ek_] := Module[
\t\t(* build all edges from the incident faces and edges info in <ek> *)
\t\t{cek, nrv, edges, i, edge, faces, ii, jj, al},
\t\tcek = ek;
\t\tnrv = Length[ cek ];
\t\tedges = {};
\t\tFor[ i=1, i < nrv,
\t\t\tFor[ j=1, j<= Length[ cek[[i,2]] ],
\t\t\t\tedge = cek[[i,2]][[j]];
\t\t\t\tfaces = cek[[i,3]];
\t\t\t\tFor[ ii=i+1, ii <= nrv,
\t\t\t\t\tal = Length[ cek[[ii,2]] ];
\t\t\t\t\tFor[ jj=1, jj <= al,
\t\t\t\t\t\tIf[ cek[[ii,2]][[ jj]] == edge,
\t\t\t\t\t\t\tIf[ Length[Intersection[ faces, cek[[ii,3]] ] ] == 2,
\t\t\t\t\t\t\t\tedges = Append[ edges, {cek[[i,1]], cek[[ii,1]] }];
\t\t\t\t\t\t\t\tcek = Delete[ cek, {ii,2,jj} ];
\t\t\t\t\t\t\t\tjj = 4;
\t\t\t\t\t\t\t  ii = nrv+1;
\t\t\t\t\t\t\t];
\t\t\t\t\t\t];
\t\t\t\t\t\tjj++;
\t\t\t\t\t];
\t\t\t\t\tii++;
\t\t\t\t];
\t\t\t j++;
\t\t];
\ti++;
];
edges
]\t\t
\t\t

showallkvoronoi[ points_] := Module[
\t\t(* Shows all n-1 Voronoi diagrams, defined by $n$ points,  simultanuously \
*)
      
\t\t{n, akv, alllines, k, edreg, edgerel, lines, i, b, e},
\t\tn = Length[ points ];
\t\takv = allkvoronoi[ points ];
\t\talllines = {Thickness[.01]};
\t\tFor[k=1, k <n,
\t\t\tedreg = edgesandregions[ akv, k];
\t\t\tedgerel = buildedges[edreg];
\t\t\tlines = {};
\t\t\tFor[ i = 1, i <= Length[ edgerel],
\t\t\t\tb  = edgerel[[i,1]];
\t\t\t\te = edgerel[[i,2]];
\t\t\t\tlines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
\t\t\t  i++;
\t\t\t];
\t\t\talllines  = Append[alllines, Join[ {Hue[N[k/n]]}, lines]];
\t\t\tk++;
\t\t];
\t\tShow[Graphics[alllines],PlotRange -> All, AspectRatio -> Automatic ,
      Axes -> False]
\t\t]\t\t

showallkvoronoi[ points_, range_] := Module[
\t\t(* Shows all n-1 Voronoi diagrams, defined by $n$ points,  simultanuously \
*)
      
\t\t{n, akv, alllines, k, edreg, edgerel, lines, i, b, e},
\t\tn = Length[ points ];
\t\takv = allkvoronoi[ points ];
\t\talllines = {Thickness[.004]};
\t\tFor[k=1, k <n,
\t\t\tedreg = edgesandregions[ akv, k];
\t\t\tedgerel = buildedges[edreg];
\t\t\tlines = {};
\t\t\tFor[ i = 1, i <= Length[ edgerel],
\t\t\t\tb  = edgerel[[i,1]];
\t\t\t\te = edgerel[[i,2]];
\t\t\t\tlines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
\t\t\t  i++;
\t\t\t];
\t\t\talllines  = Append[alllines, Join[ {Hue[N[k/n]]}, lines]];
\t\t\tk++;
\t\t];
\t\tShow[Graphics[alllines],PlotRange -> range, AspectRatio -> Automatic ,
      Axes -> True]
\t\t]\t\t

showkvoronoi[points_, k_Integer] := Module[
\t\t(* Shows the k-th order Voronoi diagram *)
\t\t{n, akv, alllines, edreg, edgerel, i, lines},
\t\tn = Length[ points ];
\t  akv = allkvoronoi[ points ];
\t  alllines = {Thickness[.002]};
\t\tedreg = edgesandregions[ akv, k];
\t\tedgerel = buildedges[edreg];
\t\tlines = {};
\t\tFor[ i = 1, i <= Length[ edgerel],
\t\t\t\tb  = edgerel[[i,1]];
\t\t\t\te = edgerel[[i,2]];
\t\t\t\tlines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
\t\t\t  i++;
\t\t];
\t\talllines = Join[{Thickness[.01], Hue[N[k/n]]}, lines];
\t\tShow[Graphics[alllines], PlotRange -> All,AspectRatio -> Automatic]
\t]\t

showsomekvoronoi[ points_,l_List] := Module[
\t\t(* Shows the combined Voronoi diagrams for the orders present in <l> *)
\t\t{n, akv, alllines, k, edreg, edgerel, lines, i, b, e,j},
\t\tn = Length[ points ];
\t\takv = allkvoronoi[ points ];
\t\talllines = {Thickness[.01]};
\t\tFor[j=1, j <=Length[l],
\t\t\tk = l[[j]];
\t\t\tedreg = edgesandregions[ akv, k];
\t\t\tedgerel = buildedges[edreg];
\t\t\tlines = {};
\t\t\tFor[ i = 1, i <= Length[ edgerel],
\t\t\t\tb  = edgerel[[i,1]];
\t\t\t\te = edgerel[[i,2]];
\t\t\t\tlines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
\t\t\t  i++;
\t\t\t];
\t\t\talllines  = Append[alllines, Join[ {Hue[N[k/n]]}, lines]];
\t\t\tj++;
\t\t];
\t\tShow[Graphics[alllines],PlotRange -> All, AspectRatio -> Automatic ,
      Axes -> True]
\t\t]\t\t

showsomekvoronoi[ points_,l_List, range_] := Module[
\t\t(* Shows the combined Voronoi diagrams for the orders present in <l> *)
\t\t{n, akv, alllines, k, edreg, edgerel, lines, i, b, e,j},
\t\tn = Length[ points ];
\t\takv = allkvoronoi[ points ];
\t\talllines = {Thickness[.01]};
\t\tFor[j=1, j <=Length[l],
\t\t\tk = l[[j]];
\t\t\tedreg = edgesandregions[ akv, k];
\t\t\tedgerel = buildedges[edreg];
\t\t\tlines = {};
\t\t\tFor[ i = 1, i <= Length[ edgerel],
\t\t\t\tb  = edgerel[[i,1]];
\t\t\t\te = edgerel[[i,2]];
\t\t\t\tlines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
\t\t\t  i++;
\t\t\t];
\t\t\talllines  = Append[alllines, Join[ {Hue[N[k/n]]}, lines]];
\t\t\tj++;
\t\t];
\t\tShow[Graphics[alllines],PlotRange -> range, AspectRatio -> Automatic ,
      Axes -> True]
\t\t]\t\t

voronoiflist[points_] := Module[
\t\t(* Returns all Voronoi faces defined by <points> *)
\t\t{n, akv, allfaces, i, faces, ear, j},
\t\tn = Length[points];
\t\takv = allkvoronoi[points];
\t\tallfaces = {{Range[n]}};
\t\tFor[ i=n-1, i >= 1,
\t\t\tfaces = {};
\t\t\tear = edgesandregions[akv,i];
\t\t\tFor[j=1, j <= Length[ear],
\t\t\t\tfaces = Join[faces, ear[[j,3]]];
\t\t\t\tj++;
\t\t\t];
\t\t\tallfaces = Append[allfaces, Union[faces]];
\t\t\ti--;
\t\t];
\t\tAppend[allfaces,{{}}]
\t\t]

randomnvor[n_Integer] := Module[
\t\t(* computes the euler chararacteristic for a set of <n> random points *)
\t\t{r, mp, flist, fvector, avector, euler},
\t\tr := Random[Integer,{0,1000}];
\t\tmp = Table[{r,r},{i,n}];
\t\tPrint[\"points = \",mp];
\t\tflist = voronoiflist[mp];
\t\tfvector=FVector[flist];
\t\tPrint[\"fvector = \", fvector];
\t\tavector=Table[(-1)^(i),{i,Length[fvector]}];
\t   euler = avector.fvector;
\t\tPrint[\"\[CurlyKappa] = \", euler];
\t\teuler
\t\t]

cvector[points_List] := Module[
\t\t(* Returns the number of circles of order $i$ for i=0,...,n-3 *)
\t\t{n, nr, cvec, akv, i},
\t\tn = Length[points];
\t\tnr = Binomial[n,3];
\t\tcvec = Table[0, {i, n-2}];
\t  akv = allkvoronoi[points];
\t\tFor[i=1, i <= nr,
\t\t\tcvec[[Length[ akv[[i,4]] ] + 1]]++;
\t\t\ti++;
\t\t];
\t\tcvec
]

 vvector[points_List] := Module[
\t\t(* Returns the number of vertices of the $k$-th order VD  for k=1,...,n-1 \
*)
      
\t\t{cvec,l,i},
\t\tcvec = Append[Prepend[ cvector[points],0],0];
\t\tl = Length[cvec];
\t\tTable[cvec[[i]] + cvec[[i+1]], {i,l-1}]
]

finfvector[points_] := Module[
\t\t{akv, finfvec, k, edreg, occur, vertices, lv, valencies, nrunbound, j},
\t\takv = allkvoronoi[ points ];
\t\tn= Length[points];
\t\tfinfvec = {};
\t\tFor[k=1, k <n,
\t\t\tedreg = edgesandregions[akv,k];
\t\t\toccur =Sort[Flatten[ buildedges[ edreg ] ] ];
\t\t\tvertices = Union[occur];
\t\t\tlv = Length[vertices];
\t\t\tvalencies = Table[Count[occur,vertices[[i]] ],{i,lv}];
\t\t\tnrunbound = 0;
\t\t\tFor[ j = 1, j <= lv,
\t\t\t\tnrunbound = nrunbound + (3-valencies[[j]]);
\t\t\t\tj++;
\t\t\t];
\t\t\t(* added later, RCL *)
\t\t\tIf[ lv == 0,
\t\t\t  nrunbound = 3;
\t\t\t];
\t\t\tfinfvec = Append[finfvec, nrunbound];
\t\t\tk++;
\t\t];
\t\tfinfvec
\t\t]\t\t

fvector[points_] := FVector[voronoiflist[points]]

evector[points_] := Module[
\t\t(* Returns the number of edges in all $k$-
        th order VD's  Note that the number of unbounded edges equals the 
          number of unbounded regions *)
\t\t{akv, finfvec, k, edreg, occur, vertices, lv, valencies, nrunbound, j},
\t\takv = allkvoronoi[ points ];
\t\tn= Length[points];
\t\tevec  = {};
\t\tfinfvec = {};
\t\tFor[k=1, k <n,
\t\t\tedreg = edgesandregions[akv,k];
\t\t\tbedges =  buildedges[ edreg ];
\t\t\toccur =Sort[Flatten[ bedges ] ];
\t\t\tvertices = Union[occur];
\t\t\tlv = Length[vertices];
\t\t\tvalencies = Table[Count[occur,vertices[[i]] ],{i,lv}];
\t\t\tnrunbound = 0;
\t\t\tFor[ j = 1, j <= lv,
\t\t\t\tnrunbound = nrunbound + (3-valencies[[j]]);
\t\t\t\tj++;
\t\t\t];
\t\t\t(* added later, RCL *)
\t\t\tIf[ lv == 0,
\t\t\t  nrunbound = 3;
\t\t\t];
\t\t\tfinfvec = Append[finfvec, nrunbound+Length[bedges]];
\t\t\tk++;
\t\t];
\t\tfinfvec
\t\t]\t\t


ffpoly[points_] := Module[
\t\t{rfvec, poly, i},
\t\trfvec = Reverse[fvector[points]];
\t\tpoly = 0;
\t\tFor[ i =1, i <= Length[rfvec],
\t\t\tpoly = poly + rfvec[[i]] t^(i-1);
\t\t\ti++;
\t\t];
\t\tFactor[poly]
\t]\t

chi[points_] := Module[
\t\t{flist, fvector, avector},
\t\tflist = voronoiflist[points];
\t\tfvector = FVector[flist];
\t\tavector= Table[(-1)^i,{i,Length[fvector]}];
\t\tavector.fvector+1
\t];

totalvertices[n_] := (1/3)n(n-1)(n-2)
totaledges[n_] := (1/2)n(n-1)^2
totalregions[n_] := (1/6)n(n^2+5)

fnk[n_,k_] := 2(k+1)(n-k)-(k+1)-(n-k)+2

bounds[vlist_] := Module[
\t\t{l,i,j},
\t\tl = Length[vlist];
\t\tTable[{Min[Table[vlist[[i,j]],{i,l}]],Max[Table[vlist[[i,j]],{i,l}]]},
    \t{j,Length[vlist[[1]]]}]
]

showsomecircles[points_List,k_Integer] := Module[
\t\t(* Shows the <points> and the Voronoi circles of order <k> *)
\t\t{n,akv,circles,i,center,radius},
\t\tn=Length[points];
\t\tIf[ k>n-3, Return[Print[\"k is too big...\"]] ];
\t\takv = allkvoronoi[points];
\t\tcircles = { Hue[N[k/n] ]};
\t\tFor[ i=1, i <= Length[akv],
\t\t\tIf[Length[ akv[[i,4]] ]==k,
\t\t\t\t  center =  akv[[i,3]];
\t\t\t\t  radius = distance[ {points[[  akv[[i,2,1]] ]], center }]
\t\t\t\t\t\t// Sqrt//N;
\t\t\t\t  circles = Append[circles,Circle[center,radius]];
\t\t\t];
\t\t\ti++;
\t\t];
\t\tl1 = Table[Text[i,points[[i]] ],{i,n}];
\t\t{Show[Graphics[{circles,l1}],
\t AspectRatio -> Automatic,
\t\t\tAxes -> False\t], Length[circles]-1}
]

<<$HOME/cas/matica/kvoronoi/ihull.m

bisector[lijst_] := Module[
\t\t\t{p1,a,b,p2,c,d},
\t\t
\t\t  p1 = lijst[[1]];
\t\t  a = p1[[1]]; b = p1[[2]];
\t\t  p2 = lijst[[2]];
\t\t  c = p2[[1]]; d = p2[[2]];
\t\t  If[d==b,1/2 (a +c),
\t\t  (a-c)/(d-b) x + (1/2) ( d+b +(c c-a a)/(d-b))]
\t\t\t]

bisector[p1_, p2_] := Module[
\t\t\t{a,b,c,d},
\t\t
\t\t  a = p1[[1]]; b = p1[[2]];
\t\t  c = p2[[1]]; d = p2[[2]];
\t\t  If[d==b,1/2 (a +c),
\t\t  (a-c)/(d-b) x + (1/2) ( d+b +(c c-a a)/(d-b))]
\t\t\t]


bisectorarr[points_] := Module[
\t\t{n,k,bis,x},
\t\tn = Length[points];
\t\tk = KSubsets[Range[n],2];
\t\tbis = Table[bisector[ points[[k[[i,1]] ]],points[[k[[i,2]] ]] ],
\t\t\t{i,Length[k]}]//N;
\t\tbis
\t\t]

pairtoposition[pair_,n_] := Module[
\t\t(* computes the position of an ordered pair in the enumeration of 
        combinations of length two of [n] *)
\t\t{p1,k},
\t\tp1=pair[[1]];
\t\tSum[n-k,{k,1,p1-1}]+ pair[[2]]-p1
\t\t]

oneslope[p1_,p2_] := Module[
\t\t\t{a,b,c,d},
\t\t  a = p1[[1]]; b = p1[[2]];
\t\t  c = p2[[1]]; d = p2[[2]];
\t\t  If[d==b,1/2 (a +c)//N, (a-c)/(d-b) //N]
\t\t\t]

allslopes[points_] := Module[
\t\t{n,k,i},
\t\tn = Length[points];
\t\tk = KSubsets[Range[n],2];
\t Table[oneslope[ points[[ k[[i,1]] ]],points[[ k[[i,2]] ]] ],
\t\t\t{i,Length[k]} ]
\t\t]

buildedgesinf[ek_] := Module[
\t\t(* build all edges from the incident faces and edges info in <ek> *)
\t\t{cek, nrv, edges, i,j, edge, faces, ii, jj, al},
\t\tcek = ek;
\t\tnrv = Length[ cek ];
\t\tedges = {};
\t\tinfedges={};
\t\tFor[ i=1, i < nrv,
\t\t\tFor[ j=1, j<= Length[ cek[[i,2]] ],
\t\t\t\tedge = cek[[i,2]][[j]];
\t\t\t\tfaces = cek[[i,3]];
\t\t\t\tFor[ ii=i+1, ii <= nrv,
\t\t\t\t\tal = Length[ cek[[ii,2]] ];
\t\t\t\t\tFor[ jj=1, jj <= al,
\t\t\t\t\t\tIf[ cek[[ii,2]][[ jj]] == edge,
\t\t\t\t\t\t\tIf[ Length[Intersection[ faces, cek[[ii,3]] ] ] == 2,
\t\t\t\t\t\t\t\tedges = Append[ edges, {cek[[i,1]], cek[[ii,1]] }];
\t\t\t\t\t\t\t\tcek = Delete[ cek, {ii,2,jj} ];
\t\t\t\t\t\t\t\tjj = 5;
\t\t\t\t\t\t\t  ii = nrv+1;
\t\t\t\t\t\t\t];
\t\t\t\t\t\t];
\t\t\t\t\t\tjj++;
\t\t\t\t\t];
\t\t\t\t\tii++;
\t\t\t\t];
\t\t\t\tIf[ jj < 5, infedges=Append[infedges,{i,edge}]];
\t\t\t j++;
\t\t];
\ti++;
];
For[k=1, k <= Length[cek[[nrv,2]]],
\t\t\tinfedges = Append[infedges, {i,cek[[nrv,2,k]]} ];
\t\t\tk++;
];
{edges,infedges}
]\t\t

sign[p_,q_,r_] := Module[
\t\t{},
     Det[{{1,p[[1]],p[[2]]},
\t\t\t\t\t{1,q[[1]],q[[2]]},
\t\t\t\t\t{1,r[[1]],r[[2]]}}]//Sign
\t\t]

showkvoronoiinf[points_, k_Integer] := Module[
\t\t{lines,n, akv,edreg, edgeinfo,edgerel,i,b, e,infinfo, infedge,thiscircle,
      slopes, slope, centre,radius,ss,s1,p1label, p2label, p3label,  p1, p2, 
      p3,alllines, tpoints, cyoung, testsign},
                If[ k < 1, Return[{}]];
\t\tlines = {};
\t\tn = Length[ points ];
\t\tIf[ k >= n, Return[{}]];
\t  akv = allkvoronoi[ points ];
\t\tedreg = edgesandregions[ akv, k];
\t\tedgeinfo = buildedgesinf[edreg];
\t\tedgerel = edgeinfo[[1]];
\t\t(* Add the edges *)
\t\tFor[ i = 1, i <= Length[ edgerel],
\t\t\t\tb  = edgerel[[i,1]];
\t\t\t\te = edgerel[[i,2]];
\t\t\t\tlines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
\t\t\t  i++;
\t\t];
\t\t(* Now we add the halflines *)
\t   slopes = allslopes[points];
\t\tinfinfo = edgeinfo[[2]];
\t\tFor[i=1, i <= Length[infinfo],
\t\t\tinfedge = infinfo[[i]];
\t\t\tthiscircle=edreg[[infedge[[1]],1]] ;
\t\t\t(* determine whether <thiscircle> is young or old *)
\t  \thcircle = Length[ akv[[thiscircle,4]] ];
\t\t\tIf[ hcircle+2==k, cyoung = -1, cyoung = 1];
\t\t\tslope= slopes[[pairtoposition[infedge[[2]],n ] ]];
\t\t\tcentre= akv[[ edreg[[infedge[[1]],1]],3]];
\t\t\tradius = N[Sqrt[ distance[ {centre, points[[infedge[[2,1]] ]]} ] ] ];
\t\t\tss = (2.radius)/(Sqrt[1+slope^2]){1,slope};
\t\t\ts1 = ss+centre;
\t\t\tp1label=infedge[[2,1]];
\t\t\tp2label=infedge[[2,2]];
\t\t\tp3label= Complement[akv[[thiscircle,2]],{p1label,p2label}][[1]];
\t\t\tp1 = points[[p1label]];
\t\t\tp2 = points[[p2label]];
\t\t\tp3= points[[p3label]];
\t\t\ttestsign = sign[p1,p2, s1]  sign[ p1, p2, p3 ];
\t\t\tIf[ cyoung== testsign,
\t\t\t\t\t\t\t\t\ts1 = -ss + centre;
\t\t\t];
\t\t\tlines = Append[lines,Line[{centre,s1}]];
\t\t\ti++;
\t\t];
\t\talllines = Join[{Thickness[.005], Hue[N[k/n]]}, lines];
\t\ttpoints = Table[Text[i,points[[i]] ],{i,n}];
     Show[Graphics[{alllines,tpoints}], PlotRange -> All,
      AspectRatio -> Automatic]
\t]\t

showkvoronoiinf[points_, k_Integer,akv_ ]:= Module[
\t\t{lines,n, akv,edreg, edgeinfo,edgerel,i,b, e,infinfo, infedge,thiscircle,
      slopes, slope, centre,radius,ss,s1,p1label, p2label, p3label,  p1, p2, 
      p3,alllines, tpoints, cyoung, testsign},
\t\tlines = {};
\t\tn = Length[ points ];
\t\tIf[ k >= n, Return[{}]];
\t\tedreg = edgesandregions[ akv, k];
\t\tedgeinfo = buildedgesinf[edreg];
\t\tedgerel = edgeinfo[[1]];
\t\t(* Add the edges *)
\t\tFor[ i = 1, i <= Length[ edgerel],
\t\t\t\tb  = edgerel[[i,1]];
\t\t\t\te = edgerel[[i,2]];
\t\t\t\tlines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
\t\t\t  i++;
\t\t];
\t\t(* Now we add the halflines *)
\t   slopes = allslopes[points];
\t\tinfinfo = edgeinfo[[2]];
\t\tFor[i=1, i <= Length[infinfo],
\t\t\tinfedge = infinfo[[i]];
\t\t\tthiscircle=edreg[[infedge[[1]],1]] ;
\t\t\t(* determine whether <thiscircle> is young or old *)
\t  \thcircle = Length[ akv[[thiscircle,4]] ];
\t\t\tIf[ hcircle+2==k, cyoung = -1, cyoung = 1];
\t\t\tslope= slopes[[pairtoposition[infedge[[2]],n ] ]];
\t\t\tcentre= akv[[ edreg[[infedge[[1]],1]],3]];
\t\t\tradius = N[Sqrt[ distance[ {centre, points[[infedge[[2,1]] ]]} ] ] ];
\t\t\tss = (2.radius)/(Sqrt[1+slope^2]){1,slope};
\t\t\ts1 = ss+centre;
\t\t\tp1label=infedge[[2,1]];
\t\t\tp2label=infedge[[2,2]];
\t\t\tp3label= Complement[akv[[thiscircle,2]],{p1label,p2label}][[1]];
\t\t\tp1 = points[[p1label]];
\t\t\tp2 = points[[p2label]];
\t\t\tp3= points[[p3label]];
\t\t\ttestsign = sign[p1,p2, s1]  sign[ p1, p2, p3 ];
\t\t\tIf[ cyoung== testsign,
\t\t\t\t\t\t\t\t\ts1 = -ss + centre;
\t\t\t];
\t\t\tlines = Append[lines,Line[{centre,s1}]];
\t\t\ti++;
\t\t];
\t\talllines = Join[{Thickness[.005], Hue[N[k/n]]}, lines];
\t\ttpoints = Table[Text[i,points[[i]] ],{i,n}];
     Show[Graphics[{alllines,tpoints}], PlotRange -> All,
      AspectRatio -> Automatic]
\t]\t

\
\>", "Input",
  PageWidth->Infinity,
  InitializationCell->True,
  ShowSpecialCharacters->False,
  FormatType->InputForm]
}, Closed]],

Cell[CellGroupData[{

Cell["Second order Delauany triangulations.", "Subsection"],

Cell["We start with a list S of 7 points.", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(\(points = {{0.2675984295549247`, 
            0.8535461718417067`}, {0.1305841735968181`, 
            0.08773057333372192`}, {0.8833203358116403`, 
            0.7133112243927627`}, {0.7148603434510898`, 
            0.6061062367112292`}, {0.3088088797261596`, 
            0.3988386007048817`}, {0.86262835750783`, 
            0.19703287214655868`}, {0.7603336554833071`, 
            0.12006112332043573`}};\)\)], "Input",
  PageWidth->WindowWidth],

Cell["\<\
For these points we know how to build the 2nd order Voronoi \
diagram. More precisely, we can make a list <duos> of the cells V2(pi,pj) \
that occur in the 2nd order diagram V2(S).\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[{
    \(\(akv = allkvoronoi[points];\)\), "\[IndentingNewLine]", 
    \(\(e2\  = \ edgesandregions[akv, 2];\)\), "\[IndentingNewLine]", 
    \(duos = 
      Union[Flatten[Table[e2[\([i, 3]\)], {i, Length[e2]}], 1]]\)}], "Input",
  PageWidth->WindowWidth],

Cell["\<\
For every pair (pi,pj) of points in <duos> we make a new point \
(pi+pj)/2. In this way we associate a point (pi+pj)/2 to every cel V2(pi,pj) \
in <duos>, resulting in a list <duop>.\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[{
    \(h[duo_]\  := \((\ 
          points[\([duo[\([1]\)]]\)] + points[\([duo[\([2]\)]]\)])\)/
        2\), "\[IndentingNewLine]", 
    \(\(duop = 
        Table[\ h[duos[\([i]\)]\ ], {i, \ Length[duos]}];\)\)}], "Input",
  PageWidth->WindowWidth],

Cell["\<\
Ofcourse we can compute the ordinary Delaunay triangulation D(duop) \
for the point set <duop>. We do so and plot D(duop) together with \
V2(S).\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[{
    \(\(bit = { .01,  .02};\)\), "\[IndentingNewLine]", 
    \(\(teksten = 
        Table[\ Text[duos[\([i]\)], \ duop[\([i]\)]\ ], \ {i, \ 
            Length[duos]}];\)\), "\[IndentingNewLine]", 
    \(\(s2 = showkvoronoiinf[points, 2];\)\), "\[IndentingNewLine]", 
    \(\(s3\  = \ 
        PlanarGraphPlot[duop, 
          LabelPoints\  \[Rule] \ False];\)\), "\[IndentingNewLine]", 
    \(\(Show[s2, \ 
        Graphics[{{PointSize[ .03]}, Point\  /@ \ duop, \ teksten}], 
        s3, \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- .1\), 
              1.1}, {\(- .1\), 1.1}}, \ 
        AspectRatio\  \[Rule] \ Automatic];\)\)}], "Input",
  PageWidth->WindowWidth],

Cell["\<\
We define a graph D2(S) on the points in <duop>.  Two vertices pij \
and pkl in <duop> are connected by an edge if and only if the corresponding \
Voronoi cells  V2(pi,pj) and V2(pk,pl) have an edge in common. \
\>", "Text",
  PageWidth->WindowWidth],

Cell["\<\
In order to compute those edges, we use the following program \
duoedges[eandr_, duos_], that takes the list of edges and regions <eandr> and \
the list of cells <duos>  in the 2nd order diagram as its input and returns a \
lsit of edges between elements of <duos>.\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(duoedges[\ eandr_, \ duos_]\  := \ 
      Module[\[IndentingNewLine]{nv, \ iedges, \ edges, \ i, \ duo, \ 
          thisvertex, \ pos, \ j, \ k}, \[IndentingNewLine]nv\  = \ 
          Length[eandr]; \[IndentingNewLine]iedges = 
          Table[eandr[\([i, 3]\)], {i, 
              nv}]; \[IndentingNewLine]edges = {}; \[IndentingNewLine]For[\ 
          i = 1, \ i\  \[LessEqual] \ 
            Length[duos], \[IndentingNewLine]duo = \ 
            duos[\([i]\)]; \[IndentingNewLine]For[j = 1, \ 
            j\  \[LessEqual] \ 
              Length[iedges], \[IndentingNewLine]thisvertex\  = \ 
              iedges[\([j]\)]; \[IndentingNewLine]pos\  = \ 
              Position[\ thisvertex, \ duo]; \[IndentingNewLine]If[
              Length[pos]\  > \ 0, \[IndentingNewLine]pos = 
                pos[\([1, 1]\)]; \[IndentingNewLine]For[\ k = 1, \ 
                k\  \[LessEqual] \ 3, \[IndentingNewLine]If[\ 
                  pos\  \[NotEqual] \ k, \[IndentingNewLine]edges\  = \ 
                    Append[edges, \ {i, \(Position[duos, \ 
                            thisvertex[\([k]\)]\ ]\)[\([1, 
                            1]\)]}]\[IndentingNewLine]]; \
\[IndentingNewLine]\(\(k++\);\)\[IndentingNewLine]];\[IndentingNewLine]]; \
\[IndentingNewLine]\(\(j++\);\)\[IndentingNewLine]]; \[IndentingNewLine]\(i\
++\);\[IndentingNewLine]]; \[IndentingNewLine]Union[
          Table[Union[edges[\([i]\)]], {i, 
              Length[edges]}]]\[IndentingNewLine]]\)], "Input",
  PageWidth->WindowWidth,
  InitializationCell->True],

Cell["We perform the program and draw all edges that it returns .", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(eded = duoedges[\ e2, duos]\)], "Input",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(\(lines = 
        Table[Line[duop[\([\ eded[\([i]\)]\ ]\)]\ ], {i, 
            Length[eded]}];\)\)], "Input",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(\(s4 = \ 
        Show[Graphics[{\[IndentingNewLine]Orange, \[IndentingNewLine]lines\
\[IndentingNewLine]}], \ AspectRatio\  \[Rule] \ Automatic];\)\)], "Input",
  PageWidth->WindowWidth],

Cell["\<\
We superimpose the result on the the Delaunay triangulation D(duop) \
that we have computed above. It's not the same triangulation.\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(\(Show[s2, s3, s4, \ 
        Graphics[{{PointSize[ .05]}, Point\  /@ \ duop, \ 
            teksten}], \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- .1\), 
              1.1}, {\(- .1\), 1.1}}, \ 
        AspectRatio\  \[Rule] \ Automatic];\)\)], "Input",
  PageWidth->WindowWidth],

Cell["\<\
We group all commands. This gives a program that generates a random \
point set of <n> points and compute both triangulations for it.\
\>", "Text"],

Cell[BoxData[
    \(randon2ddelaunay[n_]\  := \ 
      Module[\[IndentingNewLine]{points, r, \ akv, \ e2, \ duos, \ bit, h, \ 
          teksten, \ s2, \ s3, \ s4}, \[IndentingNewLine]r\  := \ 
          Random[]; \[IndentingNewLine]points\  = \ 
          Table[{r, r}, {i, n}]; \[IndentingNewLine]akv = 
          allkvoronoi[points]; \[IndentingNewLine]e2\  = \ 
          edgesandregions[akv, 2]; \[IndentingNewLine]duos = 
          Union[Flatten[Table[e2[\([i, 3]\)], {i, Length[e2]}], 
              1]]; \[IndentingNewLine]h[
            duo_]\  := \((\ 
              points[\([duo[\([1]\)]]\)] + points[\([duo[\([2]\)]]\)])\)/
            2; \[IndentingNewLine]duop = 
          Table[\ h[duos[\([i]\)]\ ], {i, \ 
              Length[duos]}]; \[IndentingNewLine]bit = { .01,  .02}; \
\[IndentingNewLine]teksten = 
          Table[\ Text[duos[\([i]\)], \ duop[\([i]\)]\ ], \ {i, \ 
              Length[duos]}]; \[IndentingNewLine]s2 = 
          showkvoronoiinf[points, 2]; \[IndentingNewLine]s3\  = \ 
          PlanarGraphPlot[duop, 
            LabelPoints\  \[Rule] \ False]; \[IndentingNewLine]Show[s2, \ 
          Graphics[{{PointSize[ .03]}, Point\  /@ \ duop, \ teksten}], 
          s3, \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- .1\), 
                1.1}, {\(- .1\), 1.1}}, \ 
          AspectRatio\  \[Rule] \ Automatic]; \[IndentingNewLine]eded = 
          duoedges[\ e2, duos]; \[IndentingNewLine]lines = 
          Table[Line[duop[\([\ eded[\([i]\)]\ ]\)]\ ], {i, 
              Length[eded]}]; \[IndentingNewLine]s4 = \ 
          Show[Graphics[{\[IndentingNewLine]Orange, \[IndentingNewLine]lines\
\[IndentingNewLine]}], \ 
            AspectRatio\  \[Rule] \ Automatic]; \[IndentingNewLine]Show[s2, 
          s3, s4, \ 
          Graphics[{{PointSize[ .05]}, Point\  /@ \ duop, \ 
              teksten}], \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- \
.1\), 1.1}, {\(- .1\), 1.1}}, \ 
          AspectRatio\  \[Rule] \ Automatic];\[IndentingNewLine]]\)], "Input",\

  InitializationCell->True],

Cell["Initialize  the programs fromn the  kvoronoi package.", "Text"],

Cell["We start with a list S of 7 points.", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(\(points = {{0.2675984295549247`, 
            0.8535461718417067`}, {0.1305841735968181`, 
            0.08773057333372192`}, {0.8833203358116403`, 
            0.7133112243927627`}, {0.7148603434510898`, 
            0.6061062367112292`}, {0.3088088797261596`, 
            0.3988386007048817`}, {0.86262835750783`, 
            0.19703287214655868`}, {0.7603336554833071`, 
            0.12006112332043573`}};\)\)], "Input",
  PageWidth->WindowWidth],

Cell["\<\
For these points we know how to build the 2nd order Voronoi \
diagram. More precisely, we can make a list <duos> of the cells V2(pi,pj) \
that occur in the 2nd order diagram V2(S).\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[{
    \(\(akv = allkvoronoi[points];\)\), "\[IndentingNewLine]", 
    \(\(e2\  = \ edgesandregions[akv, 2];\)\), "\[IndentingNewLine]", 
    \(duos = 
      Union[Flatten[Table[e2[\([i, 3]\)], {i, Length[e2]}], 1]]\)}], "Input",
  PageWidth->WindowWidth],

Cell["\<\
For every pair (pi,pj) of points in <duos> we make a new point \
(pi+pj)/2. In this way we associate a point (pi+pj)/2 to every cel V2(pi,pj) \
in <duos>, resulting in a list <duop>.\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[{
    \(h[duo_]\  := \((\ 
          points[\([duo[\([1]\)]]\)] + points[\([duo[\([2]\)]]\)])\)/
        2\), "\[IndentingNewLine]", 
    \(\(duop = 
        Table[\ h[duos[\([i]\)]\ ], {i, \ Length[duos]}];\)\)}], "Input",
  PageWidth->WindowWidth],

Cell["\<\
Ofcourse we can compute the ordinary Delaunay triangulation D(duop) \
for the point set <duop>. We do so and plot D(duop) together with \
V2(S).\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[{
    \(\(bit = { .01,  .02};\)\), "\[IndentingNewLine]", 
    \(\(teksten = 
        Table[\ Text[duos[\([i]\)], \ duop[\([i]\)]\ ], \ {i, \ 
            Length[duos]}];\)\), "\[IndentingNewLine]", 
    \(\(s2 = showkvoronoiinf[points, 2];\)\), "\[IndentingNewLine]", 
    \(\(s3\  = \ 
        PlanarGraphPlot[duop, 
          LabelPoints\  \[Rule] \ False];\)\), "\[IndentingNewLine]", 
    \(\(Show[s2, \ 
        Graphics[{{PointSize[ .03]}, Point\  /@ \ duop, \ teksten}], 
        s3, \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- .1\), 
              1.1}, {\(- .1\), 1.1}}, \ 
        AspectRatio\  \[Rule] \ Automatic];\)\)}], "Input",
  PageWidth->WindowWidth],

Cell["\<\
We define a graph D2(S) on the points in <duop>.  Two vertices pij \
and pkl in <duop> are connected by an edge if and only if the corresponding \
Voronoi cells  V2(pi,pj) and V2(pk,pl) have an edge in common. \
\>", "Text",
  PageWidth->WindowWidth],

Cell["\<\
In order to compute those edges, we use the following program \
duoedges[eandr_, duos_], that takes the list of edges and regions <eandr> and \
the list of cells <duos>  in the 2nd order diagram as its input and returns a \
lsit of edges between elements of <duos>.\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(duoedges[\ eandr_, \ duos_]\  := \ 
      Module[\[IndentingNewLine]{nv, \ iedges, \ edges, \ i, \ duo, \ 
          thisvertex, \ pos, \ j, \ k}, \[IndentingNewLine]nv\  = \ 
          Length[eandr]; \[IndentingNewLine]iedges = 
          Table[eandr[\([i, 3]\)], {i, 
              nv}]; \[IndentingNewLine]edges = {}; \[IndentingNewLine]For[\ 
          i = 1, \ i\  \[LessEqual] \ 
            Length[duos], \[IndentingNewLine]duo = \ 
            duos[\([i]\)]; \[IndentingNewLine]For[j = 1, \ 
            j\  \[LessEqual] \ 
              Length[iedges], \[IndentingNewLine]thisvertex\  = \ 
              iedges[\([j]\)]; \[IndentingNewLine]pos\  = \ 
              Position[\ thisvertex, \ duo]; \[IndentingNewLine]If[
              Length[pos]\  > \ 0, \[IndentingNewLine]pos = 
                pos[\([1, 1]\)]; \[IndentingNewLine]For[\ k = 1, \ 
                k\  \[LessEqual] \ 3, \[IndentingNewLine]If[\ 
                  pos\  \[NotEqual] \ k, \[IndentingNewLine]edges\  = \ 
                    Append[edges, \ {i, \(Position[duos, \ 
                            thisvertex[\([k]\)]\ ]\)[\([1, 
                            1]\)]}]\[IndentingNewLine]]; \
\[IndentingNewLine]\(\(k++\);\)\[IndentingNewLine]];\[IndentingNewLine]]; \
\[IndentingNewLine]\(\(j++\);\)\[IndentingNewLine]]; \[IndentingNewLine]\(i\
++\);\[IndentingNewLine]]; \[IndentingNewLine]Union[
          Table[Union[edges[\([i]\)]], {i, 
              Length[edges]}]]\[IndentingNewLine]]\)], "Input",
  PageWidth->WindowWidth,
  InitializationCell->True],

Cell["We perform the program and draw all edges that it returns .", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(eded = duoedges[\ e2, duos]\)], "Input",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(\(lines = 
        Table[Line[duop[\([\ eded[\([i]\)]\ ]\)]\ ], {i, 
            Length[eded]}];\)\)], "Input",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(\(s4 = \ 
        Show[Graphics[{\[IndentingNewLine]Orange, \[IndentingNewLine]lines\
\[IndentingNewLine]}], \ AspectRatio\  \[Rule] \ Automatic];\)\)], "Input",
  PageWidth->WindowWidth],

Cell["\<\
We superimpose the result on the the Delaunay triangulation D(duop) \
that we have computed above. It's not the same triangulation.\
\>", "Text",
  PageWidth->WindowWidth],

Cell[BoxData[
    \(\(Show[s2, s3, s4, \ 
        Graphics[{{PointSize[ .05]}, Point\  /@ \ duop, \ 
            teksten}], \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- .1\), 
              1.1}, {\(- .1\), 1.1}}, \ 
        AspectRatio\  \[Rule] \ Automatic];\)\)], "Input",
  PageWidth->WindowWidth],

Cell["\<\
We group all commands. This gives a program that generates a random \
point set of <n> points and compute both triangulations for it.\
\>", "Text"],

Cell[BoxData[
    \(randon2ddelaunay[n_]\  := \ 
      Module[\[IndentingNewLine]{points, r, \ akv, \ e2, \ duos, \ bit, h, \ 
          teksten, \ s2, \ s3, \ s4}, \[IndentingNewLine]r\  := \ 
          Random[]; \[IndentingNewLine]points\  = \ 
          Table[{r, r}, {i, n}]; \[IndentingNewLine]akv = 
          allkvoronoi[points]; \[IndentingNewLine]e2\  = \ 
          edgesandregions[akv, 2]; \[IndentingNewLine]duos = 
          Union[Flatten[Table[e2[\([i, 3]\)], {i, Length[e2]}], 
              1]]; \[IndentingNewLine]h[
            duo_]\  := \((\ 
              points[\([duo[\([1]\)]]\)] + points[\([duo[\([2]\)]]\)])\)/
            2; \[IndentingNewLine]duop = 
          Table[\ h[duos[\([i]\)]\ ], {i, \ 
              Length[duos]}]; \[IndentingNewLine]bit = { .01,  .02}; \
\[IndentingNewLine]teksten = 
          Table[\ Text[duos[\([i]\)], \ duop[\([i]\)]\ ], \ {i, \ 
              Length[duos]}]; \[IndentingNewLine]s2 = 
          showkvoronoiinf[points, 2]; \[IndentingNewLine]s3\  = \ 
          PlanarGraphPlot[duop, 
            LabelPoints\  \[Rule] \ False]; \[IndentingNewLine]Show[s2, \ 
          Graphics[{{PointSize[ .03]}, Point\  /@ \ duop, \ teksten}], 
          s3, \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- .1\), 
                1.1}, {\(- .1\), 1.1}}, \ 
          AspectRatio\  \[Rule] \ Automatic]; \[IndentingNewLine]eded = 
          duoedges[\ e2, duos]; \[IndentingNewLine]lines = 
          Table[Line[duop[\([\ eded[\([i]\)]\ ]\)]\ ], {i, 
              Length[eded]}]; \[IndentingNewLine]s4 = \ 
          Show[Graphics[{\[IndentingNewLine]Orange, \[IndentingNewLine]lines\
\[IndentingNewLine]}], \ 
            AspectRatio\  \[Rule] \ Automatic]; \[IndentingNewLine]Show[s2, 
          s3, s4, \ 
          Graphics[{{PointSize[ .05]}, Point\  /@ \ duop, \ 
              teksten}], \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- \
.1\), 1.1}, {\(- .1\), 1.1}}, \ 
          AspectRatio\  \[Rule] \ Automatic];\[IndentingNewLine]]\)], "Input",\

  InitializationCell->True],

Cell[BoxData[
    \(randon2ddelaunay[10]\)], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["k-th order Delauany triangulations.", "Subsection"],

Cell[TextData[{
  "We can easily generalize these programs. First of all, note that the \
procedure",
  StyleBox[" ",
    FontWeight->"Bold"],
  StyleBox["duoedges[ eandr_, duos_] ",
    FontFamily->"Courier",
    FontWeight->"Bold"],
  "works not only for 2nd order Voronoi regions, but for k-th order Voronoi \
regions, for k=1,...,n-1. Here `works` means that ",
  StyleBox["duoedges[ eandr_, duos_]",
    FontFamily->"Courier",
    FontWeight->"Bold"],
  "\n returns an edge exactly for those pairs of k-th order Voronoi regions \
that have a k-th -order Voronoi edge in common in Vk(S)."
}], "Text"],

Cell[BoxData[
    \(order3delaunay[points_]\  := \ 
      Module[\[IndentingNewLine]{r, \ akv, \ e2, \ duos, \ bit, h, \ 
          teksten, \ s2, \ s3, \ s4}, \[IndentingNewLine]akv = 
          allkvoronoi[points]; \[IndentingNewLine]e3\  = \ 
          edgesandregions[akv, 3]; \[IndentingNewLine]kregions = 
          Union[Flatten[Table[e3[\([i, 3]\)], {i, Length[e3]}], 
              1]]; \[IndentingNewLine]h[duo_]\  := \ 
          Sum[points[\([\ duo[\([i]\)]\ ]\)], \ {i, 1, 3}]/
            3; \[IndentingNewLine]duop = 
          Table[\ h[kregions[\([i]\)]\ ], {i, \ 
              Length[kregions]}]; \[IndentingNewLine]bit = {0,  .05}; \
\[IndentingNewLine]teksten = 
          Table[\ Text[kregions[\([i]\)], \ duop[\([i]\)]\ ], \ {i, \ 
              Length[kregions]}]; \[IndentingNewLine]s2 = 
          showkvoronoiinf[points, 3]; \[IndentingNewLine]s3\  = \ 
          PlanarGraphPlot[duop, 
            LabelPoints\  \[Rule] \ False]; \[IndentingNewLine]Show[s2, \ 
          Graphics[{{PointSize[ .03]}, Point\  /@ \ duop, \ teksten}], 
          s3, \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- .1\), 
                1.1}, {\(- .1\), 1.1}}, \ 
          AspectRatio\  \[Rule] \ Automatic]; \[IndentingNewLine]eded = 
          duoedges[\ e3, kregions]; \[IndentingNewLine]lines = 
          Table[Line[duop[\([\ eded[\([i]\)]\ ]\)]\ ], {i, 
              Length[eded]}]; \[IndentingNewLine]s4 = \ 
          Show[Graphics[{\[IndentingNewLine]Orange, \[IndentingNewLine]lines\
\[IndentingNewLine]}], \ 
            AspectRatio\  \[Rule] \ Automatic]; \[IndentingNewLine]Show[s2, 
          s3, s4, \ 
          Graphics[{{PointSize[ .05]}, Point\  /@ \ duop, \ 
              teksten}], \ \[IndentingNewLine]PlotRange\  \[Rule] \ {{\(- \
.1\), 1.1}, {\(- .1\), 1.1}}, \ 
          AspectRatio\  \[Rule] \ Automatic];\[IndentingNewLine]]\)], "Input"],

Cell[TextData[{
  StyleBox["order3delaunay[points]",
    FontFamily->"Courier",
    FontWeight->"Bold"],
  " returns the Delaunay triangulation of the barycenters of the triples of \
points that occur as 3th order Voronoi region and the graph D3(S)."
}], "Text"],

Cell[BoxData[
    \(order3delaunay[points]\)], "Input"],

Cell[TextData[{
  StyleBox["orderdelaunay[points_,k_]",
    FontFamily->"Courier",
    FontWeight->"Bold"],
  " returns the graph Dk(S)."
}], "Text"],

Cell[BoxData[
    \(orderdelaunay[points_, k_]\  := \ 
      Module[\[IndentingNewLine]{n, box, \ akv, erinfo, kregions, \ 
          nrkpoints, bar, \ kpoints, \ eded, \ 
          lines\ }, \[IndentingNewLine]n\  = \ 
          Length[points]; \[IndentingNewLine]box = {{\(- .1\), 
              1.1}, {\(- .1\), 1.1}}; \[IndentingNewLine]akv = 
          allkvoronoi[points]; \[IndentingNewLine]erinfo\  = \ 
          edgesandregions[akv, k]; \[IndentingNewLine]kregions = 
          Union[Flatten[Table[erinfo[\([i, 3]\)], {i, Length[erinfo]}], 
              1]]; \[IndentingNewLine]nrkpoints\  = \ 
          Length[kregions]; \[IndentingNewLine]bar[tu_]\  := \ 
          Sum[points[\([\ tu[\([i]\)]\ ]\)], \ {i, 1, k}]/
            k; \[IndentingNewLine]kpoints = 
          Table[\ bar[kregions[\([i]\)]\ ], {i, \ 
              nrkpoints}]; \[IndentingNewLine]eded = 
          duoedges[\ erinfo, kregions]; \[IndentingNewLine]lines = 
          Table[Line[kpoints[\([\ eded[\([i]\)]\ ]\)]\ ], {i, 
              Length[eded]}]; \[IndentingNewLine]Show[
          Graphics[{\[IndentingNewLine]Hue[
                N[1 - k/n]], \[IndentingNewLine]lines\[IndentingNewLine]}], \ 
          AspectRatio\  \[Rule] \ Automatic]\[IndentingNewLine]]\)], "Input",
  InitializationCell->True],

Cell["\<\
We run it on our example configuration, and plot Dk(S) together \
with Vk(S) for each k in 1,...,n-1, where n = |S|.\
\>", "Text"],

Cell[BoxData[
    \(\(\(For[kk = 1, \ 
      kk\  \[LessEqual] \ 
        Length[points]\  - 1, \[IndentingNewLine]Print["\<ORDER \>", 
        kk]; \[IndentingNewLine]t1 = \ 
        showkvoronoiinf[points, \ kk]; \[IndentingNewLine]t2 = 
        orderdelaunay[points, kk]; \[IndentingNewLine]Show[{t1, t2}, \ 
        PlotRange\  \[Rule] \ {{\(- .1\), 1.1}, {\(- .1\), 
              1.1}}]; \[IndentingNewLine]\(kk++\);\[IndentingNewLine]]\)\(\
\[IndentingNewLine]\)\)\)], "Input"],

Cell["\<\
We can do the same on arbitrary examples if n doesn't exceed, say, \
n=25.\
\>", "Text"],

Cell[BoxData[{
    \(\(n = 25;\)\), "\[IndentingNewLine]", 
    \(\(points\  = \ Table[{Random[], Random[]}, {i, 25}];\)\)}], "Input"],

Cell[BoxData[
    \(\(\(For[kk = 1, \ 
      kk\  \[LessEqual] \ 
        Length[points]\  - 1, \[IndentingNewLine]Print["\<ORDER \>", 
        kk]; \[IndentingNewLine]t1 = \ 
        showkvoronoiinf[points, \ kk]; \[IndentingNewLine]t2 = 
        orderdelaunay[points, kk]; \[IndentingNewLine]Show[{t1, t2}, \ 
        PlotRange\  \[Rule] \ {{\(- .1\), 1.1}, {\(- .1\), 
              1.1}}]; \[IndentingNewLine]\(kk++\);\[IndentingNewLine]]\)\(\
\[IndentingNewLine]\)\)\)], "Input"]
}, Closed]]
}, Open  ]]
},
FrontEndVersion->"4.0 for X",
ScreenRectangle->{{0, 1280}, {0, 1024}},
AutoGeneratedPackage->None,
WindowSize->{753, 974},
WindowMargins->{{Automatic, 6}, {Automatic, 4}},
PrintingPageRange->{Automatic, Automatic},
PrintingOptions->{"PaperSize"->{597.562, 842.375},
"PaperOrientation"->"Portrait",
"PostScriptOutputFile":>FrontEnd`FileName[{$RootDirectory, "user1", "home6", \
"Lindenbergh"}, "delaunay3.nb.ps", CharacterEncoding -> "ISO8859-1"],
"Magnification"->1}
]


(***********************************************************************
Cached data follows.  If you edit this Notebook file directly, not using
Mathematica, you must remove the line containing CacheID at the top of 
the file.  The cache data will then be recreated when you save this file 
from within Mathematica.
***********************************************************************)

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{

Cell[CellGroupData[{
Cell[1739, 51, 82, 1, 63, "Section"],

Cell[CellGroupData[{
Cell[1846, 56, 40, 0, 45, "Subsection"],
Cell[1889, 58, 69, 0, 32, "Text"],
Cell[1961, 60, 18763, 650, 9627, "Input",
  InitializationCell->True]
}, Closed]],

Cell[CellGroupData[{
Cell[20761, 715, 59, 0, 29, "Subsection"],
Cell[20823, 717, 77, 1, 32, "Text"],
Cell[20903, 720, 473, 9, 75, "Input"],
Cell[21379, 731, 230, 5, 50, "Text"],
Cell[21612, 738, 267, 5, 59, "Input"],
Cell[21882, 745, 232, 5, 50, "Text"],
Cell[22117, 752, 262, 6, 43, "Input"],
Cell[22382, 760, 194, 5, 50, "Text"],
Cell[22579, 767, 687, 14, 107, "Input"],
Cell[23269, 783, 260, 5, 50, "Text"],
Cell[23532, 790, 315, 6, 50, "Text"],
Cell[23850, 798, 1559, 27, 411, "Input",
  InitializationCell->True],
Cell[25412, 827, 101, 1, 32, "Text"],
Cell[25516, 830, 86, 2, 27, "Input"],
Cell[25605, 834, 157, 4, 27, "Input"],
Cell[25765, 840, 207, 4, 75, "Input"],
Cell[25975, 846, 181, 4, 32, "Text"],
Cell[26159, 852, 305, 6, 43, "Input"],
Cell[26467, 860, 157, 3, 32, "Text"],
Cell[26627, 865, 2037, 37, 395, "Input",
  InitializationCell->True],
Cell[28667, 904, 69, 0, 32, "Text"],
Cell[28739, 906, 77, 1, 32, "Text"],
Cell[28819, 909, 473, 9, 75, "Input"],
Cell[29295, 920, 230, 5, 50, "Text"],
Cell[29528, 927, 267, 5, 59, "Input"],
Cell[29798, 934, 232, 5, 50, "Text"],
Cell[30033, 941, 262, 6, 43, "Input"],
Cell[30298, 949, 194, 5, 50, "Text"],
Cell[30495, 956, 687, 14, 107, "Input"],
Cell[31185, 972, 260, 5, 50, "Text"],
Cell[31448, 979, 315, 6, 50, "Text"],
Cell[31766, 987, 1559, 27, 411, "Input",
  InitializationCell->True],
Cell[33328, 1016, 101, 1, 32, "Text"],
Cell[33432, 1019, 86, 2, 27, "Input"],
Cell[33521, 1023, 157, 4, 27, "Input"],
Cell[33681, 1029, 207, 4, 75, "Input"],
Cell[33891, 1035, 181, 4, 32, "Text"],
Cell[34075, 1041, 305, 6, 43, "Input"],
Cell[34383, 1049, 157, 3, 32, "Text"],
Cell[34543, 1054, 2037, 37, 395, "Input",
  InitializationCell->True],
Cell[36583, 1093, 53, 1, 27, "Input"]
}, Closed]],

Cell[CellGroupData[{
Cell[36673, 1099, 57, 0, 29, "Subsection"],
Cell[36733, 1101, 604, 15, 68, "Text"],
Cell[37340, 1118, 1878, 32, 363, "Input"],
Cell[39221, 1152, 262, 6, 50, "Text"],
Cell[39486, 1160, 55, 1, 27, "Input"],
Cell[39544, 1163, 149, 5, 32, "Text"],
Cell[39696, 1170, 1293, 22, 283, "Input",
  InitializationCell->True],
Cell[40992, 1194, 140, 3, 32, "Text"],
Cell[41135, 1199, 484, 9, 139, "Input"],
Cell[41622, 1210, 98, 3, 32, "Text"],
Cell[41723, 1215, 134, 2, 43, "Input"],
Cell[41860, 1219, 484, 9, 139, "Input"]
}, Closed]]
}, Open  ]]
}
]
*)




(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)