13
$\begingroup$

I'd like to use mma to recreate this curve shortening flow effect. I have something that works for simple shapes (LHS), but not for more complex curves, which it causes to self-intersect (RHS gif):

Clear[f, np, p, t];
p1 = RandomReal[{-1, 1}, {7, 2}];
f1 = BSplineFunction[p1, SplineDegree -> 5, SplineClosed -> True];
np[u_, dt_] := u + dt/Norm[D[f1[t], t]] /. t -> u;
newpts = 
  Table[f1[t], {t, NestWhileList[np[#, 1/25.] &, 0, # < 1 &]}];
newPursuitList[1, 250, newpts]

What is the correct approach for this?

Code for newPursuitList function used above:

newPursuitList[k_, n_, list_] := 
  Module[{}, 
   folf[l_List, i_] := 
    newptRadiusMid[##, i] & @@@ Partition[l, 3, 1, 2];
   midpt[pc_, pd_] := (Tr /@ Thread@{pc, pd})/2;
   newptRadiusMid[pa_, pc_, pd_, time_] := 
    Module[{pnew, xnew, ynew, pab, radius, pb}, pb = midpt[pc, pd];
     incr = 1; dist = EuclideanDistance[pa, pb];
     radius = If[dist > incr, incr, dist];
     pab = {pb[[1]] - pa[[1]], pb[[2]] - pa[[2]]};
     xnew = 
      If[pab[[1]] == 0, 
       pa[[1]],(*rate*)(1/2)(*rate*)
         Sign[pab[[1]]] radius/
          Sqrt[1 + (pa[[2]] - pb[[2]])^2/(pa[[1]] - pb[[1]])^2] + 
        pa[[1]]];
     ynew = 
      If[pab[[1]] == 0, 
       pa[[2]] +(*rate*)(1/2)(*rate*)Sign[pab[[2]]] radius, 
       pa[[2]] + (xnew - pa[[1]]) #[[2]]/#[[1]] &@pab];
     pnew = {xnew, ynew}];
   abs[li_, p_] := 
    Module[{a, b, c, d, e, f, g}, 
     a = {Min@#, Max@#} &@li[[p]][[All, 1]];
     b = {Min@#, Max@#} &@li[[p]][[All, 2]];
     c = Abs[Differences@a];
     d = Abs[Differences@b];
     e = Max@{c, d}/2; f = 1.1; g = f*e;
     {{# - g, # + g} &@Mean@a, {# - g, # + g} &@Mean@b}];
   newlist = Chop /@ FoldList[folf, list, Range@n];
   plots = 
    Table[Graphics[{Line[Join[{Last@#}, #] &@newlist[[m]]]}, 
      PlotRange -> {{# - .2 Abs@# &@Min@#[[All, 1]], # + .2 Abs@# &@
            Max@#[[All, 1]]}, {# - .2 Abs@# &@
            Min@#[[All, 2]], # + .2 Abs@# &@Max@#[[All, 2]]}} &@list, 
      Frame -> False], {m, 1, n, 1}];
   ListAnimate[plots]];

Update

@Goofy's code works as desired with swirl:

enter image description here

pts1 = Catenate[
   Cases[PolarPlot[Sqrt[t], {t, 0, 5 Pi}], Line[data_] :> data, 
    Infinity]];
pts2 = Catenate[
   Cases[PolarPlot[-Sqrt[t], {t, 0, 4 Pi}], Line[data_] :> data, 
    Infinity]];
pts = Join[pts1, Reverse@pts2]/(2 Pi);
Graphics[{Polygon@pts}];

DynamicModule[{pts = remesh@Most@pts, rate = 0.0005}, 
 Graphics[{Dynamic@
    Polygon[{pts = 
       remesh[pts + 
         rate*RotateLeft@
           Divide[Transpose@{ListConvolve[{1, -2, 1}, pts[[All, 1]], 
               1], ListConvolve[{1, -2, 1}, pts[[All, 2]], 1]}, 
            ListConvolve[{1, 0, -1}, pts[[All, 1]], 1]^2 + 
             ListConvolve[{1, 0, -1}, pts[[All, 2]], 1]^2]]}]}, 
  Frame -> True, PlotRange -> 1]]
$\endgroup$
1
  • $\begingroup$ By the Gage–Hamilton–Grayson theorem,the Curve Shortening Flow should at first make a simple closed curve converge to a convex closed curve and then to a circle. In the above code, if there an example convergence to a circle? – $\endgroup$
    – cvgmt
    Jan 29 at 1:34

1 Answer 1

11
$\begingroup$

This is roughly the algorithm in the link. Like it says there, it's potentially numerically unstable. The funny curlicue example in the OP was not provided, so I can't check that. In fact, most of the time, the provided random example gives an initial curve that is self-intersecting, which I thought was to be avoided. Anyway, here's the idea in the link:

remesh[{}, ___] := {};
remesh[pp_, threshold_ : 0.0008] := Replace[
   Reap[
     Module[{oldp = Sow@First@pp},
      Do[If[SquaredEuclideanDistance[oldp, p] > threshold,
        Sow[oldp = p]],
       {p, Rest@pp}]
      ]
     ][[2, 1]],
   l_List /; Length[l] < 5 :> {}];
DynamicModule[{
  pts = remesh@Most@Table[f1[t], {t, 0., 1., 1/128}],
  rate = 0.001},
 (*Table[*)  (* use Table instead of Dynamic for GIF *)
 Graphics[{
   Dynamic@Line[{pts = remesh[
        pts + rate*RotateLeft@
           Divide[ (* compute 4*laplacian *)
            Transpose@{
              ListConvolve[{1, -2, 1}, pts[[All, 1]], 1],
              ListConvolve[{1, -2, 1}, pts[[All, 2]], 1]},
            ListConvolve[{1, 0, -1}, pts[[All, 1]], 1]^2 +
             ListConvolve[{1, 0, -1}, pts[[All, 2]], 1]^2
            ]
        ]}]
   },
  Frame -> True, PlotRange -> 1](*,
 500][[;;;;10]]*)
 ]

enter image description here

$\endgroup$
7
  • 1
    $\begingroup$ There is a “Preserve Area” button on the right side of the website. Can this be implemented? $\endgroup$
    – yode
    Jan 29 at 4:42
  • $\begingroup$ @Goofy Very interesting answer! Could you please describe the estimation of the curvature in more detail? Thanks! $\endgroup$ Jan 29 at 8:27
  • 1
    $\begingroup$ @UlrichNeumann $\kappa = d^2{\bf r}/ds^2 \approx [{\bf r}(t+h)-2{\bf r}(t)+{\bf r}(t-h)]/\Delta s^2$, $\Delta s \approx \|\Delta {\bf r}\|/2=\|{\bf r}(t+h)-{\bf r}(t-h)\|/2$. The missing factor of $0.25$ is included in the rate. $\endgroup$
    – Goofy
    Jan 30 at 21:49
  • $\begingroup$ @yode I don't see why not. It's probably why remesh() in the link adds points when a gap grows between adjacent points of the curve. I couldn't see a reason for it in this use case and stuck with a simpler answer. $\endgroup$
    – Goofy
    Jan 30 at 21:53
  • 1
    $\begingroup$ @UlrichNeumann One way to derive a formula is to interpolate and differentiate. For $2n+1$ points, you can get it with With[{n = 2}, NDSolve`FiniteDifferenceDerivative[2, Array[x, 2 n + 1, 0], Array[r, 2 n + 1, 0], "DifferenceOrder" -> 2 n - 1][[n + 1]] ] /. dx : HoldPattern[-x[j_] + x[i_]] :> (j - i) ds // Simplify, where $ds$ is the average step size. If the step sizes vary greatly from step to step, one could use dx : HoldPattern[-_x + _x] :> Norm[dx] instead of the previous rule. (Luckily, FiniteDifferenceDerivative gives a formula in terms of $\Delta x$!) $\endgroup$
    – Goofy
    Feb 1 at 11:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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