fix: fix bug in quad collision detection

Arbitrary quad collision function was using code that only works for
rectangles. This refactors the code to perform intersection checks with
the sum of two triangles.

This also fixes a bug preventing pyraydeon shapes from using native
geometry for collision detection.

Author: cbgbt

Makefile

@@ -10,11 +10,23 @@ check-fmt:
 	cargo fmt --check
 	uv --project pyraydeon run ruff format --check pyraydeon
 
-.PHONY: render-test
-render-test:
+.PHONY: rust-render-test
+rust-render-test:
 	./raydeon/check-examples.sh
+
+.PHONY: reinstall-py-venv
+reinstall-py-venv:
+	echo "Reinstalling native dependencies in virtualenv..."
+	uv --project pyraydeon run --reinstall python -c 'print("Reinstalled dependencies")'
+
+
+.PHONY: py-render-test
+py-render-test: reinstall-py-venv
 	./pyraydeon/check-examples.sh
 
+.PHONY: render-test
+render-test: rust-render-test py-render-test
+
 .PHONY: unit-test
 unit-test:
 	cargo test --locked

pyraydeon/check-examples.sh

@@ -22,9 +22,6 @@ if [ "$ALL_TOOLS_FOUND" = false ]; then
     exit 1
 fi
 
-echo "Reinstalling native dependencies in virtualenv..."
-uv --project ${SCRIPT_DIR} run --reinstall python -c 'print("Reinstalled dependencies")'
-
 for example in ${SCRIPT_DIR}/examples/*.py; do
     example_name=$(basename "$example" .py)
 

pyraydeon/examples/py_cubes.py

@@ -19,36 +19,6 @@
 )
 
 
-class Quad(CollisionGeometry):
-    def __init__(self, vertices):
-        self.vertices = vertices
-        self.plane = self.compute_plane(vertices)
-
-    def compute_plane(self, points):
-        p1, p2, p3 = points[:3]
-        normal = np.cross(p2 - p1, p3 - p1)
-        normal /= np.linalg.norm(normal)
-        return Plane(p1, normal)
-
-    def is_point_in_face(self, point):
-        edge1 = self.vertices[1] - self.vertices[0]
-        edge2 = self.vertices[3] - self.vertices[0]
-        v = point - self.vertices[0]
-        u1 = np.dot(v, edge1) / np.dot(edge1, edge1)
-        u2 = np.dot(v, edge2) / np.dot(edge2, edge2)
-        return 0 <= u1 <= 1 and 0 <= u2 <= 1
-
-    def hit_by(self, ray) -> HitData | None:
-        intersection = self.plane.hit_by(ray)
-        if intersection is not None and self.is_point_in_face(intersection.hit_point):
-            return intersection
-
-    def bounding_box(self):
-        my_min = np.minimum.reduce(self.vertices)
-        my_max = np.maximum.reduce(self.vertices)
-        return AABB3(my_min, my_max)
-
-
 class RectPrism(Geometry):
     def __init__(
         self,
@@ -117,7 +87,7 @@ def __repr__(self):
         return f"RectPrism(basis='[{self.right}, {self.up}, {self.fwd}]', dims='[{self.width}, {self.height}, {self.depth}]')"
 
     def collision_geometry(self):
-        return [Quad(self.vertices[face]) for face in self.faces]
+        return [PyQuad(self.vertices[face]) for face in self.faces]
 
     def paths(self, cam):
         edges = set(
@@ -134,6 +104,36 @@ def paths(self, cam):
         return paths
 
 
+class PyQuad(CollisionGeometry):
+    def __init__(self, vertices):
+        self.vertices = vertices
+        self.plane = self.compute_plane(vertices)
+
+    def compute_plane(self, points):
+        p1, p2, p3 = points[:3]
+        normal = np.cross(p2 - p1, p3 - p1)
+        normal /= np.linalg.norm(normal)
+        return Plane(p1, normal)
+
+    def is_point_in_face(self, point):
+        edge1 = self.vertices[1] - self.vertices[0]
+        edge2 = self.vertices[3] - self.vertices[0]
+        v = point - self.vertices[0]
+        u1 = np.dot(v, edge1) / np.dot(edge1, edge1)
+        u2 = np.dot(v, edge2) / np.dot(edge2, edge2)
+        return 0 <= u1 <= 1 and 0 <= u2 <= 1
+
+    def hit_by(self, ray) -> HitData | None:
+        intersection = self.plane.hit_by(ray)
+        if intersection is not None and self.is_point_in_face(intersection.hit_point):
+            return intersection
+
+    def bounding_box(self):
+        my_min = np.minimum.reduce(self.vertices)
+        my_max = np.maximum.reduce(self.vertices)
+        return AABB3(my_min, my_max)
+
+
 up = np.array([-1.0, 1.0, 0.0])
 up = up / np.linalg.norm(up)
 right = np.array([1.0, 1.0, 0.0])

pyraydeon/examples/py_rhombohedron.py

@@ -0,0 +1,155 @@
+"""Demonstrates custom objects with native collision geometry"""
+
+import numpy as np
+import svg
+
+from pyraydeon import (
+    Camera,
+    Geometry,
+    LineSegment3D,
+    Scene,
+    Quad,
+)
+
+
+class Rhombohedron(Geometry):
+    def __init__(self, origin, basis, dims):
+        basis = basis / np.linalg.norm(basis, axis=1, keepdims=True)
+
+        self.origin = origin
+        self.basis = basis
+        self.dims = dims
+
+        combinations = np.array(np.meshgrid([0, 1], [0, 1], [0, 1])).T.reshape(-1, 3)
+        scaled_combinations = combinations * dims
+
+        # Transform the scaled combinations using the basis
+        transformed_vertices = np.dot(scaled_combinations, basis)
+        # Shift by the origin
+        self.vertices = transformed_vertices + origin
+
+        # Make edges stand out slightly so as to not be intersected by their own faces
+        centroid = np.mean(self.vertices, axis=0)
+        vert_move_dirs = self.vertices - centroid
+        unit_move_dirs = vert_move_dirs / np.linalg.norm(
+            vert_move_dirs, axis=0, keepdims=True
+        )
+        move_vectors = unit_move_dirs * 0.0015
+        self.path_vertices = self.vertices + move_vectors
+
+        self.faces = [
+            [0, 1, 3, 2],  # Bottom face
+            [4, 5, 7, 6],  # Top face
+            [0, 1, 5, 4],  # Front face
+            [2, 3, 7, 6],  # Back face
+            [0, 2, 6, 4],  # Left face
+            [1, 3, 7, 5],  # Right face
+        ]
+        self.quads = self.compute_quads()
+
+    def __repr__(self):
+        return f"Rhomboid(origin='{self.origin}', basis='{self.basis}', dims='{self.dims}')"
+
+    def compute_quads(self):
+        quads = []
+        for face in self.faces:
+            verts = self.vertices[face]
+            origin = verts[0]
+            basis = np.array(
+                [
+                    verts[1] - origin,
+                    verts[3] - origin,
+                ]
+            )
+            dims = np.linalg.norm(basis, axis=1)
+
+            quads.append(Quad(origin, basis, dims))
+        return quads
+
+    def collision_geometry(self):
+        return [geom for quad in self.quads for geom in quad.collision_geometry()]
+
+    def paths(self, cam):
+        edges = set(
+            [
+                tuple(sorted((face[i], face[(i + 1) % len(face)])))
+                for face in self.faces
+                for i in range(len(face))
+            ]
+        )
+        paths = [
+            LineSegment3D(self.path_vertices[edge[0]], self.path_vertices[edge[1]])
+            for edge in edges
+        ]
+        return paths
+
+
+scene = Scene(
+    [
+        Rhombohedron(
+            origin=np.array([0.0, 0.0, 0.0]),
+            basis=np.array(
+                [
+                    [0.9, 0.5, 0.0],
+                    [-0.3, 1.0, 0.0],
+                    [-0.5, 0.25, -0.7],
+                ]
+            ),
+            dims=np.array([1.0, 1.0, 1.0]),
+        ),
+        Rhombohedron(
+            origin=np.array([2.0, 0.0, -2.0]),
+            basis=np.array(
+                [
+                    [-0.9, 0.5, 0.0],
+                    [0.3, 1.0, 0.5],
+                    [1.5, 0.25, -0.7],
+                ]
+            ),
+            dims=np.array([1.0, 1.0, 1.0]),
+        ),
+    ]
+)
+
+eye = np.array([0, 0.4, 5])
+focus = np.array([0, 0.4, 0])
+up = np.array([0, 1, 0])
+
+fovy = 60.0
+width = 1024
+height = 1024
+znear = 0.1
+zfar = 20.0
+
+cam = Camera.look_at(eye, focus, up).perspective(fovy, width, height, znear, zfar)
+
+paths = scene.render(cam)
+
+canvas = svg.SVG(
+    width="8in",
+    height="8in",
+    viewBox="0 0 1024 1024",
+)
+backing_rect = svg.Rect(
+    x=0,
+    y=0,
+    width="100%",
+    height="100%",
+    fill="white",
+)
+svg_lines = [
+    svg.Line(
+        x1=f"{path.p1[0]}",
+        y1=f"{path.p1[1]}",
+        x2=f"{path.p2[0]}",
+        y2=f"{path.p2[1]}",
+        stroke_width="0.7mm",
+        stroke="black",
+    )
+    for path in paths
+]
+line_group = svg.G(transform=f"translate(0, {height}) scale(1, -1)", elements=svg_lines)
+canvas.elements = [backing_rect, line_group]
+
+
+print(canvas)

pyraydeon/src/linear.rs

@@ -1,9 +1,6 @@
-use numpy::{Ix1, PyArray, PyReadonlyArray1};
+use numpy::PyReadonlyArray1;
 use pyo3::exceptions::PyIndexError;
 use pyo3::prelude::*;
-use raydeon::CollisionGeometry as _;
-
-use crate::ray::{HitData, Ray};
 
 #[derive(Debug, Copy, Clone)]
 pub struct ArbitrarySpace;
@@ -248,42 +245,8 @@ impl FloatIter {
     }
 }
 
-pywrap!(AABB3, raydeon::AABB3<ArbitrarySpace>);
-
-#[pymethods]
-impl AABB3 {
-    #[new]
-    fn new(min: PyReadonlyArray1<f64>, max: PyReadonlyArray1<f64>) -> PyResult<Self> {
-        let min = Point3::try_from(min)?;
-        let max = Point3::try_from(max)?;
-        Ok(raydeon::AABB3::new(min.0, max.0).into())
-    }
-
-    #[getter]
-    fn min<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray<f64, Ix1>> {
-        PyArray::from_slice_bound(py, &self.0.min.to_array())
-    }
-
-    #[getter]
-    fn max<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray<f64, Ix1>> {
-        PyArray::from_slice_bound(py, &self.0.max.to_array())
-    }
-
-    fn hit_by(&self, _py: Python, ray: Ray) -> Option<HitData> {
-        raydeon::shapes::AxisAlignedCuboid::from(self.0.cast_unit())
-            .hit_by(&ray.0)
-            .map(Into::into)
-    }
-
-    fn __repr__(slf: &Bound<'_, Self>) -> PyResult<String> {
-        let class_name = slf.get_type().qualname()?;
-        Ok(format!("{}<{:?}>", class_name, slf.borrow().0))
-    }
-}
-
 pub(crate) fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_class::<Vec3>()?;
     m.add_class::<Point3>()?;
-    m.add_class::<AABB3>()?;
     Ok(())
 }

pyraydeon/src/ray.rs

@@ -1,7 +1,8 @@
 use numpy::{Ix1, PyArray, PyReadonlyArray1};
 use pyo3::prelude::*;
+use raydeon::CollisionGeometry as _;
 
-use crate::linear::{Point3, Vec3};
+use crate::linear::{ArbitrarySpace, Point3, Vec3};
 
 pywrap!(Ray, raydeon::Ray);
 
@@ -56,8 +57,42 @@ impl HitData {
     }
 }
 
+pywrap!(AABB3, raydeon::AABB3<ArbitrarySpace>);
+
+#[pymethods]
+impl AABB3 {
+    #[new]
+    fn new(min: PyReadonlyArray1<f64>, max: PyReadonlyArray1<f64>) -> PyResult<Self> {
+        let min = Point3::try_from(min)?;
+        let max = Point3::try_from(max)?;
+        Ok(raydeon::AABB3::new(min.0, max.0).into())
+    }
+
+    #[getter]
+    fn min<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray<f64, Ix1>> {
+        PyArray::from_slice_bound(py, &self.0.min.to_array())
+    }
+
+    #[getter]
+    fn max<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray<f64, Ix1>> {
+        PyArray::from_slice_bound(py, &self.0.max.to_array())
+    }
+
+    fn hit_by(&self, _py: Python, ray: Ray) -> Option<HitData> {
+        raydeon::shapes::AxisAlignedCuboid::from(self.0.cast_unit())
+            .hit_by(&ray.0)
+            .map(Into::into)
+    }
+
+    fn __repr__(slf: &Bound<'_, Self>) -> PyResult<String> {
+        let class_name = slf.get_type().qualname()?;
+        Ok(format!("{}<{:?}>", class_name, slf.borrow().0))
+    }
+}
+
 pub(crate) fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_class::<Ray>()?;
     m.add_class::<HitData>()?;
+    m.add_class::<AABB3>()?;
     Ok(())
 }