From f5cfb3b7815d219518a0a29d7f48fa97221dc30d Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Mon, 18 Dec 2023 16:52:18 +0000 Subject: [PATCH 01/17] Private method for SDF internals --- python/ncollpyde/_ncollpyde.pyi | 3 ++ python/ncollpyde/main.py | 39 ++++++++++++++++----- src/interface.rs | 62 +++++++++++++++++++++++++++++++++ src/utils.rs | 13 +++++-- 4 files changed, 106 insertions(+), 11 deletions(-) diff --git a/python/ncollpyde/_ncollpyde.pyi b/python/ncollpyde/_ncollpyde.pyi index f94a936..1d60b2c 100644 --- a/python/ncollpyde/_ncollpyde.pyi +++ b/python/ncollpyde/_ncollpyde.pyi @@ -30,3 +30,6 @@ class TriMeshWrapper: def intersections_many_threaded( self, src_points: Points, tgt_points: Points ) -> Tuple[List[int], Points, List[bool]]: ... + def sdf_intersections( + self, points: Points, vectors: Points + ) -> Tuple[npt.NDArray[np.float_], npt.NDArray[np.float_]]: ... diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index 327ceea..ed2da22 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -2,7 +2,7 @@ import random import warnings from multiprocessing import cpu_count -from typing import TYPE_CHECKING, Optional, Tuple, Union +from typing import TYPE_CHECKING, Optional, Tuple, Union, List import numpy as np from numpy.typing import ArrayLike, NDArray @@ -218,6 +218,34 @@ def contains( return self._impl.contains(coords, self._interpret_threads(threads)) + def _as_points(self, points: ArrayLike) -> NDArray: + p = np.asarray(points, self.dtype) + if p.shape[1:] != (3,): + raise ValueError("Points must be Nx3 array-like") + return p + + def _validate_points(self, *points: ArrayLike) -> List[NDArray]: + """Ensure that arrays are equal-length sets of points.""" + ndim = None + out = [] + + for p_raw in points: + p = self._as_points(p_raw) + nd = p.shape[1:] + if ndim is None: + ndim = nd + elif ndim != nd: + raise ValueError("Point arrays are not the same shape") + out.append(p) + + return out + + def _sdf_intersections( + self, points: ArrayLike, vectors: ArrayLike + ) -> Tuple[NDArray, NDArray]: + p, v = self._validate_points(points, vectors) + return self._impl.sdf_intersections(p, v) + def intersections( self, src_points: ArrayLike, @@ -257,14 +285,7 @@ def intersections( float array of locations (Nx3), bool array of is_backface (N) """ - src = np.asarray(src_points, self.dtype) - tgt = np.asarray(tgt_points, self.dtype) - - if src.shape != tgt.shape: - raise ValueError("Source and target points arrays must be the same shape") - - if src.shape[1:] != (3,): - raise ValueError("Points must be Nx3 array-like") + src, tgt = self._validate_points(src_points, tgt_points) if self._interpret_threads(threads): idxs, points, bfs = self._impl.intersections_many_threaded( diff --git a/src/interface.rs b/src/interface.rs index f3a5319..7e2eeee 100644 --- a/src/interface.rs +++ b/src/interface.rs @@ -4,6 +4,7 @@ use std::iter::repeat_with; use numpy::ndarray::{Array, Zip}; use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray2}; use parry3d_f64::math::{Point, Vector}; +use parry3d_f64::query::{Ray, RayCast}; use parry3d_f64::shape::TriMesh; use pyo3::exceptions::PyRuntimeError; use pyo3::prelude::*; @@ -157,6 +158,67 @@ impl TriMeshWrapper { PyArray2::from_vec2(py, &[point_to_vec(&aabb.mins), point_to_vec(&aabb.maxs)]).unwrap() } + pub fn sdf_intersections<'py>( + &self, + py: Python<'py>, + points: PyReadonlyArray2, + vecs: PyReadonlyArray2, + ) -> (&'py PyArray1, &'py PyArray1) { + let diameter = self.mesh.local_bounding_sphere().radius() * 2.0; + + let (dists, dot_norms): (Vec, Vec) = points + .as_array() + .rows() + .into_iter() + .map(|p| Point::new(p[0], p[1], p[2])) + .zip( + vecs.as_array() + .rows() + .into_iter() + .map(|v| Vector::new(v[0], v[1], v[2]).normalize()), + ) + .map(|(p, v)| { + let ray = Ray::new(p, v); + if let Some(inter) = self.mesh.cast_local_ray_and_get_normal( + &ray, diameter, false, // unused + ) { + (inter.toi, v.dot(&inter.normal)) + } else { + (Precision::INFINITY, Precision::NAN) + } + }) + .unzip(); + ( + PyArray1::from_vec(py, dists), + PyArray1::from_vec(py, dot_norms), + ) + } + + // pub fn sdf_intersections_threaded<'py>( + // &self, + // py: Python<'py>, + // points: PyReadonlyArray2, + // vecs: PyReadonlyArray2, + // ) -> (&'py PyArray1, &'py PyArray1) { + // let diameter = self.mesh.local_bounding_sphere().radius() * 2.0; + + // Zip::from(points.as_array().rows()) + // .and(vecs.as_array().rows()) + // .par_map_collect(|point, vector| { + // let p = Point::new(point[0], point[1], point[2]); + // let v = Vector::new(vector[0], vector[1], vector[2]).normalize(); + + // let ray = Ray::new(p, v); + // if let Some(inter) = self.mesh.cast_local_ray_and_get_normal( + // &ray, diameter, false, // unused + // ) { + // (inter.toi, v.dot(&inter.normal)) + // } else { + // (Precision::INFINITY, Precision::NAN) + // } + // }) + // } + pub fn intersections_many<'py>( &self, py: Python<'py>, diff --git a/src/utils.rs b/src/utils.rs index b1e37a3..2081de6 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -1,6 +1,6 @@ use parry3d_f64::math::{Isometry, Point, Vector}; use parry3d_f64::query::{PointQuery, Ray, RayCast}; -use parry3d_f64::shape::TriMesh; +use parry3d_f64::shape::{FeatureId, TriMesh}; use rand::Rng; pub type Precision = f64; @@ -61,11 +61,20 @@ pub fn points_cross_mesh( src_point: &Point, tgt_point: &Point, ) -> Option<(Point, bool)> { + points_cross_mesh_info(mesh, src_point, tgt_point) + .map(|(inter, _, ft)| (inter, mesh.is_backface(ft))) +} + +pub fn points_cross_mesh_info( + mesh: &TriMesh, + src_point: &Point, + tgt_point: &Point, +) -> Option<(Point, Vector, FeatureId)> { let ray = Ray::new(*src_point, tgt_point - src_point); mesh.cast_local_ray_and_get_normal( &ray, 1.0, false, // unused ) - .map(|i| (ray.point_at(i.toi), mesh.is_backface(i.feature))) + .map(|i| (ray.point_at(i.toi), i.normal, i.feature)) } pub fn dist_from_mesh(mesh: &TriMesh, point: &Point, rays: Option<&[Vector]>) -> f64 { From 31e78de39f1dd470809539fc8b3d0a049228bbe7 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Thu, 21 Dec 2023 11:36:19 +0000 Subject: [PATCH 02/17] Parallelise SDF Also refactors to reduce copies into numpy arrays, and use vecs of vecs less. --- python/ncollpyde/_ncollpyde.pyi | 2 +- python/ncollpyde/main.py | 4 +- src/interface.rs | 207 ++++++++++++++++---------------- src/utils.rs | 27 +++++ tests/test_ncollpyde.py | 5 +- 5 files changed, 132 insertions(+), 113 deletions(-) diff --git a/python/ncollpyde/_ncollpyde.pyi b/python/ncollpyde/_ncollpyde.pyi index 1d60b2c..97f8763 100644 --- a/python/ncollpyde/_ncollpyde.pyi +++ b/python/ncollpyde/_ncollpyde.pyi @@ -31,5 +31,5 @@ class TriMeshWrapper: self, src_points: Points, tgt_points: Points ) -> Tuple[List[int], Points, List[bool]]: ... def sdf_intersections( - self, points: Points, vectors: Points + self, points: Points, vectors: Points, threaded: bool ) -> Tuple[npt.NDArray[np.float_], npt.NDArray[np.float_]]: ... diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index ed2da22..b3b986a 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -241,10 +241,10 @@ def _validate_points(self, *points: ArrayLike) -> List[NDArray]: return out def _sdf_intersections( - self, points: ArrayLike, vectors: ArrayLike + self, points: ArrayLike, vectors: ArrayLike, threads: Optional[bool] = None ) -> Tuple[NDArray, NDArray]: p, v = self._validate_points(points, vectors) - return self._impl.sdf_intersections(p, v) + return self._impl.sdf_intersections(p, v, self._interpret_threads(threads)) def intersections( self, diff --git a/src/interface.rs b/src/interface.rs index 7e2eeee..2af8f4f 100644 --- a/src/interface.rs +++ b/src/interface.rs @@ -1,10 +1,10 @@ use std::fmt::Debug; use std::iter::repeat_with; +use ndarray::{Array2, ArrayView1}; use numpy::ndarray::{Array, Zip}; use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray2}; use parry3d_f64::math::{Point, Vector}; -use parry3d_f64::query::{Ray, RayCast}; use parry3d_f64::shape::TriMesh; use pyo3::exceptions::PyRuntimeError; use pyo3::prelude::*; @@ -12,7 +12,10 @@ use rand::SeedableRng; use rand_pcg::Pcg64Mcg; use rayon::{prelude::*, ThreadPoolBuilder}; -use crate::utils::{dist_from_mesh, mesh_contains_point, points_cross_mesh, random_dir, Precision}; +use crate::utils::{ + aabb_diag, dist_from_mesh, mesh_contains_point, points_cross_mesh, random_dir, sdf_inner, + Precision, +}; fn vec_to_point(v: Vec) -> Point { Point::new(v[0], v[1], v[2]) @@ -22,14 +25,6 @@ fn point_to_vec(p: &Point) -> Vec { vec![p.x, p.y, p.z] } -fn vector_to_vec(v: &Vector) -> Vec { - vec![v[0], v[1], v[2]] -} - -// fn face_to_vec(f: &TriMeshFace) -> Vec { -// f.indices.iter().cloned().collect() -// } - #[pyclass] pub struct TriMeshWrapper { mesh: TriMesh, @@ -60,8 +55,7 @@ impl TriMeshWrapper { let mesh = TriMesh::new(points2, indices2); if n_rays > 0 { - let bsphere = mesh.local_bounding_sphere(); - let len = bsphere.radius() * 2.0; + let len = aabb_diag(&mesh); let mut rng = Pcg64Mcg::seed_from_u64(ray_seed); @@ -91,17 +85,17 @@ impl TriMeshWrapper { } else { None }; - if parallel { - Zip::from(points.as_array().rows()) - .par_map_collect(|v| { - dist_from_mesh(&self.mesh, &Point::new(v[0], v[1], v[2]), rays) - }) - .into_pyarray(py) + let p_arr = points.as_array(); + let zipped = Zip::from(p_arr.rows()); + let clos = + |v: ArrayView1| dist_from_mesh(&self.mesh, &Point::new(v[0], v[1], v[2]), rays); + + let collected = if parallel { + zipped.par_map_collect(clos) } else { - Zip::from(points.as_array().rows()) - .map_collect(|v| dist_from_mesh(&self.mesh, &Point::new(v[0], v[1], v[2]), rays)) - .into_pyarray(py) - } + zipped.map_collect(clos) + }; + collected.into_pyarray(py) } pub fn contains<'py>( @@ -110,52 +104,77 @@ impl TriMeshWrapper { points: PyReadonlyArray2, parallel: bool, ) -> &'py PyArray1 { - if parallel { - Zip::from(points.as_array().rows()) - .par_map_collect(|r| { - mesh_contains_point( - &self.mesh, - &Point::new(r[0], r[1], r[2]), - &self.ray_directions, - ) - }) - .into_pyarray(py) + let p_arr = points.as_array(); + let zipped = Zip::from(p_arr.rows()); + let clos = |r: ArrayView1| { + mesh_contains_point( + &self.mesh, + &Point::new(r[0], r[1], r[2]), + &self.ray_directions, + ) + }; + + let collected = if parallel { + zipped.par_map_collect(clos) } else { - Zip::from(points.as_array().rows()) - .map_collect(|r| { - mesh_contains_point( - &self.mesh, - &Point::new(r[0], r[1], r[2]), - &self.ray_directions, - ) - }) - .into_pyarray(py) - } + zipped.map_collect(clos) + }; + collected.into_pyarray(py) } pub fn points<'py>(&self, py: Python<'py>) -> &'py PyArray2 { - let vv: Vec> = self.mesh.vertices().iter().map(point_to_vec).collect(); - PyArray2::from_vec2(py, &vv).unwrap() + let vs = self.mesh.vertices(); + let v = vs + .iter() + .fold(Vec::with_capacity(vs.len() * 3), |mut out, p| { + out.push(p.x); + out.push(p.y); + out.push(p.z); + out + }); + Array2::from_shape_vec((vs.len(), 3), v) + .unwrap() + .into_pyarray(py) } pub fn faces<'py>(&self, py: Python<'py>) -> &'py PyArray2 { - let vv: Vec> = self - .mesh - .indices() - .iter() - .map(|arr| vec![arr[0], arr[1], arr[2]]) - .collect(); - PyArray2::from_vec2(py, &vv).unwrap() + let vs = self.mesh.indices(); + let v: Vec<_> = vs.iter().flatten().cloned().collect(); + Array2::from_shape_vec((vs.len(), 3), v) + .unwrap() + .into_pyarray(py) } pub fn rays<'py>(&self, py: Python<'py>) -> &'py PyArray2 { - let vv: Vec> = self.ray_directions.iter().map(vector_to_vec).collect(); - PyArray2::from_vec2(py, &vv).unwrap() + let vs = &self.ray_directions; + let v = vs + .iter() + .fold(Vec::with_capacity(vs.len() * 3), |mut out, p| { + out.push(p.x); + out.push(p.y); + out.push(p.z); + out + }); + Array2::from_shape_vec((vs.len(), 3), v) + .unwrap() + .into_pyarray(py) } pub fn aabb<'py>(&self, py: Python<'py>) -> &'py PyArray2 { let aabb = self.mesh.local_aabb(); - PyArray2::from_vec2(py, &[point_to_vec(&aabb.mins), point_to_vec(&aabb.maxs)]).unwrap() + Array2::from_shape_vec( + (2, 3), + vec![ + aabb.mins.x, + aabb.mins.y, + aabb.mins.z, + aabb.maxs.x, + aabb.maxs.y, + aabb.maxs.z, + ], + ) + .unwrap() + .into_pyarray(py) } pub fn sdf_intersections<'py>( @@ -163,61 +182,37 @@ impl TriMeshWrapper { py: Python<'py>, points: PyReadonlyArray2, vecs: PyReadonlyArray2, + threaded: bool, ) -> (&'py PyArray1, &'py PyArray1) { - let diameter = self.mesh.local_bounding_sphere().radius() * 2.0; + let diameter = aabb_diag(&self.mesh); - let (dists, dot_norms): (Vec, Vec) = points - .as_array() - .rows() - .into_iter() - .map(|p| Point::new(p[0], p[1], p[2])) - .zip( - vecs.as_array() - .rows() - .into_iter() - .map(|v| Vector::new(v[0], v[1], v[2]).normalize()), - ) - .map(|(p, v)| { - let ray = Ray::new(p, v); - if let Some(inter) = self.mesh.cast_local_ray_and_get_normal( - &ray, diameter, false, // unused - ) { - (inter.toi, v.dot(&inter.normal)) - } else { - (Precision::INFINITY, Precision::NAN) - } - }) - .unzip(); - ( - PyArray1::from_vec(py, dists), - PyArray1::from_vec(py, dot_norms), - ) - } + let n = points.shape()[0]; + + let mut dists = Array::from_elem((n,), 0.0); + let mut dot_norms = Array::from_elem((n,), 0.0); - // pub fn sdf_intersections_threaded<'py>( - // &self, - // py: Python<'py>, - // points: PyReadonlyArray2, - // vecs: PyReadonlyArray2, - // ) -> (&'py PyArray1, &'py PyArray1) { - // let diameter = self.mesh.local_bounding_sphere().radius() * 2.0; - - // Zip::from(points.as_array().rows()) - // .and(vecs.as_array().rows()) - // .par_map_collect(|point, vector| { - // let p = Point::new(point[0], point[1], point[2]); - // let v = Vector::new(vector[0], vector[1], vector[2]).normalize(); - - // let ray = Ray::new(p, v); - // if let Some(inter) = self.mesh.cast_local_ray_and_get_normal( - // &ray, diameter, false, // unused - // ) { - // (inter.toi, v.dot(&inter.normal)) - // } else { - // (Precision::INFINITY, Precision::NAN) - // } - // }) - // } + let p_arr = points.as_array(); + let v_arr = vecs.as_array(); + + let zipped = Zip::from(p_arr.rows()) + .and(v_arr.rows()) + .and(&mut dists) + .and(&mut dot_norms); + + let clos = |point, vector, dist: &mut f64, dot_norm: &mut f64| { + let (d, dn) = sdf_inner(point, vector, diameter, &self.mesh); + *dist = d; + *dot_norm = dn; + }; + + if threaded { + zipped.par_for_each(clos); + } else { + zipped.for_each(clos); + } + + (dists.into_pyarray(py), dot_norms.into_pyarray(py)) + } pub fn intersections_many<'py>( &self, diff --git a/src/utils.rs b/src/utils.rs index 2081de6..7445bef 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -1,3 +1,4 @@ +use ndarray::ArrayView1; use parry3d_f64::math::{Isometry, Point, Vector}; use parry3d_f64::query::{PointQuery, Ray, RayCast}; use parry3d_f64::shape::{FeatureId, TriMesh}; @@ -87,6 +88,32 @@ pub fn dist_from_mesh(mesh: &TriMesh, point: &Point, rays: Option<&[Vector< dist } +/// The diagonal length of the mesh's axis-aligned bounding box. +/// +/// Useful as an upper bound for ray length. +pub fn aabb_diag(mesh: &TriMesh) -> f64 { + mesh.local_aabb().extents().norm() +} + +pub fn sdf_inner( + point: ArrayView1, + vector: ArrayView1, + diameter: Precision, + mesh: &TriMesh, +) -> (Precision, Precision) { + let p = Point::new(point[0], point[1], point[2]); + let v = Vector::new(vector[0], vector[1], vector[2]).normalize(); + + let ray = Ray::new(p, v); + if let Some(inter) = mesh.cast_local_ray_and_get_normal( + &ray, diameter, false, // unused + ) { + (inter.toi, v.dot(&inter.normal)) + } else { + (Precision::INFINITY, Precision::NAN) + } +} + #[cfg(test)] mod tests { use std::fs::OpenOptions; diff --git a/tests/test_ncollpyde.py b/tests/test_ncollpyde.py index 743d143..45f7ea8 100644 --- a/tests/test_ncollpyde.py +++ b/tests/test_ncollpyde.py @@ -10,10 +10,7 @@ import numpy as np import pytest -try: - import trimesh -except ImportError: - trimesh = None +import trimesh from ncollpyde import PRECISION, Volume, configure_threadpool From 0e1fab41b05f81d9595ea9aeb8991c9a05ee391a Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Thu, 21 Dec 2023 11:56:50 +0000 Subject: [PATCH 03/17] Fix SDF function So that it is actually signed. Also some documentation. --- python/ncollpyde/main.py | 18 ++++++++++++++++++ src/interface.rs | 4 +--- src/utils.rs | 9 +++++++-- 3 files changed, 26 insertions(+), 5 deletions(-) diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index b3b986a..1a1c620 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -243,6 +243,24 @@ def _validate_points(self, *points: ArrayLike) -> List[NDArray]: def _sdf_intersections( self, points: ArrayLike, vectors: ArrayLike, threads: Optional[bool] = None ) -> Tuple[NDArray, NDArray]: + """Compute values required for signed distance field. + + :param points: Nx3 ndarray of floats + Points to calculate the distance from. + Should be within the axis-aligned bounding box of the mesh. + :param vectors: Nx3 ndarray of floats + Directions to fire rays from the given points. + :param threads: None, + Whether to parallelise the queries. If ``None`` (default), + refer to the instance's ``threads`` attribute. + :return: 2-tuple N-length np.ndarrays of floats. + The first is the distance, + which is negative if the collision is with a backface, + and infinity if there is no collision. + The second is the dot product of the vector + with the normal of the feature the ray hit, + NaN if there was no collision. + """ p, v = self._validate_points(points, vectors) return self._impl.sdf_intersections(p, v, self._interpret_threads(threads)) diff --git a/src/interface.rs b/src/interface.rs index 2af8f4f..86a71ba 100644 --- a/src/interface.rs +++ b/src/interface.rs @@ -127,9 +127,7 @@ impl TriMeshWrapper { let v = vs .iter() .fold(Vec::with_capacity(vs.len() * 3), |mut out, p| { - out.push(p.x); - out.push(p.y); - out.push(p.z); + out.extend(p.iter().cloned()); out }); Array2::from_shape_vec((vs.len(), 3), v) diff --git a/src/utils.rs b/src/utils.rs index 7445bef..6dc1470 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -79,7 +79,7 @@ pub fn points_cross_mesh_info( } pub fn dist_from_mesh(mesh: &TriMesh, point: &Point, rays: Option<&[Vector]>) -> f64 { - let mut dist = mesh.distance_to_point(&Isometry::identity(), point, true); + let mut dist = mesh.distance_to_local_point(point, true); if let Some(r) = rays { if mesh_contains_point(mesh, point, r) { dist = -dist; @@ -108,7 +108,12 @@ pub fn sdf_inner( if let Some(inter) = mesh.cast_local_ray_and_get_normal( &ray, diameter, false, // unused ) { - (inter.toi, v.dot(&inter.normal)) + let dot = v.dot(&inter.normal); + if mesh.is_backface(inter.feature) { + (inter.toi, dot) + } else { + (-inter.toi, dot) + } } else { (Precision::INFINITY, Precision::NAN) } From 042d7c95c9ab99ed6e2e7c8bcdcb407c09acd143 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Thu, 21 Dec 2023 12:06:40 +0000 Subject: [PATCH 04/17] bump dependencies --- .github/workflows/ci.yaml | 194 +++++++++++++++++++------------------- docs/requirements.txt | 2 +- requirements.txt | 16 ++-- 3 files changed, 106 insertions(+), 106 deletions(-) diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 867c726..4c80ddd 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -1,106 +1,106 @@ on: [push, pull_request] defaults: - run: - shell: bash + run: + shell: bash jobs: + lint-rust: + runs-on: ubuntu-22.04 + steps: + - uses: actions/checkout@v2 + - uses: actions-rs/toolchain@v1 + with: + toolchain: stable + components: rustfmt, clippy + - uses: Swatinem/rust-cache@v1 + - uses: actions/cache@v2 + with: + path: ~/.cache/pip + key: ${{ runner.os }}-pip-${{ hashFiles('**/requirements.txt') }} + restore-keys: | + ${{ runner.os }}-pip- + - run: make lint-rust - lint-rust: - runs-on: ubuntu-20.04 - steps: - - uses: actions/checkout@v2 - - uses: actions-rs/toolchain@v1 - with: - toolchain: stable - components: rustfmt, clippy - - uses: Swatinem/rust-cache@v1 - - uses: actions/cache@v2 - with: - path: ~/.cache/pip - key: ${{ runner.os }}-pip-${{ hashFiles('**/requirements.txt') }} - restore-keys: | - ${{ runner.os }}-pip- - - run: make lint-rust + lint-python: + runs-on: ubuntu-22.04 + steps: + - uses: actions/checkout@v2 + - uses: actions/setup-python@v2 + with: + python-version: "3.x" + - run: pip install $(grep -E '^(black|ruff|mypy|numpy)' requirements.txt) + - run: make lint-python - lint-python: - runs-on: ubuntu-20.04 - steps: - - uses: actions/checkout@v2 - - uses: actions/setup-python@v2 - with: - python-version: '3.x' - - run: pip install $(grep -E '^(black|ruff|mypy|numpy)' requirements.txt) - - run: make lint-python + test-rust: + runs-on: ubuntu-22.04 + steps: + - uses: actions/checkout@v2 + - uses: actions-rs/toolchain@v1 + with: + toolchain: stable + - uses: Swatinem/rust-cache@v1 + - run: cargo test - test-rust: - runs-on: ubuntu-20.04 - steps: - - uses: actions/checkout@v2 - - uses: actions-rs/toolchain@v1 - with: - toolchain: stable - - uses: Swatinem/rust-cache@v1 - - run: cargo test + test-python: + strategy: + fail-fast: false + matrix: + python-version: + - "3.8" + - "3.9" + - "3.10" + - "3.11" + - "3.12" + runs-on: ubuntu-22.04 + steps: + - uses: actions/checkout@v2 + - uses: actions-rs/toolchain@v1 + with: + toolchain: stable + - uses: Swatinem/rust-cache@v1 + - uses: actions/setup-python@v2 + with: + python-version: ${{ matrix.python-version }} + - uses: actions/cache@v2 + with: + path: ~/.cache/pip + key: ${{ runner.os }}-pip-${{ hashFiles('**/requirements.txt') }} + restore-keys: | + ${{ runner.os }}-pip- + - run: | + sudo apt-get install libspatialindex-dev + pip install -U pip wheel + pip install -r requirements.txt + name: Install dependencies + - run: | + mkdir -p $TGT_DIR + rm -f $TGT_DIR/*.whl + maturin build --release --interpreter python --out $TGT_DIR + pip install $TGT_DIR/*.whl + name: Install package + env: + TGT_DIR: "target/wheels/${{ matrix.python-version }}" + - run: pytest --verbose - test-python: - strategy: - fail-fast: false - matrix: - python-version: - - '3.8' - - '3.9' - - '3.10' - - '3.11' - runs-on: ubuntu-20.04 - steps: - - uses: actions/checkout@v2 - - uses: actions-rs/toolchain@v1 - with: - toolchain: stable - - uses: Swatinem/rust-cache@v1 - - uses: actions/setup-python@v2 - with: - python-version: ${{ matrix.python-version }} - - uses: actions/cache@v2 - with: - path: ~/.cache/pip - key: ${{ runner.os }}-pip-${{ hashFiles('**/requirements.txt') }} - restore-keys: | - ${{ runner.os }}-pip- - - run: | - sudo apt-get install libspatialindex-dev - pip install -U pip wheel - pip install -r requirements.txt - name: Install dependencies - - run: | - mkdir -p $TGT_DIR - rm -f $TGT_DIR/*.whl - maturin build --release --interpreter python --out $TGT_DIR - pip install $TGT_DIR/*.whl - name: Install package - env: - TGT_DIR: "target/wheels/${{ matrix.python-version }}" - - run: pytest --verbose - - deploy: - strategy: - matrix: - os: [macos-latest, windows-latest, ubuntu-20.04] - needs: [lint-rust, lint-python, test-rust, test-python] - runs-on: ${{ matrix.os }} - if: github.event_name == 'push' && contains(github.ref, 'refs/tags/v') - steps: - - uses: actions/checkout@v2 - - uses: actions-rs/toolchain@v1 - with: - toolchain: stable - - uses: actions/setup-python@v2 - with: - python-version: '3.9' - - uses: messense/maturin-action@v1 - with: - manylinux: auto - command: publish - args: -u __token__ -p ${{ secrets.MATURIN_PASSWORD }} --skip-existing --universal2 - name: Deploy wheels + deploy: + strategy: + matrix: + os: [macos-latest, windows-latest, ubuntu-20.04] + needs: [lint-rust, lint-python, test-rust, test-python] + runs-on: ${{ matrix.os }} + if: github.event_name == 'push' && contains(github.ref, 'refs/tags/v') + steps: + - uses: actions/checkout@v2 + - uses: actions-rs/toolchain@v1 + with: + toolchain: stable + - uses: actions/setup-python@v2 + with: + python-version: "3.9" + - uses: messense/maturin-action@v1 + with: + manylinux: auto + command: publish + args: -u __token__ -p ${{ secrets.MATURIN_PASSWORD }} --skip-existing --universal2 + name: Deploy wheels diff --git a/docs/requirements.txt b/docs/requirements.txt index f5988fb..5a95f62 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,2 +1,2 @@ -Sphinx==6.1.3 +Sphinx==7.2.6 sphinx-rtd-theme diff --git a/requirements.txt b/requirements.txt index cc1f9d5..39b7f0c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,14 +1,14 @@ # build -maturin==0.14.16 +maturin==1.4.0 # run -numpy==1.24.2 -trimesh[easy]==3.21.2 +numpy==1.26.2 +trimesh[easy]==4.0.7 # test meshio==5.3.4 -pytest==7.2.2 -pytest-runner==6.0.0 +pytest==7.4.3 +pytest-runner==6.0.1 # bench # pyoctree==0.2.10 @@ -16,11 +16,11 @@ pytest-benchmark==4.0.0 # develop pip -black==23.3.0 +black==23.12.0 watchdog==3.0.0 ruff -coverage==7.2.2 -mypy==1.1.1 +coverage==7.3.4 +mypy==1.7.1 # docs -r docs/requirements.txt From 0b3d2049d4d57603a799c4cf1c74d264777e649e Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Thu, 21 Dec 2023 12:11:13 +0000 Subject: [PATCH 05/17] bump python version --- .github/workflows/ci.yaml | 1 - pyproject.toml | 6 +++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 4c80ddd..e8c5a12 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -47,7 +47,6 @@ jobs: fail-fast: false matrix: python-version: - - "3.8" - "3.9" - "3.10" - "3.11" diff --git a/pyproject.toml b/pyproject.toml index 7798565..0351ef0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ name = "ncollpyde" # version = "0.19.0" description = "Point/line-mesh intersection queries in python" readme = "README.rst" -requires-python = ">=3.8" +requires-python = ">=3.9,<4" license = {file = "LICENSE"} authors = [ {name = "Chris L. Barnes", email = "chrislloydbarnes@gmail.com"}, @@ -15,10 +15,10 @@ classifiers = [ "Natural Language :: English", "Programming Language :: Rust", "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", ] dependencies = [ @@ -34,7 +34,7 @@ documentation = "https://ncollpyde.readthedocs.io/" repository = "https://github.com/clbarnes/ncollpyde/" [build-system] -requires = ["maturin==0.14", "numpy>=1.21"] +requires = ["maturin>=1.4", "numpy>=1.21"] build-backend = "maturin" [tool.maturin] From 5544551bdc390b9abc97c4fd0bace6af7e8eb2f8 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Thu, 21 Dec 2023 12:19:05 +0000 Subject: [PATCH 06/17] Use the correct normal --- src/utils.rs | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/utils.rs b/src/utils.rs index 6dc1470..cd18a28 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -1,7 +1,7 @@ use ndarray::ArrayView1; use parry3d_f64::math::{Isometry, Point, Vector}; use parry3d_f64::query::{PointQuery, Ray, RayCast}; -use parry3d_f64::shape::{FeatureId, TriMesh}; +use parry3d_f64::shape::{FeatureId, Shape, TriMesh}; use rand::Rng; pub type Precision = f64; @@ -108,7 +108,10 @@ pub fn sdf_inner( if let Some(inter) = mesh.cast_local_ray_and_get_normal( &ray, diameter, false, // unused ) { - let dot = v.dot(&inter.normal); + let normal = mesh + .feature_normal_at_point(inter.feature, &ray.point_at(inter.toi)) + .expect("Already checked collision"); + let dot = v.dot(&normal); if mesh.is_backface(inter.feature) { (inter.toi, dot) } else { From c0d7db1baaf708e4c9586338f1420a5c17f48efc Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Thu, 21 Dec 2023 12:36:24 +0000 Subject: [PATCH 07/17] add tests for SDF internals --- tests/test_ncollpyde.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/test_ncollpyde.py b/tests/test_ncollpyde.py index 45f7ea8..0a5980f 100644 --- a/tests/test_ncollpyde.py +++ b/tests/test_ncollpyde.py @@ -3,6 +3,7 @@ """Tests for `ncollpyde` package.""" from itertools import product +from math import pi, sqrt import sys import subprocess as sp import logging @@ -311,3 +312,17 @@ def test_configure_threadpool_twice(): with pytest.raises(RuntimeError): configure_threadpool(3, "prefix") configure_threadpool(3, "prefix") + + +@pytest.mark.parametrize( + ["point", "vec", "exp_dist", "exp_dot"], + [ + ([0.5, 0.5, 0.5], [1, 0, 0], 0.5, -1), + ([-0.5, 0.5, 0.5], [1, 0, 0], -0.5, -1), + ([0.75, 0.5, 0.5], [1, 1, 0], sqrt(2 * 0.25**2), -np.cos(pi / 4)), + ], +) +def test_sdf_inner(simple_volume: Volume, point, vec, exp_dist, exp_dot): + dists, dots = simple_volume._sdf_intersections([point], [vec]) + assert np.allclose(dists[0], exp_dist) + assert np.allclose(dots[0], exp_dot) From 87266724524b7503b8084ec2dfa677f170b3b0f7 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Thu, 21 Dec 2023 12:39:04 +0000 Subject: [PATCH 08/17] Improve conditional trimesh use --- python/ncollpyde/main.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index 1a1c620..915d7a0 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -7,11 +7,6 @@ import numpy as np from numpy.typing import ArrayLike, NDArray -try: - import trimesh -except ImportError: - trimesh = None - from ._ncollpyde import ( TriMeshWrapper, _index, @@ -136,7 +131,9 @@ def __init__( def _validate( self, vertices: np.ndarray, triangles: np.ndarray ) -> Tuple[NDArray[np.float64], NDArray[np.uint32]]: - if trimesh: + try: + import trimesh + tm = trimesh.Trimesh(vertices, triangles, validate=True) if not tm.is_volume: logger.info("Mesh not valid, attempting to fix") @@ -150,8 +147,7 @@ def _validate( ) return tm.vertices.astype(self.dtype), tm.faces.astype(np.uint32) - - else: + except ImportError: warnings.warn("trimesh not installed; full validation not possible") if vertices.shape[1:] != (3,): From 1b31ca12affaeb2237381641bf66ba449e30dfe2 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Thu, 21 Dec 2023 13:12:34 +0000 Subject: [PATCH 09/17] Update rtd config --- .readthedocs.yaml | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 01937d7..74ad84a 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -3,26 +3,26 @@ version: 2 sphinx: - builder: html + builder: html build: - os: "ubuntu-22.04" - tools: - python: "3.9" - apt_packages: - - curl - - build-essential - - gcc - - make - jobs: - pre_create_environment: - - curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y - pre_install: - - /bin/bash scripts/cargo_hack.sh - + os: "ubuntu-22.04" + tools: + python: "3.9" + rust: "1.70" + # apt_packages: + # - curl + # - build-essential + # - gcc + # - make + # jobs: + # pre_create_environment: + # - curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y + # pre_install: + # - /bin/bash scripts/cargo_hack.sh python: - install: - - requirements: docs/requirements.txt - - method: pip - path: . + install: + - requirements: docs/requirements.txt + - method: pip + path: . From 2e6ef5e34e4a0d510dde0892871ffbd72c49cdf0 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Wed, 10 Jan 2024 16:46:24 +0000 Subject: [PATCH 10/17] SDF internals, robust containment checks, benches - private method _sdf_intersections, which casts rays and returns the distance to a mesh and the dot product of the ray and the face normal - containment checks now may use signed distance, which is slower in real-world meshes but more robust - ray-based containment checks now may use consensus of a number of rays - different containment check strategies are benchmarked --- .gitignore | 2 + Cargo.lock | 465 ++++++++++++++++++++++++++++++-- Cargo.toml | 7 + README.rst | 2 +- benches/containment.rs | 206 ++++++++++++++ python/ncollpyde/_ncollpyde.pyi | 7 +- python/ncollpyde/main.py | 91 +++++-- src/interface.rs | 143 +++++++--- src/lib.rs | 2 +- src/utils.rs | 215 ++++++++++----- tests/conftest.py | 2 +- tests/test_bench.py | 19 ++ tests/test_ncollpyde.py | 30 ++- 13 files changed, 1040 insertions(+), 151 deletions(-) create mode 100644 benches/containment.rs diff --git a/.gitignore b/.gitignore index cf8bf8e..bc784d1 100644 --- a/.gitignore +++ b/.gitignore @@ -102,3 +102,5 @@ ENV/ .mypy_cache/ .benchmarks/ + +tmp/ diff --git a/Cargo.lock b/Cargo.lock index 425ee01..c17c328 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -11,6 +11,18 @@ dependencies = [ "memchr", ] +[[package]] +name = "anes" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4b46cbb362ab8752921c97e041f5e366ee6297bd428a31275b9fcf1e380f7299" + +[[package]] +name = "anstyle" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7079075b41f533b8c61d2a4d073c4676e1f8b249ff94a393b0595db304e0dd87" + [[package]] name = "approx" version = "0.5.1" @@ -49,6 +61,12 @@ version = "1.3.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" +[[package]] +name = "bitflags" +version = "2.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "327762f6e5a765692301e5bb513e0d9fef63be86bbc14528052b1cd3e6f03e07" + [[package]] name = "bstr" version = "0.2.17" @@ -70,6 +88,12 @@ dependencies = [ "serde", ] +[[package]] +name = "bumpalo" +version = "3.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f30e7476521f6f8af1a1c4c0b8cc94f0bee37d91763d0ca2665f299b6cd8aec" + [[package]] name = "bytemuck" version = "1.13.1" @@ -114,7 +138,7 @@ checksum = "e663fe9e3bdc9b639c25b5a6ca103c981545f680852aae2f05e7c50b88a5418b" dependencies = [ "bstr 0.2.17", "cargo_metadata", - "clap", + "clap 3.2.23", "clap-cargo", "crates-index", "difflib", @@ -149,6 +173,12 @@ dependencies = [ "serde_json", ] +[[package]] +name = "cast" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" + [[package]] name = "cc" version = "1.0.79" @@ -164,6 +194,33 @@ version = "1.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" +[[package]] +name = "ciborium" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "effd91f6c78e5a4ace8a5d3c0b6bfaec9e2baaef55f3efc00e45fb2e477ee926" +dependencies = [ + "ciborium-io", + "ciborium-ll", + "serde", +] + +[[package]] +name = "ciborium-io" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cdf919175532b369853f5d5e20b26b43112613fd6fe7aee757e35f7a44642656" + +[[package]] +name = "ciborium-ll" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "defaa24ecc093c77630e6c15e17c51f5e187bf35ee514f4e2d67baaa96dae22b" +dependencies = [ + "ciborium-io", + "half", +] + [[package]] name = "clap" version = "3.2.23" @@ -171,9 +228,9 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "71655c45cb9845d3270c9d6df84ebe72b4dad3c2ba3f7023ad47c144e4e473a5" dependencies = [ "atty", - "bitflags", + "bitflags 1.3.2", "clap_derive", - "clap_lex", + "clap_lex 0.2.4", "indexmap", "once_cell", "strsim", @@ -181,6 +238,15 @@ dependencies = [ "textwrap", ] +[[package]] +name = "clap" +version = "4.4.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bfaff671f6b22ca62406885ece523383b9b64022e341e53e009a62ebc47a45f2" +dependencies = [ + "clap_builder", +] + [[package]] name = "clap-cargo" version = "0.8.0" @@ -188,10 +254,20 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "551b6aa534ced210e29bc4ea2016bc11c74770f0a2b94b29dfc1a92ab6fc28d4" dependencies = [ "cargo_metadata", - "clap", + "clap 3.2.23", "doc-comment", ] +[[package]] +name = "clap_builder" +version = "4.4.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a216b506622bb1d316cd51328dce24e07bdff4a6128a47c7e7fad11878d5adbb" +dependencies = [ + "anstyle", + "clap_lex 0.6.0", +] + [[package]] name = "clap_derive" version = "3.2.18" @@ -214,6 +290,12 @@ dependencies = [ "os_str_bytes", ] +[[package]] +name = "clap_lex" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "702fc72eb24e5a1e48ce58027a675bc24edd52096d5397d4aea7c6dd9eca0bd1" + [[package]] name = "combine" version = "4.6.6" @@ -245,6 +327,42 @@ dependencies = [ "toml", ] +[[package]] +name = "criterion" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2b12d017a929603d80db1831cd3a24082f8137ce19c69e6447f54f5fc8d692f" +dependencies = [ + "anes", + "cast", + "ciborium", + "clap 4.4.11", + "criterion-plot", + "is-terminal", + "itertools", + "num-traits", + "once_cell", + "oorandom", + "plotters", + "rayon", + "regex", + "serde", + "serde_derive", + "serde_json", + "tinytemplate", + "walkdir", +] + +[[package]] +name = "criterion-plot" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6b50826342786a51a89e2da3a28f1c32b06e387201bc2d19791f622c673706b1" +dependencies = [ + "cast", + "itertools", +] + [[package]] name = "crossbeam-channel" version = "0.5.7" @@ -346,6 +464,16 @@ dependencies = [ "termcolor", ] +[[package]] +name = "errno" +version = "0.3.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a258e46cdc063eb8519c00b9fc845fc47bcfca4130e2f08e88665ceda8474245" +dependencies = [ + "libc", + "windows-sys 0.52.0", +] + [[package]] name = "float-cmp" version = "0.8.0" @@ -387,7 +515,7 @@ version = "0.14.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d0155506aab710a86160ddb504a480d2964d7ab5b9e62419be69e0032bc5931c" dependencies = [ - "bitflags", + "bitflags 1.3.2", "libc", "libgit2-sys", "log", @@ -409,6 +537,12 @@ dependencies = [ "regex", ] +[[package]] +name = "half" +version = "1.8.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eabb4a44450da02c90444cf74558da904edde8fb4e9035a9a6a4e15445af0bd7" + [[package]] name = "hashbrown" version = "0.12.3" @@ -439,6 +573,12 @@ dependencies = [ "libc", ] +[[package]] +name = "hermit-abi" +version = "0.3.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d77f7ec81a6d05a3abb01ab6eb7590f6083d08449fe5a1c8b1e620283546ccb7" + [[package]] name = "hex" version = "0.4.3" @@ -507,6 +647,17 @@ version = "1.0.9" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "bfa799dd5ed20a7e349f3b4639aa80d74549c81716d9ec4f994c9b5815598306" +[[package]] +name = "is-terminal" +version = "0.4.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cb0889898416213fab133e1d33a0e5858a48177452750691bde3666d0fdbaf8b" +dependencies = [ + "hermit-abi 0.3.3", + "rustix", + "windows-sys 0.48.0", +] + [[package]] name = "itertools" version = "0.10.5" @@ -531,6 +682,15 @@ dependencies = [ "libc", ] +[[package]] +name = "js-sys" +version = "0.3.66" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cee9c64da59eae3b50095c18d3e74f8b73c0b86d2792824ff01bbce68ba229ca" +dependencies = [ + "wasm-bindgen", +] + [[package]] name = "lazy_static" version = "1.4.0" @@ -539,9 +699,9 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" [[package]] name = "libc" -version = "0.2.140" +version = "0.2.151" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "99227334921fae1a979cf0bfdfcc6b3e5ce376ef57e16fb6fb3ea2ed6095f80c" +checksum = "302d7ab3130588088d277783b1e2d2e10c9e9e4a16dd9050e6ec93fb3e7048f4" [[package]] name = "libgit2-sys" @@ -589,6 +749,12 @@ dependencies = [ "vcpkg", ] +[[package]] +name = "linux-raw-sys" +version = "0.4.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c4cd1a83af159aa67994778be9070f0ae1bd732942279cabb14f86f986a21456" + [[package]] name = "lock_api" version = "0.4.9" @@ -601,12 +767,9 @@ dependencies = [ [[package]] name = "log" -version = "0.4.17" +version = "0.4.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "abb12e687cfb44aa40f41fc3978ef76448f9b6038cad6aef4259d3c095a2382e" -dependencies = [ - "cfg-if", -] +checksum = "b5e6163cb8c49088c2c36f57875e58ccd8c87c7427f7fbd50ea6710b2f3f2e8f" [[package]] name = "maplit" @@ -657,7 +820,10 @@ dependencies = [ name = "ncollpyde" version = "0.19.0" dependencies = [ + "approx", "cargo-release", + "criterion", + "log", "ndarray", "numpy", "parry3d-f64", @@ -764,6 +930,12 @@ version = "1.17.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b7e5500299e16ebb147ae15a00a942af264cf3688f47923b8fc2cd5858f23ad3" +[[package]] +name = "oorandom" +version = "11.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ab1bc2a289d34bd04a330323ac98a1b4bc82c9d9fcb1e66b63caa84da26b575" + [[package]] name = "openssl-probe" version = "0.1.5" @@ -815,7 +987,7 @@ dependencies = [ "libc", "redox_syscall", "smallvec", - "windows-sys", + "windows-sys 0.45.0", ] [[package]] @@ -826,7 +998,7 @@ checksum = "5dd9aeaccbb2bdc6b1ee9b9c7976a54523af68cb88112050b5f4b02935d5a346" dependencies = [ "approx", "arrayvec", - "bitflags", + "bitflags 1.3.2", "downcast-rs", "either", "indexmap", @@ -858,6 +1030,34 @@ version = "0.3.26" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "6ac9a59f73473f1b8d852421e59e64809f025994837ef743615c6d0c5b305160" +[[package]] +name = "plotters" +version = "0.3.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2c224ba00d7cadd4d5c660deaf2098e5e80e07846537c51f9cfa4be50c1fd45" +dependencies = [ + "num-traits", + "plotters-backend", + "plotters-svg", + "wasm-bindgen", + "web-sys", +] + +[[package]] +name = "plotters-backend" +version = "0.3.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e76628b4d3a7581389a35d5b6e2139607ad7c75b17aed325f210aa91f4a9609" + +[[package]] +name = "plotters-svg" +version = "0.3.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "38f6d39893cca0701371e3c27294f09797214b86f1fb951b89ade8ec04e2abab" +dependencies = [ + "plotters-backend", +] + [[package]] name = "ppv-lite86" version = "0.2.17" @@ -1045,7 +1245,7 @@ version = "0.2.16" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "fb5a58c1855b4b6819d59012155603f0b22ad30cad752600aadfcb695265519a" dependencies = [ - "bitflags", + "bitflags 1.3.2", ] [[package]] @@ -1094,6 +1294,19 @@ version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" +[[package]] +name = "rustix" +version = "0.38.28" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72e572a5e8ca657d7366229cdde4bd14c4eb5499a9573d4d366fe1b599daa316" +dependencies = [ + "bitflags 2.4.1", + "errno", + "libc", + "linux-raw-sys", + "windows-sys 0.52.0", +] + [[package]] name = "ryu" version = "1.0.13" @@ -1339,6 +1552,16 @@ dependencies = [ "time-core", ] +[[package]] +name = "tinytemplate" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "be4d6b5f19ff7664e8c98d03e2139cb510db9b0a60b55f8e8709b689d939b6bc" +dependencies = [ + "serde", + "serde_json", +] + [[package]] name = "tinyvec" version = "1.6.0" @@ -1447,6 +1670,70 @@ version = "0.11.0+wasi-snapshot-preview1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423" +[[package]] +name = "wasm-bindgen" +version = "0.2.89" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ed0d4f68a3015cc185aff4db9506a015f4b96f95303897bfa23f846db54064e" +dependencies = [ + "cfg-if", + "wasm-bindgen-macro", +] + +[[package]] +name = "wasm-bindgen-backend" +version = "0.2.89" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1b56f625e64f3a1084ded111c4d5f477df9f8c92df113852fa5a374dbda78826" +dependencies = [ + "bumpalo", + "log", + "once_cell", + "proc-macro2", + "quote", + "syn 2.0.11", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-macro" +version = "0.2.89" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0162dbf37223cd2afce98f3d0785506dcb8d266223983e4b5b525859e6e182b2" +dependencies = [ + "quote", + "wasm-bindgen-macro-support", +] + +[[package]] +name = "wasm-bindgen-macro-support" +version = "0.2.89" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0eb82fcb7930ae6219a7ecfd55b217f5f0893484b7a13022ebb2b2bf20b5283" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.11", + "wasm-bindgen-backend", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-shared" +version = "0.2.89" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ab9b36309365056cd639da3134bf87fa8f3d86008abf99e612384a6eecd459f" + +[[package]] +name = "web-sys" +version = "0.3.66" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "50c24a44ec86bb68fbecd1b3efed7e85ea5621b39b35ef2766b66cd984f8010f" +dependencies = [ + "js-sys", + "wasm-bindgen", +] + [[package]] name = "wide" version = "0.7.8" @@ -1494,7 +1781,25 @@ version = "0.45.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "75283be5efb2831d37ea142365f009c02ec203cd29a3ebecbc093d52315b66d0" dependencies = [ - "windows-targets", + "windows-targets 0.42.2", +] + +[[package]] +name = "windows-sys" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "677d2418bec65e3338edb076e806bc1ec15693c5d0104683f2efe857f61056a9" +dependencies = [ + "windows-targets 0.48.5", +] + +[[package]] +name = "windows-sys" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "282be5f36a8ce781fad8c8ae18fa3f9beff57ec1b52cb3de0789201425d9a33d" +dependencies = [ + "windows-targets 0.52.0", ] [[package]] @@ -1503,13 +1808,43 @@ version = "0.42.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8e5180c00cd44c9b1c88adb3693291f1cd93605ded80c250a75d472756b4d071" dependencies = [ - "windows_aarch64_gnullvm", - "windows_aarch64_msvc", - "windows_i686_gnu", - "windows_i686_msvc", - "windows_x86_64_gnu", - "windows_x86_64_gnullvm", - "windows_x86_64_msvc", + "windows_aarch64_gnullvm 0.42.2", + "windows_aarch64_msvc 0.42.2", + "windows_i686_gnu 0.42.2", + "windows_i686_msvc 0.42.2", + "windows_x86_64_gnu 0.42.2", + "windows_x86_64_gnullvm 0.42.2", + "windows_x86_64_msvc 0.42.2", +] + +[[package]] +name = "windows-targets" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a2fa6e2155d7247be68c096456083145c183cbbbc2764150dda45a87197940c" +dependencies = [ + "windows_aarch64_gnullvm 0.48.5", + "windows_aarch64_msvc 0.48.5", + "windows_i686_gnu 0.48.5", + "windows_i686_msvc 0.48.5", + "windows_x86_64_gnu 0.48.5", + "windows_x86_64_gnullvm 0.48.5", + "windows_x86_64_msvc 0.48.5", +] + +[[package]] +name = "windows-targets" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8a18201040b24831fbb9e4eb208f8892e1f50a37feb53cc7ff887feb8f50e7cd" +dependencies = [ + "windows_aarch64_gnullvm 0.52.0", + "windows_aarch64_msvc 0.52.0", + "windows_i686_gnu 0.52.0", + "windows_i686_msvc 0.52.0", + "windows_x86_64_gnu 0.52.0", + "windows_x86_64_gnullvm 0.52.0", + "windows_x86_64_msvc 0.52.0", ] [[package]] @@ -1518,38 +1853,122 @@ version = "0.42.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "597a5118570b68bc08d8d59125332c54f1ba9d9adeedeef5b99b02ba2b0698f8" +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b38e32f0abccf9987a4e3079dfb67dcd799fb61361e53e2882c3cbaf0d905d8" + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cb7764e35d4db8a7921e09562a0304bf2f93e0a51bfccee0bd0bb0b666b015ea" + [[package]] name = "windows_aarch64_msvc" version = "0.42.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e08e8864a60f06ef0d0ff4ba04124db8b0fb3be5776a5cd47641e942e58c4d43" +[[package]] +name = "windows_aarch64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc35310971f3b2dbbf3f0690a219f40e2d9afcf64f9ab7cc1be722937c26b4bc" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bbaa0368d4f1d2aaefc55b6fcfee13f41544ddf36801e793edbbfd7d7df075ef" + [[package]] name = "windows_i686_gnu" version = "0.42.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c61d927d8da41da96a81f029489353e68739737d3beca43145c8afec9a31a84f" +[[package]] +name = "windows_i686_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a75915e7def60c94dcef72200b9a8e58e5091744960da64ec734a6c6e9b3743e" + +[[package]] +name = "windows_i686_gnu" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a28637cb1fa3560a16915793afb20081aba2c92ee8af57b4d5f28e4b3e7df313" + [[package]] name = "windows_i686_msvc" version = "0.42.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "44d840b6ec649f480a41c8d80f9c65108b92d89345dd94027bfe06ac444d1060" +[[package]] +name = "windows_i686_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f55c233f70c4b27f66c523580f78f1004e8b5a8b659e05a4eb49d4166cca406" + +[[package]] +name = "windows_i686_msvc" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ffe5e8e31046ce6230cc7215707b816e339ff4d4d67c65dffa206fd0f7aa7b9a" + [[package]] name = "windows_x86_64_gnu" version = "0.42.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8de912b8b8feb55c064867cf047dda097f92d51efad5b491dfb98f6bbb70cb36" +[[package]] +name = "windows_x86_64_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53d40abd2583d23e4718fddf1ebec84dbff8381c07cae67ff7768bbf19c6718e" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3d6fa32db2bc4a2f5abeacf2b69f7992cd09dca97498da74a151a3132c26befd" + [[package]] name = "windows_x86_64_gnullvm" version = "0.42.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "26d41b46a36d453748aedef1486d5c7a85db22e56aff34643984ea85514e94a3" +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b7b52767868a23d5bab768e390dc5f5c55825b6d30b86c844ff2dc7414044cc" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1a657e1e9d3f514745a572a6846d3c7aa7dbe1658c056ed9c3344c4109a6949e" + [[package]] name = "windows_x86_64_msvc" version = "0.42.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9aec5da331524158c6d1a4ac0ab1541149c0b9505fde06423b02f5ef0106b9f0" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed94fce61571a4006852b7389a063ab983c02eb1bb37b47f8272ce92d06d9538" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dff9641d1cd4be8d1a070daf9e3773c5f67e78b4d9d42263020c057706765c04" diff --git a/Cargo.toml b/Cargo.toml index 52613da..40200c8 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -10,6 +10,7 @@ rayon = "1.7" rand = "0.8" rand_pcg = "0.3" numpy = "0.18" +log = "0.4.20" # ndarray is a dependency of numpy. # There is currently no way to directly specify sub-dependency features; @@ -23,6 +24,8 @@ features = ["rayon"] [dev-dependencies] stl_io = "0.6" cargo-release = "0.20" +criterion = "0.5.1" +approx = "0.5.1" [lib] name = "ncollpyde" @@ -31,3 +34,7 @@ crate-type = ["cdylib"] [package.metadata.release] disable-publish = true no-dev-version = true + +[[bench]] +name = "containment" +harness = false diff --git a/README.rst b/README.rst index fc8ff13..c5a4467 100644 --- a/README.rst +++ b/README.rst @@ -105,7 +105,7 @@ Known issues * Performance gains for multi-threaded queries are underwhelming, especially for ray intersections: see `this issue `_ * Very rare false positives for containment - * Due to a `bug in the underlying library `_ + * Due to a `bug in the underlying library `_ * Only happens when the point is outside the mesh and fires a ray which touches a single edge or vertex of the mesh. * Also affects ``is_backface`` result for ray intersection checks * manylinux-compatible wheels are built on CI but not necessarily in your local environment. Always allow CI to deploy the wheels. diff --git a/benches/containment.rs b/benches/containment.rs new file mode 100644 index 0000000..18c0d56 --- /dev/null +++ b/benches/containment.rs @@ -0,0 +1,206 @@ +use std::{fs::OpenOptions, iter::repeat_with, path::PathBuf}; + +use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion}; +use log::trace; +use parry3d_f64::{ + math::{Point, Vector}, + na::{Point3, Vector3}, + query::{PointQuery, Ray, RayCast}, + shape::{TriMesh, TriMeshFlags}, +}; +use rand::Rng; +use rand_pcg::Pcg64; +use stl_io::read_stl; + +type Precision = f64; + +const SEED: u128 = 1991; + +pub fn mesh_contains_point_ray( + mesh: &TriMesh, + point: &Point, + ray_direction: &Vector, +) -> bool { + let intersection_opt = mesh.cast_local_ray_and_get_normal( + &Ray::new(*point, *ray_direction), + 1.0, + false, // unused + ); + + if let Some(intersection) = intersection_opt { + mesh.is_backface(intersection.feature) + } else { + false + } +} + +pub fn mesh_contains_point_rays_simple( + mesh: &TriMesh, + point: &Point, + ray_directions: &[Vector], +) -> bool { + if ray_directions.is_empty() { + false + } else { + ray_directions + .iter() + .all(|v| mesh_contains_point_ray(mesh, point, v)) + } +} + +pub fn mesh_contains_point_rays( + mesh: &TriMesh, + point: &Point, + ray_directions: &[Vector], +) -> bool { + if !mesh.local_aabb().contains_local_point(point) { + trace!("Returning false, not in AABB"); + return false; + } + + // check whether point is on boundary + if mesh.contains_local_point(point) { + trace!("Returning true, point is on boundary"); + return true; + } + + mesh_contains_point_rays_simple(mesh, point, ray_directions) +} + +pub fn mesh_contains_point_oriented(mesh: &TriMesh, point: &Point) -> bool { + mesh.pseudo_normals().expect("Mesh orientation not checked"); + mesh.contains_local_point(point) +} + +pub fn mesh_contains_point_oriented_aabb(mesh: &TriMesh, point: &Point) -> bool { + mesh.pseudo_normals().expect("Mesh orientation not checked"); + mesh.local_aabb().contains_local_point(point) && mesh.contains_local_point(point) +} + +fn read_mesh(name: &'static str) -> TriMesh { + let mut stl_path = PathBuf::from(env!("CARGO_MANIFEST_DIR")) + .canonicalize() + .expect("couldn't resolve"); + stl_path.push("meshes"); + stl_path.push(name); + + let mut f = OpenOptions::new() + .read(true) + .open(stl_path) + .expect("Couldn't open file"); + let io_obj = read_stl(&mut f).expect("Couldn't parse STL"); + io_obj.validate().expect("Mesh is invalid"); + + TriMesh::new( + io_obj + .vertices + .iter() + .map(|v| Point::new(v[0] as Precision, v[1] as Precision, v[2] as Precision)) + .collect(), + io_obj + .faces + .iter() + .map(|t| { + [ + t.vertices[0] as u32, + t.vertices[1] as u32, + t.vertices[2] as u32, + ] + }) + .collect(), + ) +} + +fn orient(mesh: &TriMesh) -> TriMesh { + let mut m2 = mesh.clone(); + m2.set_flags(TriMeshFlags::ORIENTED).unwrap(); + m2 +} + +fn points_for_mesh(m: &TriMesh, n: usize) -> Vec> { + let aabb = m.local_aabb(); + let range = aabb.maxs - aabb.mins; + let frac = range * 0.1; + let mins = aabb.mins - frac; + let new_range = range + 2.0 * frac; + + let mut rng = Pcg64::new(0, SEED); + + repeat_with(|| { + Point3::from([ + rng.gen::() * new_range.x + mins.x, + rng.gen::() * new_range.y + mins.y, + rng.gen::() * new_range.z + mins.z, + ]) + }) + .take(n) + .collect() +} + +fn rays_for_mesh(mesh: &TriMesh, n: usize) -> Vec> { + let diameter = mesh.local_bounding_sphere().radius() * 2.0; + let mut rng = Pcg64::new(0, SEED); + + repeat_with(|| { + Vector3::from([rng.gen::(), rng.gen::(), rng.gen::()]).normalize() * diameter + }) + .take(n) + .collect() +} + +pub fn containment(c: &mut Criterion) { + let mut group = c.benchmark_group("containment"); + let max_rays = 100; + for name in &[ + "cube.stl", + "teapot.stl", + // "SEZ_right.stl", + ] { + let mesh_raw = read_mesh(name); + let points = points_for_mesh(&mesh_raw, 100); + let rays = rays_for_mesh(&mesh_raw, max_rays); + group.throughput(criterion::Throughput::Elements(points.len() as u64)); + + for n_rays in vec![0, 1, 2, 3, 4, 5] { + group.bench_with_input( + BenchmarkId::new(*name, format!("ray{n_rays}")), + &points, + |b, i| { + b.iter(|| { + i.iter().for_each(|p| { + mesh_contains_point_rays_simple( + &mesh_raw, + black_box(p), + &rays[..n_rays], + ); + }) + }) + }, + ); + } + + let mesh_oriented = orient(&mesh_raw); + group.bench_with_input(BenchmarkId::new(*name, "pseudonormal"), &points, |b, i| { + b.iter(|| { + i.iter().for_each(|p| { + mesh_contains_point_oriented(&mesh_oriented, black_box(p)); + }) + }) + }); + + group.bench_with_input( + BenchmarkId::new(*name, "pseudonormal_plus_aabb"), + &points, + |b, i| { + b.iter(|| { + i.iter().for_each(|p| { + mesh_contains_point_oriented_aabb(&mesh_oriented, black_box(p)); + }) + }) + }, + ); + } +} + +criterion_group!(benches, containment); +criterion_main!(benches); diff --git a/python/ncollpyde/_ncollpyde.pyi b/python/ncollpyde/_ncollpyde.pyi index 97f8763..2ac2874 100644 --- a/python/ncollpyde/_ncollpyde.pyi +++ b/python/ncollpyde/_ncollpyde.pyi @@ -16,7 +16,9 @@ class TriMeshWrapper: def __init__( self, points: Points, indices: Indices, n_rays: int, ray_seed: int ): ... - def contains(self, points: Points, parallel: bool) -> npt.NDArray[np.bool_]: ... + def contains( + self, points: Points, n_rays: Optional[int], consensus: int, parallel: bool + ) -> npt.NDArray[np.bool_]: ... def distance( self, points: Points, signed: bool, parallel: bool ) -> npt.NDArray[np.float64]: ... @@ -27,6 +29,9 @@ class TriMeshWrapper: def intersections_many( self, src_points: Points, tgt_points: Points ) -> Tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... + def intersections_many_threaded2( + self, src_points: Points, tgt_points: Points + ) -> Tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... def intersections_many_threaded( self, src_points: Points, tgt_points: Points ) -> Tuple[List[int], Points, List[bool]]: ... diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index 915d7a0..553435c 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -126,12 +126,16 @@ def __init__( ) ray_seed = random.randrange(0, 2**64) - self._impl = TriMeshWrapper(vert, tri, int(n_rays), ray_seed) + self.n_rays = int(n_rays) + inner_rays = 0 if self.n_rays < 0 else self.n_rays + + self._impl = TriMeshWrapper(vert, tri, inner_rays, ray_seed) def _validate( self, vertices: np.ndarray, triangles: np.ndarray ) -> Tuple[NDArray[np.float64], NDArray[np.uint32]]: try: + # todo: may not be necessary now parry can do some topology checks import trimesh tm = trimesh.Trimesh(vertices, triangles, validate=True) @@ -162,8 +166,12 @@ def _validate( return vertices, triangles def __contains__(self, item: ArrayLike) -> bool: - """Check whether a single point is in the volume.""" - return self.contains(np.asarray([item]), False)[0] + """Check whether a single point is in the volume. + + Uses the slower, more robust signed distance strategy. + For more control over the strategy, see the ``self.contains()`` method. + """ + return self.contains(np.asarray([item]), -1)[0] def _interpret_threads(self, threads: Optional[Union[int, bool]]) -> bool: return interpret_threads(threads, self.threads) @@ -172,6 +180,7 @@ def distance( self, coords: ArrayLike, signed: bool = True, + *, threads: Optional[bool] = None, ) -> np.ndarray: """Check the distance from the volume to multiple points (as a Px3 array-like). @@ -179,8 +188,6 @@ def distance( Distances are reported to the boundary of the volume. By default, if the point is inside the volume, the distance will be reported as negative. - ``signed=False`` is faster but cannot distinguish - between internal and external points. :param coords: :param signed: bool, default True. @@ -198,21 +205,59 @@ def distance( return self._impl.distance(coords, signed, self._interpret_threads(threads)) def contains( - self, coords: ArrayLike, threads: Optional[bool] = None + self, + coords: ArrayLike, + n_rays: Optional[int] = None, + consensus: Optional[int] = None, + *, + threads: Optional[bool] = None, ) -> NDArray[np.bool_]: """Check whether multiple points (as a Px3 array-like) are in the volume. :param coords: + :param n_rays: Optional[int] + If None, use the maximum rays defined on construction. + If < 0, use signed distance strategy + (more robust, but slower for many meshes). + Otherwise, use this many meshes, up to the maximum defined on construction. + :param consensus: Optional[int] + If using ray casting strategy, how many rays need to hit a backface + in order to call the point internal? + Rays from an external point may erroneously report hitting a backface + if they skim an edge (see https://github.com/clbarnes/ncollpyde/issues/3 ) + If None, will be set to ``n_rays // 2 + 1`` (i.e. >50%). + Ignored if using signed distance strategy. :param threads: None, Whether to parallelise the queries. If ``None`` (default), refer to the instance's ``threads`` attribute. :return: np.ndarray of bools, whether each point is inside the volume """ - coords = np.asarray(coords, self.dtype) - if coords.shape[1:] != (3,): - raise ValueError("Coords is not a Nx3 array-like") + coords = self._as_points(coords) + if n_rays is None: + n_rays = self.n_rays + elif n_rays < 0: + n_rays = None + elif n_rays > self.n_rays: + logger.warning( + "Requested %s rays, using the maximum of %s", n_rays, self.n_rays + ) + n_rays = self.n_rays - return self._impl.contains(coords, self._interpret_threads(threads)) + if n_rays is None: + consensus = 1 + else: + if consensus is None: + consensus = n_rays // 2 + 1 + elif consensus > n_rays: + raise ValueError( + "Requested consensus of %s rays but only casting %s", + consensus, + n_rays, + ) + + return self._impl.contains( + coords, n_rays, consensus, self._interpret_threads(threads) + ) def _as_points(self, points: ArrayLike) -> NDArray: p = np.asarray(points, self.dtype) @@ -237,7 +282,7 @@ def _validate_points(self, *points: ArrayLike) -> List[NDArray]: return out def _sdf_intersections( - self, points: ArrayLike, vectors: ArrayLike, threads: Optional[bool] = None + self, points: ArrayLike, vectors: ArrayLike, *, threads: Optional[bool] = None ) -> Tuple[NDArray, NDArray]: """Compute values required for signed distance field. @@ -246,6 +291,7 @@ def _sdf_intersections( Should be within the axis-aligned bounding box of the mesh. :param vectors: Nx3 ndarray of floats Directions to fire rays from the given points. + Need not be normalized. :param threads: None, Whether to parallelise the queries. If ``None`` (default), refer to the instance's ``threads`` attribute. @@ -253,8 +299,8 @@ def _sdf_intersections( The first is the distance, which is negative if the collision is with a backface, and infinity if there is no collision. - The second is the dot product of the vector - with the normal of the feature the ray hit, + The second is the absolute dot product of the normalized vector + with the unit normal of the feature the ray hit, NaN if there was no collision. """ p, v = self._validate_points(points, vectors) @@ -264,6 +310,7 @@ def intersections( self, src_points: ArrayLike, tgt_points: ArrayLike, + *, threads: Optional[bool] = None, ) -> Tuple[NDArray[np.uint64], NDArray[np.float64], NDArray[np.bool_]]: """Get intersections between line segments and volume. @@ -302,14 +349,7 @@ def intersections( src, tgt = self._validate_points(src_points, tgt_points) if self._interpret_threads(threads): - idxs, points, bfs = self._impl.intersections_many_threaded( - src.tolist(), tgt.tolist() - ) - return ( - np.array(idxs, INDEX), - np.array(points, self.dtype), - np.array(bfs, bool), - ) + return self._impl.intersections_many_threaded2(src, tgt) else: return self._impl.intersections_many(src, tgt) @@ -379,3 +419,12 @@ def rays(self) -> NDArray[np.float64]: and are the length of the diameter of the volume's bounding sphere. """ return self._impl.rays() + + +def points_around_vol(vol: Volume, n: int, pad: float = 0.2, seed=1991): + ext = vol.extents + ranges = ext[1] - ext[0] + to_pad = ranges * pad + + rng = np.random.default_rng(seed) + return rng.uniform(ext[0] - to_pad, ext[1] + to_pad, (n, 3)) diff --git a/src/interface.rs b/src/interface.rs index 86a71ba..d1b4c56 100644 --- a/src/interface.rs +++ b/src/interface.rs @@ -6,15 +6,15 @@ use numpy::ndarray::{Array, Zip}; use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray2}; use parry3d_f64::math::{Point, Vector}; use parry3d_f64::shape::TriMesh; -use pyo3::exceptions::PyRuntimeError; +use pyo3::exceptions::{PyRuntimeError, PyValueError}; use pyo3::prelude::*; use rand::SeedableRng; use rand_pcg::Pcg64Mcg; use rayon::{prelude::*, ThreadPoolBuilder}; use crate::utils::{ - aabb_diag, dist_from_mesh, mesh_contains_point, points_cross_mesh, random_dir, sdf_inner, - Precision, + aabb_diag, dist_from_mesh, mesh_contains_point, mesh_contains_point_oriented, + points_cross_mesh, random_dir, sdf_inner, Precision, FLAGS, }; fn vec_to_point(v: Vec) -> Point { @@ -39,7 +39,7 @@ impl TriMeshWrapper { indices: PyReadonlyArray2, n_rays: usize, ray_seed: u64, - ) -> Self { + ) -> PyResult { let points2 = points .as_array() .rows() @@ -52,24 +52,27 @@ impl TriMeshWrapper { .into_iter() .map(|v| [v[0], v[1], v[2]]) .collect(); - let mesh = TriMesh::new(points2, indices2); + let mut mesh = TriMesh::new(points2, indices2); + + mesh.set_flags(FLAGS) + .map_err(|e| PyValueError::new_err(format!("Invalid mesh topology: {e}")))?; if n_rays > 0 { let len = aabb_diag(&mesh); let mut rng = Pcg64Mcg::seed_from_u64(ray_seed); - Self { + Ok(Self { mesh, ray_directions: repeat_with(|| random_dir(&mut rng, len)) .take(n_rays) .collect(), - } + }) } else { - Self { + Ok(Self { mesh, ray_directions: Vec::default(), - } + }) } } @@ -80,15 +83,10 @@ impl TriMeshWrapper { signed: bool, parallel: bool, ) -> &'py PyArray1 { - let rays = if signed { - Some(&self.ray_directions[..]) - } else { - None - }; let p_arr = points.as_array(); let zipped = Zip::from(p_arr.rows()); let clos = - |v: ArrayView1| dist_from_mesh(&self.mesh, &Point::new(v[0], v[1], v[2]), rays); + |v: ArrayView1| dist_from_mesh(&self.mesh, &Point::new(v[0], v[1], v[2]), signed); let collected = if parallel { zipped.par_map_collect(clos) @@ -102,23 +100,38 @@ impl TriMeshWrapper { &self, py: Python<'py>, points: PyReadonlyArray2, + n_rays: Option, + consensus: usize, parallel: bool, ) -> &'py PyArray1 { let p_arr = points.as_array(); let zipped = Zip::from(p_arr.rows()); - let clos = |r: ArrayView1| { - mesh_contains_point( - &self.mesh, - &Point::new(r[0], r[1], r[2]), - &self.ray_directions, - ) - }; - let collected = if parallel { - zipped.par_map_collect(clos) + let collected = if let Some(n) = n_rays { + // use ray casting + let rays = &self.ray_directions[..n]; + let clos = |r: ArrayView1| { + mesh_contains_point(&self.mesh, &Point::new(r[0], r[1], r[2]), rays, consensus) + }; + + if parallel { + zipped.par_map_collect(clos) + } else { + zipped.map_collect(clos) + } } else { - zipped.map_collect(clos) + // use pseudonormals + let clos = |r: ArrayView1| { + mesh_contains_point_oriented(&self.mesh, &Point::new(r[0], r[1], r[2])) + }; + + if parallel { + zipped.par_map_collect(clos) + } else { + zipped.map_collect(clos) + } }; + collected.into_pyarray(py) } @@ -225,37 +238,34 @@ impl TriMeshWrapper { let mut idxs = Vec::default(); let mut intersections = Vec::default(); let mut is_backface = Vec::default(); - let mut count = 0; - for (idx, point, is_bf) in src_points + src_points .as_array() .rows() .into_iter() .zip(tgt_points.as_array().rows().into_iter()) .zip(0_u64..) - .filter_map(|((src, tgt), i)| { - points_cross_mesh( + .for_each(|((src, tgt), i)| { + if let Some((pt, is_bf)) = points_cross_mesh( &self.mesh, &Point::new(src[0], src[1], src[2]), &Point::new(tgt[0], tgt[1], tgt[2]), - ) - .map(|o| (i, o.0, o.1)) - }) - { - idxs.push(idx); - intersections.extend(point.iter().cloned()); - is_backface.push(is_bf); - count += 1; - } + ) { + idxs.push(i); + intersections.extend(pt.iter()); + is_backface.push(is_bf); + } + }); ( PyArray1::from_vec(py, idxs), - Array::from_shape_vec((count, 3), intersections) + Array::from_shape_vec((is_backface.len(), 3), intersections) .unwrap() .into_pyarray(py), PyArray1::from_vec(py, is_backface), ) } + #[deprecated] pub fn intersections_many_threaded( &self, src_points: Vec>, @@ -273,6 +283,61 @@ impl TriMeshWrapper { (idxs, intersections, is_backface) } + + pub fn intersections_many_threaded2<'py>( + &self, + py: Python<'py>, + src_points: PyReadonlyArray2, + tgt_points: PyReadonlyArray2, + ) -> ( + &'py PyArray1, + &'py PyArray2, + &'py PyArray1, + ) { + let src_arr = src_points.as_array(); + let tgt_arr = tgt_points.as_array(); + let zipped = Zip::indexed(src_arr.rows()).and(tgt_arr.rows()); + + let (out_idxs, out_pts, out_is_bf) = zipped + .into_par_iter() + // calculate the output rows, discarding non-intersections + .filter_map(|(idx, src, tgt)| { + points_cross_mesh( + &self.mesh, + &Point::new(src[0], src[1], src[2]), + &Point::new(tgt[0], tgt[1], tgt[2]), + ) + .map(|o| (idx as u64, o.0, o.1)) + }) + // convert chunks of results into flat vecs + .fold( + || (vec![], vec![], vec![]), + |(mut idx, mut pts, mut is_backface), (i, pt, is_bf)| { + idx.push(i); + pts.extend(pt.iter()); + is_backface.push(is_bf); + (idx, pts, is_backface) + }, + ) + // concatenate the chunked flat vecs + .reduce( + || (vec![], vec![], vec![]), + |(mut idxs, mut pts, mut is_backfaces), (mut idx, mut pt, mut is_bf)| { + idxs.append(&mut idx); + pts.append(&mut pt); + is_backfaces.append(&mut is_bf); + (idxs, pts, is_backfaces) + }, + ); + + ( + PyArray1::from_vec(py, out_idxs), + Array::from_shape_vec((out_is_bf.len(), 3), out_pts) + .unwrap() + .into_pyarray(py), + PyArray1::from_vec(py, out_is_bf), + ) + } } #[pymodule] diff --git a/src/lib.rs b/src/lib.rs index 14cbf3e..0567ded 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,4 @@ -mod utils; +pub mod utils; mod interface; diff --git a/src/utils.rs b/src/utils.rs index cd18a28..8b74549 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -1,11 +1,20 @@ use ndarray::ArrayView1; -use parry3d_f64::math::{Isometry, Point, Vector}; +use parry3d_f64::math::{Point, Vector}; +use parry3d_f64::na::{distance, Unit}; use parry3d_f64::query::{PointQuery, Ray, RayCast}; -use parry3d_f64::shape::{FeatureId, Shape, TriMesh}; +use parry3d_f64::shape::{FeatureId, Shape, TriMesh, TriMeshFlags}; use rand::Rng; pub type Precision = f64; +pub const FLAGS: TriMeshFlags = TriMeshFlags::empty() + .union(TriMeshFlags::ORIENTED) + .union(TriMeshFlags::DELETE_BAD_TOPOLOGY_TRIANGLES) + .union(TriMeshFlags::HALF_EDGE_TOPOLOGY) + .union(TriMeshFlags::MERGE_DUPLICATE_VERTICES) + .union(TriMeshFlags::DELETE_DEGENERATE_TRIANGLES) + .union(TriMeshFlags::DELETE_DUPLICATE_TRIANGLES); + pub fn random_dir(rng: &mut R, length: Precision) -> Vector { let unscaled: Vector = [ rng.gen::() - 0.5, @@ -38,23 +47,47 @@ pub fn mesh_contains_point( mesh: &TriMesh, point: &Point, ray_directions: &[Vector], + consensus: usize, ) -> bool { if !mesh.local_aabb().contains_local_point(point) { return false; } - // check whether point is on boundary - if mesh.contains_point(&Isometry::identity(), point) { - return true; + // this previously checked if point was on boundary, + // now it does a full (slow) containment check + // if mesh.contains_local_point(point) { + // return true; + // } + + // ray_directions + // .iter() + // .filter(|v| mesh_contains_point_ray(mesh, point, v)) + // .count() + // >= consensus + + let mut inside_remaining = consensus; + let mut remaining = ray_directions.len(); + for contains in ray_directions + .iter() + .map(|v| mesh_contains_point_ray(mesh, point, v)) + { + remaining -= 1; + if contains { + if inside_remaining == 1 { + // early return if we've met consensus + return true; + } + inside_remaining -= 1; + } else if remaining < inside_remaining { + // early return if there aren't enough rays left to reach consensus + return false; + } } + return false; +} - if ray_directions.is_empty() { - false - } else { - ray_directions - .iter() - .all(|v| mesh_contains_point_ray(mesh, point, v)) - } +pub fn mesh_contains_point_oriented(mesh: &TriMesh, point: &Point) -> bool { + mesh.local_aabb().contains_local_point(point) && mesh.contains_local_point(point) } pub fn points_cross_mesh( @@ -78,14 +111,14 @@ pub fn points_cross_mesh_info( .map(|i| (ray.point_at(i.toi), i.normal, i.feature)) } -pub fn dist_from_mesh(mesh: &TriMesh, point: &Point, rays: Option<&[Vector]>) -> f64 { - let mut dist = mesh.distance_to_local_point(point, true); - if let Some(r) = rays { - if mesh_contains_point(mesh, point, r) { - dist = -dist; - } +pub fn dist_from_mesh(mesh: &TriMesh, point: &Point, signed: bool) -> f64 { + let pp = mesh.project_local_point(point, true); + let dist = distance(&pp.point, point); + if signed && pp.is_inside { + -dist + } else { + dist } - dist } /// The diagonal length of the mesh's axis-aligned bounding box. @@ -95,23 +128,17 @@ pub fn aabb_diag(mesh: &TriMesh) -> f64 { mesh.local_aabb().extents().norm() } -pub fn sdf_inner( - point: ArrayView1, - vector: ArrayView1, - diameter: Precision, +pub fn ray_toi_dot( + p: Point, + v: Unit>, + length: Precision, mesh: &TriMesh, ) -> (Precision, Precision) { - let p = Point::new(point[0], point[1], point[2]); - let v = Vector::new(vector[0], vector[1], vector[2]).normalize(); - - let ray = Ray::new(p, v); + let ray = Ray::new(p, *v); if let Some(inter) = mesh.cast_local_ray_and_get_normal( - &ray, diameter, false, // unused + &ray, length, true, // unused ) { - let normal = mesh - .feature_normal_at_point(inter.feature, &ray.point_at(inter.toi)) - .expect("Already checked collision"); - let dot = v.dot(&normal); + let dot = v.dot(&inter.normal).abs(); if mesh.is_backface(inter.feature) { (inter.toi, dot) } else { @@ -122,8 +149,26 @@ pub fn sdf_inner( } } +/// Find the distance to a mesh boundary from a point along a particular direction. +/// +/// Returns a tuple of the distance and the absolute dot product between the vector and the face normal. +/// The distance is positive if the intersection was with a backface, +/// negative if the intersection was with an external face. +pub fn sdf_inner( + point: ArrayView1, + vector: ArrayView1, + diameter: Precision, + mesh: &TriMesh, +) -> (Precision, Precision) { + let p = Point::new(point[0], point[1], point[2]); + let v = Unit::new_normalize(Vector::new(vector[0], vector[1], vector[2])); + ray_toi_dot(p, v, diameter, mesh) +} + #[cfg(test)] mod tests { + use approx::assert_relative_eq; + use std::f64::consts::TAU; use std::fs::OpenOptions; use std::path::PathBuf; use stl_io::read_stl; @@ -146,7 +191,7 @@ mod tests { let io_obj = read_stl(&mut f).expect("Couldn't parse STL"); io_obj.validate().expect("Mesh is invalid"); - TriMesh::new( + let mut tm = TriMesh::new( io_obj .vertices .iter() @@ -163,7 +208,9 @@ mod tests { ] }) .collect(), - ) + ); + tm.set_flags(FLAGS).unwrap(); + tm } fn cube() -> TriMesh { @@ -172,17 +219,17 @@ mod tests { #[test] fn corner_contains() { - assert!(cube().contains_point(&Isometry::identity(), &Point::new(0.0, 0.0, 0.0))) + assert!(cube().contains_local_point(&Point::new(0.0, 0.0, 0.0))) } #[test] fn edge_contains() { - assert!(cube().contains_point(&Isometry::identity(), &Point::new(0.5, 0.0, 0.0))) + assert!(cube().contains_local_point(&Point::new(0.5, 0.0, 0.0))) } #[test] fn face_contains() { - assert!(cube().contains_point(&Isometry::identity(), &Point::new(0.5, 0.5, 0.0))) + assert!(cube().contains_local_point(&Point::new(0.5, 0.5, 0.0))) } fn assert_ray(p: [Precision; 3], v: [Precision; 3], is_inside: bool) { @@ -309,41 +356,85 @@ mod tests { assert!(get_cross([1.5, 1.5, 1.5], [2.0, 2.0, 2.0]).is_none()); } - fn axis_rays() -> Vec> { - vec![ - Vector::new(1.0, 0.0, 0.0), - Vector::new(0.0, 1.0, 0.0), - Vector::new(0.0, 0.0, 1.0), - ] - } - - fn assert_dist( - mesh: &TriMesh, - point: &Point, - rays: Option<&[Vector]>, - expected: Precision, - ) { - assert_eq!(dist_from_mesh(mesh, point, rays), expected) + fn assert_dist(mesh: &TriMesh, point: &Point, signed: bool, expected: Precision) { + assert_eq!(dist_from_mesh(mesh, point, signed), expected) } #[test] fn distance_signed() { - let rays = axis_rays(); let cube = cube(); - assert_dist(&cube, &Point::new(1.0, 1.0, 1.0), Some(&rays), 0.0); - assert_dist(&cube, &Point::new(0.5, 0.5, 0.5), Some(&rays), -0.5); - assert_dist(&cube, &Point::new(2.0, 1.0, 1.0), Some(&rays), 1.0); + assert_dist(&cube, &Point::new(1.0, 1.0, 1.0), true, 0.0); + assert_dist(&cube, &Point::new(0.5, 0.5, 0.5), true, -0.5); + assert_dist(&cube, &Point::new(2.0, 1.0, 1.0), true, 1.0); let three: Precision = 3.0; - assert_dist(&cube, &Point::new(2.0, 2.0, 2.0), Some(&rays), three.sqrt()); + assert_dist(&cube, &Point::new(2.0, 2.0, 2.0), true, three.sqrt()); } #[test] fn distance_unsigned() { let cube = cube(); - assert_dist(&cube, &Point::new(1.0, 1.0, 1.0), None, 0.0); - assert_dist(&cube, &Point::new(0.5, 0.5, 0.5), None, 0.5); - assert_dist(&cube, &Point::new(2.0, 1.0, 1.0), None, 1.0); - let three: Precision = 3.0; - assert_dist(&cube, &Point::new(2.0, 2.0, 2.0), None, three.sqrt()); + assert_dist(&cube, &Point::new(0.5, 0.5, 0.5), false, 0.5); + } + + #[test] + fn containment_with_psnorms() { + let cube = cube(); + assert!( + !cube.contains_local_point(&[1.5, 0.5, 0.5].into()), + "containment check failed for outside" + ); + assert!( + cube.contains_local_point(&[0.5, 0.5, 0.5].into()), + "containment check failed for center" + ); + assert!( + cube.contains_local_point(&[0.0, 0.0, 0.0].into()), + "containment check failed for vertex" + ); + assert!( + cube.contains_local_point(&[0.5, 0.0, 0.0].into()), + "containment check failed for edge" + ); + assert!( + cube.contains_local_point(&[0.5, 0.5, 0.0].into()), + "containment check failed for face" + ); + } + + fn get_sdf_results( + point: &[Precision; 3], + direction: &[Precision; 3], + mesh: &TriMesh, + ) -> (Precision, Precision) { + let length = aabb_diag(mesh); + let v = Unit::new_normalize(Vector::new(direction[0], direction[1], direction[2])); + ray_toi_dot(point.clone().into(), v, length, mesh) + } + + #[test] + fn sdf_miss() { + let vol = cube(); + let (dist, dot) = get_sdf_results(&[2.0, 2.0, 2.0], &[1.0, 1.0, 1.0], &vol); + assert_eq!(dist, Precision::INFINITY); + assert!(dot.is_nan()); + } + + #[test] + fn sdf_direct() { + let vol = cube(); + let res = get_sdf_results(&[-0.5, 0.5, 0.5], &[1.0, 0.0, 0.0], &vol); + assert_eq!(res, (-0.5, 1.0)); + } + + #[test] + fn sdf_diag() { + let vol = cube(); + let offset = -0.1; + let angle = TAU / 8.; // 45deg + let exp_dist = offset / angle.cos(); + let exp_dot = (TAU / 8.).cos(); + let (dist, dot) = get_sdf_results(&[offset, 0.5, 0.5], &[1.0, 1.0, 0.0], &vol); + assert_relative_eq!(dist, exp_dist); + assert_relative_eq!(dot, exp_dot); } } diff --git a/tests/conftest.py b/tests/conftest.py index 5d94a13..9cd3fb2 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -17,7 +17,7 @@ def mesh(): @pytest.fixture def volume(mesh): - return Volume.from_meshio(meshio, validate=True) + return Volume.from_meshio(mesh, validate=True) @pytest.fixture diff --git a/tests/test_bench.py b/tests/test_bench.py index 387f2b4..4dc5d1f 100644 --- a/tests/test_bench.py +++ b/tests/test_bench.py @@ -5,6 +5,7 @@ import pytest from ncollpyde import Volume +from ncollpyde.main import points_around_vol SAMPLES_PER_DIM = 10 PADDING = 0.2 @@ -15,6 +16,7 @@ INTERSECTION_SERIAL = "intersection serial" CONTAINS_PARALLEL = "containment parallel" INTERSECTION_PARALLEL = "intersection parallel" +INTERSECTION_PARALLEL_IMPL = "intersection parallel implementation" DISTANCE = "distance" @@ -260,6 +262,23 @@ def test_ncollpyde_intersection(mesh, benchmark, threads): benchmark(vol.intersections, starts, stops) +@pytest.mark.benchmark(group=INTERSECTION_PARALLEL_IMPL) +@pytest.mark.parametrize( + ("method_name",), + [("intersections_many_threaded",), ("intersections_many_threaded2",)], +) +def test_ncollpyde_intersection_impls(mesh, benchmark, method_name): + n_edges = 1_000 + + vol = Volume.from_meshio(mesh, threads=True) + + starts = points_around_vol(vol, n_edges, seed=1991) + stops = points_around_vol(vol, n_edges, seed=1992) + + fn = getattr(vol._impl, method_name) + benchmark(fn, starts, stops) + + @pytest.mark.benchmark(group=DISTANCE) @pytest.mark.parametrize("signed", [True, False]) def test_ncollpyde_distance(mesh, benchmark, signed): diff --git a/tests/test_ncollpyde.py b/tests/test_ncollpyde.py index 0a5980f..6205289 100644 --- a/tests/test_ncollpyde.py +++ b/tests/test_ncollpyde.py @@ -7,6 +7,7 @@ import sys import subprocess as sp import logging +from ncollpyde.main import points_around_vol import numpy as np import pytest @@ -49,6 +50,13 @@ def test_many(mesh, threads): assert np.array_equal(vol.contains(points, threads=threads), expected) +def test_contains_results(volume: Volume): + pts = points_around_vol(volume, 1000, 0.1) + ray = volume.contains(pts, n_rays=3, consensus=3, threads=False) + psnorms = volume.contains(pts, n_rays=-1, threads=False) + assert np.allclose(ray, psnorms) + + def test_0_rays(mesh): vol = Volume.from_meshio(mesh, n_rays=0) points = [p for p, _ in points_expected] @@ -233,7 +241,7 @@ def test_intersections_threads(simple_volume, threads): [1.5, 0.5, 1.5], ] - idxs, _, _ = simple_volume.intersections(sources, targets, threads) + idxs, _, _ = simple_volume.intersections(sources, targets, threads=threads) assert len(idxs) == 1 assert idxs[0] == 1 @@ -281,7 +289,7 @@ def test_near_miss(simple_volume: Volume, steps, angle): ([0.5, 0.5, 0.5], -0.5), ], ) -def test_distance_unsigned(simple_volume, point, expected, signed): +def test_distance(simple_volume, point, expected, signed): if not signed: expected = np.abs(expected) assert np.allclose( @@ -326,3 +334,21 @@ def test_sdf_inner(simple_volume: Volume, point, vec, exp_dist, exp_dot): dists, dots = simple_volume._sdf_intersections([point], [vec]) assert np.allclose(dists[0], exp_dist) assert np.allclose(dots[0], exp_dot) + + +def assert_intersection_results(test, ref): + test_dict = {idx: (tuple(pt), is_bf) for idx, pt, is_bf in zip(*test)} + ref_dict = {idx: (tuple(pt), is_bf) for idx, pt, is_bf in zip(*ref)} + assert test_dict == ref_dict + + +def test_intersections_impls(volume: Volume): + n_edges = 1000 + starts = points_around_vol(volume, n_edges, seed=1991) + stops = points_around_vol(volume, n_edges, seed=1992) + + serial = volume.intersections(starts, stops, threads=False) + par = volume.intersections(starts, stops, threads=True) + assert_intersection_results(serial, par) + par2 = volume._impl.intersections_many_threaded2(starts, stops) + assert_intersection_results(serial, par2) From 6138656561f38584a9c8f386f49d78ec28201c10 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Wed, 10 Jan 2024 17:02:41 +0000 Subject: [PATCH 11/17] Fix python tests for SDF internals --- src/utils.rs | 2 +- tests/test_ncollpyde.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/utils.rs b/src/utils.rs index 8b74549..5594583 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -2,7 +2,7 @@ use ndarray::ArrayView1; use parry3d_f64::math::{Point, Vector}; use parry3d_f64::na::{distance, Unit}; use parry3d_f64::query::{PointQuery, Ray, RayCast}; -use parry3d_f64::shape::{FeatureId, Shape, TriMesh, TriMeshFlags}; +use parry3d_f64::shape::{FeatureId, TriMesh, TriMeshFlags}; use rand::Rng; pub type Precision = f64; diff --git a/tests/test_ncollpyde.py b/tests/test_ncollpyde.py index 6205289..2130375 100644 --- a/tests/test_ncollpyde.py +++ b/tests/test_ncollpyde.py @@ -325,9 +325,9 @@ def test_configure_threadpool_twice(): @pytest.mark.parametrize( ["point", "vec", "exp_dist", "exp_dot"], [ - ([0.5, 0.5, 0.5], [1, 0, 0], 0.5, -1), - ([-0.5, 0.5, 0.5], [1, 0, 0], -0.5, -1), - ([0.75, 0.5, 0.5], [1, 1, 0], sqrt(2 * 0.25**2), -np.cos(pi / 4)), + ([0.5, 0.5, 0.5], [1, 0, 0], 0.5, 1), + ([-0.5, 0.5, 0.5], [1, 0, 0], -0.5, 1), + ([0.75, 0.5, 0.5], [1, 1, 0], sqrt(2 * 0.25**2), np.cos(pi / 4)), ], ) def test_sdf_inner(simple_volume: Volume, point, vec, exp_dist, exp_dot): From f76b0ae20129ba756deb4831863cd36268f3ebe3 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Wed, 10 Jan 2024 17:26:47 +0000 Subject: [PATCH 12/17] Bump pyo3 version switch to pseudonormal containment check is now n_rays=0 rather than None, due to https://github.com/PyO3/pyo3/issues/3735 --- Cargo.lock | 50 ++++++++++++++++++++------------- Cargo.toml | 4 +-- python/ncollpyde/_ncollpyde.pyi | 2 +- python/ncollpyde/main.py | 5 ++-- src/interface.rs | 41 ++++++++------------------- src/utils.rs | 2 +- tests/test_bench.py | 2 +- tests/test_ncollpyde.py | 6 ---- 8 files changed, 49 insertions(+), 63 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index c17c328..c7f66c4 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -393,7 +393,7 @@ dependencies = [ "autocfg", "cfg-if", "crossbeam-utils", - "memoffset", + "memoffset 0.8.0", "scopeguard", ] @@ -643,9 +643,9 @@ dependencies = [ [[package]] name = "indoc" -version = "1.0.9" +version = "2.0.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bfa799dd5ed20a7e349f3b4639aa80d74549c81716d9ec4f994c9b5815598306" +checksum = "1e186cfbae8084e513daff4240b4797e342f988cecda4fb6c939150f96315fd8" [[package]] name = "is-terminal" @@ -801,6 +801,15 @@ dependencies = [ "autocfg", ] +[[package]] +name = "memoffset" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a634b1c61a95585bd15607c6ab0c4e5b226e695ff2800ba0cdccddf208c406c" +dependencies = [ + "autocfg", +] + [[package]] name = "nalgebra" version = "0.32.2" @@ -911,9 +920,9 @@ dependencies = [ [[package]] name = "numpy" -version = "0.18.0" +version = "0.20.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "96b0fee4571867d318651c24f4a570c3f18408cf95f16ccb576b3ce85496a46e" +checksum = "bef41cbb417ea83b30525259e30ccef6af39b31c240bda578889494c5392d331" dependencies = [ "libc", "ndarray", @@ -1099,14 +1108,14 @@ dependencies = [ [[package]] name = "pyo3" -version = "0.18.2" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cfb848f80438f926a9ebddf0a539ed6065434fd7aae03a89312a9821f81b8501" +checksum = "9a89dc7a5850d0e983be1ec2a463a171d20990487c3cfcd68b5363f1ee3d6fe0" dependencies = [ "cfg-if", "indoc", "libc", - "memoffset", + "memoffset 0.9.0", "parking_lot", "pyo3-build-config", "pyo3-ffi", @@ -1116,9 +1125,9 @@ dependencies = [ [[package]] name = "pyo3-build-config" -version = "0.18.2" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "98a42e7f42e917ce6664c832d5eee481ad514c98250c49e0b03b20593e2c7ed0" +checksum = "07426f0d8fe5a601f26293f300afd1a7b1ed5e78b2a705870c5f30893c5163be" dependencies = [ "once_cell", "target-lexicon", @@ -1126,9 +1135,9 @@ dependencies = [ [[package]] name = "pyo3-ffi" -version = "0.18.2" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a0707f0ab26826fe4ccd59b69106e9df5e12d097457c7b8f9c0fd1d2743eec4d" +checksum = "dbb7dec17e17766b46bca4f1a4215a85006b4c2ecde122076c562dd058da6cf1" dependencies = [ "libc", "pyo3-build-config", @@ -1136,25 +1145,26 @@ dependencies = [ [[package]] name = "pyo3-macros" -version = "0.18.2" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "978d18e61465ecd389e1f235ff5a467146dc4e3c3968b90d274fe73a5dd4a438" +checksum = "05f738b4e40d50b5711957f142878cfa0f28e054aa0ebdfc3fd137a843f74ed3" dependencies = [ "proc-macro2", "pyo3-macros-backend", "quote", - "syn 1.0.109", + "syn 2.0.11", ] [[package]] name = "pyo3-macros-backend" -version = "0.18.2" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8e0e1128f85ce3fca66e435e08aa2089a2689c1c48ce97803e13f63124058462" +checksum = "0fc910d4851847827daf9d6cdd4a823fbdaab5b8818325c5e97a86da79e8881f" dependencies = [ + "heck", "proc-macro2", "quote", - "syn 1.0.109", + "syn 2.0.11", ] [[package]] @@ -1627,9 +1637,9 @@ dependencies = [ [[package]] name = "unindent" -version = "0.1.11" +version = "0.2.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e1766d682d402817b5ac4490b3c3002d91dfa0d22812f341609f97b08757359c" +checksum = "c7de7d73e1754487cb58364ee906a499937a0dfabd86bcb980fa99ec8c8fa2ce" [[package]] name = "url" diff --git a/Cargo.toml b/Cargo.toml index 40200c8..07f38c0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,12 +4,12 @@ version = "0.19.0" edition = "2018" [dependencies] -pyo3 = { version = "0.18", features = ["extension-module", "abi3-py38"] } +pyo3 = { version = "0.20", features = ["extension-module", "abi3-py38"] } parry3d-f64 = { version = "0.13", features = ["dim3", "f64", "enhanced-determinism"] } rayon = "1.7" rand = "0.8" rand_pcg = "0.3" -numpy = "0.18" +numpy = "0.20" log = "0.4.20" # ndarray is a dependency of numpy. diff --git a/python/ncollpyde/_ncollpyde.pyi b/python/ncollpyde/_ncollpyde.pyi index 2ac2874..4d0f231 100644 --- a/python/ncollpyde/_ncollpyde.pyi +++ b/python/ncollpyde/_ncollpyde.pyi @@ -17,7 +17,7 @@ class TriMeshWrapper: self, points: Points, indices: Indices, n_rays: int, ray_seed: int ): ... def contains( - self, points: Points, n_rays: Optional[int], consensus: int, parallel: bool + self, points: Points, n_rays: int, consensus: int, parallel: bool ) -> npt.NDArray[np.bool_]: ... def distance( self, points: Points, signed: bool, parallel: bool diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index 553435c..74662a3 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -217,7 +217,7 @@ def contains( :param coords: :param n_rays: Optional[int] If None, use the maximum rays defined on construction. - If < 0, use signed distance strategy + If < 1, use signed distance strategy (more robust, but slower for many meshes). Otherwise, use this many meshes, up to the maximum defined on construction. :param consensus: Optional[int] @@ -235,7 +235,7 @@ def contains( coords = self._as_points(coords) if n_rays is None: n_rays = self.n_rays - elif n_rays < 0: + elif n_rays < 1: n_rays = None elif n_rays > self.n_rays: logger.warning( @@ -245,6 +245,7 @@ def contains( if n_rays is None: consensus = 1 + n_rays = 0 else: if consensus is None: consensus = n_rays // 2 + 1 diff --git a/src/interface.rs b/src/interface.rs index d1b4c56..addc5b2 100644 --- a/src/interface.rs +++ b/src/interface.rs @@ -1,4 +1,3 @@ -use std::fmt::Debug; use std::iter::repeat_with; use ndarray::{Array2, ArrayView1}; @@ -17,13 +16,13 @@ use crate::utils::{ points_cross_mesh, random_dir, sdf_inner, Precision, FLAGS, }; -fn vec_to_point(v: Vec) -> Point { - Point::new(v[0], v[1], v[2]) -} +// fn vec_to_point(v: Vec) -> Point { +// Point::new(v[0], v[1], v[2]) +// } -fn point_to_vec(p: &Point) -> Vec { - vec![p.x, p.y, p.z] -} +// fn point_to_vec(p: &Point) -> Vec { +// vec![p.x, p.y, p.z] +// } #[pyclass] pub struct TriMeshWrapper { @@ -96,20 +95,21 @@ impl TriMeshWrapper { collected.into_pyarray(py) } + // #[pyo3(signature = (points, n_rays=None, consensus=None, parallel=None))] pub fn contains<'py>( &self, py: Python<'py>, points: PyReadonlyArray2, - n_rays: Option, + n_rays: usize, consensus: usize, parallel: bool, ) -> &'py PyArray1 { let p_arr = points.as_array(); let zipped = Zip::from(p_arr.rows()); - let collected = if let Some(n) = n_rays { + let collected = if n_rays >= 1 { // use ray casting - let rays = &self.ray_directions[..n]; + let rays = &self.ray_directions[..n_rays]; let clos = |r: ArrayView1| { mesh_contains_point(&self.mesh, &Point::new(r[0], r[1], r[2]), rays, consensus) }; @@ -242,7 +242,7 @@ impl TriMeshWrapper { .as_array() .rows() .into_iter() - .zip(tgt_points.as_array().rows().into_iter()) + .zip(tgt_points.as_array().rows()) .zip(0_u64..) .for_each(|((src, tgt), i)| { if let Some((pt, is_bf)) = points_cross_mesh( @@ -265,25 +265,6 @@ impl TriMeshWrapper { ) } - #[deprecated] - pub fn intersections_many_threaded( - &self, - src_points: Vec>, - tgt_points: Vec>, - ) -> (Vec, Vec>, Vec) { - let (idxs, (intersections, is_backface)) = src_points - .into_par_iter() - .zip(tgt_points.into_par_iter()) - .enumerate() - .filter_map(|(i, (src, tgt))| { - points_cross_mesh(&self.mesh, &vec_to_point(src), &vec_to_point(tgt)) - .map(|o| (i as u64, (point_to_vec(&o.0), o.1))) - }) - .unzip(); - - (idxs, intersections, is_backface) - } - pub fn intersections_many_threaded2<'py>( &self, py: Python<'py>, diff --git a/src/utils.rs b/src/utils.rs index 5594583..00db46b 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -83,7 +83,7 @@ pub fn mesh_contains_point( return false; } } - return false; + false } pub fn mesh_contains_point_oriented(mesh: &TriMesh, point: &Point) -> bool { diff --git a/tests/test_bench.py b/tests/test_bench.py index 4dc5d1f..9f26011 100644 --- a/tests/test_bench.py +++ b/tests/test_bench.py @@ -221,7 +221,7 @@ def test_trimesh_contains(trimesh_volume, sample_points, expected, benchmark): @pytest.mark.benchmark(group=CONTAINS_SERIAL) -@pytest.mark.parametrize("n_rays", [0, 1, 2, 4, 8, 16]) +@pytest.mark.parametrize("n_rays", [1, 2, 4, 8, 16]) def test_ncollpyde_contains(mesh, n_rays, sample_points, expected, benchmark): ncollpyde_volume = Volume.from_meshio(mesh, n_rays=n_rays) ncollpyde_volume.threads = False diff --git a/tests/test_ncollpyde.py b/tests/test_ncollpyde.py index 2130375..c030506 100644 --- a/tests/test_ncollpyde.py +++ b/tests/test_ncollpyde.py @@ -57,12 +57,6 @@ def test_contains_results(volume: Volume): assert np.allclose(ray, psnorms) -def test_0_rays(mesh): - vol = Volume.from_meshio(mesh, n_rays=0) - points = [p for p, _ in points_expected] - assert np.array_equal(vol.contains(points), [False] * len(points)) - - def test_no_validation(mesh): triangles = mesh.cells_dict["triangle"] Volume(mesh.points, triangles, True) From aa9191d1db88472e64aee019d03955188eed0a38 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Wed, 10 Jan 2024 17:32:01 +0000 Subject: [PATCH 13/17] Remove deprecated intersections_many_threaded --- python/ncollpyde/_ncollpyde.pyi | 3 --- tests/test_bench.py | 5 ++++- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/python/ncollpyde/_ncollpyde.pyi b/python/ncollpyde/_ncollpyde.pyi index 4d0f231..3889bfc 100644 --- a/python/ncollpyde/_ncollpyde.pyi +++ b/python/ncollpyde/_ncollpyde.pyi @@ -32,9 +32,6 @@ class TriMeshWrapper: def intersections_many_threaded2( self, src_points: Points, tgt_points: Points ) -> Tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... - def intersections_many_threaded( - self, src_points: Points, tgt_points: Points - ) -> Tuple[List[int], Points, List[bool]]: ... def sdf_intersections( self, points: Points, vectors: Points, threaded: bool ) -> Tuple[npt.NDArray[np.float_], npt.NDArray[np.float_]]: ... diff --git a/tests/test_bench.py b/tests/test_bench.py index 9f26011..3b7deda 100644 --- a/tests/test_bench.py +++ b/tests/test_bench.py @@ -265,7 +265,10 @@ def test_ncollpyde_intersection(mesh, benchmark, threads): @pytest.mark.benchmark(group=INTERSECTION_PARALLEL_IMPL) @pytest.mark.parametrize( ("method_name",), - [("intersections_many_threaded",), ("intersections_many_threaded2",)], + [ + # ("intersections_many_threaded",), + ("intersections_many_threaded2",), + ], ) def test_ncollpyde_intersection_impls(mesh, benchmark, method_name): n_edges = 1_000 From 14089a3bb0fb27475807ef84bf8eb336c76ef5d8 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Wed, 10 Jan 2024 17:33:32 +0000 Subject: [PATCH 14/17] remove unused typings --- python/ncollpyde/_ncollpyde.pyi | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/ncollpyde/_ncollpyde.pyi b/python/ncollpyde/_ncollpyde.pyi index 3889bfc..8c03f41 100644 --- a/python/ncollpyde/_ncollpyde.pyi +++ b/python/ncollpyde/_ncollpyde.pyi @@ -1,4 +1,4 @@ -from typing import List, Tuple, Optional +from typing import Optional import numpy as np import numpy.typing as npt @@ -28,10 +28,10 @@ class TriMeshWrapper: def aabb(self) -> Points: ... def intersections_many( self, src_points: Points, tgt_points: Points - ) -> Tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... + ) -> tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... def intersections_many_threaded2( self, src_points: Points, tgt_points: Points - ) -> Tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... + ) -> tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... def sdf_intersections( self, points: Points, vectors: Points, threaded: bool - ) -> Tuple[npt.NDArray[np.float_], npt.NDArray[np.float_]]: ... + ) -> tuple[npt.NDArray[np.float_], npt.NDArray[np.float_]]: ... From aa031221fac0e107c446319a44b6ec28625cfd3f Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Wed, 10 Jan 2024 19:05:00 +0000 Subject: [PATCH 15/17] rename intersections_many_threaded2 --- python/ncollpyde/_ncollpyde.pyi | 2 +- python/ncollpyde/main.py | 2 +- src/interface.rs | 2 +- tests/test_bench.py | 4 ++-- tests/test_ncollpyde.py | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/python/ncollpyde/_ncollpyde.pyi b/python/ncollpyde/_ncollpyde.pyi index 8c03f41..3fe119f 100644 --- a/python/ncollpyde/_ncollpyde.pyi +++ b/python/ncollpyde/_ncollpyde.pyi @@ -29,7 +29,7 @@ class TriMeshWrapper: def intersections_many( self, src_points: Points, tgt_points: Points ) -> tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... - def intersections_many_threaded2( + def intersections_many_threaded( self, src_points: Points, tgt_points: Points ) -> tuple[npt.NDArray[np.uint64], Points, npt.NDArray[np.bool_]]: ... def sdf_intersections( diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index 74662a3..54efe1e 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -350,7 +350,7 @@ def intersections( src, tgt = self._validate_points(src_points, tgt_points) if self._interpret_threads(threads): - return self._impl.intersections_many_threaded2(src, tgt) + return self._impl.intersections_many_threaded(src, tgt) else: return self._impl.intersections_many(src, tgt) diff --git a/src/interface.rs b/src/interface.rs index addc5b2..81d39f8 100644 --- a/src/interface.rs +++ b/src/interface.rs @@ -265,7 +265,7 @@ impl TriMeshWrapper { ) } - pub fn intersections_many_threaded2<'py>( + pub fn intersections_many_threaded<'py>( &self, py: Python<'py>, src_points: PyReadonlyArray2, diff --git a/tests/test_bench.py b/tests/test_bench.py index 3b7deda..2171c34 100644 --- a/tests/test_bench.py +++ b/tests/test_bench.py @@ -266,8 +266,8 @@ def test_ncollpyde_intersection(mesh, benchmark, threads): @pytest.mark.parametrize( ("method_name",), [ - # ("intersections_many_threaded",), - ("intersections_many_threaded2",), + ("intersections_many_threaded",), + # ("intersections_many_threaded2",), ], ) def test_ncollpyde_intersection_impls(mesh, benchmark, method_name): diff --git a/tests/test_ncollpyde.py b/tests/test_ncollpyde.py index c030506..d4b9f3a 100644 --- a/tests/test_ncollpyde.py +++ b/tests/test_ncollpyde.py @@ -344,5 +344,5 @@ def test_intersections_impls(volume: Volume): serial = volume.intersections(starts, stops, threads=False) par = volume.intersections(starts, stops, threads=True) assert_intersection_results(serial, par) - par2 = volume._impl.intersections_many_threaded2(starts, stops) + par2 = volume._impl.intersections_many_threaded(starts, stops) assert_intersection_results(serial, par2) From 9a09434315bfb9c0b236395884c921af637334f8 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Tue, 16 Jan 2024 17:43:06 +0000 Subject: [PATCH 16/17] Update validation routines --- python/ncollpyde/__init__.py | 4 +-- python/ncollpyde/_ncollpyde.pyi | 7 +++- python/ncollpyde/main.py | 64 +++++++++++++++++++++++++++++---- src/interface.rs | 14 ++++++-- tests/conftest.py | 8 ++--- tests/test_ncollpyde.py | 54 ++++++++++++++-------------- 6 files changed, 107 insertions(+), 44 deletions(-) diff --git a/python/ncollpyde/__init__.py b/python/ncollpyde/__init__.py index 96614cc..04df5b1 100644 --- a/python/ncollpyde/__init__.py +++ b/python/ncollpyde/__init__.py @@ -11,7 +11,7 @@ from .main import INDEX # noqa: F401 from .main import N_CPUS # noqa: F401 from .main import PRECISION # noqa: F401 -from .main import Volume # noqa: F401 +from .main import Volume, Validation from .main import configure_threadpool # noqa: F401 from ._ncollpyde import n_threads # noqa: F401 from ._ncollpyde import _version @@ -19,4 +19,4 @@ __version__ = _version() __version_info__ = tuple(int(n) for n in __version__.split("-")[0].split(".")) -__all__ = ["Volume"] +__all__ = ["Volume", "Validation"] diff --git a/python/ncollpyde/_ncollpyde.pyi b/python/ncollpyde/_ncollpyde.pyi index 3fe119f..a58cd90 100644 --- a/python/ncollpyde/_ncollpyde.pyi +++ b/python/ncollpyde/_ncollpyde.pyi @@ -14,7 +14,12 @@ Indices = npt.NDArray[np.uint32] class TriMeshWrapper: def __init__( - self, points: Points, indices: Indices, n_rays: int, ray_seed: int + self, + points: Points, + indices: Indices, + n_rays: int, + ray_seed: int, + validate: int, ): ... def contains( self, points: Points, n_rays: int, consensus: int, parallel: bool diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index 54efe1e..1cc9ca8 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -3,6 +3,7 @@ import warnings from multiprocessing import cpu_count from typing import TYPE_CHECKING, Optional, Tuple, Union, List +from enum import IntFlag, auto import numpy as np from numpy.typing import ArrayLike, NDArray @@ -29,6 +30,41 @@ INDEX = np.dtype(_index()) +class Validation(IntFlag): + """Enum representing the different validations which can be applied. + + Combine with `|`. + Must contain ORIENTED (see `Validation.minimum()`). + """ + HALF_EDGE_TOPOLOGY = 1 + CONNECTED_COMPONENTS = 2 + DELETE_BAD_TOPOLOGY_TRIANGLES = 4 + ORIENTED = 8 + MERGE_DUPLICATE_VERTICES = 16 + MERGE_DEGENERATE_TRIANGLES = 32 + MERGE_DUPLICATE_TRIANGLES = 64 + + @classmethod + def minimum(cls): + return cls.ORIENTED + + @classmethod + def default(cls): + return cls.all() + + @classmethod + def all(cls): + return ( + cls.HALF_EDGE_TOPOLOGY + | cls.CONNECTED_COMPONENTS + | cls.DELETE_BAD_TOPOLOGY_TRIANGLES + | cls.ORIENTED + | cls.MERGE_DUPLICATE_VERTICES + | cls.MERGE_DEGENERATE_TRIANGLES + | cls.MERGE_DUPLICATE_TRIANGLES + ) + + def configure_threadpool(n_threads: Optional[int], name_prefix: Optional[str]): """Configure the thread pool used for parallelisation. @@ -85,7 +121,7 @@ def __init__( self, vertices: ArrayLike, triangles: ArrayLike, - validate=False, + validate: Validation = Validation.default(), threads: Optional[bool] = None, n_rays=DEFAULT_RAYS, ray_seed=DEFAULT_SEED, @@ -114,10 +150,24 @@ def __init__( """ vert = np.asarray(vertices, self.dtype) if len(vert) > np.iinfo(INDEX).max: - raise ValueError(f"Cannot represent {len(vert)} vertices with {INDEX}") + raise ValueError( + f"Cannot represent {len(vert)} vertices with {INDEX}" + ) tri = np.asarray(triangles, INDEX) - if validate: - vert, tri = self._validate(vert, tri) + if isinstance(validate, bool): + warnings.warn( + "`validate: bool` should be replaced by a Validation enum " + "controlling validation by this library. " + "The previous behaviour of validating using the " + "external trimesh library may be removed in future.", + DeprecationWarning, + ) + if validate: + vert, tri = self._validate(vert, tri) + validate = Validation.all() + else: + validate = Validation.minimum() + self.threads = self._interpret_threads(threads) if ray_seed is None: logger.warning( @@ -129,7 +179,7 @@ def __init__( self.n_rays = int(n_rays) inner_rays = 0 if self.n_rays < 0 else self.n_rays - self._impl = TriMeshWrapper(vert, tri, inner_rays, ray_seed) + self._impl = TriMeshWrapper(vert, tri, inner_rays, ray_seed, validate) def _validate( self, vertices: np.ndarray, triangles: np.ndarray @@ -358,7 +408,7 @@ def intersections( def from_meshio( cls, mesh: "meshio.Mesh", - validate=False, + validate: Validation = Validation.default(), threads=None, n_rays=DEFAULT_RAYS, ray_seed=DEFAULT_SEED, @@ -367,7 +417,7 @@ def from_meshio( Convenience function for instantiating a Volume from a meshio Mesh. :param mesh: meshio Mesh whose only cells are triangles. - :param validate: as passed to ``__init__``, defaults to False + :param validate: as passed to ``__init__``, defaults to Validation.default() :param threads: as passed to ``__init__``, defaults to None :param n_rays: as passed to ``__init__``, defaults to 3 :param ray_seed: as passed to ``__init__``, defaults to None (random) diff --git a/src/interface.rs b/src/interface.rs index 81d39f8..5fee118 100644 --- a/src/interface.rs +++ b/src/interface.rs @@ -4,7 +4,7 @@ use ndarray::{Array2, ArrayView1}; use numpy::ndarray::{Array, Zip}; use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray2}; use parry3d_f64::math::{Point, Vector}; -use parry3d_f64::shape::TriMesh; +use parry3d_f64::shape::{TriMesh, TriMeshFlags}; use pyo3::exceptions::{PyRuntimeError, PyValueError}; use pyo3::prelude::*; use rand::SeedableRng; @@ -13,7 +13,7 @@ use rayon::{prelude::*, ThreadPoolBuilder}; use crate::utils::{ aabb_diag, dist_from_mesh, mesh_contains_point, mesh_contains_point_oriented, - points_cross_mesh, random_dir, sdf_inner, Precision, FLAGS, + points_cross_mesh, random_dir, sdf_inner, Precision, }; // fn vec_to_point(v: Vec) -> Point { @@ -38,6 +38,7 @@ impl TriMeshWrapper { indices: PyReadonlyArray2, n_rays: usize, ray_seed: u64, + validate: u8, ) -> PyResult { let points2 = points .as_array() @@ -53,7 +54,14 @@ impl TriMeshWrapper { .collect(); let mut mesh = TriMesh::new(points2, indices2); - mesh.set_flags(FLAGS) + let flags = TriMeshFlags::from_bits(validate) + .ok_or(PyValueError::new_err("Invalid `validate` enum"))?; + if !flags.contains(TriMeshFlags::ORIENTED) { + return Err(PyValueError::new_err( + "`validate` enum must contain ORIENTED", + )); + } + mesh.set_flags(flags) .map_err(|e| PyValueError::new_err(format!("Invalid mesh topology: {e}")))?; if n_rays > 0 { diff --git a/tests/conftest.py b/tests/conftest.py index 9cd3fb2..365ca99 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -3,7 +3,7 @@ import meshio import pytest -from ncollpyde import Volume +from ncollpyde import Volume, Validation test_dir = Path(__file__).resolve().parent project_dir = test_dir.parent @@ -17,7 +17,7 @@ def mesh(): @pytest.fixture def volume(mesh): - return Volume.from_meshio(mesh, validate=True) + return Volume.from_meshio(mesh, validate=Validation.all()) @pytest.fixture @@ -27,11 +27,11 @@ def simple_mesh(): @pytest.fixture def simple_volume(simple_mesh): - return Volume.from_meshio(simple_mesh, validate=True) + return Volume.from_meshio(simple_mesh, validate=Validation.all()) @pytest.fixture def sez_right(): return Volume.from_meshio( - meshio.read(str(mesh_dir / "SEZ_right.stl")), validate=True + meshio.read(str(mesh_dir / "SEZ_right.stl")), validate=Validation.all() ) diff --git a/tests/test_ncollpyde.py b/tests/test_ncollpyde.py index d4b9f3a..c8c0455 100644 --- a/tests/test_ncollpyde.py +++ b/tests/test_ncollpyde.py @@ -57,44 +57,44 @@ def test_contains_results(volume: Volume): assert np.allclose(ray, psnorms) -def test_no_validation(mesh): - triangles = mesh.cells_dict["triangle"] - Volume(mesh.points, triangles, True) +# def test_no_validation(mesh): +# triangles = mesh.cells_dict["triangle"] +# Volume(mesh.points, triangles, True) -@pytest.mark.skipif(not trimesh, reason="Requires trimesh") -def test_can_repair_hole(mesh): - triangles = mesh.cells_dict["triangle"] - triangles = triangles[:-1] - Volume(mesh.points, triangles, True) +# @pytest.mark.skipif(not trimesh, reason="Requires trimesh") +# def test_can_repair_hole(mesh): +# triangles = mesh.cells_dict["triangle"] +# triangles = triangles[:-1] +# Volume(mesh.points, triangles, True) -@pytest.mark.skipif(not trimesh, reason="Requires trimesh") -def test_can_repair_inversion(mesh): - triangles = mesh.cells_dict["triangle"] - triangles[-1] = triangles[-1, ::-1] - Volume(mesh.points, triangles, True) +# @pytest.mark.skipif(not trimesh, reason="Requires trimesh") +# def test_can_repair_inversion(mesh): +# triangles = mesh.cells_dict["triangle"] +# triangles[-1] = triangles[-1, ::-1] +# Volume(mesh.points, triangles, True) -@pytest.mark.skipif(not trimesh, reason="Requires trimesh") -def test_can_repair_inversions(mesh): - triangles = mesh.cells_dict["triangle"] - triangles = triangles[:, ::-1] - Volume(mesh.points, triangles, True) +# @pytest.mark.skipif(not trimesh, reason="Requires trimesh") +# def test_can_repair_inversions(mesh): +# triangles = mesh.cells_dict["triangle"] +# triangles = triangles[:, ::-1] +# Volume(mesh.points, triangles, True) -@pytest.mark.skipif(not trimesh, reason="Requires trimesh") -def test_inversions_repaired(simple_mesh): - center = [0.5, 0.5, 0.5] +# @pytest.mark.skipif(not trimesh, reason="Requires trimesh") +# def test_inversions_repaired(simple_mesh): +# center = [0.5, 0.5, 0.5] - orig_points = simple_mesh.points - orig_triangles = simple_mesh.cells_dict["triangle"] - assert center in Volume(orig_points, orig_triangles) +# orig_points = simple_mesh.points +# orig_triangles = simple_mesh.cells_dict["triangle"] +# assert center in Volume(orig_points, orig_triangles) - inv_triangles = orig_triangles[:, ::-1] - assert center not in Volume(orig_points, inv_triangles) +# inv_triangles = orig_triangles[:, ::-1] +# assert center not in Volume(orig_points, inv_triangles) - assert center in Volume(orig_points, inv_triangles, validate=True) +# assert center in Volume(orig_points, inv_triangles, validate=True) def test_points(mesh): From 743ae4565b199adc420b8b5edcfb24ca86c02ef3 Mon Sep 17 00:00:00 2001 From: Chris Barnes Date: Tue, 16 Jan 2024 17:45:40 +0000 Subject: [PATCH 17/17] fmt --- python/ncollpyde/main.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/python/ncollpyde/main.py b/python/ncollpyde/main.py index 1cc9ca8..d812aaa 100644 --- a/python/ncollpyde/main.py +++ b/python/ncollpyde/main.py @@ -3,7 +3,7 @@ import warnings from multiprocessing import cpu_count from typing import TYPE_CHECKING, Optional, Tuple, Union, List -from enum import IntFlag, auto +from enum import IntFlag import numpy as np from numpy.typing import ArrayLike, NDArray @@ -36,6 +36,7 @@ class Validation(IntFlag): Combine with `|`. Must contain ORIENTED (see `Validation.minimum()`). """ + HALF_EDGE_TOPOLOGY = 1 CONNECTED_COMPONENTS = 2 DELETE_BAD_TOPOLOGY_TRIANGLES = 4 @@ -150,9 +151,7 @@ def __init__( """ vert = np.asarray(vertices, self.dtype) if len(vert) > np.iinfo(INDEX).max: - raise ValueError( - f"Cannot represent {len(vert)} vertices with {INDEX}" - ) + raise ValueError(f"Cannot represent {len(vert)} vertices with {INDEX}") tri = np.asarray(triangles, INDEX) if isinstance(validate, bool): warnings.warn(