From a55753ffb13c9d7e666613b4fd7b9bc4e3896c66 Mon Sep 17 00:00:00 2001 From: "Sean P. Kelly" Date: Tue, 19 Nov 2024 00:03:07 -0800 Subject: [PATCH] feat(pyraydeon): add numpy support --- Cargo.toml | 3 +- pyraydeon/Cargo.toml | 1 + pyraydeon/examples/py_cubes.py | 221 +++++++++++++++++++++++ pyraydeon/examples/py_cubes_expected.svg | 1 + pyraydeon/examples/triangles.py | 4 +- pyraydeon/pyproject.toml | 6 +- pyraydeon/src/linear.rs | 97 +++++++++- pyraydeon/src/ray.rs | 27 ++- pyraydeon/src/scene.rs | 48 ++++- pyraydeon/src/shapes/primitive.rs | 31 +++- 10 files changed, 401 insertions(+), 38 deletions(-) create mode 100644 pyraydeon/examples/py_cubes.py create mode 100644 pyraydeon/examples/py_cubes_expected.svg diff --git a/Cargo.toml b/Cargo.toml index 3a67b45..7a493a3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -16,7 +16,8 @@ env_logger = "0.11" euclid = "0.22" float-cmp = "0.5" log = "0.4" -pyo3 = "0.23" +pyo3 = "0.22" rayon = "1.2" +numpy = "0.22" svg = "0.18" tracing = "0.1" diff --git a/pyraydeon/Cargo.toml b/pyraydeon/Cargo.toml index 59b2e01..aeeee4f 100644 --- a/pyraydeon/Cargo.toml +++ b/pyraydeon/Cargo.toml @@ -11,3 +11,4 @@ crate-type = ["cdylib"] [dependencies] pyo3 = { workspace = true, features = ["extension-module"] } raydeon.workspace = true +numpy.workspace = true diff --git a/pyraydeon/examples/py_cubes.py b/pyraydeon/examples/py_cubes.py new file mode 100644 index 0000000..c1596cf --- /dev/null +++ b/pyraydeon/examples/py_cubes.py @@ -0,0 +1,221 @@ +"""Demonstrates custom object and numpy support + +This example is really slow, and a good case for finding more ways to express +geometry as the sum of native parts. +""" + +import numpy as np +import svg + +from pyraydeon import AABB3, Camera, Geometry, HitData, LineSegment3D, Scene + + +class RectPrism(Geometry): + def __init__( + self, + origin: np.ndarray, + right: np.ndarray, + width: float, + up: np.ndarray, + height: float, + depth: float, + ): + up = up / np.linalg.norm(up) + right = right / np.linalg.norm(right) + fwd = np.cross(up, right) + fwd = fwd / np.linalg.norm(fwd) + + self.origin = origin + self.right = right + self.up = up + self.fwd = fwd + + self.width = width + self.height = height + self.depth = depth + + self.vertices = np.array( + [ + origin, + origin + right * width, + origin + right * width + fwd * depth, + origin + fwd * depth, + origin + up * height, + origin + right * width + up * height, + origin + right * width + fwd * depth + up * height, + origin + fwd * depth + up * height, + ] + ) + + # Make edges stand out slightly so as to not be intersected by their own faces + origin = origin - (right * 0.0015) - (up * 0.0015) - (fwd * 0.0015) + width = width + 0.003 + height = height + 0.003 + depth = depth + 0.003 + self.path_vertices = np.array( + [ + origin, + origin + right * width, + origin + right * width + fwd * depth, + origin + fwd * depth, + origin + up * height, + origin + right * width + up * height, + origin + right * width + fwd * depth + up * height, + origin + fwd * depth + up * height, + ] + ) + + self.faces = [ + [0, 1, 2, 3], + [4, 5, 6, 7], + [0, 1, 5, 4], + [1, 2, 6, 5], + [2, 3, 7, 6], + [3, 0, 4, 7], + ] + self.planes = [self.compute_plane(self.vertices[face]) for face in self.faces] + + def __repr__(self): + return f"RectPrism(basis='[{self.right}, {self.up}, {self.fwd}]', dims='[{self.width}, {self.height}, {self.depth}]')" + + def compute_plane(self, points): + p1, p2, p3 = points[:3] + normal = np.cross(p2 - p1, p3 - p1) + normal /= np.linalg.norm(normal) + d = -np.dot(normal, p1) + return normal, d + + def ray_intersects_plane(self, ray, plane) -> HitData | None: + normal, d = plane + + denom = np.dot(normal, ray.dir) + if abs(denom) < 1e-6: + return None + t = -(np.dot(normal, ray.point) + d) / denom + return HitData(ray.point + t * ray.dir, t) if t >= 0 else None + + def is_point_in_face(self, point, face): + face_vertices = self.vertices[face] + edge1 = face_vertices[1] - face_vertices[0] + edge2 = face_vertices[3] - face_vertices[0] + v = point - face_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: + if not self.bounding_box().hit_by(ray): + return None + for face, plane in zip(self.faces, self.planes): + intersection = self.ray_intersects_plane(ray, plane) + if intersection is not None and self.is_point_in_face( + intersection.hit_point, face + ): + 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) + + 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 + + +up = np.array([-1.0, 1.0, 0.0]) +up = up / np.linalg.norm(up) +right = np.array([1.0, 1.0, 0.0]) +right = right / np.linalg.norm(right) + + +r = RectPrism( + origin=np.array([0.0, 0.0, 0.0]), + right=right, + width=1.0, + up=up, + height=1.0, + depth=1.0, +) + +scene = Scene( + [ + RectPrism( + origin=np.array([0.0, 0.0, 0.0]), + right=right, + width=1.0, + up=up, + height=1.0, + depth=1.0, + ), + RectPrism( + origin=np.array([0.0, 0.0, 1.25]), + right=right, + width=1.0, + up=up, + height=1.0, + depth=1.0, + ), + RectPrism( + origin=right * 1.1, + right=right, + width=1.0, + up=up, + height=1.0, + depth=1.0, + ), + ] +) + +eye = np.array([0.25, 3, 6]) +focus = np.array([0, 0, 0]) +up = np.array([0, 1, 0]) + +fovy = 60.0 +width = 1024 +height = 1024 +znear = 0.1 +zfar = 10.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.x}", + y1=f"{path.p1.y}", + x2=f"{path.p2.x}", + y2=f"{path.p2.y}", + 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) diff --git a/pyraydeon/examples/py_cubes_expected.svg b/pyraydeon/examples/py_cubes_expected.svg new file mode 100644 index 0000000..32d1829 --- /dev/null +++ b/pyraydeon/examples/py_cubes_expected.svg @@ -0,0 +1 @@ + diff --git a/pyraydeon/examples/triangles.py b/pyraydeon/examples/triangles.py index 152877f..7f978c4 100644 --- a/pyraydeon/examples/triangles.py +++ b/pyraydeon/examples/triangles.py @@ -3,7 +3,7 @@ from pyraydeon import Camera, Point3, Scene, Tri, Vec3, Geometry -class CustomObject(Geometry): +class CustomTriangle(Geometry): def __init__(self, p1, p2, p3): self.tri = Tri(p1, p2, p3) @@ -19,7 +19,7 @@ def bounding_box(self): scene = Scene( [ - CustomObject( + CustomTriangle( Point3(0, 0, 0), Point3(0, 0, 1), Point3(1, 0, 1), diff --git a/pyraydeon/pyproject.toml b/pyraydeon/pyproject.toml index 789958a..96f47f8 100644 --- a/pyraydeon/pyproject.toml +++ b/pyraydeon/pyproject.toml @@ -4,9 +4,11 @@ build-backend = "maturin" [project] name = "pyraydeon" -version = "0.1.0-alpha" +version = "0.1.0-alpha3" requires-python = ">=3.9" -dependencies = [] +dependencies = [ + "numpy", +] classifiers = [ "Programming Language :: Rust", "Programming Language :: Python :: Implementation :: CPython", diff --git a/pyraydeon/src/linear.rs b/pyraydeon/src/linear.rs index 846d73f..ff97263 100644 --- a/pyraydeon/src/linear.rs +++ b/pyraydeon/src/linear.rs @@ -1,4 +1,9 @@ +use numpy::PyArray1; +use pyo3::exceptions::PyIndexError; use pyo3::prelude::*; +use raydeon::Shape; + +use crate::ray::{HitData, Ray}; #[derive(Debug, Copy, Clone)] pub struct ArbitrarySpace; @@ -12,10 +17,23 @@ impl Vec3 { raydeon::Vec3::new(x, y, z).into() } + fn __len__(&self) -> usize { + 3 + } + fn as_point(slf: PyRef<'_, Self>) -> PyResult> { Py::new(slf.py(), Point3(slf.0.to_point())) } + fn __getitem__(&self, idx: usize) -> PyResult { + match idx { + 0 => Ok(self.0.x), + 1 => Ok(self.0.y), + 2 => Ok(self.0.z), + _ => Err(PyIndexError::new_err(format!("'{idx}' is out of bounds"))), + } + } + #[getter] fn x(&self) -> f64 { self.0.x @@ -61,6 +79,17 @@ impl Vec3 { } } +impl TryFrom<&Bound<'_, PyAny>> for Vec3 { + type Error = PyErr; + + fn try_from(value: &Bound<'_, PyAny>) -> Result { + let x = value.get_item(0)?.extract()?; + let y = value.get_item(1)?.extract()?; + let z = value.get_item(2)?.extract()?; + Ok(Self::new(x, y, z)) + } +} + pywrap!(Point3, raydeon::Point3); #[pymethods] @@ -70,6 +99,19 @@ impl Point3 { raydeon::Point3::new(x, y, z).into() } + fn __len__(&self) -> usize { + 3 + } + + fn __getitem__(&self, idx: usize) -> PyResult { + match idx { + 0 => Ok(self.0.x), + 1 => Ok(self.0.y), + 2 => Ok(self.0.z), + _ => Err(PyIndexError::new_err(format!("'{idx}' is out of bounds"))), + } + } + #[getter] fn x(&self) -> f64 { self.0.x @@ -102,6 +144,17 @@ impl Point3 { } } +impl TryFrom<&Bound<'_, PyAny>> for Point3 { + type Error = PyErr; + + fn try_from(value: &Bound<'_, PyAny>) -> Result { + let x = value.get_item(0)?.extract()?; + let y = value.get_item(1)?.extract()?; + let z = value.get_item(2)?.extract()?; + Ok(Self::new(x, y, z)) + } +} + pywrap!(Point2, raydeon::Point2); #[pymethods] @@ -121,6 +174,18 @@ impl Point2 { self.0.y } + fn __len__(&self) -> usize { + 2 + } + + fn __getitem__(&self, idx: usize) -> PyResult { + match idx { + 0 => Ok(self.0.x), + 1 => Ok(self.0.y), + _ => Err(PyIndexError::new_err(format!("'{idx}' is out of bounds"))), + } + } + fn __iter__(slf: PyRef<'_, Self>) -> PyResult> { let iter = FloatIter { iter: Box::new([slf.0.x, slf.0.y].into_iter()), @@ -134,6 +199,16 @@ impl Point2 { } } +impl TryFrom<&Bound<'_, PyAny>> for Point2 { + type Error = PyErr; + + fn try_from(value: &Bound<'_, PyAny>) -> Result { + let x = value.get_item(0)?.extract()?; + let y = value.get_item(1)?.extract()?; + Ok(Self::new(x, y)) + } +} + #[pyclass] struct FloatIter { iter: Box + Send + Sync>, @@ -155,18 +230,28 @@ pywrap!(AABB3, raydeon::AABB3); #[pymethods] impl AABB3 { #[new] - fn new(min: &Point3, max: &Point3) -> Self { - raydeon::AABB3::new(min.0, max.0).into() + fn new(min: &Bound<'_, PyAny>, max: &Bound<'_, PyAny>) -> PyResult { + let min = Point3::try_from(min)?; + let max = Point3::try_from(max)?; + Ok(raydeon::AABB3::new(min.0, max.0).into()) } #[getter] - fn min(&self) -> Point3 { - self.0.min.into() + fn min<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let min = [self.0.min.x, self.0.min.y, self.0.min.z]; + PyArray1::from_slice_bound(py, &min) } #[getter] - fn max(&self) -> Point3 { - self.0.max.into() + fn max<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let max = [self.0.max.x, self.0.max.y, self.0.max.z]; + PyArray1::from_slice_bound(py, &max) + } + + fn hit_by(&self, ray: Ray) -> Option { + raydeon::shapes::RectPrism::from(self.0.cast_unit()) + .hit_by(&ray.0) + .map(Into::into) } fn __repr__(slf: &Bound<'_, Self>) -> PyResult { diff --git a/pyraydeon/src/ray.rs b/pyraydeon/src/ray.rs index 5ef9176..49d5a7d 100644 --- a/pyraydeon/src/ray.rs +++ b/pyraydeon/src/ray.rs @@ -1,3 +1,4 @@ +use numpy::PyArray1; use pyo3::prelude::*; use crate::linear::{Point3, Vec3}; @@ -7,18 +8,22 @@ pywrap!(Ray, raydeon::Ray); #[pymethods] impl Ray { #[new] - fn new(point: &Point3, dir: &Vec3) -> Self { - raydeon::Ray::new(point.0.cast_unit(), dir.0.cast_unit()).into() + fn new(point: &Bound<'_, PyAny>, dir: &Bound<'_, PyAny>) -> PyResult { + let point = Point3::try_from(point)?; + let dir = Vec3::try_from(dir)?; + Ok(raydeon::Ray::new(point.0.cast_unit(), dir.0.cast_unit()).into()) } #[getter] - fn point(&self) -> Point3 { - self.0.point.cast_unit().into() + fn point<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let point = [self.0.point.x, self.0.point.y, self.0.point.z]; + PyArray1::from_slice_bound(py, &point) } #[getter] - fn dir(&self) -> Vec3 { - self.0.dir.cast_unit().into() + fn dir<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let dir = [self.0.dir.x, self.0.dir.y, self.0.dir.z]; + PyArray1::from_slice_bound(py, &dir) } fn __repr__(slf: &Bound<'_, Self>) -> PyResult { @@ -32,13 +37,15 @@ pywrap!(HitData, raydeon::HitData); #[pymethods] impl HitData { #[new] - fn new(hit_point: &Point3, dist_to: f64) -> Self { - raydeon::HitData::new(hit_point.0.cast_unit(), dist_to).into() + fn new(hit_point: &Bound<'_, PyAny>, dist_to: f64) -> PyResult { + let hit_point = Point3::try_from(hit_point)?; + Ok(raydeon::HitData::new(hit_point.0.cast_unit(), dist_to).into()) } #[getter] - fn hit_point(&self) -> Point3 { - self.0.hit_point.cast_unit().into() + fn hit_point<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let hp = [self.0.hit_point.x, self.0.hit_point.y, self.0.hit_point.z]; + PyArray1::from_slice_bound(py, &hp) } #[getter] diff --git a/pyraydeon/src/scene.rs b/pyraydeon/src/scene.rs index ba2756e..2627c8c 100644 --- a/pyraydeon/src/scene.rs +++ b/pyraydeon/src/scene.rs @@ -1,5 +1,6 @@ use std::sync::Arc; +use numpy::PyArray1; use pyo3::prelude::*; use raydeon::WorldSpace; @@ -11,8 +12,15 @@ pywrap!(Camera, raydeon::Camera); #[pymethods] impl Camera { #[staticmethod] - fn look_at(eye: &Point3, center: &Vec3, up: &Vec3) -> LookingCamera { - raydeon::Camera::look_at(eye.cast_unit(), center.cast_unit(), up.cast_unit()).into() + fn look_at( + eye: &Bound<'_, PyAny>, + center: &Bound<'_, PyAny>, + up: &Bound<'_, PyAny>, + ) -> PyResult { + let eye = Point3::try_from(eye)?; + let center = Vec3::try_from(center)?; + let up = Vec3::try_from(up)?; + Ok(raydeon::Camera::look_at(eye.cast_unit(), center.cast_unit(), up.cast_unit()).into()) } fn __repr__(slf: &Bound<'_, Self>) -> PyResult { @@ -77,14 +85,23 @@ pywrap!(LineSegment2D, raydeon::path::LineSegment2D); #[pymethods] impl LineSegment2D { + #[new] + fn new(p1: &Bound<'_, PyAny>, p2: &Bound<'_, PyAny>) -> PyResult { + let p1 = Point2::try_from(p1)?; + let p2 = Point2::try_from(p2)?; + Ok(raydeon::path::LineSegment2D::new(p1.cast_unit(), p2.cast_unit()).into()) + } + #[getter] - fn p1(&self) -> Point2 { - self.p1.cast_unit().into() + fn p1<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let p1 = [self.0.p1.x, self.0.p1.y]; + PyArray1::from_slice_bound(py, &p1) } #[getter] - fn p2(&self) -> Point2 { - self.p2.cast_unit().into() + fn p2<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let p2 = [self.0.p2.x, self.0.p2.y]; + PyArray1::from_slice_bound(py, &p2) } fn __repr__(slf: &Bound<'_, Self>) -> PyResult { @@ -97,14 +114,23 @@ pywrap!(LineSegment3D, raydeon::path::LineSegment3D); #[pymethods] impl LineSegment3D { + #[new] + fn new(p1: &Bound<'_, PyAny>, p2: &Bound<'_, PyAny>) -> PyResult { + let p1 = Point3::try_from(p1)?; + let p2 = Point3::try_from(p2)?; + Ok(raydeon::path::LineSegment3D::new(p1.cast_unit(), p2.cast_unit()).into()) + } + #[getter] - fn p1(&self) -> Point3 { - self.p1.cast_unit().into() + fn p1<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let p1 = [self.0.p1.x, self.0.p1.y, self.0.p1.z]; + PyArray1::from_slice_bound(py, &p1) } #[getter] - fn p2(&self) -> Point3 { - self.p2.cast_unit().into() + fn p2<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + let p2 = [self.0.p2.x, self.0.p2.y, self.0.p2.z]; + PyArray1::from_slice_bound(py, &p2) } fn __repr__(slf: &Bound<'_, Self>) -> PyResult { @@ -117,5 +143,7 @@ pub(crate) fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_class::()?; // `LookingCamera` remains "private" m.add_class::()?; + m.add_class::()?; + m.add_class::()?; Ok(()) } diff --git a/pyraydeon/src/shapes/primitive.rs b/pyraydeon/src/shapes/primitive.rs index aea88d1..0778115 100644 --- a/pyraydeon/src/shapes/primitive.rs +++ b/pyraydeon/src/shapes/primitive.rs @@ -1,5 +1,5 @@ use super::Geometry; -use crate::linear::Point3; +use crate::linear::{Point3, Vec3}; use pyo3::prelude::*; use raydeon::WorldSpace; use std::sync::Arc; @@ -25,14 +25,22 @@ impl From> for RectPrism { impl RectPrism { #[new] #[pyo3(signature = (min, max, tag=0))] - fn new(min: &Point3, max: &Point3, tag: usize) -> (Self, Geometry) { + fn new( + min: &Bound<'_, PyAny>, + max: &Bound<'_, PyAny>, + tag: usize, + ) -> PyResult<(Self, Geometry)> { + let min: Vec3 = min.try_into()?; + let max: Vec3 = max.try_into()?; + let shape = Arc::new(raydeon::shapes::RectPrism::tagged( - min.to_vector().cast_unit(), - max.to_vector().cast_unit(), + min.cast_unit(), + max.cast_unit(), tag, )); let geom = Geometry::native(Arc::clone(&shape) as Arc>); - (Self(shape), geom) + + Ok((Self(shape), geom)) } } @@ -57,7 +65,16 @@ impl From> for Tri { impl Tri { #[new] #[pyo3(signature = (p1, p2, p3, tag=0))] - fn new(p1: &Point3, p2: &Point3, p3: &Point3, tag: usize) -> (Self, Geometry) { + fn new( + p1: &Bound<'_, PyAny>, + p2: &Bound<'_, PyAny>, + p3: &Bound<'_, PyAny>, + tag: usize, + ) -> PyResult<(Self, Geometry)> { + let p1: Point3 = p1.try_into()?; + let p2: Point3 = p2.try_into()?; + let p3: Point3 = p3.try_into()?; + let shape = Arc::new(raydeon::shapes::Triangle::tagged( p1.cast_unit(), p2.cast_unit(), @@ -65,6 +82,6 @@ impl Tri { tag, )); let geom = Geometry::native(Arc::clone(&shape) as Arc>); - (Self(shape), geom) + Ok((Self(shape), geom)) } }