"""Constructive Solid Geometry (CSG) for boolean mesh operations.
Provides CSG shape nodes that generate meshes and a combiner node that performs
boolean operations (union, subtract, intersect) on child CSG shapes using a
BSP-tree algorithm.
Usage::
from simvx.core.csg import CSGBox3D, CSGSphere3D, CSGCombiner3D, CSGOperation
combiner = CSGCombiner3D(operation=CSGOperation.SUBTRACT)
combiner.add_child(CSGBox3D(size=(2, 2, 2)))
combiner.add_child(CSGSphere3D(radius=1.2))
mesh = combiner.get_mesh() # Box with sphere carved out
"""
from __future__ import annotations
from __future__ import annotations
import logging
from collections.abc import Sequence
from enum import IntEnum
import numpy as np
from .descriptors import Property
from .graphics.mesh import Mesh
from .math.types import Vec3
from .nodes_3d.node3d import Node3D
log = logging.getLogger(__name__)
__all__ = [
"CSGOperation",
"CSGShape3D",
"CSGBox3D",
"CSGSphere3D",
"CSGCylinder3D",
"CSGCombiner3D",
]
# Tolerance for plane-side classification
_EPS = 1e-5
[docs]
class CSGOperation(IntEnum):
"""Boolean operation type for CSG combining."""
UNION = 0
SUBTRACT = 1
INTERSECT = 2
# ============================================================================
# BSP Polygon / Plane / Tree
# ============================================================================
class _Polygon:
"""A convex polygon with vertices (Nx3 float32) and a plane normal/w."""
__slots__ = ("vertices", "normal", "w")
def __init__(self, vertices: np.ndarray, normal: np.ndarray | None = None, w: float | None = None) -> None:
self.vertices = np.asarray(vertices, dtype=np.float32)
if normal is not None and w is not None:
self.normal = np.asarray(normal, dtype=np.float32)
self.w = float(w)
else:
self._compute_plane()
def _compute_plane(self) -> None:
"""Compute plane normal and distance from first three vertices."""
v = self.vertices
if len(v) < 3:
self.normal = np.array([0, 0, 1], dtype=np.float32)
self.w = 0.0
return
a = v[1] - v[0]
b = v[2] - v[0]
n = np.cross(a, b)
length = np.linalg.norm(n)
if length < 1e-10:
self.normal = np.array([0, 0, 1], dtype=np.float32)
self.w = 0.0
else:
self.normal = n / length
self.w = float(np.dot(self.normal, v[0]))
def flip(self) -> _Polygon:
"""Return a new polygon with reversed winding and flipped normal."""
return _Polygon(self.vertices[::-1].copy(), -self.normal.copy(), -self.w)
def clone(self) -> _Polygon:
"""Return a deep copy."""
return _Polygon(self.vertices.copy(), self.normal.copy(), self.w)
# Side classification constants
_COPLANAR = 0
_FRONT = 1
_BACK = 2
_SPANNING = 3
def _classify_polygon(plane_normal: np.ndarray, plane_w: float, poly: _Polygon) -> int:
"""Classify a polygon relative to a splitting plane."""
dots = poly.vertices @ plane_normal - plane_w
front = np.any(dots > _EPS)
back = np.any(dots < -_EPS)
if front and back:
return _SPANNING
if front:
return _FRONT
if back:
return _BACK
return _COPLANAR
def _split_polygon(
plane_normal: np.ndarray,
plane_w: float,
poly: _Polygon,
coplanar_front: list[_Polygon],
coplanar_back: list[_Polygon],
front: list[_Polygon],
back: list[_Polygon],
) -> None:
"""Split a polygon by a plane, distributing pieces to the four output lists."""
dots = poly.vertices @ plane_normal - plane_w
types = np.where(dots > _EPS, _FRONT, np.where(dots < -_EPS, _BACK, _COPLANAR))
poly_type = 0
for t in types:
poly_type |= t
if poly_type == _COPLANAR:
if np.dot(plane_normal, poly.normal) > 0:
coplanar_front.append(poly)
else:
coplanar_back.append(poly)
elif poly_type == _FRONT:
front.append(poly)
elif poly_type == _BACK:
back.append(poly)
else:
# Spanning — split the polygon
f_verts: list[np.ndarray] = []
b_verts: list[np.ndarray] = []
n = len(poly.vertices)
for i in range(n):
j = (i + 1) % n
ti, tj = int(types[i]), int(types[j])
vi, vj = poly.vertices[i], poly.vertices[j]
if ti != _BACK:
f_verts.append(vi)
if ti != _FRONT:
b_verts.append(vi)
if (ti | tj) == _SPANNING:
# Compute intersection point
denom = np.dot(plane_normal, vj - vi)
if abs(denom) > 1e-10:
t_val = (plane_w - np.dot(plane_normal, vi)) / denom
t_val = max(0.0, min(1.0, t_val))
v_interp = vi + t_val * (vj - vi)
else:
v_interp = vi.copy()
f_verts.append(v_interp)
b_verts.append(v_interp.copy())
if len(f_verts) >= 3:
front.append(_Polygon(np.array(f_verts, dtype=np.float32), poly.normal.copy(), poly.w))
if len(b_verts) >= 3:
back.append(_Polygon(np.array(b_verts, dtype=np.float32), poly.normal.copy(), poly.w))
class _BSPNode:
"""BSP tree node. Each node holds a splitting plane and polygon lists."""
__slots__ = ("plane_normal", "plane_w", "polygons", "front", "back")
def __init__(self, polygons: list[_Polygon] | None = None) -> None:
self.plane_normal: np.ndarray | None = None
self.plane_w: float = 0.0
self.polygons: list[_Polygon] = []
self.front: _BSPNode | None = None
self.back: _BSPNode | None = None
if polygons:
self.build(polygons)
def clone(self) -> _BSPNode:
node = _BSPNode()
if self.plane_normal is not None:
node.plane_normal = self.plane_normal.copy()
node.plane_w = self.plane_w
node.polygons = [p.clone() for p in self.polygons]
node.front = self.front.clone() if self.front else None
node.back = self.back.clone() if self.back else None
return node
def invert(self) -> None:
"""Flip all polygons and swap front/back subtrees (inside-out)."""
self.polygons = [p.flip() for p in self.polygons]
if self.plane_normal is not None:
self.plane_normal = -self.plane_normal
self.plane_w = -self.plane_w
if self.front:
self.front.invert()
if self.back:
self.back.invert()
self.front, self.back = self.back, self.front
def clip_polygons(self, polygons: list[_Polygon]) -> list[_Polygon]:
"""Recursively clip polygons against this BSP tree."""
if self.plane_normal is None:
return list(polygons)
front_list: list[_Polygon] = []
back_list: list[_Polygon] = []
for p in polygons:
_split_polygon(self.plane_normal, self.plane_w, p, front_list, back_list, front_list, back_list)
if self.front:
front_list = self.front.clip_polygons(front_list)
if self.back:
back_list = self.back.clip_polygons(back_list)
else:
back_list = []
return front_list + back_list
def clip_to(self, bsp: _BSPNode) -> None:
"""Clip all polygons in this tree against another BSP tree."""
self.polygons = bsp.clip_polygons(self.polygons)
if self.front:
self.front.clip_to(bsp)
if self.back:
self.back.clip_to(bsp)
def all_polygons(self) -> list[_Polygon]:
"""Collect all polygons from this tree."""
result = list(self.polygons)
if self.front:
result.extend(self.front.all_polygons())
if self.back:
result.extend(self.back.all_polygons())
return result
def build(self, polygons: list[_Polygon]) -> None:
"""Build the BSP tree from a list of polygons (iterative to avoid stack overflow)."""
if not polygons:
return
# Use an iterative worklist to avoid deep recursion on large meshes
stack: list[tuple[_BSPNode, list[_Polygon]]] = [(self, polygons)]
while stack:
node, polys = stack.pop()
if not polys:
continue
if node.plane_normal is None:
node.plane_normal = polys[0].normal.copy()
node.plane_w = polys[0].w
front_list: list[_Polygon] = []
back_list: list[_Polygon] = []
for p in polys:
_split_polygon(node.plane_normal, node.plane_w, p, node.polygons, node.polygons, front_list, back_list)
if front_list:
if node.front is None:
node.front = _BSPNode()
stack.append((node.front, front_list))
if back_list:
if node.back is None:
node.back = _BSPNode()
stack.append((node.back, back_list))
# ============================================================================
# BSP Boolean Operations
# ============================================================================
def _csg_union(a: _BSPNode, b: _BSPNode) -> list[_Polygon]:
"""Compute union of two BSP trees, return polygon list."""
a = a.clone()
b = b.clone()
a.clip_to(b)
b.clip_to(a)
b.invert()
b.clip_to(a)
b.invert()
a.build(b.all_polygons())
return a.all_polygons()
def _csg_subtract(a: _BSPNode, b: _BSPNode) -> list[_Polygon]:
"""Compute A - B (subtract B from A), return polygon list."""
a = a.clone()
b = b.clone()
a.invert()
a.clip_to(b)
b.clip_to(a)
b.invert()
b.clip_to(a)
b.invert()
a.build(b.all_polygons())
a.invert()
return a.all_polygons()
def _csg_intersect(a: _BSPNode, b: _BSPNode) -> list[_Polygon]:
"""Compute intersection of two BSP trees, return polygon list."""
a = a.clone()
b = b.clone()
a.invert()
b.clip_to(a)
b.invert()
a.clip_to(b)
b.clip_to(a)
a.build(b.all_polygons())
a.invert()
return a.all_polygons()
# ============================================================================
# Mesh <-> Polygon conversion
# ============================================================================
def _mesh_to_polygons(mesh: Mesh) -> list[_Polygon]:
"""Convert a Mesh into a list of triangle _Polygons.
Uses stored mesh normals when available to determine correct face orientation,
since winding order may differ from the normal direction (e.g. Vulkan CW convention).
"""
positions = mesh.positions
has_normals = mesh.normals is not None and len(mesh.normals) == len(positions)
if mesh.indices is not None:
tris = mesh.indices.reshape(-1, 3)
else:
tris = np.arange(len(positions), dtype=np.uint32).reshape(-1, 3)
polygons: list[_Polygon] = []
for tri in tris:
verts = positions[tri].copy()
edge1 = verts[1] - verts[0]
edge2 = verts[2] - verts[0]
cross_n = np.cross(edge1, edge2)
cross_len = np.linalg.norm(cross_n)
if cross_len < 1e-10:
continue # Skip degenerate triangles
cross_n /= cross_len
if has_normals:
# Use the mesh's stored normal as the authoritative direction
mesh_normal = mesh.normals[tri[0]].copy()
mesh_normal_len = np.linalg.norm(mesh_normal)
if mesh_normal_len > 1e-10:
mesh_normal /= mesh_normal_len
# If winding-derived normal disagrees with mesh normal, flip winding
if np.dot(cross_n, mesh_normal) < 0:
verts = verts[::-1].copy()
cross_n = -cross_n
w = float(np.dot(cross_n, verts[0]))
polygons.append(_Polygon(verts, cross_n, w))
return polygons
def _polygons_to_mesh(polygons: list[_Polygon]) -> Mesh:
"""Convert a list of _Polygons back into a Mesh with normals."""
if not polygons:
# Return a minimal empty-ish mesh (single degenerate triangle)
positions = np.zeros((3, 3), dtype=np.float32)
normals = np.zeros((3, 3), dtype=np.float32)
normals[:, 2] = 1.0
indices = np.array([0, 1, 2], dtype=np.uint32)
return Mesh(positions, indices, normals)
all_positions: list[np.ndarray] = []
all_normals: list[np.ndarray] = []
all_indices: list[int] = []
vert_offset = 0
for poly in polygons:
n_verts = len(poly.vertices)
if n_verts < 3:
continue
# Fan triangulation from first vertex (CCW winding, right-hand rule)
for i in range(1, n_verts - 1):
all_indices.extend([vert_offset, vert_offset + i, vert_offset + i + 1])
all_positions.append(poly.vertices)
all_normals.append(np.tile(poly.normal, (n_verts, 1)))
vert_offset += n_verts
if not all_positions:
positions = np.zeros((3, 3), dtype=np.float32)
normals = np.zeros((3, 3), dtype=np.float32)
normals[:, 2] = 1.0
return Mesh(positions, np.array([0, 1, 2], dtype=np.uint32), normals)
positions = np.vstack(all_positions).astype(np.float32)
normals = np.vstack(all_normals).astype(np.float32)
indices = np.array(all_indices, dtype=np.uint32)
return Mesh(positions, indices, normals)
[docs]
def csg_combine(mesh_a: Mesh, mesh_b: Mesh, operation: CSGOperation) -> Mesh:
"""Perform a CSG boolean operation on two meshes.
Args:
mesh_a: First operand mesh.
mesh_b: Second operand mesh.
operation: The boolean operation to apply.
Returns:
A new ``Mesh`` with the result of the boolean operation.
"""
polys_a = _mesh_to_polygons(mesh_a)
polys_b = _mesh_to_polygons(mesh_b)
if not polys_a or not polys_b:
return mesh_a if polys_a else mesh_b
bsp_a = _BSPNode(polys_a)
bsp_b = _BSPNode(polys_b)
if operation == CSGOperation.UNION:
result_polys = _csg_union(bsp_a, bsp_b)
elif operation == CSGOperation.SUBTRACT:
result_polys = _csg_subtract(bsp_a, bsp_b)
elif operation == CSGOperation.INTERSECT:
result_polys = _csg_intersect(bsp_a, bsp_b)
else:
raise ValueError(f"Unknown CSG operation: {operation}")
return _polygons_to_mesh(result_polys)
# ============================================================================
# CSG Shape Nodes
# ============================================================================
[docs]
class CSGShape3D(Node3D):
"""Base class for CSG shape nodes.
Holds a computed mesh and an operation type used when combined by a
``CSGCombiner3D`` parent. Subclasses override ``_build_mesh()`` to
generate their primitive geometry.
"""
operation = Property(CSGOperation.UNION, enum=[CSGOperation.UNION, CSGOperation.SUBTRACT, CSGOperation.INTERSECT],
hint="Boolean operation when combined with siblings", group="CSG")
def __init__(self, operation: CSGOperation = CSGOperation.UNION, **kwargs):
super().__init__(**kwargs)
self.operation = operation
self._cached_mesh: Mesh | None = None
self._dirty: bool = True
def _build_mesh(self) -> Mesh:
"""Generate the primitive mesh. Override in subclasses."""
raise NotImplementedError
def _invalidate_csg(self) -> None:
"""Mark the cached mesh as stale."""
self._dirty = True
self._cached_mesh = None
[docs]
def get_mesh(self) -> Mesh:
"""Return the (possibly cached) mesh for this shape."""
if self._dirty or self._cached_mesh is None:
self._cached_mesh = self._build_mesh()
self._dirty = False
return self._cached_mesh
[docs]
class CSGBox3D(CSGShape3D):
"""CSG box primitive. Generates an axis-aligned box mesh.
Attributes:
size: Box dimensions as ``Vec3`` (width, height, depth).
"""
size = Property(Vec3(1, 1, 1), hint="Box dimensions (width, height, depth)", group="CSG")
def __init__(self, size: Vec3 | Sequence[float] | None = None, **kwargs):
super().__init__(**kwargs)
if size is not None:
self.size = Vec3(size)
self._invalidate_csg()
def _build_mesh(self) -> Mesh:
from .surface_tool import create_box
return create_box(self.size)
[docs]
class CSGSphere3D(CSGShape3D):
"""CSG sphere primitive. Generates a UV sphere mesh.
Attributes:
radius: Sphere radius.
rings: Number of horizontal rings (latitude divisions).
sectors: Number of vertical sectors (longitude divisions).
"""
radius = Property(1.0, range=(0.001, 100.0), hint="Sphere radius", group="CSG")
rings = Property(16, range=(3, 128), hint="Horizontal ring count", group="CSG")
sectors = Property(16, range=(3, 128), hint="Vertical sector count", group="CSG")
def __init__(
self,
radius: float = 1.0,
rings: int = 16,
sectors: int = 16,
**kwargs,
):
super().__init__(**kwargs)
self.radius = radius
self.rings = rings
self.sectors = sectors
self._invalidate_csg()
def _build_mesh(self) -> Mesh:
from .surface_tool import create_sphere
return create_sphere(self.radius, self.rings, self.sectors)
[docs]
class CSGCylinder3D(CSGShape3D):
"""CSG cylinder primitive. Generates a cylinder mesh along the Y axis.
Attributes:
radius: Cylinder radius.
height: Total height.
segments: Number of radial segments.
"""
radius = Property(0.5, range=(0.001, 100.0), hint="Cylinder radius", group="CSG")
height = Property(1.0, range=(0.001, 100.0), hint="Cylinder height", group="CSG")
segments = Property(16, range=(3, 128), hint="Radial segment count", group="CSG")
def __init__(
self,
radius: float = 0.5,
height: float = 1.0,
segments: int = 16,
**kwargs,
):
super().__init__(**kwargs)
self.radius = radius
self.height = height
self.segments = segments
self._invalidate_csg()
def _build_mesh(self) -> Mesh:
from .surface_tool import create_cylinder
return create_cylinder(self.radius, self.height, self.segments)
[docs]
class CSGCombiner3D(CSGShape3D):
"""Combines child CSG shapes using boolean operations.
The first child is the base shape. Each subsequent child is combined with
the accumulated result using that child's ``operation`` property.
Usage::
combiner = CSGCombiner3D()
combiner.add_child(CSGBox3D(size=(2, 2, 2)))
combiner.add_child(CSGSphere3D(radius=1.2, operation=CSGOperation.SUBTRACT))
result = combiner.get_mesh()
"""
def __init__(self, **kwargs):
super().__init__(**kwargs)
[docs]
def get_mesh(self) -> Mesh:
"""Combine all child CSG shapes and return the resulting mesh."""
if not self._dirty and self._cached_mesh is not None:
return self._cached_mesh
csg_children = [c for c in self.children if isinstance(c, CSGShape3D)]
if not csg_children:
# No children — return a degenerate mesh
self._cached_mesh = _polygons_to_mesh([])
self._dirty = False
return self._cached_mesh
result_mesh = csg_children[0].get_mesh()
for child in csg_children[1:]:
child_mesh = child.get_mesh()
result_mesh = csg_combine(result_mesh, child_mesh, child.operation)
self._cached_mesh = result_mesh
self._dirty = False
return self._cached_mesh
def _build_mesh(self) -> Mesh:
"""Not used directly — get_mesh() handles combining."""
return self.get_mesh()