Curvature

Accuracy

The following benchmark compares the accuracy of the principal curvatures for three meshes (a sphere, a sphere with improved quality mesh, and a torus) as computed by three mymesh methods (QuadFit(), CubicFit(), AnalyticalCurvature()) and the methods of two other libraries: VTK (via pyvista) and trimesh. For the unit sphere, both principal curvatures are theoretically equal to 1, and for the torus with minor radius 1, the maximum principal curvature is theoretically equal to 1 while the minimum principal curvature varies spatially.

The results highlight two points:

1. the sensitivity of mesh-based curvature to the quality of the mesh (as evidenced by the more accurate results for the smoothed sphere mesh), and

2. that the methods implemented in mymesh are comparable or superior to the alternatives that were tested (at least in these test cases).

The best performing method in all cases was the implicit function-based curvature (AnalyticalCurvature()), however this is only possible when an implicit function representation of the surface exists. The cubic surface fitting method (CubicFit(), [GI04]) is generally superior to other mesh-based approaches, but may be more susceptible to errors due to the low-quality elements.

See also: Curvature, curvature, Curvature Analysis.

import mymesh
from mymesh import curvature, implicit
import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import trimesh
import subprocess

def test_curvature(M, func=None):
    # Function to test different methods of curvature calculation

    # mymesh - quadratic fitting
    k1_quadratic, k2_quadratic = curvature.QuadFit(M.NodeCoords,
                                            M.NodeConn,
                                            M.NodeNeighbors,
                                            M.NodeNormals)

    # mymesh - cubic fitting
    k1_cubic, k2_cubic = curvature.CubicFit(M.NodeCoords,
                                            M.NodeConn,
                                            M.NodeNeighbors,
                                            M.NodeNormals)


    # vtk (via pyvista)
    k1_vtk = M.to_pyvista().extract_surface().curvature(curv_type='maximum')
    k2_vtk = M.to_pyvista().extract_surface().curvature(curv_type='minimum')
    # trimesh
    r = 0.05
    M_trimesh = trimesh.base.Trimesh(M.NodeCoords, M.NodeConn)
    K_trimesh = trimesh.curvature.discrete_gaussian_curvature_measure(M_trimesh, M_trimesh.vertices, r) / trimesh.curvature.sphere_ball_intersection(1, r)
    H_trimesh = trimesh.curvature.discrete_mean_curvature_measure (M_trimesh, M_trimesh.vertices, r) / trimesh.curvature.sphere_ball_intersection(1, r)
    k1_trimesh = H_trimesh + np.sqrt(np.maximum(H_trimesh**2-K_trimesh,0))
    k2_trimesh = H_trimesh - np.sqrt(np.maximum(H_trimesh**2-K_trimesh,0))

    if func is not None:
        # mymesh - implicit
        k1_implicit, k2_implicit, K_implicit, H_implicit = curvature.AnalyticalCurvature(func, M.NodeCoords)

        return (k1_quadratic, k2_quadratic), \
                (k1_cubic, k2_cubic), \
                (k1_vtk, k2_vtk), \
                (k1_trimesh, k2_trimesh), \
                (k1_implicit, k2_implicit)
    subprocess.run('echo "debug g"', shell=True)
    return (k1_quadratic, k2_quadratic), \
                (k1_cubic, k2_cubic), \
                (k1_vtk, k2_vtk), \
                (k1_trimesh, k2_trimesh)

(Source code)