Generators

The voids.generators sub-package provides functions to create synthetic and mesh-based pore networks for controlled experiments and validation.

---

Network Generators

voids.generators.network

sample_depth

sample_depth(net)

Infer slab thickness for quasi-2D network models.

Parameters:

Name	Type	Description	Default
net	Network	Network whose sample geometry is inspected.	required

Returns:

Type	Description
float	Effective depth (thickness) used to convert 2D areas to volumes.

Notes

If sample.lengths['z'] exists, it is used directly. Otherwise depth is inferred from

depth = bulk_volume / (Lx * Ly).

This inference assumes orthogonal sample axes and consistent sample metadata.

Source code in src/voids/generators/network.py

def sample_depth(net: Network) -> float:
    """Infer slab thickness for quasi-2D network models.

    Parameters
    ----------
    net :
        Network whose sample geometry is inspected.

    Returns
    -------
    float
        Effective depth (thickness) used to convert 2D areas to volumes.

    Notes
    -----
    If ``sample.lengths['z']`` exists, it is used directly. Otherwise depth is
    inferred from

    ``depth = bulk_volume / (Lx * Ly)``.

    This inference assumes orthogonal sample axes and consistent sample
    metadata.
    """

    if "z" in net.sample.lengths:
        depth = float(net.sample.lengths["z"])
    else:
        lx = float(net.sample.length_for_axis("x"))
        ly = float(net.sample.length_for_axis("y"))
        bulk = float(net.sample.resolved_bulk_volume())
        depth = float(bulk / max(lx * ly, 1.0e-30))
    if depth <= 0:
        raise ValueError("sample depth must be positive")
    return depth

update_network_geometry_from_radii

update_network_geometry_from_radii(
    net,
    *,
    pore_radius,
    throat_radius,
    shape_factor=DEFAULT_G_REF,
    pore_length_fraction=0.45,
    min_core_fraction=0.05,
)

Recompute 3D pore/throat geometric fields from prescribed radii.

Parameters:

Name	Type	Description	Default
net	Network	Target network to update in place.	required
pore_radius	ndarray	Pore inscribed radii, shape (Np,).	required
throat_radius	ndarray	Throat inscribed radii, shape (Nt,).	required
shape_factor	float	Shape factor assigned uniformly to pores and throats. For circular conduits, the common value is 1 / (4 * pi).	DEFAULT_G_REF
pore_length_fraction	float	Fraction of center-to-center throat length used as an upper bound for each pore-body contribution.	0.45
min_core_fraction	float	Lower bound for throat core length as fraction of direct length.	0.05

Returns:

Type	Description
None	The network is modified in place.

Notes

Geometry model: - A = pi r^2 - P = 2 pi r - pore volume V_p = (4/3) pi r^3 - throat core volume V_t = A_t * L_core

Here L_core is computed from center distance minus pore-body lengths and clipped by min_core_fraction.

Scientific caveat

This is a synthetic geometric closure used for controlled studies; it is not a direct pore-space reconstruction from imaging.

Source code in src/voids/generators/network.py

def update_network_geometry_from_radii(
    net: Network,
    *,
    pore_radius: np.ndarray,
    throat_radius: np.ndarray,
    shape_factor: float = DEFAULT_G_REF,
    pore_length_fraction: float = 0.45,
    min_core_fraction: float = 0.05,
) -> None:
    """Recompute 3D pore/throat geometric fields from prescribed radii.

    Parameters
    ----------
    net :
        Target network to update in place.
    pore_radius :
        Pore inscribed radii, shape ``(Np,)``.
    throat_radius :
        Throat inscribed radii, shape ``(Nt,)``.
    shape_factor :
        Shape factor assigned uniformly to pores and throats. For circular
        conduits, the common value is ``1 / (4 * pi)``.
    pore_length_fraction :
        Fraction of center-to-center throat length used as an upper bound for
        each pore-body contribution.
    min_core_fraction :
        Lower bound for throat core length as fraction of direct length.

    Returns
    -------
    None
        The network is modified in place.

    Notes
    -----
    Geometry model:
    - ``A = pi r^2``
    - ``P = 2 pi r``
    - pore volume ``V_p = (4/3) pi r^3``
    - throat core volume ``V_t = A_t * L_core``

    Here ``L_core`` is computed from center distance minus pore-body lengths and
    clipped by ``min_core_fraction``.

    Scientific caveat
    -----------------
    This is a synthetic geometric closure used for controlled studies; it is not
    a direct pore-space reconstruction from imaging.
    """

    _validate_geometry_update_controls(
        shape_factor=shape_factor,
        pore_length_fraction=pore_length_fraction,
        min_core_fraction=min_core_fraction,
    )
    pore_r = _as_float_vector(pore_radius, expected_size=net.Np, name="pore_radius")
    throat_r = _as_float_vector(throat_radius, expected_size=net.Nt, name="throat_radius")
    if np.any(pore_r <= 0) or np.any(throat_r <= 0):
        raise ValueError("pore_radius and throat_radius must be strictly positive")

    pore_area = np.pi * pore_r**2
    pore_perimeter = 2.0 * np.pi * pore_r
    pore_volume = (4.0 / 3.0) * np.pi * pore_r**3

    net.pore["radius_inscribed"] = pore_r
    net.pore["diameter_inscribed"] = 2.0 * pore_r
    net.pore["area"] = pore_area
    net.pore["perimeter"] = pore_perimeter
    net.pore["volume"] = pore_volume
    net.pore["shape_factor"] = np.full(net.Np, float(shape_factor), dtype=float)

    conns = np.asarray(net.throat_conns, dtype=int)
    c0 = net.pore_coords[conns[:, 0]]
    c1 = net.pore_coords[conns[:, 1]]
    direct_length = np.linalg.norm(c1 - c0, axis=1)

    p1_length = np.minimum(pore_r[conns[:, 0]], float(pore_length_fraction) * direct_length)
    p2_length = np.minimum(pore_r[conns[:, 1]], float(pore_length_fraction) * direct_length)
    core_length = np.maximum(
        direct_length - p1_length - p2_length,
        float(min_core_fraction) * direct_length,
    )

    throat_area = np.pi * throat_r**2
    throat_perimeter = 2.0 * np.pi * throat_r
    throat_volume = throat_area * core_length

    net.throat["radius_inscribed"] = throat_r
    net.throat["diameter_inscribed"] = 2.0 * throat_r
    net.throat["area"] = throat_area
    net.throat["perimeter"] = throat_perimeter
    net.throat["shape_factor"] = np.full(net.Nt, float(shape_factor), dtype=float)
    net.throat["length"] = direct_length
    net.throat["direct_length"] = direct_length
    net.throat["pore1_length"] = p1_length
    net.throat["core_length"] = core_length
    net.throat["pore2_length"] = p2_length
    net.throat["volume"] = throat_volume

update_network_geometry_2d

update_network_geometry_2d(
    net,
    *,
    pore_radius,
    throat_radius,
    depth=None,
    shape_factor=DEFAULT_G_REF,
    pore_length_fraction=0.45,
    min_core_fraction=0.05,
)

Recompute geometry fields for 2D slab-like network models.

Parameters:

Name	Type	Description	Default
net	Network	Target network modified in place.	required
pore_radius	ndarray	Pore radii, shape (Np,).	required
throat_radius	ndarray	Throat radii, shape (Nt,).	required
depth	float | None	Slab thickness. If omitted, inferred via :func:sample_depth.	None
shape_factor	float	Same geometric controls as :func:update_network_geometry_from_radii.	DEFAULT_G_REF
pore_length_fraction	float	Same geometric controls as :func:update_network_geometry_from_radii.	DEFAULT_G_REF
min_core_fraction	float	Same geometric controls as :func:update_network_geometry_from_radii.	DEFAULT_G_REF

Returns:

Type	Description
None	The network is updated in place.

Notes

In quasi-2D mode, pore volume is approximated as V_p = A_p * depth with A_p = pi r^2. Throat volume uses the same A_t * L_core model used in the 3D helper.

Source code in src/voids/generators/network.py

def update_network_geometry_2d(
    net: Network,
    *,
    pore_radius: np.ndarray,
    throat_radius: np.ndarray,
    depth: float | None = None,
    shape_factor: float = DEFAULT_G_REF,
    pore_length_fraction: float = 0.45,
    min_core_fraction: float = 0.05,
) -> None:
    """Recompute geometry fields for 2D slab-like network models.

    Parameters
    ----------
    net :
        Target network modified in place.
    pore_radius :
        Pore radii, shape ``(Np,)``.
    throat_radius :
        Throat radii, shape ``(Nt,)``.
    depth :
        Slab thickness. If omitted, inferred via :func:`sample_depth`.
    shape_factor, pore_length_fraction, min_core_fraction :
        Same geometric controls as :func:`update_network_geometry_from_radii`.

    Returns
    -------
    None
        The network is updated in place.

    Notes
    -----
    In quasi-2D mode, pore volume is approximated as ``V_p = A_p * depth`` with
    ``A_p = pi r^2``. Throat volume uses the same ``A_t * L_core`` model used in
    the 3D helper.
    """

    _validate_geometry_update_controls(
        shape_factor=shape_factor,
        pore_length_fraction=pore_length_fraction,
        min_core_fraction=min_core_fraction,
    )
    pore_r = _as_float_vector(pore_radius, expected_size=net.Np, name="pore_radius")
    throat_r = _as_float_vector(throat_radius, expected_size=net.Nt, name="throat_radius")
    if np.any(pore_r <= 0) or np.any(throat_r <= 0):
        raise ValueError("pore_radius and throat_radius must be strictly positive")

    slab_depth = float(sample_depth(net) if depth is None else depth)
    if slab_depth <= 0:
        raise ValueError("depth must be positive")

    pore_area = np.pi * pore_r**2
    pore_perimeter = 2.0 * np.pi * pore_r
    pore_volume = pore_area * slab_depth

    net.pore["radius_inscribed"] = pore_r
    net.pore["diameter_inscribed"] = 2.0 * pore_r
    net.pore["area"] = pore_area
    net.pore["perimeter"] = pore_perimeter
    net.pore["volume"] = pore_volume
    net.pore["shape_factor"] = np.full(net.Np, float(shape_factor), dtype=float)

    conns = np.asarray(net.throat_conns, dtype=int)
    c0 = net.pore_coords[conns[:, 0]]
    c1 = net.pore_coords[conns[:, 1]]
    direct_length = np.linalg.norm(c1 - c0, axis=1)

    p1_length = np.minimum(pore_r[conns[:, 0]], float(pore_length_fraction) * direct_length)
    p2_length = np.minimum(pore_r[conns[:, 1]], float(pore_length_fraction) * direct_length)
    core_length = np.maximum(
        direct_length - p1_length - p2_length,
        float(min_core_fraction) * direct_length,
    )

    throat_area = np.pi * throat_r**2
    throat_perimeter = 2.0 * np.pi * throat_r
    throat_volume = throat_area * core_length

    net.throat["radius_inscribed"] = throat_r
    net.throat["diameter_inscribed"] = 2.0 * throat_r
    net.throat["area"] = throat_area
    net.throat["perimeter"] = throat_perimeter
    net.throat["shape_factor"] = np.full(net.Nt, float(shape_factor), dtype=float)
    net.throat["length"] = direct_length
    net.throat["direct_length"] = direct_length
    net.throat["pore1_length"] = p1_length
    net.throat["core_length"] = core_length
    net.throat["pore2_length"] = p2_length
    net.throat["volume"] = throat_volume

insert_vug_superpore

insert_vug_superpore(
    net,
    *,
    radii_xyz,
    center=None,
    shape_factor=DEFAULT_G_REF,
    connector_neighbor_weight=0.35,
    connector_vug_weight=0.1,
    connector_min_scale=1.1,
    connector_max_scale=2.8,
    pore_length_fraction=0.4,
    pore_length_radius_scale=0.8,
    min_core_fraction=0.02,
)

Insert one 3D vug as a super-pore replacing an ellipsoidal region.

Parameters:

Name	Type	Description	Default
net	Network	Source network.	required
radii_xyz	tuple[float, float, float]	Semi-axes of the ellipsoidal replacement region.	required
center	ndarray | tuple[float, float, float] | None	Optional vug center coordinate. Defaults to bounding-box center.	None
shape_factor	float	Shape factor assigned to inserted pore/throats.	DEFAULT_G_REF
connector_neighbor_weight	float	Linear weights used to estimate connector throat radii from neighboring pore radii and vug inscribed radius.	0.35
connector_vug_weight	float	Linear weights used to estimate connector throat radii from neighboring pore radii and vug inscribed radius.	0.35
connector_min_scale	float	Lower/upper clipping multipliers against median existing throat radius.	1.1
connector_max_scale	float	Lower/upper clipping multipliers against median existing throat radius.	1.1
pore_length_fraction	float	Controls for partitioning connector throat length into pore-side and core contributions.	0.4
pore_length_radius_scale	float	Controls for partitioning connector throat length into pore-side and core contributions.	0.4
min_core_fraction	float	Controls for partitioning connector throat length into pore-side and core contributions.	0.4

Returns:

Type	Description
tuple[Network, dict[str, object]]	(net_vug, metadata) where net_vug is the transformed network and metadata summarizes removed pores, boundary neighbors, and masks.

Raises:

Type	Description
ValueError	If radii or controls are invalid.
RuntimeError	If the selected region yields no usable interface connections.
KeyError	If required pore radius fields are absent.

Algorithm summary

1. Identify pores inside an ellipsoid.
2. Build induced subnetwork excluding those pores.
3. Add one new vug pore at center.
4. Connect vug pore to boundary neighbors of removed region.
5. Rebuild/extend geometry fields and labels.

Scientific caveat

The connector geometry is a heuristic closure model. It preserves topological continuity and plausible scale relations, but is not a unique physical derivation from first principles or direct imaging.

Source code in src/voids/generators/network.py

def insert_vug_superpore(
    net: Network,
    *,
    radii_xyz: tuple[float, float, float],
    center: np.ndarray | tuple[float, float, float] | None = None,
    shape_factor: float = DEFAULT_G_REF,
    connector_neighbor_weight: float = 0.35,
    connector_vug_weight: float = 0.10,
    connector_min_scale: float = 1.10,
    connector_max_scale: float = 2.80,
    pore_length_fraction: float = 0.40,
    pore_length_radius_scale: float = 0.80,
    min_core_fraction: float = 0.02,
) -> tuple[Network, dict[str, object]]:
    """Insert one 3D vug as a super-pore replacing an ellipsoidal region.

    Parameters
    ----------
    net :
        Source network.
    radii_xyz :
        Semi-axes of the ellipsoidal replacement region.
    center :
        Optional vug center coordinate. Defaults to bounding-box center.
    shape_factor :
        Shape factor assigned to inserted pore/throats.
    connector_neighbor_weight, connector_vug_weight :
        Linear weights used to estimate connector throat radii from neighboring
        pore radii and vug inscribed radius.
    connector_min_scale, connector_max_scale :
        Lower/upper clipping multipliers against median existing throat radius.
    pore_length_fraction, pore_length_radius_scale, min_core_fraction :
        Controls for partitioning connector throat length into pore-side and core
        contributions.

    Returns
    -------
    tuple[Network, dict[str, object]]
        ``(net_vug, metadata)`` where ``net_vug`` is the transformed network and
        ``metadata`` summarizes removed pores, boundary neighbors, and masks.

    Raises
    ------
    ValueError
        If radii or controls are invalid.
    RuntimeError
        If the selected region yields no usable interface connections.
    KeyError
        If required pore radius fields are absent.

    Algorithm summary
    -----------------
    1. Identify pores inside an ellipsoid.
    2. Build induced subnetwork excluding those pores.
    3. Add one new vug pore at ``center``.
    4. Connect vug pore to boundary neighbors of removed region.
    5. Rebuild/extend geometry fields and labels.

    Scientific caveat
    -----------------
    The connector geometry is a heuristic closure model. It preserves
    topological continuity and plausible scale relations, but is not a unique
    physical derivation from first principles or direct imaging.
    """

    rx, ry, rz = (float(radii_xyz[0]), float(radii_xyz[1]), float(radii_xyz[2]))
    if min(rx, ry, rz) <= 0:
        raise ValueError("All radii_xyz values must be positive")
    if shape_factor <= 0:
        raise ValueError("shape_factor must be positive")

    base = net.copy()
    coords = np.asarray(base.pore_coords, dtype=float)
    if center is None:
        center_arr = 0.5 * (coords.min(axis=0) + coords.max(axis=0))
    else:
        center_arr = np.asarray(center, dtype=float)
    if center_arr.shape != (3,):
        raise ValueError("center must have shape (3,)")

    inside = _ellipsoid_mask(coords, center=center_arr, radii_xyz=(rx, ry, rz))
    if not np.any(inside):
        nearest = int(np.argmin(np.linalg.norm(coords - center_arr[None, :], axis=1)))
        inside[nearest] = True

    conns = np.asarray(base.throat_conns, dtype=int)
    ci = inside[conns[:, 0]]
    cj = inside[conns[:, 1]]
    outside_neighbors = np.unique(np.concatenate([conns[ci & (~cj), 1], conns[(~ci) & cj, 0]]))
    if outside_neighbors.size == 0:
        raise RuntimeError("Vug insertion produced zero interface neighbors")

    subnet, kept_old_idx, _ = induced_subnetwork(base, ~inside)
    old_to_new = -np.ones(base.Np, dtype=int)
    old_to_new[kept_old_idx] = np.arange(kept_old_idx.size, dtype=int)
    boundary_new = old_to_new[outside_neighbors]
    boundary_new = np.unique(boundary_new[boundary_new >= 0])
    if boundary_new.size == 0:
        raise RuntimeError("No boundary pores survived after induced-subnetwork reduction")
    if "radius_inscribed" not in subnet.pore:
        raise KeyError("net.pore['radius_inscribed'] is required for super-pore insertion")

    r_eq = _equivalent_radius_3d((rx, ry, rz))
    r_ins = float(min(rx, ry, rz))

    vug_idx = subnet.Np
    new_coords = np.vstack([subnet.pore_coords, center_arr[None, :]])
    new_conns = np.column_stack([np.full(boundary_new.size, vug_idx, dtype=int), boundary_new])

    net_vug = subnet.copy()
    net_vug.pore_coords = new_coords
    net_vug.throat_conns = np.vstack([subnet.throat_conns, new_conns])

    pore_append = {
        "radius_inscribed": np.array([r_ins], dtype=float),
        "diameter_inscribed": np.array([2.0 * r_ins], dtype=float),
        "area": np.array([np.pi * r_eq**2], dtype=float),
        "perimeter": np.array([2.0 * np.pi * r_eq], dtype=float),
        "shape_factor": np.array([float(shape_factor)], dtype=float),
        "volume": np.array([(4.0 / 3.0) * np.pi * rx * ry * rz], dtype=float),
    }
    _extend_entity_fields(net_vug.pore, n_before=subnet.Np, n_append=1, append_fields=pore_append)

    boundary_coords = net_vug.pore_coords[boundary_new]
    direct_length = np.linalg.norm(boundary_coords - center_arr[None, :], axis=1)
    direct_length = np.maximum(direct_length, 1.0e-12)

    throat_r_median = _median_throat_radius(subnet, fallback=max(0.1 * r_ins, 1.0e-12))
    neigh_r = np.asarray(net_vug.pore["radius_inscribed"], dtype=float)[boundary_new]
    conn_radius = np.clip(
        connector_neighbor_weight * neigh_r + connector_vug_weight * r_ins,
        connector_min_scale * throat_r_median,
        connector_max_scale * throat_r_median,
    )

    p1_length = np.minimum(pore_length_fraction * direct_length, pore_length_radius_scale * r_ins)
    p2_length = np.minimum(pore_length_fraction * direct_length, pore_length_radius_scale * neigh_r)
    core_length = np.maximum(
        direct_length - p1_length - p2_length,
        min_core_fraction * direct_length,
    )

    conn_area = np.pi * conn_radius**2
    conn_perimeter = 2.0 * np.pi * conn_radius
    conn_volume = conn_area * core_length

    throat_append = {
        "radius_inscribed": conn_radius,
        "diameter_inscribed": 2.0 * conn_radius,
        "area": conn_area,
        "perimeter": conn_perimeter,
        "shape_factor": np.full(boundary_new.size, float(shape_factor), dtype=float),
        "length": direct_length,
        "direct_length": direct_length,
        "pore1_length": p1_length,
        "core_length": core_length,
        "pore2_length": p2_length,
        "volume": conn_volume,
    }
    _extend_entity_fields(
        net_vug.throat,
        n_before=subnet.Nt,
        n_append=boundary_new.size,
        append_fields=throat_append,
    )

    _extend_labels(net_vug.pore_labels, n_before=subnet.Np, n_append=1)
    vug_label = np.zeros(net_vug.Np, dtype=bool)
    vug_label[vug_idx] = True
    net_vug.pore_labels["vug"] = vug_label

    _extend_labels(net_vug.throat_labels, n_before=subnet.Nt, n_append=boundary_new.size)
    vug_conn = np.zeros(net_vug.Nt, dtype=bool)
    vug_conn[subnet.Nt :] = True
    net_vug.throat_labels["vug_connection"] = vug_conn

    net_vug.extra = {
        **net_vug.extra,
        "vug_radii_xyz_m": (rx, ry, rz),
        "vug_equivalent_radius_m": float(r_eq),
        "vug_removed_pores": int(inside.sum()),
        "vug_boundary_neighbors": int(boundary_new.size),
    }

    metadata: dict[str, object] = {
        "inside_mask_original": inside,
        "outside_neighbors_original": outside_neighbors,
        "removed_pores": int(inside.sum()),
        "boundary_neighbors": int(boundary_new.size),
        "equivalent_radius_m": float(r_eq),
    }
    return net_vug, metadata

insert_vug_superpore_2d

insert_vug_superpore_2d(
    net,
    *,
    radii_xy,
    center_xy=None,
    depth=None,
    shape_factor=DEFAULT_G_REF,
    connector_neighbor_weight=0.4,
    connector_vug_weight=0.12,
    connector_min_scale=1.05,
    connector_max_scale=3.0,
    pore_length_fraction=0.42,
    pore_length_radius_scale=0.85,
    min_core_fraction=0.02,
)

Insert one 2D vug as a super-pore replacing an elliptical region.

Parameters:

Name	Type	Description	Default
net	Network	Source network (typically quasi-2D mesh network).	required
radii_xy	tuple[float, float]	Semi-axes of the elliptical replacement region in the XY plane.	required
center_xy	tuple[float, float] | None	Optional XY center. Defaults to XY bounding-box center.	None
depth	float | None	Optional slab thickness; inferred from sample metadata when omitted.	None
shape_factor	float	Shape factor assigned to inserted pore/throats.	DEFAULT_G_REF
connector_neighbor_weight	float	Weights used to estimate connector throat radii.	0.4
connector_vug_weight	float	Weights used to estimate connector throat radii.	0.4
connector_min_scale	float	Connector radius clipping multipliers against median throat radius.	1.05
connector_max_scale	float	Connector radius clipping multipliers against median throat radius.	1.05
pore_length_fraction	float	Controls for connector throat length partitioning.	0.42
pore_length_radius_scale	float	Controls for connector throat length partitioning.	0.42
min_core_fraction	float	Controls for connector throat length partitioning.	0.42

Returns:

Type	Description
tuple[Network, dict[str, object]]	Updated network and diagnostic metadata.

Notes

The inserted 2D vug pore uses: - equivalent radius from area matching, - pore volume pi * rx * ry * depth, - connector geometry from the same heuristic strategy used in 3D.

Source code in src/voids/generators/network.py

def insert_vug_superpore_2d(
    net: Network,
    *,
    radii_xy: tuple[float, float],
    center_xy: tuple[float, float] | None = None,
    depth: float | None = None,
    shape_factor: float = DEFAULT_G_REF,
    connector_neighbor_weight: float = 0.40,
    connector_vug_weight: float = 0.12,
    connector_min_scale: float = 1.05,
    connector_max_scale: float = 3.00,
    pore_length_fraction: float = 0.42,
    pore_length_radius_scale: float = 0.85,
    min_core_fraction: float = 0.02,
) -> tuple[Network, dict[str, object]]:
    """Insert one 2D vug as a super-pore replacing an elliptical region.

    Parameters
    ----------
    net :
        Source network (typically quasi-2D mesh network).
    radii_xy :
        Semi-axes of the elliptical replacement region in the XY plane.
    center_xy :
        Optional XY center. Defaults to XY bounding-box center.
    depth :
        Optional slab thickness; inferred from sample metadata when omitted.
    shape_factor :
        Shape factor assigned to inserted pore/throats.
    connector_neighbor_weight, connector_vug_weight :
        Weights used to estimate connector throat radii.
    connector_min_scale, connector_max_scale :
        Connector radius clipping multipliers against median throat radius.
    pore_length_fraction, pore_length_radius_scale, min_core_fraction :
        Controls for connector throat length partitioning.

    Returns
    -------
    tuple[Network, dict[str, object]]
        Updated network and diagnostic metadata.

    Notes
    -----
    The inserted 2D vug pore uses:
    - equivalent radius from area matching,
    - pore volume ``pi * rx * ry * depth``,
    - connector geometry from the same heuristic strategy used in 3D.
    """

    rx, ry = (float(radii_xy[0]), float(radii_xy[1]))
    if min(rx, ry) <= 0:
        raise ValueError("All radii_xy values must be positive")
    if shape_factor <= 0:
        raise ValueError("shape_factor must be positive")

    base = net.copy()
    coords = np.asarray(base.pore_coords, dtype=float)
    if center_xy is None:
        center_xy = (
            float(0.5 * (coords[:, 0].min() + coords[:, 0].max())),
            float(0.5 * (coords[:, 1].min() + coords[:, 1].max())),
        )

    inside = _ellipse_mask_2d(coords, center_xy=center_xy, radii_xy=(rx, ry))
    if not np.any(inside):
        nearest = int(
            np.argmin((coords[:, 0] - center_xy[0]) ** 2 + (coords[:, 1] - center_xy[1]) ** 2)
        )
        inside[nearest] = True

    conns = np.asarray(base.throat_conns, dtype=int)
    ci = inside[conns[:, 0]]
    cj = inside[conns[:, 1]]
    outside_neighbors = np.unique(np.concatenate([conns[ci & (~cj), 1], conns[(~ci) & cj, 0]]))
    if outside_neighbors.size == 0:
        raise RuntimeError("Vug insertion produced zero interface neighbors")

    subnet, kept_old_idx, _ = induced_subnetwork(base, ~inside)
    old_to_new = -np.ones(base.Np, dtype=int)
    old_to_new[kept_old_idx] = np.arange(kept_old_idx.size, dtype=int)
    boundary_new = old_to_new[outside_neighbors]
    boundary_new = np.unique(boundary_new[boundary_new >= 0])
    if boundary_new.size == 0:
        raise RuntimeError("No interface pores remained after induced-subnetwork reduction")
    if "radius_inscribed" not in subnet.pore:
        raise KeyError("net.pore['radius_inscribed'] is required for super-pore insertion")

    r_eq = _equivalent_radius_2d((rx, ry))
    r_ins = float(min(rx, ry))
    slab_depth = float(sample_depth(subnet) if depth is None else depth)
    if slab_depth <= 0:
        raise ValueError("depth must be positive")

    vug_idx = subnet.Np
    zref = float(np.mean(subnet.pore_coords[:, 2]))
    center3 = np.array([center_xy[0], center_xy[1], zref], dtype=float)
    new_coords = np.vstack([subnet.pore_coords, center3[None, :]])
    new_conns = np.column_stack([np.full(boundary_new.size, vug_idx, dtype=int), boundary_new])

    net_vug = subnet.copy()
    net_vug.pore_coords = new_coords
    net_vug.throat_conns = np.vstack([subnet.throat_conns, new_conns])

    pore_append = {
        "radius_inscribed": np.array([r_ins], dtype=float),
        "diameter_inscribed": np.array([2.0 * r_ins], dtype=float),
        "area": np.array([np.pi * r_eq**2], dtype=float),
        "perimeter": np.array([2.0 * np.pi * r_eq], dtype=float),
        "shape_factor": np.array([float(shape_factor)], dtype=float),
        "volume": np.array([np.pi * rx * ry * slab_depth], dtype=float),
    }
    _extend_entity_fields(net_vug.pore, n_before=subnet.Np, n_append=1, append_fields=pore_append)

    boundary_coords = net_vug.pore_coords[boundary_new]
    direct_length = np.linalg.norm(boundary_coords - center3[None, :], axis=1)
    direct_length = np.maximum(direct_length, 1.0e-12)

    throat_r_median = _median_throat_radius(subnet, fallback=max(0.1 * r_ins, 1.0e-12))
    neigh_r = np.asarray(net_vug.pore["radius_inscribed"], dtype=float)[boundary_new]
    conn_radius = np.clip(
        connector_neighbor_weight * neigh_r + connector_vug_weight * r_ins,
        connector_min_scale * throat_r_median,
        connector_max_scale * throat_r_median,
    )

    p1_length = np.minimum(pore_length_fraction * direct_length, pore_length_radius_scale * r_ins)
    p2_length = np.minimum(pore_length_fraction * direct_length, pore_length_radius_scale * neigh_r)
    core_length = np.maximum(
        direct_length - p1_length - p2_length,
        min_core_fraction * direct_length,
    )

    conn_area = np.pi * conn_radius**2
    conn_perimeter = 2.0 * np.pi * conn_radius
    conn_volume = conn_area * core_length

    throat_append = {
        "radius_inscribed": conn_radius,
        "diameter_inscribed": 2.0 * conn_radius,
        "area": conn_area,
        "perimeter": conn_perimeter,
        "shape_factor": np.full(boundary_new.size, float(shape_factor), dtype=float),
        "length": direct_length,
        "direct_length": direct_length,
        "pore1_length": p1_length,
        "core_length": core_length,
        "pore2_length": p2_length,
        "volume": conn_volume,
    }
    _extend_entity_fields(
        net_vug.throat,
        n_before=subnet.Nt,
        n_append=boundary_new.size,
        append_fields=throat_append,
    )

    _extend_labels(net_vug.pore_labels, n_before=subnet.Np, n_append=1)
    vug_label = np.zeros(net_vug.Np, dtype=bool)
    vug_label[vug_idx] = True
    net_vug.pore_labels["vug"] = vug_label

    _extend_labels(net_vug.throat_labels, n_before=subnet.Nt, n_append=boundary_new.size)
    vug_conn = np.zeros(net_vug.Nt, dtype=bool)
    vug_conn[subnet.Nt :] = True
    net_vug.throat_labels["vug_connection"] = vug_conn

    net_vug.extra = {
        **net_vug.extra,
        "vug_radii_xy_m": (rx, ry),
        "vug_equivalent_radius_m": float(r_eq),
        "vug_removed_pores": int(inside.sum()),
        "vug_boundary_neighbors": int(boundary_new.size),
    }

    metadata: dict[str, object] = {
        "inside_mask_original": inside,
        "removed_pores": int(inside.sum()),
        "boundary_neighbors": int(boundary_new.size),
        "equivalent_radius_m": float(r_eq),
        "center_xy": center_xy,
    }
    return net_vug, metadata

---

Porous Image Generators

voids.generators.porous_image

generate_spanning_blobs_matrix

generate_spanning_blobs_matrix(
    *,
    shape,
    porosity,
    blobiness,
    axis_index,
    seed_start,
    max_tries,
)

Generate a percolating porous matrix using PoreSpy's blobs model.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	Image shape in voxels. Supports 2D and 3D.	required
porosity	float	Target void fraction passed to porespy.generators.blobs.	required
blobiness	float	Correlation-length control used by PoreSpy; larger values generally produce smoother/larger features.	required
axis_index	int	Axis used for percolation (spanning) acceptance.	required
seed_start	int	Initial random seed. Subsequent attempts use seed_start + i.	required
max_tries	int	Maximum number of seeds to test.	required

Returns:

Type	Description
tuple[ndarray, int]	(matrix, seed_used) where matrix is boolean with True as void and seed_used is the accepted seed.

Raises:

Type	Description
ValueError	If input controls are outside valid ranges.
RuntimeError	If no percolating realization is found in max_tries attempts.

Scientific interpretation

The acceptance criterion is topological percolation only. It ensures a connected pathway exists, but does not guarantee a target hydraulic conductance or morphological realism for specific rock classes.

Source code in src/voids/generators/porous_image.py

def generate_spanning_blobs_matrix(
    *,
    shape: Sequence[int],
    porosity: float,
    blobiness: float,
    axis_index: int,
    seed_start: int,
    max_tries: int,
) -> tuple[np.ndarray, int]:
    """Generate a percolating porous matrix using PoreSpy's `blobs` model.

    Parameters
    ----------
    shape :
        Image shape in voxels. Supports 2D and 3D.
    porosity :
        Target void fraction passed to ``porespy.generators.blobs``.
    blobiness :
        Correlation-length control used by PoreSpy; larger values generally
        produce smoother/larger features.
    axis_index :
        Axis used for percolation (spanning) acceptance.
    seed_start :
        Initial random seed. Subsequent attempts use ``seed_start + i``.
    max_tries :
        Maximum number of seeds to test.

    Returns
    -------
    tuple[numpy.ndarray, int]
        ``(matrix, seed_used)`` where ``matrix`` is boolean with ``True`` as
        void and ``seed_used`` is the accepted seed.

    Raises
    ------
    ValueError
        If input controls are outside valid ranges.
    RuntimeError
        If no percolating realization is found in ``max_tries`` attempts.

    Scientific interpretation
    -------------------------
    The acceptance criterion is topological percolation only. It ensures a
    connected pathway exists, but does not guarantee a target hydraulic
    conductance or morphological realism for specific rock classes.
    """

    dims = normalize_shape(shape, allowed_ndim=(2, 3))
    axis = validate_axis_index(axis_index=axis_index, ndim=len(dims))
    if not (0.0 < porosity < 1.0):
        raise ValueError("porosity must be in (0, 1)")
    if blobiness <= 0:
        raise ValueError("blobiness must be positive")
    if max_tries < 1:
        raise ValueError("max_tries must be >= 1")

    for i in range(max_tries):
        seed = int(seed_start + i)
        matrix = np.asarray(
            ps.generators.blobs(
                shape=dims,
                porosity=float(porosity),
                blobiness=float(blobiness),
                seed=seed,
            ),
            dtype=bool,
        )
        if has_spanning_cluster(matrix, axis_index=axis):
            return matrix, seed
    raise RuntimeError(
        f"Could not generate spanning blobs matrix for seed_start={int(seed_start)} "
        f"after {int(max_tries)} trials"
    )

generate_connected_matrix

generate_connected_matrix(
    *,
    shape,
    porosity,
    blobiness,
    axis_index,
    seed_start,
    max_tries,
    show_progress=False,
    progress_desc=None,
)

Backward-compatible wrapper for notebook 08 API signatures.

Parameters:

Name	Type	Description	Default
shape	tuple[int, int, int]	Same physical meaning as :func:generate_spanning_blobs_matrix.	required
porosity	tuple[int, int, int]	Same physical meaning as :func:generate_spanning_blobs_matrix.	required
blobiness	tuple[int, int, int]	Same physical meaning as :func:generate_spanning_blobs_matrix.	required
axis_index	tuple[int, int, int]	Same physical meaning as :func:generate_spanning_blobs_matrix.	required
seed_start	tuple[int, int, int]	Same physical meaning as :func:generate_spanning_blobs_matrix.	required
max_tries	tuple[int, int, int]	Same physical meaning as :func:generate_spanning_blobs_matrix.	required
show_progress	bool	Retained for compatibility with notebook code. Currently ignored by the packaged implementation.	False
progress_desc	bool	Retained for compatibility with notebook code. Currently ignored by the packaged implementation.	False

Returns:

Type	Description
tuple[ndarray, int]	Same as :func:generate_spanning_blobs_matrix.

Notes

This wrapper exists to reduce migration cost for notebook scripts while preserving deterministic behavior in the core implementation.

Source code in src/voids/generators/porous_image.py

def generate_connected_matrix(
    *,
    shape: tuple[int, int, int],
    porosity: float,
    blobiness: float,
    axis_index: int,
    seed_start: int,
    max_tries: int,
    show_progress: bool = False,
    progress_desc: str | None = None,
) -> tuple[np.ndarray, int]:
    """Backward-compatible wrapper for notebook `08` API signatures.

    Parameters
    ----------
    shape, porosity, blobiness, axis_index, seed_start, max_tries :
        Same physical meaning as :func:`generate_spanning_blobs_matrix`.
    show_progress, progress_desc :
        Retained for compatibility with notebook code. Currently ignored by the
        packaged implementation.

    Returns
    -------
    tuple[numpy.ndarray, int]
        Same as :func:`generate_spanning_blobs_matrix`.

    Notes
    -----
    This wrapper exists to reduce migration cost for notebook scripts while
    preserving deterministic behavior in the core implementation.
    """

    del show_progress, progress_desc
    return generate_spanning_blobs_matrix(
        shape=shape,
        porosity=porosity,
        blobiness=blobiness,
        axis_index=axis_index,
        seed_start=seed_start,
        max_tries=max_tries,
    )

estimate_voronoi_ncells_for_porosity_2d

estimate_voronoi_ncells_for_porosity_2d(
    shape,
    porosity,
    *,
    intercept=0.08,
    slope=0.000322,
    reference_shape=(180, 180),
    min_ncells=40,
)

Estimate Voronoi-cell count needed for a target 2D void fraction.

Parameters:

Name	Type	Description	Default
shape	tuple[int, int]	Output image shape (nx, ny).	required
porosity	float	Target void porosity in (0, 1).	required
intercept	float	Empirical linear model coefficients for phi_void ~= intercept + slope * ncells at reference_shape.	0.08
slope	float	Empirical linear model coefficients for phi_void ~= intercept + slope * ncells at reference_shape.	0.08
reference_shape	tuple[int, int]	Calibration image shape for the linear model.	(180, 180)
min_ncells	int	Lower bound applied to the returned estimate.	40

Returns:

Type	Description
int	Estimated ncells for porespy.generators.voronoi_edges.

Notes

The default calibration was measured in notebook 09 for reference_shape=(180, 180) and r=0 with void = ~voronoi_edges.

Assumptions and caveats

The relation is empirical, not universal. Different resolutions, filtering, or post-processing can shift the mapping significantly.

Source code in src/voids/generators/porous_image.py

def estimate_voronoi_ncells_for_porosity_2d(
    shape: tuple[int, int],
    porosity: float,
    *,
    intercept: float = 0.080,
    slope: float = 3.22e-4,
    reference_shape: tuple[int, int] = (180, 180),
    min_ncells: int = 40,
) -> int:
    """Estimate Voronoi-cell count needed for a target 2D void fraction.

    Parameters
    ----------
    shape :
        Output image shape ``(nx, ny)``.
    porosity :
        Target void porosity in ``(0, 1)``.
    intercept, slope :
        Empirical linear model coefficients for
        ``phi_void ~= intercept + slope * ncells`` at ``reference_shape``.
    reference_shape :
        Calibration image shape for the linear model.
    min_ncells :
        Lower bound applied to the returned estimate.

    Returns
    -------
    int
        Estimated ``ncells`` for ``porespy.generators.voronoi_edges``.

    Notes
    -----
    The default calibration was measured in notebook `09` for
    ``reference_shape=(180, 180)`` and ``r=0`` with ``void = ~voronoi_edges``.

    Assumptions and caveats
    -----------------------
    The relation is empirical, not universal. Different resolutions, filtering,
    or post-processing can shift the mapping significantly.
    """

    dims = normalize_shape(shape, allowed_ndim=(2,))
    if not (0.0 < porosity < 1.0):
        raise ValueError("porosity must be in (0, 1)")
    if slope <= 0:
        raise ValueError("slope must be positive")
    if min_ncells < 1:
        raise ValueError("min_ncells must be >= 1")

    ref = normalize_shape(reference_shape, allowed_ndim=(2,))
    area = float(dims[0] * dims[1])
    area_ref = float(ref[0] * ref[1])
    ncells_ref = (float(porosity) - float(intercept)) / float(slope)
    ncells_scaled = ncells_ref * (area / area_ref)
    return int(max(int(min_ncells), round(ncells_scaled)))

generate_spanning_voronoi_matrix_2d

generate_spanning_voronoi_matrix_2d(
    *,
    shape,
    porosity,
    axis_index,
    seed_start,
    max_tries,
    edge_radius_vox=0,
    target_tol=0.003,
    ncells_step=10,
    search_half_window=70,
    min_ncells=40,
)

Generate a percolating 2D matrix from Voronoi-edge microstructures.

Parameters:

Name	Type	Description	Default
shape	tuple[int, int]	Output image shape (nx, ny).	required
porosity	float	Target void porosity.	required
axis_index	int	Axis used for percolation acceptance.	required
seed_start	int	Seed-search controls.	required
max_tries	int	Seed-search controls.	required
edge_radius_vox	int	Edge thickening radius passed to voronoi_edges(..., r=...).	0
target_tol	float	Relative acceptance tolerance on porosity mismatch.	0.003
ncells_step	int	Search controls around the estimated ncells value.	10
search_half_window	int	Search controls around the estimated ncells value.	10
min_ncells	int	Lower bound for candidate ncells.	40

Returns:

Type	Description
tuple[ndarray, int]	(matrix, seed_used) with boolean matrix encoded as True=void.

Raises:

Type	Description
ValueError	If controls are invalid.
RuntimeError	If no spanning realization is found.

Scientific interpretation

The generator targets low-porosity, edge-dominated connectivity. This is useful for sensitivity studies but remains a synthetic morphology model, not a physically faithful reconstruction of any specific rock sample.

Source code in src/voids/generators/porous_image.py

def generate_spanning_voronoi_matrix_2d(
    *,
    shape: tuple[int, int],
    porosity: float,
    axis_index: int,
    seed_start: int,
    max_tries: int,
    edge_radius_vox: int = 0,
    target_tol: float = 0.003,
    ncells_step: int = 10,
    search_half_window: int = 70,
    min_ncells: int = 40,
) -> tuple[np.ndarray, int]:
    """Generate a percolating 2D matrix from Voronoi-edge microstructures.

    Parameters
    ----------
    shape :
        Output image shape ``(nx, ny)``.
    porosity :
        Target void porosity.
    axis_index :
        Axis used for percolation acceptance.
    seed_start, max_tries :
        Seed-search controls.
    edge_radius_vox :
        Edge thickening radius passed to ``voronoi_edges(..., r=...)``.
    target_tol :
        Relative acceptance tolerance on porosity mismatch.
    ncells_step, search_half_window :
        Search controls around the estimated ``ncells`` value.
    min_ncells :
        Lower bound for candidate ``ncells``.

    Returns
    -------
    tuple[numpy.ndarray, int]
        ``(matrix, seed_used)`` with boolean matrix encoded as ``True=void``.

    Raises
    ------
    ValueError
        If controls are invalid.
    RuntimeError
        If no spanning realization is found.

    Scientific interpretation
    -------------------------
    The generator targets low-porosity, edge-dominated connectivity. This is
    useful for sensitivity studies but remains a synthetic morphology model, not
    a physically faithful reconstruction of any specific rock sample.
    """

    dims = normalize_shape(shape, allowed_ndim=(2,))
    dims2 = (int(dims[0]), int(dims[1]))
    axis = validate_axis_index(axis_index=axis_index, ndim=2)
    if not (0.0 < porosity < 1.0):
        raise ValueError("porosity must be in (0, 1)")
    if max_tries < 1:
        raise ValueError("max_tries must be >= 1")
    if edge_radius_vox < 0:
        raise ValueError("edge_radius_vox must be >= 0")
    if target_tol < 0:
        raise ValueError("target_tol must be >= 0")
    if ncells_step < 1:
        raise ValueError("ncells_step must be >= 1")
    if search_half_window < 0:
        raise ValueError("search_half_window must be >= 0")
    if min_ncells < 1:
        raise ValueError("min_ncells must be >= 1")

    ncells_guess = estimate_voronoi_ncells_for_porosity_2d(
        dims2,
        float(porosity),
        min_ncells=min_ncells,
    )
    ncells_candidates = sorted(
        {
            max(int(min_ncells), ncells_guess + delta)
            for delta in range(-search_half_window, search_half_window + 1, ncells_step)
        }
    )

    best_matrix: np.ndarray | None = None
    best_seed: int | None = None
    best_error = np.inf
    for i in range(max_tries):
        seed = int(seed_start + i)
        for ncells in ncells_candidates:
            matrix = np.asarray(
                ps.generators.voronoi_edges(
                    shape=dims2,
                    ncells=int(ncells),
                    r=int(edge_radius_vox),
                    seed=seed,
                ),
                dtype=bool,
            )
            void = ~matrix
            if not has_spanning_cluster(void, axis_index=axis):
                continue

            err = abs(float(void.mean()) - float(porosity))
            if err < best_error:
                best_matrix = void
                best_seed = seed
                best_error = err
            if err <= target_tol:
                return void, seed

    if best_matrix is not None and best_seed is not None:
        return best_matrix, best_seed

    raise RuntimeError(
        "Could not generate spanning Voronoi matrix "
        f"for seed_start={int(seed_start)} after {int(max_tries)} trials"
    )

generate_spanning_matrix_2d

generate_spanning_matrix_2d(
    *,
    shape,
    porosity,
    axis_index,
    generator_name,
    seed_start,
    max_tries,
    blobs_blobiness=1.8,
    blobs_fallback_porosity_levels=(),
    voronoi_edge_radius_vox=0,
    voronoi_target_tol=0.003,
    voronoi_ncells_step=10,
    voronoi_search_half_window=70,
    voronoi_min_ncells=40,
)

Dispatch 2D matrix generation across supported synthetic families.

Parameters:

Name	Type	Description	Default
shape	tuple[int, int]	Base generation/percolation controls.	required
porosity	tuple[int, int]	Base generation/percolation controls.	required
axis_index	tuple[int, int]	Base generation/percolation controls.	required
seed_start	tuple[int, int]	Base generation/percolation controls.	required
max_tries	tuple[int, int]	Base generation/percolation controls.	required
generator_name	str	Generator family selector: "voronoi_edges" or "blobs".	required
blobs_blobiness	float	Controls for blobs mode, including optional fallback porosity trials if the requested porosity fails to percolate.	1.8
blobs_fallback_porosity_levels	float	Controls for blobs mode, including optional fallback porosity trials if the requested porosity fails to percolate.	1.8
voronoi_edge_radius_vox	int	Edge-thickening radius for Voronoi mode.	0
voronoi_target_tol	float	Porosity mismatch tolerance for Voronoi mode.	0.003
voronoi_ncells_step	int	Candidate-search controls for Voronoi mode.	10
voronoi_search_half_window	int	Candidate-search controls for Voronoi mode.	10
voronoi_min_ncells	int	Candidate-search controls for Voronoi mode.	10

Returns:

Type	Description
tuple[ndarray, int, float]	(matrix, seed_used, porosity_used) where porosity_used is: - requested/fallback porosity parameter for blobs; - achieved void porosity for voronoi_edges.

Raises:

Type	Description
ValueError	If generator_name is unsupported or controls are invalid.
RuntimeError	If a spanning realization cannot be found.

Source code in src/voids/generators/porous_image.py

def generate_spanning_matrix_2d(
    *,
    shape: tuple[int, int],
    porosity: float,
    axis_index: int,
    generator_name: str,
    seed_start: int,
    max_tries: int,
    blobs_blobiness: float = 1.8,
    blobs_fallback_porosity_levels: Sequence[float] = (),
    voronoi_edge_radius_vox: int = 0,
    voronoi_target_tol: float = 0.003,
    voronoi_ncells_step: int = 10,
    voronoi_search_half_window: int = 70,
    voronoi_min_ncells: int = 40,
) -> tuple[np.ndarray, int, float]:
    """Dispatch 2D matrix generation across supported synthetic families.

    Parameters
    ----------
    shape, porosity, axis_index, seed_start, max_tries :
        Base generation/percolation controls.
    generator_name :
        Generator family selector: ``"voronoi_edges"`` or ``"blobs"``.
    blobs_blobiness, blobs_fallback_porosity_levels :
        Controls for ``blobs`` mode, including optional fallback porosity trials
        if the requested porosity fails to percolate.
    voronoi_edge_radius_vox :
        Edge-thickening radius for Voronoi mode.
    voronoi_target_tol :
        Porosity mismatch tolerance for Voronoi mode.
    voronoi_ncells_step, voronoi_search_half_window, voronoi_min_ncells :
        Candidate-search controls for Voronoi mode.

    Returns
    -------
    tuple[numpy.ndarray, int, float]
        ``(matrix, seed_used, porosity_used)`` where ``porosity_used`` is:
        - requested/fallback porosity parameter for ``blobs``;
        - achieved void porosity for ``voronoi_edges``.

    Raises
    ------
    ValueError
        If ``generator_name`` is unsupported or controls are invalid.
    RuntimeError
        If a spanning realization cannot be found.
    """

    generator = str(generator_name).strip().lower()
    if generator == "voronoi_edges":
        matrix, seed = generate_spanning_voronoi_matrix_2d(
            shape=shape,
            porosity=float(porosity),
            axis_index=axis_index,
            seed_start=seed_start,
            max_tries=max_tries,
            edge_radius_vox=voronoi_edge_radius_vox,
            target_tol=voronoi_target_tol,
            ncells_step=voronoi_ncells_step,
            search_half_window=voronoi_search_half_window,
            min_ncells=voronoi_min_ncells,
        )
        return matrix, seed, float(matrix.mean())

    if generator == "blobs":
        if not (0.0 < porosity < 1.0):
            raise ValueError("porosity must be in (0, 1)")
        porosity_trials = [float(porosity)] + [
            float(v) for v in blobs_fallback_porosity_levels if float(v) > float(porosity)
        ]
        porosity_trials = list(dict.fromkeys(porosity_trials))
        last_error: Exception | None = None
        for porosity_try in porosity_trials:
            try:
                matrix, seed = generate_spanning_blobs_matrix(
                    shape=shape,
                    porosity=porosity_try,
                    blobiness=blobs_blobiness,
                    axis_index=axis_index,
                    seed_start=seed_start,
                    max_tries=max_tries,
                )
            except RuntimeError as exc:
                last_error = exc
                continue
            return matrix, seed, float(porosity_try)

        if last_error is not None:
            raise RuntimeError(
                "Could not generate spanning blobs matrix for requested or fallback porosities"
            ) from last_error
        raise RuntimeError("No porosity trial values were available for blobs generation")

    raise ValueError(f"Unsupported generator_name: {generator_name}")

generate_connected_matrix_2d

generate_connected_matrix_2d(
    *,
    shape,
    porosity,
    axis_index,
    generator_name,
    seed_start,
    max_tries,
    show_progress=False,
    progress_desc=None,
)

Backward-compatible wrapper for notebook 09 matrix API.

Parameters:

Name	Type	Description	Default
shape	tuple[int, int]	Same as :func:generate_spanning_matrix_2d.	required
porosity	tuple[int, int]	Same as :func:generate_spanning_matrix_2d.	required
axis_index	tuple[int, int]	Same as :func:generate_spanning_matrix_2d.	required
generator_name	tuple[int, int]	Same as :func:generate_spanning_matrix_2d.	required
seed_start	tuple[int, int]	Same as :func:generate_spanning_matrix_2d.	required
max_tries	tuple[int, int]	Same as :func:generate_spanning_matrix_2d.	required
show_progress	bool	Retained for compatibility. Currently ignored.	False
progress_desc	bool	Retained for compatibility. Currently ignored.	False

Returns:

Type	Description
tuple[ndarray, int, float]	Same as :func:generate_spanning_matrix_2d.

Source code in src/voids/generators/porous_image.py

def generate_connected_matrix_2d(
    *,
    shape: tuple[int, int],
    porosity: float,
    axis_index: int,
    generator_name: str,
    seed_start: int,
    max_tries: int,
    show_progress: bool = False,
    progress_desc: str | None = None,
) -> tuple[np.ndarray, int, float]:
    """Backward-compatible wrapper for notebook `09` matrix API.

    Parameters
    ----------
    shape, porosity, axis_index, generator_name, seed_start, max_tries :
        Same as :func:`generate_spanning_matrix_2d`.
    show_progress, progress_desc :
        Retained for compatibility. Currently ignored.

    Returns
    -------
    tuple[numpy.ndarray, int, float]
        Same as :func:`generate_spanning_matrix_2d`.
    """

    del show_progress, progress_desc
    return generate_spanning_matrix_2d(
        shape=shape,
        porosity=porosity,
        axis_index=axis_index,
        generator_name=generator_name,
        seed_start=seed_start,
        max_tries=max_tries,
    )

insert_ellipsoidal_vug

insert_ellipsoidal_vug(
    matrix_void, *, radii_vox, center=None
)

Insert an axis-aligned ellipsoidal vug into a 3D binary void image.

Parameters:

Name	Type	Description	Default
matrix_void	ndarray	Input 3D boolean array where True denotes void.	required
radii_vox	tuple[int, int, int]	Ellipsoid semi-axes (rx, ry, rz) in voxels.	required
center	tuple[int, int, int] | None	Ellipsoid center index. Defaults to the image center.	None

Returns:

Type	Description
tuple[ndarray, ndarray]	(updated_void, inserted_mask) where inserted_mask identifies the ellipsoidal support.

Raises:

Type	Description
ValueError	If dimensionality or radii are invalid.

Notes

The operation is a boolean union: pre-existing void voxels remain void.

Source code in src/voids/generators/porous_image.py

def insert_ellipsoidal_vug(
    matrix_void: np.ndarray,
    *,
    radii_vox: tuple[int, int, int],
    center: tuple[int, int, int] | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Insert an axis-aligned ellipsoidal vug into a 3D binary void image.

    Parameters
    ----------
    matrix_void :
        Input 3D boolean array where ``True`` denotes void.
    radii_vox :
        Ellipsoid semi-axes ``(rx, ry, rz)`` in voxels.
    center :
        Ellipsoid center index. Defaults to the image center.

    Returns
    -------
    tuple[numpy.ndarray, numpy.ndarray]
        ``(updated_void, inserted_mask)`` where ``inserted_mask`` identifies the
        ellipsoidal support.

    Raises
    ------
    ValueError
        If dimensionality or radii are invalid.

    Notes
    -----
    The operation is a boolean union: pre-existing void voxels remain void.
    """

    arr = np.asarray(matrix_void, dtype=bool)
    if arr.ndim != 3:
        raise ValueError("matrix_void must be a 3D array")

    out = arr.copy()
    nx, ny, nz = out.shape
    if center is None:
        cx, cy, cz = nx // 2, ny // 2, nz // 2
    else:
        if len(center) != 3:
            raise ValueError("center must have length 3")
        cx, cy, cz = center
    rx, ry, rz = (float(radii_vox[0]), float(radii_vox[1]), float(radii_vox[2]))
    if min(rx, ry, rz) <= 0:
        raise ValueError("All ellipsoid radii must be positive")

    x = np.arange(nx, dtype=float) - float(cx)
    y = np.arange(ny, dtype=float) - float(cy)
    z = np.arange(nz, dtype=float) - float(cz)
    xx, yy, zz = np.meshgrid(x, y, z, indexing="ij")
    ellipsoid_mask = (xx / rx) ** 2 + (yy / ry) ** 2 + (zz / rz) ** 2 <= 1.0
    out[ellipsoid_mask] = True
    return out, ellipsoid_mask

insert_spherical_vug

insert_spherical_vug(
    matrix_void, *, radius_vox, center=None
)

Insert a spherical vug into a 3D void image.

Parameters:

Name	Type	Description	Default
matrix_void	ndarray	Input 3D binary void mask.	required
radius_vox	int	Sphere radius in voxels.	required
center	tuple[int, int, int] | None	Optional center index. Defaults to image center.	None

Returns:

Type	Description
tuple[ndarray, ndarray]	Updated void image and spherical support mask.

Source code in src/voids/generators/porous_image.py

def insert_spherical_vug(
    matrix_void: np.ndarray,
    *,
    radius_vox: int,
    center: tuple[int, int, int] | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Insert a spherical vug into a 3D void image.

    Parameters
    ----------
    matrix_void :
        Input 3D binary void mask.
    radius_vox :
        Sphere radius in voxels.
    center :
        Optional center index. Defaults to image center.

    Returns
    -------
    tuple[numpy.ndarray, numpy.ndarray]
        Updated void image and spherical support mask.
    """

    radius = int(radius_vox)
    if radius <= 0:
        raise ValueError("radius_vox must be positive")
    return insert_ellipsoidal_vug(
        matrix_void,
        radii_vox=(radius, radius, radius),
        center=center,
    )

insert_elliptical_vug_2d

insert_elliptical_vug_2d(
    matrix_void, *, radii_vox, center=None
)

Insert an axis-aligned elliptical vug into a 2D binary void image.

Parameters:

Name	Type	Description	Default
matrix_void	ndarray	Input 2D boolean array where True denotes void.	required
radii_vox	tuple[int, int]	Ellipse semi-axes (rx, ry) in voxels.	required
center	tuple[int, int] | None	Optional center index. Defaults to image center.	None

Returns:

Type	Description
tuple[ndarray, ndarray]	Updated void image and ellipse mask.

Source code in src/voids/generators/porous_image.py

def insert_elliptical_vug_2d(
    matrix_void: np.ndarray,
    *,
    radii_vox: tuple[int, int],
    center: tuple[int, int] | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Insert an axis-aligned elliptical vug into a 2D binary void image.

    Parameters
    ----------
    matrix_void :
        Input 2D boolean array where ``True`` denotes void.
    radii_vox :
        Ellipse semi-axes ``(rx, ry)`` in voxels.
    center :
        Optional center index. Defaults to image center.

    Returns
    -------
    tuple[numpy.ndarray, numpy.ndarray]
        Updated void image and ellipse mask.
    """

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

    out = arr.copy()
    nx, ny = out.shape
    if center is None:
        cx, cy = nx // 2, ny // 2
    else:
        if len(center) != 2:
            raise ValueError("center must have length 2")
        cx, cy = center
    rx, ry = (float(radii_vox[0]), float(radii_vox[1]))
    if min(rx, ry) <= 0:
        raise ValueError("All ellipse radii must be positive")

    x = np.arange(nx, dtype=float) - float(cx)
    y = np.arange(ny, dtype=float) - float(cy)
    xx, yy = np.meshgrid(x, y, indexing="ij")
    ellipse_mask = (xx / rx) ** 2 + (yy / ry) ** 2 <= 1.0
    out[ellipse_mask] = True
    return out, ellipse_mask

insert_circular_vug_2d

insert_circular_vug_2d(
    matrix_void, *, radius_vox, center=None
)

Insert a circular vug into a 2D binary void image.

Parameters:

Name	Type	Description	Default
matrix_void	ndarray	Input 2D void mask.	required
radius_vox	int	Circle radius in voxels.	required
center	tuple[int, int] | None	Optional center index.	None

Returns:

Type	Description
tuple[ndarray, ndarray]	Updated image and inserted circular mask.

Source code in src/voids/generators/porous_image.py

def insert_circular_vug_2d(
    matrix_void: np.ndarray,
    *,
    radius_vox: int,
    center: tuple[int, int] | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Insert a circular vug into a 2D binary void image.

    Parameters
    ----------
    matrix_void :
        Input 2D void mask.
    radius_vox :
        Circle radius in voxels.
    center :
        Optional center index.

    Returns
    -------
    tuple[numpy.ndarray, numpy.ndarray]
        Updated image and inserted circular mask.
    """

    radius = int(radius_vox)
    if radius <= 0:
        raise ValueError("radius_vox must be positive")
    return insert_elliptical_vug_2d(
        matrix_void,
        radii_vox=(radius, radius),
        center=center,
    )

make_synthetic_grayscale

make_synthetic_grayscale(
    binary_void,
    *,
    seed,
    void_mean=65.0,
    solid_mean=185.0,
    noise_std=8.0,
    clip_min=0.0,
    clip_max=255.0,
)

Generate synthetic grayscale contrast from binary void/solid phases.

Parameters:

Name	Type	Description	Default
binary_void	ndarray	2D or 3D binary mask where True is void.	required
seed	int	Seed for reproducible random noise.	required
void_mean	float	Mean gray levels assigned to void and solid phases before noise.	65.0
solid_mean	float	Mean gray levels assigned to void and solid phases before noise.	65.0
noise_std	float	Standard deviation of additive Gaussian noise.	8.0
clip_min	float	Output clipping range.	0.0
clip_max	float	Output clipping range.	0.0

Returns:

Type	Description
ndarray	Floating-point grayscale image with same shape as binary_void.

Scientific interpretation

This is a synthetic observation model useful for controlled algorithm benchmarking. It does not represent scanner-specific physics (beam hardening, ring artifacts, partial-volume effects, etc.).

Source code in src/voids/generators/porous_image.py

def make_synthetic_grayscale(
    binary_void: np.ndarray,
    *,
    seed: int,
    void_mean: float = 65.0,
    solid_mean: float = 185.0,
    noise_std: float = 8.0,
    clip_min: float = 0.0,
    clip_max: float = 255.0,
) -> np.ndarray:
    """Generate synthetic grayscale contrast from binary void/solid phases.

    Parameters
    ----------
    binary_void :
        2D or 3D binary mask where ``True`` is void.
    seed :
        Seed for reproducible random noise.
    void_mean, solid_mean :
        Mean gray levels assigned to void and solid phases before noise.
    noise_std :
        Standard deviation of additive Gaussian noise.
    clip_min, clip_max :
        Output clipping range.

    Returns
    -------
    numpy.ndarray
        Floating-point grayscale image with same shape as ``binary_void``.

    Scientific interpretation
    -------------------------
    This is a synthetic observation model useful for controlled algorithm
    benchmarking. It does not represent scanner-specific physics (beam hardening,
    ring artifacts, partial-volume effects, etc.).
    """

    phase = np.asarray(binary_void, dtype=bool)
    if phase.ndim not in {2, 3}:
        raise ValueError("binary_void must be 2D or 3D")
    if noise_std < 0:
        raise ValueError("noise_std must be non-negative")
    if clip_max <= clip_min:
        raise ValueError("clip_max must be larger than clip_min")

    rng = np.random.default_rng(seed)
    base = np.where(phase, float(void_mean), float(solid_mean))
    noise = rng.normal(loc=0.0, scale=float(noise_std), size=phase.shape)
    gray = np.clip(base + noise, float(clip_min), float(clip_max))
    return gray.astype(float)

make_synthetic_grayscale_2d

make_synthetic_grayscale_2d(
    binary_void,
    seed,
    *,
    void_mean=70.0,
    solid_mean=185.0,
    noise_std=8.0,
    clip_min=0.0,
    clip_max=255.0,
)

Backward-compatible 2D wrapper around :func:make_synthetic_grayscale.

Parameters:

Name	Type	Description	Default
binary_void	ndarray	2D binary void mask.	required
seed	int	Random seed for noise realization.	required
void_mean	float	Same meaning as in :func:make_synthetic_grayscale.	70.0
solid_mean	float	Same meaning as in :func:make_synthetic_grayscale.	70.0
noise_std	float	Same meaning as in :func:make_synthetic_grayscale.	70.0
clip_min	float	Same meaning as in :func:make_synthetic_grayscale.	70.0
clip_max	float	Same meaning as in :func:make_synthetic_grayscale.	70.0

Returns:

Type	Description
ndarray	2D floating-point grayscale image.

Source code in src/voids/generators/porous_image.py

def make_synthetic_grayscale_2d(
    binary_void: np.ndarray,
    seed: int,
    *,
    void_mean: float = 70.0,
    solid_mean: float = 185.0,
    noise_std: float = 8.0,
    clip_min: float = 0.0,
    clip_max: float = 255.0,
) -> np.ndarray:
    """Backward-compatible 2D wrapper around :func:`make_synthetic_grayscale`.

    Parameters
    ----------
    binary_void :
        2D binary void mask.
    seed :
        Random seed for noise realization.
    void_mean, solid_mean, noise_std, clip_min, clip_max :
        Same meaning as in :func:`make_synthetic_grayscale`.

    Returns
    -------
    numpy.ndarray
        2D floating-point grayscale image.
    """

    phase = np.asarray(binary_void, dtype=bool)
    if phase.ndim != 2:
        raise ValueError("binary_void must be 2D for make_synthetic_grayscale_2d")
    return make_synthetic_grayscale(
        phase,
        seed=seed,
        void_mean=void_mean,
        solid_mean=solid_mean,
        noise_std=noise_std,
        clip_min=clip_min,
        clip_max=clip_max,
    )

---

Vug Templates

voids.generators.vug_templates

format_radius_token

format_radius_token(value)

Return a stable filename-safe token for radius values.

Source code in src/voids/generators/vug_templates.py

def format_radius_token(value: float) -> str:
    """Return a stable filename-safe token for radius values."""

    return f"{value:.2f}".replace(".", "p")

equivalent_radius_2d

equivalent_radius_2d(radii_xy)

Return the area-equivalent circular radius for an ellipse.

Source code in src/voids/generators/vug_templates.py

def equivalent_radius_2d(radii_xy: tuple[float, float]) -> float:
    """Return the area-equivalent circular radius for an ellipse."""

    rx, ry = radii_xy
    if min(rx, ry) <= 0:
        raise ValueError("Ellipse radii must be positive")
    return float((rx * ry) ** 0.5)

equivalent_radius_3d

equivalent_radius_3d(radii_xyz)

Return the volume-equivalent spherical radius for an ellipsoid.

Source code in src/voids/generators/vug_templates.py

def equivalent_radius_3d(radii_xyz: tuple[float, float, float]) -> float:
    """Return the volume-equivalent spherical radius for an ellipsoid."""

    rx, ry, rz = radii_xyz
    if min(rx, ry, rz) <= 0:
        raise ValueError("Ellipsoid radii must be positive")
    return float((rx * ry * rz) ** (1.0 / 3.0))

match_ellipse_to_circle

match_ellipse_to_circle(
    radius_vox, *, aspect, search_window
)

Match integer ellipse radii (a, b) to circular area r^2.

The optimization prioritizes area error and then aspect-ratio error.

Source code in src/voids/generators/vug_templates.py

def match_ellipse_to_circle(
    radius_vox: int,
    *,
    aspect: float,
    search_window: int,
) -> tuple[int, int]:
    """
    Match integer ellipse radii `(a, b)` to circular area `r^2`.

    The optimization prioritizes area error and then aspect-ratio error.
    """

    _validate_aspect_and_window(aspect=aspect, search_window=search_window)
    r = int(radius_vox)
    if r <= 0:
        raise ValueError("radius_vox must be positive")

    target = float(r**2)
    b_real = r / (aspect**0.5)
    a_real = aspect * b_real
    a0 = int(round(a_real))
    b0 = int(round(b_real))

    best_key: tuple[float, float, float, float] | None = None
    best_tuple: tuple[int, int] | None = None

    b_min = max(1, b0 - search_window)
    b_max = max(b_min, b0 + search_window)
    for b in range(b_min, b_max + 1):
        a_target = target / float(b)
        a_center = int(round(a_target))
        a_min = max(1, min(a0, a_center) - search_window)
        a_max = max(a_min, max(a0, a_center) + search_window)
        for a in range(a_min, a_max + 1):
            area = float(a * b)
            area_err = abs(area - target) / target
            asp = float(a / b)
            asp_err = abs(asp - aspect) / aspect
            center_err = abs(a - a_real) + abs(b - b_real)
            key = (area_err + 0.10 * asp_err, area_err, asp_err, center_err)
            if best_key is None or key < best_key:
                best_key = key
                best_tuple = (int(a), int(b))

    if best_tuple is None:  # pragma: no cover - finite search always yields candidates
        raise RuntimeError("Could not find matched ellipse radii")
    return best_tuple

match_ellipsoid_to_sphere

match_ellipsoid_to_sphere(
    radius_vox, *, aspect, search_window
)

Match integer ellipsoid radii (a, b, b) to spherical volume r^3.

The optimization prioritizes volume error and then aspect-ratio error.

Source code in src/voids/generators/vug_templates.py

def match_ellipsoid_to_sphere(
    radius_vox: int,
    *,
    aspect: float,
    search_window: int,
) -> tuple[int, int, int]:
    """
    Match integer ellipsoid radii `(a, b, b)` to spherical volume `r^3`.

    The optimization prioritizes volume error and then aspect-ratio error.
    """

    _validate_aspect_and_window(aspect=aspect, search_window=search_window)
    r = int(radius_vox)
    if r <= 0:
        raise ValueError("radius_vox must be positive")

    target = float(r**3)
    b_real = r / (aspect ** (1.0 / 3.0))
    a_real = aspect * b_real
    a0 = int(round(a_real))
    b0 = int(round(b_real))

    best_key: tuple[float, float, float, float] | None = None
    best_tuple: tuple[int, int, int] | None = None

    b_min = max(1, b0 - search_window)
    b_max = max(b_min, b0 + search_window)
    for b in range(b_min, b_max + 1):
        a_target = target / float(b * b)
        a_center = int(round(a_target))
        a_min = max(1, min(a0, a_center) - search_window)
        a_max = max(a_min, max(a0, a_center) + search_window)
        for a in range(a_min, a_max + 1):
            vol = float(a * b * b)
            vol_err = abs(vol - target) / target
            asp = float(a / b)
            asp_err = abs(asp - aspect) / aspect
            center_err = abs(a - a_real) + abs(b - b_real)
            key = (vol_err + 0.12 * asp_err, vol_err, asp_err, center_err)
            if best_key is None or key < best_key:
                best_key = key
                best_tuple = (int(a), int(b), int(b))

    if best_tuple is None:  # pragma: no cover - finite search always yields candidates
        raise RuntimeError("Could not find matched ellipsoid radii")
    return best_tuple

build_image_vug_radii_2d

build_image_vug_radii_2d(
    circle_radii_vox, *, aspect, search_window
)

Build area-matched 2D anisotropic radii from circular base radii.

Returns:

Type	Description
tuple	(flow_radii, orth_radii, match_report) where match_report contains (config_index, flow_rel_error, orth_rel_error).

Source code in src/voids/generators/vug_templates.py

def build_image_vug_radii_2d(
    circle_radii_vox: Sequence[int],
    *,
    aspect: float,
    search_window: int,
) -> tuple[list[tuple[int, int]], list[tuple[int, int]], list[tuple[int, float, float]]]:
    """
    Build area-matched 2D anisotropic radii from circular base radii.

    Returns
    -------
    tuple
        ``(flow_radii, orth_radii, match_report)`` where ``match_report`` contains
        ``(config_index, flow_rel_error, orth_rel_error)``.
    """

    flow_radii: list[tuple[int, int]] = []
    orth_radii: list[tuple[int, int]] = []
    report: list[tuple[int, float, float]] = []

    for idx, radius in enumerate(circle_radii_vox, start=1):
        r = int(radius)
        if r <= 0:
            raise ValueError("All circle_radii_vox values must be positive")
        flow = match_ellipse_to_circle(r, aspect=aspect, search_window=search_window)
        orth = (flow[1], flow[0])
        target = float(r**2)
        flow_err = abs(float(flow[0] * flow[1]) - target) / target
        orth_err = abs(float(orth[0] * orth[1]) - target) / target
        flow_radii.append(flow)
        orth_radii.append(orth)
        report.append((idx, flow_err, orth_err))

    return flow_radii, orth_radii, report

build_image_vug_radii_3d

build_image_vug_radii_3d(
    sphere_radii_vox, *, aspect, search_window
)

Build volume-matched 3D anisotropic radii from spherical base radii.

Returns:

Type	Description
tuple	(flow_radii, orth_radii, match_report) where match_report contains (config_index, flow_rel_error, orth_rel_error).

Source code in src/voids/generators/vug_templates.py

def build_image_vug_radii_3d(
    sphere_radii_vox: Sequence[int],
    *,
    aspect: float,
    search_window: int,
) -> tuple[
    list[tuple[int, int, int]],
    list[tuple[int, int, int]],
    list[tuple[int, float, float]],
]:
    """
    Build volume-matched 3D anisotropic radii from spherical base radii.

    Returns
    -------
    tuple
        ``(flow_radii, orth_radii, match_report)`` where ``match_report`` contains
        ``(config_index, flow_rel_error, orth_rel_error)``.
    """

    flow_radii: list[tuple[int, int, int]] = []
    orth_radii: list[tuple[int, int, int]] = []
    report: list[tuple[int, float, float]] = []

    for idx, radius in enumerate(sphere_radii_vox, start=1):
        r = int(radius)
        if r <= 0:
            raise ValueError("All sphere_radii_vox values must be positive")
        flow = match_ellipsoid_to_sphere(r, aspect=aspect, search_window=search_window)
        orth = (flow[1], flow[2], flow[0])
        target = float(r**3)
        flow_err = abs(float(flow[0] * flow[1] * flow[2]) - target) / target
        orth_err = abs(float(orth[0] * orth[1] * orth[2]) - target) / target
        flow_radii.append(flow)
        orth_radii.append(orth)
        report.append((idx, flow_err, orth_err))

    return flow_radii, orth_radii, report

build_lattice_vug_templates_2d

build_lattice_vug_templates_2d(
    *, equiv_radii_spacing, spacing_m, aspect
)

Build 2D lattice vug templates with matched equivalent radius per config.

Returns:

Type	Description
tuple	(templates, match_report), where match_report contains (config_index, circle_err, flow_err, orth_err) relative to target area.

Source code in src/voids/generators/vug_templates.py

def build_lattice_vug_templates_2d(
    *,
    equiv_radii_spacing: Sequence[float],
    spacing_m: float,
    aspect: float,
) -> tuple[list[dict[str, object]], list[tuple[int, float, float, float]]]:
    """
    Build 2D lattice vug templates with matched equivalent radius per config.

    Returns
    -------
    tuple
        ``(templates, match_report)``, where ``match_report`` contains
        ``(config_index, circle_err, flow_err, orth_err)`` relative to target area.
    """

    if spacing_m <= 0:
        raise ValueError("spacing_m must be positive")
    if aspect <= 1.0:
        raise ValueError("aspect must be > 1.0")

    templates: list[dict[str, object]] = []
    report: list[tuple[int, float, float, float]] = []
    aspect_root = aspect**0.5

    for idx, req_over_spacing in enumerate(equiv_radii_spacing, start=1):
        req = float(req_over_spacing)
        if req <= 0:
            raise ValueError("All equiv_radii_spacing values must be positive")
        req_m = req * spacing_m
        target = float(req_m**2)

        circle = (req_m, req_m)
        b = req_m / aspect_root
        a = aspect * b
        flow = (a, b)
        orth = (b, a)

        circle_err = abs(float(circle[0] * circle[1]) - target) / target
        flow_err = abs(float(flow[0] * flow[1]) - target) / target
        orth_err = abs(float(orth[0] * orth[1]) - target) / target

        token = format_radius_token(req)
        templates.extend(
            [
                {
                    "case": f"circle_cfg{idx}_req{token}",
                    "family": "circular",
                    "orientation": "isotropic",
                    "config_index": idx,
                    "radii_xy_m": circle,
                    "r_eq_spacing": req,
                    "target_equivalent_radius_m": req_m,
                    "template_area_rel_error": circle_err,
                },
                {
                    "case": f"ellipse_flow_cfg{idx}_req{token}",
                    "family": "elliptical",
                    "orientation": "flow_stretched",
                    "config_index": idx,
                    "radii_xy_m": flow,
                    "r_eq_spacing": req,
                    "target_equivalent_radius_m": req_m,
                    "template_area_rel_error": flow_err,
                },
                {
                    "case": f"ellipse_orth_cfg{idx}_req{token}",
                    "family": "elliptical",
                    "orientation": "orthogonal_stretched",
                    "config_index": idx,
                    "radii_xy_m": orth,
                    "r_eq_spacing": req,
                    "target_equivalent_radius_m": req_m,
                    "template_area_rel_error": orth_err,
                },
            ]
        )
        report.append((idx, circle_err, flow_err, orth_err))

    return templates, report

build_lattice_vug_templates_3d

build_lattice_vug_templates_3d(
    *, equiv_radii_spacing, spacing_m, aspect
)

Build 3D lattice vug templates with matched equivalent radius per config.

Returns:

Type	Description
tuple	(templates, match_report), where match_report contains (config_index, sphere_err, flow_err, orth_err) relative to target volume.

Source code in src/voids/generators/vug_templates.py

def build_lattice_vug_templates_3d(
    *,
    equiv_radii_spacing: Sequence[float],
    spacing_m: float,
    aspect: float,
) -> tuple[list[dict[str, object]], list[tuple[int, float, float, float]]]:
    """
    Build 3D lattice vug templates with matched equivalent radius per config.

    Returns
    -------
    tuple
        ``(templates, match_report)``, where ``match_report`` contains
        ``(config_index, sphere_err, flow_err, orth_err)`` relative to target
        volume.
    """

    if spacing_m <= 0:
        raise ValueError("spacing_m must be positive")
    if aspect <= 1.0:
        raise ValueError("aspect must be > 1.0")

    templates: list[dict[str, object]] = []
    report: list[tuple[int, float, float, float]] = []
    aspect_root = aspect ** (1.0 / 3.0)

    for idx, req_over_spacing in enumerate(equiv_radii_spacing, start=1):
        req = float(req_over_spacing)
        if req <= 0:
            raise ValueError("All equiv_radii_spacing values must be positive")
        req_m = req * spacing_m
        target = float(req_m**3)

        sphere = (req_m, req_m, req_m)
        b = req_m / aspect_root
        a = aspect * b
        flow = (a, b, b)
        orth = (b, b, a)

        sphere_err = abs(float(sphere[0] * sphere[1] * sphere[2]) - target) / target
        flow_err = abs(float(flow[0] * flow[1] * flow[2]) - target) / target
        orth_err = abs(float(orth[0] * orth[1] * orth[2]) - target) / target

        token = format_radius_token(req)
        templates.extend(
            [
                {
                    "case": f"sphere_cfg{idx}_req{token}",
                    "family": "spherical",
                    "orientation": "isotropic",
                    "config_index": idx,
                    "radii_xyz_m": sphere,
                    "r_eq_spacing": req,
                    "target_equivalent_radius_m": req_m,
                    "template_volume_rel_error": sphere_err,
                },
                {
                    "case": f"ellipsoid_flow_cfg{idx}_req{token}",
                    "family": "ellipsoidal",
                    "orientation": "flow_stretched",
                    "config_index": idx,
                    "radii_xyz_m": flow,
                    "r_eq_spacing": req,
                    "target_equivalent_radius_m": req_m,
                    "template_volume_rel_error": flow_err,
                },
                {
                    "case": f"ellipsoid_orth_cfg{idx}_req{token}",
                    "family": "ellipsoidal",
                    "orientation": "orthogonal_stretched",
                    "config_index": idx,
                    "radii_xyz_m": orth,
                    "r_eq_spacing": req,
                    "target_equivalent_radius_m": req_m,
                    "template_volume_rel_error": orth_err,
                },
            ]
        )
        report.append((idx, sphere_err, flow_err, orth_err))

    return templates, report