# unbounded Voronoi iteration with finite and non-equal point radii

I would like to create a cool differential line growth animation based on this great tutorial.

The only missing part is the relax node, but I couldn't find any pointers on the web how to implement this.

Point Relax geometry node - it moves points with overlapping radii away from each other, optionally on a surface. The relaxation algorithm is a generalization of Lloyd’s algorithm, (Voronoi iteration), with support added for finite and non-equal point radii, as well as unbounded 2D & 3D.

Needs["NumericalDifferentialEquationAnalysis"]

periodicSpline[\[Xi]_] := Module[{}, a0 + a1 (\[Xi] - \[Xi]j) + a2 (\[Xi] - \[Xi]j)^2 + a3 (\[Xi] - \[Xi]j)^3]

periodicSplineCoeffCalc[pts_] :=
Module[{x, xDD, nCurve, coefs, qAhlberg, uAhlberg, sAhlberg,
tAhlberg, vAhlberg},
nCurve = Length[pts];
x = Join[pts, pts[[{1, 2}]]];
coefs = {};
xDD = ConstantArray[0, nCurve + 1];
qAhlberg = uAhlberg = sAhlberg = ConstantArray[0, nCurve + 1];
sAhlberg[[(0 + 1) + 0]] = 1;
tAhlberg = vAhlberg = ConstantArray[0, nCurve + 1];
tAhlberg[[(0 + 1) + nCurve]] = 1;
Do[
qAhlberg[[(0 + 1) + n]] = -(1/(4 + qAhlberg[[(0 + 1) + n - 1]]));
uAhlberg[[(0 + 1) + n]] =
qAhlberg[[(n + 1)]] (uAhlberg[[(0 + 1) + n - 1]] -
6 nCurve nCurve (x[[(0 + 1) + n + 1]] - 2 x[[(0 + 1) + n]] +
x[[(0 + 1) + n - 1]]));
sAhlberg[[(0 + 1) + n]] =
qAhlberg[[(n + 1)]] sAhlberg[[(0 + 1) + n - 1]];
, {n, 1, nCurve, 1}];
Do[
tAhlberg[[(0 + 1) + n]] =
qAhlberg[[(0 + 1) + n]] tAhlberg[[(0 + 1) + n + 1]] +
sAhlberg[[(0 + 1) + n]];
vAhlberg[[(0 + 1) + n]] =
qAhlberg[[(0 + 1) + n]] vAhlberg[[(0 + 1) + n + 1]] +
uAhlberg[[(0 + 1) + n]];
, {n, nCurve - 1, 0, -1}];
xDD[[(0 + 1) +
nCurve]] = (6 nCurve nCurve (x[[(0 + 1) + 1]] -
2 x[[(0 + 1) + nCurve]] + x[[(0 + 1) + nCurve - 1]]) -
vAhlberg[[(0 + 1) + 1]] - vAhlberg[[(0 + 1) + nCurve - 1]])/(4 +
tAhlberg[[(0 + 1) + 1]] + tAhlberg[[(0 + 1) + nCurve - 1]]);
Do[
xDD[[(0 + 1) + n]] =
tAhlberg[[(0 + 1) + n]] xDD[[(0 + 1) + nCurve]] +
vAhlberg[[(0 + 1) + n]];
, {n, 0, nCurve - 1, 1}];
Do[
AppendTo[
coefs, {x[[(0 + 1) + n]],
nCurve (x[[(0 + 1) + n]] - x[[(0 + 1) + n - 1]]) + (
2 xDD[[(0 + 1) + n]] + xDD[[(0 + 1) + n - 1]])/(6 nCurve),
1/2 xDD[[(0 + 1) + n]],
nCurve (xDD[[(0 + 1) + n]] - xDD[[(0 + 1) + n - 1]])/6}];
, {n, 1, nCurve, 1}];
coefs
]

resampleList[pts_, resampleDist_] :=
Module[{dist, nSeg, n, remain = 0, it = 0},
nSeg = ConstantArray[0, Length[pts]];
dist = ConstantArray[0, Quotient[Fold[Plus, pts], resampleDist]];
For[i = 0, i < Length[pts], i++,
n = Quotient[pts[[i + 1]] + remain, resampleDist];
nSeg[[i + 1]] = n;
For[j = 0, j < n, j++, dist[[it + 1 + j]] = resampleDist - remain];
For[j = 1, j < n, j++,
dist[[it + 1 + j]] = dist[[it + 1 + j - 1]] + resampleDist];
remain = Mod[pts[[i + 1]] + remain, resampleDist];
it += n;
];
{nSeg, dist}
]

calcCurveParam[crv_, range_, segDist_] :=
Module[{a, b, err, numer = 1, denom = 1, totalLength, l},
a = range[[1]];
b = range[[2]];
totalLength =
Total[(cGQ ((b - a)/2 crv) /. \[Xi] -> ((b - a)/2 xGQ + (a + b)/
2)) /. MapThread[{xGQ -> #1, cGQ -> #2} &,
err = segDist;
If[segDist > totalLength, -1,
While[Abs[err] > totalLength 10^-6,
b = numer/2^denom range[[2]];
l = Total[(cGQ ((b - a)/2 crv) /. \[Xi] -> ((b - a)/2 xGQ + (
a + b)/2)) /.
MapThread[{xGQ -> #1, cGQ -> #2} &,
err = segDist - l;
numer = If[err > 0, 2 numer + 1, 2 numer - 1];
denom++;
];
N[b]
]
]

calcCubicSpline[pts_] :=
Module[{splineCoeffX, splineCoeffY, curveCubicSplineEqus},
splineCoeffX = {periodicSplineCoeffCalc[pts[[All, 1]]],
Array[#/Length[pts] &, Length[pts]]};
splineCoeffY = {periodicSplineCoeffCalc[pts[[All, 2]]],
Array[#/Length[pts] &, Length[pts]]};
curveCubicSplineEqus = {Partition[
Array[#/Length[pts] &, Length[pts] + 1, 0], 2, 1],
ComplexExpand[Norm[#]] & /@
Transpose[{(D[
periodicSpline[\[Xi]], {\[Xi], 1}] /. {a0 -> #[[1]][[1]],
a1 -> #[[1]][[2]], a2 -> #[[1]][[3]],
a3 -> #[[1]][[4]], \[Xi]j -> #[[2]]}) & /@
splineCoeffX], (D[
periodicSpline[\[Xi]], {\[Xi], 1}] /. {a0 -> #[[1]][[1]],
a1 -> #[[1]][[2]], a2 -> #[[1]][[3]],
a3 -> #[[1]][[4]], \[Xi]j -> #[[2]]}) & /@
{splineCoeffX, splineCoeffY, curveCubicSplineEqus}
]

resampleSpline[pts_, crv_, resampleDist_] :=
Module[{splineCoeffX, splineCoeffY, curveCubicSplineEqus,
lengthSegGQ, nSeg, dist, it, newPts},
{splineCoeffX, splineCoeffY, curveCubicSplineEqus} = crv;
lengthSegGQ =
Total[(cGQ (Apply[Subtract, Reverse[#1]]/
2 #2) /. \[Xi] -> (Apply[Subtract, Reverse[#1]]/2 xGQ +
Apply[Plus, #1]/2)) /.
MapThread[{xGQ -> #1, cGQ -> #2} &,
Transpose[
GaussianQuadratureWeights[8, -1, 1, 16]]]] &, {Partition[
Array[#/Length[pts] &, Length[pts] + 1, 0], 2, 1],
ComplexExpand[Norm[#]] & /@
Transpose[{(D[
periodicSpline[\[Xi]], {\[Xi], 1}] /. {a0 -> #[[1]][[1]],
a1 -> #[[1]][[2]], a2 -> #[[1]][[3]],
a3 -> #[[1]][[4]], \[Xi]j -> #[[2]]}) & /@
splineCoeffX], (D[
periodicSpline[\[Xi]], {\[Xi], 1}] /. {a0 -> #[[1]][[1]],
a1 -> #[[1]][[2]], a2 -> #[[1]][[3]],
a3 -> #[[1]][[4]], \[Xi]j -> #[[2]]}) & /@
{nSeg, dist} = resampleList[lengthSegGQ, resampleDist];
it = 0;
newPts = ConstantArray[0, Length[dist]];
For[i = 0, i < Length[pts], i++,
n = nSeg[[i + 1]];
For[j = 0, j < n, j++,
newPts[[
it + j +
1]] = {((periodicSpline[\[Xi]] /. {a0 -> #[[1]][[1]],
a1 -> #[[1]][[2]], a2 -> #[[1]][[3]],
a3 -> #[[1]][[4]], \[Xi]j -> #[[2]]}) & /@
i + 1]], ((periodicSpline[\[Xi]] /. {a0 -> #[[1]][[1]],
a1 -> #[[1]][[2]], a2 -> #[[1]][[3]],
a3 -> #[[1]][[4]], \[Xi]j -> #[[2]]}) & /@
MapThread[{#1, #2} &, splineCoeffY])[[i + 1]]} /. \[Xi] ->
calcCurveParam[curveCubicSplineEqus[[2, i + 1]],
curveCubicSplineEqus[[1, i + 1]], dist[[it + j + 1]]]
];
it += n;
];
If[Mod[Total[lengthSegGQ], resampleDist]/resampleDist < 1/3, newPts,
Join[pts[[{1}]], newPts]]
]

curvePts = {{4, 4}, {3, 5}, {2, 4}, {3, 3}};
curveSpline = calcCubicSpline[curvePts];
newPts = resampleSpline[curvePts, curveSpline, 0.12];

ParametricPlot[{periodicSpline[\[Xi]] /. {a0 -> #1[[1]],
a1 -> #1[[2]], a2 -> #1[[3]],
a3 -> #1[[4]], \[Xi]j -> #3[[2]]},
periodicSpline[\[Xi]] /. {a0 -> #2[[1]], a1 -> #2[[2]],
a2 -> #2[[3]],
a3 -> #2[[4]], \[Xi]j -> #3[[2]]}}, {\[Xi], #3[[1]], #3[[2]]},
Axes -> False] &, {curveSpline[[1]][[1]], curveSpline[[2]][[1]],
Partition[Prepend[curveSpline[[1]][[2]], 0], 2, 1]}],
Graphics[{{Red, PointSize[Large],
Point[curvePts]}, {PointSize[Medium], Point[newPts]},
Circle[#, 0.07] & /@ newPts}]}, PlotRange -> All]


There are some great Voronoi iteration / Lloyd algorithm solutions, I wonder if they can be altered to fit for this relaxation problem.

Lloyd's algorithm by J. M. & Stippling Drawing by Silvia Hao

• That is quite a lot of code given that there is actually no question in your post. So what do you expect from us? Commented Feb 19, 2020 at 19:06
• The code is there to demonstrate that I've done some work, not firing questions from the hip - usually this is the first thing that got asked. The question is about how can one relax the points - "look for any points that, if they were spheres with the specified radii, would be overlapping, and attempts to move them to nearby locations that reduce the overlap". Commented Feb 19, 2020 at 19:38
• "The code is there to demonstrate that I've done some work, not firing questions from the hip - usually this is the first thing that got asked. " Right. But we ask users also to be concise and to post only minimal examples. If I understood it correctly, the actual question here could be phrased much shorter like "How to push apart those points in a point cloud that are closer to each other than a given tolerance?" Personally, I usually do not have the time to read through all details and to guess what the actual question is. I have already wasted too a big chunk of my life on that. Commented Feb 19, 2020 at 20:04

Hmm. Nearest would do a good job of finding those point pairs. Then pushing each point with into a direction which proportion to and opposing the sum of those vectors pointing from it to the other neaby points would maybe do want you want.

So let's generate a bunch of points. Then, for each point p, we find each other point q in distance r or less and push p in direction u = p-q with a given stepsize t. Doing all at once looks like this:

pts = RandomPoint[Ball[], 100];
r = 0.2;
t = 0.5;
u = pts - Total[Nearest[pts, pts, {∞, r}], {2}];
newpts = pts + t u;


So in one line it reads as follows:

newpts2 = (1 + t) pts - t Total[Nearest[pts, pts, {∞, r}], {2}];
`