"""Routing voxelizations: matrix-backed 3D grids with world/index conversions, morphology and components."""
import warnings as _warnings
from collections.abc import Callable
from functools import cached_property
from itertools import product
from math import ceil
import networkx as nx
import numpy as np
from numpy.typing import NDArray
from scipy.ndimage import binary_dilation as ndimage_binary_dilation
from scipy.ndimage import label
from scipy.spatial import KDTree
from skimage.morphology import (
ball,
binary_closing,
binary_dilation,
binary_erosion,
binary_opening,
)
from volmdlr.core import BoundingBox, Primitive3D
from volmdlr.discrete_representation import DECIMALS, MatrixBasedVoxelization
from volmdlr.model import VolumeModel
from volmdlr.shapes import Solid
import routing.core.functions as fun
from routing.core._types import Bounds, Coord3D, Index3D
from routing.core.functions_compiled import bresenham_3d_c
from routing.core.geometry import Box, RoutingShape, VoxelFrame
from routing.core.verbose import verror
_ALIGNMENT_TOLERANCE = 1e-6 # below this, two grid offsets are the same place
# 13 offsets covering 26-connectivity, one representative per {o, -o} pair.
# Keeps only lexicographically positive offsets, so each undirected edge is visited once.
_OFFSETS_26 = [(i, j, k) for i, j, k in product((-1, 0, 1), repeat=3) if (i, j, k) > (0, 0, 0)]
[docs]
class RoutingVoxelization(MatrixBasedVoxelization): # noqa: PLR0904
"""
A class representing a 3D boolean grid of voxels, placed in space by an origin and a voxel size.
Carries the world/index conversions, the morphology operators and the connected components
that every routing stage needs. ``True`` marks the voxels the grid owns: matter for a
voxelized shape, free space for a carved room. Which of the two it means comes from the
factory that built it, not from the type.
:param matrix: The 3D numpy boolean array where True marks an owned voxel.
:param origin: The centre of voxel (0, 0, 0) in global coordinates.
:param voxel_size: The edge length of a voxel, in the same unit as the primitives.
:param name: The name to identify this voxelization. Defaults to ''.
:param primitives: The CAD view; ``None`` defers it to first access, ``[]`` states there is none.
"""
_non_serializable_attributes = ["voxel_matrix", "matrix", "min_grid_center", "_primitives"] # noqa: RUF012
_MAX_DISPLAYED_VOXELS = 10000
def __init__(
self,
matrix: np.ndarray,
origin: np.ndarray | Coord3D,
voxel_size: float,
name: str = "",
primitives: list[Primitive3D] | None = None,
):
"""Wrap a boolean ``matrix`` anchored at ``origin`` with edge length ``voxel_size``."""
super().__init__(
voxel_matrix=matrix,
min_grid_center=tuple(float(coordinate) for coordinate in np.asarray(origin)),
voxel_size=voxel_size,
name=name,
)
self._primitives = primitives
def __repr__(self) -> str:
"""One-line summary with shape, voxel count and voxel size."""
return f"{type(self).__name__}(shape={self.shape}, voxels={self.voxel_count}, vs={self.voxel_size:.4f})"
[docs]
def copy(self, deep: bool = True, memo: dict | None = None) -> "RoutingVoxelization":
"""
Deep copy of the voxel matrix; ``origin`` and ``voxel_size`` are shared (immutable).
:param deep: If False, perform a shallow copy. If True, perform a deep copy.
:param memo: A dict that keep track of references.
"""
return self.filled_with(self.matrix.copy())
[docs]
@classmethod
def from_routing_shapes(
cls,
routing_shapes: list[RoutingShape],
frame: VoxelFrame,
dilation_voxels: int,
) -> "RoutingVoxelization":
"""Let every shape write itself into one shared ``frame`` grid, each touching only the cells it occupies."""
volume_matrix = np.zeros(frame.shape, dtype=bool)
for routing_shape in routing_shapes:
routing_shape.voxelize_into(volume_matrix, frame, dilation_voxels)
return cls(volume_matrix, frame.origin, frame.voxel_size)
[docs]
@classmethod
def from_voxels(cls, voxels: MatrixBasedVoxelization, build_cad: bool = False) -> "RoutingVoxelization":
"""
Adopt an existing voxelization, keeping its origin, its voxel size, its matter and its name.
:param voxels: The voxelization to adopt.
:param build_cad: Rebuild the CAD view now; otherwise the result states it has none, rather
than deferring an expensive mesh build to the first attribute access.
"""
primitives = voxels.volmdlr_primitives() if build_cad else []
fun.set_primitives_name(primitives, "voxels")
return cls(voxels.matrix, voxels.min_grid_center, voxels.voxel_size, voxels.name, primitives)
[docs]
@staticmethod
def covering_frame(
bounds: np.ndarray,
voxel_size: float,
dilation_voxels: int,
origin: np.ndarray | None = None,
) -> VoxelFrame:
"""
Frame of the grid covering ``bounds``, padded to survive a ``dilation_voxels`` dilation.
:param bounds: World points the grid must contain, as any array of ``(x, y, z)`` rows.
:param voxel_size: The edge length of a voxel.
:param dilation_voxels: The dilation the padding must leave room for.
:param origin: Anchor the grid on an existing one; ``None`` snaps a fresh origin onto whole
voxel steps, so repeated round-trips do not accumulate sub-voxel drift.
"""
corners = np.asarray(bounds, dtype=float).reshape(-1, 3)
padding = (dilation_voxels + 2) * voxel_size
if origin is None:
origin = np.floor(corners.min(axis=0) / voxel_size) * voxel_size - padding
origin = np.asarray(origin, dtype=float)
grid_shape = np.ceil((corners.max(axis=0) + padding - origin) / voxel_size).astype(int) + 1
return VoxelFrame(
tuple(float(coordinate) for coordinate in origin), voxel_size, tuple(int(side) for side in grid_shape)
)
[docs]
@classmethod
def enclosing(
cls,
routing_shapes: list[RoutingShape],
voxel_size: float,
dilation_voxels: int,
) -> "RoutingVoxelization":
"""Voxelize ``routing_shapes`` on a fresh grid auto-sized to enclose them."""
frame = cls.covering_frame(
np.array([shape.bounding_box for shape in routing_shapes]), voxel_size, dilation_voxels
)
return cls.from_routing_shapes(routing_shapes, frame, dilation_voxels)
@property
def primitives(self) -> list[Primitive3D]:
"""CAD primitives reconstructed from voxels, computed on first access."""
if self._primitives is None:
self._primitives = self.volmdlr_primitives()
return self._primitives
@property
def frame(self) -> VoxelFrame:
"""The chart this grid is sampled on."""
return VoxelFrame.of(self)
@property
def origin(self) -> np.ndarray:
"""Grid origin as ndarray (alias for ``min_grid_center``)."""
return np.asarray(self.min_grid_center, dtype=float)
@property
def voxel_origin(self) -> Coord3D:
"""Grid origin as tuple (alias for ``min_grid_center``)."""
return self.min_grid_center
@property
def shape(self) -> tuple:
"""Voxel-grid dimensions (i, j, k)."""
return self.matrix.shape
@property
def voxel_count(self) -> int:
"""Number of occupied voxels."""
return int(self.matrix.sum())
@property
def indices(self) -> np.ndarray:
""":attr:`voxel_indices` transposed to ``(3, n_voxels)``."""
return self.voxel_indices.T
@property
def coords(self) -> np.ndarray:
""":meth:`world_coords` transposed to ``(3, n_voxels)``."""
return self.world_coords().T
@staticmethod
def _homogeneous(indices: np.ndarray) -> np.ndarray:
"""Append a row of ones to ``(3, n_voxels)`` indices."""
return np.concatenate((indices, np.ones((1, indices.shape[1]))), axis=0)
@cached_property
def border_indices(self) -> np.ndarray:
"""Indices of the border voxels — those not fully surrounded by matter, as an ``(N, 3)`` array."""
return np.argwhere(self.border().matrix)
@cached_property
def border_centers(self) -> np.ndarray:
"""World centres of the border voxels, as an ``(N, 3)`` array."""
return self.indices_to_world(self.border_indices)
@cached_property
def homogeneous_indices(self) -> np.ndarray:
""":attr:`indices` with a row of ones appended, as a ``(4, n_voxels)`` array."""
return self._homogeneous(self.indices)
[docs]
def world_coords(self) -> np.ndarray:
"""World coordinates of the occupied voxel centres, as an ``(n_voxels, 3)`` array."""
return self.indices_to_world(self.voxel_indices)
[docs]
def indices_to_world(self, indices: np.ndarray) -> np.ndarray:
"""Convert integer or fractional voxel indices, ``(3,)`` or ``(N, 3)``, to world coordinates."""
return self.frame.centers(indices)
[docs]
def world_to_indices(self, points: np.ndarray) -> np.ndarray:
"""Convert world coordinates, ``(3,)`` or ``(N, 3)``, to the indices of the voxels containing them."""
return self.frame.nearest_indices(points)
[docs]
def world_to_index(self, point: Coord3D | np.ndarray) -> Index3D:
"""Index of the single voxel containing ``point``, as a plain int tuple."""
return self.frame.nearest_index(point)
[docs]
def clip_coord(self, coordinates: np.ndarray) -> np.ndarray:
"""Clip a point onto the closest voxel centre inside the grid."""
clipped = np.clip(self.world_to_indices(coordinates), 0, np.array(self.matrix.shape) - 1)
return self.indices_to_world(clipped)
[docs]
def homogeneous_indices_in_box(self, lower: np.ndarray, upper: np.ndarray) -> np.ndarray:
""":attr:`homogeneous_indices` restricted to the ``[lower, upper)`` index box, computed on that window only."""
lower = np.maximum(lower, 0)
if np.any(np.minimum(upper, self.matrix.shape) <= lower):
return np.ones((4, 0))
window = fun.clipped_grid_window(lower, upper, self.matrix.shape)
return self._homogeneous(np.argwhere(self.matrix[window]).T + lower[:, None])
[docs]
def to_matrix_based_voxelization(self) -> MatrixBasedVoxelization:
"""Downcast to the plain volmdlr voxelization, dropping every routing-side attribute."""
return MatrixBasedVoxelization(
voxel_matrix=self.matrix, voxel_size=self.voxel_size, min_grid_center=self.min_grid_center
)
[docs]
def to_graph(self) -> nx.Graph:
"""26-connectivity graph over the occupied voxels."""
matrix = self.matrix
shape = self.matrix.shape
coordinates = np.argwhere(self.matrix)
graph = nx.Graph()
for coordinate in coordinates:
key = (int(coordinate[0]), int(coordinate[1]), int(coordinate[2]))
graph.add_node(key, voxel_index=np.asarray(coordinate, dtype=float))
for coordinate in coordinates:
i, j, k = int(coordinate[0]), int(coordinate[1]), int(coordinate[2])
for di, dj, dk in _OFFSETS_26:
ni, nj, nk = i + di, j + dj, k + dk
if 0 <= ni < shape[0] and 0 <= nj < shape[1] and 0 <= nk < shape[2] and matrix[ni, nj, nk]:
graph.add_edge((i, j, k), (ni, nj, nk))
return graph
@staticmethod
def _morphed(matrix: np.ndarray, morphology_function: Callable, radius: int) -> np.ndarray:
"""Matrix-level morphology op: returns ``matrix`` unchanged if ``radius <= 0``."""
if radius <= 0:
return matrix
return morphology_function(matrix, ball(radius))
[docs]
def filled_with(self, matrix: np.ndarray) -> "RoutingVoxelization":
"""Get the same grid under the same name, carrying ``matrix`` as its matter."""
return RoutingVoxelization(matrix, self.origin, self.voxel_size, self.name)
def _volume_morphology(self, morphology_function: Callable, radius: int) -> "RoutingVoxelization":
new_matrix = self._morphed(self.matrix, morphology_function, radius)
if new_matrix is self.matrix:
return self
return self.filled_with(new_matrix)
[docs]
def apply_opening(self, radius: int) -> "RoutingVoxelization":
"""Erosion then dilation. Removes small protrusions and noise."""
return self._volume_morphology(binary_opening, radius)
[docs]
def apply_closing(self, radius: int) -> "RoutingVoxelization":
"""Dilation then erosion. Fills small holes and gaps."""
return self._volume_morphology(binary_closing, radius)
[docs]
def apply_dilation(self, radius: int) -> "RoutingVoxelization":
"""Expand all occupied regions outward by ``radius`` voxels."""
return self._volume_morphology(binary_dilation, radius)
[docs]
def apply_erosion(self, radius: int) -> "RoutingVoxelization":
"""Shrink all occupied regions inward by ``radius`` voxels."""
return self._volume_morphology(binary_erosion, radius)
[docs]
def cleaned(
self,
opening_radius: int = 0,
closing_radius: int = 0,
erosion_radius: int = 0,
) -> "RoutingVoxelization":
"""Apply opening, then closing, then erosion; ``self`` unchanged when all radii are 0 (no allocation)."""
matrix = self._morphed(self.matrix, binary_opening, opening_radius)
matrix = self._morphed(matrix, binary_closing, closing_radius)
matrix = self._morphed(matrix, binary_erosion, erosion_radius)
if matrix is self.matrix:
return self
return self.filled_with(matrix)
def _holds(self, point: np.ndarray) -> bool:
"""Whether ``point`` falls inside the grid, on an occupied voxel."""
return self.indices_are_active(self.world_to_indices(point))
[docs]
def indices_are_active(self, indices: list[Index3D] | np.ndarray) -> bool:
"""Whether every index in ``indices`` is in-grid and occupied."""
candidates = np.asarray(indices)
if candidates.size == 0 or np.any(candidates < 0) or np.any(candidates >= self.matrix.shape):
return False
return bool(self.matrix[tuple(candidates.T)].all())
[docs]
def crossed_indices(self, start_position: np.ndarray, end_position: np.ndarray) -> np.ndarray:
"""Return the voxel indices crossed by the ``[start, end]`` segment (Bresenham), dropping out-of-grid ones."""
indices = bresenham_3d_c(self.world_to_index(start_position), self.world_to_index(end_position))
in_grid = np.all((indices >= 0) & (indices < self.matrix.shape), axis=1)
return indices[in_grid]
[docs]
def sample_segment(
self, start_position: np.ndarray, end_position: np.ndarray, samples_per_voxel: float = 2.0
) -> np.ndarray:
"""Sample the ``[start, end]`` segment at ``samples_per_voxel`` density (>= 2 samples, endpoints included)."""
length = float(np.linalg.norm(end_position - start_position))
sample_count = max(2, int(np.ceil(length / self.voxel_size * samples_per_voxel)))
fractions = np.arange(sample_count, dtype=float) / (sample_count - 1)
return start_position + fractions[:, None] * (end_position - start_position)
[docs]
def biggest_component(self) -> "RoutingVoxelization":
"""Get the connected component holding the most voxels."""
parts = self.extract_connected_matrix_based_voxelizations()
return parts[int(np.argmax([len(part) for part in parts]))]
[docs]
def nearest_active_center(self, point: np.ndarray, tolerance_voxels: int | None) -> np.ndarray:
"""
Centre of the occupied voxel nearest to ``point``.
:param point: The world position to search from.
:param tolerance_voxels: Search budget: the nearest centre within a ``±tolerance_voxels``
window, raising beyond it. ``None`` removes the budget — the window grows until the
nearest occupied voxel of the whole grid is found (raising only on an empty grid).
"""
if tolerance_voxels is not None:
return self._nearest_center_within(point, tolerance_voxels)
window_radius = 2
largest_radius = max(self.matrix.shape)
while True:
try:
nearest = self._nearest_center_within(point, window_radius)
except ValueError:
if window_radius >= largest_radius:
raise
window_radius = min(window_radius * 2, largest_radius)
continue
# A corner hit inside the window may still be beaten by an unseen voxel just outside
# it: one pass at the radius certain to contain every closer voxel settles the argmin.
certain_radius = int(np.ceil(float(np.linalg.norm(nearest - point)) / self.voxel_size)) + 1
if certain_radius <= window_radius:
return nearest
return self._nearest_center_within(point, certain_radius)
def _nearest_center_within(self, point: np.ndarray, tolerance_voxels: int) -> np.ndarray:
"""Centre of the occupied voxel nearest to ``point`` within a ``±tolerance_voxels`` window; raise when none."""
index = np.array(self.world_to_index(point))
window_low = np.clip(index - tolerance_voxels, 0, self.matrix.shape)
window = fun.clipped_grid_window(window_low, index + tolerance_voxels + 1, self.matrix.shape)
window_indices = np.argwhere(self.matrix[window]) + window_low
if window_indices.size == 0:
raise ValueError(f"Point {point.tolist()} is farther than {tolerance_voxels} voxels from any active voxel.")
centers = self.indices_to_world(window_indices)
return centers[np.linalg.norm(centers - point, axis=1).argmin()]
[docs]
def component_containing(self, point: np.ndarray, tolerance_voxels: int = 2) -> "RoutingVoxelization":
"""
Return the connected component carrying the occupied voxel nearest to ``point``.
:param point: The world position to seed the search from.
:param tolerance_voxels: How far ``point`` may penetrate the carved material (voxel-quantization
slack); beyond that it is genuinely unroutable and this raises rather than selecting a wrong domain.
"""
point = np.asarray(point, dtype=float)
seed = self.nearest_active_center(point, tolerance_voxels)
for part in self.extract_connected_matrix_based_voxelizations():
if part._holds(seed):
return part
raise RuntimeError(f"Active voxel at {seed.tolist()} missing from every extracted component.")
[docs]
def component_carrying(self, point: np.ndarray) -> "RoutingVoxelization":
"""Get the connected component occupied at ``point`` exactly; the whole voxelization, logged, when none is.
The strict sibling of :meth:`component_containing`: no quantization slack, and a point on no
component falls back to the whole voxelization instead of raising.
"""
for part in self.extract_connected_matrix_based_voxelizations():
index = np.array(part.world_to_index(point))
if np.all(index >= 0) and np.all(index < part.matrix.shape) and bool(part.matrix[tuple(index)]):
return part
verror("Chosen point of Pipe's path not in voxels. Taking whole voxelization.")
return self
[docs]
def offset_in(self, parent: "RoutingVoxelization") -> np.ndarray:
"""Integer index offset of this grid's origin inside ``parent``."""
return np.round((self.origin - np.asarray(parent.min_grid_center)) / self.voxel_size).astype(int)
[docs]
def nesting_in(self, coarser: "RoutingVoxelization") -> tuple[int, np.ndarray]:
"""
Get the two numbers mapping a ``coarser`` voxel onto the block of this grid's voxels it covers.
Together they are the conversion: the block covering coarse voxel ``(i, j, k)`` starts at
``offset + (i, j, k) * ratio`` here and spans ``ratio`` voxels per axis. The multi-scale
counterpart of :meth:`offset_in`, which needs only the offset because both grids share a
voxel size.
A coarser voxel must cover a whole block and sit centred on it, so ratio and offset both
have to land on integers; anything else means the two grids were not built from one another
and the caller is about to read the wrong voxels.
:param coarser: The grid this one nests into, of a larger voxel size.
:return: The size ratio, and the index here at which ``coarser``'s first voxel starts.
"""
ratio = coarser.voxel_size / self.voxel_size
integer_ratio = round(ratio)
if abs(ratio - integer_ratio) > _ALIGNMENT_TOLERANCE:
raise ValueError(f"Voxel size ratio must be an integer, got {ratio}")
centring = (integer_ratio - 1) / 2 * self.voxel_size
misalignment = coarser.origin - (self.origin + centring)
offset = np.round(misalignment / self.voxel_size).astype(int)
if np.any(np.abs(misalignment - offset * self.voxel_size) > _ALIGNMENT_TOLERANCE):
raise ValueError("Voxelizations are not properly aligned")
return integer_ratio, offset
[docs]
def restrict_field(self, parent: "RoutingVoxelization", field: NDArray[np.float64]) -> NDArray[np.float64]:
"""
Sample a field defined on ``parent``'s grid over this sub-grid: parent values on occupied voxels, ``inf`` elsewhere.
``self`` must lie on ``parent``'s grid — the relation :meth:`intersection` with ``parent`` produces (same
voxel size, integer offset). Occupied voxels of ``self`` are then a subset of ``parent``'s, so each maps to a
finite ``field`` entry; the ``inf`` fill only reaches unoccupied voxels, which carry no value by convention.
"""
offset = self.offset_in(parent)
return np.where(self.matrix, field[fun.grid_window(offset, offset + self.matrix.shape)], np.inf)
[docs]
def set_distances_direct(self, distances_array: np.ndarray) -> np.ndarray:
"""Scatter one distance per occupied voxel into a full grid, ``inf`` on the unoccupied ones."""
result = np.full(self.matrix.shape, np.inf)
result[self.indices[0], self.indices[1], self.indices[2]] = distances_array
return result
# RELATIONS BETWEEN TWO VOXELIZATIONS
def _contact_distance(self, other: "RoutingVoxelization") -> float:
"""Distance under which a voxel centre of each grid counts as touching the other."""
return 0.5 * (self.voxel_size + other.voxel_size) + 1e-6
[docs]
def is_adjacent_to(self, other: "RoutingVoxelization") -> bool:
"""Whether the two voxelizations touch: some pair of border voxels lies within one voxel step."""
if not self.bounding_box.is_intersecting(other.bounding_box):
return False
positions_1, positions_2 = self.border_centers, other.border_centers
if len(positions_1) == 0 or len(positions_2) == 0:
return False
distance_threshold = self._contact_distance(other)
mask_1, mask_2 = self._proximity_masks(positions_1, positions_2, distance_threshold)
candidates_1, candidates_2 = positions_1[mask_1], positions_2[mask_2]
if len(candidates_1) == 0 or len(candidates_2) == 0:
return False
if len(candidates_2) < len(candidates_1):
candidates_1, candidates_2 = candidates_2, candidates_1
tree = KDTree(candidates_2)
return any(tree.query_ball_point(position, distance_threshold) for position in candidates_1)
@staticmethod
def _proximity_masks(
positions_1: np.ndarray,
positions_2: np.ndarray,
distance_threshold: float,
margin_factor: float = 1.01,
) -> tuple[np.ndarray, np.ndarray]:
"""
Masks keeping, on each side, the positions inside the other side's bounding box widened by the threshold.
A cheap pre-filter: anything outside that box cannot be within ``distance_threshold`` of the
other set, so the KD-tree only ever sees the plausible candidates.
"""
if positions_1.size == 0 or positions_2.size == 0:
return np.zeros(len(positions_1), dtype=bool), np.zeros(len(positions_2), dtype=bool)
margin = distance_threshold * margin_factor
box_1_min, box_1_max = positions_1.min(axis=0) - margin, positions_1.max(axis=0) + margin
box_2_min, box_2_max = positions_2.min(axis=0) - margin, positions_2.max(axis=0) + margin
return (
np.all((positions_1 >= box_2_min) & (positions_1 <= box_2_max), axis=1),
np.all((positions_2 >= box_1_min) & (positions_2 <= box_1_max), axis=1),
)
@staticmethod
def _position_pairs_within(
positions_a: np.ndarray, positions_b: np.ndarray, distance_threshold: float
) -> np.ndarray:
"""
Every pair of positions, one from each set, closer than ``distance_threshold``.
:return: Array of shape ``(K, 6)`` where each row is ``[x1, y1, z1, x2, y2, z2]``.
"""
tree = KDTree(positions_b)
pairs = tree.query_ball_point(positions_a, distance_threshold)
query_indices = []
target_indices = []
for query_index, nearby_indices in enumerate(pairs):
if nearby_indices:
query_indices.extend([query_index] * len(nearby_indices))
target_indices.extend(nearby_indices)
if not query_indices:
return np.empty((0, 6), dtype=np.float64)
return np.hstack([positions_a[np.array(query_indices)], positions_b[np.array(target_indices)]])
[docs]
def labelled_components(self) -> tuple[np.ndarray, int]:
"""Label each voxel with the index of its connected component, ``-1`` off matter, and count the components.
The 6-connectivity and the numbering of :meth:`extract_connected_matrix_based_voxelizations`,
without materialising one voxelization per component.
"""
labels, component_count = label(self.matrix)
return labels.astype(np.int32) - 1, int(component_count)
[docs]
def touching_border_pairs(self, other: "RoutingVoxelization") -> np.ndarray:
"""
Pairs of border voxels of the two grids lying within one voxel step of each other.
The predicate of :meth:`is_adjacent_to`, evaluated for every pair at once: grids whose bounding
boxes are apart touch nowhere, and otherwise each side's border is first reduced to the shell the
other side's border can reach (:meth:`_border_facing`), so the distance search only ever sees the
contact, whatever the size of the grids.
:return: Array of shape ``(K, 2)``: rows of ``self.border_indices`` and of ``other.border_indices``.
"""
if not self.bounding_box.is_intersecting(other.bounding_box):
return np.empty((0, 2), dtype=int)
threshold = self._contact_distance(other)
coarse, fine = (self, other) if self.voxel_size >= other.voxel_size else (other, self)
coarse_rows = np.flatnonzero(coarse._border_facing(fine.border_centers, threshold))
if not coarse_rows.size:
return np.empty((0, 2), dtype=int)
fine_rows = np.flatnonzero(fine._border_facing(coarse.border_centers[coarse_rows], threshold))
if not fine_rows.size:
return np.empty((0, 2), dtype=int)
pairs = KDTree(coarse.border_centers[coarse_rows]).sparse_distance_matrix(
KDTree(fine.border_centers[fine_rows]), threshold, output_type="ndarray"
)
coarse_to_fine = np.stack([coarse_rows[pairs["i"]], fine_rows[pairs["j"]]], axis=1)
return coarse_to_fine if coarse is self else coarse_to_fine[:, ::-1]
def _border_facing(self, centres: np.ndarray, threshold: float) -> np.ndarray:
"""Mask over :attr:`border_indices` of the border voxels that may lie within ``threshold`` of one of ``centres``.
The centres are rasterised on this grid and their cells widened by the number of steps
``threshold`` spans, so no pair is lost and the distance search after it stays local.
"""
reach = int(np.floor(threshold / self.voxel_size + 0.5))
shape = np.asarray(self.shape)
padded_shape = shape + 2 * reach
cells = self.frame.nearest_indices(centres) + reach
cells = cells[np.all((cells >= 0) & (cells < padded_shape), axis=1)]
marked = np.unique(np.ravel_multi_index(tuple(cells.T), padded_shape))
window = tuple(slice(reach, reach + side) for side in self.shape)
if len(marked) * (2 * reach + 1) ** 3 > reach * int(np.prod(padded_shape)):
# Dense marks on a small grid: dilating the grid costs less than widening each mark.
marks = np.zeros(padded_shape, dtype=bool)
marks[np.unravel_index(marked, padded_shape)] = True
marks = ndimage_binary_dilation(marks, structure=np.ones((3, 3, 3), dtype=bool), iterations=reach)
return marks[window][tuple(self.border_indices.T)]
cells = np.stack(np.unravel_index(marked, padded_shape), axis=1) - reach
reachable = np.zeros(self.shape, dtype=bool)
for offset in product(range(-reach, reach + 1), repeat=3):
shifted = cells + offset
inside = np.all((shifted >= 0) & (shifted < shape), axis=1)
reachable[tuple(shifted[inside].T)] = True
return reachable[tuple(self.border_indices.T)]
[docs]
def subtract_coarser(self, coarser: "RoutingVoxelization") -> "RoutingVoxelization":
"""
Clear every voxel covered by an occupied voxel of ``coarser`` — a multi-scale boolean difference.
Each occupied coarse voxel clears the whole block of this grid's voxels it covers,
following the block mapping of :meth:`nesting_in`.
"""
ratio, offset = self.nesting_in(coarser)
result_matrix = self.matrix.copy()
for coarse_index in np.argwhere(coarser.matrix):
block_start = np.maximum(offset + coarse_index * ratio, 0)
block_end = np.minimum(offset + (coarse_index + 1) * ratio, self.matrix.shape)
result_matrix[fun.grid_window(block_start, block_end)] = False
return RoutingVoxelization(result_matrix, self.min_grid_center, self.voxel_size)
[docs]
def component_labels(self, components: list["RoutingVoxelization"]) -> np.ndarray:
"""
Label each voxel of this grid with the index of the component covering it, ``-1`` elsewhere.
:param components: Voxelizations lying on this grid, typically its connected components.
"""
labels = np.full(self.matrix.shape, -1, dtype=np.int32)
for component_index, component in enumerate(components):
offset = component.offset_in(self)
window = fun.grid_window(offset, offset + np.array(component.matrix.shape))
labels[window][component.matrix] = component_index
return labels
# INDEX BOXES
[docs]
def maximal_boxes(self, min_fill: float = 1.0) -> list[Bounds]:
"""
Greedy cover of the occupied voxels by maximal index boxes, each box carved out as it is taken.
Every seed is grown along the six axis orders and the largest box wins, so a cluttered room
yields fewer, fatter boxes — which is the point: a :class:`SmartGrid` allocates its node
container over the whole box, so a box that is all matter wastes nothing.
:param min_fill: Smallest share of free matter in every slab a box grows by: 1 keeps the boxes all
matter; below, a box takes in some empty space and the room splits into fewer boxes. A box never
takes a cell of another, empty or not, so the boxes never overlap.
"""
orders = [("x", "y", "z"), ("x", "z", "y"), ("y", "x", "z"), ("y", "z", "x"), ("z", "x", "y"), ("z", "y", "x")]
remaining = self.matrix.copy()
unclaimed = np.ones(self.matrix.shape, dtype=bool)
boxes = []
while np.any(remaining):
seed = np.argwhere(remaining)[0]
candidates = [self._grow_box_from_seed(remaining, unclaimed, seed, order, min_fill) for order in orders]
best = max(candidates, key=self._box_volume)
remaining[fun.box_window(best)] = False
unclaimed[fun.box_window(best)] = False
boxes.append(best)
return boxes
@staticmethod
def _grow_box_from_seed(
remaining: np.ndarray, unclaimed: np.ndarray, seed: np.ndarray, order: tuple[str, str, str], min_fill: float
) -> Bounds:
"""Grow the maximal box around ``seed``, expanding one axis at a time in ``order`` (:meth:`maximal_boxes`)."""
seed_x, seed_y, seed_z = (int(coordinate) for coordinate in seed)
bounds: Bounds = {"x": (seed_x, seed_x), "y": (seed_y, seed_y), "z": (seed_z, seed_z)}
for axis_name in order:
bounds[axis_name] = RoutingVoxelization._expand_box_along_axis(
remaining, unclaimed, bounds, "xyz".index(axis_name), min_fill
)
return bounds
@staticmethod
def _expand_box_along_axis(
remaining: np.ndarray, unclaimed: np.ndarray, bounds: Bounds, axis: int, min_fill: float
) -> tuple[int, int]:
"""Expand the box along ``axis`` while every swept slab is unclaimed and ``min_fill`` free matter, the other axes fixed."""
current_min, current_max = bounds["xyz"[axis]]
while current_max + 1 < remaining.shape[axis]:
slab = RoutingVoxelization._slab(bounds, axis, current_max + 1)
if not RoutingVoxelization._slab_joins(remaining, unclaimed, slab, min_fill):
break
current_max += 1
while current_min - 1 >= 0:
slab = RoutingVoxelization._slab(bounds, axis, current_min - 1)
if not RoutingVoxelization._slab_joins(remaining, unclaimed, slab, min_fill):
break
current_min -= 1
return current_min, current_max
@staticmethod
def _slab(bounds: Bounds, axis: int, index: int) -> tuple[slice, slice, slice]:
"""Index window of the one-cell-thick slab of the box at ``index`` along ``axis``."""
window = list(fun.box_window(bounds))
window[axis] = slice(index, index + 1)
return tuple(window)
@staticmethod
def _slab_joins(
remaining: np.ndarray, unclaimed: np.ndarray, slab: tuple[slice, slice, slice], min_fill: float
) -> bool:
"""Whether a box may grow by ``slab``: no cell of it taken by another box, at least ``min_fill`` free matter."""
return bool(np.all(unclaimed[slab])) and remaining[slab].mean() >= min_fill
@staticmethod
def _box_volume(bounds: Bounds) -> int:
"""Count the voxels the inclusive box spans."""
return (
(bounds["x"][1] - bounds["x"][0] + 1)
* (bounds["y"][1] - bounds["y"][0] + 1)
* (bounds["z"][1] - bounds["z"][0] + 1)
)
[docs]
def sub_voxelization(self, bounds: Bounds, copy: bool = True) -> "RoutingVoxelization":
"""
Extract the voxelization spanning the inclusive index box ``bounds``, re-anchored on its own origin.
:param bounds: Per-axis inclusive ``(lowest, highest)`` indices.
:param copy: Copy the extracted block. ``False`` shares memory with this grid's matrix —
cheaper for a short-lived window, but it keeps the whole parent matrix alive.
"""
window = fun.box_window(bounds)
sub_matrix = self.matrix[window].copy() if copy else self.matrix[window]
lowest_corner = np.array([bounds[axis][0] for axis in "xyz"])
sub_origin = np.round(self.origin + lowest_corner * self.voxel_size, DECIMALS)
return RoutingVoxelization(sub_matrix, sub_origin, self.voxel_size)
[docs]
def refined_box(self, bounds: Bounds, factor: int) -> "RoutingVoxelization":
"""
Build a fully-occupied voxelization, ``factor`` times finer, covering the inclusive index box ``bounds``.
:param bounds: Per-axis inclusive ``(lowest, highest)`` indices on this grid.
:param factor: How many fine voxels a voxel of this grid spans per axis.
"""
box_offset = np.array([bounds[axis][0] for axis in "xyz"])
box_origin = self.min_grid_center + box_offset * self.voxel_size
box_shape = np.array([bounds[axis][1] - bounds[axis][0] + 1 for axis in "xyz"])
refined_voxel_size = round(self.voxel_size / factor, 6)
refined_origin = np.round(box_origin - (((factor - 1) / 2) * refined_voxel_size), 6)
return RoutingVoxelization(
matrix=np.ones(box_shape * factor, dtype=bool),
origin=tuple(refined_origin.tolist()),
voxel_size=refined_voxel_size,
)
[docs]
def within_bounding_box(self, bounding_box: BoundingBox) -> "RoutingVoxelization | None":
"""Keep the occupied voxels whose centres fall inside ``bounding_box``, cropped to them; ``None`` when empty."""
if len(self) == 0:
return None
centers = self.get_voxel_centers_array()
inside = (
(centers[:, 0] >= bounding_box.xmin)
& (centers[:, 0] <= bounding_box.xmax)
& (centers[:, 1] >= bounding_box.ymin)
& (centers[:, 1] <= bounding_box.ymax)
& (centers[:, 2] >= bounding_box.zmin)
& (centers[:, 2] <= bounding_box.zmax)
)
inside_indices = self.voxel_indices[inside]
if len(inside_indices) == 0:
return None
filtered_matrix = np.zeros_like(self.matrix)
filtered_matrix[tuple(inside_indices.T)] = True
return RoutingVoxelization(filtered_matrix, self.min_grid_center, self.voxel_size).crop_matrix()
[docs]
def voxel_boxes(
self,
centers: np.ndarray | None = None,
scale: float = 1.0,
color: tuple[float, float, float] = (0.5, 0.5, 0.5),
alpha: float = 0.8,
decimation: int = 1,
) -> list[Solid]:
"""
Render voxels as cube solids of side ``scale * voxel_size`` — the inspection view of the grid.
Distinct from the inherited :meth:`volmdlr_primitives`, which returns one surface mesh for
the whole volume: here each voxel stays a cube, so the discretization itself is visible.
:param centers: World centres to draw; ``None`` draws this volume's own occupied voxels.
:param scale: Cube side as a multiple of the voxel size, to make markers stand out.
:param decimation: Keep one centre every ``decimation``, to stay under a viewable count.
"""
drawn = self.world_coords() if centers is None else np.asarray(centers, dtype=float).reshape(-1, 3)
return Box.cubes_at(drawn[::decimation], scale * self.voxel_size, color, alpha)
[docs]
def display(
self,
title: str = "",
overlay_prims: list[Primitive3D] | None = None,
decimation: int = 0,
color: tuple[float, float, float] = (0.5, 0.5, 0.5),
alpha: float = 0.4,
) -> None:
"""Visualize this volume as small box primitives, with optional overlay."""
if not decimation:
decimation = max(ceil(self.voxel_count / self._MAX_DISPLAYED_VOXELS), 1)
if title:
print(f"\n--- {title} [{self.voxel_count} voxels] ---")
boxes = self.voxel_boxes(color=color, alpha=alpha, decimation=decimation)
print(f" {self.voxel_count} voxels -> {len(boxes)} boxes (decimation={decimation})")
overlay = list(overlay_prims) if overlay_prims else []
VolumeModel(overlay + boxes).display_3d()
[docs]
class VoxelVolume(RoutingVoxelization):
"""Deprecated: use :class:`RoutingVoxelization`, which merges this class with ``RoomVoxelization``."""
def __init__(self, *args, **kwargs):
"""Emit a DeprecationWarning and forward the parent's argument layout untouched."""
_warnings.warn("VoxelVolume is deprecated, use RoutingVoxelization instead.", DeprecationWarning, stacklevel=2)
super().__init__(*args, **kwargs)
[docs]
class RoomVoxelization(RoutingVoxelization):
"""
Deprecated: use :class:`RoutingVoxelization`, which merges this class with ``VoxelVolume``.
Construction shim only. This class keeps its historical argument order, which no longer matches
the positional layout ``MatrixBasedVoxelization`` rebuilds instances with, so operations that
reconstruct — ``extract_connected_matrix_based_voxelizations``, ``crop_matrix``, the boolean
operators — must be called on a :class:`RoutingVoxelization` instead.
"""
def __init__(
self,
primitives: list[Primitive3D],
voxel_matrix: NDArray[np.bool_],
voxel_size: float,
voxel_origin: Coord3D,
name: str = "",
):
"""Emit a DeprecationWarning and forward, keeping this class's historical argument order."""
_warnings.warn(
"RoomVoxelization is deprecated, use RoutingVoxelization instead.", DeprecationWarning, stacklevel=2
)
super().__init__(voxel_matrix, voxel_origin, voxel_size, name, primitives)
