diff --git a/src/bonsai/bonsai/tool/cad.py b/src/bonsai/bonsai/tool/cad.py index 957c8339fb..2a1f5f9c69 100644 --- a/src/bonsai/bonsai/tool/cad.py +++ b/src/bonsai/bonsai/tool/cad.py @@ -206,6 +206,237 @@ class Cad: """ return geometry.intersect_line_plane(v1, v2, plane_co, plane_no) + @classmethod + def obb_world_clip_planes( + cls, + center: Vector, + axes: tuple[Vector, Vector, Vector], + half_extents: Vector, + ) -> tuple[tuple[float, float, float, float], ...]: + """Return the 6 inward world clip planes of an oriented bounding box. + + Each plane is a 4-tuple ``(a, b, c, d)`` for the equation + ``a*x + b*y + c*z + d``; a point is KEPT when the value is ``>= 0`` + for every plane, matching ``RegionView3D.clip_planes`` semantics. + Return order is ``(+x, -x, +y, -y, +z, -z)`` where ``+x`` is the face + on the positive side of ``axes[0]``. ``axes`` are assumed orthonormal. + """ + cx, cy, cz = center.x, center.y, center.z + planes: list[tuple[float, float, float, float]] = [] + for i in range(3): + ux, uy, uz = axes[i].x, axes[i].y, axes[i].z + h = float(half_extents[i]) + px, py, pz = cx + h * ux, cy + h * uy, cz + h * uz + nx, ny, nz = -ux, -uy, -uz + planes.append((nx, ny, nz, -(nx * px + ny * py + nz * pz))) + px, py, pz = cx - h * ux, cy - h * uy, cz - h * uz + planes.append((ux, uy, uz, -(ux * px + uy * py + uz * pz))) + return tuple(planes) + + @classmethod + def obb_clip_planes_from_matrix( + cls, + matrix_world: Matrix, + expand: float = 0.0, + expand_rel: float = 0.0, + ) -> tuple[tuple[float, float, float, float], ...]: + """Return the 6 inward world clip planes for the unit cube under ``matrix_world``. + + The implicit box is ``[-1, +1]^3`` in object-local space, so the + host's ``matrix_world`` translation is the world centre, its + rotation orients the box axes, and each column's magnitude is the + world half-extent along that local axis. ``expand`` (absolute + world units) and ``expand_rel`` (fraction of each axis's + half-extent) both add an outward margin — callers that visualise + the box with overlapping geometry (e.g. an empty CUBE display + sharing edges with the clip planes) pass non-zero values so the + box's own wireframe sits safely INSIDE the clip volume. Use the + relative form when the box is rendered at varying scales, since + the depth-buffer precision needed to keep an edge unclipped grows + with world-coordinate magnitude. + """ + world_center = matrix_world.col[3].xyz + linear = matrix_world.to_3x3() + world_axes = [] + world_half_list = [] + for i in range(3): + v = linear.col[i].copy() + length = v.length + if length > 0.0: + world_axes.append(v / length) + else: + world_axes.append(Vector((0.0, 0.0, 0.0))) + world_half_list.append(length + expand + length * expand_rel) + return cls.obb_world_clip_planes( + world_center, + (world_axes[0], world_axes[1], world_axes[2]), + Vector(world_half_list), + ) + + @classmethod + def point_is_inside_clip_planes( + cls, + planes: tuple[tuple[float, float, float, float], ...], + point: Vector, + eps: float = 1e-6, + ) -> bool: + """True iff ``point`` is on the kept side of every plane (inclusive).""" + x, y, z = point.x, point.y, point.z + for a, b, c, d in planes: + if a * x + b * y + c * z + d < -eps: + return False + return True + + @classmethod + def newell_normal(cls, points: Sequence) -> Vector: + """Newell's-method normal for a (possibly non-planar) 3D polygon ring. + + Robust for thin / near-degenerate rings where a two-edge cross + product would be unstable. + """ + nx = ny = nz = 0.0 + n = len(points) + for i in range(n): + cur = points[i] + nxt = points[(i + 1) % n] + nx += (cur[1] - nxt[1]) * (cur[2] + nxt[2]) + ny += (cur[2] - nxt[2]) * (cur[0] + nxt[0]) + nz += (cur[0] - nxt[0]) * (cur[1] + nxt[1]) + return Vector((nx, ny, nz)) + + @classmethod + def plane_basis(cls, points: Sequence) -> tuple[Vector, Vector]: + """Return an orthonormal ``(u, v)`` basis for the ring's best-fit plane.""" + normal = cls.newell_normal(points) + if normal.length < 1e-12: + normal = Vector((0.0, 0.0, 1.0)) + normal = normal.normalized() + ref = Vector((1.0, 0.0, 0.0)) + if abs(normal.x) > 0.9: + ref = Vector((0.0, 1.0, 0.0)) + u = normal.cross(ref) + if u.length < 1e-12: + ref = Vector((0.0, 0.0, 1.0)) + u = normal.cross(ref) + u = u.normalized() + v = normal.cross(u).normalized() + return u, v + + @classmethod + def tessellate_ring_planar(cls, polyline_list: list[list]) -> list[tuple[int, int, int]]: + """Triangulate ``[outer, *inners]`` 3D coord rings in their own plane. + + Projects every ring onto the outer ring's best-fit plane and + returns ``(i, j, k)`` index triples into the flat + ``outer + inners[0] + inners[1] + ...`` vertex list. Falls + back to a shapely constrained Delaunay triangulation when + ``mathutils.geometry.tessellate_polygon`` silently leaves ring + vertices unused (its known failure mode on complex concave + polygons-with-holes). + """ + from mathutils.geometry import tessellate_polygon + + if not polyline_list or not polyline_list[0]: + return [] + outer = polyline_list[0] + u, v = cls.plane_basis(outer) + origin = Vector(outer[0]) + + def _project_xy(ring): + return [((Vector(co) - origin).dot(u), (Vector(co) - origin).dot(v)) for co in ring] + + projected_xy = [_project_xy(ring) for ring in polyline_list] + projected = [[Vector((x, y, 0.0)) for x, y in ring] for ring in projected_xy] + triangles = tessellate_polygon(projected) + + n_total = sum(len(r) for r in projected_xy) + used = {i for tri in triangles for i in tri} + if triangles and len(used) >= n_total: + return triangles + + fallback = cls._tessellate_via_shapely(projected_xy) + return fallback if fallback else triangles + + @classmethod + def _tessellate_via_shapely(cls, projected_xy: list[list[tuple[float, float]]]) -> list[tuple[int, int, int]]: + """Constrained-Delaunay fallback for :meth:`tessellate_ring_planar`. + + Honours the polygon's boundary AND holes. Returns ``[]`` when + shapely is unavailable or the polygon can't be cleaned via + ``buffer(0)``. + """ + try: + from shapely.geometry import Polygon + except Exception: + return [] + outer = projected_xy[0] + inners = projected_xy[1:] + if len(outer) < 3: + return [] + try: + poly = Polygon(outer, inners) + poly = poly if poly.is_valid else poly.buffer(0) + if poly.is_empty: + return [] + except Exception: + return [] + + flat = list(outer) + for r in inners: + flat.extend(r) + + def _key(x, y): + return (round(x, 6), round(y, 6)) + + index_of: dict[tuple[float, float], int] = {} + for idx, (x, y) in enumerate(flat): + index_of.setdefault(_key(x, y), idx) + + try: + from shapely import constrained_delaunay_triangles + + res = constrained_delaunay_triangles(poly) + tri_geoms = list(getattr(res, "geoms", []) or []) + except Exception: + try: + from shapely.ops import triangulate + + tri_geoms = [t for t in triangulate(poly) if poly.contains(t.representative_point())] + except Exception: + return [] + + out: list[tuple[int, int, int]] = [] + for t in tri_geoms: + coords = list(t.exterior.coords)[:-1] + if len(coords) != 3: + continue + idxs = [index_of.get(_key(x, y)) for x, y in coords] + if any(i is None for i in idxs): + continue + out.append(tuple(idxs)) + return out + + @classmethod + def corners_might_cross_clip_planes( + cls, + planes: tuple[tuple[float, float, float, float], ...], + corners: Sequence[Vector], + ) -> bool: + """Conservative reject test: True if ``corners`` might cross the clip volume. + + Returns False only when at least one plane has ALL corners on its + rejected side — meaning the convex hull of ``corners`` is fully + outside the clip volume and a per-mesh bisect can be skipped. + Returns True otherwise (possibly with false positives — never + false negatives), so callers always cap any object that actually + crosses the box. ``corners`` is typically the 8 world-space corners + of an object's bound box. + """ + for a, b, c, d in planes: + if all(a * v.x + b * v.y + c * v.z + d < 0.0 for v in corners): + return False + return True + def intersect_edge_plane_v2(v1, v2, plane_co, plane_no, eps=1e-9): """ Numpy version of intersect_edge_plane diff --git a/src/bonsai/test/tool/test_cad.py b/src/bonsai/test/tool/test_cad.py index 2e84c71718..51491e259d 100644 --- a/src/bonsai/test/tool/test_cad.py +++ b/src/bonsai/test/tool/test_cad.py @@ -16,7 +16,9 @@ # You should have received a copy of the GNU General Public License # along with Bonsai. If not, see . -from mathutils import Vector +import math + +from mathutils import Matrix, Vector from bonsai.tool.cad import Cad as subject from test.bim.bootstrap import NewFile @@ -88,3 +90,87 @@ class TestClosestPoints(NewFile): edge1 = (V(0, 0, 0), V(0, 0, 0)) edge2 = (V(1, 0, 1), V(2, 0, 2)) assert subject.closest_points(edge1, edge2)[0] == (edge1[0], edge2[0]) + + +class TestObbWorldClipPlanes(NewFile): + def test_unit_box_at_origin_returns_axis_aligned_planes(self): + planes = subject.obb_world_clip_planes( + V(0, 0, 0), + (V(1, 0, 0), V(0, 1, 0), V(0, 0, 1)), + V(1, 1, 1), + ) + assert planes[0] == (-1.0, 0.0, 0.0, 1.0) + assert planes[1] == (1.0, 0.0, 0.0, 1.0) + assert planes[2] == (0.0, -1.0, 0.0, 1.0) + assert planes[3] == (0.0, 1.0, 0.0, 1.0) + assert planes[4] == (0.0, 0.0, -1.0, 1.0) + assert planes[5] == (0.0, 0.0, 1.0, 1.0) + + def test_center_is_inside_all_planes(self): + center = V(5, -3, 2) + planes = subject.obb_world_clip_planes( + center, + (V(1, 0, 0), V(0, 1, 0), V(0, 0, 1)), + V(2, 1, 0.5), + ) + assert subject.point_is_inside_clip_planes(planes, center) + + def test_point_just_outside_positive_x_face_rejected(self): + planes = subject.obb_world_clip_planes( + V(0, 0, 0), + (V(1, 0, 0), V(0, 1, 0), V(0, 0, 1)), + V(1, 1, 1), + ) + assert subject.point_is_inside_clip_planes(planes, V(0.5, 0, 0)) + assert not subject.point_is_inside_clip_planes(planes, V(1.5, 0, 0)) + + def test_rotated_obb_clips_along_rotated_axes(self): + s = math.sin(math.radians(45)) + planes = subject.obb_world_clip_planes( + V(0, 0, 0), + (V(s, s, 0), V(-s, s, 0), V(0, 0, 1)), + V(1, 1, 1), + ) + assert subject.point_is_inside_clip_planes(planes, V(1.2, 0, 0)) + assert not subject.point_is_inside_clip_planes(planes, V(1.42, 0, 0)) + + def test_zero_extent_axis_does_not_raise(self): + planes = subject.obb_world_clip_planes( + V(0, 0, 0), + (V(1, 0, 0), V(0, 1, 0), V(0, 0, 1)), + V(1, 1, 0), + ) + assert subject.point_is_inside_clip_planes(planes, V(0, 0, 0)) + + +class TestObbClipPlanesFromMatrix(NewFile): + def test_identity_matches_unit_box(self): + planes = subject.obb_clip_planes_from_matrix(Matrix.Identity(4)) + assert subject.point_is_inside_clip_planes(planes, V(0, 0, 0)) + assert not subject.point_is_inside_clip_planes(planes, V(2, 0, 0)) + assert not subject.point_is_inside_clip_planes(planes, V(0, -2, 0)) + + def test_translated_host_shifts_clip_region(self): + translated = Matrix.Translation(V(10, 0, 0)) + planes = subject.obb_clip_planes_from_matrix(translated) + assert not subject.point_is_inside_clip_planes(planes, V(0, 0, 0)) + assert subject.point_is_inside_clip_planes(planes, V(10, 0, 0)) + + def test_z_rotation_rotates_box(self): + rot = Matrix.Rotation(math.radians(45), 4, "Z") + planes = subject.obb_clip_planes_from_matrix(rot) + assert subject.point_is_inside_clip_planes(planes, V(1.2, 0, 0)) + assert not subject.point_is_inside_clip_planes(planes, V(1.42, 0, 0)) + + def test_host_scale_scales_box_extents(self): + scaled = Matrix.Diagonal((2.0, 2.0, 2.0, 1.0)) + planes = subject.obb_clip_planes_from_matrix(scaled) + assert subject.point_is_inside_clip_planes(planes, V(1.9, 0, 0)) + assert not subject.point_is_inside_clip_planes(planes, V(2.1, 0, 0)) + + def test_non_uniform_scale_axis_independent(self): + scaled = Matrix.Diagonal((3.0, 1.0, 1.0, 1.0)) + planes = subject.obb_clip_planes_from_matrix(scaled) + assert subject.point_is_inside_clip_planes(planes, V(2.9, 0, 0)) + assert not subject.point_is_inside_clip_planes(planes, V(3.1, 0, 0)) + assert not subject.point_is_inside_clip_planes(planes, V(0, 1.1, 0))