Contour Tools

Contour extraction and processing utilities.

Module Reference

Tools for creating and manipulating contours.

class physiotwin4d.contour_tools.ContourTools(log_level=20)

Tools for creating and manipulating contours.

__init__(log_level=20)

Initialize ContourTools.

Parameters:

log_level (int | str) – Logging level (default: logging.INFO)

apply_anatomy_color(mesh, anatomy_names)

Attach a structure’s USDAnatomyTools color in-place.

Sets, as WorkflowConvertImageToVTK._annotate() does, so geometry from here colors the same way in Paraview, PyVista, and the USD exporter:

• field_data['AnatomyColor'] — RGB float32 color.

• cell_data['Color'] — RGBA uint8 solid color (n_cells × 4).

Parameters:

• mesh (DataSet) – Surface or volume mesh to annotate.

• anatomy_names (Sequence[str]) – Names tried in order, most specific first, e.g. an organ name followed by its anatomy group. USDAnatomyTools carries overrides for some organs (myocardium) but not others (left_ventricle), so a group name is the usual second entry. Falls back to 'other' when none resolves.

Return type:

None

extract_contours(labelmap_image, smoothing_iterations=10, smoothing_scale=1.0, surface_reduction_rate=0.0, taubin_iterations=20)

Make contours from a labelmap image.

Every label boundary is emitted, including the internal boundaries between adjacent labels, so the result is a multi-material surface and is not watertight: an edge where three labels meet is shared by three faces. Use extract_watertight_surface() when a single label’s closed surface is needed.

Two passes keep the result off the voxel block edges the labelmap is drawn on. An anisotropic labelmap is first resampled onto an isotropic grid of its finest pitch, the way extract_label_surfaces() does, so that a boundary between two thick slices lands between them instead of terracing at one of them. The contour is then Taubin-smoothed, which the surface net’s own constrained smoothing cannot substitute for: that one may not move a point more than about a voxel, so it rounds the blocks without removing them. Taubin does not shrink the surface, and the surface net shares its points between neighboring labels, so smoothing the mesh as one moves a shared point once and the labels stay in contact.

Parameters:

• labelmap_image (itk.image) – The labelmap image to create contours from

• smoothing_iterations (int) – Surface-net smoothing iterations.

• smoothing_scale (float) – Surface-net smoothing scale.

• surface_reduction_rate (float) – Fraction of triangles to remove afterwards (0.0 disables).

• taubin_iterations (int) – Taubin smoothing iterations applied after reduction, so they act on evenly sized triangles (0 disables).

Returns:

The contours as a PyVista PolyData object

Return type:

pv.PolyData

extract_watertight_surface(mask_image, smoothing_iterations=10, gaussian_sigma_mm=0.5, surface_reduction_rate=0.0, anatomy_names=None)

Extract one binary mask’s closed, outward-oriented surface.

extract_contours() cannot produce a watertight surface: its surface nets pinch where a mask self-touches across a voxel diagonal, leaving edges shared by four faces, and they leave the mask open where it reaches the image border. Isocontouring a continuous field cannot pinch, so this pads the mask with one voxel of background, blurs it, and runs marching cubes at the half-way isovalue instead.

Reduction goes through remesh_and_smooth_surface(), whose ACVD remeshing keeps a watertight input watertight where the VTK decimators do not. That is a property of the remesher rather than a guarantee, so a reduced result is still checked and a warning logged if it degrades.

Parameters:

• mask_image (image) – Binary mask holding the single structure to extract.

• smoothing_iterations (int) – Taubin smoothing iterations (0 disables).

• gaussian_sigma_mm (float) – Blur applied before isocontouring, in millimeters, so it is independent of the voxel pitch.

• surface_reduction_rate (float) – Fraction of triangles to remove (0.0 disables).

• anatomy_names (Optional[Sequence[str]]) – Names passed to apply_anatomy_color(), most specific first. None leaves the surface uncolored.

Return type:

PolyData

Returns:

The structure’s surface, with outward normals.

extract_label_surfaces(labelmap_image, isotropic_spacing_mm=None, distance_sigma_mm=None, smoothing_iterations=30)

Extract every label’s surface, smooth and conforming with its neighbors.

extract_watertight_surface(), run per label, traces each label’s own voxel block edges: on anisotropic data the result terraces at the slice pitch, and neighboring labels are contoured independently, so their shared wall is meshed twice and the two copies do not match. This extracts all labels together instead:

1. The labelmap is resampled onto an isotropic grid with ITK’s label-aware Gaussian interpolator, which votes over the labels in a physical-space kernel rather than picking a nearest voxel. That is what removes the terracing: a boundary between two slices lands between them instead of on one of them.

2. Every label, plus the background, gets a signed distance map on that grid, and each voxel is assigned to the label whose map is smallest. The assignment is a partition, so no gap or overlap can arise.

3. Label L’s surface is the zero level of D_L minus the smallest of the other maps. On a wall between L and M that field is the negation of M’s, so marching cubes puts identical vertices on both surfaces and the two meet exactly.

4. The surfaces are merged, which welds those coincident vertices, then Taubin-smoothed as one mesh. Smoothing therefore moves a shared vertex once and the surfaces stay in contact.

Structures thinner than the interpolation kernel lose volume, coronary arteries most of all; the fraction of each label’s voxel volume that the surface encloses is logged.

Parameters:

• labelmap_image (image) – Multi-label image; every non-zero label present is extracted. A binary mask yields a single surface.

• isotropic_spacing_mm (Optional[float]) – Edge length of the isotropic grid the surfaces are contoured on, which sets both their smoothness and their triangle count. None uses the labelmap’s finest spacing.

• distance_sigma_mm (Optional[float]) – Blur applied to the distance maps, which is what takes the voxel facets out of the contoured surface. None uses the isotropic spacing; raising it smooths further and thins the smallest structures.

• smoothing_iterations (int) – Taubin smoothing iterations (0 disables).

Return type:

dict[int, PolyData]

Returns:

Label id → that label’s closed, outward-oriented surface. A label too small to survive the isotropic grid is left out, so the mapping is empty when the labelmap holds no non-zero label and may be missing labels that it does hold.

static is_watertight(surface)

Report whether every edge of surface is shared by exactly two faces.

A surface with no faces has no edge that fails the test, so it is reported as not watertight rather than vacuously watertight.

Return type:

bool

extract_tetrahedra(mask_image, element_size_mm=None, anatomy_names=None)

Build a tetrahedral mesh filling one binary mask.

Every retained voxel becomes a hexahedron, which VTK then splits into six tetrahedra sharing the hexahedra’s points, so the result is a conforming mesh whose boundary is the voxel staircase rather than the smooth surface extract_label_surfaces() returns. Pass the result through trim_tetrahedra_to_surface() to relax that staircase onto the surface.

Parameters:

• mask_image (image) – Binary mask holding the single structure to fill.

• element_size_mm (Optional[float]) – Edge length of the isotropic voxels the mask is resampled to before meshing, which is the resulting element size. None meshes the mask’s own voxels, so on anisotropic data the elements inherit that anisotropy. A size above the thinnest part of the structure drops that part.

• anatomy_names (Optional[Sequence[str]]) – Names passed to apply_anatomy_color(), most specific first. None leaves the mesh uncolored.

Return type:

UnstructuredGrid

Returns:

The structure’s tetrahedral mesh, empty if the mask is empty or element_size_mm is too coarse to keep any of it.

trim_tetrahedra_to_surface(tetrahedra, surface, iterations=5, relaxation=0.6, min_scaled_jacobian=0.1)

Relax a tetrahedral mesh onto surface, keeping every cell whole.

extract_tetrahedra() meshes voxels, so its boundary is a staircase that both protrudes through the smooth surface extract_label_surfaces() builds from the same mask and falls short of it elsewhere. Cutting the mesh at the surface with clip_surface would follow it exactly but shatters the boundary tetrahedra into slivers (about a tenth of the cells drop below a scaled Jacobian of 0.1), and neither VTK nor any current dependency can repair those.

So nothing is cut. A crinkle clip drops the cells that lie entirely outside while leaving every surviving cell intact, then the mesh is relaxed: each sweep moves every point part of the way toward the average of its neighbors, and every boundary point instead toward its closest point on surface. The interior smoothing is what makes room for the boundary to reach the surface – projecting the boundary alone flattens the cells behind it, and the quality bound below then undoes the move, which is why one projection pass leaves the staircase in place.

A move that would still wreck a cell is backed off: every point of a cell below min_scaled_jacobian has that sweep’s step halved, repeatedly, until the whole mesh clears the bound.

Parameters:

• tetrahedra (UnstructuredGrid) – Volume mesh to relax; its cell and field data survive.

• surface (PolyData) – Closed surface to relax onto, in the same frame.

• iterations (int) – Relaxation sweeps. The boundary reaches the surface in the first few and then stops moving: on the Duke heart labels, sweeps beyond the default leave the mean boundary-to-surface distance and the worst cell quality where they already were, and only cost time.

• relaxation (float) – Fraction of the way to its target a point moves per sweep. 1.0 moves the whole way and oscillates.

• min_scaled_jacobian (float) – Cell-quality bound every tetrahedron must meet after each sweep. 0.0 only rules out flattened and inverted cells; the default also rules out slivers.

Return type:

UnstructuredGrid

Returns:

The relaxed mesh.

remesh_and_smooth_surface(surface, surface_reduction_rate=0.0, smoothing_iterations=0)

Optionally remesh then smooth a surface (no-op when disabled).

Reduction is isotropic remeshing (ACVD, via pyacvd) rather than decimation: the surface is re-tiled with uniform, well-shaped triangles at the requested resolution. decimate_pro reaches the same triangle count but leaves a watertight input non-watertight; ACVD does not.

Remeshing rebuilds the topology and so discards cell data, exactly as decimate_pro did: per-cell boundary_labels (needed for anatomy splitting downstream) are transferred back onto the new cells from their nearest original cell so anatomy materials still apply. Uniform triangles cannot represent a label patch smaller than one of them, though, so such a patch is absorbed by its neighbours and its label pair disappears – a warning names the pairs lost. decimate_pro kept those patches by being non-uniform, which is the trade being made here. Smoothing uses non-shrinking Taubin smoothing, which only moves points and therefore preserves cells and their labels. It is told to move non-manifold points too, since on a multi-material surface every edge where three labels meet is non-manifold and VTK pins those points otherwise; on a manifold surface the setting has nothing to act on.

Parameters:

• surface (PolyData) – Input surface.

• surface_reduction_rate (float) – Fraction of triangles to remove (0.0 disables).

• smoothing_iterations (int) – Taubin smoothing iterations (0 disables).

Return type:

PolyData

Returns:

The remeshed and/or smoothed surface.

static extract_surface(mesh)

Extract the surface of a mesh.

Parameters:

mesh (DataSet) – Input mesh (PolyData is returned unchanged; any other DataSet is passed through extract_surface).

Returns:

The surface of the mesh.

Return type:

pv.PolyData

static transform_contours(contours, tfm, with_deformation_magnitude=False)

Transform contours using a given transform.

Parameters:

tfm (itk.Transform) – The transform to use

Returns:

The transformed contours with deformation magnitude

Return type:

pv.PolyData

merge_meshes(meshes)

Merge multiple fixed meshes into a single mesh.

Returns:

Merged mesh

Return type:

pv.PolyData

static create_reference_image(mesh, spatial_resolution=0.5, buffer_factor=0.25, ptype=itk.F)

Create a reference image from a mesh.

Return type:

Image

static create_mask_from_mesh(mesh, reference_image)

Return type:

itk.Image

create_labelmap_from_meshes(meshes, reference_image)

Create a labelmap from a list of meshes.

Return type:

itk.Image

static sample_mesh_faces(mesh, max_spacing)

Return mesh points supplemented by samples across the triangle faces.

Rasterizing vertices alone leaves gaps between them on meshes that are coarse relative to the voxel size, which makes a distance map built from them ripple. Adding barycentric samples dense enough that consecutive samples are closer than max_spacing closes those gaps.

Parameters:

• mesh (DataSet) – Source mesh; its surface is triangulated if needed.

• max_spacing (float) – Target spacing between samples, in mm.

Return type:

ndarray

Returns:

(n, 3) array of sample points, starting with the mesh’s own points.

create_distance_map(mesh, reference_image, squared_distance=False, negative_inside=True, zero_inside=False, norm_to_max_distance=0.0, sample_faces=True)

Compute a distance map of a mesh on the reference image’s grid.

Parameters:

• mesh (pv.DataSet | pv.UnstructuredGrid) – Mesh whose surface the distances are measured to.

• reference_image (itk.Image) – Image defining the output grid.

• squared_distance (bool) – Sign-preserving square of the result. Default: False

• negative_inside (bool) – Keep the signed output. Default: True

• zero_inside (bool) – Clip negative values to zero before anything else. Default: False

• norm_to_max_distance (float) – If non-zero, divide by this value and clip to [-1, 1]. Default: 0.0 (distances stay in mm)

• sample_faces (bool) – Rasterize samples across the triangle faces as well as the vertices, so that coarse meshes do not leave gaps in the rasterized surface. Default: True

Return type:

itk.Image

Returns:

ITK image of distances on the reference grid.

static create_deformation_field(points, point_displacements, reference_image, blur_sigma=2.5, ptype=itk.D)

Create a displacement map from model points and displacements.

Return type:

Image

static save_surfaces(surfaces, output_dir, prefix='')

Save each named surface to its own VTP file.

Parameters:

• surfaces (dict[str, PolyData]) – Mapping of name → surface (e.g. the 'surfaces' value from WorkflowConvertImageToVTK.process()).

• output_dir (str) – Directory to write files into (created if absent).

• prefix (str) – Optional filename prefix. Each file is named {prefix}_{name}.vtp (or {name}.vtp when prefix is empty).

Return type:

dict[str, str]

Returns:

Mapping of name → absolute path of the saved file.

static save_combined_surfaces(surfaces, output_filename)

Merge all named surfaces into a single VTP file.

The merged mesh retains per-cell Color (RGBA uint8) from each surface’s annotation, enabling colour-by-anatomy rendering in Paraview, PyVista, etc.

It also gains a per-cell SegmentationLabelIds (int32) array, which carries each cell’s originating label ID so structure identity survives the merge. Downstream, ConvertVTKToUSD splits on this array when given mask_ids, giving one prim (and one anatomy material) per structure. A surface whose field_data['SegmentationLabelIds'] does not hold exactly one ID has no per-cell attribution — that is the case for the per-group surfaces of WorkflowConvertImageToVTK, which are contoured from a merged binary mask — so its cells are tagged 0. Pass the per-label surfaces ('label_surfaces') to get real IDs.

Per-object field_data is not preserved: it is per-object, so a single merged mesh cannot carry one value per input surface. The remaining keys set by WorkflowConvertImageToVTK._annotate() are therefore lost:

• AnatomyGroup — group name, e.g. 'heart'.

• SegmentationLabelNames — structure names within the group.

• AnatomyColor — RGB float color (survives indirectly as the per-cell Color array).

Use save_surfaces() instead when structure names must be recoverable from the saved files.

Parameters:

• surfaces (dict[str, PolyData]) – Mapping of name → surface.

• output_filename (str) – Path of the VTP file to write, including its directory. Any missing parent directories are created.

Return type:

str

Returns:

Path to the saved VTP file.

Raises:

ValueError – If surfaces is empty.

Navigation

Transform Tools | Utility Modules | 4D Image Conversion