The programs.

<< mima.m

circlecenter[xyz_List] := 
      Module[
  
{a, d, 
          e}, 
a = ca[xyz]; 
If[
          a == 
            0, Print["Dividing by zero"];]; 
d = 
          cd[xyz]; 
e = 
          ce[xyz]; 
{Limit[-d/(2  a),  
            t  ->  0], 
          Limit[-e/(2  a),  
            t  ->  0]}
]

unbedges[ alledges_]  :=  
      Module[
{nrpos,  i,  
          uedges}, 
uedges = {}; 
nrpos  =  Table[Length[ Position[alledges, alledges[[i]], 1]], {i, 1,  
              Length[alledges]}]; 
For[i = 1,  
          i  ≤  Length[ alledges], 
If[ 
            nrpos[[i]]  ==  
              1, 
uedges  =  
                Append[ uedges,  
                  alledges[[i]]];
]; 
i++
]; 
uedges
]

ctoe[ 
        c_]  := 
{Sort[{c[[1]], 
            c[[2]]}], 
Sort[{c[[2]], 
            c[[3]]}], 
Sort[{c[[3]], 
            c[[1]]}]}

checkorient[pp_]  :=  
      Module[
{ld,  s}, 
ld  =  
          linedet[pp]; 
s  =  
          psign[ld]; 
If[ s  ==  0, 
          Return[False] ]; 
If[s  ==  1,  
          Return[ -1] ]; 
1
]

\(emptycircles[polpp_List] := 
      Module[
 (*\ 
          gets\ all\ voronoi\ circles\ for\ a\ list\ of\ polynomial\ sites\ *) 
{l, rr, lc, circlelist, thiscircle, i, cp, \ 
          changeor, \ ocircle, \ otherpoints, j, \ empty, \ h, \ 
          tocheck}, 
l = 
          Length[polpp]; 
rr = 
          Range[l]; 

 (*\ 
          lc = \ #  of\ different\ circles\ *) 
lc = 
          Binomial[l, 
            3]; 
circlelist = {}; 
thiscircle = {1, 2, 
            3}; 

For[i = 1, 
          i ≤ lc, 

 (*\ 
            orient\ the\ circle\ clockwise\ *) 
cp\  = \ 
            polpp[\([\ thiscircle\ ]\)]; 
changeor\  = \ 
            checkorient[\ cp\ ]; 
If[\ 
            changeor\  == \ \(-1\), \ 
ocircle\  = \ 
              thiscircle[\([\ {1, 3, 
                    2}\ ]\)], 
ocircle\  = \ 
              thiscircle
]; 

cp\  = \ 
            polpp[\([\ ocircle\ ]\)]; 
otherpoints = 
            Complement[rr, 
              ocircle]; 
 (*\ \ \(Print[
                ocircle, \ \ otherpoints];\)\ \ *) 
j\  = 
            1; 
empty\  = \ 
            True; 

While[\((\ 
              j ≤ \((l - 3)\)\  && \ 
                empty\  == \ True)\), 
h = 
              otherpoints[\([j]\)]; 
tocheck\  = \ 
              Append[\ ocircle, h]; 
If[\ 
              psign[\ circledet[\ polpp[\([\ tocheck\ ]\)]\ ]\ ]\  == \ 
                1, 
\(j++\), 
empty\  = \ 
                False
];
]; 

If[\ 
            empty\  == \ 
              True, \ 
\(circlelist\ \  = \ 
                Append[\ circlelist, 
                  ocircle\ ];\)
]; 

 (*\ 
            Use\ the\ undocumented\ function\ `NextKSubset'\ *) 
 (*\ \(Print["thiscircle ", \ 
                thiscircle];\)\ *) 
thiscircle = 
            NextKSubset[rr, 
              thiscircle]; 
\(i++\);
]; 
circlelist
]\)

firstcoefpoly[function_] := 
      Module[
  
{depth = 0, coef = 0}, 
If[
          function == 0, Return[0]]; 
While[
          coef == 0, depth++; 
coef = 
            Coefficient[function, t, 
              depth];]; 
coef]

psign[exp_] := \(Module[
\ \ {}, 
If[\ NumberQ[exp], \ 
          Return[\ Sign[exp]\ ]
]; 
If[
          PolynomialQ[exp], 
          Return[\ Sign[\ 
              firstcoefpoly[
                exp]]]]; 
False
		]\)

n2[p_]  :=  p[[1]]^2  +  p[[2]]^2

linedet[pp_]  :=  
      Module[
{mat, i}, 
mat  =  
          Table[Prepend[ pp[[i]], 1], {i, 3}]; 
Det[
          mat]
]

circledet[pp_]  :=  
      Module[
{mat, i}, 
mat  =  
          Table[Join[ pp[[i]],  {n2[  pp[[i]]  ], 1} ], {i, 1, 
              4}]; 
Det[mat]
]

bilineeq[{u_, v_}] := 
      Module[
{rmat, b, a, 
          pp}, 
rmat  =  {{0, -1}, {1, 
              0}}; 
b  =  (1/2) (u + 
                v); 
a  =  
          b  + rmat . (u - b); 
pp = {b, 
            a, {x, y}}; 
mat = 
          Table[Prepend[pp[[i]], 1], {i, 
              3}]; 
Det[mat]
]

Example: computing the type.

The sign of the matrix <linedet> is 
	0          p3 on line(p1,p2)
	1          p3 left of line(p1,p2)
          -1         p3 right  of line(p1,p2)

linedet[ {{1, 1}, {3, 3}, {2, 2}}]
linedet[{{1, 1}, {3, 3}, {0, 2}}]
linedet[{{1, 1}, {3, 3}, {2, 0}}]

The sign of the matrix <circledet> is 
	0	p4 on circle(p1,p2,p3)
	-1       p4 in clockwise oriented circle(p1,p2,p3)
	1          p4 out clockwise oriented circle(p1,p2,p3)

circledet[{{-1, 0}, {0, 1}, {1, 0}, {Cos[Pi/12], 
          Sin[Pi/12]}}]
circledet[{{-1, 0}, {0, 1}, {1, 
          0}, { .1,  .1}}]
circledet[{{-1, 0}, {0, 1}, {1, 0}, {3, 0}}]

The orientation of the circle is essential. Therefore the three points defining the circle, should be in the right cyclic order. This can be achieved as follows: pick two points p_i and p_j. These two points determine a directed line l_ij. Determine s = signum[ linedet[{pi,pj,pk}] ]. If

	 s = -1            pi, pj, pk
	 s = 1              pi, pk, pj
	 s = 0              no general position

We define some polynomial sites.

p3  =  {2 t,  2 t^3 +  t^4};
p8  =  {-2 t  -  t^3,  2 t};
p11 =   {-2 t^2  - 2 t^4,  -3 t  +  2 t^2 - 
            t^3 };
p16  =  {t^2 - 2  t^3 + 
            t^5, -t^4};
sites  =  {p3, p8, p11, p16};

We determine the orientation of the circle defined by p3, p8 and p11.

ld = linedet[{ p3,  p8,  p11}]
psign[ld]

checkorient[{ p3,  p8,  p11}]

In this case, -1 means False, so the right orientation is p3, p11, p8. Now we can check if p16 in circle(p3, p11, p8) for small t.

cdet = circledet[{p3, p11, p8, p16}]
psign[cdet]

Indeed, p16 is in the circle(p3, p11, p8). The program emptycircles[ polpp_List ] returns all empty circles of a list S of polynomial sites. By definition, this is the type of the Voronoi diagram V(S).

emptycircles[ sites]

Let's try a bigger example.

punten = {{0, (t^5)}, {0, 2 t + (t^4) + 2 (t^5)}, {2 t, 
            2 (t^3) + (t^4)}, {(t^2), 
            2 t + 2 (t^2) - (t^3)}, {-2 (t^4), -(t^5)}, {3 (t^5), 
            (t^2) + 2 (t^3) + 3 (t^5)}, {2 t + 3 (t^2) - 2 (t^3), 
            2 t - 2 (t^2) - (t^3)}, {-2 t - (t^3), 
            2 t}, {-2 t + (t^3), -2 t - (t^3)}, {t + 2 (t^2) + 
              2 (t^3), -t + 3 (t^3)}, {-2 (t^2) - 
              2 (t^4), -3 t + 2 (t^2) - (t^3) + (t^4)}, {-2 t - (t^4), 
            3 (t^2) + 3 (t^3) - (t^5)}, {-2 (t^2) - (t^4), -2 (t^2) + 
              3 (t^3)}, {-t + 3 (t^2) - 2 (t^3) + 3 (t^4), 
            t + 3 (t^2) - 3 (t^5)}, {-t - 2 (t^5), -3 (t^3) + 
              2 (t^5)}, {(t^2) - 2 (t^3) + (t^5), -(t^4)}, {3 (t^3) + 
              2 (t^5), -(t^4)}, {3 (t^3) + 3 (t^4) + 2 (t^5), 
            2 t - 2 (t^3) - (t^4)}, {(t^2) + (t^3) + 3 (t^5), 
            0}, {3 (t^2) + 2 (t^3) + 3 (t^5), 
            (t^2) + 3 (t^3) + 2 (t^4)}};

type20 = emptycircles[punten ]

From every Voronoi vertex (= circle centre), three edges leave.

ctoe[{1, 5, 6}]

alledges = 
      Flatten[Table[ctoe[type20[[i]]], {i, Length[type20]}], 
        1]

From this edge list, we can easily distill the unbounded edges.

Table[
      Length[ Position[alledges, alledges[[i]], 1]], {i, 1,  
        Length[alledges]}]

unbedges[ alledges]

These are edges between polynomial sites on the `combinatorial convex hull'.

We can compute the actual centers of the Voronoi circles.

cccc = 
      Table[circlecenter[ punten[[ type20[[i]] ]] ], {i, 
          Length[type20]}]

By orienting the combinatorial convex hull in the rigth way and by analysing the (limit) direction, this gives us all information we need in order to construct the one-skeleton + multiplicities.

Example: computing a circle center.

u =  {0, -t};
v = {-t, 0};
w  =  {t - t^2, -2 t};

Computing the circle center, using the `circular' determinants.

circlecenter[{u, v, w}]

Alternatively, we construct equations in t for the bisectors

eq1 = bilineeq[{u, v}]
eq1  =  Expand[eq1/t]
sol = Solve[eq1 == 0, y]
s1  =  sol[[1, 1, 2]]

eq2 = bilineeq[{u, w}]
eq2 = Expand[eq2/t]
sol  =  Solve[eq2 == 0, y]
s2 = sol[[1, 1, 2]]

And intersect them for finding the circle center.

sol2 = Solve[s1 == s2, x]
s3  =  sol2[[1, 1, 2]]

sx  =  s3  /.  { t  ->  0}
sy  =  s1  /.  {x  ->  sx}