# Minimization problem and Voronoi Mesh

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?

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.