Minimization problem and Voronoi Mesh

Question

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,
  AllTrue[
   Table[Not[RegionEqual[
     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?

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[
   Table[{ReplacePart[
      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]]}],
     ReplacePart[
      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, 
      n}]
    },
   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}]];

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

Any comments or questions are welcome.