I have a random convex mesh $Q$, made of $n$ polygons, and I want to test how close it is to a Voronoi tessellation. In other words, I'm looking for generators (seeds) $\{(x_i,y_i)\}_{1\leq i\leq n}$.

This requires a bit of math, therefore I'll present some of the ideas in Chapter 2.6 of Spatial Tessellations, by A. Shewhart and S. Wilks, so that you have a bit of context.

In order to check if the mesh forms a Voronoi Tessellation, we need to guarantee that

  1. Each generator is in the associated Voronoi region.
  2. An edge in $Q$ should be on the perpendicular bisector of the two side generators.

Regarding 1, let $e$ be an edge of $Q$ shared by two polygons $q_i$ and $q_j$. Then, for some $a,b\in\mathbb{R}$, the line containing $e$ can be expressed by the equation $$ ax+by=1 $$ Suppose that $q_i$ and the origin lie in the same side of $e$. Then we get $$ ax_i+by_i>0\,\text{ and }\,ax_j+by_j<0. $$ Collecting the inequalities for all edges we get a system of linear inequalities, denoted by $$ A\mathbf{x}>0. $$ For 2, the line containing $e$ should contain the midpoint of $(x_i,y_i)$ and $(x_j,y_j)$. Hence we get $$ a\frac{x_i+x_j}{2}+b\frac{y_i+y_j}{2}=1. $$ Furthermore, since the line connecting $(x_i, y_i)$ and $(x_j, y_j)$ should be perpendicular to $e$, we get $$ a(y_i-y_j)-b(x_i-x_j)=0. $$ We get similar equations for each edge. Collecting them all, we obtain a system of linear equations, which we denote by $$ B\mathbf{x}=\mathbf{c}. $$ Now, instead of searching for the exact solution of the previous equation, I merely want to introduce a certain error factor, in order to characterise the "closedness" to a Voronoi tessellation. Therefore, together with $A\mathbf{x}>0$, I want to use Mathematica to solve the problem $$ \min_{\mathbf{x}}\| B\mathbf{x}-\mathbf{c} \|^2. $$ How do I do this?

My main problem is in defining the equations for the corresponding seeds. Regarding $A\mathbf{x}>0$, I could simply use RegionIntersection with each polygon to simply force the points to be inside them. But how do I define $B$ and $\mathbf{c}$? It doesn't seem obvious.

In the end, I want something like

Minimize[{Dot[B, {Join[Table[x[i], {i, n}], Table[y[i], {i, n}]]}] - c,
  Dot[A, {Join[Table[x[i], {i, n}], Table[y[i], {i, n}]]}] > 0},
 {Join[Table[x[i], {i, n}], Table[y[i], {i, n}]]}]

where the condition $A\mathbf{x}>0$ could be replaced in the following manner

Minimize[{Dot[B, {Join[Table[x[i], {i, n}], Table[y[i], {i, n}]]}] - c,
     RegionIntersection[MeshPrimitives[Q, 2][[i]], 
      Point[{x[[i]], y[[i]]}]], EmptyRegion[2]]], {i, n}], TrueQ]},
 {Join[Table[x[i], {i, n}], Table[y[i], {i, n}]]}]

For $B$ and $\mathbf{c}$, I have access to the edges that share a polygon, but how do I make the correct association with each $(x_i,y_i)$?

Any ideas?


1 Answer 1


This is my solution

abf = Function[l, Module[{x1, y1, x2, y2},
    x1 = l[[1, 1, 1]];
    y1 = l[[1, 1, 2]];
    x2 = l[[1, 2, 1]];
    y2 = l[[1, 2, 2]];
    Solve[as x1 + bs y1 == 1 && as x2 + bs y2 == 1, {as, bs}]

n = MeshPrimitives[mesh, 2] // Length;
shre0 = Complement[MeshPrimitives[mesh, 1], 
   MeshPrimitives[BoundaryMesh[mesh], 1]];
edgn = Length[shre0];
cents = RegionCentroid[MeshPrimitives[mesh, 2]];
shre = Table[Line[SortBy[shre0[[i, 1]], Norm]], {i, edgn}];
pol = Table[
   Append[MeshPrimitives[mesh, 2][[i]][[1]], 
    MeshPrimitives[mesh, 2][[i]][[1, 1]]], {i, n}];
polin0 = Table[
   Table[Line[{pol[[j, i]], pol[[j, i + 1]]}], {i, 
     Length[pol[[j]]] - 1}], {j, n}];
polin = Table[
   Table[Line[SortBy[polin0[[j, i, 1]], Norm]], {i, 
     Length[polin0[[j]]]}], {j, n}];
lsp = Table[Intersection[shre, polin[[i]]], {i, Length[polin]}];
lor = {};
For[i = 1, i <= n, i++,
  For[j = i + 1, j <= n, j++,
   regg = Intersection[lsp[[i]], lsp[[j]]];
   If[regg =!= {},
    lor = 
     Append[lor, {{i, j}, {abf[regg[[1]]][[1, 1, 2]], 
        abf[regg[[1]]][[1, 2, 2]]}, regg[[1]]}]]
abi = Transpose[lor][[1]];
ab = Transpose[lor][[2]];
matB = Flatten[
      ConstantArray[0, 2 n], {2*abi[[i, 1]] - 1 -> ab[[i, 1]], 
       2*abi[[i, 1]] -> ab[[i, 2]], 2*abi[[i, 2]] - 1 -> ab[[i, 1]], 
       2*abi[[i, 2]] -> ab[[i, 2]]}],
      ConstantArray[0, 2 n], {2*abi[[i, 1]] - 1 -> -ab[[i, 2]], 
       2*abi[[i, 1]] -> ab[[i, 1]], 2*abi[[i, 2]] - 1 -> ab[[i, 2]], 
       2*abi[[i, 2]] -> -ab[[i, 1]]}]},
    {i, 1, edgn}], 1];
cB = Flatten[ConstantArray[{2, 0}, edgn]];

crns = Transpose[
   Table[MeshPrimitives[mesh, 0][[i, 1]], {i, 
     Length[MeshPrimitives[mesh, 0]]}]];
cr00 = Min[crns[[1]]];
cr10 = Max[crns[[1]]];
cr01 = Min[crns[[2]]];
cr11 = Max[crns[[2]]];

minV = FindMinimum[{Norm[
      Dot[matB, Flatten[Table[{xs[i], ys[i]}, {i, n}]]] - cB]^2,
    Table[cr00 <= xs[i] <= cr10 && cr01 <= ys[i] <= cr11 &&

      RegionMember[MeshPrimitives[mesh, 2][[i]], {xs[i], ys[i]}], {i, 
   Flatten[Table[{xs[i], ys[i]}, {i, n}]]];
error = minV[[1]];
ptt = Table[{minV[[2, 2 i - 1, 2]], minV[[2, 2 i, 2]]}, {i, n}];

Then, testing for a Voronoi mesh we get

n = 5; mesh = VoronoiMesh[RandomReal[1, {n, 2}]];

enter image description here

Just as an example, for a random convex mesh, we get

enter image description here

Any comments or questions are welcome.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.