diff --git a/ncollpyde/__init__.py b/ncollpyde/__init__.py index 3c793bb..287e1e4 100644 --- a/ncollpyde/__init__.py +++ b/ncollpyde/__init__.py @@ -9,10 +9,9 @@ from .main import DEFAULT_SEED # noqa: F401 from .main import DEFAULT_THREADS # noqa: F401 from .main import INDEX # noqa: F401 -from .main import N_CPUS # noqa: F401 +from .main import N_THREADS # noqa: F401 from .main import PRECISION # noqa: F401 from .main import Volume # noqa: F401 -from .ncollpyde import n_threads # noqa: F401 from .ncollpyde import _version __version__ = _version() diff --git a/ncollpyde/main.py b/ncollpyde/main.py index 0e5f932..8b63cb6 100644 --- a/ncollpyde/main.py +++ b/ncollpyde/main.py @@ -1,7 +1,6 @@ import logging import random import warnings -from multiprocessing import cpu_count from typing import TYPE_CHECKING, Optional, Tuple, Union import numpy as np @@ -12,7 +11,7 @@ except ImportError: trimesh = None -from .ncollpyde import TriMeshWrapper, _index, _precision +from .ncollpyde import TriMeshWrapper, _index, _precision, n_threads if TYPE_CHECKING: import meshio @@ -20,7 +19,7 @@ logger = logging.getLogger(__name__) -N_CPUS = cpu_count() +N_THREADS = n_threads() DEFAULT_THREADS = True DEFAULT_RAYS = 3 DEFAULT_SEED = 1991 @@ -85,6 +84,10 @@ def __init__( - one ray reports that the point is external. :param ray_seed: int >=0 (default {DEFAULT_SEED}), used for generating rays. If None, use a random seed. + :param n_rays_inside: optional int (default None), used for ray consensus. + Number of rays which must hit mesh backfaces + for point to be considered internal. + If None, use n_rays. """ vert = np.asarray(vertices, self.dtype) if len(vert) > np.iinfo(INDEX).max: @@ -100,7 +103,18 @@ def __init__( ) ray_seed = random.randrange(0, 2**64) - self._impl = TriMeshWrapper(vert, tri, int(n_rays), ray_seed) + n_rays = max(int(n_rays), 0) + if n_rays < 1: + logger.warning( + "<1 ray used; points will only be considered internal " + "if they lie on the mesh shell" + ) + if not n_rays % 2: + logger.warning( + "Even number of rays used; odd numbers are preferred to break ties" + ) + + self._impl = TriMeshWrapper(vert, tri, int(n_rays), ray_seed, None) def _validate( self, vertices: np.ndarray, triangles: np.ndarray diff --git a/ncollpyde/ncollpyde.pyi b/ncollpyde/ncollpyde.pyi index 785837a..18876b1 100644 --- a/ncollpyde/ncollpyde.pyi +++ b/ncollpyde/ncollpyde.pyi @@ -1,4 +1,4 @@ -from typing import List, Tuple +from typing import List, Optional, Tuple import numpy as np import numpy.typing as npt @@ -13,7 +13,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, + n_rays_inside: Optional[int], ): ... def contains(self, points: Points, parallel: bool) -> npt.NDArray[np.bool_]: ... def distance( diff --git a/src/interface.rs b/src/interface.rs index 598ec47..892012e 100644 --- a/src/interface.rs +++ b/src/interface.rs @@ -33,6 +33,7 @@ fn vector_to_vec(v: &Vector) -> Vec pub struct TriMeshWrapper { mesh: TriMesh, ray_directions: Vec>, + n_rays_inside: usize, } #[cfg(not(test))] @@ -44,6 +45,7 @@ impl TriMeshWrapper { indices: PyReadonlyArray2, n_rays: usize, ray_seed: u64, + n_rays_inside: Option, ) -> Self { let points2 = points .as_array() @@ -59,6 +61,8 @@ impl TriMeshWrapper { .collect(); let mesh = TriMesh::new(points2, indices2); + let actual_n_rays_inside = n_rays_inside.unwrap_or(n_rays); + if n_rays > 0 { let bsphere = mesh.local_bounding_sphere(); let len = bsphere.radius() * 2.0; @@ -70,11 +74,13 @@ impl TriMeshWrapper { ray_directions: repeat_with(|| random_dir(&mut rng, len)) .take(n_rays) .collect(), + n_rays_inside: actual_n_rays_inside, } } else { Self { mesh, ray_directions: Vec::default(), + n_rays_inside: actual_n_rays_inside, } } } @@ -94,12 +100,24 @@ impl TriMeshWrapper { 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) + dist_from_mesh( + &self.mesh, + &Point::new(v[0], v[1], v[2]), + rays, + self.n_rays_inside, + ) }) .into_pyarray(py) } else { Zip::from(points.as_array().rows()) - .map_collect(|v| dist_from_mesh(&self.mesh, &Point::new(v[0], v[1], v[2]), rays)) + .map_collect(|v| { + dist_from_mesh( + &self.mesh, + &Point::new(v[0], v[1], v[2]), + rays, + self.n_rays_inside, + ) + }) .into_pyarray(py) } } @@ -117,6 +135,7 @@ impl TriMeshWrapper { &self.mesh, &Point::new(r[0], r[1], r[2]), &self.ray_directions, + self.n_rays_inside, ) }) .into_pyarray(py) @@ -127,6 +146,7 @@ impl TriMeshWrapper { &self.mesh, &Point::new(r[0], r[1], r[2]), &self.ray_directions, + self.n_rays_inside, ) }) .into_pyarray(py) diff --git a/src/utils.rs b/src/utils.rs index b1e37a3..ab58f2f 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -5,6 +5,7 @@ use rand::Rng; pub type Precision = f64; +/// Random vector with given length pub fn random_dir(rng: &mut R, length: Precision) -> Vector { let unscaled: Vector = [ rng.gen::() - 0.5, @@ -15,6 +16,61 @@ pub fn random_dir(rng: &mut R, length: Precision) -> Vector { unscaled.normalize() * length } +// /// 3 orthonormal vectors describing a random basis +// pub fn random_basis(rng: &mut R) -> Vec> { +// let rand = random_dir(rng, 1.0); +// let rand2 = random_dir(rng, 1.0); +// let cross1 = rand.cross(&rand2); +// let cross2 = rand.cross(&cross1); + +// vec![rand, cross1, cross2] +// } + +// /// 6 vectors made up of a random orthonormal basis and each of those vectors * -1 +// /// If `interleave`, arranged as A,-A,B,-B,C,-C; otherwise A,B,C,-A,-B,-C. +// pub fn random_compass(rng: &mut R, interleave: bool) -> Vec> { +// let pos = random_basis(rng); + +// if interleave { +// vec![ +// pos[0], +// -pos[0], +// pos[1], +// -pos[1], +// pos[2], +// -pos[2], +// ] +// } else { +// vec![ +// pos[0], +// pos[1], +// pos[2], +// -pos[0], +// -pos[1], +// -pos[2], +// ] +// } +// } + +// /// Any number of random vectors, where 0-6, 6-12, 12-18 etc. are as produced by `random_compass`. +// pub fn random_rays(rng: &mut R, n: usize, interleave: bool) -> Vec> { +// let mut out = Vec::with_capacity(n); +// while n > 0 { +// let mut vecs = random_compass(rng, interleave); + +// if n > vecs.len() { +// out.append(&mut vecs); +// } else { +// out.append(&mut vecs.drain(..n).collect()) +// } +// } +// out +// } + +// pub fn random_rays_with_length(rng: &mut R, n: usize, interleave: bool, length: Precision) -> Vec> { +// random_rays(rng, n, interleave).into_iter().map(|v| v * length).collect() +// } + pub fn mesh_contains_point_ray( mesh: &TriMesh, point: &Point, @@ -37,6 +93,7 @@ pub fn mesh_contains_point( mesh: &TriMesh, point: &Point, ray_directions: &[Vector], + n_rays_inside: usize, ) -> bool { if !mesh.local_aabb().contains_local_point(point) { return false; @@ -47,13 +104,40 @@ pub fn mesh_contains_point( return true; } + if n_rays_inside == 0 { + return true; + } + if ray_directions.is_empty() { - false - } else { - ray_directions - .iter() - .all(|v| mesh_contains_point_ray(mesh, point, v)) + return false; } + + if n_rays_inside > ray_directions.len() { + return false; + } + + // counting both hits and misses allows short-circuiting in both directions + let n_rays_outside = ray_directions.len() - n_rays_inside; + let mut in_count = 0; + let mut out_count = 0; + + for ray in ray_directions { + let is_inside = mesh_contains_point_ray(mesh, point, ray); + if is_inside { + in_count += 1; + if in_count >= n_rays_inside { + return true; + } + } else { + out_count += 1; + if out_count >= n_rays_outside { + return false; + } + } + } + + // probably shouldn't happen + false } pub fn points_cross_mesh( @@ -68,10 +152,15 @@ pub fn points_cross_mesh( .map(|i| (ray.point_at(i.toi), mesh.is_backface(i.feature))) } -pub fn dist_from_mesh(mesh: &TriMesh, point: &Point, rays: Option<&[Vector]>) -> f64 { +pub fn dist_from_mesh( + mesh: &TriMesh, + point: &Point, + rays: Option<&[Vector]>, + n_rays_inside: usize, +) -> f64 { let mut dist = mesh.distance_to_point(&Isometry::identity(), point, true); if let Some(r) = rays { - if mesh_contains_point(mesh, point, r) { + if mesh_contains_point(mesh, point, r, n_rays_inside) { dist = -dist; } } @@ -279,7 +368,7 @@ mod tests { rays: Option<&[Vector]>, expected: Precision, ) { - assert_eq!(dist_from_mesh(mesh, point, rays), expected) + assert_eq!(dist_from_mesh(mesh, point, rays, 2), expected) } #[test] diff --git a/tests/test_ncollpyde.py b/tests/test_ncollpyde.py index 5fea1ee..88d3b1c 100644 --- a/tests/test_ncollpyde.py +++ b/tests/test_ncollpyde.py @@ -46,10 +46,10 @@ def test_many(mesh, threads): assert np.array_equal(vol.contains(points, threads=threads), expected) -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_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):