# Curve shortening flow

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;
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*)
Sqrt[1 + (pa[[2]] - pb[[2]])^2/(pa[[1]] - pb[[1]])^2] +
pa[[1]]];
ynew =
If[pab[[1]] == 0,
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:

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

• 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? – Jan 29 at 1:34

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]]*)
]


• There is a “Preserve Area” button on the right side of the website. Can this be implemented?
– yode
Jan 29 at 4:42
• @Goofy Very interesting answer! Could you please describe the estimation of the curvature in more detail? Thanks! Jan 29 at 8:27
• @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. Jan 30 at 21:49
• @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. Jan 30 at 21:53
• @UlrichNeumann One way to derive a formula is to interpolate and differentiate. For $2n+1$ points, you can get it with With[{n = 2}, NDSolveFiniteDifferenceDerivative[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$!) Feb 1 at 11:31