3
$\begingroup$

I'm trying to calculate the mean curvature over the surface of a discrete mesh that I generated from the computational software SE (Surface Evolver) and I want to make sure I can calculate the same curvature energy I got from SE.

In SE, I used the built-in function star_perp_sq_mean_curvature to calculate the curvature energy. There is a software manual https://kenbrakke.com/evolver/html/quants.htm#star_perp_sq_mean_curvature however, I think a lot of details are under the hood.

What I want to understand is,

  1. this 30 yo software calculates the curvature energy in a flash second while my code takes mins
  2. when the surface of the mesh is pretty smooth, the value I calculate is not much different from that of SE but once the mesh has some locally sharp patches, the calculation becomes quite off from SE

I think there are a lot of ways to approximate the mean curvature, for instance, the method here How to speed up estimation of Mean and Gaussian curvatures on triangular meshes? works pretty well but not as fast as SE does and the numbers are not the same.

Maybe I will have to open up the hood to completely mimic the calculation, but I wonder in case, if anyone knows how to code the mean curvature calculation that gives the same number with SE, or if there is a way to calculate the mean curvature of a mesh with ~10^4 vertices in a second with a good approximation level.

$\endgroup$
2
  • $\begingroup$ Is this a question about Wolfram Mathematica? Because there is no code or anything similar ... You have to at least provide a small example of a mesh and the expected result. $\endgroup$
    – Domen
    Feb 12 at 18:36
  • 3
    $\begingroup$ The 30 year old software was coded in portable C and was designed by people who knew what they were doing. $\endgroup$ Feb 12 at 18:52

1 Answer 1

7
$\begingroup$

Not sure what the Surface Evolver does, but the following should be a quite quick an appropriate approximation of the squared mean curvature integral. Note that I am using code from this post: https://mathematica.stackexchange.com/a/158356/38178.

(*Assuming that S is a MeshRegion that represents a _closed_ triangle surface with consistent orientation of triangles.*)
pts = MeshCoordinates[S];
flist = MeshCells[S, 2, "Multicells" -> True][[1, 1]];
pat = Flatten[getLaplacianCombinatorics[flist], 1];
A = LaplaceBeltrami[pts, Flatten[flist], pat];

Block[{i, j, k, U1, U2, U3, V1, V2, V3},
  {i, j, k} = Transpose[flist];
  {U1, U2, U3} = Transpose[pts[[j]] - pts[[i]]];
  {V1, V2, V3} = Transpose[pts[[k]] - pts[[i]]];
  
  triangleareanormals = 
   0.5 Transpose[{U2 V3 - U3 V2, U3 V1 - U1 V3, U1 V2 - U2 V1}];
  
  triangleareas = Sqrt[Total[triangleareanormals^2, {2}]];
];


Av = With[{spopt = SystemOptions["SparseArrayOptions"],
    vals = ConstantArray[(1/3.), 3 Length[flist]],
    pat = Transpose[{
       Flatten[flist],
       Flatten[Transpose[ConstantArray[Range[Length[flist]], 3]]]
       }]
    }, 
    Internal`WithLocalSettings[
      SetSystemOptions["SparseArrayOptions" ->{"TreatRepeatedEntries" -> Total}], 
      SparseArray[Rule[pat, vals], {Length[pts], Length[flist]}, 0.],
      SetSystemOptions[spopt]
     ]
   ];



vertexareas = Av . triangleareas;

vertexinvareas = 1./vertexareas;

vertexnormals = (Av.triangleareanormals) vertexinvareas;

A = LaplaceBeltrami[pts, Flatten[flist], pat];

H = 0.5 ((A.pts) vertexinvareas);

Hproj = vertexnormals Total[vertexnormals H, {2}] / Total[vertexnormals^2, {2}];

HSquareIntegral = Total[H^2, {2}] . vertexareas;
HprojSquareIntegral = Total[Hproj^2, {2}] . vertexareas;

Here A is the weak formulation of the Laplace-Betrami operator, pts are the vertex coordinates and A.pts is exactly the derivative of the the surface area with respect to all the vertex coordinates. H is therefore an approximation of the $L^2$-gradient of surface area; hence one of many possible approximations of the mean curvature vector field.

Surface Evolver seems to project this onto the $L^2$-gradient of the volume functional which is typically a good approximation for vertex normals. I think this just the area-weigted normals, but I could be wrong.

This runs on my machine for a mesh with 1.5 million triangles in about 0.85 seconds. It does some superfluous work though (one does not really have to compute the Laplace-Beltrami operator), and coding this in a low-level language like C++ and using parallelization, one can get this below a quarter of a second or less.

Note that no discretization is perfect. (Some are way worse than others, though.) Almost all have defects when applied to unsmooth meshes. So keeping your meshes in a good state (for example, by remeshing) is pivotal.

$\endgroup$
3
  • $\begingroup$ wow I got really impressed about the speed! Thanks a lot! However, I am having some trouble connecting this to the continuum curvature energy Could you help me a bit with this? I thought HprojSquareIntegral in this code is supposed to approximate the area integral of mean curvature square So it should output 8Pi for a sphere but it seems like to be about 1/5 of 8Pi? Am I missing something or totally misunderstanding? Thanks again $\endgroup$
    – physgj
    Feb 14 at 0:36
  • $\begingroup$ Thank you for pointing this out! The sphere has area $4 \pi$ and mean curvature $1$. So the energy should be $4 \pi \approx 12.5664$. My code deviates from this for two reasons: 1.) I forgot to divide my H by 0.5. 2.) In the definition of Av I divided by .3 instead of 3. ^^ Now for a nice sphere mesh with 82000 triangles, I get a relative error of 0.0000726882 for the energy. $\endgroup$ Feb 14 at 1:49
  • $\begingroup$ Oh yeah, 4Pi! I multiplied the conventional prefactor for the curvature energy (Helfrich). The bare integral has to be 4Pi. Thanks for fixing it! $\endgroup$
    – physgj
    Feb 14 at 21:02

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.