unbounded Voronoi iteration with finite and non-equal point radii

Question

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} &, 
      Transpose[Chop[GaussianQuadratureWeights[4, -1, 1, 16]]]]];
  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} &, 
        Transpose[Chop[GaussianQuadratureWeights[4, -1, 1, 16]]]]];
    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]]}) & /@ 
        MapThread[{#1, #2} &, 
         splineCoeffX], (D[
            periodicSpline[\[Xi]], {\[Xi], 1}] /. {a0 -> #[[1]][[1]], 
            a1 -> #[[1]][[2]], a2 -> #[[1]][[3]], 
            a3 -> #[[1]][[4]], \[Xi]j -> #[[2]]}) & /@ 
        MapThread[{#1, #2} &, splineCoeffY]}]};
  {splineCoeffX, splineCoeffY, curveCubicSplineEqus}
  ]

resampleSpline[pts_, crv_, resampleDist_] := 
 Module[{splineCoeffX, splineCoeffY, curveCubicSplineEqus, 
   lengthSegGQ, nSeg, dist, it, newPts},
  {splineCoeffX, splineCoeffY, curveCubicSplineEqus} = crv;
  lengthSegGQ = 
   MapThread[
    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]]}) & /@ 
         MapThread[{#1, #2} &, 
          splineCoeffX], (D[
             periodicSpline[\[Xi]], {\[Xi], 1}] /. {a0 -> #[[1]][[1]],
              a1 -> #[[1]][[2]], a2 -> #[[1]][[3]], 
             a3 -> #[[1]][[4]], \[Xi]j -> #[[2]]}) & /@ 
         MapThread[{#1, #2} &, splineCoeffY]}]}];
  {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]]}) & /@ 
          MapThread[{#1, #2} &, splineCoeffX])[[
        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];

Show[{MapThread[
   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

Answer 1

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}];