parry3d/transformation/vhacd/
vhacd.rs

1// Rust port, with modifications, of https://github.com/kmammou/v-hacd/blob/master/src/VHACD_Lib/src/VHACD.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 crate::math::{Point, Real, Vector, DIM};
20use crate::transformation::vhacd::VHACDParameters;
21use crate::transformation::voxelization::{VoxelSet, VoxelizedVolume};
22use alloc::sync::Arc;
23use alloc::vec::Vec;
24
25#[cfg(feature = "dim2")]
26type ConvexHull = Vec<Point<Real>>;
27#[cfg(feature = "dim3")]
28type ConvexHull = (Vec<Point<Real>>, Vec<[u32; 3]>);
29
30#[derive(Copy, Clone, Debug)]
31pub(crate) struct CutPlane {
32    pub abc: Vector<Real>,
33    pub d: Real,
34    pub axis: u8,
35    pub index: u32,
36}
37
38/// Approximate convex decomposition using the VHACD algorithm.
39pub struct VHACD {
40    // raycast_mesh: Option<RaycastMesh>,
41    voxel_parts: Vec<VoxelSet>,
42    volume_ch0: Real,
43    max_concavity: Real,
44}
45
46impl VHACD {
47    /// Decompose the given polyline (in 2D) or triangle mesh (in 3D).
48    ///
49    /// # Parameters
50    /// * `params` - The parameters for the VHACD algorithm execution.
51    /// * `points` - The vertex buffer of the polyline (in 2D) or triangle mesh (in 3D).
52    /// * `indices` - The index buffer of the polyline (in 2D) or triangle mesh (in 3D).
53    /// * `keep_voxel_to_primitives_map` - If set to `true` then a map between the voxels
54    ///   computed during the decomposition, and the primitives (triangle or segment) they
55    ///   intersect will be computed. This is required in order to compute the convex-hulls
56    ///   using the original polyline/trimesh primitives afterwards (otherwise the convex
57    ///   hulls resulting from the convex decomposition will use the voxels vertices).
58    pub fn decompose(
59        params: &VHACDParameters,
60        points: &[Point<Real>],
61        indices: &[[u32; DIM]],
62        keep_voxel_to_primitives_map: bool,
63    ) -> Self {
64        // if params.project_hull_vertices || params.fill_mode == FillMode::RAYCAST_FILL {
65        //     self.raycast_mesh =
66        //         RaycastMesh::create_raycast_mesh(num_points, points, num_triangles, triangles);
67        // }
68
69        let voxelized = VoxelizedVolume::voxelize(
70            points,
71            indices,
72            params.resolution,
73            params.fill_mode,
74            keep_voxel_to_primitives_map,
75            // &self.raycast_mesh,
76        );
77
78        let mut result = Self::from_voxels(params, voxelized.into());
79
80        let primitive_classes = Arc::new(result.classify_primitives(indices.len()));
81        for part in &mut result.voxel_parts {
82            part.primitive_classes = primitive_classes.clone();
83        }
84
85        result
86    }
87
88    /// Perform an approximate convex decomposition of a set of voxels.
89    pub fn from_voxels(params: &VHACDParameters, voxels: VoxelSet) -> Self {
90        let mut result = Self {
91            // raycast_mesh: None,
92            voxel_parts: Vec::new(),
93            volume_ch0: 0.0,
94            max_concavity: -Real::MAX,
95        };
96
97        result.do_compute_acd(params, voxels);
98        result
99    }
100
101    /// The almost-convex voxelized parts computed by the VHACD algorithm.
102    pub fn voxel_parts(&self) -> &[VoxelSet] {
103        &self.voxel_parts
104    }
105
106    #[cfg(feature = "dim2")]
107    fn compute_preferred_cutting_direction(eigenvalues: &Vector<Real>) -> (Vector<Real>, Real) {
108        let vx = eigenvalues.y * eigenvalues.y;
109        let vy = eigenvalues.x * eigenvalues.x;
110
111        if vx < vy {
112            let e = eigenvalues.y * eigenvalues.y;
113            let dir = Vector::x();
114
115            if e == 0.0 {
116                (dir, 0.0)
117            } else {
118                (dir, 1.0 - vx / e)
119            }
120        } else {
121            let e = eigenvalues.x * eigenvalues.x;
122            let dir = Vector::y();
123
124            if e == 0.0 {
125                (dir, 0.0)
126            } else {
127                (dir, 1.0 - vy / e)
128            }
129        }
130    }
131
132    #[cfg(feature = "dim3")]
133    fn compute_preferred_cutting_direction(eigenvalues: &Vector<Real>) -> (Vector<Real>, Real) {
134        let vx = (eigenvalues.y - eigenvalues.z) * (eigenvalues.y - eigenvalues.z);
135        let vy = (eigenvalues.x - eigenvalues.z) * (eigenvalues.x - eigenvalues.z);
136        let vz = (eigenvalues.x - eigenvalues.y) * (eigenvalues.x - eigenvalues.y);
137
138        if vx < vy && vx < vz {
139            let e = eigenvalues.y * eigenvalues.y + eigenvalues.z * eigenvalues.z;
140            let dir = Vector::x();
141
142            if e == 0.0 {
143                (dir, 0.0)
144            } else {
145                (dir, 1.0 - vx / e)
146            }
147        } else if vy < vx && vy < vz {
148            let e = eigenvalues.x * eigenvalues.x + eigenvalues.z * eigenvalues.z;
149            let dir = Vector::y();
150
151            if e == 0.0 {
152                (dir, 0.0)
153            } else {
154                (dir, 1.0 - vy / e)
155            }
156        } else {
157            let e = eigenvalues.x * eigenvalues.x + eigenvalues.y * eigenvalues.y;
158            let dir = Vector::z();
159
160            if e == 0.0 {
161                (dir, 0.0)
162            } else {
163                (dir, 1.0 - vz / e)
164            }
165        }
166    }
167
168    fn refine_axes_aligned_clipping_planes(
169        vset: &VoxelSet,
170        best_plane: &CutPlane,
171        downsampling: u32,
172        planes: &mut Vec<CutPlane>,
173    ) {
174        let min_v = vset.min_bb_voxels();
175        let max_v = vset.max_bb_voxels();
176
177        let best_id = best_plane.axis as usize;
178        let i0 = min_v[best_id].max(best_plane.index.saturating_sub(downsampling));
179        let i1 = max_v[best_id].min(best_plane.index + downsampling);
180
181        for i in i0..=i1 {
182            let plane = CutPlane {
183                abc: Vector::ith(best_id, 1.0),
184                axis: best_plane.axis,
185                d: -(vset.origin[best_id] + (i as Real + 0.5) * vset.scale),
186                index: i,
187            };
188            planes.push(plane);
189        }
190    }
191
192    // Returns the best plane, and the min concavity.
193    fn compute_best_clipping_plane(
194        &self,
195        input_voxels: &VoxelSet,
196        input_voxels_ch: &ConvexHull,
197        planes: &[CutPlane],
198        preferred_cutting_direction: &Vector<Real>,
199        w: Real,
200        alpha: Real,
201        beta: Real,
202        convex_hull_downsampling: u32,
203        params: &VHACDParameters,
204    ) -> (CutPlane, Real) {
205        let mut best_plane = planes[0];
206        let mut min_concavity = Real::MAX;
207        let mut i_best = -1;
208        let mut min_total = Real::MAX;
209
210        let mut left_ch;
211        let mut right_ch;
212        let mut left_ch_pts = Vec::new();
213        let mut right_ch_pts = Vec::new();
214        let mut left_voxels = VoxelSet::new();
215        let mut right_voxels = VoxelSet::new();
216        let mut on_surface_voxels = VoxelSet::new();
217
218        input_voxels.select_on_surface(&mut on_surface_voxels);
219
220        for (x, plane) in planes.iter().enumerate() {
221            // Compute convex hulls.
222            if params.convex_hull_approximation {
223                right_ch_pts.clear();
224                left_ch_pts.clear();
225
226                on_surface_voxels.intersect(
227                    plane,
228                    &mut right_ch_pts,
229                    &mut left_ch_pts,
230                    convex_hull_downsampling * 32,
231                );
232
233                clip_mesh(
234                    #[cfg(feature = "dim2")]
235                    input_voxels_ch,
236                    #[cfg(feature = "dim3")]
237                    &input_voxels_ch.0,
238                    plane,
239                    &mut right_ch_pts,
240                    &mut left_ch_pts,
241                );
242                right_ch = convex_hull(&right_ch_pts);
243                left_ch = convex_hull(&left_ch_pts);
244            } else {
245                on_surface_voxels.clip(plane, &mut right_voxels, &mut left_voxels);
246                right_ch = right_voxels.compute_convex_hull(convex_hull_downsampling);
247                left_ch = left_voxels.compute_convex_hull(convex_hull_downsampling);
248            }
249
250            let volume_left_ch = compute_volume(&left_ch);
251            let volume_right_ch = compute_volume(&right_ch);
252
253            // compute clipped volumes
254            let (volume_left, volume_right) = input_voxels.compute_clipped_volumes(plane);
255            let concavity_left = compute_concavity(volume_left, volume_left_ch, self.volume_ch0);
256            let concavity_right = compute_concavity(volume_right, volume_right_ch, self.volume_ch0);
257            let concavity = concavity_left + concavity_right;
258
259            // compute cost
260            let balance = alpha * (volume_left - volume_right).abs() / self.volume_ch0;
261            let d = w * plane.abc.dot(preferred_cutting_direction);
262            let symmetry = beta * d;
263            let total = concavity + balance + symmetry;
264
265            if total < min_total || (total == min_total && (x as i32) < i_best) {
266                min_concavity = concavity;
267                best_plane = *plane;
268                min_total = total;
269                i_best = x as i32;
270            }
271        }
272
273        (best_plane, min_concavity)
274    }
275
276    fn process_primitive_set(
277        &mut self,
278        params: &VHACDParameters,
279        first_iteration: bool,
280        parts: &mut Vec<VoxelSet>,
281        temp: &mut Vec<VoxelSet>,
282        mut voxels: VoxelSet,
283    ) {
284        let volume = voxels.compute_volume(); // Compute the volume for this primitive set
285        voxels.compute_bb(); // Compute the bounding box for this primitive set.
286        let voxels_convex_hull = voxels.compute_convex_hull(params.convex_hull_downsampling); // Generate the convex hull for this primitive set.
287
288        // Compute the volume of the convex hull
289        let volume_ch = compute_volume(&voxels_convex_hull);
290
291        // If this is the first iteration, store the volume of the base
292        if first_iteration {
293            self.volume_ch0 = volume_ch;
294        }
295
296        // Compute the concavity of this volume
297        let concavity = compute_concavity(volume, volume_ch, self.volume_ch0);
298
299        // Compute the volume error.
300        if concavity > params.concavity {
301            let eigenvalues = voxels.compute_principal_axes();
302            let (preferred_cutting_direction, w) =
303                Self::compute_preferred_cutting_direction(&eigenvalues);
304
305            let mut planes = Vec::new();
306            voxels.compute_axes_aligned_clipping_planes(params.plane_downsampling, &mut planes);
307
308            let (mut best_plane, mut min_concavity) = self.compute_best_clipping_plane(
309                &voxels,
310                &voxels_convex_hull,
311                &planes,
312                &preferred_cutting_direction,
313                w,
314                concavity * params.alpha,
315                concavity * params.beta,
316                params.convex_hull_downsampling,
317                params,
318            );
319
320            if params.plane_downsampling > 1 || params.convex_hull_downsampling > 1 {
321                let mut planes_ref = Vec::new();
322
323                Self::refine_axes_aligned_clipping_planes(
324                    &voxels,
325                    &best_plane,
326                    params.plane_downsampling,
327                    &mut planes_ref,
328                );
329
330                let best = self.compute_best_clipping_plane(
331                    &voxels,
332                    &voxels_convex_hull,
333                    &planes_ref,
334                    &preferred_cutting_direction,
335                    w,
336                    concavity * params.alpha,
337                    concavity * params.beta,
338                    1, // convex_hull_downsampling = 1
339                    params,
340                );
341
342                best_plane = best.0;
343                min_concavity = best.1;
344            }
345
346            if min_concavity > self.max_concavity {
347                self.max_concavity = min_concavity;
348            }
349
350            let mut best_left = VoxelSet::new();
351            let mut best_right = VoxelSet::new();
352
353            voxels.clip(&best_plane, &mut best_right, &mut best_left);
354
355            temp.push(best_left);
356            temp.push(best_right);
357        } else {
358            parts.push(voxels);
359        }
360    }
361
362    fn do_compute_acd(&mut self, params: &VHACDParameters, mut voxels: VoxelSet) {
363        let intersections = voxels.intersections.clone();
364        let mut input_parts = Vec::new();
365        let mut parts = Vec::new();
366        let mut temp = Vec::new();
367        input_parts.push(core::mem::take(&mut voxels));
368
369        let mut first_iteration = true;
370        self.volume_ch0 = 1.0;
371
372        // Compute the decomposition depth based on the number of convex hulls being requested.
373        let mut hull_count = 2;
374        let mut depth = 1;
375
376        while params.max_convex_hulls > hull_count {
377            depth += 1;
378            hull_count *= 2;
379        }
380
381        // We must always increment the decomposition depth one higher than the maximum number of hulls requested.
382        // The reason for this is as follows.
383        // Say, for example, the user requests 32 convex hulls exactly.  This would be a decomposition depth of 5.
384        // However, when we do that, we do *not* necessarily get 32 hulls as a result.  This is because, during
385        // the recursive descent of the binary tree, one or more of the leaf nodes may have no concavity and
386        // will not be split.  So, in this way, even with a decomposition depth of 5, you can produce fewer than
387        // 32 hulls.  So, in this case, we would set the decomposition depth to 6 (producing up to as high as 64 convex
388        // hulls). Then, the merge step which combines over-described hulls down to the user requested amount, we will end
389        // up getting exactly 32 convex hulls as a result. We could just allow the artist to directly control the
390        // decomposition depth directly, but this would be a bit too complex and the preference is simply to let them
391        // specify how many hulls they want and derive the solution from that.
392        depth += 1;
393
394        for _ in 0..depth {
395            if input_parts.is_empty() {
396                break;
397            }
398
399            for input_part in input_parts.drain(..) {
400                self.process_primitive_set(
401                    params,
402                    first_iteration,
403                    &mut parts,
404                    &mut temp,
405                    input_part,
406                );
407                first_iteration = false;
408            }
409
410            core::mem::swap(&mut input_parts, &mut temp);
411            // Note that temp is already clear because our previous for
412            // loop used `drain`. However we call `clear` here explicitly
413            // to make sure it still works if we remove the `drain` in the
414            // future.
415            temp.clear();
416        }
417
418        parts.append(&mut input_parts);
419        self.voxel_parts = parts;
420
421        for part in &mut self.voxel_parts {
422            part.intersections = intersections.clone();
423        }
424    }
425
426    // Returns a vector such that `result[i]` gives the index of the voxelized convex part that
427    // intersects it.
428    //
429    // If multiple convex parts intersect the same primitive, then `result[i]` is set to `u32::MAX`.
430    // This is used to avoid some useless triangle/segment cutting when computing the exact convex hull
431    // of a voxelized convex part.
432    fn classify_primitives(&self, num_primitives: usize) -> Vec<u32> {
433        if num_primitives == 0 {
434            return Vec::new();
435        }
436
437        const NO_CLASS: u32 = u32::MAX - 1;
438        const MULTICLASS: u32 = u32::MAX;
439
440        let mut primitive_classes = Vec::new();
441        primitive_classes.resize(num_primitives, NO_CLASS);
442
443        for (ipart, part) in self.voxel_parts.iter().enumerate() {
444            for voxel in &part.voxels {
445                let range = voxel.intersections_range.0..voxel.intersections_range.1;
446                for inter in &part.intersections[range] {
447                    let class = &mut primitive_classes[*inter as usize];
448                    if *class == NO_CLASS {
449                        *class = ipart as u32;
450                    } else if *class != ipart as u32 {
451                        *class = MULTICLASS;
452                    }
453                }
454            }
455        }
456
457        primitive_classes
458    }
459
460    /// Compute the intersections between the voxelized convex part of this decomposition,
461    /// and all the primitives from the original decomposed polyline/trimesh,
462    ///
463    /// This will panic if `keep_voxel_to_primitives_map` was set to `false` when initializing
464    /// `self`.
465    pub fn compute_primitive_intersections(
466        &self,
467        points: &[Point<Real>],
468        indices: &[[u32; DIM]],
469    ) -> Vec<Vec<Point<Real>>> {
470        self.voxel_parts
471            .iter()
472            .map(|part| part.compute_primitive_intersections(points, indices))
473            .collect()
474    }
475
476    /// Compute the convex-hulls of the parts computed by this approximate convex-decomposition,
477    /// taking into account the primitives from the original polyline/trimesh being decomposed.
478    ///
479    /// This will panic if `keep_voxel_to_primitives_map` was set to `false` when initializing
480    /// `self`.
481    #[cfg(feature = "dim2")]
482    pub fn compute_exact_convex_hulls(
483        &self,
484        points: &[Point<Real>],
485        indices: &[[u32; DIM]],
486    ) -> Vec<Vec<Point<Real>>> {
487        self.voxel_parts
488            .iter()
489            .map(|part| part.compute_exact_convex_hull(points, indices))
490            .collect()
491    }
492
493    /// Compute the convex-hulls of the parts computed by this approximate convex-decomposition,
494    /// taking into account the primitives from the original polyline/trimesh being decomposed.
495    ///
496    /// This will panic if `keep_voxel_to_primitives_map` was set to `false` when initializing
497    /// `self`.
498    #[cfg(feature = "dim3")]
499    pub fn compute_exact_convex_hulls(
500        &self,
501        points: &[Point<Real>],
502        indices: &[[u32; DIM]],
503    ) -> Vec<(Vec<Point<Real>>, Vec<[u32; DIM]>)> {
504        self.voxel_parts
505            .iter()
506            .map(|part| part.compute_exact_convex_hull(points, indices))
507            .collect()
508    }
509
510    /// Compute the convex hulls of the voxelized approximately-convex parts
511    /// computed by `self` on the voxelized model.
512    ///
513    /// Use `compute_exact_convex_hulls` instead if the original polyline/trimesh geometry
514    /// needs to be taken into account.
515    #[cfg(feature = "dim2")]
516    pub fn compute_convex_hulls(&self, downsampling: u32) -> Vec<Vec<Point<Real>>> {
517        let downsampling = downsampling.max(1);
518        self.voxel_parts
519            .iter()
520            .map(|part| part.compute_convex_hull(downsampling))
521            .collect()
522    }
523
524    /// Compute the convex hulls of the voxelized approximately-convex parts
525    /// computed by `self` on the voxelized model.
526    ///
527    /// Use `compute_exact_convex_hulls` instead if the original polyline/trimesh geometry
528    /// needs to be taken into account.
529    #[cfg(feature = "dim3")]
530    pub fn compute_convex_hulls(
531        &self,
532        downsampling: u32,
533    ) -> Vec<(Vec<Point<Real>>, Vec<[u32; DIM]>)> {
534        let downsampling = downsampling.max(1);
535        self.voxel_parts
536            .iter()
537            .map(|part| part.compute_convex_hull(downsampling))
538            .collect()
539    }
540}
541
542fn compute_concavity(volume: Real, volume_ch: Real, volume0: Real) -> Real {
543    (volume_ch - volume).abs() / volume0
544}
545
546fn clip_mesh(
547    points: &[Point<Real>],
548    plane: &CutPlane,
549    positive_part: &mut Vec<Point<Real>>,
550    negative_part: &mut Vec<Point<Real>>,
551) {
552    for pt in points {
553        let d = plane.abc.dot(&pt.coords) + plane.d;
554
555        if d > 0.0 {
556            positive_part.push(*pt);
557        } else if d < 0.0 {
558            negative_part.push(*pt);
559        } else {
560            positive_part.push(*pt);
561            negative_part.push(*pt);
562        }
563    }
564}
565
566#[cfg(feature = "dim2")]
567fn convex_hull(vertices: &[Point<Real>]) -> Vec<Point<Real>> {
568    if vertices.len() > 1 {
569        crate::transformation::convex_hull(vertices)
570    } else {
571        Vec::new()
572    }
573}
574
575#[cfg(feature = "dim3")]
576fn convex_hull(vertices: &[Point<Real>]) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
577    if vertices.len() > 2 {
578        crate::transformation::convex_hull(vertices)
579    } else {
580        (Vec::new(), Vec::new())
581    }
582}
583
584#[cfg(feature = "dim2")]
585fn compute_volume(polygon: &[Point<Real>]) -> Real {
586    if !polygon.is_empty() {
587        crate::mass_properties::details::convex_polygon_area_and_center_of_mass(polygon).0
588    } else {
589        0.0
590    }
591}
592
593#[cfg(feature = "dim3")]
594fn compute_volume(mesh: &(Vec<Point<Real>>, Vec<[u32; DIM]>)) -> Real {
595    if !mesh.0.is_empty() {
596        crate::mass_properties::details::trimesh_signed_volume_and_center_of_mass(&mesh.0, &mesh.1)
597            .0
598    } else {
599        0.0
600    }
601}