Measuring surface curvature

[marcomusy]

[forwarded message]
Can I calculate the principal curvatures of the mesh and their directions (see the image)? Also, is there another method to measure curvature besides using mesh.compute_curvature(method=0), where the measurement is independent of quality or fineness of the mesh?

Thanks for the help. I am very happy with the vedo package.

[marcomusy]

I don't think there is a method that is completely independent of the mesh precision, but you can make this dependence negligible if the mesh is sufficiently high res (your mesh looks good in this sense).
Here is an example where you can choose the level of averaging (first part of the code). To run the code, grab the latest dev version

pip install -U git+https://github.com/marcomusy/vedo.git

then run

import numpy as np
from vedo import Mesh, Plotter, show, Points, Point, Arrows
from vedo import project_point_on_variety, progressbar, print_histogram

msh = Mesh("data/smoothened-mesh.stl").scale(0.01)
msh.compute_normals()
adlist = msh.compute_adjacency()

# Example of computing curvature at a point ##################
pid = 11131  # 11131 5759
pts = msh.points
ids = msh.find_adjacent_vertices(pid, depth=8)
bpts = msh.points[ids]
p = pts[pid]
q, res = project_point_on_variety(
    p, bpts, degree=2, compute_surface=True, compute_curvature=True
)
show(Points(pts[ids], r=5), Point(p), res[0], msh.clone().alpha(0.2)).close()


##############################################################
# Example of computing curvature at all points in the mesh 
# by custom fitting a quadratic surface
curvs_g = []
curvs_m = []
# Compute curvature at the point pid
for i in progressbar(range(msh.npoints)):
    ids = msh.find_adjacent_vertices(i, depth=3, adjacency_list=adlist)
    bpts = msh.points[ids]
    _, res = project_point_on_variety(pts[i], bpts, degree=2, compute_curvature=True)
    curvs_g.append(res[4])
    curvs_m.append(res[5])
msh.pointdata["Gauss_Curvature"] = curvs_g
msh.pointdata["Mean_Curvature"]  = curvs_m

print_histogram(curvs_g, bins=12, title="Gaussian", logscale=0, vrange=(-2, 2))
print_histogram(curvs_m, bins=12, title="Mean ", logscale=0, vrange=(0, 2))

msh.lighting("glossy")
msh.cmap('coolwarm', "Gauss_Curvature", vmin=-2, vmax=2).add_scalarbar()
show(msh, "Gauss Curvature").close()
msh.cmap("Reds", "Mean_Curvature", vmin=0, vmax=2).add_scalarbar()
show(msh, "Mean Curvature").close()


############################################################## default examples
msh = Mesh("data/smoothened-mesh.stl").scale(0.01)
msh.compute_normals()
msh1 = msh.clone().compute_curvature(method=0)
arr = msh1.pointdata["Gauss_Curvature"]
print_histogram(arr, bins=12, title="Gaussian Curvature", logscale=0, vrange=(-2, 2))
msh1.cmap("coolwarm", "Gauss_Curvature", vmin=-2, vmax=2).add_scalarbar()

msh2 = msh.clone().compute_curvature(method=1)
arr = msh2.pointdata["Mean_Curvature"]
print_histogram(arr, bins=12, title="Mean Curvature", logscale=0, vrange=(-2, 2))
msh2.cmap("coolwarm_r", "Mean_Curvature", vmin=-2, vmax=2).add_scalarbar()

msh3 = msh1.clone().lighting("off")
vgrad = msh3.gradient("Gauss_Curvature")
# normalize the gradient vector
vgrad /= np.linalg.norm(vgrad, axis=1)[:, np.newaxis] * 10
arrows = Arrows(msh3.points, msh3.points+vgrad, c="black")
msh3.pointdata.select("Gauss_Curvature")
isolines = msh3.isolines(n=10, vmin=-1, vmax=1).c("black")

plt = Plotter(N=4, axes=1)
plt.at(0).show(msh.lw(1), "Original shape")
plt.at(1).show(msh1, "Gaussian curvature")
plt.at(2).show(msh2, "Mean curvature")
plt.at(3).show(msh3, isolines, arrows, "Gradient")
plt.interactive().close()

Hope this helps

[marcomusy]

Btw, I just fixed a small bug - so please update vedo (eg. with git pull).

I also created a new example, as I think this might interesting for others too (@jo-mueller )

[jo-mueller]

Super cool! @marcomusy Maybe one comment for the above - am I right in the assumption that the depth keyword in this command

ids = msh.find_adjacent_vertices(i, depth=3, adjacency_list=adlist)

returns all 3-hop neighbors, so to say, for a queried vertex? I would find it a bit cleaner to search for neighbors by a given radius. Depending on how the mesh is constructed (i.e., if it has scarce vertices in flat areas), this could influence the outcome in strange ways? I guess my "safe route" to take at this is to try to smooth/decimate the mesh into triangles of regular size and then calculating curvature in a neighborhood specifiec by radius rather than by neighborhood.

Then again, maybe taking far-away neighbors for scarce vertices in flat areas is actually really desirable? 🤔 I could see this being legitimate, too.

[marcomusy]

Hi Johannes, yes that would be cleaner but also computationally more expensive because you then need to compute geodesics...
I would need to quantify how much it takes to compute it. But it's true, it would make sense at least to have a max_radius keyword.

PS: in fact maybe specifying a desired nr of nieghbors is even a better and cheaper option.

[jo-mueller]

because you then need to compute geodesics...

Yosh, that's true. At least in geometrically intricate regions I wouldn't count much on the euclidian distance to be a good measure of neighborhood. Maybe a combination of both would be reasonable - if that's what you are suggesting? I.e., obtain a group of neighboring vertices and cap these vertices by a max_radius (geodesic/euclidian) distance?

[marcomusy]

Thinking about it I believe that specifying a desired number of points is probably the best option because tipically meshes that are locally more dense of vertices in some regions are so because higher precision is needed in those locations.
So in a sense it takes automatically in account the complexity of the surface...