Curve shortening flow

Question

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:

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

Answer 1

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