Image Processing

The voids.image sub-package provides utilities for segmented image processing, connectivity analysis, and pore network extraction used in vug sensitivity studies.

---

Maximal-Ball Extraction

voids.image.maximal_ball

MaximalBallSettings dataclass

User-facing controls for the native maximal-ball extraction stages.

These settings mirror the main Imperial pnextract controls closely enough for staged verification work, while keeping Python names explicit and readable. The current implementation covers the maximal-ball candidate stage and the initial overlap suppression stage. Hierarchy construction and voxel growth are planned follow-on steps.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class MaximalBallSettings:
    """User-facing controls for the native maximal-ball extraction stages.

    These settings mirror the main Imperial `pnextract` controls closely enough
    for staged verification work, while keeping Python names explicit and
    readable. The current implementation covers the maximal-ball candidate stage
    and the initial overlap suppression stage. Hierarchy construction and voxel
    growth are planned follow-on steps.
    """

    minimal_pore_radius_voxels: float | None = None
    clip_radius_fraction_streamwise: float = 0.05
    clip_radius_fraction_transverse: float = 0.98
    medial_surface_mid_radius_fraction: float = 0.7
    medial_surface_noise_voxels: float | None = None
    hierarchy_length_factor: float = 0.6
    hierarchy_radius_factor: float = 1.1
    radius_smoothing_iterations: int = 3
    retention_radius_factor: float = 0.15
    retention_radius_offset_voxels: float | None = None
    radius_field_mode: str = "half_voxel"
    candidate_selection_mode: str = "threshold_all"

ResolvedMaximalBallSettings dataclass

Concrete maximal-ball settings after Imperial-style default resolution.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class ResolvedMaximalBallSettings:
    """Concrete maximal-ball settings after Imperial-style default resolution."""

    minimal_pore_radius_voxels: float
    clip_radius_fraction_streamwise: float
    clip_radius_fraction_transverse: float
    medial_surface_mid_radius_fraction: float
    medial_surface_noise_voxels: float
    hierarchy_length_factor: float
    hierarchy_radius_factor: float
    radius_smoothing_iterations: int
    retention_radius_factor: float
    retention_radius_offset_voxels: float
    radius_field_mode: str
    candidate_selection_mode: str

MaximalBallCandidates dataclass

Candidate and retained maximal-ball data on voxel centers.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class MaximalBallCandidates:
    """Candidate and retained maximal-ball data on voxel centers."""

    center_indices: np.ndarray
    radii_voxels: np.ndarray
    candidate_mask: np.ndarray
    retained_mask: np.ndarray
    distance_map: np.ndarray
    settings: ResolvedMaximalBallSettings

    @property
    def retained_center_indices(self) -> np.ndarray:
        """Return retained maximal-ball centers in descending-radius order."""

        return cast(np.ndarray, self.center_indices[self.retained_mask])

    @property
    def retained_radii_voxels(self) -> np.ndarray:
        """Return retained maximal-ball radii in descending-radius order."""

        return cast(np.ndarray, self.radii_voxels[self.retained_mask])

retained_center_indices property

retained_center_indices

Return retained maximal-ball centers in descending-radius order.

retained_radii_voxels property

retained_radii_voxels

Return retained maximal-ball radii in descending-radius order.

MaximalBallHierarchy dataclass

Parent-child hierarchy over retained maximal-ball candidates.

The hierarchy is stored on the retained-ball order of :class:MaximalBallCandidates, which is already sorted by descending radius. Each ball points either to itself, if it is a root/master ball, or to the index of its parent ball in the same retained ordering.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class MaximalBallHierarchy:
    """Parent-child hierarchy over retained maximal-ball candidates.

    The hierarchy is stored on the retained-ball order of
    :class:`MaximalBallCandidates`, which is already sorted by descending
    radius. Each ball points either to itself, if it is a root/master ball, or
    to the index of its parent ball in the same retained ordering.
    """

    center_indices: np.ndarray
    center_coordinates: np.ndarray
    radii_voxels: np.ndarray
    parent_indices: np.ndarray
    master_indices: np.ndarray
    hierarchy_levels: np.ndarray
    distance_map: np.ndarray
    settings: ResolvedMaximalBallSettings

    @property
    def root_mask(self) -> np.ndarray:
        """Return a boolean mask of root/master balls."""

        return cast(
            np.ndarray,
            self.parent_indices == np.arange(self.parent_indices.size, dtype=np.int64),
        )

root_mask property

root_mask

Return a boolean mask of root/master balls.

MaximalBallVoxelRegions dataclass

Voxel ownership assignment grown from maximal-ball hierarchy roots.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class MaximalBallVoxelRegions:
    """Voxel ownership assignment grown from maximal-ball hierarchy roots."""

    label_image: np.ndarray
    root_ball_indices: np.ndarray
    root_labels: np.ndarray
    root_center_indices: np.ndarray
    root_radii_voxels: np.ndarray
    root_of_ball_index: np.ndarray
    unassigned_label: int

    @property
    def assigned_void_mask(self) -> np.ndarray:
        """Return a mask of voxels assigned to some pore/root region."""

        return self.label_image >= 0

assigned_void_mask property

assigned_void_mask

Return a mask of voxels assigned to some pore/root region.

MaximalBallRegionAdjacency dataclass

Region-wise geometric summaries derived from voxel ownership labels.

Notes

The fields here are deliberately close to the intermediate quantities that the Imperial extractor builds before CNM export:

• per-region occupied voxel counts
• per-region exposed face counts
• region-to-region interface face counts
• interface centroids in voxel-index coordinates
• boundary-face contact counts on each sample side

This is still an intermediate voxel-geometry product, not yet a final voids.Network.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class MaximalBallRegionAdjacency:
    """Region-wise geometric summaries derived from voxel ownership labels.

    Notes
    -----
    The fields here are deliberately close to the intermediate quantities that
    the Imperial extractor builds before CNM export:

    - per-region occupied voxel counts
    - per-region exposed face counts
    - region-to-region interface face counts
    - interface centroids in voxel-index coordinates
    - boundary-face contact counts on each sample side

    This is still an intermediate voxel-geometry product, not yet a final
    ``voids.Network``.
    """

    region_labels: np.ndarray
    region_volume_voxels: np.ndarray
    region_surface_face_counts: np.ndarray
    throat_region_pairs: np.ndarray
    throat_face_counts: np.ndarray
    throat_axis_face_balance: np.ndarray
    throat_centroid_indices: np.ndarray
    throat_max_touch_radius_side1_voxels: np.ndarray
    throat_max_touch_radius_side2_voxels: np.ndarray
    throat_max_touch_index_side1: np.ndarray
    throat_max_touch_index_side2: np.ndarray
    boundary_face_counts: np.ndarray

MaximalBallExtractionResult dataclass

Staged native maximal-ball extraction outputs before CNM assembly.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class MaximalBallExtractionResult:
    """Staged native maximal-ball extraction outputs before CNM assembly."""

    candidates: MaximalBallCandidates
    hierarchy: MaximalBallHierarchy
    voxel_regions: MaximalBallVoxelRegions
    region_adjacency: MaximalBallRegionAdjacency

MaximalBallNetworkDictResult dataclass

PoreSpy-style network mapping assembled from native maximal-ball regions.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class MaximalBallNetworkDictResult:
    """PoreSpy-style network mapping assembled from native maximal-ball regions."""

    network_dict: dict[str, np.ndarray]
    extraction: MaximalBallExtractionResult

MaximalBallExtractionDiagnostics dataclass

Compact diagnostics for step-by-step maximal-ball extraction comparison.

Source code in src/voids/image/maximal_ball.py

@dataclass(slots=True)
class MaximalBallExtractionDiagnostics:
    """Compact diagnostics for step-by-step maximal-ball extraction comparison."""

    retained_ball_count: int
    root_region_count: int
    occupied_region_count: int
    assigned_void_fraction: float
    unassigned_void_voxel_count: int
    zero_throat_region_count: int
    internal_zero_throat_region_count: int
    boundary_zero_throat_region_count: int
    throat_touch_radius_side1_mean_voxels: float
    throat_touch_radius_side2_mean_voxels: float
    throat_refined_support_radius_side1_mean_voxels: float
    throat_refined_support_radius_side2_mean_voxels: float

compute_void_distance_map

compute_void_distance_map(
    void_phase_mask,
    *,
    backend="auto",
    edt_parallel_threads=None,
)

Compute the void-space Euclidean distance map.

Parameters:

Name	Type	Description	Default
void_phase_mask	ndarray	Boolean array where True marks void voxels.	required
backend	str	Distance-transform backend. "auto" prefers the optional edt package for 3D arrays when available, otherwise falls back to SciPy. Explicit options are "scipy" and "edt".	'auto'
edt_parallel_threads	int | None	Number of worker threads to use when the optional edt backend is active. When omitted, voids first checks VOIDS_EDT_THREADS and otherwise uses one thread for stable default behavior.	None

Source code in src/voids/image/maximal_ball.py

def compute_void_distance_map(
    void_phase_mask: np.ndarray,
    *,
    backend: str = "auto",
    edt_parallel_threads: int | None = None,
) -> np.ndarray:
    """Compute the void-space Euclidean distance map.

    Parameters
    ----------
    void_phase_mask :
        Boolean array where ``True`` marks void voxels.
    backend :
        Distance-transform backend. ``"auto"`` prefers the optional `edt`
        package for 3D arrays when available, otherwise falls back to SciPy.
        Explicit options are ``"scipy"`` and ``"edt"``.
    edt_parallel_threads :
        Number of worker threads to use when the optional `edt` backend is
        active. When omitted, `voids` first checks ``VOIDS_EDT_THREADS`` and
        otherwise uses one thread for stable default behavior.
    """

    mask = np.asarray(void_phase_mask, dtype=bool)
    if mask.ndim not in {2, 3}:
        raise ValueError("void_phase_mask must be a 2D or 3D boolean array")

    normalized_backend = str(backend).strip().lower()
    if normalized_backend not in {"auto", "scipy", "edt"}:
        raise ValueError("backend must be one of {'auto', 'scipy', 'edt'}")

    use_fast_edt = normalized_backend == "edt" or (
        normalized_backend == "auto" and mask.ndim == 3 and fast_edt is not None
    )
    if use_fast_edt:
        if fast_edt is None:
            raise ImportError(
                "backend='edt' requested, but the optional 'edt' package is unavailable"
            )
        return np.asarray(
            fast_edt.edt(
                mask,
                black_border=True,
                parallel=_resolve_edt_parallel_threads(edt_parallel_threads),
            ),
            dtype=float,
        )
    return np.asarray(ndi.distance_transform_edt(mask), dtype=float)

compute_maximal_ball_radius_field

compute_maximal_ball_radius_field(
    void_phase_mask,
    *,
    backend="auto",
    edt_parallel_threads=None,
    mode="half_voxel",
)

Compute the radius field used by the maximal-ball extractor.

Parameters:

Name	Type	Description	Default
void_phase_mask	ndarray	Boolean array where True marks void voxels.	required
backend	str	Distance-transform backend passed through to :func:compute_void_distance_map.	'auto'
edt_parallel_threads	int | None	Number of worker threads to use when the optional edt backend is active.	None
mode	str	Radius-field convention. "half_voxel" returns the nearest non-void Euclidean distance minus half a voxel. "edt" returns the plain Euclidean distance map unchanged. "imperial_pnextract" is accepted as a backward-compatible alias for "half_voxel".	'half_voxel'

Source code in src/voids/image/maximal_ball.py

def compute_maximal_ball_radius_field(
    void_phase_mask: np.ndarray,
    *,
    backend: str = "auto",
    edt_parallel_threads: int | None = None,
    mode: str = "half_voxel",
) -> np.ndarray:
    """Compute the radius field used by the maximal-ball extractor.

    Parameters
    ----------
    void_phase_mask :
        Boolean array where ``True`` marks void voxels.
    backend :
        Distance-transform backend passed through to
        :func:`compute_void_distance_map`.
    edt_parallel_threads :
        Number of worker threads to use when the optional `edt` backend is
        active.
    mode :
        Radius-field convention. ``"half_voxel"`` returns the nearest non-void
        Euclidean distance minus half a voxel. ``"edt"`` returns the plain
        Euclidean distance map unchanged. ``"imperial_pnextract"`` is accepted
        as a backward-compatible alias for ``"half_voxel"``.
    """

    normalized_mode = _normalize_radius_field_mode(mode)

    distance_map = compute_void_distance_map(
        void_phase_mask,
        backend=backend,
        edt_parallel_threads=edt_parallel_threads,
    )
    if normalized_mode == "edt":
        return distance_map
    radius_field = np.asarray(distance_map, dtype=float) - 0.5
    return np.where(radius_field > 0.0, radius_field, 0.0)

smooth_radius_field_local_relaxation

smooth_radius_field_local_relaxation(
    radius_field, void_phase_mask, *, iterations
)

Smooth a radius field with a compact local relaxation stencil.

Source code in src/voids/image/maximal_ball.py

def smooth_radius_field_local_relaxation(
    radius_field: np.ndarray,
    void_phase_mask: np.ndarray,
    *,
    iterations: int,
) -> np.ndarray:
    """Smooth a radius field with a compact local relaxation stencil."""

    if iterations < 0:
        raise ValueError("iterations must be nonnegative")

    mask = np.asarray(void_phase_mask, dtype=bool)
    smoothed_radius = np.asarray(radius_field, dtype=float).copy()
    if smoothed_radius.shape != mask.shape:
        raise ValueError("radius_field and void_phase_mask must have the same shape")
    if smoothed_radius.ndim not in {2, 3}:
        raise ValueError("radius_field must be a 2D or 3D array")
    if iterations == 0:
        return smoothed_radius

    neighborhood_kernel = np.ones((3,) * smoothed_radius.ndim, dtype=float)
    mask_float = mask.astype(float, copy=False)
    for _ in range(iterations):
        neighbor_count = ndi.convolve(
            mask_float,
            neighborhood_kernel,
            mode="constant",
            cval=0.0,
        )
        radius_sum = ndi.convolve(
            np.where(mask, smoothed_radius, 0.0),
            neighborhood_kernel,
            mode="constant",
            cval=0.0,
        )
        radius_delta = np.zeros_like(smoothed_radius, dtype=float)
        radius_delta[mask] = (
            4.0 * radius_sum[mask] / (3.0 * neighbor_count[mask] + 27.0) - smoothed_radius[mask]
        )
        radius_delta_sum = ndi.convolve(
            np.where(mask, radius_delta, 0.0),
            neighborhood_kernel,
            mode="constant",
            cval=0.0,
        )
        local_update = np.zeros_like(smoothed_radius, dtype=float)
        local_update[mask] = 0.02 * (
            radius_delta[mask] - 1.98 * radius_delta_sum[mask] / (neighbor_count[mask] + 27.0)
        )
        smoothed_radius[mask] += np.clip(local_update[mask], -0.005, 0.01)
    return smoothed_radius

resolve_maximal_ball_settings

resolve_maximal_ball_settings(distance_map, settings=None)

Resolve Imperial-style default settings from a distance map.

The Imperial code derives several defaults from the average void-space radius. We mirror that default logic here so staged comparisons use the same parameter semantics even before the full extractor is implemented.

Source code in src/voids/image/maximal_ball.py

def resolve_maximal_ball_settings(
    distance_map: np.ndarray,
    settings: MaximalBallSettings | None = None,
) -> ResolvedMaximalBallSettings:
    """Resolve Imperial-style default settings from a distance map.

    The Imperial code derives several defaults from the average void-space
    radius. We mirror that default logic here so staged comparisons use the
    same parameter semantics even before the full extractor is implemented.
    """

    raw_settings = settings or MaximalBallSettings()
    positive_radii = np.asarray(distance_map, dtype=float)
    positive_radii = positive_radii[positive_radii > 0.0]
    average_radius = float(positive_radii.mean()) if positive_radii.size else 0.0

    default_minimal_radius = min(1.25, 0.25 * average_radius) + 0.5
    minimal_pore_radius_voxels = (
        default_minimal_radius
        if raw_settings.minimal_pore_radius_voxels is None
        else float(raw_settings.minimal_pore_radius_voxels)
    )
    if minimal_pore_radius_voxels <= 0.0:
        raise ValueError("minimal_pore_radius_voxels must be positive")

    medial_surface_noise_voxels = (
        abs(minimal_pore_radius_voxels) + 1.0
        if raw_settings.medial_surface_noise_voxels is None
        else float(raw_settings.medial_surface_noise_voxels)
    )
    retention_radius_offset_voxels = (
        abs(minimal_pore_radius_voxels)
        if raw_settings.retention_radius_offset_voxels is None
        else float(raw_settings.retention_radius_offset_voxels)
    )
    if medial_surface_noise_voxels <= 0.0:
        raise ValueError("medial_surface_noise_voxels must be positive")
    if retention_radius_offset_voxels <= 0.0:
        raise ValueError("retention_radius_offset_voxels must be positive")
    if raw_settings.radius_smoothing_iterations < 0:
        raise ValueError("radius_smoothing_iterations must be nonnegative")
    radius_field_mode = _normalize_radius_field_mode(raw_settings.radius_field_mode)
    candidate_selection_mode = str(raw_settings.candidate_selection_mode).strip().lower()
    if candidate_selection_mode not in {"threshold_all", "local_maxima"}:
        raise ValueError(
            "candidate_selection_mode must be one of {'threshold_all', 'local_maxima'}"
        )

    return ResolvedMaximalBallSettings(
        minimal_pore_radius_voxels=float(minimal_pore_radius_voxels),
        clip_radius_fraction_streamwise=float(raw_settings.clip_radius_fraction_streamwise),
        clip_radius_fraction_transverse=float(raw_settings.clip_radius_fraction_transverse),
        medial_surface_mid_radius_fraction=float(raw_settings.medial_surface_mid_radius_fraction),
        medial_surface_noise_voxels=float(medial_surface_noise_voxels),
        hierarchy_length_factor=float(raw_settings.hierarchy_length_factor),
        hierarchy_radius_factor=float(raw_settings.hierarchy_radius_factor),
        radius_smoothing_iterations=int(raw_settings.radius_smoothing_iterations),
        retention_radius_factor=float(raw_settings.retention_radius_factor),
        retention_radius_offset_voxels=float(retention_radius_offset_voxels),
        radius_field_mode=radius_field_mode,
        candidate_selection_mode=candidate_selection_mode,
    )

clip_distance_map_to_domain_boundaries

clip_distance_map_to_domain_boundaries(
    distance_map, *, settings
)

Apply the Imperial-style boundary clipping heuristic to a distance map.

Source code in src/voids/image/maximal_ball.py

def clip_distance_map_to_domain_boundaries(
    distance_map: np.ndarray,
    *,
    settings: ResolvedMaximalBallSettings,
) -> np.ndarray:
    """Apply the Imperial-style boundary clipping heuristic to a distance map."""

    clipped_distance_map = np.asarray(distance_map, dtype=float).copy()
    if clipped_distance_map.ndim not in {2, 3}:
        raise ValueError("distance_map must be a 2D or 3D array")

    for axis_index, axis_size in enumerate(clipped_distance_map.shape):
        voxel_coordinates = np.arange(axis_size, dtype=float)
        boundary_distance = np.minimum(
            voxel_coordinates + 2.0,
            axis_size - voxel_coordinates + 1.0,
        )
        broadcast_shape = [1] * clipped_distance_map.ndim
        broadcast_shape[axis_index] = axis_size
        boundary_distance = boundary_distance.reshape(broadcast_shape)

        if axis_index == 0:
            clip_fraction = settings.clip_radius_fraction_streamwise
            radius_floor = 0.1
        else:
            clip_fraction = settings.clip_radius_fraction_transverse
            radius_floor = 0.01

        needs_clipping = boundary_distance < clipped_distance_map
        blended_radius = (
            1.0 - clip_fraction
        ) * clipped_distance_map + clip_fraction * boundary_distance
        clipped_distance_map = np.where(
            needs_clipping,
            np.maximum(blended_radius, radius_floor),
            clipped_distance_map,
        )
    return clipped_distance_map

find_maximal_ball_candidates

find_maximal_ball_candidates(
    distance_map,
    *,
    minimal_radius_voxels,
    footprint=None,
    selection_mode="local_maxima",
)

Find maximal-ball candidates from a radius field.

Source code in src/voids/image/maximal_ball.py

def find_maximal_ball_candidates(
    distance_map: np.ndarray,
    *,
    minimal_radius_voxels: float,
    footprint: np.ndarray | None = None,
    selection_mode: str = "local_maxima",
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Find maximal-ball candidates from a radius field."""

    if minimal_radius_voxels <= 0.0:
        raise ValueError("minimal_radius_voxels must be positive")

    working_distance_map = np.asarray(distance_map, dtype=float)
    if working_distance_map.ndim not in {2, 3}:
        raise ValueError("distance_map must be a 2D or 3D array")

    normalized_selection_mode = str(selection_mode).strip().lower()
    if normalized_selection_mode not in {"threshold_all", "local_maxima"}:
        raise ValueError("selection_mode must be one of {'threshold_all', 'local_maxima'}")

    candidate_mask = (working_distance_map > 0.0) & (working_distance_map >= minimal_radius_voxels)
    if normalized_selection_mode == "local_maxima":
        if footprint is None:
            footprint = np.ones((3,) * working_distance_map.ndim, dtype=bool)
        local_maxima = ndi.maximum_filter(
            working_distance_map,
            footprint=footprint,
            mode="nearest",
        )
        candidate_mask &= np.isclose(working_distance_map, local_maxima)

    center_indices = np.argwhere(candidate_mask)
    candidate_radii = working_distance_map[candidate_mask]
    if center_indices.size == 0:
        empty_centers = np.zeros((0, working_distance_map.ndim), dtype=np.int64)
        empty_radii = np.zeros(0, dtype=float)
        return empty_centers, empty_radii, candidate_mask

    descending_order = np.argsort(-candidate_radii, kind="stable")
    return (
        center_indices[descending_order].astype(np.int64, copy=False),
        candidate_radii[descending_order].astype(float, copy=False),
        candidate_mask,
    )

suppress_overlapping_maximal_balls

suppress_overlapping_maximal_balls(
    center_indices, radii_voxels, *, settings
)

Retain descending-radius maximal-ball candidates after overlap suppression.

Source code in src/voids/image/maximal_ball.py

def suppress_overlapping_maximal_balls(
    center_indices: np.ndarray,
    radii_voxels: np.ndarray,
    *,
    settings: ResolvedMaximalBallSettings,
) -> np.ndarray:
    """Retain descending-radius maximal-ball candidates after overlap suppression."""

    centers = np.asarray(center_indices, dtype=np.int64)
    radii = np.asarray(radii_voxels, dtype=float)
    if centers.ndim != 2:
        raise ValueError("center_indices must have shape (count, ndim)")
    if radii.ndim != 1 or radii.shape[0] != centers.shape[0]:
        raise ValueError("radii_voxels must have shape (count,) matching center_indices")
    return np.asarray(
        _suppress_overlapping_maximal_balls_spatial_compiled(
            centers=centers,
            radii=radii,
            retention_radius_factor=float(settings.retention_radius_factor),
            retention_radius_offset_voxels=float(settings.retention_radius_offset_voxels),
            medial_surface_noise_voxels=float(settings.medial_surface_noise_voxels),
            cell_size_voxels=max(
                4.0, min(12.0, 8.0 * float(settings.retention_radius_offset_voxels))
            ),
        ),
        dtype=bool,
    )

refine_retained_ball_coordinates

refine_retained_ball_coordinates(maximal_ball_candidates)

Apply Imperial-style uphill refinements to retained maximal balls.

Returns:

Type	Description
tuple[ndarray, ndarray, ndarray]	Refined integer voxel indices, refined floating-point center coordinates, and refined radii.

Source code in src/voids/image/maximal_ball.py

def refine_retained_ball_coordinates(
    maximal_ball_candidates: MaximalBallCandidates,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Apply Imperial-style uphill refinements to retained maximal balls.

    Returns
    -------
    tuple[np.ndarray, np.ndarray, np.ndarray]
        Refined integer voxel indices, refined floating-point center
        coordinates, and refined radii.
    """

    retained_center_indices = np.asarray(
        maximal_ball_candidates.retained_center_indices,
        dtype=np.int64,
    )
    retained_radii_voxels = np.asarray(
        maximal_ball_candidates.retained_radii_voxels,
        dtype=float,
    )
    retained_count = retained_radii_voxels.size
    if retained_count == 0:
        ndim = maximal_ball_candidates.center_indices.shape[1]
        return (
            np.zeros((0, ndim), dtype=np.int64),
            np.zeros((0, ndim), dtype=float),
            np.zeros(0, dtype=float),
        )

    distance_map = np.asarray(maximal_ball_candidates.distance_map, dtype=float)
    if distance_map.ndim == 3 and retained_center_indices.shape[1] == 3:
        return cast(
            tuple[np.ndarray, np.ndarray, np.ndarray],
            _refine_retained_ball_coordinates_compiled_3d(
                distance_map,
                retained_center_indices,
                retained_radii_voxels,
            ),
        )

    refined_integer_indices = retained_center_indices.copy()
    refined_center_coordinates = retained_center_indices.astype(float) - 0.5
    refined_radii_voxels = retained_radii_voxels.copy()

    for ball_index in range(retained_count):
        refined_center_coordinates[ball_index], refined_radii_voxels[ball_index] = (
            _refine_ball_center_subvoxel(
                distance_map,
                refined_integer_indices[ball_index],
                displacement_limit=0.49,
                radius_gain_factor=0.95,
            )
        )

    occupied_center_lookup = {
        tuple(int(value) for value in integer_index) for integer_index in refined_integer_indices
    }
    for ball_index in range(retained_count):
        original_index_tuple = tuple(int(value) for value in refined_integer_indices[ball_index])
        occupied_center_lookup.discard(original_index_tuple)
        relocated_index, relocated_radius = _refine_ball_center_relocation(
            distance_map,
            refined_integer_indices[ball_index],
            occupied_center_lookup,
        )
        refined_integer_indices[ball_index] = relocated_index
        occupied_center_lookup.add(tuple(int(value) for value in relocated_index))
        refined_radii_voxels[ball_index] = relocated_radius

    for ball_index in range(retained_count):
        refined_center_coordinates[ball_index], refined_radii_voxels[ball_index] = (
            _refine_ball_center_subvoxel(
                distance_map,
                refined_integer_indices[ball_index],
                displacement_limit=0.49,
                radius_gain_factor=0.95,
            )
        )

    return refined_integer_indices, refined_center_coordinates, refined_radii_voxels

refine_ball_from_seed_index

refine_ball_from_seed_index(distance_map, seed_index)

Refine one voxel-centered ball seed with the Imperial uphill sequence.

Source code in src/voids/image/maximal_ball.py

def refine_ball_from_seed_index(
    distance_map: np.ndarray,
    seed_index: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, float]:
    """Refine one voxel-centered ball seed with the Imperial uphill sequence."""

    working_distance_map = np.asarray(distance_map, dtype=float)
    integer_seed_index = np.asarray(seed_index, dtype=np.int64)
    refined_center_coordinate, refined_radius_voxels = _refine_ball_center_subvoxel(
        working_distance_map,
        integer_seed_index,
        displacement_limit=0.49,
        radius_gain_factor=0.95,
    )
    relocated_index, relocated_radius_voxels = _refine_ball_center_relocation(
        working_distance_map,
        integer_seed_index,
        occupied_center_lookup=set(),
    )
    refined_center_coordinate, refined_radius_voxels = _refine_ball_center_subvoxel(
        working_distance_map,
        relocated_index,
        displacement_limit=0.49,
        radius_gain_factor=0.95,
    )
    return (
        relocated_index,
        refined_center_coordinate,
        max(
            float(refined_radius_voxels),
            float(relocated_radius_voxels),
        ),
    )

build_maximal_ball_hierarchy

build_maximal_ball_hierarchy(maximal_ball_candidates)

Build an Imperial-inspired hierarchy over retained maximal balls.

Notes

This stage mirrors the main geometric ideas in the Imperial parent competition logic:

• only retained maximal balls participate
• nearby balls interact only when their midpoint is supported by the void distance map
• smaller balls preferentially attach to larger nearby balls
• nearby master balls can also merge into a higher-level hierarchy

This is still a staged native implementation. The downstream voxel-growth and throat-construction stages are not yet included here.

Source code in src/voids/image/maximal_ball.py

def build_maximal_ball_hierarchy(
    maximal_ball_candidates: MaximalBallCandidates,
) -> MaximalBallHierarchy:
    """Build an Imperial-inspired hierarchy over retained maximal balls.

    Notes
    -----
    This stage mirrors the main geometric ideas in the Imperial parent
    competition logic:

    - only retained maximal balls participate
    - nearby balls interact only when their midpoint is supported by the void
      distance map
    - smaller balls preferentially attach to larger nearby balls
    - nearby master balls can also merge into a higher-level hierarchy

    This is still a staged native implementation. The downstream voxel-growth
    and throat-construction stages are not yet included here.
    """

    retained_center_indices, retained_center_coordinates, retained_radii_voxels = (
        refine_retained_ball_coordinates(maximal_ball_candidates)
    )
    retained_count = retained_radii_voxels.size
    parent_indices = np.arange(retained_count, dtype=np.int64)
    if retained_count == 0:
        return MaximalBallHierarchy(
            center_indices=retained_center_indices,
            center_coordinates=retained_center_coordinates,
            radii_voxels=retained_radii_voxels,
            parent_indices=parent_indices,
            master_indices=parent_indices.copy(),
            hierarchy_levels=np.zeros(0, dtype=np.int64),
            distance_map=np.asarray(maximal_ball_candidates.distance_map, dtype=float),
            settings=maximal_ball_candidates.settings,
        )

    settings = maximal_ball_candidates.settings
    retained_center_coordinates_float = retained_center_coordinates.astype(float, copy=False)
    center_tree = cKDTree(retained_center_coordinates_float)
    distance_map = np.asarray(maximal_ball_candidates.distance_map, dtype=float)
    neighbor_search_radii_voxels = (
        settings.hierarchy_length_factor * retained_radii_voxels
        + 2.0 * settings.medial_surface_noise_voxels
        + 2.0
    )
    nearby_ball_indices_by_ball = center_tree.query_ball_point(
        retained_center_coordinates_float,
        r=neighbor_search_radii_voxels,
        workers=-1,
        return_sorted=False,
    )

    for first_ball_index in range(retained_count):
        first_center_coordinate = retained_center_coordinates[first_ball_index]
        first_radius_voxels = float(retained_radii_voxels[first_ball_index])
        for second_ball_index in nearby_ball_indices_by_ball[first_ball_index]:
            if second_ball_index <= first_ball_index:
                continue
            second_center_coordinate = retained_center_coordinates[second_ball_index]
            second_radius_voxels = float(retained_radii_voxels[second_ball_index])
            if not _pair_has_supported_midpoint(
                first_center_coordinate,
                first_radius_voxels,
                second_center_coordinate,
                second_radius_voxels,
                distance_map=distance_map,
                settings=settings,
            ):
                continue

            larger_ball_index = (
                first_ball_index
                if first_radius_voxels >= second_radius_voxels
                else second_ball_index
            )
            smaller_ball_index = (
                second_ball_index if larger_ball_index == first_ball_index else first_ball_index
            )
            _assign_parent_if_allowed(
                parent_indices,
                smaller_ball_index,
                larger_ball_index,
                radii_voxels=retained_radii_voxels,
            )

            first_root_index = _find_root_index(parent_indices, first_ball_index)
            second_root_index = _find_root_index(parent_indices, second_ball_index)
            if first_root_index == second_root_index:
                continue

            first_root_radius_voxels = float(retained_radii_voxels[first_root_index])
            second_root_radius_voxels = float(retained_radii_voxels[second_root_index])
            first_root_center_index = retained_center_coordinates[first_root_index]
            second_root_center_index = retained_center_coordinates[second_root_index]
            root_distance_voxels = float(
                np.linalg.norm(
                    first_root_center_index.astype(float) - second_root_center_index.astype(float)
                )
            )
            average_root_radius_voxels = 0.5 * (
                first_root_radius_voxels + second_root_radius_voxels
            )
            merge_threshold_voxels = np.sqrt(settings.hierarchy_length_factor) * (
                average_root_radius_voxels + 2.0 * settings.medial_surface_noise_voxels
            )
            if root_distance_voxels > merge_threshold_voxels:
                continue

            larger_root_index = (
                first_root_index
                if first_root_radius_voxels >= second_root_radius_voxels
                else second_root_index
            )
            smaller_root_index = (
                second_root_index if larger_root_index == first_root_index else first_root_index
            )
            if retained_radii_voxels[smaller_root_index] < (
                settings.hierarchy_radius_factor * retained_radii_voxels[smaller_ball_index]
                + settings.medial_surface_noise_voxels
            ):
                _assign_parent_if_allowed(
                    parent_indices,
                    smaller_root_index,
                    larger_root_index,
                    radii_voxels=retained_radii_voxels,
                )

    master_indices = np.array(
        [_find_root_index(parent_indices, ball_index) for ball_index in range(retained_count)],
        dtype=np.int64,
    )
    hierarchy_levels = np.zeros(retained_count, dtype=np.int64)
    for ball_index in range(retained_count):
        current_index = ball_index
        while parent_indices[current_index] != current_index:
            hierarchy_levels[ball_index] += 1
            current_index = int(parent_indices[current_index])

    return MaximalBallHierarchy(
        center_indices=retained_center_indices,
        center_coordinates=retained_center_coordinates,
        radii_voxels=retained_radii_voxels,
        parent_indices=parent_indices,
        master_indices=master_indices,
        hierarchy_levels=hierarchy_levels,
        distance_map=distance_map,
        settings=settings,
    )

initialize_root_region_labels

initialize_root_region_labels(
    void_phase_mask,
    maximal_ball_hierarchy,
    *,
    unassigned_label=-1,
)

Seed voxel-region labels from hierarchy root balls.

Notes

This stage mirrors the first pore-element seeding in the Imperial code: each root/master maximal ball defines an initial pore region, and each retained non-root ball maps to the region of its hierarchy root.

Source code in src/voids/image/maximal_ball.py

def initialize_root_region_labels(
    void_phase_mask: np.ndarray,
    maximal_ball_hierarchy: MaximalBallHierarchy,
    *,
    unassigned_label: int = -1,
) -> MaximalBallVoxelRegions:
    """Seed voxel-region labels from hierarchy root balls.

    Notes
    -----
    This stage mirrors the first pore-element seeding in the Imperial code:
    each root/master maximal ball defines an initial pore region, and each
    retained non-root ball maps to the region of its hierarchy root.
    """

    mask = np.asarray(void_phase_mask, dtype=bool)
    if mask.shape != maximal_ball_hierarchy.distance_map.shape:
        raise ValueError("void_phase_mask must match the hierarchy distance-map shape")

    root_ball_indices = np.flatnonzero(maximal_ball_hierarchy.root_mask).astype(np.int64)
    root_labels = np.arange(root_ball_indices.size, dtype=np.int64)
    label_image = np.full(
        mask.shape,
        unassigned_label,
        dtype=_label_dtype_for_region_count(int(root_ball_indices.size)),
    )
    if root_ball_indices.size == 0:
        return MaximalBallVoxelRegions(
            label_image=label_image,
            root_ball_indices=root_ball_indices,
            root_labels=root_labels,
            root_center_indices=np.zeros((0, mask.ndim), dtype=np.int64),
            root_radii_voxels=np.zeros(0, dtype=float),
            root_of_ball_index=np.zeros(0, dtype=np.int64),
            unassigned_label=int(unassigned_label),
        )

    root_center_indices = maximal_ball_hierarchy.center_indices[root_ball_indices]
    root_radii_voxels = maximal_ball_hierarchy.radii_voxels[root_ball_indices]
    root_lookup = {
        int(ball_index): int(label) for label, ball_index in enumerate(root_ball_indices)
    }
    root_of_ball_index = np.array(
        [root_lookup[int(root_index)] for root_index in maximal_ball_hierarchy.master_indices],
        dtype=np.int64,
    )

    for root_label, root_center_index in zip(root_labels, root_center_indices, strict=False):
        label_image[tuple(int(value) for value in root_center_index)] = int(root_label)

    return MaximalBallVoxelRegions(
        label_image=label_image,
        root_ball_indices=root_ball_indices,
        root_labels=root_labels,
        root_center_indices=root_center_indices,
        root_radii_voxels=root_radii_voxels,
        root_of_ball_index=root_of_ball_index,
        unassigned_label=int(unassigned_label),
    )

seed_root_region_ball_interiors

seed_root_region_ball_interiors(
    void_phase_mask, maximal_ball_hierarchy, voxel_regions
)

Assign small interior neighborhoods around retained balls to their root regions.

Source code in src/voids/image/maximal_ball.py

def seed_root_region_ball_interiors(
    void_phase_mask: np.ndarray,
    maximal_ball_hierarchy: MaximalBallHierarchy,
    voxel_regions: MaximalBallVoxelRegions,
) -> MaximalBallVoxelRegions:
    """Assign small interior neighborhoods around retained balls to their root regions."""

    mask = np.asarray(void_phase_mask, dtype=bool)
    label_image = np.asarray(voxel_regions.label_image).copy()
    retained_center_indices = maximal_ball_hierarchy.center_indices
    retained_radii_voxels = maximal_ball_hierarchy.radii_voxels

    for ball_index, center_index in enumerate(retained_center_indices):
        root_label = int(voxel_regions.root_of_ball_index[ball_index])
        radius_voxels = float(retained_radii_voxels[ball_index])
        seeding_radius_voxels = max(radius_voxels * 0.25 - 1.0, 1.001)
        radius_squared = seeding_radius_voxels * seeding_radius_voxels

        lower_bounds = np.maximum(
            np.floor(center_index - seeding_radius_voxels).astype(np.int64),
            0,
        )
        upper_bounds = np.minimum(
            np.ceil(center_index + seeding_radius_voxels).astype(np.int64) + 1,
            np.asarray(mask.shape, dtype=np.int64),
        )
        index_slices = tuple(
            slice(int(lower), int(upper)) for lower, upper in zip(lower_bounds, upper_bounds)
        )
        candidate_offsets = np.indices(
            tuple(int(upper - lower) for lower, upper in zip(lower_bounds, upper_bounds))
        )
        for axis_index, lower_bound in enumerate(lower_bounds):
            candidate_offsets[axis_index] = (
                candidate_offsets[axis_index] + int(lower_bound) - int(center_index[axis_index])
            )
        candidate_distance_squared = np.sum(candidate_offsets.astype(float) ** 2, axis=0)
        local_mask = candidate_distance_squared <= radius_squared
        local_void_mask = mask[index_slices]
        assignable_mask = local_mask & local_void_mask
        local_labels = label_image[index_slices]
        local_labels[assignable_mask] = root_label
        label_image[index_slices] = local_labels

    return MaximalBallVoxelRegions(
        label_image=label_image,
        root_ball_indices=voxel_regions.root_ball_indices,
        root_labels=voxel_regions.root_labels,
        root_center_indices=voxel_regions.root_center_indices,
        root_radii_voxels=voxel_regions.root_radii_voxels,
        root_of_ball_index=voxel_regions.root_of_ball_index,
        unassigned_label=voxel_regions.unassigned_label,
    )

grow_root_regions_by_radius

grow_root_regions_by_radius(
    void_phase_mask,
    distance_map,
    voxel_regions,
    *,
    minimum_supporting_neighbors,
    radius_support_mode=None,
    require_strictly_larger_radius=None,
    iterations=1,
    void_indices=None,
)

Grow root regions across unassigned void voxels using local radius rules.

Source code in src/voids/image/maximal_ball.py

def grow_root_regions_by_radius(
    void_phase_mask: np.ndarray,
    distance_map: np.ndarray,
    voxel_regions: MaximalBallVoxelRegions,
    *,
    minimum_supporting_neighbors: int,
    radius_support_mode: str | None = None,
    require_strictly_larger_radius: bool | None = None,
    iterations: int = 1,
    void_indices: np.ndarray | None = None,
) -> MaximalBallVoxelRegions:
    """Grow root regions across unassigned void voxels using local radius rules."""

    if minimum_supporting_neighbors < 1:
        raise ValueError("minimum_supporting_neighbors must be at least 1")
    if iterations < 1:
        raise ValueError("iterations must be at least 1")

    mask = np.asarray(void_phase_mask, dtype=bool)
    working_distance_map = np.asarray(distance_map, dtype=float)
    label_image = np.asarray(voxel_regions.label_image).copy()
    if mask.shape != label_image.shape or mask.shape != working_distance_map.shape:
        raise ValueError("void_phase_mask, distance_map, and voxel_regions.label_image must match")

    normalized_radius_support_mode = _normalize_radius_support_mode(
        radius_support_mode=radius_support_mode,
        require_strictly_larger_radius=require_strictly_larger_radius,
    )
    if void_indices is None:
        resolved_void_indices = np.argwhere(mask).astype(np.int64, copy=False)
    else:
        resolved_void_indices = np.asarray(void_indices, dtype=np.int64)
    image_shape = np.asarray(mask.shape, dtype=np.int64)
    radius_support_mode_code = _encode_radius_support_mode(normalized_radius_support_mode)

    for _ in range(iterations):
        previous_labels = label_image.copy()
        previous_labels_flat = previous_labels.reshape(-1)
        label_image_flat = label_image.reshape(-1)
        working_distance_map_flat = working_distance_map.reshape(-1)
        if mask.ndim == 2:
            changed_any_voxel = _grow_root_regions_by_radius_compiled_2d(
                previous_labels_flat,
                label_image_flat,
                working_distance_map_flat,
                resolved_void_indices,
                shape0=int(image_shape[0]),
                shape1=int(image_shape[1]),
                unassigned_label=int(voxel_regions.unassigned_label),
                minimum_supporting_neighbors=int(minimum_supporting_neighbors),
                radius_support_mode_code=int(radius_support_mode_code),
            )
        else:
            changed_any_voxel = _grow_root_regions_by_radius_compiled_3d(
                previous_labels_flat,
                label_image_flat,
                working_distance_map_flat,
                resolved_void_indices,
                shape0=int(image_shape[0]),
                shape1=int(image_shape[1]),
                shape2=int(image_shape[2]),
                unassigned_label=int(voxel_regions.unassigned_label),
                minimum_supporting_neighbors=int(minimum_supporting_neighbors),
                radius_support_mode_code=int(radius_support_mode_code),
            )
        if not changed_any_voxel:
            break

    return MaximalBallVoxelRegions(
        label_image=label_image,
        root_ball_indices=voxel_regions.root_ball_indices,
        root_labels=voxel_regions.root_labels,
        root_center_indices=voxel_regions.root_center_indices,
        root_radii_voxels=voxel_regions.root_radii_voxels,
        root_of_ball_index=voxel_regions.root_of_ball_index,
        unassigned_label=voxel_regions.unassigned_label,
    )

reassign_region_boundary_voxels_by_majority

reassign_region_boundary_voxels_by_majority(
    void_phase_mask,
    distance_map,
    voxel_regions,
    *,
    radius_support_mode="any",
    iterations=1,
    void_indices=None,
)

Reassign weakly supported labeled voxels using a neighbor majority rule.

This mirrors the Imperial medianElem stage conceptually: if a labeled voxel is more exposed to different neighboring pore labels than to its own label, it may be reassigned to the strongest competing neighbor label.

Source code in src/voids/image/maximal_ball.py

def reassign_region_boundary_voxels_by_majority(
    void_phase_mask: np.ndarray,
    distance_map: np.ndarray,
    voxel_regions: MaximalBallVoxelRegions,
    *,
    radius_support_mode: str = "any",
    iterations: int = 1,
    void_indices: np.ndarray | None = None,
) -> MaximalBallVoxelRegions:
    """Reassign weakly supported labeled voxels using a neighbor majority rule.

    This mirrors the Imperial `medianElem` stage conceptually: if a labeled
    voxel is more exposed to different neighboring pore labels than to its own
    label, it may be reassigned to the strongest competing neighbor label.
    """

    if iterations < 1:
        raise ValueError("iterations must be at least 1")

    mask = np.asarray(void_phase_mask, dtype=bool)
    working_distance_map = np.asarray(distance_map, dtype=float)
    label_image = np.asarray(voxel_regions.label_image).copy()
    if mask.shape != label_image.shape or mask.shape != working_distance_map.shape:
        raise ValueError("void_phase_mask, distance_map, and voxel_regions.label_image must match")

    normalized_radius_support_mode = _normalize_radius_support_mode(
        radius_support_mode=radius_support_mode,
        require_strictly_larger_radius=None,
    )
    if void_indices is None:
        resolved_void_indices = np.argwhere(mask).astype(np.int64, copy=False)
    else:
        resolved_void_indices = np.asarray(void_indices, dtype=np.int64)
    image_shape = np.asarray(mask.shape, dtype=np.int64)
    radius_support_mode_code = _encode_radius_support_mode(normalized_radius_support_mode)

    for _ in range(iterations):
        previous_labels = label_image.copy()
        previous_labels_flat = previous_labels.reshape(-1)
        label_image_flat = label_image.reshape(-1)
        working_distance_map_flat = working_distance_map.reshape(-1)
        if mask.ndim == 2:
            changed_any_voxel = _reassign_region_boundary_voxels_by_majority_compiled_2d(
                previous_labels_flat,
                label_image_flat,
                working_distance_map_flat,
                resolved_void_indices,
                shape0=int(image_shape[0]),
                shape1=int(image_shape[1]),
                radius_support_mode_code=int(radius_support_mode_code),
            )
        else:
            changed_any_voxel = _reassign_region_boundary_voxels_by_majority_compiled_3d(
                previous_labels_flat,
                label_image_flat,
                working_distance_map_flat,
                resolved_void_indices,
                shape0=int(image_shape[0]),
                shape1=int(image_shape[1]),
                shape2=int(image_shape[2]),
                radius_support_mode_code=int(radius_support_mode_code),
            )
        if not changed_any_voxel:
            break

    return MaximalBallVoxelRegions(
        label_image=label_image,
        root_ball_indices=voxel_regions.root_ball_indices,
        root_labels=voxel_regions.root_labels,
        root_center_indices=voxel_regions.root_center_indices,
        root_radii_voxels=voxel_regions.root_radii_voxels,
        root_of_ball_index=voxel_regions.root_of_ball_index,
        unassigned_label=voxel_regions.unassigned_label,
    )

retreat_mixed_region_boundary_voxels

retreat_mixed_region_boundary_voxels(
    void_phase_mask, voxel_regions, *, void_indices=None
)

Retreat mixed boundary voxels back to the unassigned state.

This mirrors the Imperial retreatPoresMedian stage: labeled voxels that touch both same-label and different-label neighbors are temporarily removed so later growth passes can rebuild cleaner interfaces.

Source code in src/voids/image/maximal_ball.py

def retreat_mixed_region_boundary_voxels(
    void_phase_mask: np.ndarray,
    voxel_regions: MaximalBallVoxelRegions,
    *,
    void_indices: np.ndarray | None = None,
) -> MaximalBallVoxelRegions:
    """Retreat mixed boundary voxels back to the unassigned state.

    This mirrors the Imperial `retreatPoresMedian` stage: labeled voxels that
    touch both same-label and different-label neighbors are temporarily removed
    so later growth passes can rebuild cleaner interfaces.
    """

    mask = np.asarray(void_phase_mask, dtype=bool)
    label_image = np.asarray(voxel_regions.label_image).copy()
    if mask.shape != label_image.shape:
        raise ValueError("void_phase_mask and voxel_regions.label_image must match")

    image_shape = np.asarray(mask.shape, dtype=np.int64)
    previous_labels = label_image.copy()
    if void_indices is None:
        resolved_void_indices = np.argwhere(mask).astype(np.int64, copy=False)
    else:
        resolved_void_indices = np.asarray(void_indices, dtype=np.int64)
    previous_labels_flat = previous_labels.reshape(-1)
    label_image_flat = label_image.reshape(-1)
    if mask.ndim == 2:
        _retreat_mixed_region_boundary_voxels_compiled_2d(
            previous_labels_flat,
            label_image_flat,
            resolved_void_indices,
            shape0=int(image_shape[0]),
            shape1=int(image_shape[1]),
            unassigned_label=int(voxel_regions.unassigned_label),
        )
    else:
        _retreat_mixed_region_boundary_voxels_compiled_3d(
            previous_labels_flat,
            label_image_flat,
            resolved_void_indices,
            shape0=int(image_shape[0]),
            shape1=int(image_shape[1]),
            shape2=int(image_shape[2]),
            unassigned_label=int(voxel_regions.unassigned_label),
        )

    return MaximalBallVoxelRegions(
        label_image=label_image,
        root_ball_indices=voxel_regions.root_ball_indices,
        root_labels=voxel_regions.root_labels,
        root_center_indices=voxel_regions.root_center_indices,
        root_radii_voxels=voxel_regions.root_radii_voxels,
        root_of_ball_index=voxel_regions.root_of_ball_index,
        unassigned_label=voxel_regions.unassigned_label,
    )

stamp_retained_ball_centers_to_root_labels

stamp_retained_ball_centers_to_root_labels(
    voxel_regions, maximal_ball_hierarchy
)

Restore retained-ball centers to their hierarchy-root region labels.

Source code in src/voids/image/maximal_ball.py

def stamp_retained_ball_centers_to_root_labels(
    voxel_regions: MaximalBallVoxelRegions,
    maximal_ball_hierarchy: MaximalBallHierarchy,
) -> MaximalBallVoxelRegions:
    """Restore retained-ball centers to their hierarchy-root region labels."""

    label_image = np.asarray(voxel_regions.label_image).copy()
    for ball_index, center_index in enumerate(maximal_ball_hierarchy.center_indices):
        root_label = int(voxel_regions.root_of_ball_index[ball_index])
        center_index_tuple = tuple(int(value) for value in center_index)
        label_image[center_index_tuple] = root_label

    return MaximalBallVoxelRegions(
        label_image=label_image,
        root_ball_indices=voxel_regions.root_ball_indices,
        root_labels=voxel_regions.root_labels,
        root_center_indices=voxel_regions.root_center_indices,
        root_radii_voxels=voxel_regions.root_radii_voxels,
        root_of_ball_index=voxel_regions.root_of_ball_index,
        unassigned_label=voxel_regions.unassigned_label,
    )

grow_root_regions_by_neighbor_priority

grow_root_regions_by_neighbor_priority(
    void_phase_mask,
    voxel_regions,
    *,
    iterations=1,
    void_indices=None,
)

Grow unassigned voxels by direct neighbor propagation in sweep order.

This mirrors the late Imperial growPores / growPores_X2 stages more closely than the earlier radius-aware majority passes.

Source code in src/voids/image/maximal_ball.py

def grow_root_regions_by_neighbor_priority(
    void_phase_mask: np.ndarray,
    voxel_regions: MaximalBallVoxelRegions,
    *,
    iterations: int = 1,
    void_indices: np.ndarray | None = None,
) -> MaximalBallVoxelRegions:
    """Grow unassigned voxels by direct neighbor propagation in sweep order.

    This mirrors the late Imperial `growPores` / `growPores_X2` stages more
    closely than the earlier radius-aware majority passes.
    """

    if iterations < 1:
        raise ValueError("iterations must be at least 1")

    mask = np.asarray(void_phase_mask, dtype=bool)
    label_image = np.asarray(voxel_regions.label_image).copy()
    if mask.shape != label_image.shape:
        raise ValueError("void_phase_mask and voxel_regions.label_image must match")

    image_shape = np.asarray(mask.shape, dtype=np.int64)
    if void_indices is None:
        forward_indices = np.argwhere(mask).astype(np.int64, copy=False)
    else:
        forward_indices = np.asarray(void_indices, dtype=np.int64)
    backward_indices = forward_indices[::-1].copy()
    label_image_flat = label_image.reshape(-1)
    if mask.ndim == 2:
        _grow_root_regions_by_neighbor_priority_compiled_2d(
            label_image_flat,
            forward_indices,
            backward_indices,
            shape0=int(image_shape[0]),
            shape1=int(image_shape[1]),
            unassigned_label=int(voxel_regions.unassigned_label),
            iterations=int(iterations),
        )
    else:
        _grow_root_regions_by_neighbor_priority_compiled_3d(
            label_image_flat,
            forward_indices,
            backward_indices,
            shape0=int(image_shape[0]),
            shape1=int(image_shape[1]),
            shape2=int(image_shape[2]),
            unassigned_label=int(voxel_regions.unassigned_label),
            iterations=int(iterations),
        )

    return MaximalBallVoxelRegions(
        label_image=label_image,
        root_ball_indices=voxel_regions.root_ball_indices,
        root_labels=voxel_regions.root_labels,
        root_center_indices=voxel_regions.root_center_indices,
        root_radii_voxels=voxel_regions.root_radii_voxels,
        root_of_ball_index=voxel_regions.root_of_ball_index,
        unassigned_label=voxel_regions.unassigned_label,
    )

assign_voxel_regions_from_hierarchy

assign_voxel_regions_from_hierarchy(
    void_phase_mask, maximal_ball_hierarchy
)

Assign voxel ownership using an Imperial-inspired staged growth schedule.

Source code in src/voids/image/maximal_ball.py

def assign_voxel_regions_from_hierarchy(
    void_phase_mask: np.ndarray,
    maximal_ball_hierarchy: MaximalBallHierarchy,
) -> MaximalBallVoxelRegions:
    """Assign voxel ownership using an Imperial-inspired staged growth schedule."""

    mask = np.asarray(void_phase_mask, dtype=bool)
    void_indices = np.argwhere(mask).astype(np.int64, copy=False)
    voxel_regions = initialize_root_region_labels(
        mask,
        maximal_ball_hierarchy,
    )
    voxel_regions = seed_root_region_ball_interiors(
        mask,
        maximal_ball_hierarchy,
        voxel_regions,
    )
    initial_growth_schedule = [
        (3, "strictly_larger", 3),
        (2, "strictly_larger", 3),
        (2, "greater_or_equal", 4),
    ]
    for minimum_supporting_neighbors, radius_support_mode, iterations in initial_growth_schedule:
        voxel_regions = grow_root_regions_by_radius(
            mask,
            maximal_ball_hierarchy.distance_map,
            voxel_regions,
            minimum_supporting_neighbors=minimum_supporting_neighbors,
            radius_support_mode=radius_support_mode,
            iterations=iterations,
            void_indices=void_indices,
        )

    voxel_regions = reassign_region_boundary_voxels_by_majority(
        mask,
        maximal_ball_hierarchy.distance_map,
        voxel_regions,
        radius_support_mode="any",
        iterations=2,
        void_indices=void_indices,
    )
    voxel_regions = retreat_mixed_region_boundary_voxels(
        mask,
        voxel_regions,
        void_indices=void_indices,
    )
    voxel_regions = stamp_retained_ball_centers_to_root_labels(
        voxel_regions,
        maximal_ball_hierarchy,
    )

    late_growth_schedule = [
        (2, "greater_or_equal", 6),
        (2, "any", 4),
        (1, "any", 4),
    ]
    for minimum_supporting_neighbors, radius_support_mode, iterations in late_growth_schedule:
        voxel_regions = grow_root_regions_by_radius(
            mask,
            maximal_ball_hierarchy.distance_map,
            voxel_regions,
            minimum_supporting_neighbors=minimum_supporting_neighbors,
            radius_support_mode=radius_support_mode,
            iterations=iterations,
            void_indices=void_indices,
        )

    voxel_regions = reassign_region_boundary_voxels_by_majority(
        mask,
        maximal_ball_hierarchy.distance_map,
        voxel_regions,
        radius_support_mode="any",
        iterations=2,
        void_indices=void_indices,
    )
    voxel_regions = grow_root_regions_by_neighbor_priority(
        mask,
        voxel_regions,
        iterations=max(16, 2 * max(mask.shape)),
        void_indices=void_indices,
    )
    return voxel_regions

measure_region_adjacency

measure_region_adjacency(
    void_phase_mask, voxel_regions, *, distance_map=None
)

Measure pore-region volumes, interfaces, and boundary contacts.

Parameters:

Name	Type	Description	Default
void_phase_mask	ndarray	Boolean void-domain mask used for extraction.	required
voxel_regions	MaximalBallVoxelRegions	Labeled pore/root ownership image.	required

Returns:

Type	Description
MaximalBallRegionAdjacency	Region-wise voxel volumes and region-pair interface measurements.

Notes

This stage converts the voxel partition into the basic discrete geometry we need for a native pnextract-like network assembly:

• region voxel counts become pore-region volumes
• region-pair contact faces become throat candidates
• boundary-face contacts expose inlet/outlet touching regions

The centroid coordinates are reported in voxel-index units, using face midpoint locations such as i + 0.5 along the axis normal to the interface.

Source code in src/voids/image/maximal_ball.py

def measure_region_adjacency(
    void_phase_mask: np.ndarray,
    voxel_regions: MaximalBallVoxelRegions,
    *,
    distance_map: np.ndarray | None = None,
) -> MaximalBallRegionAdjacency:
    """Measure pore-region volumes, interfaces, and boundary contacts.

    Parameters
    ----------
    void_phase_mask :
        Boolean void-domain mask used for extraction.
    voxel_regions :
        Labeled pore/root ownership image.

    Returns
    -------
    MaximalBallRegionAdjacency
        Region-wise voxel volumes and region-pair interface measurements.

    Notes
    -----
    This stage converts the voxel partition into the basic discrete geometry we
    need for a native `pnextract`-like network assembly:

    - region voxel counts become pore-region volumes
    - region-pair contact faces become throat candidates
    - boundary-face contacts expose inlet/outlet touching regions

    The centroid coordinates are reported in voxel-index units, using face
    midpoint locations such as ``i + 0.5`` along the axis normal to the
    interface.
    """

    mask = np.asarray(void_phase_mask, dtype=bool)
    label_image = np.asarray(voxel_regions.label_image)
    if mask.shape != label_image.shape:
        raise ValueError("void_phase_mask and voxel_regions.label_image must match")
    working_distance_map: np.ndarray | None = None
    if distance_map is not None:
        working_distance_map = np.asarray(distance_map, dtype=float)
        if working_distance_map.shape != mask.shape:
            raise ValueError(
                "distance_map must match void_phase_mask and voxel_regions.label_image"
            )

    region_labels = np.asarray(voxel_regions.root_labels, dtype=np.int64)
    region_count = int(region_labels.size)
    region_volume_voxels = np.zeros(region_count, dtype=np.int64)
    region_surface_face_counts = np.zeros(region_count, dtype=np.int64)
    boundary_face_counts = np.zeros((region_count, 2 * mask.ndim), dtype=np.int64)

    assigned_void_mask = mask & (label_image >= 0)
    if np.any(label_image[assigned_void_mask] >= region_count):
        raise ValueError("voxel region labels must be contiguous root labels starting at zero")

    if np.any(assigned_void_mask):
        region_volume_voxels += np.bincount(
            label_image[assigned_void_mask],
            minlength=region_count,
        ).astype(np.int64, copy=False)

    pair_face_counts: dict[tuple[int, int], int] = {}
    pair_axis_face_balance: dict[tuple[int, int], np.ndarray] = {}
    pair_centroid_sums: dict[tuple[int, int], np.ndarray] = {}
    pair_max_touch_radius_side1: dict[tuple[int, int], float] = {}
    pair_max_touch_radius_side2: dict[tuple[int, int], float] = {}
    pair_max_touch_index_side1: dict[tuple[int, int], np.ndarray] = {}
    pair_max_touch_index_side2: dict[tuple[int, int], np.ndarray] = {}

    image_shape = np.asarray(mask.shape, dtype=np.int64)
    for axis_index in range(mask.ndim):
        lower_boundary_selector: list[int | slice] = [slice(None)] * mask.ndim
        lower_boundary_selector[axis_index] = 0
        lower_boundary_void_mask = assigned_void_mask[tuple(lower_boundary_selector)]
        lower_boundary_labels = label_image[tuple(lower_boundary_selector)]
        if np.any(lower_boundary_void_mask):
            lower_counts = np.bincount(
                lower_boundary_labels[lower_boundary_void_mask],
                minlength=region_count,
            )
            boundary_face_counts[:, 2 * axis_index] += lower_counts.astype(np.int64, copy=False)
            region_surface_face_counts += lower_counts.astype(np.int64, copy=False)

        upper_boundary_selector: list[int | slice] = [slice(None)] * mask.ndim
        upper_boundary_selector[axis_index] = int(image_shape[axis_index] - 1)
        upper_boundary_void_mask = assigned_void_mask[tuple(upper_boundary_selector)]
        upper_boundary_labels = label_image[tuple(upper_boundary_selector)]
        if np.any(upper_boundary_void_mask):
            upper_counts = np.bincount(
                upper_boundary_labels[upper_boundary_void_mask],
                minlength=region_count,
            )
            boundary_face_counts[:, 2 * axis_index + 1] += upper_counts.astype(
                np.int64,
                copy=False,
            )
            region_surface_face_counts += upper_counts.astype(np.int64, copy=False)

        lower_slices = [slice(None)] * mask.ndim
        upper_slices = [slice(None)] * mask.ndim
        lower_slices[axis_index] = slice(0, -1)
        upper_slices[axis_index] = slice(1, None)
        lower_slices_tuple = tuple(lower_slices)
        upper_slices_tuple = tuple(upper_slices)

        lower_assigned_mask = assigned_void_mask[lower_slices_tuple]
        upper_assigned_mask = assigned_void_mask[upper_slices_tuple]
        lower_labels = label_image[lower_slices_tuple]
        upper_labels = label_image[upper_slices_tuple]

        differing_assignment_mask = lower_assigned_mask & (
            (~upper_assigned_mask) | (lower_labels != upper_labels)
        )
        if np.any(differing_assignment_mask):
            lower_surface_counts = np.bincount(
                lower_labels[differing_assignment_mask],
                minlength=region_count,
            )
            region_surface_face_counts += lower_surface_counts.astype(np.int64, copy=False)

        differing_assignment_mask = upper_assigned_mask & (
            (~lower_assigned_mask) | (lower_labels != upper_labels)
        )
        if np.any(differing_assignment_mask):
            upper_surface_counts = np.bincount(
                upper_labels[differing_assignment_mask],
                minlength=region_count,
            )
            region_surface_face_counts += upper_surface_counts.astype(np.int64, copy=False)

        shared_interface_mask = (
            lower_assigned_mask
            & upper_assigned_mask
            & (lower_labels >= 0)
            & (upper_labels >= 0)
            & (lower_labels != upper_labels)
        )
        if not np.any(shared_interface_mask):
            continue

        lower_interface_labels = lower_labels[shared_interface_mask]
        upper_interface_labels = upper_labels[shared_interface_mask]
        interface_indices = np.argwhere(shared_interface_mask).astype(float, copy=False)
        face_midpoint_indices = interface_indices.copy()
        face_midpoint_indices[:, axis_index] += 0.5
        lower_touch_radii = (
            working_distance_map[lower_slices_tuple][shared_interface_mask]
            if working_distance_map is not None
            else None
        )
        upper_touch_radii = (
            working_distance_map[upper_slices_tuple][shared_interface_mask]
            if working_distance_map is not None
            else None
        )

        for face_index in range(interface_indices.shape[0]):
            first_label = int(lower_interface_labels[face_index])
            second_label = int(upper_interface_labels[face_index])
            lower_voxel_index = interface_indices[face_index].astype(np.int64, copy=False)
            upper_voxel_index = lower_voxel_index.copy()
            upper_voxel_index[axis_index] += 1
            if first_label < second_label:
                ordered_pair = (first_label, second_label)
                orientation_sign = 1.0
                side1_touch_radius = (
                    float(lower_touch_radii[face_index])
                    if lower_touch_radii is not None
                    else float("nan")
                )
                side2_touch_radius = (
                    float(upper_touch_radii[face_index])
                    if upper_touch_radii is not None
                    else float("nan")
                )
                side1_touch_index = lower_voxel_index
                side2_touch_index = upper_voxel_index
            else:
                ordered_pair = (second_label, first_label)
                orientation_sign = -1.0
                side1_touch_radius = (
                    float(upper_touch_radii[face_index])
                    if upper_touch_radii is not None
                    else float("nan")
                )
                side2_touch_radius = (
                    float(lower_touch_radii[face_index])
                    if lower_touch_radii is not None
                    else float("nan")
                )
                side1_touch_index = upper_voxel_index
                side2_touch_index = lower_voxel_index

            pair_face_counts[ordered_pair] = pair_face_counts.get(ordered_pair, 0) + 1
            if ordered_pair not in pair_axis_face_balance:
                pair_axis_face_balance[ordered_pair] = np.zeros(mask.ndim, dtype=float)
            if ordered_pair not in pair_centroid_sums:
                pair_centroid_sums[ordered_pair] = np.zeros(mask.ndim, dtype=float)
            pair_axis_face_balance[ordered_pair][axis_index] += orientation_sign
            pair_centroid_sums[ordered_pair] += face_midpoint_indices[face_index]
            if np.isfinite(side1_touch_radius):
                previous_side1_radius = pair_max_touch_radius_side1.get(ordered_pair, -np.inf)
                if side1_touch_radius > previous_side1_radius:
                    pair_max_touch_radius_side1[ordered_pair] = side1_touch_radius
                    pair_max_touch_index_side1[ordered_pair] = side1_touch_index.copy()
            if np.isfinite(side2_touch_radius):
                previous_side2_radius = pair_max_touch_radius_side2.get(ordered_pair, -np.inf)
                if side2_touch_radius > previous_side2_radius:
                    pair_max_touch_radius_side2[ordered_pair] = side2_touch_radius
                    pair_max_touch_index_side2[ordered_pair] = side2_touch_index.copy()

    occupied_region_mask = region_volume_voxels > 0
    occupied_region_indices = np.flatnonzero(occupied_region_mask).astype(np.int64)
    if occupied_region_indices.size != region_count:
        compact_label_of_region = np.full(region_count, -1, dtype=np.int64)
        compact_label_of_region[occupied_region_indices] = np.arange(
            occupied_region_indices.size,
            dtype=np.int64,
        )
        region_labels = region_labels[occupied_region_indices]
        region_volume_voxels = region_volume_voxels[occupied_region_indices]
        region_surface_face_counts = region_surface_face_counts[occupied_region_indices]
        boundary_face_counts = boundary_face_counts[occupied_region_indices]
        remapped_pair_face_counts: dict[tuple[int, int], int] = {}
        remapped_pair_axis_face_balance: dict[tuple[int, int], np.ndarray] = {}
        remapped_pair_centroid_sums: dict[tuple[int, int], np.ndarray] = {}
        remapped_pair_max_touch_radius_side1: dict[tuple[int, int], float] = {}
        remapped_pair_max_touch_radius_side2: dict[tuple[int, int], float] = {}
        remapped_pair_max_touch_index_side1: dict[tuple[int, int], np.ndarray] = {}
        remapped_pair_max_touch_index_side2: dict[tuple[int, int], np.ndarray] = {}
        for ordered_pair, face_count in pair_face_counts.items():
            first_region = int(compact_label_of_region[ordered_pair[0]])
            second_region = int(compact_label_of_region[ordered_pair[1]])
            if first_region < 0 or second_region < 0:
                continue
            remapped_pair = (first_region, second_region)
            remapped_pair_face_counts[remapped_pair] = face_count
            remapped_pair_axis_face_balance[remapped_pair] = pair_axis_face_balance[ordered_pair]
            remapped_pair_centroid_sums[remapped_pair] = pair_centroid_sums[ordered_pair]
            if ordered_pair in pair_max_touch_radius_side1:
                remapped_pair_max_touch_radius_side1[remapped_pair] = pair_max_touch_radius_side1[
                    ordered_pair
                ]
                remapped_pair_max_touch_index_side1[remapped_pair] = pair_max_touch_index_side1[
                    ordered_pair
                ].copy()
            if ordered_pair in pair_max_touch_radius_side2:
                remapped_pair_max_touch_radius_side2[remapped_pair] = pair_max_touch_radius_side2[
                    ordered_pair
                ]
                remapped_pair_max_touch_index_side2[remapped_pair] = pair_max_touch_index_side2[
                    ordered_pair
                ].copy()
        pair_face_counts = remapped_pair_face_counts
        pair_axis_face_balance = remapped_pair_axis_face_balance
        pair_centroid_sums = remapped_pair_centroid_sums
        pair_max_touch_radius_side1 = remapped_pair_max_touch_radius_side1
        pair_max_touch_radius_side2 = remapped_pair_max_touch_radius_side2
        pair_max_touch_index_side1 = remapped_pair_max_touch_index_side1
        pair_max_touch_index_side2 = remapped_pair_max_touch_index_side2

    ordered_pairs = sorted(pair_face_counts)
    throat_region_pairs = np.asarray(ordered_pairs, dtype=np.int64)
    throat_count = len(ordered_pairs)
    throat_face_counts = np.zeros(throat_count, dtype=np.int64)
    throat_axis_face_balance = np.zeros((throat_count, mask.ndim), dtype=float)
    throat_centroid_indices = np.zeros((throat_count, mask.ndim), dtype=float)
    throat_max_touch_radius_side1_voxels = np.full(throat_count, np.nan, dtype=float)
    throat_max_touch_radius_side2_voxels = np.full(throat_count, np.nan, dtype=float)
    throat_max_touch_index_side1 = np.full((throat_count, mask.ndim), -1, dtype=np.int64)
    throat_max_touch_index_side2 = np.full((throat_count, mask.ndim), -1, dtype=np.int64)
    for throat_index, ordered_pair in enumerate(ordered_pairs):
        face_count = int(pair_face_counts[ordered_pair])
        throat_face_counts[throat_index] = face_count
        throat_axis_face_balance[throat_index] = pair_axis_face_balance[ordered_pair]
        throat_centroid_indices[throat_index] = pair_centroid_sums[ordered_pair] / max(
            face_count, 1
        )
        if ordered_pair in pair_max_touch_radius_side1:
            throat_max_touch_radius_side1_voxels[throat_index] = pair_max_touch_radius_side1[
                ordered_pair
            ]
            throat_max_touch_index_side1[throat_index] = pair_max_touch_index_side1[ordered_pair]
        if ordered_pair in pair_max_touch_radius_side2:
            throat_max_touch_radius_side2_voxels[throat_index] = pair_max_touch_radius_side2[
                ordered_pair
            ]
            throat_max_touch_index_side2[throat_index] = pair_max_touch_index_side2[ordered_pair]

    return MaximalBallRegionAdjacency(
        region_labels=region_labels,
        region_volume_voxels=region_volume_voxels,
        region_surface_face_counts=region_surface_face_counts,
        throat_region_pairs=throat_region_pairs,
        throat_face_counts=throat_face_counts,
        throat_axis_face_balance=throat_axis_face_balance,
        throat_centroid_indices=throat_centroid_indices,
        throat_max_touch_radius_side1_voxels=throat_max_touch_radius_side1_voxels,
        throat_max_touch_radius_side2_voxels=throat_max_touch_radius_side2_voxels,
        throat_max_touch_index_side1=throat_max_touch_index_side1,
        throat_max_touch_index_side2=throat_max_touch_index_side2,
        boundary_face_counts=boundary_face_counts,
    )

extract_maximal_ball_regions

extract_maximal_ball_regions(
    void_phase_mask,
    *,
    distance_map_backend="auto",
    edt_parallel_threads=None,
    settings=None,
    apply_boundary_clipping=True,
)

Run the staged native maximal-ball pipeline up to region adjacency.

This is the current highest-level native extraction entry point that stays independent of PoreSpy network generation. It stops at voxel-region and interface geometry because the final pore/throat-to-Network assembly is still under active implementation.

Source code in src/voids/image/maximal_ball.py

def extract_maximal_ball_regions(
    void_phase_mask: np.ndarray,
    *,
    distance_map_backend: str = "auto",
    edt_parallel_threads: int | None = None,
    settings: MaximalBallSettings | None = None,
    apply_boundary_clipping: bool = True,
) -> MaximalBallExtractionResult:
    """Run the staged native maximal-ball pipeline up to region adjacency.

    This is the current highest-level native extraction entry point that stays
    independent of PoreSpy network generation. It stops at voxel-region and
    interface geometry because the final pore/throat-to-`Network` assembly is
    still under active implementation.
    """

    candidates = extract_maximal_ball_candidates(
        void_phase_mask,
        distance_map_backend=distance_map_backend,
        edt_parallel_threads=edt_parallel_threads,
        settings=settings,
        apply_boundary_clipping=apply_boundary_clipping,
    )
    hierarchy = build_maximal_ball_hierarchy(candidates)
    voxel_regions = assign_voxel_regions_from_hierarchy(void_phase_mask, hierarchy)
    region_adjacency = measure_region_adjacency(
        void_phase_mask,
        voxel_regions,
        distance_map=hierarchy.distance_map,
    )
    return MaximalBallExtractionResult(
        candidates=candidates,
        hierarchy=hierarchy,
        voxel_regions=voxel_regions,
        region_adjacency=region_adjacency,
    )

summarize_maximal_ball_extraction_diagnostics

summarize_maximal_ball_extraction_diagnostics(
    void_phase_mask, extraction_result
)

Summarize intermediate maximal-ball extraction behavior for comparison work.

Source code in src/voids/image/maximal_ball.py

def summarize_maximal_ball_extraction_diagnostics(
    void_phase_mask: np.ndarray,
    extraction_result: MaximalBallExtractionResult,
) -> MaximalBallExtractionDiagnostics:
    """Summarize intermediate maximal-ball extraction behavior for comparison work."""

    mask = np.asarray(void_phase_mask, dtype=bool)
    label_image = np.asarray(extraction_result.voxel_regions.label_image)
    if mask.shape != label_image.shape:
        raise ValueError("void_phase_mask must match extraction_result voxel-region labels")

    assigned_void_mask = mask & (label_image >= 0)
    unassigned_void_voxel_count = int(np.count_nonzero(mask & (label_image < 0)))
    region_adjacency = extraction_result.region_adjacency
    (
        refined_support_radius_side1_voxels,
        refined_support_radius_side2_voxels,
    ) = _select_interface_supporting_ball_radii(extraction_result)
    occupied_region_count = int(region_adjacency.region_volume_voxels.size)
    neighbor_counts = np.zeros(occupied_region_count, dtype=np.int64)
    throat_region_pairs = np.asarray(region_adjacency.throat_region_pairs, dtype=np.int64)
    if throat_region_pairs.size:
        np.add.at(neighbor_counts, throat_region_pairs[:, 0], 1)
        np.add.at(neighbor_counts, throat_region_pairs[:, 1], 1)
    zero_throat_region_mask = neighbor_counts == 0
    boundary_face_counts = np.asarray(region_adjacency.boundary_face_counts, dtype=np.int64)
    boundary_touch_mask = np.sum(boundary_face_counts, axis=1) > 0
    assigned_void_count = int(np.count_nonzero(assigned_void_mask))
    void_voxel_count = max(int(np.count_nonzero(mask)), 1)

    return MaximalBallExtractionDiagnostics(
        retained_ball_count=int(extraction_result.candidates.retained_center_indices.shape[0]),
        root_region_count=int(extraction_result.voxel_regions.root_labels.size),
        occupied_region_count=occupied_region_count,
        assigned_void_fraction=float(assigned_void_count / void_voxel_count),
        unassigned_void_voxel_count=unassigned_void_voxel_count,
        zero_throat_region_count=int(np.count_nonzero(zero_throat_region_mask)),
        internal_zero_throat_region_count=int(
            np.count_nonzero(zero_throat_region_mask & ~boundary_touch_mask)
        ),
        boundary_zero_throat_region_count=int(
            np.count_nonzero(zero_throat_region_mask & boundary_touch_mask)
        ),
        throat_touch_radius_side1_mean_voxels=float(
            np.nanmean(region_adjacency.throat_max_touch_radius_side1_voxels)
            if region_adjacency.throat_max_touch_radius_side1_voxels.size
            else np.nan
        ),
        throat_touch_radius_side2_mean_voxels=float(
            np.nanmean(region_adjacency.throat_max_touch_radius_side2_voxels)
            if region_adjacency.throat_max_touch_radius_side2_voxels.size
            else np.nan
        ),
        throat_refined_support_radius_side1_mean_voxels=float(
            np.nanmean(refined_support_radius_side1_voxels)
            if refined_support_radius_side1_voxels.size
            else np.nan
        ),
        throat_refined_support_radius_side2_mean_voxels=float(
            np.nanmean(refined_support_radius_side2_voxels)
            if refined_support_radius_side2_voxels.size
            else np.nan
        ),
    )

build_network_dict_from_maximal_ball_regions

build_network_dict_from_maximal_ball_regions(
    extraction_result,
    *,
    voxel_size,
    axis_names=("x", "y", "z"),
    flow_boundary_mode="direct",
    boundary_axis=None,
    boundary_length_epsilon=1e-300,
    boundary_radius_scale=1.1,
    throat_area_mode="face_count",
    throat_shape_factor_radius_mode="inscribed",
    throat_anchor_mode="second_side",
)

Assemble a PoreSpy-style network mapping from maximal-ball regions.

Parameters:

Name	Type	Description	Default
extraction_result	MaximalBallExtractionResult	Native maximal-ball extraction outputs through the region-adjacency stage.	required
voxel_size	float	Physical edge length of one voxel.	required
axis_names	tuple[str, ...]	Axis labels associated with the image dimensions. Only the first ndim entries are used.	('x', 'y', 'z')

Notes

This builder intentionally uses explicit, readable geometric rules rather than hidden heuristics:

• pore coordinates are the root maximal-ball centers
• pore volumes are the labeled region voxel counts
• throat areas are the counted interface faces
• throat centroids are the mean interface-face midpoints
• conduit lengths are derived from pore-center to interface-centroid distances with a minimum half-voxel regularization

The current implementation now follows the Imperial export logic more closely in three places:

• throat radii are taken from interface-supporting maximal balls
• throat and pore shape factors follow the Imperial export repair and throat-area-weighted pore averaging logic
• region volumes are redistributed between pores and throats using the same area-partition rule used by the Imperial CNM writer

This is still not full pnextract parity, because the upstream voxel ownership and throat-surface ball construction remain a native approximation.

Source code in src/voids/image/maximal_ball.py

def build_network_dict_from_maximal_ball_regions(
    extraction_result: MaximalBallExtractionResult,
    *,
    voxel_size: float,
    axis_names: tuple[str, ...] = ("x", "y", "z"),
    flow_boundary_mode: str = "direct",
    boundary_axis: str | None = None,
    boundary_length_epsilon: float = 1.0e-300,
    boundary_radius_scale: float = 1.1,
    throat_area_mode: str = "face_count",
    throat_shape_factor_radius_mode: str = "inscribed",
    throat_anchor_mode: str = "second_side",
) -> dict[str, np.ndarray]:
    """Assemble a PoreSpy-style network mapping from maximal-ball regions.

    Parameters
    ----------
    extraction_result :
        Native maximal-ball extraction outputs through the region-adjacency
        stage.
    voxel_size :
        Physical edge length of one voxel.
    axis_names :
        Axis labels associated with the image dimensions. Only the first
        ``ndim`` entries are used.

    Notes
    -----
    This builder intentionally uses explicit, readable geometric rules rather
    than hidden heuristics:

    - pore coordinates are the root maximal-ball centers
    - pore volumes are the labeled region voxel counts
    - throat areas are the counted interface faces
    - throat centroids are the mean interface-face midpoints
    - conduit lengths are derived from pore-center to interface-centroid
      distances with a minimum half-voxel regularization

    The current implementation now follows the Imperial export logic more
    closely in three places:

    - throat radii are taken from interface-supporting maximal balls
    - throat and pore shape factors follow the Imperial export repair and
      throat-area-weighted pore averaging logic
    - region volumes are redistributed between pores and throats using the same
      area-partition rule used by the Imperial CNM writer

    This is still not full `pnextract` parity, because the upstream voxel
    ownership and throat-surface ball construction remain a native
    approximation.
    """

    if voxel_size <= 0.0:
        raise ValueError("voxel_size must be positive")
    normalized_flow_boundary_mode = _resolve_flow_boundary_mode(flow_boundary_mode)
    normalized_throat_area_mode = _resolve_throat_area_mode(throat_area_mode)
    normalized_throat_shape_factor_radius_mode = _resolve_throat_shape_factor_radius_mode(
        throat_shape_factor_radius_mode
    )
    normalized_throat_anchor_mode = _resolve_throat_anchor_mode(throat_anchor_mode)
    if boundary_length_epsilon <= 0.0:
        raise ValueError("boundary_length_epsilon must be positive")
    if boundary_radius_scale <= 0.0:
        raise ValueError("boundary_radius_scale must be positive")

    hierarchy = extraction_result.hierarchy
    voxel_regions = extraction_result.voxel_regions
    region_adjacency = extraction_result.region_adjacency
    image_ndim = (
        int(hierarchy.center_indices.shape[1])
        if hierarchy.center_indices.size
        else int(voxel_regions.label_image.ndim)
    )
    if image_ndim not in {2, 3}:
        raise ValueError("maximal-ball network assembly supports only 2D or 3D images")
    if len(axis_names) < image_ndim:
        raise ValueError("axis_names must provide at least one label per image dimension")

    active_axis_names = axis_names[:image_ndim]
    if boundary_axis is not None and boundary_axis not in active_axis_names:
        raise ValueError("boundary_axis must be one of the active axis names")
    if normalized_flow_boundary_mode == "external_reservoir" and boundary_axis is None:
        boundary_axis = active_axis_names[0]
    retained_region_labels = np.asarray(region_adjacency.region_labels, dtype=np.int64)
    region_count = int(retained_region_labels.size)
    all_root_center_indices = np.asarray(voxel_regions.root_center_indices, dtype=float)
    all_root_radii_voxels = np.asarray(voxel_regions.root_radii_voxels, dtype=float)

    if all_root_center_indices.ndim != 2 or all_root_center_indices.shape[1] != image_ndim:
        raise ValueError("root_center_indices must have shape (n_roots, image_ndim)")
    if all_root_radii_voxels.shape != (all_root_center_indices.shape[0],):
        raise ValueError("root_radii_voxels must align with root_center_indices")
    if np.any(retained_region_labels < 0) or np.any(
        retained_region_labels >= all_root_center_indices.shape[0]
    ):
        raise ValueError("region_adjacency.region_labels must index the root-region arrays")

    all_root_center_coordinates = np.asarray(
        hierarchy.center_coordinates[voxel_regions.root_ball_indices],
        dtype=float,
    )
    if all_root_center_coordinates.ndim != 2 or all_root_center_coordinates.shape[1] != image_ndim:
        raise ValueError("root center coordinates must have shape (n_roots, image_ndim)")

    root_center_coordinates = all_root_center_coordinates[retained_region_labels]
    root_radii_voxels = all_root_radii_voxels[retained_region_labels]

    pore_coords = root_center_coordinates * float(voxel_size)
    if pore_coords.shape[1] == 2:
        pore_coords = np.column_stack([pore_coords, np.zeros(region_count, dtype=float)])
    pore_radius = root_radii_voxels * float(voxel_size)
    pore_area = np.pi * pore_radius**2
    pore_volume = (
        np.asarray(region_adjacency.region_volume_voxels, dtype=float) * float(voxel_size) ** 3
    )
    pore_surface_area = (
        np.asarray(region_adjacency.region_surface_face_counts, dtype=float)
        * float(voxel_size) ** 2
    )

    throat_region_pairs = np.asarray(region_adjacency.throat_region_pairs, dtype=np.int64)
    if throat_region_pairs.size == 0:
        throat_region_pairs = np.zeros((0, 2), dtype=np.int64)
    throat_count = int(throat_region_pairs.shape[0])
    throat_face_counts = np.asarray(region_adjacency.throat_face_counts, dtype=float)
    axis_face_balance = np.asarray(region_adjacency.throat_axis_face_balance, dtype=float)
    if normalized_throat_area_mode == "vector_magnitude":
        throat_cross_area_face_counts = np.linalg.norm(axis_face_balance, axis=1)
        throat_cross_area_face_counts = np.where(
            throat_cross_area_face_counts > 0.0,
            throat_cross_area_face_counts,
            np.maximum(throat_face_counts, 0.1),
        )
    else:
        throat_cross_area_face_counts = throat_face_counts
    throat_area = throat_cross_area_face_counts * float(voxel_size) ** 2
    throat_centroid_indices = np.asarray(region_adjacency.throat_centroid_indices, dtype=float)
    if throat_centroid_indices.shape[0] != throat_count:
        raise ValueError("throat_centroid_indices must align with throat_region_pairs")
    throat_centroid_coords = throat_centroid_indices * float(voxel_size)
    if throat_centroid_coords.shape[1] == 2:
        throat_centroid_coords = np.column_stack(
            [throat_centroid_coords, np.zeros(throat_count, dtype=float)]
        )

    (
        first_side_center_indices,
        first_side_radius_voxels,
        second_side_center_indices,
        second_side_radius_voxels,
    ) = _select_interface_supporting_ball_data(extraction_result)
    pore_radius_by_region = root_radii_voxels
    equivalent_interface_radius = np.sqrt(np.maximum(throat_area, 0.0) / np.pi) / float(voxel_size)
    throat_support_radius_voxels = np.empty(throat_count, dtype=float)
    throat_shape_factor_radius_voxels = np.empty(throat_count, dtype=float)
    for throat_index, region_pair in enumerate(throat_region_pairs):
        first_region_label = int(region_pair[0])
        second_region_label = int(region_pair[1])
        first_radius_voxels = float(first_side_radius_voxels[throat_index])
        second_radius_voxels = float(second_side_radius_voxels[throat_index])
        if not np.isfinite(first_radius_voxels):
            first_radius_voxels = min(
                float(pore_radius_by_region[first_region_label]),
                float(equivalent_interface_radius[throat_index]),
            )
        if not np.isfinite(second_radius_voxels):
            second_radius_voxels = min(
                float(pore_radius_by_region[second_region_label]),
                float(equivalent_interface_radius[throat_index]),
            )
        throat_support_radius_voxels[throat_index] = min(
            0.5 * (first_radius_voxels + second_radius_voxels),
            float(pore_radius_by_region[first_region_label]),
            float(pore_radius_by_region[second_region_label]),
        )
        surface_ball_radius_voxels = min(
            max(second_radius_voxels, 0.5),
            float(pore_radius_by_region[first_region_label]),
            float(pore_radius_by_region[second_region_label]),
        )
        if normalized_throat_shape_factor_radius_mode == "surface_ball":
            throat_shape_factor_radius_voxels[throat_index] = surface_ball_radius_voxels
        else:
            throat_shape_factor_radius_voxels[throat_index] = throat_support_radius_voxels[
                throat_index
            ]
    throat_radius = throat_support_radius_voxels * float(voxel_size)
    throat_shape_factor_radius = throat_shape_factor_radius_voxels * float(voxel_size)
    first_side_radius_physical = np.where(
        np.isfinite(first_side_radius_voxels),
        first_side_radius_voxels * float(voxel_size),
        throat_radius,
    )
    second_side_radius_physical = np.where(
        np.isfinite(second_side_radius_voxels),
        second_side_radius_voxels * float(voxel_size),
        throat_radius,
    )

    minimum_pore_segment_length = 0.5 * float(voxel_size)
    minimum_total_length = 3.01 * float(voxel_size)
    minimum_throat_core_length = 1.0 * float(voxel_size)
    if throat_count:
        first_pore_coordinates = pore_coords[throat_region_pairs[:, 0]]
        second_pore_coordinates = pore_coords[throat_region_pairs[:, 1]]
        if normalized_throat_anchor_mode == "second_side":
            throat_anchor_center_indices = second_side_center_indices.copy()
        else:
            throat_anchor_center_indices = np.where(
                (first_side_radius_voxels >= second_side_radius_voxels)[:, np.newaxis],
                first_side_center_indices,
                second_side_center_indices,
            )
        invalid_anchor_mask = ~np.isfinite(throat_anchor_center_indices).all(axis=1)
        throat_anchor_center_indices[invalid_anchor_mask] = throat_centroid_indices[
            invalid_anchor_mask
        ]
        throat_anchor_coordinates = throat_anchor_center_indices * float(voxel_size)
        if throat_anchor_coordinates.shape[1] == 2:
            throat_anchor_coordinates = np.column_stack(
                [throat_anchor_coordinates, np.zeros(throat_count, dtype=float)]
            )
        pore1_to_anchor_length = np.linalg.norm(
            throat_anchor_coordinates - first_pore_coordinates,
            axis=1,
        )
        pore2_to_anchor_length = np.linalg.norm(
            second_pore_coordinates - throat_anchor_coordinates,
            axis=1,
        )
        throat_total_length = np.maximum(
            pore1_to_anchor_length + pore2_to_anchor_length,
            minimum_total_length,
        )
        pore1_length = np.maximum(0.67 * pore1_to_anchor_length, minimum_pore_segment_length)
        pore2_length = np.maximum(0.67 * pore2_to_anchor_length, minimum_pore_segment_length)
        core_length = np.maximum(
            throat_total_length - pore1_length - pore2_length,
            minimum_throat_core_length,
        )
        throat_total_length = pore1_length + core_length + pore2_length
    else:
        pore1_length = np.zeros(0, dtype=float)
        pore2_length = np.zeros(0, dtype=float)
        core_length = np.zeros(0, dtype=float)
        throat_total_length = np.zeros(0, dtype=float)

    direct_boundary_label_masks: dict[str, tuple[np.ndarray, np.ndarray]] = {}
    connected_boundary_label_masks: dict[str, tuple[np.ndarray, np.ndarray]] = {}
    boundary_touch_radii_voxels = _max_boundary_touch_radii_by_side(
        voxel_regions.label_image,
        retained_region_labels,
        hierarchy.distance_map,
    )
    pore_boundary = np.zeros(region_count, dtype=bool)
    for axis_index, axis_name in enumerate(active_axis_names):
        lower_face_count = np.asarray(
            region_adjacency.boundary_face_counts[:, 2 * axis_index],
            dtype=np.int64,
        )
        upper_face_count = np.asarray(
            region_adjacency.boundary_face_counts[:, 2 * axis_index + 1],
            dtype=np.int64,
        )
        lower_contact = lower_face_count > 0
        upper_contact = upper_face_count > 0
        lower_contact, upper_contact = _resolve_axis_boundary_label_overlap(
            lower_contact,
            upper_contact,
            lower_face_count=lower_face_count,
            upper_face_count=upper_face_count,
            pore_axis_coordinates=pore_coords[:, axis_index],
            sample_axis_length=voxel_regions.label_image.shape[axis_index] * float(voxel_size),
        )
        direct_boundary_label_masks[axis_name] = (lower_contact, upper_contact)
        connected_boundary_label_masks[axis_name] = (lower_contact.copy(), upper_contact.copy())
        pore_boundary |= lower_contact | upper_contact

    if normalized_flow_boundary_mode == "external_reservoir" and boundary_axis is not None:
        boundary_axis_index = active_axis_names.index(boundary_axis)
        lower_contact, upper_contact = direct_boundary_label_masks[boundary_axis]
        lower_face_count = np.asarray(
            region_adjacency.boundary_face_counts[:, 2 * boundary_axis_index],
            dtype=np.int64,
        )
        upper_face_count = np.asarray(
            region_adjacency.boundary_face_counts[:, 2 * boundary_axis_index + 1],
            dtype=np.int64,
        )
        sample_axis_length = voxel_regions.label_image.shape[boundary_axis_index] * float(
            voxel_size
        )

        helper_coordinates: list[np.ndarray] = []
        helper_radii: list[float] = []
        helper_areas: list[float] = []
        helper_volumes: list[float] = []
        helper_surface_areas: list[float] = []
        boundary_throat_connections: list[tuple[int, int]] = []
        boundary_throat_area: list[float] = []
        boundary_throat_radius: list[float] = []
        boundary_throat_shape_factor_radius: list[float] = []
        boundary_throat_pore1_length: list[float] = []
        boundary_throat_core_length: list[float] = []
        boundary_throat_pore2_length: list[float] = []
        boundary_throat_centroids: list[np.ndarray] = []
        boundary_throat_face_counts: list[float] = []
        boundary_axis_face_balance: list[np.ndarray] = []
        boundary_support_radius_side1: list[float] = []
        boundary_support_radius_side2: list[float] = []
        inlet_helper_mask: list[bool] = []
        outlet_helper_mask: list[bool] = []

        def append_boundary_helper_pores(
            contact_mask: np.ndarray,
            face_counts: np.ndarray,
            side: str,
        ) -> None:
            boundary_coordinate = 0.0 if side == "lower" else sample_axis_length
            contact_indices = np.flatnonzero(contact_mask)
            for pore_index in contact_indices:
                helper_index = region_count + len(helper_coordinates)
                physical_coordinate = pore_coords[pore_index]
                helper_coordinate = physical_coordinate.copy()
                helper_coordinate[boundary_axis_index] = boundary_coordinate
                face_count = max(int(face_counts[pore_index]), 1)
                contact_area = float(face_count) * float(voxel_size) ** 2
                boundary_touch_radius_voxels = boundary_touch_radii_voxels[
                    pore_index,
                    2 * boundary_axis_index + (0 if side == "lower" else 1),
                ]
                if np.isfinite(boundary_touch_radius_voxels):
                    contact_radius = min(
                        float(pore_radius[pore_index]),
                        float(boundary_touch_radius_voxels) * float(voxel_size),
                    )
                else:
                    contact_radius = min(
                        float(pore_radius[pore_index]),
                        float(np.sqrt(contact_area / np.pi)),
                    )
                contact_radius = max(contact_radius, 0.5 * float(voxel_size))
                helper_radius = boundary_radius_scale * contact_radius
                center_to_boundary_length = abs(
                    float(physical_coordinate[boundary_axis_index] - boundary_coordinate)
                )
                total_boundary_length = max(
                    center_to_boundary_length,
                    3.01 * float(voxel_size),
                )
                internal_pore_length = max(0.67 * center_to_boundary_length, 0.0)
                boundary_core_length = max(
                    total_boundary_length - boundary_length_epsilon - internal_pore_length,
                    float(voxel_size),
                )

                helper_coordinates.append(helper_coordinate)
                helper_radii.append(helper_radius)
                helper_areas.append(np.pi * helper_radius**2)
                helper_volumes.append(0.0)
                helper_surface_areas.append(0.0)
                boundary_throat_area.append(contact_area)
                boundary_throat_radius.append(contact_radius)
                boundary_throat_shape_factor_radius.append(contact_radius)
                boundary_throat_centroids.append(helper_coordinate.copy())
                boundary_throat_face_counts.append(float(face_count))
                axis_balance = np.zeros(image_ndim, dtype=float)
                axis_balance[boundary_axis_index] = -face_count if side == "lower" else face_count
                boundary_axis_face_balance.append(axis_balance)
                if side == "lower":
                    boundary_throat_connections.append((helper_index, int(pore_index)))
                    boundary_throat_pore1_length.append(boundary_length_epsilon)
                    boundary_throat_pore2_length.append(internal_pore_length)
                    inlet_helper_mask.append(True)
                    outlet_helper_mask.append(False)
                else:
                    boundary_throat_connections.append((int(pore_index), helper_index))
                    boundary_throat_pore1_length.append(internal_pore_length)
                    boundary_throat_pore2_length.append(boundary_length_epsilon)
                    inlet_helper_mask.append(False)
                    outlet_helper_mask.append(True)
                boundary_throat_core_length.append(boundary_core_length)
                boundary_support_radius_side1.append(contact_radius)
                boundary_support_radius_side2.append(contact_radius)

        append_boundary_helper_pores(lower_contact, lower_face_count, "lower")
        append_boundary_helper_pores(upper_contact, upper_face_count, "upper")

        helper_count = len(helper_coordinates)
        if helper_count:
            pore_coords = np.vstack([pore_coords, np.asarray(helper_coordinates, dtype=float)])
            pore_radius = np.concatenate([pore_radius, np.asarray(helper_radii, dtype=float)])
            pore_area = np.concatenate([pore_area, np.asarray(helper_areas, dtype=float)])
            pore_volume = np.concatenate([pore_volume, np.asarray(helper_volumes, dtype=float)])
            pore_surface_area = np.concatenate(
                [pore_surface_area, np.asarray(helper_surface_areas, dtype=float)]
            )
            throat_region_pairs = np.vstack(
                [
                    throat_region_pairs,
                    np.asarray(boundary_throat_connections, dtype=np.int64),
                ]
            )
            throat_area = np.concatenate(
                [throat_area, np.asarray(boundary_throat_area, dtype=float)]
            )
            throat_radius = np.concatenate(
                [throat_radius, np.asarray(boundary_throat_radius, dtype=float)]
            )
            throat_shape_factor_radius = np.concatenate(
                [
                    throat_shape_factor_radius,
                    np.asarray(boundary_throat_shape_factor_radius, dtype=float),
                ]
            )
            pore1_length = np.concatenate(
                [pore1_length, np.asarray(boundary_throat_pore1_length, dtype=float)]
            )
            core_length = np.concatenate(
                [core_length, np.asarray(boundary_throat_core_length, dtype=float)]
            )
            pore2_length = np.concatenate(
                [pore2_length, np.asarray(boundary_throat_pore2_length, dtype=float)]
            )
            added_total_lengths = (
                pore1_length[-helper_count:]
                + core_length[-helper_count:]
                + pore2_length[-helper_count:]
            )
            throat_total_length = np.concatenate([throat_total_length, added_total_lengths])
            throat_centroid_coords = np.vstack(
                [throat_centroid_coords, np.asarray(boundary_throat_centroids, dtype=float)]
            )
            throat_face_counts = np.concatenate(
                [throat_face_counts, np.asarray(boundary_throat_face_counts, dtype=float)]
            )
            axis_face_balance = np.vstack(
                [
                    axis_face_balance,
                    np.asarray(boundary_axis_face_balance, dtype=float),
                ]
            )
            first_side_radius_physical = np.concatenate(
                [
                    first_side_radius_physical,
                    np.asarray(boundary_support_radius_side1, dtype=float),
                ]
            )
            second_side_radius_physical = np.concatenate(
                [
                    second_side_radius_physical,
                    np.asarray(boundary_support_radius_side2, dtype=float),
                ]
            )
            pore_boundary = np.concatenate([pore_boundary, np.ones(helper_count, dtype=bool)])
            for axis_name in active_axis_names:
                lower_mask, upper_mask = direct_boundary_label_masks[axis_name]
                direct_boundary_label_masks[axis_name] = (
                    np.concatenate([lower_mask, np.zeros(helper_count, dtype=bool)]),
                    np.concatenate([upper_mask, np.zeros(helper_count, dtype=bool)]),
                )
                connected_lower, connected_upper = connected_boundary_label_masks[axis_name]
                connected_boundary_label_masks[axis_name] = (
                    np.concatenate([connected_lower, np.zeros(helper_count, dtype=bool)]),
                    np.concatenate([connected_upper, np.zeros(helper_count, dtype=bool)]),
                )
            inlet_label = np.concatenate(
                [
                    np.zeros(region_count, dtype=bool),
                    np.asarray(inlet_helper_mask, dtype=bool),
                ]
            )
            outlet_label = np.concatenate(
                [
                    np.zeros(region_count, dtype=bool),
                    np.asarray(outlet_helper_mask, dtype=bool),
                ]
            )
            direct_boundary_label_masks[boundary_axis] = (inlet_label, outlet_label)
            region_count = int(pore_coords.shape[0])
            throat_count = int(throat_region_pairs.shape[0])

    pore_data: dict[str, np.ndarray] = {
        "radius_inscribed": pore_radius.copy(),
        "area": pore_area.copy(),
        "shape_factor": np.full(region_count, _CIRCULAR_SHAPE_FACTOR, dtype=float),
        "volume": pore_volume.copy(),
        "region_volume": pore_volume.copy(),
        "surface_area": pore_surface_area,
    }
    throat_data: dict[str, np.ndarray] = {
        "radius_inscribed": throat_radius.copy(),
        "shape_factor_radius": throat_shape_factor_radius.copy(),
        "area": throat_area.copy(),
        "shape_factor": np.full(throat_count, _CIRCULAR_SHAPE_FACTOR, dtype=float),
        "volume": np.zeros(throat_count, dtype=float),
        "total_length": throat_total_length,
        "pore1_length": pore1_length,
        "core_length": core_length,
        "pore2_length": pore2_length,
        "centroid": throat_centroid_coords,
        "face_count": throat_face_counts.astype(np.int64, copy=False),
        "axis_face_balance": axis_face_balance,
        "supporting_radius_side1": first_side_radius_physical,
        "supporting_radius_side2": second_side_radius_physical,
    }

    from voids.io.porespy import _apply_imperial_export_geometry_repairs, _derive_missing_geometry

    geometry_repair_summary = _apply_imperial_export_geometry_repairs(
        pore_data,
        throat_data,
        throat_region_pairs,
        num_pores=region_count,
        random_seed=1001,
    )
    _redistribute_region_volumes_like_imperial_export(
        pore_data,
        throat_data,
        throat_region_pairs,
    )
    _derive_missing_geometry(pore_data, throat_data)
    boundary_face_count = np.asarray(
        region_adjacency.boundary_face_counts.sum(axis=1),
        dtype=np.int64,
    )
    if boundary_face_count.size < region_count:
        boundary_face_count = np.concatenate(
            [
                boundary_face_count,
                np.zeros(region_count - boundary_face_count.size, dtype=np.int64),
            ]
        )

    network_dict: dict[str, np.ndarray] = {
        "pore.coords": pore_coords,
        "throat.conns": throat_region_pairs,
        "pore.radius_inscribed": pore_data["radius_inscribed"],
        "pore.area": pore_data["area"],
        "pore.shape_factor": pore_data["shape_factor"],
        "pore.volume": pore_data["volume"],
        "pore.region_volume": pore_data["region_volume"],
        "pore.surface_area": pore_data["surface_area"],
        "throat.radius_inscribed": throat_data["radius_inscribed"],
        "throat.shape_factor_radius": throat_data["shape_factor_radius"],
        "throat.cross_sectional_area": throat_data["area"],
        "throat.shape_factor": throat_data["shape_factor"],
        "throat.volume": throat_data["volume"],
        "throat.total_length": throat_data["length"],
        "throat.conduit_lengths.pore1": throat_data["pore1_length"],
        "throat.conduit_lengths.throat": throat_data["core_length"],
        "throat.conduit_lengths.pore2": throat_data["pore2_length"],
        "throat.centroid": throat_data["centroid"],
        "throat.face_count": throat_data["face_count"],
        "throat.axis_face_balance": throat_data["axis_face_balance"],
        "throat.supporting_radius_side1": throat_data["supporting_radius_side1"],
        "throat.supporting_radius_side2": throat_data["supporting_radius_side2"],
        "pore.boundary_face_count": boundary_face_count,
        "throat.geometry_repairs_mode": np.full(
            throat_count, geometry_repair_summary["mode"], dtype=object
        ),
    }

    for axis_name in active_axis_names:
        lower_contact, upper_contact = direct_boundary_label_masks[axis_name]
        network_dict[f"pore.inlet_{axis_name}min"] = lower_contact
        network_dict[f"pore.outlet_{axis_name}max"] = upper_contact
        connected_lower, connected_upper = connected_boundary_label_masks[axis_name]
        network_dict[f"pore.boundary_connected_inlet_{axis_name}min"] = connected_lower
        network_dict[f"pore.boundary_connected_outlet_{axis_name}max"] = connected_upper
    network_dict["pore.boundary"] = pore_boundary
    return network_dict

extract_maximal_ball_network_dict

extract_maximal_ball_network_dict(
    void_phase_mask,
    *,
    voxel_size,
    distance_map_backend="auto",
    edt_parallel_threads=None,
    settings=None,
    apply_boundary_clipping=True,
    axis_names=("x", "y", "z"),
    flow_boundary_mode="direct",
    boundary_axis=None,
    boundary_length_epsilon=1e-300,
    boundary_radius_scale=1.1,
    throat_area_mode="face_count",
    throat_shape_factor_radius_mode="inscribed",
    throat_anchor_mode="second_side",
)

Run the staged native maximal-ball path and assemble a network mapping.

Source code in src/voids/image/maximal_ball.py

def extract_maximal_ball_network_dict(
    void_phase_mask: np.ndarray,
    *,
    voxel_size: float,
    distance_map_backend: str = "auto",
    edt_parallel_threads: int | None = None,
    settings: MaximalBallSettings | None = None,
    apply_boundary_clipping: bool = True,
    axis_names: tuple[str, ...] = ("x", "y", "z"),
    flow_boundary_mode: str = "direct",
    boundary_axis: str | None = None,
    boundary_length_epsilon: float = 1.0e-300,
    boundary_radius_scale: float = 1.1,
    throat_area_mode: str = "face_count",
    throat_shape_factor_radius_mode: str = "inscribed",
    throat_anchor_mode: str = "second_side",
) -> MaximalBallNetworkDictResult:
    """Run the staged native maximal-ball path and assemble a network mapping."""

    extraction_result = extract_maximal_ball_regions(
        void_phase_mask,
        distance_map_backend=distance_map_backend,
        edt_parallel_threads=edt_parallel_threads,
        settings=settings,
        apply_boundary_clipping=apply_boundary_clipping,
    )
    network_dict = build_network_dict_from_maximal_ball_regions(
        extraction_result,
        voxel_size=voxel_size,
        axis_names=axis_names,
        flow_boundary_mode=flow_boundary_mode,
        boundary_axis=boundary_axis,
        boundary_length_epsilon=boundary_length_epsilon,
        boundary_radius_scale=boundary_radius_scale,
        throat_area_mode=throat_area_mode,
        throat_shape_factor_radius_mode=throat_shape_factor_radius_mode,
        throat_anchor_mode=throat_anchor_mode,
    )
    return MaximalBallNetworkDictResult(
        network_dict=network_dict,
        extraction=extraction_result,
    )

extract_maximal_ball_candidates

extract_maximal_ball_candidates(
    void_phase_mask,
    *,
    distance_map_backend="auto",
    edt_parallel_threads=None,
    settings=None,
    apply_boundary_clipping=True,
)

Compute and suppress maximal-ball candidates for a void-phase image.

Source code in src/voids/image/maximal_ball.py

def extract_maximal_ball_candidates(
    void_phase_mask: np.ndarray,
    *,
    distance_map_backend: str = "auto",
    edt_parallel_threads: int | None = None,
    settings: MaximalBallSettings | None = None,
    apply_boundary_clipping: bool = True,
) -> MaximalBallCandidates:
    """Compute and suppress maximal-ball candidates for a void-phase image."""

    raw_settings = settings or MaximalBallSettings()
    distance_map = compute_maximal_ball_radius_field(
        void_phase_mask,
        backend=distance_map_backend,
        edt_parallel_threads=edt_parallel_threads,
        mode=raw_settings.radius_field_mode,
    )
    resolved_settings = resolve_maximal_ball_settings(distance_map, raw_settings)
    working_distance_map = (
        clip_distance_map_to_domain_boundaries(distance_map, settings=resolved_settings)
        if apply_boundary_clipping
        else np.asarray(distance_map, dtype=float)
    )
    center_indices, radii_voxels, candidate_mask = find_maximal_ball_candidates(
        working_distance_map,
        minimal_radius_voxels=resolved_settings.minimal_pore_radius_voxels,
        selection_mode=resolved_settings.candidate_selection_mode,
    )
    retained_mask = suppress_overlapping_maximal_balls(
        center_indices,
        radii_voxels,
        settings=resolved_settings,
    )
    return MaximalBallCandidates(
        center_indices=center_indices,
        radii_voxels=radii_voxels,
        candidate_mask=candidate_mask,
        retained_mask=retained_mask,
        distance_map=working_distance_map,
        settings=resolved_settings,
    )

---

PREGO Extraction

voids.image.prego

PregoSettings dataclass

Controls for seed-based Pore Region Growing segmentation.

The defaults use r_max=4 and Gaussian smoothing sigma=0.4. PREGO is currently implemented for a single active pore phase, with nonzero input voxels treated as void.

Source code in src/voids/image/prego.py

@dataclass(slots=True)
class PregoSettings:
    """Controls for seed-based Pore Region Growing segmentation.

    The defaults use ``r_max=4`` and Gaussian smoothing ``sigma=0.4``. PREGO is
    currently implemented for a single active pore phase, with nonzero input
    voxels treated as void.
    """

    r_max: int = 4
    sigma: float = 0.4
    peak_footprint: str = "sphere"
    distance_map_backend: str = "auto"
    edt_parallel_threads: int | None = None
    cleanup_unassigned: bool = True
    growth_mode: str = "level_queue"

PregoSegmentationResult dataclass

Intermediate PREGO segmentation data before network construction.

Source code in src/voids/image/prego.py

@dataclass(slots=True)
class PregoSegmentationResult:
    """Intermediate PREGO segmentation data before network construction."""

    im: np.ndarray
    distance_map: np.ndarray
    peaks: np.ndarray
    regions: np.ndarray
    seed_indices: np.ndarray
    seed_radii_voxels: np.ndarray
    seed_activation_levels: np.ndarray
    settings: PregoSettings

PregoNetworkDictResult dataclass

PoreSpy-style network mapping assembled from PREGO regions.

Source code in src/voids/image/prego.py

@dataclass(slots=True)
class PregoNetworkDictResult:
    """PoreSpy-style network mapping assembled from PREGO regions."""

    network_dict: dict[str, object]
    segmentation: PregoSegmentationResult

snow_seed_points

snow_seed_points(
    void_phase_mask,
    *,
    distance_map=None,
    r_max=4,
    sigma=0.4,
    peak_footprint="sphere",
    peaks=None,
    distance_map_backend="auto",
    edt_parallel_threads=None,
    porespy_module=ps,
)

Find PREGO seed points using the peak-filtering stages of SNOW.

Returns:

Type	Description
tuple	(seed_indices, seed_labels, distance_map) where seed_labels has one labeled voxel per seed and labels are ordered by descending seed radius.

Source code in src/voids/image/prego.py

def snow_seed_points(
    void_phase_mask: np.ndarray,
    *,
    distance_map: np.ndarray | None = None,
    r_max: int = 4,
    sigma: float = 0.4,
    peak_footprint: str = "sphere",
    peaks: np.ndarray | None = None,
    distance_map_backend: str = "auto",
    edt_parallel_threads: int | None = None,
    porespy_module: Any = ps,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Find PREGO seed points using the peak-filtering stages of SNOW.

    Returns
    -------
    tuple
        ``(seed_indices, seed_labels, distance_map)`` where ``seed_labels`` has
        one labeled voxel per seed and labels are ordered by descending seed
        radius.
    """

    mask = np.asarray(void_phase_mask, dtype=bool)
    if mask.ndim not in {2, 3}:
        raise ValueError("void_phase_mask must be a 2D or 3D array")
    if distance_map is None:
        dt = compute_void_distance_map(
            mask,
            backend=distance_map_backend,
            edt_parallel_threads=edt_parallel_threads,
        )
    else:
        dt = np.asarray(distance_map, dtype=float)
        if dt.shape != mask.shape:
            raise ValueError("distance_map must match void_phase_mask")

    if peaks is None:
        if sigma > 0:
            peak_distance_map = ndi.gaussian_filter(input=dt, sigma=float(sigma)) * mask
        else:
            peak_distance_map = dt.copy()
        normalized_footprint = str(peak_footprint).strip().lower()
        if normalized_footprint == "sphere":
            peak_mask = porespy_module.filters.find_peaks(dt=peak_distance_map, r_max=int(r_max))
        elif normalized_footprint in {"cube", "box"}:
            peak_distance_for_filter = peak_distance_map + 2.0 * (~mask)
            peak_max = ndi.maximum_filter(
                peak_distance_for_filter,
                size=(2 * int(r_max) + 1,) * mask.ndim,
            )
            peak_mask = (peak_distance_map == peak_max) & mask & (peak_distance_map > 0)
        else:
            raise ValueError("peak_footprint must be 'cube' or 'sphere'")
        peak_mask = porespy_module.filters.trim_saddle_points(peaks=peak_mask, dt=dt)
        peak_mask = porespy_module.filters.trim_nearby_peaks(peaks=peak_mask, dt=dt)
        peak_labels, _ = ndi.label(peak_mask > 0, structure=_connectivity_structure(mask.ndim))
    else:
        supplied_peaks = np.asarray(peaks)
        if supplied_peaks.shape != mask.shape:
            raise ValueError("peaks must match void_phase_mask")
        if supplied_peaks.dtype == bool:
            peak_labels, _ = ndi.label(
                supplied_peaks & mask,
                structure=_connectivity_structure(mask.ndim),
            )
        else:
            supplied_labels = np.asarray(supplied_peaks)
            max_label = int(np.max(supplied_labels)) if supplied_labels.size else 0
            peak_dtype = _prego_label_dtype(max_label=max_label, shape=mask.shape)
            peak_labels = supplied_labels.astype(peak_dtype, copy=False) * mask

    seed_indices, seed_labels = _reduce_peak_labels_to_seed_points(peak_labels, dt)
    if seed_indices.size == 0 and np.any(mask):
        fallback_index = np.asarray(
            np.unravel_index(int(np.argmax(dt)), dt.shape),
            dtype=np.int64,
        )
        seed_indices = fallback_index.reshape(1, mask.ndim)
        seed_labels = np.zeros(mask.shape, dtype=np.int64)
        seed_labels[tuple(int(value) for value in fallback_index)] = 1
    return seed_indices, seed_labels, dt

prego_partitioning

prego_partitioning(
    im,
    *,
    settings=None,
    distance_map=None,
    peaks=None,
    porespy_module=ps,
)

Partition a binary pore image with PREGO-style region growing.

Notes

The default growth_mode="level_queue" follows the paper's delayed seed activation and level-by-level FIFO queue before the final expansion of unassigned foreground voxels. growth_mode="fast" remains available as a faster approximation that stamps non-overlapping seed spheres in descending radius order before the same final FIFO fill.

Source code in src/voids/image/prego.py

def prego_partitioning(
    im: np.ndarray,
    *,
    settings: PregoSettings | None = None,
    distance_map: np.ndarray | None = None,
    peaks: np.ndarray | None = None,
    porespy_module: Any = ps,
) -> PregoSegmentationResult:
    """Partition a binary pore image with PREGO-style region growing.

    Notes
    -----
    The default ``growth_mode="level_queue"`` follows the paper's delayed seed
    activation and level-by-level FIFO queue before the final expansion of
    unassigned foreground voxels. ``growth_mode="fast"`` remains available as
    a faster approximation that stamps non-overlapping seed spheres in
    descending radius order before the same final FIFO fill.
    """

    resolved_settings = settings or PregoSettings()
    mask = np.asarray(im, dtype=bool)
    if mask.ndim not in {2, 3}:
        raise ValueError("im must be a 2D or 3D binary image")
    growth_mode = _normalize_growth_mode(resolved_settings.growth_mode)

    seed_indices, seed_label_map, dt = snow_seed_points(
        mask,
        distance_map=distance_map,
        r_max=resolved_settings.r_max,
        sigma=resolved_settings.sigma,
        peak_footprint=resolved_settings.peak_footprint,
        peaks=peaks,
        distance_map_backend=resolved_settings.distance_map_backend,
        edt_parallel_threads=resolved_settings.edt_parallel_threads,
        porespy_module=porespy_module,
    )
    seed_radii = (
        dt[tuple(seed_indices.T)].astype(float, copy=False)
        if seed_indices.size
        else np.zeros(0, dtype=float)
    )
    activation_levels = _seed_activation_levels(seed_radii)
    labels = np.zeros(
        mask.shape,
        dtype=_prego_label_dtype(max_label=seed_indices.shape[0], shape=mask.shape),
    )
    if seed_indices.size:
        for label, seed_index in enumerate(seed_indices, start=1):
            labels[tuple(int(value) for value in seed_index)] = label
        if growth_mode == "fast":
            if mask.ndim == 2:
                _stamp_seed_spheres_2d(
                    mask, labels, seed_label_map, seed_indices, seed_radii, 1e-12
                )
            else:
                _stamp_seed_spheres_3d(
                    mask, labels, seed_label_map, seed_indices, seed_radii, 1e-12
                )
        elif mask.ndim == 2:
            _level_queue_grow_regions_2d(
                mask,
                labels,
                seed_label_map,
                seed_indices,
                seed_radii,
                activation_levels,
                1e-12,
            )
        else:
            _level_queue_grow_regions_3d(
                mask,
                labels,
                seed_label_map,
                seed_indices,
                seed_radii,
                activation_levels,
                1e-12,
            )
        if mask.ndim == 2:
            _fifo_fill_regions_2d(mask, labels, seed_label_map)
        else:
            _fifo_fill_regions_3d(mask, labels, seed_label_map)

    if resolved_settings.cleanup_unassigned:
        labels = labels * mask
    return PregoSegmentationResult(
        im=mask,
        distance_map=dt,
        peaks=seed_label_map,
        regions=labels,
        seed_indices=seed_indices,
        seed_radii_voxels=seed_radii,
        seed_activation_levels=activation_levels,
        settings=resolved_settings,
    )

extract_prego_network_dict

extract_prego_network_dict(
    im,
    *,
    settings=None,
    distance_map=None,
    peaks=None,
    porespy_module=ps,
    regions_to_network_kwargs=None,
)

Run PREGO segmentation and convert regions to a PoreSpy network dict.

Source code in src/voids/image/prego.py

def extract_prego_network_dict(
    im: np.ndarray,
    *,
    settings: PregoSettings | None = None,
    distance_map: np.ndarray | None = None,
    peaks: np.ndarray | None = None,
    porespy_module: Any = ps,
    regions_to_network_kwargs: dict[str, object] | None = None,
) -> PregoNetworkDictResult:
    """Run PREGO segmentation and convert regions to a PoreSpy network dict."""

    segmentation = prego_partitioning(
        im,
        settings=settings,
        distance_map=distance_map,
        peaks=peaks,
        porespy_module=porespy_module,
    )
    kwargs = dict(regions_to_network_kwargs or {})
    if _regions_have_interfaces(segmentation.regions):
        network_dict = dict(
            porespy_module.networks.regions_to_network(segmentation.regions, **kwargs)
        )
    else:
        network_dict = _network_dict_without_interfaces(
            segmentation.regions,
            segmentation.distance_map,
        )
    return PregoNetworkDictResult(network_dict=network_dict, segmentation=segmentation)

---

Porosity Maps

For conceptual background, block_shape interpretation, and synthetic verification context, see Porosity Maps.

voids.image.porosity

PorosityMap dataclass

Cell-wise porosity field for continuum or external-solver workflows.

Parameters:

Name	Type	Description	Default
values	ndarray	Two- or three-dimensional porosity field with entries in [0, 1]. Values are interpreted as cell averages.	required
cell_size	float | Sequence[float]	Scalar or per-axis cell side length in physical units.	1.0
origin	Sequence[float] | None	Physical coordinate of the lower corner of the first cell.	None
units	dict[str, str]	Unit metadata for reporting and serialization.	(lambda: {'length': 'm'})()
metadata	dict[str, Any]	Additional JSON-serializable provenance and calibration metadata.	dict()

Source code in src/voids/image/porosity.py

@dataclass(slots=True)
class PorosityMap:
    """Cell-wise porosity field for continuum or external-solver workflows.

    Parameters
    ----------
    values :
        Two- or three-dimensional porosity field with entries in ``[0, 1]``.
        Values are interpreted as cell averages.
    cell_size :
        Scalar or per-axis cell side length in physical units.
    origin :
        Physical coordinate of the lower corner of the first cell.
    units :
        Unit metadata for reporting and serialization.
    metadata :
        Additional JSON-serializable provenance and calibration metadata.
    """

    values: np.ndarray
    cell_size: float | Sequence[float] = 1.0
    origin: Sequence[float] | None = None
    units: dict[str, str] = field(default_factory=lambda: {"length": "m"})
    metadata: dict[str, Any] = field(default_factory=dict)

    def __post_init__(self) -> None:
        """Normalize and validate porosity-map fields."""

        values = _validate_porosity_values(self.values)
        self.values = values
        self.cell_size = _normalize_positive_tuple(
            self.cell_size,
            ndim=values.ndim,
            name="cell_size",
        )
        self.origin = _normalize_origin(self.origin, ndim=values.ndim)
        self.units = {str(k): str(v) for k, v in self.units.items()}
        self.metadata = dict(self.metadata)

    @property
    def ndim(self) -> int:
        """Return the spatial dimensionality of the porosity map."""

        return int(self.values.ndim)

    @property
    def shape(self) -> tuple[int, ...]:
        """Return the cell-count shape of the porosity map."""

        return tuple(int(v) for v in self.values.shape)

    @property
    def mean_porosity(self) -> float:
        """Return the arithmetic mean porosity over all cells."""

        return float(np.mean(self.values))

    @property
    def cell_volume(self) -> float:
        """Return the physical volume represented by one porosity-map cell."""

        return float(np.prod(np.asarray(self.cell_size, dtype=float)))

    @property
    def bulk_volume(self) -> float:
        """Return the physical bulk volume represented by the map."""

        return float(np.prod(np.asarray(self.shape, dtype=float)) * self.cell_volume)

    @property
    def void_volume(self) -> float:
        """Return the pore volume implied by the cell-average porosity field."""

        return float(np.sum(self.values) * self.cell_volume)

    def to_metadata(self) -> dict[str, Any]:
        """Serialize map metadata without the porosity array."""

        return {
            "schema_version": _SCHEMA_VERSION,
            "shape": self.shape,
            "cell_size": self.cell_size,
            "origin": self.origin,
            "units": self.units,
            "metadata": self.metadata,
        }

    @classmethod
    def from_metadata(cls, values: np.ndarray, metadata: dict[str, Any]) -> "PorosityMap":
        """Reconstruct a porosity map from values and serialized metadata."""

        return cls(
            values=values,
            cell_size=metadata.get("cell_size", 1.0),
            origin=metadata.get("origin"),
            units={str(k): str(v) for k, v in (metadata.get("units") or {}).items()},
            metadata=dict(metadata.get("metadata") or {}),
        )

ndim property

ndim

Return the spatial dimensionality of the porosity map.

shape property

shape

Return the cell-count shape of the porosity map.

mean_porosity property

mean_porosity

Return the arithmetic mean porosity over all cells.

cell_volume property

cell_volume

Return the physical volume represented by one porosity-map cell.

bulk_volume property

bulk_volume

Return the physical bulk volume represented by the map.

void_volume property

void_volume

Return the pore volume implied by the cell-average porosity field.

to_metadata

to_metadata()

Serialize map metadata without the porosity array.

Source code in src/voids/image/porosity.py

def to_metadata(self) -> dict[str, Any]:
    """Serialize map metadata without the porosity array."""

    return {
        "schema_version": _SCHEMA_VERSION,
        "shape": self.shape,
        "cell_size": self.cell_size,
        "origin": self.origin,
        "units": self.units,
        "metadata": self.metadata,
    }

from_metadata classmethod

from_metadata(values, metadata)

Reconstruct a porosity map from values and serialized metadata.

Source code in src/voids/image/porosity.py

@classmethod
def from_metadata(cls, values: np.ndarray, metadata: dict[str, Any]) -> "PorosityMap":
    """Reconstruct a porosity map from values and serialized metadata."""

    return cls(
        values=values,
        cell_size=metadata.get("cell_size", 1.0),
        origin=metadata.get("origin"),
        units={str(k): str(v) for k, v in (metadata.get("units") or {}).items()},
        metadata=dict(metadata.get("metadata") or {}),
    )

PermeabilityMap dataclass

Cell-wise permeability field paired with a continuum porosity map.

Parameters:

Name	Type	Description	Default
values	ndarray	Two- or three-dimensional permeability field. Values are interpreted in the square of the length unit recorded in units.	required
cell_size	float | Sequence[float]	Scalar or per-axis cell side length in physical units.	1.0
origin	Sequence[float] | None	Physical coordinate of the lower corner of the first cell.	None
units	dict[str, str]	Unit metadata for reporting and serialization.	(lambda: {'length': 'm', 'permeability': 'm^2'})()
metadata	dict[str, Any]	Additional JSON-serializable provenance and closure metadata.	dict()

Source code in src/voids/image/porosity.py

@dataclass(slots=True)
class PermeabilityMap:
    """Cell-wise permeability field paired with a continuum porosity map.

    Parameters
    ----------
    values :
        Two- or three-dimensional permeability field. Values are interpreted in
        the square of the length unit recorded in ``units``.
    cell_size :
        Scalar or per-axis cell side length in physical units.
    origin :
        Physical coordinate of the lower corner of the first cell.
    units :
        Unit metadata for reporting and serialization.
    metadata :
        Additional JSON-serializable provenance and closure metadata.
    """

    values: np.ndarray
    cell_size: float | Sequence[float] = 1.0
    origin: Sequence[float] | None = None
    units: dict[str, str] = field(default_factory=lambda: {"length": "m", "permeability": "m^2"})
    metadata: dict[str, Any] = field(default_factory=dict)

    def __post_init__(self) -> None:
        """Normalize and validate permeability-map fields."""

        values = _validate_permeability_values(self.values)
        self.values = values
        self.cell_size = _normalize_positive_tuple(
            self.cell_size,
            ndim=values.ndim,
            name="cell_size",
        )
        self.origin = _normalize_origin(self.origin, ndim=values.ndim)
        self.units = {str(k): str(v) for k, v in self.units.items()}
        self.metadata = dict(self.metadata)

    @property
    def ndim(self) -> int:
        """Return the spatial dimensionality of the permeability map."""

        return int(self.values.ndim)

    @property
    def shape(self) -> tuple[int, ...]:
        """Return the cell-count shape of the permeability map."""

        return tuple(int(v) for v in self.values.shape)

    @property
    def finite_mean_permeability(self) -> float:
        """Return the arithmetic mean over finite permeability values."""

        finite = self.values[np.isfinite(self.values)]
        if finite.size == 0:
            return float("nan")
        return float(np.mean(finite))

    @property
    def inverse_values(self) -> np.ndarray:
        """Return the inverse permeability field with reciprocal endpoint limits."""

        with np.errstate(divide="ignore", invalid="ignore"):
            return np.asarray(np.reciprocal(self.values), dtype=float)

    def to_metadata(self) -> dict[str, Any]:
        """Serialize map metadata without the permeability array."""

        return {
            "schema_version": _PERMEABILITY_SCHEMA_VERSION,
            "shape": self.shape,
            "cell_size": self.cell_size,
            "origin": self.origin,
            "units": self.units,
            "metadata": self.metadata,
        }

    @classmethod
    def from_metadata(cls, values: np.ndarray, metadata: dict[str, Any]) -> "PermeabilityMap":
        """Reconstruct a permeability map from values and serialized metadata."""

        return cls(
            values=values,
            cell_size=metadata.get("cell_size", 1.0),
            origin=metadata.get("origin"),
            units={str(k): str(v) for k, v in (metadata.get("units") or {}).items()},
            metadata=dict(metadata.get("metadata") or {}),
        )

ndim property

ndim

Return the spatial dimensionality of the permeability map.

shape property

shape

Return the cell-count shape of the permeability map.

finite_mean_permeability property

finite_mean_permeability

Return the arithmetic mean over finite permeability values.

inverse_values property

inverse_values

Return the inverse permeability field with reciprocal endpoint limits.

to_metadata

to_metadata()

Serialize map metadata without the permeability array.

Source code in src/voids/image/porosity.py

def to_metadata(self) -> dict[str, Any]:
    """Serialize map metadata without the permeability array."""

    return {
        "schema_version": _PERMEABILITY_SCHEMA_VERSION,
        "shape": self.shape,
        "cell_size": self.cell_size,
        "origin": self.origin,
        "units": self.units,
        "metadata": self.metadata,
    }

from_metadata classmethod

from_metadata(values, metadata)

Reconstruct a permeability map from values and serialized metadata.

Source code in src/voids/image/porosity.py

@classmethod
def from_metadata(cls, values: np.ndarray, metadata: dict[str, Any]) -> "PermeabilityMap":
    """Reconstruct a permeability map from values and serialized metadata."""

    return cls(
        values=values,
        cell_size=metadata.get("cell_size", 1.0),
        origin=metadata.get("origin"),
        units={str(k): str(v) for k, v in (metadata.get("units") or {}).items()},
        metadata=dict(metadata.get("metadata") or {}),
    )

calibrated_porosity_from_grayscale

calibrated_porosity_from_grayscale(
    grayscale,
    *,
    solid_gray,
    pore_gray,
    background_porosity=0.0,
    clip=True,
)

Map grayscale intensity to porosity with a two-point linear calibration.

Parameters:

Name	Type	Description	Default
grayscale	ndarray	Two- or three-dimensional grayscale image.	required
solid_gray	float	Grayscale value assigned to the background or unresolved-microporosity porosity.	required
pore_gray	float	Grayscale value assigned to porosity equal to one.	required
background_porosity	float	Porosity assigned at solid_gray. This represents unresolved microporosity when nonzero.	0.0
clip	bool	Whether to clip extrapolated values to [background_porosity, 1].	True

Returns:

Type	Description
ndarray	Voxel-wise porosity field.

Source code in src/voids/image/porosity.py

def calibrated_porosity_from_grayscale(
    grayscale: np.ndarray,
    *,
    solid_gray: float,
    pore_gray: float,
    background_porosity: float = 0.0,
    clip: bool = True,
) -> np.ndarray:
    """Map grayscale intensity to porosity with a two-point linear calibration.

    Parameters
    ----------
    grayscale :
        Two- or three-dimensional grayscale image.
    solid_gray :
        Grayscale value assigned to the background or unresolved-microporosity
        porosity.
    pore_gray :
        Grayscale value assigned to porosity equal to one.
    background_porosity :
        Porosity assigned at ``solid_gray``. This represents unresolved
        microporosity when nonzero.
    clip :
        Whether to clip extrapolated values to
        ``[background_porosity, 1]``.

    Returns
    -------
    numpy.ndarray
        Voxel-wise porosity field.
    """

    arr = np.asarray(grayscale, dtype=float)
    normalize_shape(arr.shape, allowed_ndim=(2, 3))
    if not np.all(np.isfinite(arr)):
        raise ValueError("grayscale values must be finite")

    solid = float(solid_gray)
    pore = float(pore_gray)
    if not np.isfinite(solid) or not np.isfinite(pore):
        raise ValueError("solid_gray and pore_gray must be finite")
    if solid == pore:
        raise ValueError("solid_gray and pore_gray must differ")

    background = float(background_porosity)
    if not np.isfinite(background) or background < 0.0 or background > 1.0:
        raise ValueError("background_porosity must lie in [0, 1]")

    phi = background + (1.0 - background) * (arr - solid) / (pore - solid)
    if clip:
        phi = np.clip(phi, background, 1.0)
    return np.asarray(phi, dtype=float)

kozeny_carman_permeability

kozeny_carman_permeability(
    porosity,
    *,
    characteristic_length,
    kozeny_constant=180.0,
    solid_permeability=0.0,
    free_flow_permeability=np.inf,
    max_permeability=None,
)

Estimate permeability from porosity with a Kozeny-Carman closure.

Parameters:

Name	Type	Description	Default
porosity	ndarray	Two- or three-dimensional porosity field with entries in [0, 1].	required
characteristic_length	float	Characteristic length d in the same physical units used by the field. The returned permeability is in d**2 units.	required
kozeny_constant	float	Denominator constant. The default 180 is the common packed-sphere Kozeny-Carman value and should be treated as a calibration parameter for image-derived continuum fields.	180.0
solid_permeability	float	Permeability assigned at porosity == 0.	0.0
free_flow_permeability	float	Permeability assigned at porosity == 1. The mathematical limit is infinity; use a finite value when a downstream solver requires a cap.	inf
max_permeability	float | None	Optional finite cap applied after evaluating the closure.	None

Returns:

Type	Description
ndarray	Permeability field computed as k = d**2 * phi**3 / (C * (1 - phi)**2) for 0 < phi < 1.

Source code in src/voids/image/porosity.py

def kozeny_carman_permeability(
    porosity: np.ndarray,
    *,
    characteristic_length: float,
    kozeny_constant: float = 180.0,
    solid_permeability: float = 0.0,
    free_flow_permeability: float = np.inf,
    max_permeability: float | None = None,
) -> np.ndarray:
    """Estimate permeability from porosity with a Kozeny-Carman closure.

    Parameters
    ----------
    porosity :
        Two- or three-dimensional porosity field with entries in ``[0, 1]``.
    characteristic_length :
        Characteristic length ``d`` in the same physical units used by the field.
        The returned permeability is in ``d**2`` units.
    kozeny_constant :
        Denominator constant. The default ``180`` is the common packed-sphere
        Kozeny-Carman value and should be treated as a calibration parameter for
        image-derived continuum fields.
    solid_permeability :
        Permeability assigned at ``porosity == 0``.
    free_flow_permeability :
        Permeability assigned at ``porosity == 1``. The mathematical limit is
        infinity; use a finite value when a downstream solver requires a cap.
    max_permeability :
        Optional finite cap applied after evaluating the closure.

    Returns
    -------
    numpy.ndarray
        Permeability field computed as
        ``k = d**2 * phi**3 / (C * (1 - phi)**2)`` for ``0 < phi < 1``.
    """

    phi = _validate_porosity_values(porosity)
    length = _normalize_positive_scalar(characteristic_length, name="characteristic_length")
    constant = _normalize_positive_scalar(kozeny_constant, name="kozeny_constant")
    solid = _normalize_nonnegative_scalar(solid_permeability, name="solid_permeability")
    free = _normalize_nonnegative_scalar(
        free_flow_permeability,
        name="free_flow_permeability",
        allow_inf=True,
    )
    cap = (
        None
        if max_permeability is None
        else _normalize_positive_scalar(max_permeability, name="max_permeability")
    )

    permeability = np.empty_like(phi, dtype=float)
    solid_mask = phi <= 0.0
    free_mask = phi >= 1.0
    interior = ~(solid_mask | free_mask)

    permeability[solid_mask] = solid
    permeability[free_mask] = free
    interior_phi = phi[interior]
    permeability[interior] = length**2 * interior_phi**3 / (constant * (1.0 - interior_phi) ** 2)
    if cap is not None:
        permeability = np.minimum(permeability, cap)
    return permeability

kozeny_carman_inverse_permeability

kozeny_carman_inverse_permeability(
    porosity,
    *,
    characteristic_length,
    kozeny_constant=180.0,
    solid_inverse_permeability=np.inf,
    free_flow_inverse_permeability=0.0,
    max_inverse_permeability=None,
)

Estimate inverse permeability from porosity with a Kozeny-Carman closure.

Parameters:

Name	Type	Description	Default
porosity	ndarray	Two- or three-dimensional porosity field with entries in [0, 1].	required
characteristic_length	float	Characteristic length d in the same physical units used by the field.	required
kozeny_constant	float	Numerator constant. The default 180 is the common packed-sphere Kozeny-Carman value and should be treated as a calibration parameter for image-derived continuum fields.	180.0
solid_inverse_permeability	float	Inverse permeability assigned at porosity == 0. The mathematical limit is infinity.	inf
free_flow_inverse_permeability	float	Inverse permeability assigned at porosity == 1. The mathematical limit is zero.	0.0
max_inverse_permeability	float | None	Optional finite cap applied after evaluating the closure.	None

Returns:

Type	Description
ndarray	Inverse permeability field computed as k_inv = C * (1 - phi)**2 / (d**2 * phi**3) for 0 < phi < 1.

Source code in src/voids/image/porosity.py

def kozeny_carman_inverse_permeability(
    porosity: np.ndarray,
    *,
    characteristic_length: float,
    kozeny_constant: float = 180.0,
    solid_inverse_permeability: float = np.inf,
    free_flow_inverse_permeability: float = 0.0,
    max_inverse_permeability: float | None = None,
) -> np.ndarray:
    """Estimate inverse permeability from porosity with a Kozeny-Carman closure.

    Parameters
    ----------
    porosity :
        Two- or three-dimensional porosity field with entries in ``[0, 1]``.
    characteristic_length :
        Characteristic length ``d`` in the same physical units used by the field.
    kozeny_constant :
        Numerator constant. The default ``180`` is the common packed-sphere
        Kozeny-Carman value and should be treated as a calibration parameter for
        image-derived continuum fields.
    solid_inverse_permeability :
        Inverse permeability assigned at ``porosity == 0``. The mathematical
        limit is infinity.
    free_flow_inverse_permeability :
        Inverse permeability assigned at ``porosity == 1``. The mathematical
        limit is zero.
    max_inverse_permeability :
        Optional finite cap applied after evaluating the closure.

    Returns
    -------
    numpy.ndarray
        Inverse permeability field computed as
        ``k_inv = C * (1 - phi)**2 / (d**2 * phi**3)`` for ``0 < phi < 1``.
    """

    phi = _validate_porosity_values(porosity)
    length = _normalize_positive_scalar(characteristic_length, name="characteristic_length")
    constant = _normalize_positive_scalar(kozeny_constant, name="kozeny_constant")
    solid = _normalize_nonnegative_scalar(
        solid_inverse_permeability,
        name="solid_inverse_permeability",
        allow_inf=True,
    )
    free = _normalize_nonnegative_scalar(
        free_flow_inverse_permeability,
        name="free_flow_inverse_permeability",
    )
    cap = (
        None
        if max_inverse_permeability is None
        else _normalize_positive_scalar(max_inverse_permeability, name="max_inverse_permeability")
    )

    inverse_permeability = np.empty_like(phi, dtype=float)
    solid_mask = phi <= 0.0
    free_mask = phi >= 1.0
    interior = ~(solid_mask | free_mask)

    inverse_permeability[solid_mask] = solid
    inverse_permeability[free_mask] = free
    interior_phi = phi[interior]
    inverse_permeability[interior] = (
        constant * (1.0 - interior_phi) ** 2 / (length**2 * interior_phi**3)
    )
    if cap is not None:
        inverse_permeability = np.minimum(inverse_permeability, cap)
    return inverse_permeability

permeability_map_from_porosity

permeability_map_from_porosity(
    porosity_map,
    *,
    characteristic_length,
    kozeny_constant=180.0,
    solid_permeability=0.0,
    free_flow_permeability=np.inf,
    max_permeability=None,
    metadata=None,
)

Generate a Kozeny-Carman permeability map from a porosity map.

Parameters:

Name	Type	Description	Default
porosity_map	PorosityMap	Porosity map supplying the cell-wise phi field and spatial metadata.	required
characteristic_length	float		required
kozeny_constant	float		required
solid_permeability	float		required
free_flow_permeability	float	Closure parameters passed to :func:kozeny_carman_permeability.	inf
max_permeability	float	Closure parameters passed to :func:kozeny_carman_permeability.	inf
metadata	dict[str, Any] | None	Optional extra metadata merged into the generated map metadata.	None

Returns:

Type	Description
PermeabilityMap	Permeability field paired with the porosity-map grid metadata.

Source code in src/voids/image/porosity.py

def permeability_map_from_porosity(
    porosity_map: PorosityMap,
    *,
    characteristic_length: float,
    kozeny_constant: float = 180.0,
    solid_permeability: float = 0.0,
    free_flow_permeability: float = np.inf,
    max_permeability: float | None = None,
    metadata: dict[str, Any] | None = None,
) -> PermeabilityMap:
    """Generate a Kozeny-Carman permeability map from a porosity map.

    Parameters
    ----------
    porosity_map :
        Porosity map supplying the cell-wise ``phi`` field and spatial metadata.
    characteristic_length, kozeny_constant, solid_permeability,
    free_flow_permeability, max_permeability :
        Closure parameters passed to :func:`kozeny_carman_permeability`.
    metadata :
        Optional extra metadata merged into the generated map metadata.

    Returns
    -------
    PermeabilityMap
        Permeability field paired with the porosity-map grid metadata.
    """

    permeability = kozeny_carman_permeability(
        porosity_map.values,
        characteristic_length=characteristic_length,
        kozeny_constant=kozeny_constant,
        solid_permeability=solid_permeability,
        free_flow_permeability=free_flow_permeability,
        max_permeability=max_permeability,
    )
    length_unit = porosity_map.units.get("length", "m")
    map_metadata: dict[str, Any] = {
        "source_kind": "kozeny_carman_permeability",
        "porosity_source_metadata": porosity_map.metadata,
        "characteristic_length": float(characteristic_length),
        "kozeny_constant": float(kozeny_constant),
        "solid_permeability": float(solid_permeability),
        "free_flow_permeability": float(free_flow_permeability),
        "max_permeability": None if max_permeability is None else float(max_permeability),
    }
    if metadata:
        map_metadata.update(metadata)

    return PermeabilityMap(
        values=permeability,
        cell_size=porosity_map.cell_size,
        origin=porosity_map.origin,
        units={**porosity_map.units, "permeability": f"{length_unit}^2"},
        metadata=map_metadata,
    )

porosity_map_from_binary

porosity_map_from_binary(
    image,
    *,
    block_shape=None,
    voxel_size=1.0,
    image_is_void=True,
    strict=True,
    origin=None,
    units=None,
    metadata=None,
)

Compute a local porosity map from a segmented binary image.

Parameters:

Name	Type	Description	Default
image	ndarray	Binary 2D or 3D image. By default, True or 1 denotes void.	required
block_shape	Sequence[int] | None	Number of fine image voxels in each porosity-map cell. When omitted, each image voxel becomes one porosity-map cell.	None
voxel_size	float | Sequence[float]	Scalar or per-axis fine-image voxel spacing in physical units.	1.0
image_is_void	bool	If False, nonzero image values are interpreted as solid and are inverted before averaging.	True
strict	bool	If True, image shape must be exactly divisible by block_shape. If False, trailing partial blocks are trimmed.	True
origin	Sequence[float] | None	Optional physical and provenance metadata.	None
units	Sequence[float] | None	Optional physical and provenance metadata.	None
metadata	Sequence[float] | None	Optional physical and provenance metadata.	None

Returns:

Type	Description
PorosityMap	Cell-average local porosity map.

Source code in src/voids/image/porosity.py

def porosity_map_from_binary(
    image: np.ndarray,
    *,
    block_shape: Sequence[int] | None = None,
    voxel_size: float | Sequence[float] = 1.0,
    image_is_void: bool = True,
    strict: bool = True,
    origin: Sequence[float] | None = None,
    units: dict[str, str] | None = None,
    metadata: dict[str, Any] | None = None,
) -> PorosityMap:
    """Compute a local porosity map from a segmented binary image.

    Parameters
    ----------
    image :
        Binary 2D or 3D image. By default, ``True`` or ``1`` denotes void.
    block_shape :
        Number of fine image voxels in each porosity-map cell. When omitted,
        each image voxel becomes one porosity-map cell.
    voxel_size :
        Scalar or per-axis fine-image voxel spacing in physical units.
    image_is_void :
        If ``False``, nonzero image values are interpreted as solid and are
        inverted before averaging.
    strict :
        If ``True``, image shape must be exactly divisible by ``block_shape``.
        If ``False``, trailing partial blocks are trimmed.
    origin, units, metadata :
        Optional physical and provenance metadata.

    Returns
    -------
    PorosityMap
        Cell-average local porosity map.
    """

    void_mask = _as_binary_mask(image)
    if not image_is_void:
        void_mask = ~void_mask
    block = _normalize_block_shape(block_shape, ndim=void_mask.ndim)
    values, trimmed_shape = _block_mean(void_mask.astype(float), block_shape=block, strict=strict)
    voxel = _normalize_positive_tuple(voxel_size, ndim=void_mask.ndim, name="voxel_size")
    cell_size = tuple(v * b for v, b in zip(voxel, block, strict=True))

    map_metadata: dict[str, Any] = {
        "source_kind": "binary_void_fraction",
        "fine_shape": tuple(int(v) for v in void_mask.shape),
        "trimmed_shape": trimmed_shape,
        "block_shape": block,
        "image_is_void": bool(image_is_void),
        "strict": bool(strict),
    }
    if metadata:
        map_metadata.update(metadata)

    return PorosityMap(
        values=values,
        cell_size=cell_size,
        origin=origin,
        units=units or {"length": "m"},
        metadata=map_metadata,
    )

porosity_map_from_grayscale

porosity_map_from_grayscale(
    grayscale,
    *,
    solid_gray,
    pore_gray,
    background_porosity=0.0,
    block_shape=None,
    voxel_size=1.0,
    clip=True,
    strict=True,
    origin=None,
    units=None,
    metadata=None,
)

Compute a local porosity map from a calibrated grayscale image.

Parameters:

Name	Type	Description	Default
grayscale	ndarray	Two- or three-dimensional grayscale image.	required
solid_gray	float	Calibration parameters passed to :func:calibrated_porosity_from_grayscale.	required
pore_gray	float	Calibration parameters passed to :func:calibrated_porosity_from_grayscale.	required
background_porosity	float	Calibration parameters passed to :func:calibrated_porosity_from_grayscale.	required
clip	float	Calibration parameters passed to :func:calibrated_porosity_from_grayscale.	required
block_shape	Sequence[int] | None	Number of fine image voxels in each porosity-map cell. When omitted, each image voxel becomes one porosity-map cell.	None
voxel_size	float | Sequence[float]	Scalar or per-axis fine-image voxel spacing in physical units.	1.0
strict	bool	If True, image shape must be exactly divisible by block_shape. If False, trailing partial blocks are trimmed.	True
origin	Sequence[float] | None	Optional physical and provenance metadata.	None
units	Sequence[float] | None	Optional physical and provenance metadata.	None
metadata	Sequence[float] | None	Optional physical and provenance metadata.	None

Returns:

Type	Description
PorosityMap	Cell-average local porosity map.

Source code in src/voids/image/porosity.py

def porosity_map_from_grayscale(
    grayscale: np.ndarray,
    *,
    solid_gray: float,
    pore_gray: float,
    background_porosity: float = 0.0,
    block_shape: Sequence[int] | None = None,
    voxel_size: float | Sequence[float] = 1.0,
    clip: bool = True,
    strict: bool = True,
    origin: Sequence[float] | None = None,
    units: dict[str, str] | None = None,
    metadata: dict[str, Any] | None = None,
) -> PorosityMap:
    """Compute a local porosity map from a calibrated grayscale image.

    Parameters
    ----------
    grayscale :
        Two- or three-dimensional grayscale image.
    solid_gray, pore_gray, background_porosity, clip :
        Calibration parameters passed to
        :func:`calibrated_porosity_from_grayscale`.
    block_shape :
        Number of fine image voxels in each porosity-map cell. When omitted,
        each image voxel becomes one porosity-map cell.
    voxel_size :
        Scalar or per-axis fine-image voxel spacing in physical units.
    strict :
        If ``True``, image shape must be exactly divisible by ``block_shape``.
        If ``False``, trailing partial blocks are trimmed.
    origin, units, metadata :
        Optional physical and provenance metadata.

    Returns
    -------
    PorosityMap
        Cell-average local porosity map.
    """

    voxel_porosity = calibrated_porosity_from_grayscale(
        grayscale,
        solid_gray=solid_gray,
        pore_gray=pore_gray,
        background_porosity=background_porosity,
        clip=clip,
    )
    block = _normalize_block_shape(block_shape, ndim=voxel_porosity.ndim)
    values, trimmed_shape = _block_mean(voxel_porosity, block_shape=block, strict=strict)
    voxel = _normalize_positive_tuple(voxel_size, ndim=voxel_porosity.ndim, name="voxel_size")
    cell_size = tuple(v * b for v, b in zip(voxel, block, strict=True))

    map_metadata: dict[str, Any] = {
        "source_kind": "grayscale_linear_calibration",
        "fine_shape": tuple(int(v) for v in voxel_porosity.shape),
        "trimmed_shape": trimmed_shape,
        "block_shape": block,
        "solid_gray": float(solid_gray),
        "pore_gray": float(pore_gray),
        "background_porosity": float(background_porosity),
        "clip": bool(clip),
        "strict": bool(strict),
    }
    if metadata:
        map_metadata.update(metadata)

    return PorosityMap(
        values=values,
        cell_size=cell_size,
        origin=origin,
        units=units or {"length": "m"},
        metadata=map_metadata,
    )

save_porosity_map_hdf5

save_porosity_map_hdf5(porosity_map, path)

Write a porosity map to a compact HDF5 interchange file.

Parameters:

Name	Type	Description	Default
porosity_map	PorosityMap	Porosity map to serialize.	required
path	str | Path	Destination HDF5 file. Parent directories must already exist.	required

Notes

The stored array is named /porosity and contains cell-average porosity values. Root attributes store schema and calibration metadata as JSON.

Source code in src/voids/image/porosity.py

def save_porosity_map_hdf5(porosity_map: PorosityMap, path: str | Path) -> None:
    """Write a porosity map to a compact HDF5 interchange file.

    Parameters
    ----------
    porosity_map :
        Porosity map to serialize.
    path :
        Destination HDF5 file. Parent directories must already exist.

    Notes
    -----
    The stored array is named ``/porosity`` and contains cell-average porosity
    values. Root attributes store schema and calibration metadata as JSON.
    """

    path = Path(path)
    with h5py.File(path, "w") as f:
        f.attrs["schema_version"] = _SCHEMA_VERSION
        f.create_dataset("porosity", data=porosity_map.values)
        _write_json_attr(f, "metadata", porosity_map.to_metadata())

load_porosity_map_hdf5

load_porosity_map_hdf5(path)

Load a porosity map written by :func:save_porosity_map_hdf5.

Source code in src/voids/image/porosity.py

def load_porosity_map_hdf5(path: str | Path) -> PorosityMap:
    """Load a porosity map written by :func:`save_porosity_map_hdf5`."""

    path = Path(path)
    with h5py.File(path, "r") as f:
        schema_version = f.attrs.get("schema_version")
        if isinstance(schema_version, bytes):
            schema_version = schema_version.decode("utf-8")
        if schema_version != _SCHEMA_VERSION:
            raise ValueError(f"Unsupported porosity-map schema version {schema_version!r}")
        values = f["porosity"][()]
        metadata = _read_json_attr(f, "metadata", {})
    return PorosityMap.from_metadata(values, metadata)

save_permeability_map_hdf5

save_permeability_map_hdf5(permeability_map, path)

Write a permeability map to a compact HDF5 interchange file.

Parameters:

Name	Type	Description	Default
permeability_map	PermeabilityMap	Permeability map to serialize.	required
path	str | Path	Destination HDF5 file. Parent directories must already exist.	required

Notes

The stored array is named /permeability and contains cell-wise permeability values. Root attributes store schema and closure metadata as JSON.

Source code in src/voids/image/porosity.py

def save_permeability_map_hdf5(permeability_map: PermeabilityMap, path: str | Path) -> None:
    """Write a permeability map to a compact HDF5 interchange file.

    Parameters
    ----------
    permeability_map :
        Permeability map to serialize.
    path :
        Destination HDF5 file. Parent directories must already exist.

    Notes
    -----
    The stored array is named ``/permeability`` and contains cell-wise
    permeability values. Root attributes store schema and closure metadata as
    JSON.
    """

    path = Path(path)
    with h5py.File(path, "w") as f:
        f.attrs["schema_version"] = _PERMEABILITY_SCHEMA_VERSION
        f.create_dataset("permeability", data=permeability_map.values)
        _write_json_attr(f, "metadata", permeability_map.to_metadata())

load_permeability_map_hdf5

load_permeability_map_hdf5(path)

Load a permeability map written by :func:save_permeability_map_hdf5.

Source code in src/voids/image/porosity.py

def load_permeability_map_hdf5(path: str | Path) -> PermeabilityMap:
    """Load a permeability map written by :func:`save_permeability_map_hdf5`."""

    path = Path(path)
    with h5py.File(path, "r") as f:
        schema_version = f.attrs.get("schema_version")
        if isinstance(schema_version, bytes):
            schema_version = schema_version.decode("utf-8")
        if schema_version != _PERMEABILITY_SCHEMA_VERSION:
            raise ValueError(f"Unsupported permeability-map schema version {schema_version!r}")
        values = f["permeability"][()]
        metadata = _read_json_attr(f, "metadata", {})
    return PermeabilityMap.from_metadata(values, metadata)

---

Morphometry

The morphometry helpers compute local-thickness diameter maps for binary 2D/3D phase images. This API operates on an explicit phase mask supplied by the caller, requires isotropic voxel spacing, returns diameter-valued maps in the requested physical units, and summarizes values over the selected phase only.

For the scientific definition, radius-to-diameter conversion, backend method choices, and verification notes, see Local Thickness Morphometry.

voids.image.morphometry

LocalThicknessSummary dataclass

Summary statistics for a local-thickness diameter map.

Attributes:

Name	Type	Description
label	str	Human-readable phase label, for example "bone" or "marrow".
voxel_count	int	Number of phase voxels used for the statistics.
mean, std, p10, p50, p90, max		Summary statistics of the local-thickness diameter over phase voxels.
units	str	Physical units of the returned diameter values. Use "voxel" when voxel_size=1 is a dimensionless grid spacing.
method	str	PoreSpy local-thickness method used to generate the map.
voxel_size	float	Isotropic voxel edge length in units.

Notes

BoneJ-style trabecular thickness and separation are diameter quantities: each phase voxel is assigned the diameter of the largest sphere that fits inside the phase and contains that voxel. PoreSpy reports the corresponding local radius field, so voids multiplies by two and by voxel_size.

Source code in src/voids/image/morphometry.py

@dataclass(slots=True)
class LocalThicknessSummary:
    """Summary statistics for a local-thickness diameter map.

    Attributes
    ----------
    label :
        Human-readable phase label, for example ``"bone"`` or ``"marrow"``.
    voxel_count :
        Number of phase voxels used for the statistics.
    mean, std, p10, p50, p90, max :
        Summary statistics of the local-thickness diameter over phase voxels.
    units :
        Physical units of the returned diameter values. Use ``"voxel"`` when
        ``voxel_size=1`` is a dimensionless grid spacing.
    method :
        PoreSpy local-thickness method used to generate the map.
    voxel_size :
        Isotropic voxel edge length in ``units``.

    Notes
    -----
    BoneJ-style trabecular thickness and separation are diameter quantities:
    each phase voxel is assigned the diameter of the largest sphere that fits
    inside the phase and contains that voxel. PoreSpy reports the corresponding
    local radius field, so `voids` multiplies by two and by ``voxel_size``.
    """

    label: str
    voxel_count: int
    mean: float
    std: float
    p10: float
    p50: float
    p90: float
    max: float
    units: str
    method: str
    voxel_size: float

    def as_dict(self) -> dict[str, int | float | str]:
        """Return the summary as a JSON-serializable dictionary."""

        return asdict(self)

as_dict

as_dict()

Return the summary as a JSON-serializable dictionary.

Source code in src/voids/image/morphometry.py

def as_dict(self) -> dict[str, int | float | str]:
    """Return the summary as a JSON-serializable dictionary."""

    return asdict(self)

LocalThicknessResult dataclass

Local-thickness diameter map and summary for one binary phase.

Source code in src/voids/image/morphometry.py

@dataclass(slots=True)
class LocalThicknessResult:
    """Local-thickness diameter map and summary for one binary phase."""

    thickness_map: np.ndarray
    summary: LocalThicknessSummary

local_thickness_map

local_thickness_map(
    phase_mask,
    *,
    voxel_size=1.0,
    method="dt",
    sizes=64,
    smooth=True,
    approx=False,
    distance_map=None,
)

Compute a BoneJ-style local-thickness diameter map.

Parameters:

Name	Type	Description	Default
phase_mask	ndarray	Boolean or binary 2D/3D image where True marks the phase of interest.	required
voxel_size	float | Sequence[float]	Isotropic voxel edge length. A scalar is preferred; a sequence is accepted only when all entries are equal.	1.0
method	str	PoreSpy local-thickness backend. Common choices are "dt" for a practical distance-transform/opening approach and "imj" for a closer ImageJ-style sphere-insertion workflow.	'dt'
sizes	int | Sequence[float] | None	Radius sampling control forwarded to PoreSpy for "dt" and "conv" methods. None uses all unique distance-map values and can be expensive for large volumes.	64
smooth	bool	Controls forwarded to PoreSpy. approx is only used by "imj".	True
approx	bool	Controls forwarded to PoreSpy. approx is only used by "imj".	True
distance_map	ndarray | None	Optional Euclidean distance map in voxel units for phase_mask.	None

Returns:

Type	Description
ndarray	Local-thickness diameter map in units of voxel_size. Voxels outside phase_mask are zero.

Notes

The public BoneJ description defines local thickness as the diameter of the largest sphere contained in the object and containing the point. PoreSpy's local-thickness filters return the corresponding local radius field; this function converts that radius to a diameter in physical units.

Source code in src/voids/image/morphometry.py

def local_thickness_map(
    phase_mask: np.ndarray,
    *,
    voxel_size: float | Sequence[float] = 1.0,
    method: str = "dt",
    sizes: int | Sequence[float] | None = 64,
    smooth: bool = True,
    approx: bool = False,
    distance_map: np.ndarray | None = None,
) -> np.ndarray:
    """Compute a BoneJ-style local-thickness diameter map.

    Parameters
    ----------
    phase_mask :
        Boolean or binary 2D/3D image where ``True`` marks the phase of
        interest.
    voxel_size :
        Isotropic voxel edge length. A scalar is preferred; a sequence is
        accepted only when all entries are equal.
    method :
        PoreSpy local-thickness backend. Common choices are ``"dt"`` for a
        practical distance-transform/opening approach and ``"imj"`` for a
        closer ImageJ-style sphere-insertion workflow.
    sizes :
        Radius sampling control forwarded to PoreSpy for ``"dt"`` and
        ``"conv"`` methods. ``None`` uses all unique distance-map values and can
        be expensive for large volumes.
    smooth, approx :
        Controls forwarded to PoreSpy. ``approx`` is only used by ``"imj"``.
    distance_map :
        Optional Euclidean distance map in voxel units for ``phase_mask``.

    Returns
    -------
    numpy.ndarray
        Local-thickness diameter map in units of ``voxel_size``. Voxels outside
        ``phase_mask`` are zero.

    Notes
    -----
    The public BoneJ description defines local thickness as the diameter of the
    largest sphere contained in the object and containing the point. PoreSpy's
    local-thickness filters return the corresponding local radius field; this
    function converts that radius to a diameter in physical units.
    """

    mask = _binary_phase_mask(phase_mask)
    spacing = _isotropic_voxel_size(voxel_size, ndim=mask.ndim)
    normalized_method = str(method).strip().lower()

    if not bool(np.any(mask)):
        return np.zeros(mask.shape, dtype=float)

    dt = None if distance_map is None else _validate_distance_map(distance_map, shape=mask.shape)
    radius_map = _porespy_local_thickness(
        mask,
        dt=dt,
        method=normalized_method,
        smooth=bool(smooth),
        approx=bool(approx),
        sizes=sizes,
    )
    thickness = 2.0 * spacing * np.asarray(radius_map, dtype=float)
    return np.where(mask, thickness, 0.0)

summarize_local_thickness_map

summarize_local_thickness_map(
    thickness_map,
    phase_mask,
    *,
    label="phase",
    units="voxel",
    method="unknown",
    voxel_size=1.0,
)

Summarize local-thickness values over phase voxels.

Source code in src/voids/image/morphometry.py

def summarize_local_thickness_map(
    thickness_map: np.ndarray,
    phase_mask: np.ndarray,
    *,
    label: str = "phase",
    units: str = "voxel",
    method: str = "unknown",
    voxel_size: float = 1.0,
) -> LocalThicknessSummary:
    """Summarize local-thickness values over phase voxels."""

    mask = _binary_phase_mask(phase_mask)
    thickness = np.asarray(thickness_map, dtype=float)
    if thickness.shape != mask.shape:
        raise ValueError("thickness_map must have the same shape as phase_mask")
    if not bool(np.all(np.isfinite(thickness))):
        raise ValueError("thickness_map entries must be finite")
    if bool(np.any(thickness < 0.0)):
        raise ValueError("thickness_map entries must be nonnegative")

    values = np.asarray(thickness[mask], dtype=float)
    if values.size == 0:
        nan = float("nan")
        return LocalThicknessSummary(
            label=str(label),
            voxel_count=0,
            mean=nan,
            std=nan,
            p10=nan,
            p50=nan,
            p90=nan,
            max=nan,
            units=str(units),
            method=str(method),
            voxel_size=float(voxel_size),
        )

    return LocalThicknessSummary(
        label=str(label),
        voxel_count=int(values.size),
        mean=float(np.mean(values)),
        std=float(np.std(values)),
        p10=float(np.percentile(values, 10.0)),
        p50=float(np.percentile(values, 50.0)),
        p90=float(np.percentile(values, 90.0)),
        max=float(np.max(values)),
        units=str(units),
        method=str(method),
        voxel_size=float(voxel_size),
    )

local_thickness_analysis

local_thickness_analysis(
    phase_mask,
    *,
    voxel_size=1.0,
    units="voxel",
    label="phase",
    method="dt",
    sizes=64,
    smooth=True,
    approx=False,
    distance_map=None,
)

Compute a local-thickness diameter map and summary for one phase.

Source code in src/voids/image/morphometry.py

def local_thickness_analysis(
    phase_mask: np.ndarray,
    *,
    voxel_size: float | Sequence[float] = 1.0,
    units: str = "voxel",
    label: str = "phase",
    method: str = "dt",
    sizes: int | Sequence[float] | None = 64,
    smooth: bool = True,
    approx: bool = False,
    distance_map: np.ndarray | None = None,
) -> LocalThicknessResult:
    """Compute a local-thickness diameter map and summary for one phase."""

    mask = _binary_phase_mask(phase_mask)
    spacing = _isotropic_voxel_size(voxel_size, ndim=mask.ndim)
    thickness = local_thickness_map(
        mask,
        voxel_size=spacing,
        method=method,
        sizes=sizes,
        smooth=smooth,
        approx=approx,
        distance_map=distance_map,
    )
    summary = summarize_local_thickness_map(
        thickness,
        mask,
        label=label,
        units=units,
        method=method,
        voxel_size=spacing,
    )
    return LocalThicknessResult(thickness_map=thickness, summary=summary)

---

Network Extraction

voids.image.network_extraction

NetworkExtractionResult dataclass

Store outputs of an image-to-network extraction workflow.

Attributes:

Name	Type	Description
image	ndarray	Input phase image used for extraction.
voxel_size	float	Physical voxel edge length.
axis_lengths	dict[str, float]	Sample lengths by axis.
axis_areas	dict[str, float]	Cross-sectional areas normal to each axis.
flow_axis	str	Axis used for spanning subnetwork pruning.
network_dict	dict[str, object]	Intermediate extracted network mapping before conversion to :class:voids.core.network.Network.
sample	SampleGeometry	Sample geometry attached to the network.
provenance	Provenance	Extraction provenance metadata.
net_full	Network	Full imported network before spanning pruning.
net	Network	Axis-spanning subnetwork.
pore_indices	ndarray	Indices of retained pores in net_full.
throat_mask	ndarray	Mask of retained throats in net_full.
backend	str	Extraction backend identifier (currently "porespy").
backend_version	str | None	Backend version string when available.

Source code in src/voids/image/network_extraction.py

@dataclass(slots=True)
class NetworkExtractionResult:
    """Store outputs of an image-to-network extraction workflow.

    Attributes
    ----------
    image :
        Input phase image used for extraction.
    voxel_size :
        Physical voxel edge length.
    axis_lengths :
        Sample lengths by axis.
    axis_areas :
        Cross-sectional areas normal to each axis.
    flow_axis :
        Axis used for spanning subnetwork pruning.
    network_dict :
        Intermediate extracted network mapping before conversion to
        :class:`voids.core.network.Network`.
    sample :
        Sample geometry attached to the network.
    provenance :
        Extraction provenance metadata.
    net_full :
        Full imported network before spanning pruning.
    net :
        Axis-spanning subnetwork.
    pore_indices :
        Indices of retained pores in ``net_full``.
    throat_mask :
        Mask of retained throats in ``net_full``.
    backend :
        Extraction backend identifier (currently ``"porespy"``).
    backend_version :
        Backend version string when available.
    """

    image: np.ndarray
    voxel_size: float
    axis_lengths: dict[str, float]
    axis_areas: dict[str, float]
    flow_axis: str
    network_dict: dict[str, object]
    sample: SampleGeometry
    provenance: Provenance
    net_full: Network
    net: Network
    pore_indices: np.ndarray
    throat_mask: np.ndarray
    backend: str
    backend_version: str | None

NetworkConstructionResult dataclass

Store outputs of a general network-construction workflow.

This result type covers both image-based extraction backends and imported external-network backends such as Imperial CNM text files.

Source code in src/voids/image/network_extraction.py

@dataclass(slots=True)
class NetworkConstructionResult:
    """Store outputs of a general network-construction workflow.

    This result type covers both image-based extraction backends and imported
    external-network backends such as Imperial CNM text files.
    """

    backend: str
    flow_axis: str
    sample: SampleGeometry
    provenance: Provenance
    net_full: Network
    net: Network
    image: np.ndarray | None = None
    voxel_size: float | None = None
    axis_lengths: dict[str, float] | None = None
    axis_areas: dict[str, float] | None = None
    network_dict: dict[str, object] | None = None
    pore_indices: np.ndarray | None = None
    throat_mask: np.ndarray | None = None
    backend_version: str | None = None
    backend_details: dict[str, object] = field(default_factory=dict)

infer_sample_axes

infer_sample_axes(
    shape, *, voxel_size, axis_names=("x", "y", "z")
)

Infer per-axis counts, lengths, areas, and the longest flow axis.

Parameters:

Name	Type	Description	Default
shape	tuple[int, ...]	Image shape in voxel counts.	required
voxel_size	float	Edge length of one voxel in the target length unit.	required
axis_names	tuple[str, ...]	Axis labels mapped onto the image shape.	('x', 'y', 'z')

Returns:

Type	Description
tuple	(axis_counts, axis_lengths, axis_areas, flow_axis).

Source code in src/voids/image/network_extraction.py

def infer_sample_axes(
    shape: tuple[int, ...],
    *,
    voxel_size: float,
    axis_names: tuple[str, ...] = ("x", "y", "z"),
) -> tuple[dict[str, int], dict[str, float], dict[str, float], str]:
    """Infer per-axis counts, lengths, areas, and the longest flow axis.

    Parameters
    ----------
    shape :
        Image shape in voxel counts.
    voxel_size :
        Edge length of one voxel in the target length unit.
    axis_names :
        Axis labels mapped onto the image shape.

    Returns
    -------
    tuple
        ``(axis_counts, axis_lengths, axis_areas, flow_axis)``.
    """

    if voxel_size <= 0:
        raise ValueError("voxel_size must be positive")
    if len(shape) not in {2, 3}:
        raise ValueError("shape must have length 2 or 3")
    if len(axis_names) < len(shape):
        raise ValueError("axis_names must cover every image dimension")

    active_axes = axis_names[: len(shape)]
    axis_counts = {ax: int(n) for ax, n in zip(active_axes, shape)}
    axis_lengths = {ax: count * float(voxel_size) for ax, count in axis_counts.items()}
    axis_areas: dict[str, float] = {}
    for ax in active_axes:
        others = [other for other in active_axes if other != ax]
        area = float(voxel_size) ** max(len(others), 1)
        for other in others:
            area *= axis_counts[other]
        axis_areas[ax] = area
    flow_axis = max(axis_lengths, key=lambda ax: axis_lengths[ax])
    return axis_counts, axis_lengths, axis_areas, flow_axis

extract_spanning_pore_network

extract_spanning_pore_network(
    phases,
    *,
    voxel_size,
    backend="porespy",
    flow_axis=None,
    length_unit="m",
    pressure_unit="Pa",
    extraction_kwargs=None,
    provenance_notes=None,
    strict=True,
    geometry_repairs="imperial_export",
    repair_seed=0,
)

Extract, import, and prune an axis-spanning pore network from an image.

Parameters:

Name	Type	Description	Default
phases	ndarray	Binary or integer-labeled phase image where nonzero values are active phases passed to the extraction backend.	required
voxel_size	float	Edge length of one voxel in the declared length_unit.	required
backend	str	Image-to-network extraction backend. Currently supported values are "porespy", "snow2", "porespy_snow2", the calibrated approximation aliases "porespy_imperial", "imperial_snow2", and "snow2_imperial", the native PREGO backend "prego", plus the native maximal-ball aliases "maximal_ball", "native_maximal_ball", and "maxball".	'porespy'
flow_axis	str | None	Requested spanning axis. When omitted, the longest image axis is used.	None
length_unit	str	Units stored in resulting :class:SampleGeometry.	'm'
pressure_unit	str	Units stored in resulting :class:SampleGeometry.	'm'
extraction_kwargs	dict[str, object] | None	Keyword arguments forwarded to the extraction backend call. For the Imperial-calibrated snow2 aliases, user-supplied values override the built-in defaults sigma=1.0, r_max=4, and boundary_width=1. For the PREGO backend, supported keys are settings or prego_settings, distance_map, peaks, and regions_to_network_kwargs. PoreSpy-style backends ("porespy", the snow2 aliases, and "prego") also accept post-import transport keys flow_boundary_mode, boundary_axis, boundary_length_epsilon, boundary_radius_scale, and transport_geometry. For the native maximal-ball backend, supported keys are distance_map_backend, edt_parallel_threads, apply_boundary_clipping, flow_boundary_mode, boundary_axis, boundary_length_epsilon, boundary_radius_scale, throat_area_mode, throat_shape_factor_radius_mode, throat_anchor_mode, and either settings or maximal_ball_settings.	None
provenance_notes	dict[str, object] | None	Optional extra provenance metadata attached to the resulting network.	None
strict	bool	Forwarded to :func:voids.io.porespy.from_porespy.	True
geometry_repairs	str | None	Optional importer preprocessing mode. The default "imperial_export" applies the Imperial College export-style shape-factor repair heuristics during the PoreSpy-to-voids conversion.	'imperial_export'
repair_seed	int | None	Seed for any stochastic repair branch when geometry_repairs is not None.	0

Returns:

Type	Description
NetworkExtractionResult	Full and pruned networks together with intermediate metadata.

Notes

Current implementation uses PoreSpy's snow2 backend and normalizes accepted return styles into a standard network mapping before import. The calibrated Imperial-style aliases still use snow2, but start from a benchmark-tuned parameter profile that is closer to the committed pnextract reference cases than the plain default.

Source code in src/voids/image/network_extraction.py

def extract_spanning_pore_network(
    phases: np.ndarray,
    *,
    voxel_size: float,
    backend: str = "porespy",
    flow_axis: str | None = None,
    length_unit: str = "m",
    pressure_unit: str = "Pa",
    extraction_kwargs: dict[str, object] | None = None,
    provenance_notes: dict[str, object] | None = None,
    strict: bool = True,
    geometry_repairs: str | None = "imperial_export",
    repair_seed: int | None = 0,
) -> NetworkExtractionResult:
    """Extract, import, and prune an axis-spanning pore network from an image.

    Parameters
    ----------
    phases :
        Binary or integer-labeled phase image where nonzero values are active
        phases passed to the extraction backend.
    voxel_size :
        Edge length of one voxel in the declared ``length_unit``.
    backend :
        Image-to-network extraction backend. Currently supported values are
        ``"porespy"``, ``"snow2"``, ``"porespy_snow2"``, the calibrated
        approximation aliases ``"porespy_imperial"``, ``"imperial_snow2"``,
        and ``"snow2_imperial"``, the native PREGO backend ``"prego"``,
        plus the native maximal-ball aliases
        ``"maximal_ball"``, ``"native_maximal_ball"``, and ``"maxball"``.
    flow_axis :
        Requested spanning axis. When omitted, the longest image axis is used.
    length_unit, pressure_unit :
        Units stored in resulting :class:`SampleGeometry`.
    extraction_kwargs :
        Keyword arguments forwarded to the extraction backend call. For the
        Imperial-calibrated `snow2` aliases, user-supplied values override the
        built-in defaults ``sigma=1.0``, ``r_max=4``, and ``boundary_width=1``.
        For the PREGO backend, supported keys are ``settings`` or
        ``prego_settings``, ``distance_map``, ``peaks``, and
        ``regions_to_network_kwargs``. PoreSpy-style backends
        (``"porespy"``, the `snow2` aliases, and ``"prego"``) also accept
        post-import transport keys ``flow_boundary_mode``, ``boundary_axis``,
        ``boundary_length_epsilon``, ``boundary_radius_scale``, and
        ``transport_geometry``. For the native maximal-ball backend, supported
        keys are
        ``distance_map_backend``, ``edt_parallel_threads``,
        ``apply_boundary_clipping``,
        ``flow_boundary_mode``, ``boundary_axis``,
        ``boundary_length_epsilon``, ``boundary_radius_scale``,
        ``throat_area_mode``, ``throat_shape_factor_radius_mode``,
        ``throat_anchor_mode``, and either ``settings`` or
        ``maximal_ball_settings``.
    provenance_notes :
        Optional extra provenance metadata attached to the resulting network.
    strict :
        Forwarded to :func:`voids.io.porespy.from_porespy`.
    geometry_repairs :
        Optional importer preprocessing mode. The default
        ``"imperial_export"`` applies the Imperial College export-style
        shape-factor repair heuristics during the PoreSpy-to-``voids``
        conversion.
    repair_seed :
        Seed for any stochastic repair branch when ``geometry_repairs`` is not
        ``None``.

    Returns
    -------
    NetworkExtractionResult
        Full and pruned networks together with intermediate metadata.

    Notes
    -----
    Current implementation uses PoreSpy's ``snow2`` backend and normalizes
    accepted return styles into a standard network mapping before import. The
    calibrated Imperial-style aliases still use `snow2`, but start from a
    benchmark-tuned parameter profile that is closer to the committed
    `pnextract` reference cases than the plain default.
    """

    arr = np.asarray(phases, dtype=int)
    if arr.ndim not in {2, 3}:
        raise ValueError("phases must be a 2D or 3D integer image")

    _, axis_lengths, axis_areas, inferred_axis = infer_sample_axes(arr.shape, voxel_size=voxel_size)
    selected_axis = inferred_axis if flow_axis is None else flow_axis
    if selected_axis not in axis_lengths:
        raise ValueError(f"flow_axis '{selected_axis}' is not compatible with shape {arr.shape}")

    backend_normalized = _normalize_extraction_backend(backend)
    extraction_kwargs_for_backend = dict(extraction_kwargs or {})
    flow_boundary_mode = "direct"
    boundary_axis = selected_axis
    boundary_length_epsilon = 1.0e-300
    boundary_radius_scale = 1.1
    transport_geometry: str | None = None
    if backend_normalized in _PORESPY_STYLE_IMAGE_BACKENDS:
        flow_boundary_mode = _resolve_flow_boundary_mode(
            str(extraction_kwargs_for_backend.pop("flow_boundary_mode", "direct"))
        )
        boundary_axis = str(extraction_kwargs_for_backend.pop("boundary_axis", selected_axis))
        if boundary_axis not in axis_lengths:
            raise ValueError(
                f"boundary_axis '{boundary_axis}' is not compatible with shape {arr.shape}"
            )
        boundary_length_epsilon = float(
            cast(
                float | int | str,
                extraction_kwargs_for_backend.pop("boundary_length_epsilon", 1.0e-300),
            )
        )
        boundary_radius_scale = float(
            cast(float | int | str, extraction_kwargs_for_backend.pop("boundary_radius_scale", 1.1))
        )
        transport_geometry = _resolve_transport_geometry(
            extraction_kwargs_for_backend.pop("transport_geometry", None)
        )
    network_dict = _extract_network_dict(
        arr,
        backend=backend_normalized,
        voxel_size=float(voxel_size),
        extraction_kwargs=extraction_kwargs_for_backend or None,
        flow_axis=selected_axis,
    )
    importer_geometry_repairs = geometry_repairs
    if backend_normalized != "native_maximal_ball":
        network_dict = scale_porespy_geometry(network_dict, voxel_size=voxel_size)
        network_dict = ensure_cartesian_boundary_labels(network_dict, axes=(selected_axis,))
    else:
        importer_geometry_repairs = None

    shape_2d_or_3d = tuple(int(n) for n in arr.shape)
    bulk_shape: tuple[int, int, int] = (
        shape_2d_or_3d[0],
        shape_2d_or_3d[1],
        shape_2d_or_3d[2] if arr.ndim == 3 else 1,
    )
    sample = SampleGeometry(
        voxel_size=float(voxel_size),
        bulk_shape_voxels=bulk_shape,
        lengths=axis_lengths,
        cross_sections=axis_areas,
        units={"length": length_unit, "pressure": pressure_unit},
    )
    if backend_normalized in {"native_maximal_ball", "prego"}:
        source_version: str | None = _voids_version
    else:
        source_version = getattr(ps, "__version__", None)
    provenance = Provenance(
        source_kind="image_extraction",
        source_version=source_version,
        extraction_method=backend_normalized,
        random_seed=repair_seed if geometry_repairs is not None else None,
        user_notes=dict(provenance_notes or {}),
    )
    net_full = from_porespy(
        network_dict,
        sample=sample,
        provenance=provenance,
        strict=strict,
        geometry_repairs=importer_geometry_repairs,
        repair_seed=repair_seed,
    )
    if (
        backend_normalized in _PORESPY_STYLE_IMAGE_BACKENDS
        and flow_boundary_mode == "external_reservoir"
    ):
        net_full = _add_external_reservoirs_to_network(
            net_full,
            axis=boundary_axis,
            axis_length=axis_lengths[boundary_axis],
            voxel_size=float(voxel_size),
            boundary_length_epsilon=boundary_length_epsilon,
            boundary_radius_scale=boundary_radius_scale,
        )
    if transport_geometry == "pyramids_and_cuboids":
        net_full = _assign_pyramids_and_cuboids_transport_geometry(
            net_full,
            voxel_size=float(voxel_size),
        )
    net, pore_indices, throat_mask = spanning_subnetwork(net_full, axis=selected_axis)
    return NetworkExtractionResult(
        image=arr,
        voxel_size=float(voxel_size),
        axis_lengths=axis_lengths,
        axis_areas=axis_areas,
        flow_axis=selected_axis,
        network_dict=network_dict,
        sample=sample,
        provenance=provenance,
        net_full=net_full,
        net=net,
        pore_indices=pore_indices,
        throat_mask=throat_mask,
        backend=backend_normalized,
        backend_version=source_version,
    )

construct_spanning_network

construct_spanning_network(
    *,
    backend,
    phases=None,
    voxel_size=None,
    pnflow_cnm_prefix=None,
    pnflow_solver_box_compat=False,
    flow_axis=None,
    length_unit="m",
    pressure_unit="Pa",
    extraction_kwargs=None,
    provenance_notes=None,
    strict=True,
    geometry_repairs="imperial_export",
    repair_seed=0,
)

Construct a pore network from an image backend or imported CNM files.

Parameters:

Name	Type	Description	Default
backend	str	Construction backend identifier. Supported values include the existing image-extraction aliases "porespy", "snow2", "porespy_snow2", "prego", the native maximal-ball aliases "maximal_ball", "native_maximal_ball", "maxball", and the imported-network aliases "pnflow_cnm", "imperial_cnm", and "pnextract_cnm".	required
phases	ndarray | None	Required for image-based backends and forwarded to :func:extract_spanning_pore_network.	None
voxel_size	ndarray | None	Required for image-based backends and forwarded to :func:extract_spanning_pore_network.	None
pnflow_cnm_prefix	str | Path | None	Required for the Imperial CNM backend. This is the shared path prefix before the *_node*.dat and *_link*.dat suffixes.	None
pnflow_solver_box_compat	bool	If True and backend selects the Imperial CNM path, reproduce the checked-in pnflow solver-box preprocessing quirk so the imported network matches Imperial single-phase benchmark behavior. Leave this False for a generic CNM import.	False
flow_axis	str | None		None
length_unit	str | None		None
pressure_unit	str | None		None
extraction_kwargs	str | None		None
provenance_notes	str | None		None
strict	bool	Forwarded to the selected backend where applicable.	True
geometry_repairs	bool	Forwarded to the selected backend where applicable.	True
repair_seed	bool	Forwarded to the selected backend where applicable.	True

Returns:

Type	Description
NetworkConstructionResult	Unified network-construction result.

Source code in src/voids/image/network_extraction.py

def construct_spanning_network(
    *,
    backend: str,
    phases: np.ndarray | None = None,
    voxel_size: float | None = None,
    pnflow_cnm_prefix: str | Path | None = None,
    pnflow_solver_box_compat: bool = False,
    flow_axis: str | None = None,
    length_unit: str = "m",
    pressure_unit: str = "Pa",
    extraction_kwargs: dict[str, object] | None = None,
    provenance_notes: dict[str, object] | None = None,
    strict: bool = True,
    geometry_repairs: str | None = "imperial_export",
    repair_seed: int | None = 0,
) -> NetworkConstructionResult:
    """Construct a pore network from an image backend or imported CNM files.

    Parameters
    ----------
    backend :
        Construction backend identifier. Supported values include the existing
        image-extraction aliases ``"porespy"``, ``"snow2"``,
        ``"porespy_snow2"``, ``"prego"``, the native maximal-ball aliases
        ``"maximal_ball"``, ``"native_maximal_ball"``, ``"maxball"``, and the imported-network aliases
        ``"pnflow_cnm"``, ``"imperial_cnm"``, and ``"pnextract_cnm"``.
    phases, voxel_size :
        Required for image-based backends and forwarded to
        :func:`extract_spanning_pore_network`.
    pnflow_cnm_prefix :
        Required for the Imperial CNM backend. This is the shared path prefix
        before the ``*_node*.dat`` and ``*_link*.dat`` suffixes.
    pnflow_solver_box_compat :
        If ``True`` and ``backend`` selects the Imperial CNM path, reproduce
        the checked-in `pnflow` solver-box preprocessing quirk so the imported
        network matches Imperial single-phase benchmark behavior. Leave this
        ``False`` for a generic CNM import.
    flow_axis, length_unit, pressure_unit, extraction_kwargs, provenance_notes,
    strict, geometry_repairs, repair_seed :
        Forwarded to the selected backend where applicable.

    Returns
    -------
    NetworkConstructionResult
        Unified network-construction result.
    """

    backend_normalized = _normalize_construction_backend(backend)
    if backend_normalized in {
        "porespy_snow2",
        "porespy_snow2_imperial",
        "prego",
        "native_maximal_ball",
    }:
        if phases is None:
            raise ValueError("phases is required for the image-extraction backends")
        if voxel_size is None:
            raise ValueError("voxel_size is required for the image-extraction backends")
        extracted = extract_spanning_pore_network(
            phases,
            voxel_size=float(voxel_size),
            backend=backend_normalized,
            flow_axis=flow_axis,
            length_unit=length_unit,
            pressure_unit=pressure_unit,
            extraction_kwargs=extraction_kwargs,
            provenance_notes=provenance_notes,
            strict=strict,
            geometry_repairs=geometry_repairs,
            repair_seed=repair_seed,
        )
        return _construction_result_from_extraction(extracted)

    if pnflow_cnm_prefix is None:
        raise ValueError("pnflow_cnm_prefix is required for backend='pnflow_cnm'")
    selected_axis = "x" if flow_axis is None else str(flow_axis)
    if selected_axis != "x":
        raise ValueError("Imperial CNM construction currently supports only flow_axis='x'")

    imported = load_pnflow_cnm(
        pnflow_cnm_prefix,
        boundary_axis=selected_axis,
        length_unit=length_unit,
        pressure_unit=pressure_unit,
        pnflow_solver_box_compat=pnflow_solver_box_compat,
    )
    net = imported.net.copy()
    net.provenance = _merge_provenance_notes(net.provenance, provenance_notes)
    axis_lengths = dict(imported.box_lengths)
    axis_areas = {
        "x": imported.box_lengths["y"] * imported.box_lengths["z"],
        "y": imported.box_lengths["x"] * imported.box_lengths["z"],
        "z": imported.box_lengths["x"] * imported.box_lengths["y"],
    }
    return NetworkConstructionResult(
        backend=backend_normalized,
        flow_axis=selected_axis,
        sample=net.sample,
        provenance=net.provenance,
        net_full=net,
        net=net,
        image=None,
        voxel_size=None,
        axis_lengths=axis_lengths,
        axis_areas=axis_areas,
        network_dict=None,
        pore_indices=np.arange(net.Np, dtype=np.int64),
        throat_mask=np.ones(net.Nt, dtype=bool),
        backend_version=None,
        backend_details={
            "pnflow_cnm_prefix": str(Path(pnflow_cnm_prefix)),
            "n_physical_pores": int(imported.n_physical_pores),
            "n_boundary_mirror_pores": int(imported.n_boundary_mirror_pores),
        },
    )

---

Segmentation

voids.image.segmentation

VolumeCropResult dataclass

Store cylindrical-support cropping outputs from a grayscale volume.

Attributes:

Name	Type	Description
raw	ndarray	Original grayscale volume as float array.
specimen_mask	ndarray	Slice-wise support mask after hole filling.
common_mask	ndarray	Per-pixel intersection of support masks over all slices.
crop_bounds_yx	tuple[int, int, int, int]	Maximal common rectangle bounds (y0, y1, x0, x1).
cropped	ndarray	Cropped grayscale volume containing the common inscribed rectangle.

Source code in src/voids/image/segmentation.py

@dataclass(slots=True)
class VolumeCropResult:
    """Store cylindrical-support cropping outputs from a grayscale volume.

    Attributes
    ----------
    raw :
        Original grayscale volume as float array.
    specimen_mask :
        Slice-wise support mask after hole filling.
    common_mask :
        Per-pixel intersection of support masks over all slices.
    crop_bounds_yx :
        Maximal common rectangle bounds ``(y0, y1, x0, x1)``.
    cropped :
        Cropped grayscale volume containing the common inscribed rectangle.
    """

    raw: np.ndarray
    specimen_mask: np.ndarray
    common_mask: np.ndarray
    crop_bounds_yx: tuple[int, int, int, int]
    cropped: np.ndarray

GrayscaleSegmentationResult dataclass

Store grayscale preprocessing and binary segmentation outputs.

Attributes:

Name	Type	Description
crop	VolumeCropResult	Cylindrical-support crop outputs.
threshold	float	Threshold used for binarization.
binary	ndarray	Segmented binary volume encoded as void=1, solid=0.
void_phase	str	Phase polarity used for thresholding ("dark" or "bright").
threshold_method	str	Automatic method used when threshold was not explicitly supplied.

Source code in src/voids/image/segmentation.py

@dataclass(slots=True)
class GrayscaleSegmentationResult:
    """Store grayscale preprocessing and binary segmentation outputs.

    Attributes
    ----------
    crop :
        Cylindrical-support crop outputs.
    threshold :
        Threshold used for binarization.
    binary :
        Segmented binary volume encoded as ``void=1``, ``solid=0``.
    void_phase :
        Phase polarity used for thresholding (``"dark"`` or ``"bright"``).
    threshold_method :
        Automatic method used when threshold was not explicitly supplied.
    """

    crop: VolumeCropResult
    threshold: float
    binary: np.ndarray
    void_phase: str
    threshold_method: str

largest_true_rectangle

largest_true_rectangle(mask2d)

Return maximal-area axis-aligned rectangle fully contained in a mask.

Parameters:

Name	Type	Description	Default
mask2d	ndarray	Two-dimensional boolean support mask.	required

Returns:

Type	Description
tuple[int, int, int, int]	Rectangle bounds (y0, y1, x0, x1) in NumPy slicing convention.

Raises:

Type	Description
ValueError	If mask2d is not 2D or contains no True pixels.

Source code in src/voids/image/segmentation.py

def largest_true_rectangle(mask2d: np.ndarray) -> tuple[int, int, int, int]:
    """Return maximal-area axis-aligned rectangle fully contained in a mask.

    Parameters
    ----------
    mask2d :
        Two-dimensional boolean support mask.

    Returns
    -------
    tuple[int, int, int, int]
        Rectangle bounds ``(y0, y1, x0, x1)`` in NumPy slicing convention.

    Raises
    ------
    ValueError
        If ``mask2d`` is not 2D or contains no ``True`` pixels.
    """

    mask = np.asarray(mask2d, dtype=bool)
    if mask.ndim != 2:
        raise ValueError("mask2d must be a 2D boolean array")

    heights = [0] * mask.shape[1]
    best_area = 0
    best_bounds: tuple[int, int, int, int] | None = None
    for y in range(mask.shape[0]):
        for x in range(mask.shape[1]):
            heights[x] = heights[x] + 1 if mask[y, x] else 0
        stack: list[int] = []
        x = 0
        while x <= mask.shape[1]:
            cur = heights[x] if x < mask.shape[1] else 0
            if not stack or cur >= heights[stack[-1]]:
                stack.append(x)
                x += 1
            else:
                top = stack.pop()
                height = heights[top]
                left = stack[-1] + 1 if stack else 0
                width = x - left
                area = height * width
                if area > best_area:
                    best_area = area
                    best_bounds = (y + 1 - height, y + 1, left, x)
    if best_bounds is None:
        raise ValueError("mask2d does not contain any True pixels")
    return best_bounds

crop_nonzero_cylindrical_volume

crop_nonzero_cylindrical_volume(
    raw,
    *,
    background_value=0.0,
    show_progress=False,
    progress_desc=None,
)

Crop cylindrical specimen support to a common rectangular field of view.

Parameters:

Name	Type	Description	Default
raw	ndarray	Raw 3D grayscale volume.	required
background_value	float	Voxels strictly above this value are interpreted as specimen support before hole filling.	0.0
show_progress	bool	Whether to show progress bars for slice-wise operations.	False
progress_desc	str | None	Optional progress description string.	None

Returns:

Type	Description
VolumeCropResult	Structured crop result with masks, bounds, and cropped volume.

Source code in src/voids/image/segmentation.py

def crop_nonzero_cylindrical_volume(
    raw: np.ndarray,
    *,
    background_value: float = 0.0,
    show_progress: bool = False,
    progress_desc: str | None = None,
) -> VolumeCropResult:
    """Crop cylindrical specimen support to a common rectangular field of view.

    Parameters
    ----------
    raw :
        Raw 3D grayscale volume.
    background_value :
        Voxels strictly above this value are interpreted as specimen support
        before hole filling.
    show_progress :
        Whether to show progress bars for slice-wise operations.
    progress_desc :
        Optional progress description string.

    Returns
    -------
    VolumeCropResult
        Structured crop result with masks, bounds, and cropped volume.
    """

    arr = np.asarray(raw, dtype=float)
    if arr.ndim != 3:
        raise ValueError("raw must be a 3D grayscale volume")

    specimen_mask = np.zeros_like(arr, dtype=bool)
    iterator = _progress_iter(
        range(arr.shape[0]),
        show_progress=show_progress,
        desc=progress_desc or "Filling support mask slices",
        total=int(arr.shape[0]),
    )
    for i in iterator:
        specimen_mask[i] = ndi.binary_fill_holes(arr[i] > background_value)

    common_mask = np.asarray(specimen_mask.all(axis=0), dtype=bool)
    crop_bounds_yx = largest_true_rectangle(common_mask)
    y0, y1, x0, x1 = crop_bounds_yx
    cropped = arr[:, y0:y1, x0:x1]
    return VolumeCropResult(
        raw=arr,
        specimen_mask=specimen_mask,
        common_mask=common_mask,
        crop_bounds_yx=crop_bounds_yx,
        cropped=cropped,
    )

binarize_grayscale_volume

binarize_grayscale_volume(
    cropped,
    *,
    threshold=None,
    method="otsu",
    void_phase="dark",
)

Segment grayscale volume into binary void/solid phases.

Parameters:

Name	Type	Description	Default
cropped	ndarray	Cropped 3D grayscale volume.	required
threshold	float | None	Explicit threshold; when omitted, an automatic threshold is computed.	None
method	str	Automatic threshold method name. Supported values are "otsu", "li", "yen", "isodata", and "triangle".	'otsu'
void_phase	str	Which side of threshold corresponds to void: "dark" or "bright".	'dark'

Returns:

Type	Description
tuple[ndarray, float]	(binary, threshold_used) where binary is integer encoded as void=1 and solid=0.

Source code in src/voids/image/segmentation.py

def binarize_grayscale_volume(
    cropped: np.ndarray,
    *,
    threshold: float | None = None,
    method: str = "otsu",
    void_phase: str = "dark",
) -> tuple[np.ndarray, float]:
    """Segment grayscale volume into binary void/solid phases.

    Parameters
    ----------
    cropped :
        Cropped 3D grayscale volume.
    threshold :
        Explicit threshold; when omitted, an automatic threshold is computed.
    method :
        Automatic threshold method name. Supported values are ``"otsu"``,
        ``"li"``, ``"yen"``, ``"isodata"``, and ``"triangle"``.
    void_phase :
        Which side of threshold corresponds to void: ``"dark"`` or
        ``"bright"``.

    Returns
    -------
    tuple[numpy.ndarray, float]
        ``(binary, threshold_used)`` where ``binary`` is integer encoded as
        ``void=1`` and ``solid=0``.
    """

    arr = np.asarray(cropped, dtype=float)
    if arr.ndim != 3:
        raise ValueError("cropped must be a 3D grayscale volume")
    if void_phase not in {"dark", "bright"}:
        raise ValueError("void_phase must be either 'dark' or 'bright'")

    if threshold is None:
        if method not in _THRESHOLD_METHODS:
            raise ValueError(f"Unsupported threshold method '{method}'")
        threshold = float(_THRESHOLD_METHODS[method](arr))
    else:
        threshold = float(threshold)

    if void_phase == "dark":
        binary = (arr < threshold).astype(int)
    else:
        binary = (arr > threshold).astype(int)
    return binary, threshold

preprocess_grayscale_cylindrical_volume

preprocess_grayscale_cylindrical_volume(
    raw,
    *,
    background_value=0.0,
    threshold=None,
    threshold_method="otsu",
    void_phase="dark",
    show_progress=False,
    progress_desc=None,
)

Run cylindrical crop and grayscale segmentation in one workflow call.

Parameters:

Name	Type	Description	Default
raw	ndarray	Raw 3D grayscale specimen volume.	required
background_value	float	Background/support discriminator for cropping.	0.0
threshold	float | None	Explicit segmentation threshold.	None
threshold_method	str	Method used when threshold is omitted.	'otsu'
void_phase	str	Phase polarity selector for thresholding.	'dark'
show_progress	bool	Whether to request progress reporting.	False
progress_desc	str | None	Optional progress message.	None

Returns:

Type	Description
GrayscaleSegmentationResult	Crop metadata plus segmented binary volume.

Source code in src/voids/image/segmentation.py

def preprocess_grayscale_cylindrical_volume(
    raw: np.ndarray,
    *,
    background_value: float = 0.0,
    threshold: float | None = None,
    threshold_method: str = "otsu",
    void_phase: str = "dark",
    show_progress: bool = False,
    progress_desc: str | None = None,
) -> GrayscaleSegmentationResult:
    """Run cylindrical crop and grayscale segmentation in one workflow call.

    Parameters
    ----------
    raw :
        Raw 3D grayscale specimen volume.
    background_value :
        Background/support discriminator for cropping.
    threshold :
        Explicit segmentation threshold.
    threshold_method :
        Method used when ``threshold`` is omitted.
    void_phase :
        Phase polarity selector for thresholding.
    show_progress :
        Whether to request progress reporting.
    progress_desc :
        Optional progress message.

    Returns
    -------
    GrayscaleSegmentationResult
        Crop metadata plus segmented binary volume.
    """

    crop = crop_nonzero_cylindrical_volume(
        raw,
        background_value=background_value,
        show_progress=show_progress,
        progress_desc=progress_desc,
    )
    binary, used_threshold = binarize_grayscale_volume(
        crop.cropped,
        threshold=threshold,
        method=threshold_method,
        void_phase=void_phase,
    )
    return GrayscaleSegmentationResult(
        crop=crop,
        threshold=used_threshold,
        binary=binary,
        void_phase=void_phase,
        threshold_method=threshold_method,
    )

binarize_2d_with_voids

binarize_2d_with_voids(
    gray2d,
    *,
    threshold=None,
    method="otsu",
    void_phase="dark",
)

Segment a 2D grayscale image using the same thresholding policy as 3D.

Parameters:

Name	Type	Description	Default
gray2d	ndarray	Two-dimensional grayscale image.	required
threshold	float | None	Explicit threshold value.	None
method	str	Automatic threshold method when threshold is omitted.	'otsu'
void_phase	str	Which side of threshold corresponds to void.	'dark'

Returns:

Type	Description
tuple[ndarray, float]	(binary2d, threshold_used) with binary image encoded as integers in {0, 1}.

Source code in src/voids/image/segmentation.py

def binarize_2d_with_voids(
    gray2d: np.ndarray,
    *,
    threshold: float | None = None,
    method: str = "otsu",
    void_phase: str = "dark",
) -> tuple[np.ndarray, float]:
    """Segment a 2D grayscale image using the same thresholding policy as 3D.

    Parameters
    ----------
    gray2d :
        Two-dimensional grayscale image.
    threshold :
        Explicit threshold value.
    method :
        Automatic threshold method when ``threshold`` is omitted.
    void_phase :
        Which side of threshold corresponds to void.

    Returns
    -------
    tuple[numpy.ndarray, float]
        ``(binary2d, threshold_used)`` with binary image encoded as integers
        in ``{0, 1}``.
    """

    gray = np.asarray(gray2d, dtype=float)
    if gray.ndim != 2:
        raise ValueError("gray2d must be a 2D array")
    seg3d, threshold_used = binarize_grayscale_volume(
        gray[None, :, :],
        threshold=threshold,
        method=method,
        void_phase=void_phase,
    )
    return np.asarray(seg3d[0], dtype=int), float(threshold_used)

---

Image Connectivity

voids.image.connectivity

has_spanning_cluster

has_spanning_cluster(void_mask, axis_index)

Test whether void space percolates from one boundary to the opposite.

Parameters:

Name	Type	Description	Default
void_mask	ndarray	Binary image where True denotes void phase and False denotes solid phase. Supported dimensionalities are 2D and 3D.	required
axis_index	int	Flow/percolation axis index. For a 3D array (nx, ny, nz), values are typically 0 (x), 1 (y), or 2 (z).	required

Returns:

Type	Description
bool	True if at least one connected void component intersects both opposite boundary planes along axis_index.

Raises:

Type	Description
ValueError	If void_mask is not 2D/3D or if axis_index is invalid.

Notes

Connectivity is computed via :func:scipy.ndimage.label with its default structuring element (face-connected in 3D, edge-connected in 2D). The criterion is geometric percolation only; it does not assess hydraulic conductance magnitude.

Assumptions and limitations

• Boundaries are interpreted as the first and last index along the target axis.
• Periodic boundaries are not considered.
• Very thin connections may percolate topologically even if they are not representative of realistic transport in a given experiment.

Source code in src/voids/image/connectivity.py

def has_spanning_cluster(void_mask: np.ndarray, axis_index: int) -> bool:
    """Test whether void space percolates from one boundary to the opposite.

    Parameters
    ----------
    void_mask :
        Binary image where ``True`` denotes void phase and ``False`` denotes
        solid phase. Supported dimensionalities are 2D and 3D.
    axis_index :
        Flow/percolation axis index. For a 3D array ``(nx, ny, nz)``, values
        are typically ``0`` (x), ``1`` (y), or ``2`` (z).

    Returns
    -------
    bool
        ``True`` if at least one connected void component intersects both
        opposite boundary planes along ``axis_index``.

    Raises
    ------
    ValueError
        If ``void_mask`` is not 2D/3D or if ``axis_index`` is invalid.

    Notes
    -----
    Connectivity is computed via :func:`scipy.ndimage.label` with its default
    structuring element (face-connected in 3D, edge-connected in 2D). The
    criterion is geometric percolation only; it does not assess hydraulic
    conductance magnitude.

    Assumptions and limitations
    ---------------------------
    - Boundaries are interpreted as the first and last index along the target
      axis.
    - Periodic boundaries are not considered.
    - Very thin connections may percolate topologically even if they are not
      representative of realistic transport in a given experiment.
    """

    mask = np.asarray(void_mask, dtype=bool)
    if mask.ndim not in {2, 3}:
        raise ValueError("void_mask must be 2D or 3D")
    axis = validate_axis_index(axis_index=axis_index, ndim=mask.ndim)

    labels, n_labels = ndi.label(mask)
    if n_labels == 0:
        return False

    inlet_labels = np.unique(np.take(labels, indices=0, axis=axis))
    outlet_labels = np.unique(np.take(labels, indices=-1, axis=axis))
    inlet_labels = inlet_labels[inlet_labels > 0]
    outlet_labels = outlet_labels[outlet_labels > 0]
    return bool(np.intersect1d(inlet_labels, outlet_labels).size > 0)

has_spanning_cluster_2d

has_spanning_cluster_2d(void_mask, axis_index)

2D-specialized wrapper for axis-spanning connectivity checks.

Parameters:

Name	Type	Description	Default
void_mask	ndarray	Two-dimensional binary void mask.	required
axis_index	int	Axis along which percolation is tested (0 or 1).	required

Returns:

Type	Description
bool	True when at least one 2D connected void component spans the selected axis.

Raises:

Type	Description
ValueError	If void_mask is not 2D.

Notes

This function exists for notebook/API compatibility and delegates to :func:has_spanning_cluster after enforcing 2D inputs.

Source code in src/voids/image/connectivity.py

def has_spanning_cluster_2d(void_mask: np.ndarray, axis_index: int) -> bool:
    """2D-specialized wrapper for axis-spanning connectivity checks.

    Parameters
    ----------
    void_mask :
        Two-dimensional binary void mask.
    axis_index :
        Axis along which percolation is tested (``0`` or ``1``).

    Returns
    -------
    bool
        ``True`` when at least one 2D connected void component spans the
        selected axis.

    Raises
    ------
    ValueError
        If ``void_mask`` is not 2D.

    Notes
    -----
    This function exists for notebook/API compatibility and delegates to
    :func:`has_spanning_cluster` after enforcing 2D inputs.
    """

    mask = np.asarray(void_mask, dtype=bool)
    if mask.ndim != 2:
        raise ValueError("void_mask must be 2D for has_spanning_cluster_2d")
    return has_spanning_cluster(mask, axis_index=axis_index)