As announced before, here my take on the mean curvature flow for surfaces. The code is rather lengthy and I tried to recycle as much as possible from this post about finding minimal surfaces (solving Plateau's problem). Please find the code at the end of this post.
Background
Mean curvature flow is the $L^2$-gradient flow of the area functional on the space of immersed surfaces. For a time-dependent immersion $f \colon \varSigma \times I \to \mathbb{R}^3$ of a two-dimensional manifold $\varSigma$, the governing partial differential equation is
$$\partial_t f(x,t) = \operatorname{dim}(\varSigma) \, H_f (x,t),$$
where $H_f(x,t)$ is the mean curvature of the surface $f(\varSigma, t)$ at point $f(x,t)$. Note that I understand $H_f$ as a vector-valued function $H_f \colon \varSigma \times I \to \mathbb{R}^3$; it is defined as the trace of the second fundamental form $I\!I_f$ with respect to the Riemannian metric on $\varSigma$ induced by $f$ via pullback of the Euclidean metric along $f$:
$$H_f \colon= \tfrac{1}{\operatorname{dim}(\varSigma)} \operatorname{tr}_f (I\!I_f).$$
The mean curvature can also be written as
$$H_f(x,t) = \tfrac{1}{\operatorname{dim}(\varSigma)} \Delta_{f(\cdot,t)} \,f(x,t),$$
where $\Delta_{f(\cdot,t)}$ denotes the Laplace-Beltrami operator of the surface $f(\varSigma,t)$. This way, the PDE looks a lot like the heat flow PDE
$$\partial_t f - \Delta_{f} \,f = 0,$$
but one has to take into account that $\Delta_{f(\cdot,t)}$ depends on time as well as on $f$, so it is a nonlinear system of PDEs with space- and time-dependent coefficients.
Usually, one considers the mean curvature flow for surfaces without boundary or for Dirichlet boundary conditions. Since we also want to smooth the boundary of surfaces, we apply the curve shortening flow (the 1D-analogon of the mean curvature flow) to the boundary curve $\gamma \colon \partial \varSigma \times I \to \mathbb{R^3}$ and couple these flows in the following way:
$$\begin{aligned}
\partial_t f -\Delta_f \, f &= 0, \quad \text{on $\varSigma \setminus \partial \varSigma$,}\\
\partial_t \gamma - \Delta_\gamma \, \gamma &= 0, \quad \text{on $\partial \varSigma$,}\\
f|_{\partial \varSigma \times I} &= \gamma,
\end{aligned}$$
where $\Delta_\gamma \, \gamma$ equals the curvature vector $\kappa_\gamma$ of $\gamma$.
Like heat flow, mean curvature flow has the strong tendency to remove high freqency oscilliations from the surface while moving the bulk of the surface rather slowly.
That renders the flow rather inefficient for minimizing area.
But here it is an advantange because that's precisely what we need.
Example
n = 100000;
pts = RandomReal[{-1, 1}, {n, 3}];
vals = Dot[Sin[3 pts]^2, ConstantArray[1., 3]] + RandomVariate[NormalDistribution[0, .005], n];
data = Join[pts, Partition[vals, 1], 2];
pl = ListContourPlot3D[data, Contours -> {1.5},
PerformanceGoal -> "Quality",
Mesh -> None, ContourStyle -> Directive[EdgeForm[Thin]],
MaxPlotPoints -> 50
];
R = RepairMesh[DiscretizeGraphics[pl],
{"TinyComponents", "TinyFaces", "IsolatedVertices", "SingularVertices", "DanglingEdges", "TJunctionEdges"},
PerformanceGoal -> "Quality",
MeshCellStyle -> {{2, All} -> Directive[Darker@Orange, Specularity[White, 30]]}
]
Let's apply 5 steps of mean curvature flow with stepzise 0.00125 and theta-value 0.8:
S = MeanCurvatureFlow[R, 5, 0.00125, 0.8]
Here a direct comparison:
Show[R, S]
Usage Notes
Finding good step sizes is usually quite a mess. The integrators for the PDE require something like stepsize ~ minimal triangle diameter of the current mesh. As a rule of thumb, one should determine the stepsize as a multiple of
ρ = Min[PropertyValue[{R, 1}, MeshCellMeasure]];
If the Min is too small, Mean might also do.
Moreover, mean curvature flow is known to develop singularities within finite time. Remember: Mean curvature flow is the $L^2$-gradient flow of area. That means that a closed, connected surface will inevitably shrink to a point. With the boundary components following a curve shortening flow, they also try to collapse to points. So the interior of the face and its boundary components struggle both for minimality which leads to some intricate interplay for large time horizons. Moreover, bottleneck regions have the tendency to collapse to lines (with a faster rate than the overall collapse to a point) and this is what happens with the ears of the Stanford bunny (thanks to chris for pointing me to this):
R = ExampleData[{"Geometry3D", "StanfordBunny"}, "MeshRegion"];
ρ = Min[PropertyValue[{R, 1}, MeshCellMeasure]];
NestList[GraphDiffusionFlow[#, 1, ρ, 0.8] &, R, 4]
This is a well-known (and feared) issue in geometry processing. Somewhat more desired behavior can be obtained by shrinking the time horizon by a factor of 100:
NestList[MeanCurvatureFlow[#, 1, ρ/100, 0.8] &, R, 5]
Moreover, replacing the Laplace-Betrami operator by the graph Laplacian of the underlying edge graph of the mesh leads to a flow with seemingly better longtime behavior. This is also called Laplacian smoothing. It basically equivalent to successivly averaging vertex positions with the positions of the direct neighbor vertices (with a special treatment of boundary vertices). This is very similar to kglr's method, yet, the averaging stencil is chosen by connectivity and not by distance.
NestList[GraphDiffusionFlow[#, 25, 0.125, 0.8] &, R, 4]
Code Dump
This is the code for assempling mass matrices and discrete Laplace-Beltrami operators for the surface and its boundary curves.
Block[{xx, x, PP, P, UU, U, VV, V, f, Df, u, Du, v, Dv, g, integrand, quadraturepoints, quadratureweights},
xx = Table[Compile`GetElement[x, i], {i, 1, 1}];
PP = Table[Compile`GetElement[P, i, j], {i, 1, 2}, {j, 1, 3}];
UU = Table[Compile`GetElement[U, i], {i, 1, 2}];
VV = Table[Compile`GetElement[V, i], {i, 1, 2}];
(*local affine parameterization of the curve with respect to the unit interval*)
f = x \[Function] PP[[1]] + x[[1]] (PP[[2]] - PP[[1]]);
Df = x \[Function] Evaluate[D[f[xx], {xx}]];
(*the Riemannian pullback metric with respect to f*)
g = x \[Function] Evaluate[Df[xx]\[Transpose].Df[xx]];
(*two affine functions u and v and their derivatives*)
u = x \[Function] UU[[1]] + x[[1]] (UU[[2]] - UU[[1]]);
Du = x \[Function] Evaluate[D[u[xx], {xx}]];
v = x \[Function] VV[[1]] + x[[1]] (VV[[2]] - VV[[1]]);
Dv = x \[Function] Evaluate[D[v[xx], {xx}]];
integrand = x \[Function] Evaluate[D[D[v[xx] u[xx] Sqrt[Abs[Det[g[xx]]]], {UU}, {VV}]]];
(*since the integrand is quadratic over each edge,we use a two-
point Gauss quadrature rule (for the standard triangle)*)
{quadraturepoints, quadratureweights} = Most[NIntegrate`GaussRuleData[2, $MachinePrecision]];
quadraturepoints = Partition[quadraturepoints, 1];
getCurveMass =
With[{code = N[quadratureweights.Map[integrand, quadraturepoints]]},
Compile[{{P, _Real, 2}}, code, CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"]];
integrand = x \[Function] Evaluate[D[D[Dv[xx].Inverse[g[xx]].Du[xx] Sqrt[Abs[Det[g[xx]]]], {UU}, {VV}]]];
(*since the integrand is constant over each edge,we use a one-
point Gauss quadrature rule (for the standard triangle)*)
quadraturepoints = {{1/2}};
quadratureweights = {1};
getCurveLaplaceBeltrami =
With[{code = Together@N[quadratureweights.Map[integrand, quadraturepoints]]},
Compile[{{P, _Real, 2}}, code, CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"
]
]
];
getCurveLaplacianCombinatorics =
Quiet[Module[{ff},
With[{code = Flatten[Table[Table[{ff[[i]], ff[[j]]}, {i, 1, 2}], {j, 1, 2}], 1]},
Compile[{{ff, _Integer, 1}}, code,
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"]]]];
CurveLaplaceBeltrami[pts_, flist_, pat_] :=
With[{
spopt = SystemOptions["SparseArrayOptions"],
vals = Flatten[getCurveLaplaceBeltrami[Partition[pts[[flist]], 2]]]
},
Internal`WithLocalSettings[
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
SparseArray[Rule[pat, vals], {Length[pts], Length[pts]}, 0.],
SetSystemOptions[spopt]]];
CurveMassMatrix[pts_, flist_, pat_] :=
With[{
spopt = SystemOptions["SparseArrayOptions"],
vals = Flatten[getCurveMass[Partition[pts[[flist]], 2]]]
},
Internal`WithLocalSettings[
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
SparseArray[Rule[pat, vals], {Length[pts], Length[pts]}, 0.],
SetSystemOptions[spopt]]];
Block[{xx, x, PP, P, UU, U, VV, V, f, Df, u, Du, v, Dv, g, integranf, integrand, quadraturepoints, quadratureweights},
xx = Table[Compile`GetElement[x, i], {i, 1, 2}];
PP = Table[Compile`GetElement[P, i, j], {i, 1, 3}, {j, 1, 3}];
UU = Table[Compile`GetElement[U, i], {i, 1, 3}];
VV = Table[Compile`GetElement[V, i], {i, 1, 3}];
(*local affine parameterization of the surface with respect to the \
"standard triangle"*)
f = x \[Function] PP[[1]] + x[[1]] (PP[[2]] - PP[[1]]) + x[[2]] (PP[[3]] - PP[[1]]);
Df = x \[Function] Evaluate[D[f[xx], {xx}]];
(*the Riemannian pullback metric with respect to f*)
g = x \[Function] Evaluate[Df[xx]\[Transpose].Df[xx]];
(*two affine functions u and v and their derivatives*)
u = x \[Function] UU[[1]] + x[[1]] (UU[[2]] - UU[[1]]) + x[[2]] (UU[[3]] - UU[[1]]);
Du = x \[Function] Evaluate[D[u[xx], {xx}]];
v = x \[Function] VV[[1]] + x[[1]] (VV[[2]] - VV[[1]]) + x[[2]] (VV[[3]] - VV[[1]]);
Dv = x \[Function] Evaluate[D[v[xx], {xx}]];
integrand = x \[Function] Evaluate[D[D[v[xx] u[xx] Sqrt[Abs[Det[g[xx]]]], {UU}, {VV}]]];
(*since the integrand is quadratic over each triangle,
we use a three-point Gauss quadrature rule (for the standard triangle)*)
quadraturepoints = {{0, 1/2}, {1/2, 0}, {1/2, 1/2}};
quadratureweights = {1/6, 1/6, 1/6};
getSurfaceMass =
With[{code = N[quadratureweights.Map[integrand, quadraturepoints]]},
Compile[{{P, _Real, 2}}, code, CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"]];
integrand = x \[Function] Evaluate[D[D[Dv[xx].Inverse[g[xx]].Du[xx] Sqrt[Abs[Det[g[xx]]]], {UU}, {VV}]]];
(*since the integrand is constant over each triangle,we use a one-
point Gauss quadrature rule (for the standard triangle)*)
quadraturepoints = {{1/3, 1/3}};
quadratureweights = {1/2};
getSurfaceLaplaceBeltrami =
With[{code = N[quadratureweights.Map[integrand, quadraturepoints]]},
Compile[{{P, _Real, 2}}, code, CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"]]];
getSurfaceLaplacianCombinatorics =
Quiet[Module[{ff},
With[{code = Flatten[Table[Table[{ff[[i]], ff[[j]]}, {i, 1, 3}], {j, 1, 3}], 1]},
Compile[{{ff, _Integer, 1}}, code, CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"]]]];
SurfaceLaplaceBeltrami[pts_, flist_, pat_] :=
With[{
spopt = SystemOptions["SparseArrayOptions"],
vals = Flatten[getSurfaceLaplaceBeltrami[Partition[pts[[flist]], 3]]]
},
Internal`WithLocalSettings[
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
SparseArray[Rule[pat, vals], {Length[pts], Length[pts]}, 0.],
SetSystemOptions[spopt]]];
SurfaceMassMatrix[pts_, flist_, pat_] :=
With[{spopt = SystemOptions["SparseArrayOptions"], vals = Flatten[getSurfaceMass[Partition[pts[[flist]], 3]]]},
Internal`WithLocalSettings[
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
SparseArray[Rule[pat, vals], {Length[pts], Length[pts]}, 0.], SetSystemOptions[spopt]]];
And this is the actual code for the mean curvature flow. This implements a semi-implicit $\theta$-method for integrating the flow; θ = 0.5 resemples the Crank-Nicolson scheme while θ = 1. has an implicit-Euler flavor. Note however that the integration method is not fully implicit. On the one hand, θ = 1. need not be stable (it throws usually a lot of numerical errors). On the other hand, values of θ too close to 0.5 will lead to spikes oscillating in time (a notorious behavior of the Crank-Nicolson scheme for not-so-smooth data). A good trade-off can be obtained with values of θ between 0.6 and 0.8
MeanCurvatureFlow::infy =
"Division by zero detected in computation of `1`. Flow is getting singular. Aborting the flow in step `2`.";
MeanCurvatureFlow[R_MeshRegion, steps_, stepsize_, θ_] :=
Module[{bedges, belist, faces, flist, pts, bpat, bplist, pat, a, m, aplus, aminus, τ},
τ = stepsize;
bedges = MeshCells[R, 1, "Multicells" -> True][[1, 1,
Random`Private`PositionsOf[Length /@ R["ConnectivityMatrix"[1, 2]]["AdjacencyLists"], 1]]];
belist = Flatten[bedges];
faces = MeshCells[R, 2, "Multicells" -> True][[1, 1]];
flist = Flatten[faces];
pts = MeshCoordinates[R];
bpat = If[Length[bedges] > 0, Flatten[getCurveLaplacianCombinatorics[bedges], 1], {}];
bplist = Sort[DeleteDuplicates[belist]];
pat = Flatten[getSurfaceLaplacianCombinatorics[faces], 1];
Do[
Check[
a = SurfaceLaplaceBeltrami[pts, flist, pat],
Message[MeanCurvatureFlow::infy, SurfaceLaplaceBeltrami, i];
Break[],
Power::infy
];
Check[
m = SurfaceMassMatrix[pts, flist, pat],
Message[MeanCurvatureFlow::infy, SurfaceMassMatrix, i];
Break[],
Power::infy
];
If[Length[bpat] > 0,
Check[
a[[bplist]] = CurveLaplaceBeltrami[pts, belist, bpat][[bplist]],
Message[MeanCurvatureFlow::infy, CurveLaplaceBeltrami, i];
Break[],
Power::infy
];
Check[
m[[bplist]] = CurveMassMatrix[pts, belist, bpat][[bplist]],
Message[MeanCurvatureFlow::infy, CurveMassMatrix, i];
Break[],
Power::infy
];
];
aplus = m + (θ τ) a;
aminus = m + ((1. - θ) τ) a;
pts = LinearSolve[aplus, aminus.pts];
,
{i, 1, steps}];
MeshRegion[pts, Polygon[faces]]
]
Addendum: Laplacian Smoothing
Using the graph Laplacian of the triangle mesh leads to an algorithm with similar smoothing behavoir which is also 1.) faster (since we have to factorize only one matrix), 2.) easier to implement, and 3.) probably more robust:
GraphDiffusionFlow[R_MeshRegion, steps_, stepsize_, θ_] :=
Module[{n, belist, pts, bplist, a, m, aplus, aminus, τ, edges, bedges, solve},
τ = stepsize;
n = MeshCellCount[R, 0];
edges = MeshCells[R, 1, "Multicells" -> True][[1, 1]];
a = GraphLaplacian[n, edges];
m = IdentityMatrix[Length[a], SparseArray];
belist = Random`Private`PositionsOf[Length /@ R["ConnectivityMatrix"[1, 2]]["AdjacencyLists"], 1];
If[Length[belist] > 0,
bedges = edges[[belist]];
bplist = Sort[DeleteDuplicates[Join @@ bedges]];
a[[bplist]] = GraphLaplacian[n, bedges][[bplist]];
bedges =.;
m[[bplist]] = IdentityMatrix[n, SparseArray][[bplist]];
bplist =.;
];
aplus = m + (τ θ) a;
aminus = m - (τ (1 - θ)) a;
pts = MeshCoordinates[R];
solve = LinearSolve[aplus];
Do[pts = solve[aminus.pts];, {i, 1, steps}];
MeshRegion[pts, MeshCells[R, 2, "Multicells" -> True]]]
GraphLaplacian[n_Integer,
edges_: List[List[i_Integer, j_Integer] ..]] := With[{
A = SparseArray[
Rule[
Join[edges, Transpose[Transpose[edges][[{2, 1}]]]],
ConstantArray[1, 2 Length[edges]]
],
{n, n}
]},
SparseArray[DiagonalMatrix[SparseArray[Total[A]]] - A]
]
Usage example:
T = GraphDiffusionFlow[R, 20, 0.25, 0.8];
Show[R, T]
