use super::{FillMode, VoxelizedVolume};
use crate::bounding_volume::AABB;
use crate::math::{Matrix, Point, Real, Vector, DIM};
use crate::transformation::vhacd::CutPlane;
use std::sync::Arc;
#[cfg(feature = "dim2")]
type ConvexHull = Vec<Point<Real>>;
#[cfg(feature = "dim3")]
type ConvexHull = (Vec<Point<Real>>, Vec<[u32; DIM]>);
#[derive(Copy, Clone, Debug)]
pub struct Voxel {
pub coords: Point<u32>,
pub is_on_surface: bool,
pub(crate) intersections_range: (usize, usize),
}
impl Default for Voxel {
fn default() -> Self {
Self {
coords: Point::origin(),
is_on_surface: false,
intersections_range: (0, 0),
}
}
}
pub struct VoxelSet {
pub origin: Point<Real>,
pub scale: Real,
pub(crate) min_bb_voxels: Point<u32>,
pub(crate) max_bb_voxels: Point<u32>,
pub(crate) voxels: Vec<Voxel>,
pub(crate) intersections: Arc<Vec<u32>>,
pub(crate) primitive_classes: Arc<Vec<u32>>,
}
impl VoxelSet {
pub fn new() -> Self {
Self {
origin: Point::origin(),
min_bb_voxels: Point::origin(),
max_bb_voxels: Vector::repeat(1).into(),
scale: 1.0,
voxels: Vec::new(),
intersections: Arc::new(Vec::new()),
primitive_classes: Arc::new(Vec::new()),
}
}
#[cfg(feature = "dim2")]
pub fn voxel_volume(&self) -> Real {
self.scale * self.scale
}
#[cfg(feature = "dim3")]
pub fn voxel_volume(&self) -> Real {
self.scale * self.scale * self.scale
}
pub fn voxelize(
points: &[Point<Real>],
indices: &[[u32; DIM]],
resolution: u32,
fill_mode: FillMode,
keep_voxel_to_primitives_map: bool,
) -> Self {
VoxelizedVolume::voxelize(
points,
indices,
resolution,
fill_mode,
keep_voxel_to_primitives_map,
)
.into()
}
pub fn min_bb_voxels(&self) -> Point<u32> {
self.min_bb_voxels
}
pub fn max_bb_voxels(&self) -> Point<u32> {
self.max_bb_voxels
}
pub fn compute_volume(&self) -> Real {
self.voxel_volume() * self.voxels.len() as Real
}
fn get_voxel_point(&self, voxel: &Voxel) -> Point<Real> {
self.get_point(na::convert(voxel.coords))
}
pub(crate) fn get_point(&self, voxel: Point<Real>) -> Point<Real> {
self.origin + voxel.coords * self.scale
}
pub fn len(&self) -> usize {
self.voxels.len()
}
pub fn voxels(&self) -> &[Voxel] {
&self.voxels
}
pub fn compute_bb(&mut self) {
let num_voxels = self.voxels.len();
if num_voxels == 0 {
return;
}
self.min_bb_voxels = self.voxels[0].coords;
self.max_bb_voxels = self.voxels[0].coords;
for p in 0..num_voxels {
self.min_bb_voxels = self.min_bb_voxels.inf(&self.voxels[p].coords);
self.max_bb_voxels = self.max_bb_voxels.sup(&self.voxels[p].coords);
}
}
#[cfg(feature = "dim2")]
pub fn compute_exact_convex_hull(
&self,
points: &[Point<Real>],
indices: &[[u32; DIM]],
) -> Vec<Point<Real>> {
self.do_compute_exact_convex_hull(points, indices)
}
#[cfg(feature = "dim3")]
pub fn compute_exact_convex_hull(
&self,
points: &[Point<Real>],
indices: &[[u32; DIM]],
) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
self.do_compute_exact_convex_hull(points, indices)
}
fn do_compute_exact_convex_hull(
&self,
points: &[Point<Real>],
indices: &[[u32; DIM]],
) -> ConvexHull {
assert!(!self.intersections.is_empty(),
"Cannot compute exact convex hull without voxel-to-primitives-map. Consider passing voxel_to_primitives_map = true to the voxelizer.");
let mut surface_points = Vec::new();
#[cfg(feature = "dim3")]
let (mut polygon, mut workspace) = (Vec::new(), Vec::new());
let mut pushed_points = vec![false; points.len()];
for voxel in self.voxels.iter().filter(|v| v.is_on_surface) {
let intersections =
&self.intersections[voxel.intersections_range.0..voxel.intersections_range.1];
for prim_id in intersections {
let ia = indices[*prim_id as usize][0] as usize;
let ib = indices[*prim_id as usize][1] as usize;
#[cfg(feature = "dim3")]
let ic = indices[*prim_id as usize][2] as usize;
let prim_class = self.primitive_classes.get(*prim_id as usize).copied();
if prim_class == Some(u32::MAX) || prim_class == None {
let aabb_center =
self.origin + voxel.coords.coords.map(|k| k as Real) * self.scale;
let aabb =
AABB::from_half_extents(aabb_center, Vector::repeat(self.scale / 2.0));
#[cfg(feature = "dim2")]
if let Some(seg) = aabb.clip_segment(&points[ia], &points[ib]) {
surface_points.push(seg.a);
surface_points.push(seg.b);
}
#[cfg(feature = "dim3")]
{
polygon.clear();
polygon.extend_from_slice(&[points[ia], points[ib], points[ic]]);
aabb.clip_polygon_with_workspace(&mut polygon, &mut workspace);
surface_points.append(&mut polygon);
}
} else {
let mut push_pt = |i: usize| {
if !pushed_points[i] {
surface_points.push(points[i]);
pushed_points[i] = true;
}
};
push_pt(ia);
push_pt(ib);
#[cfg(feature = "dim3")]
push_pt(ic);
}
}
if intersections.is_empty() {
self.map_voxel_points(voxel, |p| surface_points.push(p));
}
}
convex_hull(&surface_points)
}
pub fn compute_primitive_intersections(
&self,
points: &[Point<Real>],
indices: &[[u32; DIM]],
) -> Vec<Point<Real>> {
assert!(!self.intersections.is_empty(),
"Cannot compute primitive intersections voxel-to-primitives-map. Consider passing voxel_to_primitives_map = true to the voxelizer.");
let mut surface_points = Vec::new();
#[cfg(feature = "dim3")]
let (mut polygon, mut workspace) = (Vec::new(), Vec::new());
for voxel in self.voxels.iter().filter(|v| v.is_on_surface) {
let intersections =
&self.intersections[voxel.intersections_range.0..voxel.intersections_range.1];
for prim_id in intersections {
let aabb_center = self.origin + voxel.coords.coords.map(|k| k as Real) * self.scale;
let aabb = AABB::from_half_extents(aabb_center, Vector::repeat(self.scale / 2.0));
let pa = points[indices[*prim_id as usize][0] as usize];
let pb = points[indices[*prim_id as usize][1] as usize];
#[cfg(feature = "dim3")]
let pc = points[indices[*prim_id as usize][2] as usize];
#[cfg(feature = "dim2")]
if let Some(seg) = aabb.clip_segment(&pa, &pb) {
surface_points.push(seg.a);
surface_points.push(seg.b);
}
#[cfg(feature = "dim3")]
{
workspace.clear();
polygon.clear();
polygon.extend_from_slice(&[pa, pb, pc]);
aabb.clip_polygon_with_workspace(&mut polygon, &mut workspace);
if polygon.len() > 2 {
for i in 1..polygon.len() - 1 {
surface_points.push(polygon[0]);
surface_points.push(polygon[i]);
surface_points.push(polygon[i + 1]);
}
}
}
}
}
surface_points
}
#[cfg(feature = "dim2")]
pub fn compute_convex_hull(&self, sampling: u32) -> Vec<Point<Real>> {
let mut points = Vec::new();
for voxel in self
.voxels
.iter()
.filter(|v| v.is_on_surface)
.step_by(sampling as usize)
{
self.map_voxel_points(voxel, |p| points.push(p));
}
convex_hull(&points)
}
#[cfg(feature = "dim3")]
pub fn compute_convex_hull(&self, sampling: u32) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
let mut points = Vec::new();
for voxel in self
.voxels
.iter()
.filter(|v| v.is_on_surface)
.step_by(sampling as usize)
{
self.map_voxel_points(voxel, |p| points.push(p));
}
convex_hull(&points)
}
fn map_voxel_points(&self, voxel: &Voxel, mut f: impl FnMut(Point<Real>)) {
let ijk = voxel.coords.coords.map(|e| e as Real);
#[cfg(feature = "dim2")]
let shifts = [
Vector::new(-0.5, -0.5),
Vector::new(0.5, -0.5),
Vector::new(0.5, 0.5),
Vector::new(-0.5, 0.5),
];
#[cfg(feature = "dim3")]
let shifts = [
Vector::new(-0.5, -0.5, -0.5),
Vector::new(0.5, -0.5, -0.5),
Vector::new(0.5, 0.5, -0.5),
Vector::new(-0.5, 0.5, -0.5),
Vector::new(-0.5, -0.5, 0.5),
Vector::new(0.5, -0.5, 0.5),
Vector::new(0.5, 0.5, 0.5),
Vector::new(-0.5, 0.5, 0.5),
];
for shift in &shifts {
f(self.origin + (ijk + *shift) * self.scale)
}
}
pub(crate) fn intersect(
&self,
plane: &CutPlane,
positive_pts: &mut Vec<Point<Real>>,
negative_pts: &mut Vec<Point<Real>>,
sampling: u32,
) {
let num_voxels = self.voxels.len();
if num_voxels == 0 {
return;
}
let d0 = self.scale;
let mut sp = 0;
let mut sn = 0;
for v in 0..num_voxels {
let voxel = self.voxels[v];
let pt = self.get_voxel_point(&voxel);
let d = plane.abc.dot(&pt.coords) + plane.d;
if d >= 0.0 {
if d <= d0 {
self.map_voxel_points(&voxel, |p| positive_pts.push(p));
} else {
sp += 1;
if sp == sampling {
self.map_voxel_points(&voxel, |p| positive_pts.push(p));
sp = 0;
}
}
} else {
if -d <= d0 {
self.map_voxel_points(&voxel, |p| negative_pts.push(p));
} else {
sn += 1;
if sn == sampling {
self.map_voxel_points(&voxel, |p| negative_pts.push(p));
sn = 0;
}
}
}
}
}
pub(crate) fn compute_clipped_volumes(&self, plane: &CutPlane) -> (Real, Real) {
if self.voxels.is_empty() {
return (0.0, 0.0);
}
let mut num_positive_voxels = 0;
for voxel in &self.voxels {
let pt = self.get_voxel_point(voxel);
let d = plane.abc.dot(&pt.coords) + plane.d;
num_positive_voxels += (d >= 0.0) as usize;
}
let num_negative_voxels = self.voxels.len() - num_positive_voxels;
let positive_volume = self.voxel_volume() * (num_positive_voxels as Real);
let negative_volume = self.voxel_volume() * (num_negative_voxels as Real);
(negative_volume, positive_volume)
}
pub(crate) fn select_on_surface(&self, on_surf: &mut VoxelSet) {
on_surf.origin = self.origin;
on_surf.voxels.clear();
on_surf.scale = self.scale;
for voxel in &self.voxels {
if voxel.is_on_surface {
on_surf.voxels.push(*voxel);
}
}
}
pub(crate) fn clip(
&self,
plane: &CutPlane,
positive_part: &mut VoxelSet,
negative_part: &mut VoxelSet,
) {
let num_voxels = self.voxels.len();
if num_voxels == 0 {
return;
}
negative_part.origin = self.origin;
negative_part.voxels.clear();
negative_part.voxels.reserve(num_voxels);
negative_part.scale = self.scale;
positive_part.origin = self.origin;
positive_part.voxels.clear();
positive_part.voxels.reserve(num_voxels);
positive_part.scale = self.scale;
let d0 = self.scale;
for v in 0..num_voxels {
let mut voxel = self.voxels[v];
let pt = self.get_voxel_point(&voxel);
let d = plane.abc.dot(&pt.coords) + plane.d;
if d >= 0.0 {
if voxel.is_on_surface || d <= d0 {
voxel.is_on_surface = true;
positive_part.voxels.push(voxel);
} else {
positive_part.voxels.push(voxel);
}
} else {
if voxel.is_on_surface || -d <= d0 {
voxel.is_on_surface = true;
negative_part.voxels.push(voxel);
} else {
negative_part.voxels.push(voxel);
}
}
}
}
#[cfg(feature = "dim3")]
pub fn to_trimesh(
&self,
base_index: u32,
is_on_surface: bool,
) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
let mut vertices = Vec::new();
let mut indices = Vec::new();
for voxel in &self.voxels {
if voxel.is_on_surface == is_on_surface {
self.map_voxel_points(voxel, |p| vertices.push(p));
indices.push([base_index + 0, base_index + 2, base_index + 1]);
indices.push([base_index + 0, base_index + 3, base_index + 2]);
indices.push([base_index + 4, base_index + 5, base_index + 6]);
indices.push([base_index + 4, base_index + 6, base_index + 7]);
indices.push([base_index + 7, base_index + 6, base_index + 2]);
indices.push([base_index + 7, base_index + 2, base_index + 3]);
indices.push([base_index + 4, base_index + 1, base_index + 5]);
indices.push([base_index + 4, base_index + 0, base_index + 1]);
indices.push([base_index + 6, base_index + 5, base_index + 1]);
indices.push([base_index + 6, base_index + 1, base_index + 2]);
indices.push([base_index + 7, base_index + 0, base_index + 4]);
indices.push([base_index + 7, base_index + 3, base_index + 0]);
}
}
(vertices, indices)
}
pub(crate) fn compute_principal_axes(&self) -> Vector<Real> {
let num_voxels = self.voxels.len();
if num_voxels == 0 {
return Vector::zeros();
}
let mut center = Point::origin();
let denom = 1.0 / (num_voxels as Real);
for voxel in &self.voxels {
center += voxel.coords.map(|e| e as Real).coords * denom;
}
let mut cov_mat = Matrix::zeros();
for voxel in &self.voxels {
let xyz = voxel.coords.map(|e| e as Real) - center;
cov_mat.syger(denom, &xyz, &xyz, 1.0);
}
cov_mat.symmetric_eigenvalues()
}
}
#[cfg(feature = "dim2")]
fn convex_hull(vertices: &[Point<Real>]) -> Vec<Point<Real>> {
if vertices.len() > 1 {
crate::transformation::convex_hull(vertices)
} else {
Vec::new()
}
}
#[cfg(feature = "dim3")]
fn convex_hull(vertices: &[Point<Real>]) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
if vertices.len() > 2 {
crate::transformation::convex_hull(vertices)
} else {
(Vec::new(), Vec::new())
}
}