2
$\begingroup$

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]

relax

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

$\endgroup$
  • $\begingroup$ That is quite a lot of code given that there is actually no question in your post. So what do you expect from us? $\endgroup$ – Henrik Schumacher Feb 19 at 19:06
  • $\begingroup$ 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". $\endgroup$ – Toorop Feb 19 at 19:38
  • $\begingroup$ "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. $\endgroup$ – Henrik Schumacher Feb 19 at 20:04
1
$\begingroup$

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}];
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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