Higher order Delauany triangulations.

k-th Voronoi code.

Initialize  the programs fromn the  kvoronoi package.

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

ca[cl_] := Module[
		(* 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[
			(* computes the centre of the circle, defined by the points <xyz> *)
		{a,d,e},
		a = ca[xyz];
		If[ a==0,
			Print["Dividing by zero"];
		];
		d = cd[xyz];
		e = ce[xyz];
		{-d/(2a),-e/(2a)}
	]

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

allkvoronoi[points_List] := Module[
		(* gets all voronoi circles and the points they contain *)
		{l, rr, lc, circlelist, thiscircle, i, center, radius, otherpoints, k,h, j},
    
		l = Length[points];
		rr = Range[l];
		(* lc =  #of different circles *)
		lc = Binomial[l,3];
		circlelist = {};
		thiscircle = {1,2,3};
		(* Use the undocumented function `NextKSubset' *)
		For[ i=1 , i <= lc,
			center = circlecenter[ points[[thiscircle]] ];
			radius = distance[{points[[thiscircle]][[1]], center}];
			otherpoints = Complement[ rr, thiscircle];
			k = {};
			For[ j=1, j <= (l-3),
				h = otherpoints[[j]];
				If[ distance[{points[[h]],center}] < radius,
						k=Append[k,h];
					];
				j++;
			];
			circlelist = Append[circlelist, {i, thiscircle, N[center], k}];
			thiscircle = NextKSubset[rr,thiscircle];
			i++;
	];
	circlelist
]	

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

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

buildedges[ek_] := Module[
		(* build all edges from the incident faces and edges info in <ek> *)
		{cek, nrv, edges, i, edge, faces, ii, jj, al},
		cek = ek;
		nrv = Length[ cek ];
		edges = {};
		For[ i=1, i < nrv,
			For[ j=1, j<= Length[ cek[[i,2]] ],
				edge = cek[[i,2]][[j]];
				faces = cek[[i,3]];
				For[ ii=i+1, ii <= nrv,
					al = Length[ cek[[ii,2]] ];
					For[ jj=1, jj <= al,
						If[ cek[[ii,2]][[ jj]] == edge,
							If[ Length[Intersection[ faces, cek[[ii,3]] ] ] == 2,
								edges = Append[ edges, {cek[[i,1]], cek[[ii,1]] }];
								cek = Delete[ cek, {ii,2,jj} ];
								jj = 4;
							  ii = nrv+1;
							];
						];
						jj++;
					];
					ii++;
				];
			 j++;
		];
	i++;
];
edges
]		
		

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

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

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

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

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

voronoiflist[points_] := Module[
		(* Returns all Voronoi faces defined by <points> *)
		{n, akv, allfaces, i, faces, ear, j},
		n = Length[points];
		akv = allkvoronoi[points];
		allfaces = {{Range[n]}};
		For[ i=n-1, i >= 1,
			faces = {};
			ear = edgesandregions[akv,i];
			For[j=1, j <= Length[ear],
				faces = Join[faces, ear[[j,3]]];
				j++;
			];
			allfaces = Append[allfaces, Union[faces]];
			i--;
		];
		Append[allfaces,{{}}]
		]

randomnvor[n_Integer] := Module[
		(* computes the euler chararacteristic for a set of <n> random points *)
		{r, mp, flist, fvector, avector, euler},
		r := Random[Integer,{0,1000}];
		mp = Table[{r,r},{i,n}];
		Print["points = ",mp];
		flist = voronoiflist[mp];
		fvector=FVector[flist];
		Print["fvector = ", fvector];
		avector=Table[(-1)^(i),{i,Length[fvector]}];
	   euler = avector.fvector;
		Print["ϰ = ", euler];
		euler
		]

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

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

finfvector[points_] := Module[
		{akv, finfvec, k, edreg, occur, vertices, lv, valencies, nrunbound, j},
		akv = allkvoronoi[ points ];
		n= Length[points];
		finfvec = {};
		For[k=1, k <n,
			edreg = edgesandregions[akv,k];
			occur =Sort[Flatten[ buildedges[ edreg ] ] ];
			vertices = Union[occur];
			lv = Length[vertices];
			valencies = Table[Count[occur,vertices[[i]] ],{i,lv}];
			nrunbound = 0;
			For[ j = 1, j <= lv,
				nrunbound = nrunbound + (3-valencies[[j]]);
				j++;
			];
			(* added later, RCL *)
			If[ lv == 0,
			  nrunbound = 3;
			];
			finfvec = Append[finfvec, nrunbound];
			k++;
		];
		finfvec
		]		

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

evector[points_] := Module[
		(* 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 *)
		{akv, finfvec, k, edreg, occur, vertices, lv, valencies, nrunbound, j},
		akv = allkvoronoi[ points ];
		n= Length[points];
		evec  = {};
		finfvec = {};
		For[k=1, k <n,
			edreg = edgesandregions[akv,k];
			bedges =  buildedges[ edreg ];
			occur =Sort[Flatten[ bedges ] ];
			vertices = Union[occur];
			lv = Length[vertices];
			valencies = Table[Count[occur,vertices[[i]] ],{i,lv}];
			nrunbound = 0;
			For[ j = 1, j <= lv,
				nrunbound = nrunbound + (3-valencies[[j]]);
				j++;
			];
			(* added later, RCL *)
			If[ lv == 0,
			  nrunbound = 3;
			];
			finfvec = Append[finfvec, nrunbound+Length[bedges]];
			k++;
		];
		finfvec
		]		


ffpoly[points_] := Module[
		{rfvec, poly, i},
		rfvec = Reverse[fvector[points]];
		poly = 0;
		For[ i =1, i <= Length[rfvec],
			poly = poly + rfvec[[i]] t^(i-1);
			i++;
		];
		Factor[poly]
	]	

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

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[
		{l,i,j},
		l = Length[vlist];
		Table[{Min[Table[vlist[[i,j]],{i,l}]],Max[Table[vlist[[i,j]],{i,l}]]},
    	{j,Length[vlist[[1]]]}]
]

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

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

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

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


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

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

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

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

buildedgesinf[ek_] := Module[
		(* build all edges from the incident faces and edges info in <ek> *)
		{cek, nrv, edges, i,j, edge, faces, ii, jj, al},
		cek = ek;
		nrv = Length[ cek ];
		edges = {};
		infedges={};
		For[ i=1, i < nrv,
			For[ j=1, j<= Length[ cek[[i,2]] ],
				edge = cek[[i,2]][[j]];
				faces = cek[[i,3]];
				For[ ii=i+1, ii <= nrv,
					al = Length[ cek[[ii,2]] ];
					For[ jj=1, jj <= al,
						If[ cek[[ii,2]][[ jj]] == edge,
							If[ Length[Intersection[ faces, cek[[ii,3]] ] ] == 2,
								edges = Append[ edges, {cek[[i,1]], cek[[ii,1]] }];
								cek = Delete[ cek, {ii,2,jj} ];
								jj = 5;
							  ii = nrv+1;
							];
						];
						jj++;
					];
					ii++;
				];
				If[ jj < 5, infedges=Append[infedges,{i,edge}]];
			 j++;
		];
	i++;
];
For[k=1, k <= Length[cek[[nrv,2]]],
			infedges = Append[infedges, {i,cek[[nrv,2,k]]} ];
			k++;
];
{edges,infedges}
]		

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

showkvoronoiinf[points_, k_Integer] := Module[
		{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[{}]];
		lines = {};
		n = Length[ points ];
		If[ k >= n, Return[{}]];
	  akv = allkvoronoi[ points ];
		edreg = edgesandregions[ akv, k];
		edgeinfo = buildedgesinf[edreg];
		edgerel = edgeinfo[[1]];
		(* Add the edges *)
		For[ i = 1, i <= Length[ edgerel],
				b  = edgerel[[i,1]];
				e = edgerel[[i,2]];
				lines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
			  i++;
		];
		(* Now we add the halflines *)
	   slopes = allslopes[points];
		infinfo = edgeinfo[[2]];
		For[i=1, i <= Length[infinfo],
			infedge = infinfo[[i]];
			thiscircle=edreg[[infedge[[1]],1]] ;
			(* determine whether <thiscircle> is young or old *)
	  	hcircle = Length[ akv[[thiscircle,4]] ];
			If[ hcircle+2==k, cyoung = -1, cyoung = 1];
			slope= slopes[[pairtoposition[infedge[[2]],n ] ]];
			centre= akv[[ edreg[[infedge[[1]],1]],3]];
			radius = N[Sqrt[ distance[ {centre, points[[infedge[[2,1]] ]]} ] ] ];
			ss = (2.radius)/(Sqrt[1+slope^2]){1,slope};
			s1 = ss+centre;
			p1label=infedge[[2,1]];
			p2label=infedge[[2,2]];
			p3label= Complement[akv[[thiscircle,2]],{p1label,p2label}][[1]];
			p1 = points[[p1label]];
			p2 = points[[p2label]];
			p3= points[[p3label]];
			testsign = sign[p1,p2, s1]  sign[ p1, p2, p3 ];
			If[ cyoung== testsign,
									s1 = -ss + centre;
			];
			lines = Append[lines,Line[{centre,s1}]];
			i++;
		];
		alllines = Join[{Thickness[.005], Hue[N[k/n]]}, lines];
		tpoints = Table[Text[i,points[[i]] ],{i,n}];
     Show[Graphics[{alllines,tpoints}], PlotRange -> All,
      AspectRatio -> Automatic]
	]	

showkvoronoiinf[points_, k_Integer,akv_ ]:= Module[
		{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},
		lines = {};
		n = Length[ points ];
		If[ k >= n, Return[{}]];
		edreg = edgesandregions[ akv, k];
		edgeinfo = buildedgesinf[edreg];
		edgerel = edgeinfo[[1]];
		(* Add the edges *)
		For[ i = 1, i <= Length[ edgerel],
				b  = edgerel[[i,1]];
				e = edgerel[[i,2]];
				lines = Append[ lines,Line[ {akv[[b,3]], akv[[e,3]]}]];
			  i++;
		];
		(* Now we add the halflines *)
	   slopes = allslopes[points];
		infinfo = edgeinfo[[2]];
		For[i=1, i <= Length[infinfo],
			infedge = infinfo[[i]];
			thiscircle=edreg[[infedge[[1]],1]] ;
			(* determine whether <thiscircle> is young or old *)
	  	hcircle = Length[ akv[[thiscircle,4]] ];
			If[ hcircle+2==k, cyoung = -1, cyoung = 1];
			slope= slopes[[pairtoposition[infedge[[2]],n ] ]];
			centre= akv[[ edreg[[infedge[[1]],1]],3]];
			radius = N[Sqrt[ distance[ {centre, points[[infedge[[2,1]] ]]} ] ] ];
			ss = (2.radius)/(Sqrt[1+slope^2]){1,slope};
			s1 = ss+centre;
			p1label=infedge[[2,1]];
			p2label=infedge[[2,2]];
			p3label= Complement[akv[[thiscircle,2]],{p1label,p2label}][[1]];
			p1 = points[[p1label]];
			p2 = points[[p2label]];
			p3= points[[p3label]];
			testsign = sign[p1,p2, s1]  sign[ p1, p2, p3 ];
			If[ cyoung== testsign,
									s1 = -ss + centre;
			];
			lines = Append[lines,Line[{centre,s1}]];
			i++;
		];
		alllines = Join[{Thickness[.005], Hue[N[k/n]]}, lines];
		tpoints = Table[Text[i,points[[i]] ],{i,n}];
     Show[Graphics[{alllines,tpoints}], PlotRange -> All,
      AspectRatio -> Automatic]
	]

Second order Delauany triangulations.

We start with a list S of 7 points.

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`}};

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).

akv = allkvoronoi[points];
e2  =  edgesandregions[akv, 2];
duos = 
      Union[Flatten[Table[e2[[i, 3]], {i, Length[e2]}], 1]]

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>.

h[duo_]  := ( 
          points[[duo[[1]]]] + points[[duo[[2]]]])/
        2
duop = 
        Table[ h[duos[[i]] ], {i,  Length[duos]}];

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).

bit = { .01,  .02};
teksten = 
        Table[ Text[duos[[i]],  duop[[i]] ],  {i,  
            Length[duos]}];
s2 = showkvoronoiinf[points, 2];
s3  =  
        PlanarGraphPlot[duop, 
          LabelPoints  ->  False];
Show[s2,  
        Graphics[{{PointSize[ .03]}, Point  /@  duop,  teksten}], 
        s3,  
PlotRange  ->  {{- .1, 
              1.1}, {- .1, 1.1}},  
        AspectRatio  ->  Automatic];

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.

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>.

duoedges[ eandr_,  duos_]  :=  
      Module[
{nv,  iedges,  edges,  i,  duo,  
          thisvertex,  pos,  j,  k}, 
nv  =  
          Length[eandr]; 
iedges = 
          Table[eandr[[i, 3]], {i, 
              nv}]; 
edges = {}; 
For[ 
          i = 1,  i  ≤  
            Length[duos], 
duo =  
            duos[[i]]; 
For[j = 1,  
            j  ≤  
              Length[iedges], 
thisvertex  =  
              iedges[[j]]; 
pos  =  
              Position[ thisvertex,  duo]; 
If[
              Length[pos]  >  0, 
pos = 
                pos[[1, 1]]; 
For[ k = 1,  
                k  ≤  3, 
If[ 
                  pos  ≠  k, 
edges  =  
                    Append[edges,  {i, Position[duos,  
                            thisvertex[[k]] ][[1, 
                            1]]}]
]; 
k++;
];
]; 
j++;
]; 
i++;
]; 
Union[
          Table[Union[edges[[i]]], {i, 
              Length[edges]}]]
]

We perform the program and draw all edges that it returns .

eded = duoedges[ e2, duos]

lines = 
        Table[Line[duop[[ eded[[i]] ]] ], {i, 
            Length[eded]}];

s4 =  
        Show[Graphics[{
Orange, 
lines
}],  AspectRatio  ->  Automatic];

We superimpose the result on the the Delaunay triangulation D(duop) that we have computed above. It's not the same triangulation.

Show[s2, s3, s4,  
        Graphics[{{PointSize[ .05]}, Point  /@  duop,  
            teksten}],  
PlotRange  ->  {{- .1, 
              1.1}, {- .1, 1.1}},  
        AspectRatio  ->  Automatic];

We group all commands. This gives a program that generates a random point set of <n> points and compute both triangulations for it.

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

Initialize  the programs fromn the  kvoronoi package.

We start with a list S of 7 points.

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`}};

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).

akv = allkvoronoi[points];
e2  =  edgesandregions[akv, 2];
duos = 
      Union[Flatten[Table[e2[[i, 3]], {i, Length[e2]}], 1]]

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>.

h[duo_]  := ( 
          points[[duo[[1]]]] + points[[duo[[2]]]])/
        2
duop = 
        Table[ h[duos[[i]] ], {i,  Length[duos]}];

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).

bit = { .01,  .02};
teksten = 
        Table[ Text[duos[[i]],  duop[[i]] ],  {i,  
            Length[duos]}];
s2 = showkvoronoiinf[points, 2];
s3  =  
        PlanarGraphPlot[duop, 
          LabelPoints  ->  False];
Show[s2,  
        Graphics[{{PointSize[ .03]}, Point  /@  duop,  teksten}], 
        s3,  
PlotRange  ->  {{- .1, 
              1.1}, {- .1, 1.1}},  
        AspectRatio  ->  Automatic];

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.

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>.

duoedges[ eandr_,  duos_]  :=  
      Module[
{nv,  iedges,  edges,  i,  duo,  
          thisvertex,  pos,  j,  k}, 
nv  =  
          Length[eandr]; 
iedges = 
          Table[eandr[[i, 3]], {i, 
              nv}]; 
edges = {}; 
For[ 
          i = 1,  i  ≤  
            Length[duos], 
duo =  
            duos[[i]]; 
For[j = 1,  
            j  ≤  
              Length[iedges], 
thisvertex  =  
              iedges[[j]]; 
pos  =  
              Position[ thisvertex,  duo]; 
If[
              Length[pos]  >  0, 
pos = 
                pos[[1, 1]]; 
For[ k = 1,  
                k  ≤  3, 
If[ 
                  pos  ≠  k, 
edges  =  
                    Append[edges,  {i, Position[duos,  
                            thisvertex[[k]] ][[1, 
                            1]]}]
]; 
k++;
];
]; 
j++;
]; 
i++;
]; 
Union[
          Table[Union[edges[[i]]], {i, 
              Length[edges]}]]
]

We perform the program and draw all edges that it returns .

eded = duoedges[ e2, duos]

lines = 
        Table[Line[duop[[ eded[[i]] ]] ], {i, 
            Length[eded]}];

s4 =  
        Show[Graphics[{
Orange, 
lines
}],  AspectRatio  ->  Automatic];

We superimpose the result on the the Delaunay triangulation D(duop) that we have computed above. It's not the same triangulation.

Show[s2, s3, s4,  
        Graphics[{{PointSize[ .05]}, Point  /@  duop,  
            teksten}],  
PlotRange  ->  {{- .1, 
              1.1}, {- .1, 1.1}},  
        AspectRatio  ->  Automatic];

We group all commands. This gives a program that generates a random point set of <n> points and compute both triangulations for it.

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

randon2ddelaunay[10]

k-th order Delauany triangulations.

We can easily generalize these programs. First of all, note that the procedure duoedges[ eandr_, duos_] works not only for 2nd order Voronoi regions, but for k-th order Voronoi regions, for k=1,...,n-1. Here `works` means that duoedges[ eandr_, duos_]
 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).

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

order3delaunay[points] returns the Delaunay triangulation of the barycenters of the triples of points that occur as 3th order Voronoi region and the graph D3(S).

order3delaunay[points]

orderdelaunay[points_,k_] returns the graph Dk(S).

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

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|.

\(\(\(For[kk = 1, \ 
      kk\  ≤ \ 
        Length[points]\  - 1, 
Print["ORDER ", 
        kk]; 
t1 = \ 
        showkvoronoiinf[points, \ kk]; 
t2 = 
        orderdelaunay[points, kk]; 
Show[{t1, t2}, \ 
        PlotRange\  -> \ {{\(- .1\), 1.1}, {\(- .1\), 
              1.1}}]; 
\(kk++\);
]\)\(
\)\)\)

We can do the same on arbitrary examples if n doesn't exceed, say, n=25.

n = 25;
points  =  Table[{Random[], Random[]}, {i, 25}];

\(\(\(For[kk = 1, \ 
      kk\  ≤ \ 
        Length[points]\  - 1, 
Print["ORDER ", 
        kk]; 
t1 = \ 
        showkvoronoiinf[points, \ kk]; 
t2 = 
        orderdelaunay[points, kk]; 
Show[{t1, t2}, \ 
        PlotRange\  -> \ {{\(- .1\), 1.1}, {\(- .1\), 
              1.1}}]; 
\(kk++\);
]\)\(
\)\)\)