Particle–mesh intersections and local curvatures

This notebook is a toy example in which we analyse whether the inner mitochondrial membrane (IMM) directly beneath an ATP synthase differs in curvature from the surrounding membrane. To do so, we:

1. Find where particles intersect the membrane (hit triangles or vertices).

2. Grow disk-shaped regions on the mesh around each intersection, sized to the approximate in-plane extent of an ATP synthase.

3. Compare mean and Gaussian curvature in two regions: an inner disk encompassing the ATP synthase footprint (≤ ~9 nm) vs an outer annulus (~9–18 nm).

These radii were chosen based on the approximate radius of the ATP synthase base (~9 nm from STA, measured in ChimeraX) and a surrounding band at twice that radius.

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

Table of contents

1. Setup

2. Load the data

   • IMM mesh

   • ATP synthase motive list

   • Visualize IMM mesh

3. Particles to mesh

   • Cast rays

   • Grow regions around ray hits

   • Extract and save membrane patches: ray hits

   • Filter particles: ray hits

   • Visualize patches: ray hits

4. Mesh to particles

   • Build particle cloud

   • Grow regions around vertex hits

   • Extract and save membrane patches: vertex hits

   • Filter particles: vertex hits

   • Visualize patches: vertex hits

5. Per-patch curvature distributions

---

Data sources

All example data come from the public Chlamy dataset resources linked below.

Dataset	What it is
IMM mesh	Screened-Poisson surface of the IMM with precomputed curvatures (see compute_curvatures tutorial)
ATP synthase motive list	Particle positions and orientations from the Chlamy annotations repository

IMM. MemBrain v2 segmentation (Lamm et al., 2026) on a cryoCARE-denoised (Buchholz et al., 2018) C. reinhardtii tomogram (EMPIAR-11830); screened Poisson mesh (Kazhdan & Hoppe, 2013) from an external MeshLab / PyMeshLab workflow; curvatures (Rusinkiewicz, 2004) as in the compute_curvatures tutorial.

ATP synthase. Particle coordinates from ChlamyAnnotations, associated with (Kelley et al., 2025, Mol. Cell).

Representations

We use ``PleomorphicSurface`` as the notebook-facing API that delegates methods to the right underlying geometric representation:

• The IMM is a mesh (is_mesh=True) with curvature fields on vertices/triangles.

• ATP synthases can be represented as rays for particle→mesh queries or wrapped as an ``OrientedPointCloud`` for mesh→particle queries.

Two complementary query directions

Direction	Query	Seeds on IMM	Question
Particles → mesh	Rays from each particle toward the membrane	Hit triangles	Which membrane triangle does this particle meet along its casting direction?
Mesh → particles	Rays from mesh vertices toward the particle cloud	Hit vertices; regions expanded to triangles	Which membrane vertices are close to any particle along the surface normal?

Use the same distance cutoffs and surface_radii in both directions when you want comparable inner/outer bands.

Setup

%load_ext autoreload
%autoreload 2

import numpy as np
import pandas as pd

from cryocat.utils import geom
from cryocat.core import cryomotl
from cryocat.core import surface
from cryocat.analysis import assembly
from cryocat.analysis import visplot as vp
from IPython.display import Image

Load the data

IMM mesh

The IMM surface is loaded with ``PleomorphicSurface.read`` and method="mesh_curvatures". That reads a VTP triangle mesh and attaches principal, mean, and Gaussian curvature fields produced in the compute_curvatures notebook. Set the units at load time so distances and curvature units stay consistent.

After loading, IMM_surface.is_mesh should be True.

# General paths
input_path = "inputs/"
output_path = "outputs/"

IMM_mesh_path = input_path + "IMM_mesh_screenedPoisson_curvatures.vtp"
OMM_mesh_path = input_path + "OMM_mesh_screenedPoisson_curvatures.vtp"

# Read in mesh with curvatures and directly wrap as a PleomorphicSurface instance
IMM = assembly.PleomorphicSurface.read(
    IMM_mesh_path,
    method="mesh_curvatures",
    units="nm"
)

Loading triangle mesh with curvatures: inputs/IMM_mesh_screenedPoisson_curvatures.vtp
Loaded: 133654 vertices, 260614 faces
Coordinate units: nm
Curvatures will be in: nm^(-1) (principal/mean), nm^(-2) (Gaussian)
Curvature data loaded from file
  Principal curvatures: k1=[-4.830991e-01, 5.054894e-01], k2=[-8.593560e-01, 1.324550e-01]
  Mean curvature: [-5.168882e-01, 1.803246e-01]
  Gaussian curvature: [-8.929300e-02, 2.630900e-01]

print(f"IMM.is_mesh: {IMM.is_mesh}")

IMM.is_mesh: True

ATP synthase motive list

Particle coordinates are rescaled by the pixel size to match the mesh coordinate system. Adjust this factor for your data.

Particle normals are computed as the z-axis in the motive list frame.

atp_motl_path = input_path + "atp_bin4_trimmed_cleaned.em"

# Load particle motive list
atp_synthase_motl = cryomotl.Motl.load(atp_motl_path)

# Get particle coordinates and re-scale with respect to the mesh (this depends on your data)

pixel_size = 0.196 # in nm to match the pixel size of the meshes
atp_synthase_xyz = atp_synthase_motl.get_coordinates() * pixel_size

# Get z-axis normals from the motive list

atp_synthase_rotations = atp_synthase_motl.get_rotations()
atp_synthase_normals = atp_synthase_rotations.apply([0, 0, 1])

Visualize IMM mesh

As a sanity check, visualize the full IMM mesh colored by Gaussian curvature before running any intersection queries. The full motive list is shown here; only intersecting particles will be overlaid on the patch visualizations below.

fig = vp.plot_vtp_mesh(
    IMM_mesh_path,
    color_by="gaussian_curvature",
    colorscale_range=(-0.003, 0.003),
    opacity=1.0,
    overlay_points=atp_synthase_xyz,
    overlay_names=["ATP synthase"],
    overlay_colors=["gray"],
    overlay_marker_size=3,
    title="IMM – Gaussian curvature + ATP synthase positions",
)
#fig.show()
fig.write_image("plots/IMM_fullmesh_Gaussiancurvature+ATPsynthases.png", width=1200, height=800)

Image("plots/IMM_fullmesh_Gaussiancurvature+ATPsynthases.png")

Particles to mesh

Rays start at each ATP synthase position and shoot in the −z normal direction toward the IMM mesh. When one_hit_per_target=True, the triangle with the shortest ray source-to-mesh distance is stored.

Ray direction: In this dataset, particle orientations were first inspected in ChimeraX. Casting along −z matches how they are oriented relative to the membrane. Verify for your data.

Step	Method	Purpose
Build rays	geom.construct_rays	(N, 6) array: origin + direction per particle
Intersect	IMM.ray_intersections	Open3D ray cast against the mesh
Analyze	IMM.intersection_data	Hit table, radial regions, curvature summary
Extract	IMM.extract_region	Sub-meshes for visualization / export

# Construct rays from particles in minus z-direction

atp_synthase_rays = geom.construct_rays(
    points=atp_synthase_xyz,
    normals=atp_synthase_normals,
    reverse_direction=True
)

Cast rays

# Compute intersections from the particles to the mesh in the direction of the rays

cast_rays_on_IMM_mesh = IMM.ray_intersections(
    rays=atp_synthase_rays,
    one_hit_per_target=True, # stores the triangle with the shortest distance from the point source
)

Grow regions around ray hits

intersection_data turns raw ray triangle hits into a summary DataFrame and optionally computes curvatures for each grown region.

Parameter	Value here	Meaning
query_type	"ray"	Parse output of ray_intersections
max_distance_source_target	20	Drop hits farther than 20 nm (empirical; set after visual inspection)
source_id_name	"particle_id"	Ray / particle index
target_id_name	"triangle_id"	Intersected triangle on the mesh
surface_radii	[9.0, 18.0]	Grow regions around each intersected triangle: inner disk ≤ 9 nm; annulus 9–18 nm
include_curvatures	True	Attach mean / Gaussian curvature at each intersected triangle; summarize per region

Why the radii? ~9 nm approximates the radius of the ATP synthase base from STA. 18 nm defines an outer band for the surrounding membrane. Regions are grown in the mesh metric from the seed triangles.

Inspect ``region_summary``: the main curvature comparison table for the particles→mesh path.

# Wrapper to get intersection data

particle_to_mesh_data = IMM.intersection_data(
    result=cast_rays_on_IMM_mesh,
    query_type="ray",
    max_distance_source_target=20, # nm, same as mesh
    source_id_name="particle_id",
    target_id_name="triangle_id",
    surface_radii=[10.0, 20.0], # nm, same as mesh; what area sizes to grow on the mesh
    include_curvatures=True,
)

particle_to_mesh_hits = particle_to_mesh_data["hits"] # get intersected triangles
particle_to_mesh_data["region_summary"] # get summary of regions on the mesh, chosen with surface_radii

	region	n_triangles	mean_curvature_mean	mean_curvature_median	gaussian_curvature_mean	gaussian_curvature_median
0	hit triangles	360	-0.019621	-0.030229	0.000497	0.000641
1	r <= 10 nm	77874	-0.021095	-0.030045	0.000545	0.000630
2	r <= 20 nm	112150	-0.019982	-0.025765	0.000212	0.000365
3	10 < r <= 20 nm	34276	-0.017452	-0.016018	-0.000543	-0.000091

Extract and save membrane patches: ray hits

extract_region takes a triangle index array from the summary DataFrame and returns a new ``PleomorphicSurface`` sub-mesh with curvature fields preserved.

Region key	Output file	Content
"r <= 10 nm"	.vtp mesh with curvatures	Membrane directly under particles
"10 < r <= 20 nm"	.vtp mesh with curvatures	Surrounding annulus

# Grow regions starting from the intersected triangle
# Region within an ATP synthase particle

IMM_inner_patch = IMM.extract_region(
    particle_to_mesh_data["regions"]["r <= 10 nm"],
    element="triangles",
)

# Surrounding annulus
IMM_outer_patch = IMM.extract_region(
    particle_to_mesh_data["regions"]["10 < r <= 20 nm"],
    element="triangles",
)

# Save the extracted meshes

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

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

Saved mesh with curvatures to outputs/IMM_below10nm_patch_particlestomesh.vtp
  Format: VTP (VTK PolyData)
  Vertices: 42,192
  Faces: 77,874
  Scalar fields (8):
    - mean_curvature: [-4.362034e-01, 1.197936e-01]
    - gaussian_curvature: [-3.727752e-02, 1.177193e-01]
    - k1: [-1.668447e-01, 2.111054e-01]
    - k2: [-7.482412e-01, 5.034715e-02]
    - curvature_anisotropy: [2.241403e-04, 7.015738e-01]
    - shape_index: [-9.928864e-01, 9.974568e-01]
    - curvedness: [6.154567e-04, 5.301145e-01]
    - shape_category: [0, 8] (int codes)
  Vector fields (3): normals, principal_direction_1, principal_direction_2
  Curvature units: nm^(-1) (principal/mean), nm^(-2) (Gaussian)
Saved mesh with curvatures to outputs/IMM_10nmto20nm_patch_particlestomesh.vtp
  Format: VTP (VTK PolyData)
  Vertices: 21,268
  Faces: 34,276
  Scalar fields (8):
    - mean_curvature: [-3.980009e-01, 1.197936e-01]
    - gaussian_curvature: [-4.719297e-02, 8.088003e-02]
    - k1: [-1.435743e-01, 2.111054e-01]
    - k2: [-7.760181e-01, 6.802757e-02]
    - curvature_anisotropy: [3.295864e-05, 7.606622e-01]
    - shape_index: [-9.995737e-01, 9.967550e-01]
    - curvedness: [6.154567e-04, 5.489096e-01]
    - shape_category: [0, 8] (int codes)
  Vector fields (3): normals, principal_direction_1, principal_direction_2
  Curvature units: nm^(-1) (principal/mean), nm^(-2) (Gaussian)

{'format': 'VTP',
 'path': 'outputs/IMM_10nmto20nm_patch_particlestomesh.vtp',
 'n_vertices': 21268,
 'n_faces': 34276,
 '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.39800089762767943,
   'max': 0.11979361095028015,
   'mean': -0.018378446349099103},
  'gaussian_curvature': {'min': -0.04719296892544032,
   'max': 0.08088003119086297,
   'mean': -0.0005263805250443599},
  'k1': {'min': -0.1435742579183122,
   'max': 0.21110538547552327,
   'mean': 0.013931785486257675},
  'k2': {'min': -0.7760180805330399,
   'max': 0.06802756509271167,
   'mean': -0.05068867818445588},
  'curvature_anisotropy': {'min': 3.29586417993355e-05,
   'max': 0.7606621993125067,
   'mean': 0.06462046367071356},
  'shape_index': {'min': -0.9995736627961141,
   'max': 0.996754955995085,
   'mean': -0.2461664514640657},
  'curvedness': {'min': 0.0006154566632768757,
   'max': 0.5489095600225444,
   'mean': 0.04427894912544514}}}

Filter particles: ray hits

particle_id in the hits DataFrame is the 0-based row index into the motive list. This gives the subset of particles that actually intersected the mesh, used for the patch visualizations below and saved for downstream use.

# particle_id is the 0-based row index into atp_synthase_motl
hit_ids_particle_mesh = particle_to_mesh_hits["particle_id"].values

# Subset the motive list and coordinates to only intersecting particles
atp_synthase_motl_hits_particle_mesh = cryomotl.Motl(
    atp_synthase_motl.df.iloc[hit_ids_particle_mesh].reset_index(drop=True)
)
atp_synthase_xyz_hits_particle_mesh = atp_synthase_xyz[hit_ids_particle_mesh]

print(f"Intersecting particles: {len(hit_ids_particle_mesh)} / {len(atp_synthase_motl)}")

atp_synthase_motl_hits_particle_mesh.write_out(output_path + "atp_hits_particlestomesh.em")

Intersecting particles: 360 / 498

Visualize patches: ray hits

Each patch is shown colored by Gaussian curvature with the intersecting ATP synthase positions overlaid.

fig = vp.plot_vtp_mesh(
    output_path + "IMM_below10nm_patch_particlestomesh.vtp",
    color_by="gaussian_curvature",
    colorscale_range=(-0.003, 0.003),
    mesh_labels="r ≤ 10 nm",
    overlay_points=atp_synthase_xyz_hits_particle_mesh,
    overlay_names=["ATP synthase (intersecting)"],
    overlay_colors=["gray"],
    overlay_marker_size=3,
    title="Intersecting particles → mesh: Gaussian curvature (r ≤ 10 nm)",
)
#fig.show()
fig.write_image("plots/IMM_rbelow10nm_Gaussiancurvature+ATPsynthases.png", width=1200, height=800)

Image("plots/IMM_rbelow10nm_Gaussiancurvature+ATPsynthases.png")

fig = vp.plot_vtp_mesh(
    output_path + "IMM_10nmto20nm_patch_particlestomesh.vtp",
    color_by="gaussian_curvature",
    colorscale_range=(-0.003, 0.003),
    mesh_labels="10 < r ≤ 20 nm",
    overlay_points=atp_synthase_xyz_hits_particle_mesh,
    overlay_names=["ATP synthase (intersecting)"],
    overlay_colors=["gray"],
    overlay_marker_size=3,
    title="Intersecting particles → mesh: Gaussian curvature (10 < r ≤ 20 nm)",
)
#fig.show()
fig.write_image("plots/IMM_10nmto20nm_Gaussiancurvature+ATPsynthases.png", width=1200, height=800)

Image("plots/IMM_10nmto20nm_Gaussiancurvature+ATPsynthases.png")

Mesh to particles

The reverse path asks: which mesh vertices are close to the ATP synthase cloud when casting from the mesh toward particles? This approach is useful when particle orientations are unknown.

1. Build an ``OrientedPointCloud`` from particle coordinates (normals optional).

2. Wrap it in ``PleomorphicSurface``.

3. Call ``distance_to_pointcloud`` on the mesh with the particle cloud as target.

For a screened-Poisson mesh, vertex normals are not always consistently oriented. ``bidirectional=True`` casts along ±normal and keeps the closer hit. ``one_hit_per_target=True`` keeps one mesh vertex per particle (closest approach), analogous to one triangle per particle in the ray path.

Parameter	Value here	Meaning
max_distance	20.0	Same distance cutoff as the ray path (nm)
bidirectional	True	Try both normal directions
one_hit_per_target	True	One source vertex per target particle

Build particle cloud

# Construct an oriented point cloud from the ATP synthase positions
atp_synthase_pc = surface.OrientedPointCloud()

# Write vertices
atp_synthase_pc.vertices = atp_synthase_xyz

# Write normals (optional, the code runs just the same for unoriented data in this case)
atp_synthase_pc.normals = atp_synthase_normals  # optional but useful

# Convert to PleomorphicSurface instance
atp_synthase_surface = assembly.PleomorphicSurface(atp_synthase_pc)

# Intersect the ATP synthase point cloud starting from the mesh
# All mesh vertices are used
# If one_hit_per_target, return only the vertex that has the closest distance to the intersected points

intersect_mesh_with_pc = IMM.distance_to_pointcloud(
    target=atp_synthase_surface,
    max_distance=20.0,
    bidirectional=True,
    one_hit_per_target=True
)

Ray casting: 382/133654 mesh vertices intersect with point cloud

Grow regions around vertex hits

Same helper as above, with ``query_type=”distance_to_pointcloud”``.

Parameter	Value here	Meaning
source_id_name	"vertex_id"	Mesh vertex that intersects the cloud
target_id_name	"particle_id"	Closest particle index
surface_radii	[10.0, 20.0]	Same inner / outer radii as particles→mesh
surface_element	default	Seeds are intersecting vertices; regions expand to touching triangles for curvature summaries

Counts and curvature means may differ slightly from the particles→mesh table (different seeding: triangles vs vertices), but trends in inner vs outer bands should be comparable when distance cutoffs and radii match.

# Compute intersection data
# Regions are grown from the intersecting vertices (important to be one hit per target)

mesh_to_particles_data = IMM.intersection_data(
    result=intersect_mesh_with_pc,
    query_type="distance_to_pointcloud",
    max_distance_source_target=20.0,
    source_id_name="vertex_id",
    target_id_name="particle_id",
    surface_radii=[10.0, 20.0],
    surface_element="vertices",  # optional: default for mesh distance queries
    include_curvatures=True,
)

mesh_to_particle_hits = mesh_to_particles_data["hits"]
mesh_to_particles_data["region_summary"]

	region	n_triangles	mean_curvature_mean	mean_curvature_median	gaussian_curvature_mean	gaussian_curvature_median
0	hit triangles	2229	-0.020782	-0.033339	0.000760	0.000725
1	r <= 10 nm	79749	-0.020765	-0.029974	0.000453	0.000602
2	r <= 20 nm	115543	-0.019268	-0.025064	0.000198	0.000336
3	10 < r <= 20 nm	35794	-0.015931	-0.014321	-0.000371	-0.000034

Extract and save membrane patches: vertex hits

Extract annular patches the same way as in the particles→mesh section, indexing regions from the mesh→particles DataFrame.

# Extract sub-mesh region
# Region within an ATP synthase particle

IMM_inner_patch_meshtopc = IMM.extract_region(
    mesh_to_particles_data["regions"]["r <= 10 nm"],
    element="triangles",
)

# Surrounding annulus
IMM_outer_patch_meshtopc = IMM.extract_region(
    mesh_to_particles_data["regions"]["10 < r <= 20 nm"],
    element="triangles",
)

# Save the extracted meshes

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

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

Saved mesh with curvatures to outputs/IMM_below10nm_patch_meshtopc.vtp
  Format: VTP (VTK PolyData)
  Vertices: 43,226
  Faces: 79,749
  Scalar fields (8):
    - mean_curvature: [-4.362034e-01, 1.197936e-01]
    - gaussian_curvature: [-4.039314e-02, 1.177193e-01]
    - k1: [-1.668447e-01, 2.111054e-01]
    - k2: [-7.482412e-01, 5.034715e-02]
    - curvature_anisotropy: [1.156415e-04, 7.015738e-01]
    - shape_index: [-9.928864e-01, 9.974568e-01]
    - curvedness: [6.154567e-04, 5.301145e-01]
    - shape_category: [0, 8] (int codes)
  Vector fields (3): normals, principal_direction_1, principal_direction_2
  Curvature units: nm^(-1) (principal/mean), nm^(-2) (Gaussian)
Saved mesh with curvatures to outputs/IMM_10nmto20nm_patch_meshtopc.vtp
  Format: VTP (VTK PolyData)
  Vertices: 22,030
  Faces: 35,794
  Scalar fields (8):
    - mean_curvature: [-3.980009e-01, 1.186291e-01]
    - gaussian_curvature: [-4.719297e-02, 9.442748e-02]
    - k1: [-1.762917e-01, 2.889283e-01]
    - k2: [-7.760181e-01, 6.802757e-02]
    - curvature_anisotropy: [7.166662e-05, 7.606622e-01]
    - shape_index: [-9.925852e-01, 9.974568e-01]
    - curvedness: [6.154567e-04, 5.489096e-01]
    - shape_category: [0, 8] (int codes)
  Vector fields (3): normals, principal_direction_1, principal_direction_2
  Curvature units: nm^(-1) (principal/mean), nm^(-2) (Gaussian)

{'format': 'VTP',
 'path': 'outputs/IMM_10nmto20nm_patch_meshtopc.vtp',
 'n_vertices': 22030,
 'n_faces': 35794,
 '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.39800089762767943,
   'max': 0.11862914966857435,
   'mean': -0.01657110790131637},
  'gaussian_curvature': {'min': -0.04719296892544032,
   'max': 0.09442748006071483,
   'mean': -0.00036352574528635803},
  'k1': {'min': -0.17629172841642776,
   'max': 0.2889282678702083,
   'mean': 0.013808670975453581},
  'k2': {'min': -0.7760180805330399,
   'max': 0.06802756509271167,
   'mean': -0.046950886778086315},
  'curvature_anisotropy': {'min': 7.166661551948582e-05,
   'max': 0.7606621993125067,
   'mean': 0.0607595577535399},
  'shape_index': {'min': -0.9925852342577581,
   'max': 0.9974567930927871,
   'mean': -0.20062267591277785},
  'curvedness': {'min': 0.0006154566632768757,
   'max': 0.5489095600225444,
   'mean': 0.04257784849809986}}}

Filter particles: vertex hits

Same filtering as in the particles→mesh section, using the mesh→particles hits DataFrame.

# particle_id is the 0-based row index into atp_synthase_motl
hit_ids_mesh_particle = mesh_to_particle_hits["particle_id"].values

# Subset the motive list and coordinates to only intersecting particles
atp_synthase_motl_hits_mesh_particle = cryomotl.Motl(
    atp_synthase_motl.df.iloc[hit_ids_mesh_particle].reset_index(drop=True)
)
atp_synthase_xyz_hits_mesh_particle = atp_synthase_xyz[hit_ids_mesh_particle]

print(f"Intersecting particles: {len(hit_ids_mesh_particle)} / {len(atp_synthase_motl)}")

atp_synthase_motl_hits_mesh_particle.write_out(output_path + "atp_hits_meshtoparticles.em")

Intersecting particles: 382 / 498

Visualize patches: vertex hits

Same visualization as the particles→mesh patches. Using the same colorscale_range across both directions makes the curvature distributions directly comparable.

fig = vp.plot_vtp_mesh(
    output_path + "IMM_below10nm_patch_meshtopc.vtp",
    color_by="gaussian_curvature",
    colorscale_range=(-0.003, 0.003),
    mesh_labels="r ≤ 10 nm",
    overlay_points=atp_synthase_xyz_hits_mesh_particle,
    overlay_names=["ATP synthase (intersecting)"],
    overlay_colors=["gray"],
    overlay_marker_size=3,
    title="Mesh → particles: Gaussian curvature r ≤ 10 nm",
)
#fig.show()
fig.write_image("plots/IMM_below10nm_mesh-to-particles_Gaussiancurvature+ATPsynthases.png", width=1200, height=800)

Image("plots/IMM_below10nm_mesh-to-particles_Gaussiancurvature+ATPsynthases.png")

fig = vp.plot_vtp_mesh(
    output_path + "IMM_9nmto18nm_patch_meshtopc.vtp",
    color_by="gaussian_curvature",
    colorscale_range=(-0.003, 0.003),
    mesh_labels="10 < r ≤ 20 nm",
    overlay_points=atp_synthase_xyz_hits_mesh_particle,
    overlay_names=["ATP synthase (intersecting)"],
    overlay_colors=["gray"],
    overlay_marker_size=3,
    title="Mesh → particles: Gaussian curvature 10 < r ≤ 20 nm",
)
#fig.show()
fig.write_image("plots/IMM_10nmto20nm_mesh-to-particles-Gaussiancurvature+ATPsynthases.png", width=1200, height=800)

Image("plots/IMM_10nmto20nm_mesh-to-particles-Gaussiancurvature+ATPsynthases.png")

Per-patch curvature distributions

Beyond coloring the mesh, we can compare the distributions of curvature descriptors between the two grown patches — the region under each ATP synthase (r ≤ 10 nm) and the surrounding annulus (10 < r ≤ 20 nm). We read the per-vertex fields stored in the saved patch .vtp files and histogram them natively with plot_histogram from visplot.

We plot one figure per descriptor, overlaying the two patches:

Field	Meaning	Range
mean_curvature	(k1 + k2) / 2	signed
gaussian_curvature	k1 * k2	signed
shape_index	local shape class (cup … cap)	[-1, 1]

Histograms are normalized to a probability density so the two patches are comparable despite having different vertex counts.

# In-memory patches from the particle → mesh ray hits — curvatures preserved
patches = {
    "r ≤ 10 nm":      IMM_inner_patch,
    "10 < r ≤ 20 nm": IMM_outer_patch,
}

# Map each field name to the PleomorphicSurface accessor that returns it
get_field = {
    "mean_curvature":     lambda s: s.get_mean_curvature(),
    "gaussian_curvature": lambda s: s.get_gaussian_curvature(),
    "shape_index":        lambda s: s.get_shape_index(),
}

field_label = {
    "mean_curvature":     "Mean curvature (nm⁻¹)",
    "gaussian_curvature": "Gaussian curvature (nm⁻²)",
    "shape_index":        "Shape index",
}

def field_frame(field):
    """Wide DataFrame: one NaN-padded column of `field` values per patch."""
    cols = {name: np.asarray(get_field[field](s)) for name, s in patches.items()}
    n = max(map(len, cols.values()))
    return pd.DataFrame(
        {k: np.concatenate([v, np.full(n - len(v), np.nan)]) for k, v in cols.items()}
    )

for field in ["mean_curvature", "gaussian_curvature", "shape_index"]:
    fig = vp.plot_histogram(
        field_frame(field),
        bins=60,
        hist_norm="probability density",
        colors="Klimt",
        opacity=0.7,
    )
    fig.update_layout(
        title=f"IMM patches – {field_label[field]}",
        xaxis_title=field_label[field],
        barmode="overlay",
    )
    #fig.show()
    fig.write_image(f"plots/IMM_patches_hist_{field}.png", width=900, height=600)

Image(f"plots/IMM_patches_hist_mean_curvature.png")

Image(f"plots/IMM_patches_hist_gaussian_curvature.png")

Image(f"plots/IMM_patches_hist_shape_index.png")