Working with meshes: curvature computation#

The goal of this notebook is to compute per-vertex curvature descriptors on triangular meshes: principal curvatures (k1, k2), mean curvature, Gaussian curvature, curvature anisotropy, and principal curvature tangent directions.

Curvatures are computed with a Python implementation of the finite-difference method from:

Rusinkiewicz, S. (2004). Estimating Curvatures and Their Derivatives on Triangle Meshes. Proceedings of the 2nd International Symposium on 3D Data Processing, Visualization, and Transmission (3DPVT).

Rusinkiewicz estimates local curvature from how vertex normals change across triangle edges. In this implementation, those normal differences are fit as face-level curvature tensors, transformed and averaged onto vertices, optionally smoothed, and diagonalized to obtain k1, k2, and their tangent directions.

Data: the input and output files for this tutorial can be downloaded here.

Table of contents#

From segmentation to marching-cubes mesh
Starting from a precomputed screened-Poisson mesh
Native visualization

Both routes rely on the same operations: cleanup, simplification, optional smoothing and normal refinement, curvature computation, and .vtp export with curvature fields.

Starting point	What it represents	Main use in this notebook
Segmentation → marching cubes	Boundary isosurface of the segmented volume (contains both surfaces)	Mesh generation from an `.mrc` segmentation
Precomputed mesh	A regular triangle mesh loaded from disk; here, a screened Poisson mesh of a single surface	Curvature analysis after external reconstruction (e.g. PyMeshLab)

[1]:

%load_ext autoreload
%autoreload 2

[2]:

import os
from cryocat.core import surface
from cryocat.analysis import visplot as vp
from IPython.display import Image

/Users/desislavaglushkova/cryoCAT_versions/260513_refactor/cryoCAT/cryocat/utils/ioutils.py:378: SyntaxWarning: invalid escape sequence '\$'
  ValueError: The format file_/$A/$B.txt does not contain any sequence of \$ followed by C.
/Users/desislavaglushkova/cryoCAT_versions/260513_refactor/cryoCAT/cryocat/utils/ioutils.py:384: SyntaxWarning: invalid escape sequence '\$'
  pattern = f"\$(?:{test_letter})+"
/Users/desislavaglushkova/cryoCAT_versions/260513_refactor/cryoCAT/cryocat/utils/ioutils.py:389: SyntaxWarning: invalid escape sequence '\$'
  f"The format {filename_format} does not contain any sequence of \$ followed by {test_letter}."

[3]:

# Define input path
parent_input_folder = "inputs/"
input_IMM_path = parent_input_folder + "2140_z150to400_IMM.mrc"

[4]:

# Define output path
output_path = "outputs/"
os.makedirs(output_path, exist_ok=True)

From segmentation to marching-cubes mesh#

surface.Mesh.from_mrc converts an .mrc segmentation into a triangle mesh with marching cubes. Marching cubes traces the isosurface of the segmentation, generating a mesh that contains both segmentation surfaces, instead of a single midplane.

Variable	Value here	Meaning
`input_path`	`input_IMM_path`	Segmentation file to read.
`transpose`	`True`	Uses the same axis convention as the rest of the notebook.
`voxel_size`	`0.7840`	Coordinate spacing; vertices are generated directly in nm.
`smooth_sigma`	`2.5`	Gaussian smoothing applied to the segmentation before marching cubes.
`level`	`0.5`	Isovalue for a binary segmentation.

mesh.units is set separately after loading. This controls the units reported for curvature values.

[5]:

# Generate a mesh from the segmentation
IMM_mesh = surface.Mesh.from_mrc(
    input_path=input_IMM_path,
    transpose=True,
    pixel_size=0.7840, # in nm
    smooth_sigma=2.5,
    level=0.5,
)

# Optionally refine normals after marching cubes
IMM_mesh = IMM_mesh.refine_normals(radius_hit=3.0, n_iter=1)

Applying Gaussian smoothing (sigma=2.5)...
Smoothed field range: [0.000, 0.964]
Running marching cubes (level=0.5, step_size=1)...
Refining normals (radius=3.0, batch=2000, iter=1, mask=all)
Iteration 1/1

100%|██████████| 915/915 [00:52<00:00, 17.58it/s]

Normal refinement complete.

[6]:

# Make sure to set and keep track for curvature computation
IMM_mesh.units = "nm"

1. Clean-up and simplify the mesh#

The raw marching-cubes mesh is large and may contain duplicate vertices, degenerate triangles, or more triangles than needed for curvature analysis.

Function	Variable	Value here	Meaning
`cleanup_mesh`	`simplify_mesh`	`True`	Run topology cleanup and quadric decimation.
`cleanup_mesh`	`target_number_of_triangles`	`500000`	Target face count after simplification. This mesh contains both membrane sides.
`smooth`	`iterations`	`2`	Taubin smoothing passes.
`smooth`	`repair_nonfinite`	`True`	Remove vertices with NaN/Inf coordinates before smoothing if needed.
`smooth`	`recompute_normals`	`True`	Recompute normals after smoothing changes the geometry.
`refine_normals`	`radius_hit`	`3.0`	Neighborhood radius for normal averaging, in mesh units.
`refine_normals`	`n_iter`	`1`	Number of refinement passes.

Taubin smoothing reduces small-scale triangulation noise while limiting the shrinkage caused by ordinary Laplacian smoothing.

[7]:

# The marching-cubes mesh contains both segmentation surfaces (take into account when choosing the number of faces)

IMM_mesh_simplified = IMM_mesh.cleanup_mesh(
    simplify_mesh=True,
    target_number_of_triangles=500000
)

Recomputing normals after mesh cleanup...
Mesh after cleaning: 250905 vertices, 499999 faces

[8]:

# Taubin smoothing filter to prevent shrinkage
IMM_mesh_simplified.smooth(
    iterations=2,
    repair_nonfinite=True,
    recompute_normals=True
)
# Optionally refine normals to smooth artifacts from simplification
IMM_mesh_simplified = IMM_mesh_simplified.refine_normals(radius_hit=3.0, n_iter=1)

Refining normals (radius=3.0, batch=2000, iter=1, mask=all)
Iteration 1/1

100%|██████████| 126/126 [00:06<00:00, 20.95it/s]

Normal refinement complete.

1. Compute curvatures#

compute_curvatures stores curvature values directly on the mesh object.

Variable	Value here	Meaning
`smoothing_iterations`	`3`	Number of smoothing passes applied to curvature tensors before extracting `k1` and `k2`.
`smoothing_kernel_rings`	`2`	Neighborhood size for tensor smoothing.

Implementation outline, following Rusinkiewicz (2004):

Step	What the code does
Mesh geometry	Computes edge vectors, edge lengths, face areas, and face normals.
Vertex frames	Builds vertex normals, vertex areas, and local tangent coordinate systems.
Face tensors	Fits a curvature tensor per face from normal differences across triangle edges.
Vertex tensors	Transforms face tensors into vertex frames and averages them with area-based weights.
Optional smoothing	Smooths curvature tensors over mesh neighborhoods.
Principal values	Diagonalizes each 2×2 tensor to get `k1`, `k2`, and principal directions.

Stored field	Definition / interpretation
`k1`, `k2`	Principal curvatures, with `k1 >= k2` by convention.
`mean_curvature`	`(k1 + k2) / 2`.
`gaussian_curvature`	`k1 * k2`.
`curvature_anisotropy`	`abs(k1 - k2)`.
`principal_direction_1`, `principal_direction_2`	Tangent directions associated with the principal curvatures.

Curvature sign depends on mesh normal orientation, so compare signs only when normals are oriented consistently.

[9]:

IMM_mesh_cuvatures = IMM_mesh_simplified.compute_curvatures(smoothing_iterations=3, smoothing_kernel_rings=2)

Computing curvatures (units: nm^(-1) (principal/mean), nm^(-2) (Gaussian))

1. Save meshes without/with curvature fields#

Geometry-only meshes can be saved as .ply or .vtp. Meshes containing stored curvatures should be saved as .vtp files, so named per-vertex arrays are preserved.

Variable	Value here	Meaning
`output_path`	`...curvatures.vtp`	File to write.
`format`	`"vtp"`	Required when `include_curvatures=True`.
`include_curvatures`	`True`	Save curvature scalar fields and principal direction vectors.

[10]:

# Save mesh generated from marching cubes
IMM_mesh.save(
    output_path = output_path + "IMM_mesh_mc.ply",
    format="ply",
    include_curvatures=False # not computed at this stage
)

# Or as a .vtp file
IMM_mesh.save(
    output_path = output_path + "IMM_mesh_mc.vtp",
    format="vtp",
    include_curvatures=False # not computed at this stage
)

# Save simplified mesh to check whether it looks reasonable
IMM_mesh_simplified.save(
    output_path = output_path + "IMM_mesh_simplified_mc.ply",
    format="ply", # or as .vtp
    include_curvatures=False # not computed at this stage
)

Saved mesh to outputs/IMM_mesh_mc.ply
  Vertices: 250,905
  Faces: 499,999
Saved mesh to outputs/IMM_mesh_mc.vtp
  Format: VTP (VTK PolyData)
  Vertices: 250,905
  Faces: 499,999
Saved mesh to outputs/IMM_mesh_simplified_mc.ply
  Vertices: 250,905
  Faces: 499,999

[11]:

# Save the mesh with computed curvatures
IMM_mesh_cuvatures.save(
    output_path = output_path + "IMM_mesh_simplified_mc_curvatures.vtp",
    format="vtp",
    include_curvatures=True
)

Saved mesh with curvatures to outputs/IMM_mesh_simplified_mc_curvatures.vtp
  Format: VTP (VTK PolyData)
  Vertices: 250,905
  Faces: 499,999
  Scalar fields (8):
    - mean_curvature: [-3.619927e-01, 1.178025e-01]
    - gaussian_curvature: [-3.513244e-02, 1.262909e-01]
    - k1: [-2.969925e-01, 2.047013e-01]
    - k2: [-4.308967e-01, 1.068652e-01]
    - curvature_anisotropy: [0.000000e+00, 3.757315e-01]
    - shape_index: [-9.993041e-01, 9.990319e-01]
    - curvedness: [0.000000e+00, 3.684921e-01]
    - shape_category: [-1, 8] (int codes)
  Vector fields (3): normals, principal_direction_1, principal_direction_2
  Curvature units: nm^(-1) (principal/mean), nm^(-2) (Gaussian)

[11]:

{'format': 'VTP',
 'path': 'outputs/IMM_mesh_simplified_mc_curvatures.vtp',
 'n_vertices': 250905,
 'n_faces': 499999,
 'coordinate_units': 'nm',
 'scalar_fields': ['mean_curvature',
  'gaussian_curvature',
  'k1',
  'k2',
  'curvature_anisotropy',
  'shape_index',
  'curvedness',
  'shape_category'],
 'vector_fields': ['normals',
  'principal_direction_1',
  'principal_direction_2'],
 'statistics': {'mean_curvature': {'min': -0.361992672032641,
   'max': 0.11780253600558332,
   'mean': -0.019930810261905647},
  'gaussian_curvature': {'min': -0.03513244200188554,
   'max': 0.12629093398246224,
   'mean': -6.139591337126725e-05},
  'k1': {'min': -0.2969924680839255,
   'max': 0.20470130163895328,
   'mean': 0.010249963487517237},
  'k2': {'min': -0.4308966677286017,
   'max': 0.10686517438285724,
   'mean': -0.050111584011328536},
  'curvature_anisotropy': {'min': 0.0,
   'max': 0.37573153168160445,
   'mean': 0.06036154749884578},
  'shape_index': {'min': -0.9993041309801391,
   'max': 0.9990319452824967,
   'mean': -0.16893865115779508},
  'curvedness': {'min': 0.0,
   'max': 0.36849213726781227,
   'mean': 0.04456201347217928}}}

Starting from a precomputed screened-Poisson mesh#

The second route starts from a single-surface triangle mesh that was generated using the screened Poisson reconstruction (outside this notebook) and saved as a .ply file. This mesh samples ≈ the membrane midplane surface rather than the inner and outer segmentation boundaries.

Source for screened Poisson reconstruction: MeshLab / PyMeshLab-style workflow

Ranzuglia, G., Callieri, M., Dellepiane, M., Cignoni, P., & Scopigno, R. (2012). Efficient and Flexible Sampling with Blue Noise Properties of Triangular Meshes. IEEE Transactions on Visualization and Computer Graphics, 18(6), 914–924.

[12]:

input_IMM_mesh_screenedPoisson_path = parent_input_folder + "IMM_mesh_screenedPoisson.ply"

[13]:

IMM_mesh_screenedPoisson = surface.Mesh.read(input_path=input_IMM_mesh_screenedPoisson_path)

Loading triangle mesh: inputs/IMM_mesh_screenedPoisson.ply
Loaded: 4224564 vertices, 8420888 faces

Coordinate units not specified!
   Set with: mesh.units = 'nm'  (or 'angstrom', 'pixel', etc.)
   This determines curvature units (e.g., nm⁻¹ for nanometers)

[14]:

IMM_mesh_screenedPoisson.units = 'nm'

2. Clean-up and simplify the mesh#

Function	Variable	Value here	Meaning
`cleanup_mesh`	`simplify_mesh`	`True`	Run cleanup and quadric decimation.
`cleanup_mesh`	`target_number_of_triangles`	`300000`	Target face count for the midplane mesh.
`smooth`	`iterations`	`2`	Taubin smoothing passes.
`smooth`	`repair_nonfinite`	`True`	Remove NaN/Inf vertices before smoothing if needed.
`smooth`	`recompute_normals`	`True`	Recompute normals after smoothing.
`refine_normals`	`radius_hit`	`3.0`	Neighborhood radius for normal averaging.
`refine_normals`	`n_iter`	`1`	Number of normal-refinement passes.

[15]:

# This mesh has close to 8M faces, so we simplify it to 300k (a bit more than 2x fewer vertices than the marching cubes mesh)
IMM_mesh_screenedPoisson_simplified = IMM_mesh_screenedPoisson.cleanup_mesh(
    simplify_mesh=True,
    target_number_of_triangles=300000
)

Recomputing normals after mesh cleanup...
Mesh after cleaning: 159068 vertices, 299999 faces

[16]:

# Taubin smoothing filter to prevent shrinkage
IMM_mesh_screenedPoisson_simplified.smooth(
    iterations=2,
    repair_nonfinite=True,
    recompute_normals=True
)

# Optionally refine normals to smooth artifacts from simplification
IMM_mesh_screenedPoisson_simplified = IMM_mesh_screenedPoisson_simplified.refine_normals(radius_hit=3.0, n_iter=1)

Removed 68 non-finite vertices and 0 affected faces
Refining normals (radius=3.0, batch=2000, iter=1, mask=all)
Iteration 1/1

100%|██████████| 80/80 [00:03<00:00, 20.43it/s]

Normal refinement complete.

[17]:

IMM_mesh_screenedPoisson_simplified.save(
    output_path = output_path + "IMM_mesh_screenedPoisson_simplified.ply",
    format="ply",
    include_curvatures=False
)

Saved mesh to outputs/IMM_mesh_screenedPoisson_simplified.ply
  Vertices: 159,000
  Faces: 299,999

Remove small disconnected components#

Reconstructed meshes can contain small isolated pieces from segmentation noise, sparse point-cloud regions, or reconstruction artifacts.

Variable	Value here	Meaning
`min_component_size`	`500`	Keep connected components with at least this many vertices.

Choose this threshold conservatively and inspect the result. The goal is to remove obvious artifacts without deleting real disconnected structures.

[18]:

IMM_mesh_screenedPoisson_minuscomponents = IMM_mesh_screenedPoisson_simplified.remove_disconnected_components(
    min_component_size=800
)

Removing disconnected components smaller than 800 vertices...
Found 11 connected components
Removing 7 small components

[19]:

IMM_mesh_screenedPoisson_simplified.save(
    output_path = output_path + "IMM_mesh_screenedPoisson_simplified_minuscomponents.ply",
    format="ply",
    include_curvatures=False
)

Saved mesh to outputs/IMM_mesh_screenedPoisson_simplified_minuscomponents.ply
  Vertices: 157,519
  Faces: 297,599

Erode boundary vertices#

Open boundaries can give unreliable curvature estimates as boundary vertices have incomplete neighborhoods.

Variable	Value here	Meaning
`erosion_layers`	`2`	Remove boundary vertices, then one additional neighbor layer.
`min_edge_faces`	default `2`	Edges with fewer adjacent faces are treated as boundary edges.

This is useful for cropped or open screened-Poisson meshes. Each erosion layer removes real mesh area, so inspect the output before computing final curvatures.

[20]:

IMM_mesh_screenedPoisson_eroded = IMM_mesh_screenedPoisson_simplified.remove_boundary_vertices(erosion_layers=2)

Removing boundary vertices (erosion_layers=2, min_edge_faces=2)...
  Found 17489 boundary edges
  Layer 1: Marked 17488 boundary vertices
  Layer 2: Marked 6377 additional vertices
  Removed 23865 vertices (15.2%)
  Mesh now has 133654 vertices, 260614 faces

[21]:

IMM_mesh_screenedPoisson_eroded.save(
    output_path = output_path + "IMM_mesh_screenedPoisson_simplified_eroded.ply",
    format="ply",
    include_curvatures=False)

Saved mesh to outputs/IMM_mesh_screenedPoisson_simplified_eroded.ply
  Vertices: 133,654
  Faces: 260,614

2. Compute curvatures#

[22]:

IMM_mesh_screenedPoisson_curvatures = IMM_mesh_screenedPoisson_eroded.compute_curvatures(smoothing_iterations=3, smoothing_kernel_rings=2)

Computing curvatures (units: nm^(-1) (principal/mean), nm^(-2) (Gaussian))

2. Save meshes without/with curvature fields#

[23]:

IMM_mesh_screenedPoisson_curvatures.save(
    output_path = output_path + "IMM_mesh_screenedPoisson_curvatures.vtp",
    format="vtp",
    include_curvatures=True
)

Saved mesh with curvatures to outputs/IMM_mesh_screenedPoisson_curvatures.vtp
  Format: VTP (VTK PolyData)
  Vertices: 133,654
  Faces: 260,614
  Scalar fields (8):
    - mean_curvature: [-5.168882e-01, 1.803246e-01]
    - gaussian_curvature: [-8.929300e-02, 2.630900e-01]
    - k1: [-4.830991e-01, 5.054894e-01]
    - k2: [-8.593560e-01, 1.324550e-01]
    - curvature_anisotropy: [0.000000e+00, 7.606622e-01]
    - shape_index: [-9.995737e-01, 9.983791e-01]
    - curvedness: [0.000000e+00, 6.200464e-01]
    - shape_category: [-1, 8] (int codes)
  Vector fields (3): normals, principal_direction_1, principal_direction_2
  Curvature units: nm^(-1) (principal/mean), nm^(-2) (Gaussian)

[23]:

{'format': 'VTP',
 'path': 'outputs/IMM_mesh_screenedPoisson_curvatures.vtp',
 'n_vertices': 133654,
 'n_faces': 260614,
 'coordinate_units': 'nm',
 'scalar_fields': ['mean_curvature',
  'gaussian_curvature',
  'k1',
  'k2',
  'curvature_anisotropy',
  'shape_index',
  'curvedness',
  'shape_category'],
 'vector_fields': ['normals',
  'principal_direction_1',
  'principal_direction_2'],
 'statistics': {'mean_curvature': {'min': -0.5168881731607495,
   'max': 0.18032457309526292,
   'mean': -0.0063663109102357835},
  'gaussian_curvature': {'min': -0.08929300026791702,
   'max': 0.2630900188100896,
   'mean': 0.0001417383916525116},
  'k1': {'min': -0.48309905571458106,
   'max': 0.5054893754853778,
   'mean': 0.013146248638606896},
  'k2': {'min': -0.8593559905182536,
   'max': 0.1324549945010008,
   'mean': -0.025878870459078464},
  'curvature_anisotropy': {'min': 0.0,
   'max': 0.7606621993125067,
   'mean': 0.039025119097685364},
  'shape_index': {'min': -0.9995736627961141,
   'max': 0.9983791059761671,
   'mean': 0.09602621602739542},
  'curvedness': {'min': 0.0,
   'max': 0.6200464413889252,
   'mean': 0.03110502689852565}}}

Native visualization#

The meshes produced above can be visualized interactively in the notebook using plot_vtp_mesh from visplot.

Each mesh is rendered as a go.Mesh3d trace and colored by any stored per-vertex scalar field. The color_by parameter accepts any field saved in the .vtp file: "mean_curvature", "gaussian_curvature", "k1", "k2", "curvature_anisotropy", "shape_index", "curvedness", or "shape_category".

Field	Range	Sign	Suggested colorscale
`mean_curvature`, `gaussian_curvature`, `k1`, `k2`	signed	±	`"RdBu_r"` (symmetric)
`shape_index`	`[-1, 1]`	±	`"RdBu_r"` (symmetric)
`curvature_anisotropy`, `curvedness`	`>= 0`	unsigned	`"Viridis"`
`shape_category`	discrete codes	categorical	pass `categorical=True`

The colorscale is symmetric around zero by default ("RdBu_r"), which works well for signed curvatures such as mean/Gaussian curvature and shape_index. For unsigned fields (curvature_anisotropy, curvedness), use colorscale="Viridis". Use colorscale_range to clip outliers or to put two meshes on a shared scale.

[24]:

fig = vp.plot_vtp_mesh(
    output_path + "IMM_mesh_screenedPoisson_curvatures.vtp",
    color_by="mean_curvature",
    colorscale_range=(-0.06, 0.06),
    opacity=1.0,
    title="IMM – mean curvature",
)
#fig.show() # uncomment to render
fig.write_image("plots/IMM_fullmesh_meancurvature.png", width=1200, height=800)

[25]:

Image("plots/IMM_fullmesh_meancurvature.png")

[25]:

../../_images/tutorials_compute_curvatures_compute_curvatures_36_0.png

[26]:

fig = vp.plot_vtp_mesh(
    output_path + "IMM_mesh_screenedPoisson_curvatures.vtp",
    color_by="gaussian_curvature",
    colorscale_range=(-0.003, 0.003),
    opacity=1.0,
    title="IMM – Gaussian curvature",
)
#fig.show() # uncomment to render
fig.write_image("plots/IMM_fullmesh_Gaussiancurvature.png", width=1200, height=800)

[27]:

Image("plots/IMM_fullmesh_Gaussiancurvature.png")

[27]:

../../_images/tutorials_compute_curvatures_compute_curvatures_38_0.png

[28]:

fig = vp.plot_vtp_mesh(
    output_path + "IMM_mesh_screenedPoisson_curvatures.vtp",
    color_by="shape_index",
    colorscale_range=(-1.0, 1.0),
    opacity=1.0,
    title="IMM – shape index",
)
#fig.show() # uncomment to render
fig.write_image("plots/IMM_fullmesh_shapeindex.png", width=1200, height=800)

[30]:

Image("plots/IMM_fullmesh_shapeindex.png")

[30]:

../../_images/tutorials_compute_curvatures_compute_curvatures_40_0.png

[33]:

fig = vp.plot_vtp_mesh(
    output_path + "IMM_mesh_screenedPoisson_curvatures.vtp",
    color_by="curvedness",
    colorscale="BuPu",
    colorscale_range=(0.0, 0.06),
    opacity=1.0,
    title="IMM – curvedness",
)
#fig.show() # uncomment to render
fig.write_image("plots/IMM_fullmesh_curvedness.png", width=1200, height=800)

[34]:

Image("plots/IMM_fullmesh_curvedness.png")

[34]:

../../_images/tutorials_compute_curvatures_compute_curvatures_42_0.png

[ ]: