20
$\begingroup$

Context

Following this question (and great answer!), It would be nice to have a function which also smooths 3D contours plots once they have been done. There are various solutions which involve smoothing the data before making the contours, but here I am after a solution which operates on the 3D graphics itself.

Example

Let s say I produced the following plot

n = 10000; pts = RandomReal[{-1, 1}, {n, 3}]; 
vals = Dot[pts^2, ConstantArray[1., 3]] + RandomVariate[NormalDistribution[0, .15], n]; 
data = Join[pts, Partition[vals, 1], 2]; 
pl = ListContourPlot3D[data, Contours -> {0.5}, PerformanceGoal -> "Quality"]

Mathematica graphics

and I only have the plot at this stage (not the data).

Question

I would like to smooth this 3D contour as a direct post processing of the graphics3D file.

Ideally using a function like

  smoothCP3D[pl,GaussianFilter[#, 5] &]]

very much like the one presented here.

This problem is most likely related to mesh simplification which is a big field in itself.

$\endgroup$
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – Kuba May 4 '18 at 11:26
26
$\begingroup$

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

enter image description here

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]

enter image description here

Here a direct comparison:

Show[R, S]

enter image description here

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]

enter image description here

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]

enter image description here

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]

enter image description here

Code Dump

This is the code for assempling mass matrices and discrete Laplace-Beltrami operators for the surface and its boundary curves.

Quiet[Block[{xx, x, PP, P, UU, U, VV, V, f, Df, u, Du, v, Dv, g, integrant, quadraturepoints, quadratureweights}, 
   xx = Table[x[[i]], {i, 1, 1}];
   PP = Table[P[[i, j]], {i, 1, 2}, {j, 1, 3}];
   UU = Table[U[[i]], {i, 1, 2}];
   VV = Table[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}]];

   integrant = x \[Function] Evaluate[D[D[v[xx] u[xx] Sqrt[Abs[Det[g[xx]]]], {UU}, {VV}]]];
   (*since the integrant 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[integrant, quadraturepoints]] /. Part -> Compile`GetElement}, Compile[{{P, _Real, 2}},
      code,
      CompilationTarget -> "C",
      RuntimeAttributes -> {Listable},
      Parallelization -> True,
      RuntimeOptions -> "Speed"
      ]
     ];
   integrant = 
    x \[Function] Evaluate[D[D[Dv[xx].Inverse[g[xx]].Du[xx] Sqrt[ Abs[Det[g[xx]]]], {UU}, {VV}]]];
   (*since the integrant 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[integrant, quadraturepoints]] /. Part -> Compile`GetElement}, 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] /. Part -> Compile`GetElement}, 
     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]]];



Quiet[Block[{xx, x, PP, P, UU, U, VV, V, f, Df, u, Du, v, Dv, g, 
    integrant, quadraturepoints, quadratureweights}, 
   xx = Table[x[[i]], {i, 1, 2}];
   PP = Table[P[[i, j]], {i, 1, 3}, {j, 1, 3}];
   UU = Table[U[[i]], {i, 1, 3}];
   VV = Table[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}]];

   integrant = x \[Function] Evaluate[D[D[v[xx] u[xx] Sqrt[Abs[Det[g[xx]]]], {UU}, {VV}]]];
   (*since the integrant 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[integrant, quadraturepoints]] /. Part -> Compile`GetElement}, Compile[{{P, _Real, 2}},
      code,
      CompilationTarget -> "C",
      RuntimeAttributes -> {Listable},
      Parallelization -> True,
      RuntimeOptions -> "Speed"
      ]
     ];
   integrant = 
    x \[Function] 
     Evaluate[
      D[D[Dv[xx].Inverse[g[xx]].Du[xx] Sqrt[
          Abs[Det[g[xx]]]], {UU}, {VV}]]];
   (*since the integrant is constant over each trianle,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[integrant, quadraturepoints]] /. Part -> Compile`GetElement}, 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] /. Part -> Compile`GetElement}, 
     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]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you, chris! $\endgroup$ – Henrik Schumacher May 5 '18 at 9:51
  • $\begingroup$ Out of curiosity why is the mean curvature called a vector. It seems to be a scalar? $\endgroup$ – chris May 5 '18 at 9:54
  • 1
    $\begingroup$ Well, that's a question of nomenclatur. Scalar mean curvature only exists for oriented surfaces of codimension 1 (a $n$-dimensional surface in $m=n+1$-dimensional Euclidean space - or any other oriented Riemannian manifold of dimension $n+1$; that means, only when the surface normal is well-defined), while the mean curvature vector is defined for arbitrary codimension and orientability as the trace of the second fundamental form $I\!I$ (usually divided by $n$). $\endgroup$ – Henrik Schumacher May 5 '18 at 10:01
  • $\begingroup$ The second fundamental form at point $x$ is in general a bilinear map $I\!I_x \colon T_x\varSigma \times T_x\varSigma \to \mathbb{R}^m$. Usually, it gets contracted with the surface normal $\nu(x)$ in order to get a biliear form $ \colon T_x\varSigma \times T_x\varSigma \to \mathbb{R}$. Note that $\Delta_f f$ is also vector valued for $f$ is vector valued. So the PDE governing mean curvature flow is actually a nonlinear system of PDEs. $\endgroup$ – Henrik Schumacher May 5 '18 at 10:02
  • $\begingroup$ I am not exactly sure what you mean? Mean curvature flow can contract triangles to points or lines, so there might always occur numerical problems. I tried to make the code more robust by introducing several Checks. I would expect that the flow smoothes the boundaries a bit better than the generic averaging procedure. $\endgroup$ – Henrik Schumacher May 5 '18 at 11:33
8
$\begingroup$
SeedRandom[7]
n = 10000; pts = RandomReal[{-1, 1}, {n, 3}];
vals = Dot[pts^2, ConstantArray[1., 3]] + 
   RandomVariate[NormalDistribution[0, .15], n];
data = Join[pts, Partition[vals, 1], 2];
pl = ListContourPlot3D[data, Contours -> {0.5}, 
   ContourStyle -> Directive[Orange, Opacity[0.5], Specularity[White, 30]], 
   PerformanceGoal -> "Quality", ImageSize -> 300];

Playing with parameters m and k gives something not too far off:

m = 200; k = 10;
pl2 = pl /.  GraphicsComplex[a_, b_, VertexNormals -> vn_, c___] :> 
    Module[{nf = Nearest[a -> Automatic]}, 
       GraphicsComplex[Mean[a[[nf[#, m]]] + vn[[nf[#, m]]]/k] & /@ a, b, 
         VertexNormals -> vn, c]] /. Orange -> Green /. Opacity[.5] -> Opacity[.8];
Row[{pl, pl2, Show[pl, pl2]}]

enter image description here

With m = 20; k = 100; we get

enter image description here

$\endgroup$
  • $\begingroup$ I wish I could have given two votes. Thanks again! $\endgroup$ – chris May 5 '18 at 13:51
  • $\begingroup$ @chris, thank you indeed for a great question (without which we would have never known about Henrik's superb answer). $\endgroup$ – kglr May 5 '18 at 14:29

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.