parry3d/transformation/voxelization/
voxel_set.rs

1// Rust port, with modifications, of https://github.com/kmammou/v-hacd/blob/master/src/VHACD_Lib/src/vhacdVolume.cpp
2// By Khaled Mamou
3//
4// # License of the original C++ code:
5// > Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com)
6// > All rights reserved.
7// >
8// >
9// > Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
10// >
11// > 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
12// >
13// > 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
14// >
15// > 3. The names of the contributors may not be used to endorse or promote products derived from this software without specific prior written permission.
16// >
17// > THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
18
19use super::{FillMode, VoxelizedVolume};
20use crate::bounding_volume::Aabb;
21use crate::math::{Matrix, Point, Real, Vector, DIM};
22use crate::transformation::vhacd::CutPlane;
23use alloc::sync::Arc;
24use alloc::{vec, vec::Vec};
25
26#[cfg(feature = "dim2")]
27type ConvexHull = Vec<Point<Real>>;
28#[cfg(feature = "dim3")]
29type ConvexHull = (Vec<Point<Real>>, Vec<[u32; DIM]>);
30
31/// A voxel.
32#[derive(Copy, Clone, Debug)]
33pub struct Voxel {
34    /// The integer coordinates of the voxel as part of the voxel grid.
35    pub coords: Point<u32>,
36    /// Is this voxel on the surface of the volume (i.e. not inside of it)?
37    pub is_on_surface: bool,
38    /// Range of indices (to be looked up into the `VoxelSet` primitive map)
39    /// of the primitives intersected by this voxel.
40    pub(crate) intersections_range: (usize, usize),
41}
42
43impl Default for Voxel {
44    fn default() -> Self {
45        Self {
46            coords: Point::origin(),
47            is_on_surface: false,
48            intersections_range: (0, 0),
49        }
50    }
51}
52
53/// A sparse set of voxels.
54///
55/// It only contains voxels that are considered as "full" after a voxelization.
56pub struct VoxelSet {
57    /// The 3D origin of this voxel-set.
58    pub origin: Point<Real>,
59    /// The scale factor between the voxel integer coordinates and their
60    /// actual float world-space coordinates.
61    pub scale: Real,
62    pub(crate) min_bb_voxels: Point<u32>,
63    pub(crate) max_bb_voxels: Point<u32>,
64    pub(crate) voxels: Vec<Voxel>,
65    pub(crate) intersections: Arc<Vec<u32>>,
66    pub(crate) primitive_classes: Arc<Vec<u32>>,
67}
68
69impl Default for VoxelSet {
70    fn default() -> Self {
71        Self::new()
72    }
73}
74
75impl VoxelSet {
76    /// Creates a new empty set of voxels.
77    pub fn new() -> Self {
78        Self {
79            origin: Point::origin(),
80            min_bb_voxels: Point::origin(),
81            max_bb_voxels: Vector::repeat(1).into(),
82            scale: 1.0,
83            voxels: Vec::new(),
84            intersections: Arc::new(Vec::new()),
85            primitive_classes: Arc::new(Vec::new()),
86        }
87    }
88
89    /// The volume of a single voxel of this voxel set.
90    #[cfg(feature = "dim2")]
91    pub fn voxel_volume(&self) -> Real {
92        self.scale * self.scale
93    }
94
95    /// The volume of a single voxel of this voxel set.
96    #[cfg(feature = "dim3")]
97    pub fn voxel_volume(&self) -> Real {
98        self.scale * self.scale * self.scale
99    }
100
101    /// Voxelizes the given shape described by its boundary:
102    /// a triangle mesh (in 3D) or polyline (in 2D).
103    ///
104    /// # Parameters
105    /// * `points` - The vertex buffer of the boundary of the shape to voxelize.
106    /// * `indices` - The index buffer of the boundary of the shape to voxelize.
107    /// * `voxel_size` - The size of each voxel.
108    /// * `fill_mode` - Controls what is being voxelized.
109    /// * `keep_voxel_to_primitives_map` - If set to `true` a map between the voxels
110    ///   and the primitives (3D triangles or 2D segments) it intersects will be computed.
111    pub fn with_voxel_size(
112        points: &[Point<Real>],
113        indices: &[[u32; DIM]],
114        voxel_size: Real,
115        fill_mode: FillMode,
116        keep_voxel_to_primitives_map: bool,
117    ) -> Self {
118        VoxelizedVolume::with_voxel_size(
119            points,
120            indices,
121            voxel_size,
122            fill_mode,
123            keep_voxel_to_primitives_map,
124        )
125        .into()
126    }
127
128    /// Voxelizes the given shape described by its boundary:
129    /// a triangle mesh (in 3D) or polyline (in 2D).
130    ///
131    /// # Parameters
132    /// * `points` - The vertex buffer of the boundary of the shape to voxelize.
133    /// * `indices` - The index buffer of the boundary of the shape to voxelize.
134    /// * `resolution` - Controls the number of subdivision done along each axis. This number
135    ///   is the number of subdivisions along the axis where the input shape has the largest extent.
136    ///   The other dimensions will have a different automatically-determined resolution (in order to
137    ///   keep the voxels cubic).
138    /// * `fill_mode` - Controls what is being voxelized.
139    /// * `keep_voxel_to_primitives_map` - If set to `true` a map between the voxels
140    ///   and the primitives (3D triangles or 2D segments) it intersects will be computed.
141    pub fn voxelize(
142        points: &[Point<Real>],
143        indices: &[[u32; DIM]],
144        resolution: u32,
145        fill_mode: FillMode,
146        keep_voxel_to_primitives_map: bool,
147    ) -> Self {
148        VoxelizedVolume::voxelize(
149            points,
150            indices,
151            resolution,
152            fill_mode,
153            keep_voxel_to_primitives_map,
154        )
155        .into()
156    }
157
158    /// The minimal coordinates of the integer bounding-box of the voxels in this set.
159    pub fn min_bb_voxels(&self) -> Point<u32> {
160        self.min_bb_voxels
161    }
162
163    /// The maximal coordinates of the integer bounding-box of the voxels in this set.
164    pub fn max_bb_voxels(&self) -> Point<u32> {
165        self.max_bb_voxels
166    }
167
168    /// Computes the total volume of the voxels contained by this set.
169    pub fn compute_volume(&self) -> Real {
170        self.voxel_volume() * self.voxels.len() as Real
171    }
172
173    /// Gets the coordinates of the center of the given voxel.
174    pub fn get_voxel_point(&self, voxel: &Voxel) -> Point<Real> {
175        self.get_point(na::convert(voxel.coords))
176    }
177
178    pub(crate) fn get_point(&self, voxel: Point<Real>) -> Point<Real> {
179        self.origin + voxel.coords * self.scale
180    }
181
182    /// Does this voxel not contain any element?
183    pub fn is_empty(&self) -> bool {
184        self.len() == 0
185    }
186
187    /// The number of voxels in this set.
188    pub fn len(&self) -> usize {
189        self.voxels.len()
190    }
191
192    /// The set of voxels.
193    pub fn voxels(&self) -> &[Voxel] {
194        &self.voxels
195    }
196
197    /// Update the bounding box of this voxel set.
198    pub fn compute_bb(&mut self) {
199        let num_voxels = self.voxels.len();
200
201        if num_voxels == 0 {
202            return;
203        }
204
205        self.min_bb_voxels = self.voxels[0].coords;
206        self.max_bb_voxels = self.voxels[0].coords;
207
208        for p in 0..num_voxels {
209            self.min_bb_voxels = self.min_bb_voxels.inf(&self.voxels[p].coords);
210            self.max_bb_voxels = self.max_bb_voxels.sup(&self.voxels[p].coords);
211        }
212    }
213
214    // We have these cfg because we need to
215    // use the correct return type. We could just
216    // return ConvexHull but that would expose though
217    // the API a type alias that isn't really worth
218    // existing.
219    /// Compute the convex-hull of this voxel set after cutting each voxel
220    /// by the primitives (3D triangle or 2D segments) it intersects.
221    ///
222    /// This will panic if this `VoxelSet` was created with `keep_voxel_to_primitives_map = false`.
223    #[cfg(feature = "dim2")]
224    pub fn compute_exact_convex_hull(
225        &self,
226        points: &[Point<Real>],
227        indices: &[[u32; DIM]],
228    ) -> Vec<Point<Real>> {
229        self.do_compute_exact_convex_hull(points, indices)
230    }
231
232    /// Compute the convex-hull of this voxel set after cutting each voxel
233    /// by the primitives (3D triangle or 2D segments) it intersects.
234    ///
235    /// This will panic if this `VoxelSet` was created with `keep_voxel_to_primitives_map = false`.
236    #[cfg(feature = "dim3")]
237    pub fn compute_exact_convex_hull(
238        &self,
239        points: &[Point<Real>],
240        indices: &[[u32; DIM]],
241    ) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
242        self.do_compute_exact_convex_hull(points, indices)
243    }
244
245    fn do_compute_exact_convex_hull(
246        &self,
247        points: &[Point<Real>],
248        indices: &[[u32; DIM]],
249    ) -> ConvexHull {
250        assert!(!self.intersections.is_empty(),
251                "Cannot compute exact convex hull without voxel-to-primitives-map. Consider passing voxel_to_primitives_map = true to the voxelizer.");
252        let mut surface_points = Vec::new();
253        #[cfg(feature = "dim3")]
254        let (mut polygon, mut workspace) = (Vec::new(), Vec::new());
255        let mut pushed_points = vec![false; points.len()];
256
257        // Grab all the points.
258        for voxel in self.voxels.iter().filter(|v| v.is_on_surface) {
259            let intersections =
260                &self.intersections[voxel.intersections_range.0..voxel.intersections_range.1];
261            for prim_id in intersections {
262                let ia = indices[*prim_id as usize][0] as usize;
263                let ib = indices[*prim_id as usize][1] as usize;
264                #[cfg(feature = "dim3")]
265                let ic = indices[*prim_id as usize][2] as usize;
266
267                // If the primitives have been classified by VHACD, we know that:
268                // - A class equal to Some(u32::MAX) means that the primitives intersects multiple
269                //   convex parts, so we need to split it.
270                // - A class equal to None means that we did not compute any classes (so we
271                //   must assume that each triangle have to be split since it may intersect
272                //   multiple parts.
273                // - A class different from `None` and `Some(u32::MAX)` means that the triangle is
274                //   included in only one convex part. So instead of cutting it, just push the whole
275                //   triangle once.
276                let prim_class = self.primitive_classes.get(*prim_id as usize).copied();
277                if prim_class == Some(u32::MAX) || prim_class.is_none() {
278                    let aabb_center =
279                        self.origin + voxel.coords.coords.map(|k| k as Real) * self.scale;
280                    let aabb =
281                        Aabb::from_half_extents(aabb_center, Vector::repeat(self.scale / 2.0));
282
283                    #[cfg(feature = "dim2")]
284                    if let Some(seg) = aabb.clip_segment(&points[ia], &points[ib]) {
285                        surface_points.push(seg.a);
286                        surface_points.push(seg.b);
287                    }
288
289                    #[cfg(feature = "dim3")]
290                    {
291                        polygon.clear();
292                        polygon.extend_from_slice(&[points[ia], points[ib], points[ic]]);
293                        aabb.clip_polygon_with_workspace(&mut polygon, &mut workspace);
294                        surface_points.append(&mut polygon);
295                    }
296                } else {
297                    // We know this triangle is only contained by
298                    // one voxel set, i.e., `self`. So we don't
299                    // need to cut it.
300                    //
301                    // Because one triangle may intersect multiple voxels contained by
302                    // the same convex part, we only push vertices we have not pushed
303                    // so far in order to avoid some useless duplicate points (duplicate
304                    // points are OK as far as convex hull computation is concerned, but
305                    // they imply some redundant computations).
306                    let mut push_pt = |i: usize| {
307                        if !pushed_points[i] {
308                            surface_points.push(points[i]);
309                            pushed_points[i] = true;
310                        }
311                    };
312
313                    push_pt(ia);
314                    push_pt(ib);
315                    #[cfg(feature = "dim3")]
316                    push_pt(ic);
317                }
318            }
319
320            if intersections.is_empty() {
321                self.map_voxel_points(voxel, |p| surface_points.push(p));
322            }
323        }
324
325        // Compute the convex-hull.
326        convex_hull(&surface_points)
327    }
328
329    /// Computes the intersections between all the voxels of this voxel set,
330    /// and all the primitives (triangle or segments) it intersected (as per
331    /// the voxel-to-primitives-map computed during voxelization).
332    ///
333    /// Panics if the voxelization was performed without setting the parameter
334    /// `voxel_to_primitives_map = true`.
335    pub fn compute_primitive_intersections(
336        &self,
337        points: &[Point<Real>],
338        indices: &[[u32; DIM]],
339    ) -> Vec<Point<Real>> {
340        assert!(!self.intersections.is_empty(),
341                "Cannot compute primitive intersections voxel-to-primitives-map. Consider passing voxel_to_primitives_map = true to the voxelizer.");
342        let mut surface_points = Vec::new();
343        #[cfg(feature = "dim3")]
344        let (mut polygon, mut workspace) = (Vec::new(), Vec::new());
345
346        // Grab all the points.
347        for voxel in self.voxels.iter().filter(|v| v.is_on_surface) {
348            let intersections =
349                &self.intersections[voxel.intersections_range.0..voxel.intersections_range.1];
350            for prim_id in intersections {
351                let aabb_center = self.origin + voxel.coords.coords.map(|k| k as Real) * self.scale;
352                let aabb = Aabb::from_half_extents(aabb_center, Vector::repeat(self.scale / 2.0));
353
354                let pa = points[indices[*prim_id as usize][0] as usize];
355                let pb = points[indices[*prim_id as usize][1] as usize];
356                #[cfg(feature = "dim3")]
357                let pc = points[indices[*prim_id as usize][2] as usize];
358
359                #[cfg(feature = "dim2")]
360                if let Some(seg) = aabb.clip_segment(&pa, &pb) {
361                    surface_points.push(seg.a);
362                    surface_points.push(seg.b);
363                }
364
365                #[cfg(feature = "dim3")]
366                {
367                    workspace.clear();
368                    polygon.clear();
369                    polygon.extend_from_slice(&[pa, pb, pc]);
370                    aabb.clip_polygon_with_workspace(&mut polygon, &mut workspace);
371
372                    if polygon.len() > 2 {
373                        for i in 1..polygon.len() - 1 {
374                            surface_points.push(polygon[0]);
375                            surface_points.push(polygon[i]);
376                            surface_points.push(polygon[i + 1]);
377                        }
378                    }
379                }
380            }
381        }
382
383        surface_points
384    }
385
386    /// Compute the convex-hull of the voxels in this set.
387    ///
388    /// # Parameters
389    /// * `sampling` - The convex-hull computation will ignore `sampling` voxels at
390    ///   regular intervals. Useful to save some computation times if an exact result isn't need.
391    ///   Use `0` to make sure no voxel is being ignored.
392    #[cfg(feature = "dim2")]
393    pub fn compute_convex_hull(&self, sampling: u32) -> Vec<Point<Real>> {
394        let mut points = Vec::new();
395
396        // Grab all the points.
397        for voxel in self
398            .voxels
399            .iter()
400            .filter(|v| v.is_on_surface)
401            .step_by(sampling as usize)
402        {
403            self.map_voxel_points(voxel, |p| points.push(p));
404        }
405
406        // Compute the convex-hull.
407        convex_hull(&points)
408    }
409
410    /// Compute the convex-hull of the voxels in this set.
411    ///
412    /// # Parameters
413    /// * `sampling` - The convex-hull computation will ignore `sampling` voxels at
414    ///   regular intervals. Useful to save some computation times if an exact result isn't need.
415    ///   Use `0` to make sure no voxel is being ignored.
416    #[cfg(feature = "dim3")]
417    pub fn compute_convex_hull(&self, sampling: u32) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
418        let mut points = Vec::new();
419
420        // Grab all the points.
421        for voxel in self
422            .voxels
423            .iter()
424            .filter(|v| v.is_on_surface)
425            .step_by(sampling as usize)
426        {
427            self.map_voxel_points(voxel, |p| points.push(p));
428        }
429
430        // Compute the convex-hull.
431        convex_hull(&points)
432    }
433
434    /// Gets the vertices of the given voxel.
435    fn map_voxel_points(&self, voxel: &Voxel, mut f: impl FnMut(Point<Real>)) {
436        let ijk = voxel.coords.coords.map(|e| e as Real);
437
438        #[cfg(feature = "dim2")]
439        let shifts = [
440            Vector::new(-0.5, -0.5),
441            Vector::new(0.5, -0.5),
442            Vector::new(0.5, 0.5),
443            Vector::new(-0.5, 0.5),
444        ];
445
446        #[cfg(feature = "dim3")]
447        let shifts = [
448            Vector::new(-0.5, -0.5, -0.5),
449            Vector::new(0.5, -0.5, -0.5),
450            Vector::new(0.5, 0.5, -0.5),
451            Vector::new(-0.5, 0.5, -0.5),
452            Vector::new(-0.5, -0.5, 0.5),
453            Vector::new(0.5, -0.5, 0.5),
454            Vector::new(0.5, 0.5, 0.5),
455            Vector::new(-0.5, 0.5, 0.5),
456        ];
457
458        for shift in &shifts {
459            f(self.origin + (ijk + *shift) * self.scale)
460        }
461    }
462
463    pub(crate) fn intersect(
464        &self,
465        plane: &CutPlane,
466        positive_pts: &mut Vec<Point<Real>>,
467        negative_pts: &mut Vec<Point<Real>>,
468        sampling: u32,
469    ) {
470        let num_voxels = self.voxels.len();
471
472        if num_voxels == 0 {
473            return;
474        }
475
476        let d0 = self.scale;
477        let mut sp = 0;
478        let mut sn = 0;
479
480        for v in 0..num_voxels {
481            let voxel = self.voxels[v];
482            let pt = self.get_voxel_point(&voxel);
483            let d = plane.abc.dot(&pt.coords) + plane.d;
484
485            // if      (d >= 0.0 && d <= d0) positive_pts.push(pt);
486            // else if (d < 0.0 && -d <= d0) negative_pts.push(pt);
487
488            if d >= 0.0 {
489                if d <= d0 {
490                    self.map_voxel_points(&voxel, |p| positive_pts.push(p));
491                } else {
492                    sp += 1;
493
494                    if sp == sampling {
495                        self.map_voxel_points(&voxel, |p| positive_pts.push(p));
496                        sp = 0;
497                    }
498                }
499            } else if -d <= d0 {
500                self.map_voxel_points(&voxel, |p| negative_pts.push(p));
501            } else {
502                sn += 1;
503                if sn == sampling {
504                    self.map_voxel_points(&voxel, |p| negative_pts.push(p));
505                    sn = 0;
506                }
507            }
508        }
509    }
510
511    // Returns (negative_volume, positive_volume)
512    pub(crate) fn compute_clipped_volumes(&self, plane: &CutPlane) -> (Real, Real) {
513        if self.voxels.is_empty() {
514            return (0.0, 0.0);
515        }
516
517        let mut num_positive_voxels = 0;
518
519        for voxel in &self.voxels {
520            let pt = self.get_voxel_point(voxel);
521            let d = plane.abc.dot(&pt.coords) + plane.d;
522            num_positive_voxels += (d >= 0.0) as usize;
523        }
524
525        let num_negative_voxels = self.voxels.len() - num_positive_voxels;
526        let positive_volume = self.voxel_volume() * (num_positive_voxels as Real);
527        let negative_volume = self.voxel_volume() * (num_negative_voxels as Real);
528
529        (negative_volume, positive_volume)
530    }
531
532    // Set `on_surf` such that it contains only the voxel on surface contained by `self`.
533    pub(crate) fn select_on_surface(&self, on_surf: &mut VoxelSet) {
534        on_surf.origin = self.origin;
535        on_surf.voxels.clear();
536        on_surf.scale = self.scale;
537
538        for voxel in &self.voxels {
539            if voxel.is_on_surface {
540                on_surf.voxels.push(*voxel);
541            }
542        }
543    }
544
545    /// Splits this voxel set into two parts, depending on where the voxel center lies wrt. the given plane.
546    pub(crate) fn clip(
547        &self,
548        plane: &CutPlane,
549        positive_part: &mut VoxelSet,
550        negative_part: &mut VoxelSet,
551    ) {
552        let num_voxels = self.voxels.len();
553
554        if num_voxels == 0 {
555            return;
556        }
557
558        negative_part.origin = self.origin;
559        negative_part.voxels.clear();
560        negative_part.voxels.reserve(num_voxels);
561        negative_part.scale = self.scale;
562
563        positive_part.origin = self.origin;
564        positive_part.voxels.clear();
565        positive_part.voxels.reserve(num_voxels);
566        positive_part.scale = self.scale;
567
568        let d0 = self.scale;
569
570        for v in 0..num_voxels {
571            let mut voxel = self.voxels[v];
572            let pt = self.get_voxel_point(&voxel);
573            let d = plane.abc.dot(&pt.coords) + plane.d;
574
575            if d >= 0.0 {
576                if voxel.is_on_surface || d <= d0 {
577                    voxel.is_on_surface = true;
578                    positive_part.voxels.push(voxel);
579                } else {
580                    positive_part.voxels.push(voxel);
581                }
582            } else if voxel.is_on_surface || -d <= d0 {
583                voxel.is_on_surface = true;
584                negative_part.voxels.push(voxel);
585            } else {
586                negative_part.voxels.push(voxel);
587            }
588        }
589    }
590
591    /// Convert `self` into a mesh, including only the voxels on the surface or only the voxel
592    /// inside of the volume.
593    #[cfg(feature = "dim3")]
594    pub fn to_trimesh(
595        &self,
596        base_index: u32,
597        is_on_surface: bool,
598    ) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
599        let mut vertices = Vec::new();
600        let mut indices = Vec::new();
601
602        for voxel in &self.voxels {
603            if voxel.is_on_surface == is_on_surface {
604                self.map_voxel_points(voxel, |p| vertices.push(p));
605
606                indices.push([base_index, base_index + 2, base_index + 1]);
607                indices.push([base_index, base_index + 3, base_index + 2]);
608                indices.push([base_index + 4, base_index + 5, base_index + 6]);
609                indices.push([base_index + 4, base_index + 6, base_index + 7]);
610                indices.push([base_index + 7, base_index + 6, base_index + 2]);
611                indices.push([base_index + 7, base_index + 2, base_index + 3]);
612                indices.push([base_index + 4, base_index + 1, base_index + 5]);
613                indices.push([base_index + 4, base_index, base_index + 1]);
614                indices.push([base_index + 6, base_index + 5, base_index + 1]);
615                indices.push([base_index + 6, base_index + 1, base_index + 2]);
616                indices.push([base_index + 7, base_index, base_index + 4]);
617                indices.push([base_index + 7, base_index + 3, base_index]);
618            }
619        }
620
621        (vertices, indices)
622    }
623
624    pub(crate) fn compute_principal_axes(&self) -> Vector<Real> {
625        let num_voxels = self.voxels.len();
626        if num_voxels == 0 {
627            return Vector::zeros();
628        }
629
630        // TODO: find a way to reuse crate::utils::cov?
631        // The difficulty being that we need to iterate through the set of
632        // points twice. So passing an iterator to crate::utils::cov
633        // isn't really possible.
634        let mut center = Point::origin();
635        let denom = 1.0 / (num_voxels as Real);
636
637        for voxel in &self.voxels {
638            center += voxel.coords.map(|e| e as Real).coords * denom;
639        }
640
641        let mut cov_mat = Matrix::zeros();
642        for voxel in &self.voxels {
643            let xyz = voxel.coords.map(|e| e as Real) - center;
644            cov_mat.syger(denom, &xyz, &xyz, 1.0);
645        }
646
647        cov_mat.symmetric_eigenvalues()
648    }
649
650    pub(crate) fn compute_axes_aligned_clipping_planes(
651        &self,
652        downsampling: u32,
653        planes: &mut Vec<CutPlane>,
654    ) {
655        let min_v = self.min_bb_voxels();
656        let max_v = self.max_bb_voxels();
657
658        for dim in 0..DIM {
659            let i0 = min_v[dim];
660            let i1 = max_v[dim];
661
662            for i in (i0..=i1).step_by(downsampling as usize) {
663                let plane = CutPlane {
664                    abc: Vector::ith(dim, 1.0),
665                    axis: dim as u8,
666                    d: -(self.origin[dim] + (i as Real + 0.5) * self.scale),
667                    index: i,
668                };
669
670                planes.push(plane);
671            }
672        }
673    }
674}
675
676#[cfg(feature = "dim2")]
677fn convex_hull(vertices: &[Point<Real>]) -> Vec<Point<Real>> {
678    if vertices.len() > 1 {
679        crate::transformation::convex_hull(vertices)
680    } else {
681        Vec::new()
682    }
683}
684
685#[cfg(feature = "dim3")]
686fn convex_hull(vertices: &[Point<Real>]) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
687    if vertices.len() > 2 {
688        crate::transformation::convex_hull(vertices)
689    } else {
690        (Vec::new(), Vec::new())
691    }
692}
parry3d/transformation/voxelization/voxel_set.rs

parry3d/transformation/voxelization/
voxel_set.rs
