Sampling capsid surfaces#

This tutorial loads a segmented capsid volume and samples uniformly distributed oriented points on its inner and outer surfaces.

The workflow uses PleomorphicSurface as a wrapper. When methods are called on the wrapper (read, refine_normals, oversample, separate_surfaces, …), it forwards to the right surface implementation and, when an operation returns new geometry, wraps the result again so you stay in the same API. Under the hood, the surface is represented as one of two types:

  • Mesh: Triangle surface with vertices, faces, and normals.

  • OrientedPointCloud: Discrete points with normals.

Meshes and oriented point clouds share many methods (normal refinement, cropping, inner/outer separation, distance queries) but these can differ in their implementation. The representation can also change with certain methods. PleomorphicSurface internally tracks these changes so you do not re-wrap after every step.

The wrapper exposes ``is_mesh`` and ``is_point_cloud`` so you can check which representation you have at any point.

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

Table of contents#

  1. Setup

  2. Load segmentation

  3. Separate inner and outer surfaces

  4. Oversample mesh surfaces

  5. Oversample point clouds

  6. Save results

Workflow#

  1. Load the segmentation as a mesh-backed PleomorphicSurface.

  2. Separate inner and outer capsid surfaces by normal orientation around the centroid. (Optional)

  3. Oversample each surface to an approximately uniform oriented point cloud.

  4. Save meshes, motive lists, and point cloud objects.

Setup#

Set the path to the input segmentation MRC file. All outputs will be written to output_path.

[1]:

%load_ext autoreload
%autoreload 2

[2]:

import os
import numpy as np
import copy
from cryocat.analysis import assembly
import cryocat.analysis.visplot as vp
from IPython.display import Image

Load segmentation#

The capsid surface is built with ``PleomorphicSurface.read`` using method="mesh_from_mrc". Under the hood this runs marching cubes on the MRC volume and stores a ``Mesh`` in capsid.surface.

Marching cubes produces vertices, faces, and normals. Vertex coordinates are scaled by ``pixel_size`` during extraction:

  • If pixel_size=1, coordinates are in pixel/voxel units.

  • With a known physical spacing, e.g. pixel_size=0.5 (nm) or pixel_size=5.0 (Å) → vertices in those units.

The unit label is stored separately with .units and can be set after creating the mesh:

capsid.units = "pixel"       # or "nm", "A"

The main other choices at this stage are whether to smooth the segmentation before meshing and whether to call refine_normals() after loading.

Parameter

Typical value

Meaning

input_path

MRC file path

Segmentation volume to load.

transpose

True

Uses the package convention for reading MRC volumes. Keep this consistent across the notebook.

pixel_size

1 if no spacing is known

Coordinate spacing. If 1, output coordinates are in pixel/voxel units. If physical pixel size is known, use that value instead.

smooth_sigma

None to 2.5

Gaussian smoothing applied before marching cubes, in pixel/voxel units. Larger values reduce blocky artifacts but can also soften details.

level

0.5

Iso-value for marching cubes. For binary masks, 0.5 is usually appropriate.

 
[3]:

# Input segmentation path
input_path = "inputs/"
capsid_seg = input_path + "capsid_mask.mrc"

[ ]:

capsid = assembly.PleomorphicSurface.read(
    input_path=capsid_seg,
    method="mesh_from_mrc",
    transpose=True,
    pixel_size=1,
    smooth_sigma=2.0,
    level=0.5,
)

# Optional: after marching cubes; tune via radius_hit / n_iter on the method
capsid = capsid.refine_normals(radius_hit=3.0, n_iter=1)

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

100%|██████████| 243/243 [00:13<00:00, 18.50it/s]

Normal refinement complete.



[5]:

capsid.units = "pixel"

[6]:

# Check internal representation
print(f"capsid.is_mesh: {capsid.is_mesh}")
print(f"capsid.is_point_cloud: {capsid.is_point_cloud}")

capsid.is_mesh: True
capsid.is_point_cloud: False

Separate inner and outer surfaces#

Optional. Skip this section if the segmentation contains only a single surface.

For a “hollow” capsid segmentation, the outer surface normals point away from the capsid center and the inner surface normals point toward it. During separation, each normal is compared to the direction from the vertex toward a reference point (usually the centroid). Vertices whose normals point toward the reference are assigned to the inner surface; the rest to the outer surface.

separate_surfaces dispatches to the appropriate separation algorithm based on surface_type and returns two boolean masks without modifying the geometry.

Parameter

Default

Meaning

surface_type

'closed'

Separation strategy. 'closed' for enclosed volumes (capsids, vesicles); 'planar' for surfaces with lower curvature.

threshold_angle

90.0

Angle cutoff in degrees. Smaller values make the inner assignment stricter. Only used when surface_type='closed'.

reference_point

None

Reference point for classification. If None, the centroid is used. For asymmetric objects, provide a point known to lie inside the capsid. Only used when surface_type='closed'.

This approach works best for closed or nearly closed nested surfaces where normals are consistently oriented. If the normals are noisy or the shape is strongly open or asymmetric, inspect the result before continuing.

[7]:

inner_mask, outer_mask = capsid.separate_surfaces(surface_type='closed')

Extract the inner and outer surfaces#

apply_vertex_mask filters the surface to only the vertices where the mask is True, without modifying the original object.

Mask

Meaning

inner_mask

Vertices whose normals point toward the centroid (inner capsid shell).

outer_mask

Vertices whose normals point away from the centroid (outer capsid shell).

The result keeps the same representation as the input: if the input is a mesh, the output is a mesh containing only the selected vertices and valid faces. For this workflow, inner_surface and outer_surface are mesh-backed PleomorphicSurface objects after filtering.

[8]:

inner_surface = capsid.apply_vertex_mask(inner_mask)
outer_surface = capsid.apply_vertex_mask(outer_mask)

[9]:

print(f"inner_surface.is_mesh: {inner_surface.is_mesh}")
print(f"outer_surface.is_mesh: {outer_surface.is_mesh}")
print(f"inner_surface.is_point_cloud: {inner_surface.is_point_cloud}")
print(f"outer_surface.is_point_cloud: {outer_surface.is_point_cloud}")
print(f"Number of points in inner_surface: {len(inner_surface.surface.vertices)}")
print(f"Number of points in outer_surface: {len(outer_surface.surface.vertices)}")

inner_surface.is_mesh: True
outer_surface.is_mesh: True
inner_surface.is_point_cloud: False
outer_surface.is_point_cloud: False
Number of points in inner_surface: 191356
Number of points in outer_surface: 293086

Report mesh surface area#

If the input is a Mesh, get_surface_area() returns the total surface area in the same units as the vertex coordinates. If the input is an OrientedPointCloud, a NotImplementedError is raised — surface area is undefined without triangle connectivity.

[10]:

print(f"inner_surface.get_surface_area(): {inner_surface.get_surface_area():.0f} (pixels^2)")
print(f"outer_surface.get_surface_area(): {outer_surface.get_surface_area():.0f} (pixels^2)")

inner_surface.get_surface_area(): 135340 (pixels^2)
outer_surface.get_surface_area(): 204423 (pixels^2)

Visualize surface separation#

The inner and outer mesh surfaces are plotted as point clouds (vertices only) as a quick sanity check that the surface split is correct. The two surfaces are colored differently so any mixing or assignment errors are immediately visible. For large meshes use sample_fraction to keep the plot responsive.

To visualize the full triangulated mesh, save the .ply files first and use plot_ply_mesh — see the Save results section below.

[11]:

fig = vp.plot_points_with_normals(
    points=[inner_surface.surface.vertices, outer_surface.surface.vertices],
    labels=["Inner surface", "Outer surface"],
    sample_fraction=0.1,
    show_normals=False, # would be too dense
    marker_size=1,
    title="Capsid mesh inner vs outer surface",
)
fig.write_image("plots/capsid_mesh_inner_vs_outer_surface.png", width=1200, height=800)
#fig.show()

[12]:

Image("plots/capsid_mesh_inner_vs_outer_surface.png")

[12]:
../../_images/tutorials_surface_sampling_oversample_capsid_18_0.png

Oversample mesh surfaces#

Oversampling converts each mesh surface into a more uniformly sampled oriented point cloud using Poisson disk sampling. The code first estimates a target number of points from either point_spacing or oversample_factor. If the target is larger than the current number of mesh vertices, the mesh may be subdivided internally before sampling to give the sampler enough geometric detail.

Parameter

Typical use

Meaning

point_spacing

8.0 here

Desired approximate distance between sampled points, in the same units as the mesh vertices.

oversample_factor

optional

Alternative to point_spacing; target point count is set relative to the current vertex count.

poisson_init_factor

5

Controls the candidate pool size for Poisson disk sampling. Larger values can improve uniformity at the cost of runtime.

After mesh oversampling, the wrapped object changes representation from Mesh to OrientedPointCloud, because the output is sampled points with normals rather than connected faces.

[13]:

spacing = 8.0  # same units as mesh vertices (voxels here)

inner_surface_oversampled = inner_surface.oversample(point_spacing=spacing)
outer_surface_oversampled = outer_surface.oversample(point_spacing=spacing)

Target points: 1526 (0.01x original)
Final: 1526 points
Target points: 2302 (0.01x original)
Final: 2302 points

[14]:

print(f"inner_surface_oversampled.is_point_cloud: {inner_surface_oversampled.is_point_cloud}")
print(f"outer_surface_oversampled.is_point_cloud: {outer_surface_oversampled.is_point_cloud}")
print(f"Number of points in inner_surface_oversampled: {len(inner_surface_oversampled.surface.vertices)}")
print(f"Number of points in outer_surface_oversampled: {len(outer_surface_oversampled.surface.vertices)}")

inner_surface_oversampled.is_point_cloud: True
outer_surface_oversampled.is_point_cloud: True
Number of points in inner_surface_oversampled: 1526
Number of points in outer_surface_oversampled: 2302

Visualize oversampled point clouds#

Setting show_normals=True draws each point’s surface normal as a cone. Plotly has no native 3D arrow type, so cones are the closest built-in approximation to the orientation arrows shown in ChimeraX.

[15]:

fig = vp.plot_points_with_normals(
    points=[inner_surface_oversampled.surface.vertices,
            outer_surface_oversampled.surface.vertices],
    normals=[inner_surface_oversampled.surface.normals,
             outer_surface_oversampled.surface.normals],
    labels=["Inner (sampled)", "Outer (sampled)"],
    sample_fraction=0.15,
    show_normals=True,
    normal_scale=0.5,  # same units as spacing (voxels)
    marker_size=1,
    title="Oversampled capsid surfaces with normal orientations",
)
#fig.show()
fig.write_image("plots/capsid_mesh_inner_vs_outer_surface_with_normals.png", width=1200, height=800)

[16]:

Image("plots/capsid_mesh_inner_vs_outer_surface_with_normals.png")

[16]:
../../_images/tutorials_surface_sampling_oversample_capsid_24_0.png

Clean the sampled points#

Two functions let you clean up the oversampled point cloud:

  • ``clean_by_normals(max_angle_deg)`` — removes points whose normal deviates more than max_angle_deg from the surface’s mean normal. Good for trimming sampling outliers/spikes on a locally consistent patch.

  • ``clean_by_angle(max_angle_deg, reference_normal, signed=False)`` — removes points whose normal deviates more than max_angle_deg from a given axis. With signed=True, a normal pointing opposite the axis counts as a 180° deviation (directional selection); with signed=False (default) antiparallel counts as aligned.

[17]:

# Cleaners mutate in place; copy so the original (saved later) is untouched.
inner_deviated = copy.deepcopy(inner_surface_oversampled)
n0 = len(inner_deviated.surface.vertices)
top_cap = copy.deepcopy(inner_surface_oversampled)

# (1) Remove sampling outliers: normals deviating > 45 deg from the mean normal.
inner_deviated.clean_by_normals(max_angle_deg=45)
print(f"clean_by_normals: {n0} -> {len(inner_deviated.surface.vertices)} points")

# (2) Directional selection: keep points whose normals point roughly toward +z.
top_cap.clean_by_angle(max_angle_deg=60, reference_normal=[0, 0, 1], signed=True)
print(f"clean_by_angle (+z, signed): {n0} -> {len(top_cap.surface.vertices)} points")

Cleaned by normals (angle vs mean)
clean_by_normals: 1526 -> 662 points
Cleaned by angle (angle vs given axis)
clean_by_angle (+z, signed): 1526 -> 466 points

[18]:

fig = vp.plot_points_with_normals(
    points=[inner_surface_oversampled.surface.vertices, inner_deviated.surface.vertices],
    normals=[inner_surface_oversampled.surface.normals, inner_deviated.surface.normals],
    labels=["All sampled points", "Points with normals < 45 deg from mean normal"],
    show_normals=False,
    marker_size=2,
    title="clean_by_normals: remove outliers",
)
fig.write_image("plots/capsid_clean_by_normals.png", width=1200, height=800)
#fig.show()

[19]:

Image("plots/capsid_clean_by_normals.png")

[19]:
../../_images/tutorials_surface_sampling_oversample_capsid_28_0.png 
[20]:

fig = vp.plot_points_with_normals(
    points=[inner_surface_oversampled.surface.vertices, top_cap.surface.vertices],
    normals=[inner_surface_oversampled.surface.normals, top_cap.surface.normals],
    labels=["All sampled points", "Kept (normals ~ +z)"],
    show_normals=False,
    marker_size=2,
    title="clean_by_angle: directional selection toward +z",
)
fig.write_image("plots/capsid_clean_by_angle.png", width=1200, height=800)
#fig.show()

[21]:

Image("plots/capsid_clean_by_angle.png")

[21]:
../../_images/tutorials_surface_sampling_oversample_capsid_30_0.png

Oversample point clouds#

This section shows the equivalent workflow when the starting object is already an OrientedPointCloud rather than a mesh — for example, when loading from a motive list or converting directly from the mesh.

A point cloud has no faces and no continuous surface, so mesh-based subdivision and Poisson disk sampling cannot be used. Instead, the implementation applies a KD-tree-based greedy selection on the existing points, keeping only those that are at least point_spacing apart from any already-selected point.

Starting representation

Oversampling method

Output representation

Mesh

Poisson disk sampling (with optional internal subdivision)

OrientedPointCloud

OrientedPointCloud

KD-tree greedy spacing-based selection

OrientedPointCloud

 
[22]:

# Create a point cloud from the capsid mesh
# Or load a point cloud from a motive list or .ply file
capsid_pc = capsid.surface.to_oriented_points()

[ ]:

# Wrap the point cloud as a PleomorphicSurface
capsid_pc = assembly.PleomorphicSurface(capsid_pc)
print(f"capsid_pc.is_point_cloud: {capsid_pc.is_point_cloud}")

capsid_pc.is_point_cloud: True

[24]:

# The same approach for using a centroid to separate the surfaces can be applied to the point cloud
inner_mask_pc, outer_mask_pc = capsid_pc.separate_surfaces(surface_type='closed')

[25]:

# Generate two separate point clouds corresponding to the inner and outer surfaces
inner_surface_pc = capsid_pc.apply_vertex_mask(inner_mask_pc)
outer_surface_pc = capsid_pc.apply_vertex_mask(outer_mask_pc)

[26]:

spacing = 8.0

inner_surface_oversampled_pc = inner_surface_pc.oversample(point_spacing=spacing)
outer_surface_oversampled_pc = outer_surface_pc.oversample(point_spacing=spacing)

Point cloud sampling: 191356 points
  Current spacing: 0.51, target spacing: 8.00
Final: 1411 points
Point cloud sampling: 293086 points
  Current spacing: 0.50, target spacing: 8.00
Final: 2128 points

[27]:

print(f"inner_surface_oversampled_pc.is_point_cloud: {inner_surface_oversampled_pc.is_point_cloud}")
print(f"outer_surface_oversampled_pc.is_point_cloud: {outer_surface_oversampled_pc.is_point_cloud}")
print(f"Number of points in inner_surface_oversampled_pc: {len(inner_surface_oversampled_pc.surface.vertices)}")
print(f"Number of points in outer_surface_oversampled_pc: {len(outer_surface_oversampled_pc.surface.vertices)}")

inner_surface_oversampled_pc.is_point_cloud: True
outer_surface_oversampled_pc.is_point_cloud: True
Number of points in inner_surface_oversampled_pc: 1411
Number of points in outer_surface_oversampled_pc: 2128

Visualize oversampled oriented point clouds#

[28]:

fig = vp.plot_points_with_normals(
    points=[inner_surface_oversampled_pc.surface.vertices,
            outer_surface_oversampled_pc.surface.vertices],
    normals=[inner_surface_oversampled_pc.surface.normals,
             outer_surface_oversampled_pc.surface.normals],
    labels=["Inner point cloud (sampled)", "Outer point cloud (sampled)"],
    show_normals=True,
    sample_fraction=0.15,
    normal_scale=1.0,
    marker_size=1,
    title="Sampled inner and outer point clouds with normals",
)
#fig.show()
fig.write_image("plots/capsid_mesh_inner_vs_outer_surface_with_normals.png", width=800, height=600)

[29]:

Image("plots/capsid_mesh_inner_vs_outer_surface_with_normals.png", width=800, height=600)

[29]:
../../_images/tutorials_surface_sampling_oversample_capsid_40_0.png

Save results#

The .save method is implemented on PleomorphicSurface and dispatches to the correct output format based on the underlying representation.

Object type

Typical output

Method

Mesh

.ply / .vtp mesh file

save(..., format="ply" / "vtp")

OrientedPointCloud

Motive list .em

save(..., tomo_id=..., subtomo_ids=...)

OrientedPointCloud

Point cloud .ply

save(..., format="ply")

In this tutorial, inner_surface and outer_surface are meshes saved as .ply and .vtp files. After oversampling, the resulting OrientedPointCloud objects are exported as motive lists and .ply point clouds.

For motive-list export:

Parameter

Example

Meaning

tomo_id

1

Scalar when all points come from one tomogram, or an array for points from multiple tomograms.

subtomo_ids

np.ones(n_points)

Particle grouping per point. If all points belong to one object, pass an array of ones. Unique labels are mapped to sequential subtomo_id values.

 
[30]:

# Before saving, check the representation
print(f"inner_surface.is_mesh: {inner_surface.is_mesh}")
print(f"outer_surface.is_mesh: {outer_surface.is_mesh}")
print(f"inner_surface_oversampled.is_point_cloud: {inner_surface_oversampled.is_point_cloud}")
print(f"outer_surface_oversampled.is_point_cloud: {outer_surface_oversampled.is_point_cloud}")

inner_surface.is_mesh: True
outer_surface.is_mesh: True
inner_surface_oversampled.is_point_cloud: True
outer_surface_oversampled.is_point_cloud: True

[31]:

output_path = "outputs/"
os.makedirs(output_path, exist_ok=True)

[32]:

# Save the filtered mesh surfaces
inner_surface.save(
    output_path + "capsid_inner_surface.ply",
    format='ply',
    include_curvatures=False)

inner_surface.save(
    output_path + "capsid_inner_surface.vtp",
    format='vtp',
    include_curvatures=False)

outer_surface.save(
    output_path + "capsid_outer_surface.ply",
    format='ply',
    include_curvatures=False)

outer_surface.save(
    output_path + "capsid_outer_surface.vtp",
    format='vtp',
    include_curvatures=False)

Saved mesh to outputs/capsid_inner_surface.ply
  Vertices: 191,356
  Faces: 382,455
Saved mesh to outputs/capsid_inner_surface.vtp
  Format: VTP (VTK PolyData)
  Vertices: 191,356
  Faces: 382,455
Saved mesh to outputs/capsid_outer_surface.ply
  Vertices: 293,086
  Faces: 585,922
Saved mesh to outputs/capsid_outer_surface.vtp
  Format: VTP (VTK PolyData)
  Vertices: 293,086
  Faces: 585,922

[33]:

# Save the oversampled point clouds as motive lists
inner_surface_oversampled.save(
    output_path + "capsid_inner_surface_sampled_points.em",
    tomo_id=1,
    subtomo_ids=np.ones(len(inner_surface_oversampled.surface.vertices))
)

outer_surface_oversampled.save(
    output_path + "capsid_outer_surface_sampled_points.em",
    tomo_id=1,
    subtomo_ids=np.ones(len(outer_surface_oversampled.surface.vertices))
)

Saved 1526 points with normals to outputs/capsid_inner_surface_sampled_points.em
Saved 2302 points with normals to outputs/capsid_outer_surface_sampled_points.em

[34]:

# If the normals should point in the opposite direction, flip them
inner_surface_oversampled_fipped = inner_surface_oversampled.flip_normals()

[35]:

inner_surface_oversampled_fipped.save(
    output_path + "capsid_inner_surface_sampled_points_flipped.em",
    tomo_id=1,
    subtomo_ids=np.ones(len(inner_surface_oversampled.surface.vertices))
)

Saved 1526 points with normals to outputs/capsid_inner_surface_sampled_points_flipped.em

[36]:

# Save the oversampled point clouds as point clouds
inner_surface_oversampled.save(
    output_path + "capsid_inner_surface_sampled_points.ply",
    format='ply',
)

outer_surface_oversampled.save(
    output_path + "capsid_outer_surface_sampled_points.ply",
    format='ply',
)

Saved point cloud with normals to outputs/capsid_inner_surface_sampled_points.ply
Saved point cloud with normals to outputs/capsid_outer_surface_sampled_points.ply

Visualize mesh surfaces#

Now that the .ply files have been saved, render the full triangulated mesh with plot_ply_mesh.

[37]:

fig = vp.plot_ply_mesh(
    ["outputs/capsid_inner_surface.ply", "outputs/capsid_outer_surface.ply"],
    mesh_labels=["Inner surface", "Outer surface"],
    opacity=0.5,
    title="Capsid mesh surfaces",
)
#fig.show()
fig.write_image("plots/capsid_mesh_surfaces.png", width=1200, height=800)

[38]:

Image("plots/capsid_mesh_surfaces.png")

[38]:
../../_images/tutorials_surface_sampling_oversample_capsid_51_0.png

Motive lists#

ChimeraX screenshots of the motive lists

[39]:

print("Sampled inner surface")
Image("media/capsid_inner_surface_motl.png", width=487, height=682)

Sampled inner surface

[39]:
../../_images/tutorials_surface_sampling_oversample_capsid_53_1.png 
[40]:

print("Sampled outer surface")
Image("media/capsid_outer_surface_motl.png", width=487, height=682)

Sampled outer surface

[40]:
../../_images/tutorials_surface_sampling_oversample_capsid_54_1.png 
[41]:

print("Both surfaces")
Image("media/capsid_both_surfaces_motl.png", width=487, height=682)

Both surfaces

[41]:
../../_images/tutorials_surface_sampling_oversample_capsid_55_1.png

Access inner/outer split via object_id#

Instead of re-separating the surfaces every run, you can store the surface assignment as an object_id column in the motive list and reload it already split. Pass the labels through save(..., input_dict={"object_id": ...}), then load with read(method="point_cloud_from_motl", group_by="object_id"), which returns one wrapped PleomorphicSurface per unique object_id.

The group_by argument accepts any motl column, so the same call splits by subtomo_id, object_id, or any other identifier.

[42]:

# Label every point of the full capsid cloud by surface (reusing the earlier masks).
object_ids = np.where(inner_mask_pc, 1, 2)   # 1 = inner surface, 2 = outer surface

combined_motl_path = output_path + "capsid_points_by_object.em"
capsid_pc.save(
    combined_motl_path,
    format="em",
    tomo_id=1,
    subtomo_ids=np.arange(1, len(capsid_pc.surface.vertices) + 1),
    input_dict={"object_id": object_ids},
)

Saved 484442 points with normals to outputs/capsid_points_by_object.em

[ ]:

# Reload already split by object_id -> one PleomorphicSurface per surface
shells = assembly.PleomorphicSurface.read(
    combined_motl_path,
    method="point_cloud_from_motl",
    group_by="object_id",
)
for oid, surf in shells.items():
    print(f"object_id {oid}: {len(surf.surface.vertices)} points, "
          f"is_point_cloud={surf.is_point_cloud}")

object_id 1.0: 191356 points, is_point_cloud=True
object_id 2.0: 293086 points, is_point_cloud=True

[ ]: