mirror of
https://github.com/IfcOpenShell/IfcOpenShell.git
synced 2026-08-06 07:51:47 +00:00
Add OBB clip-plane and planar tessellation to tool.Cad
Adds geometry primitives the viewport clip-box feature needs: - obb_world_clip_planes / obb_clip_planes_from_matrix: derive the 6 inward clip planes of an oriented bounding box (or unit cube under a matrix_world) in RegionView3D.clip_planes form. expand / expand_rel margins let callers visualising the box with overlapping geometry (an empty CUBE display sharing edges with the planes) keep the box's own wireframe inside the clip volume. - point_is_inside_clip_planes / corners_might_cross_clip_planes: cheap reject tests for the per-mesh capping pass to skip the expensive bisect when an object's AABB is fully outside the box. - newell_normal / plane_basis: robust planar-ring normal for thin near-degenerate cap rings where a two-edge cross product is unstable. - tessellate_ring_planar: triangulate [outer, *inners] 3D rings in the outer ring's best-fit plane, with a shapely constrained-Delaunay fallback for the known failure mode of mathutils.tessellate_polygon on complex concave polygons-with-holes. Tests cover unit-box, translated, rotated, and scaled cases for the OBB-from-matrix builder + the rejection helpers. Generated with the assistance of an AI coding tool.
This commit is contained in:
@@ -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
|
||||
|
||||
@@ -16,7 +16,9 @@
|
||||
# You should have received a copy of the GNU General Public License
|
||||
# along with Bonsai. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
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))
|
||||
|
||||
Reference in New Issue
Block a user