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.
39///
40/// This structure holds the result of the V-HACD (Volumetric Hierarchical Approximate
41/// Convex Decomposition) algorithm, which decomposes a concave shape into multiple
42/// approximately-convex parts.
43///
44/// # Overview
45///
46/// The `VHACD` struct stores the decomposition result as a collection of voxelized parts,
47/// where each part is approximately convex. These parts can be converted to convex hulls
48/// for use in collision detection and physics simulation.
49///
50/// # Basic Workflow
51///
52/// 1. **Decompose**: Create a `VHACD` instance using [`VHACD::decompose`] or [`VHACD::from_voxels`]
53/// 2. **Access Parts**: Get the voxelized parts using [`voxel_parts`](VHACD::voxel_parts)
54/// 3. **Generate Hulls**: Compute convex hulls with [`compute_convex_hulls`](VHACD::compute_convex_hulls)
55///    or [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls)
56///
57/// # Examples
58///
59/// ## Basic Usage
60///
61/// ```no_run
62/// # #[cfg(all(feature = "dim3", feature = "f32"))] {
63/// use parry3d::math::Point;
64/// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
65///
66/// // Define a simple mesh (tetrahedron)
67/// let vertices = vec![
68///     Point::new(0.0, 0.0, 0.0),
69///     Point::new(1.0, 0.0, 0.0),
70///     Point::new(0.5, 1.0, 0.0),
71///     Point::new(0.5, 0.5, 1.0),
72/// ];
73/// let indices = vec![
74///     [0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2],
75/// ];
76///
77/// // Decompose with default parameters
78/// let decomposition = VHACD::decompose(
79///     &VHACDParameters::default(),
80///     &vertices,
81///     &indices,
82///     false,
83/// );
84///
85/// // Access the results
86/// println!("Generated {} parts", decomposition.voxel_parts().len());
87///
88/// // Get convex hulls for collision detection
89/// let hulls = decomposition.compute_convex_hulls(4);
90/// # }
91/// ```
92///
93/// ## With Custom Parameters
94///
95/// ```no_run
96/// # #[cfg(all(feature = "dim3", feature = "f32"))] {
97/// use parry3d::math::Point;
98/// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
99///
100/// # let vertices = vec![
101/// #     Point::new(0.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0),
102/// #     Point::new(0.5, 1.0, 0.0), Point::new(0.5, 0.5, 1.0),
103/// # ];
104/// # let indices = vec![[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2]];
105/// #
106/// // High-quality decomposition settings
107/// let params = VHACDParameters {
108///     resolution: 128,
109///     concavity: 0.001,
110///     max_convex_hulls: 32,
111///     ..Default::default()
112/// };
113///
114/// let decomposition = VHACD::decompose(&params, &vertices, &indices, false);
115/// # }
116/// ```
117///
118/// # See Also
119///
120/// - [`VHACDParameters`]: Configuration for the decomposition algorithm
121/// - [`compute_convex_hulls`](VHACD::compute_convex_hulls): Generate convex hulls from voxels
122/// - [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls): Generate hulls from original mesh
123/// - Module documentation: [`crate::transformation::vhacd`]
124pub struct VHACD {
125    // raycast_mesh: Option<RaycastMesh>,
126    voxel_parts: Vec<VoxelSet>,
127    volume_ch0: Real,
128    max_concavity: Real,
129}
130
131impl VHACD {
132    /// Decompose the given polyline (in 2D) or triangle mesh (in 3D).
133    ///
134    /// This is the primary method for performing approximate convex decomposition. It takes
135    /// a mesh defined by vertices and indices, voxelizes it, and decomposes it into
136    /// approximately-convex parts using the V-HACD algorithm.
137    ///
138    /// # Parameters
139    ///
140    /// * `params` - Configuration parameters controlling the decomposition process.
141    ///   See [`VHACDParameters`] for details on each parameter.
142    ///
143    /// * `points` - The vertex positions of your mesh (3D) or polyline (2D).
144    ///   Each point represents a vertex in world space.
145    ///
146    /// * `indices` - The connectivity information:
147    ///   - **3D**: Triangle indices `[u32; 3]` - each entry defines a triangle using 3 vertex indices
148    ///   - **2D**: Segment indices `[u32; 2]` - each entry defines a line segment using 2 vertex indices
149    ///
150    /// * `keep_voxel_to_primitives_map` - Whether to maintain a mapping between voxels and
151    ///   the original mesh primitives (triangles/segments) they intersect.
152    ///   - **`true`**: Enables [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls)
153    ///     which uses original mesh geometry for more accurate results. Uses more memory.
154    ///   - **`false`**: Only voxel-based hulls available via [`compute_convex_hulls`](VHACD::compute_convex_hulls).
155    ///     More memory efficient.
156    ///
157    /// # Returns
158    ///
159    /// A `VHACD` instance containing the decomposition results. Use [`voxel_parts`](VHACD::voxel_parts)
160    /// to access the raw voxelized parts, or [`compute_convex_hulls`](VHACD::compute_convex_hulls)
161    /// to generate convex hull geometry.
162    ///
163    /// # Performance
164    ///
165    /// Decomposition time depends primarily on:
166    /// - **`resolution`**: Higher = slower (cubic scaling in 3D)
167    /// - **`max_convex_hulls`**: More parts = more splits = longer time
168    /// - **Mesh complexity**: More vertices/triangles = slower voxelization
169    ///
170    /// Typical times (on modern CPU):
171    /// - Simple mesh (1K triangles, resolution 64): 100-500ms
172    /// - Complex mesh (10K triangles, resolution 128): 1-5 seconds
173    /// - Very complex (100K triangles, resolution 256): 10-60 seconds
174    ///
175    /// # Examples
176    ///
177    /// ## Basic Decomposition
178    ///
179    /// ```no_run
180    /// # #[cfg(all(feature = "dim3", feature = "f32"))] {
181    /// use parry3d::math::Point;
182    /// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
183    ///
184    /// // Simple L-shaped mesh
185    /// let vertices = vec![
186    ///     Point::new(0.0, 0.0, 0.0), Point::new(2.0, 0.0, 0.0),
187    ///     Point::new(2.0, 1.0, 0.0), Point::new(1.0, 1.0, 0.0),
188    /// ];
189    /// let indices = vec![
190    ///     [0, 1, 2], [0, 2, 3],
191    /// ];
192    ///
193    /// let decomposition = VHACD::decompose(
194    ///     &VHACDParameters::default(),
195    ///     &vertices,
196    ///     &indices,
197    ///     false, // Don't need exact hulls
198    /// );
199    ///
200    /// println!("Parts: {}", decomposition.voxel_parts().len());
201    /// # }
202    /// ```
203    ///
204    /// ## With Exact Hull Generation
205    ///
206    /// ```no_run
207    /// # #[cfg(all(feature = "dim3", feature = "f32"))] {
208    /// use parry3d::math::Point;
209    /// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
210    ///
211    /// # let vertices = vec![
212    /// #     Point::new(0.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0),
213    /// #     Point::new(0.5, 1.0, 0.0), Point::new(0.5, 0.5, 1.0),
214    /// # ];
215    /// # let indices = vec![[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2]];
216    /// #
217    /// // Enable voxel-to-primitive mapping for exact hulls
218    /// let decomposition = VHACD::decompose(
219    ///     &VHACDParameters::default(),
220    ///     &vertices,
221    ///     &indices,
222    ///     true, // <-- Enable for exact hulls
223    /// );
224    ///
225    /// // Now we can compute exact hulls using original mesh
226    /// let exact_hulls = decomposition.compute_exact_convex_hulls(&vertices, &indices);
227    /// # }
228    /// ```
229    ///
230    /// ## Custom Parameters
231    ///
232    /// ```no_run
233    /// # #[cfg(all(feature = "dim3", feature = "f32"))] {
234    /// use parry3d::math::Point;
235    /// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
236    ///
237    /// # let vertices = vec![
238    /// #     Point::new(0.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0),
239    /// #     Point::new(0.5, 1.0, 0.0), Point::new(0.5, 0.5, 1.0),
240    /// # ];
241    /// # let indices = vec![[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2]];
242    /// #
243    /// // High-quality settings for important objects
244    /// let params = VHACDParameters {
245    ///     resolution: 128,        // High detail
246    ///     concavity: 0.001,       // Tight fit
247    ///     max_convex_hulls: 32,   // Allow many parts
248    ///     ..Default::default()
249    /// };
250    ///
251    /// let decomposition = VHACD::decompose(&params, &vertices, &indices, false);
252    /// # }
253    /// ```
254    ///
255    /// # See Also
256    ///
257    /// - [`VHACDParameters`]: Detailed parameter documentation
258    /// - [`from_voxels`](VHACD::from_voxels): Decompose pre-voxelized data
259    /// - [`compute_convex_hulls`](VHACD::compute_convex_hulls): Generate voxel-based hulls
260    /// - [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls): Generate mesh-based hulls
261    pub fn decompose(
262        params: &VHACDParameters,
263        points: &[Point<Real>],
264        indices: &[[u32; DIM]],
265        keep_voxel_to_primitives_map: bool,
266    ) -> Self {
267        // if params.project_hull_vertices || params.fill_mode == FillMode::RAYCAST_FILL {
268        //     self.raycast_mesh =
269        //         RaycastMesh::create_raycast_mesh(num_points, points, num_triangles, triangles);
270        // }
271
272        let voxelized = VoxelizedVolume::voxelize(
273            points,
274            indices,
275            params.resolution,
276            params.fill_mode,
277            keep_voxel_to_primitives_map,
278            // &self.raycast_mesh,
279        );
280
281        let mut result = Self::from_voxels(params, voxelized.into());
282
283        let primitive_classes = Arc::new(result.classify_primitives(indices.len()));
284        for part in &mut result.voxel_parts {
285            part.primitive_classes = primitive_classes.clone();
286        }
287
288        result
289    }
290
291    /// Perform an approximate convex decomposition of a pre-voxelized set of voxels.
292    ///
293    /// This method allows you to decompose a shape that has already been voxelized,
294    /// bypassing the voxelization step. This is useful if you:
295    /// - Already have voxelized data from another source
296    /// - Want to decompose the same voxelization with different parameters
297    /// - Need more control over the voxelization process
298    ///
299    /// # Parameters
300    ///
301    /// * `params` - Configuration parameters for the decomposition algorithm.
302    ///   See [`VHACDParameters`] for details.
303    ///
304    /// * `voxels` - A pre-voxelized volume represented as a [`VoxelSet`].
305    ///   You can create this using [`VoxelizedVolume::voxelize`] or other voxelization methods.
306    ///
307    /// # Returns
308    ///
309    /// A `VHACD` instance containing the decomposition results.
310    ///
311    /// # Examples
312    ///
313    /// ```no_run
314    /// # #[cfg(all(feature = "dim3", feature = "f32"))] {
315    /// use parry3d::math::Point;
316    /// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
317    /// use parry3d::transformation::voxelization::{VoxelizedVolume, FillMode};
318    ///
319    /// # let vertices = vec![
320    /// #     Point::new(0.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0),
321    /// #     Point::new(0.5, 1.0, 0.0), Point::new(0.5, 0.5, 1.0),
322    /// # ];
323    /// # let indices = vec![[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2]];
324    /// #
325    /// // First, voxelize the mesh manually
326    /// let voxelized = VoxelizedVolume::voxelize(
327    ///     &vertices,
328    ///     &indices,
329    ///     64, // resolution
330    ///     FillMode::FloodFill {
331    ///         detect_cavities: false,
332    ///     },
333    ///     false, // don't keep primitive mapping
334    /// );
335    ///
336    /// // Then decompose the voxels
337    /// let decomposition = VHACD::from_voxels(
338    ///     &VHACDParameters::default(),
339    ///     voxelized.into(),
340    /// );
341    /// # }
342    /// ```
343    ///
344    /// # See Also
345    ///
346    /// - [`decompose`](VHACD::decompose): Decompose directly from mesh (includes voxelization)
347    /// - [`VoxelizedVolume`]: For manual voxelization
348    /// - [`VoxelSet`]: The voxel data structure
349    ///
350    /// [`VoxelizedVolume::voxelize`]: crate::transformation::voxelization::VoxelizedVolume::voxelize
351    /// [`VoxelSet`]: crate::transformation::voxelization::VoxelSet
352    /// [`VoxelizedVolume`]: crate::transformation::voxelization::VoxelizedVolume
353    pub fn from_voxels(params: &VHACDParameters, voxels: VoxelSet) -> Self {
354        let mut result = Self {
355            // raycast_mesh: None,
356            voxel_parts: Vec::new(),
357            volume_ch0: 0.0,
358            max_concavity: -Real::MAX,
359        };
360
361        result.do_compute_acd(params, voxels);
362        result
363    }
364
365    /// Returns the approximately-convex voxelized parts computed by the VHACD algorithm.
366    ///
367    /// Each part in the returned slice represents an approximately-convex region of the
368    /// original shape, stored as a set of voxels. These voxelized parts are the direct
369    /// result of the decomposition algorithm.
370    ///
371    /// # Returns
372    ///
373    /// A slice of [`VoxelSet`] structures, where each set represents one convex part.
374    /// The number of parts depends on the shape's complexity and the parameters used
375    /// (especially `concavity` and `max_convex_hulls`).
376    ///
377    /// # Usage
378    ///
379    /// You typically don't use the voxel parts directly for collision detection. Instead:
380    /// - Use [`compute_convex_hulls`](VHACD::compute_convex_hulls) for voxel-based convex hulls
381    /// - Use [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls) for mesh-based hulls
382    ///
383    /// However, accessing voxel parts directly is useful for:
384    /// - Debugging and visualization of the decomposition
385    /// - Custom processing of the voxelized representation
386    /// - Understanding how the algorithm divided the shape
387    ///
388    /// # Examples
389    ///
390    /// ```no_run
391    /// # #[cfg(all(feature = "dim3", feature = "f32"))] {
392    /// use parry3d::math::Point;
393    /// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
394    ///
395    /// # let vertices = vec![
396    /// #     Point::new(0.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0),
397    /// #     Point::new(0.5, 1.0, 0.0), Point::new(0.5, 0.5, 1.0),
398    /// # ];
399    /// # let indices = vec![[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2]];
400    /// #
401    /// let decomposition = VHACD::decompose(
402    ///     &VHACDParameters::default(),
403    ///     &vertices,
404    ///     &indices,
405    ///     false,
406    /// );
407    ///
408    /// let parts = decomposition.voxel_parts();
409    /// println!("Generated {} convex parts", parts.len());
410    ///
411    /// // Inspect individual parts
412    /// for (i, part) in parts.iter().enumerate() {
413    ///     println!("Part {}: {} voxels", i, part.voxels().len());
414    /// }
415    /// # }
416    /// ```
417    ///
418    /// # See Also
419    ///
420    /// - [`VoxelSet`]: The voxel data structure
421    /// - [`compute_convex_hulls`](VHACD::compute_convex_hulls): Convert to collision-ready convex hulls
422    ///
423    /// [`VoxelSet`]: crate::transformation::voxelization::VoxelSet
424    pub fn voxel_parts(&self) -> &[VoxelSet] {
425        &self.voxel_parts
426    }
427
428    #[cfg(feature = "dim2")]
429    fn compute_preferred_cutting_direction(eigenvalues: &Vector<Real>) -> (Vector<Real>, Real) {
430        let vx = eigenvalues.y * eigenvalues.y;
431        let vy = eigenvalues.x * eigenvalues.x;
432
433        if vx < vy {
434            let e = eigenvalues.y * eigenvalues.y;
435            let dir = Vector::x();
436
437            if e == 0.0 {
438                (dir, 0.0)
439            } else {
440                (dir, 1.0 - vx / e)
441            }
442        } else {
443            let e = eigenvalues.x * eigenvalues.x;
444            let dir = Vector::y();
445
446            if e == 0.0 {
447                (dir, 0.0)
448            } else {
449                (dir, 1.0 - vy / e)
450            }
451        }
452    }
453
454    #[cfg(feature = "dim3")]
455    fn compute_preferred_cutting_direction(eigenvalues: &Vector<Real>) -> (Vector<Real>, Real) {
456        let vx = (eigenvalues.y - eigenvalues.z) * (eigenvalues.y - eigenvalues.z);
457        let vy = (eigenvalues.x - eigenvalues.z) * (eigenvalues.x - eigenvalues.z);
458        let vz = (eigenvalues.x - eigenvalues.y) * (eigenvalues.x - eigenvalues.y);
459
460        if vx < vy && vx < vz {
461            let e = eigenvalues.y * eigenvalues.y + eigenvalues.z * eigenvalues.z;
462            let dir = Vector::x();
463
464            if e == 0.0 {
465                (dir, 0.0)
466            } else {
467                (dir, 1.0 - vx / e)
468            }
469        } else if vy < vx && vy < vz {
470            let e = eigenvalues.x * eigenvalues.x + eigenvalues.z * eigenvalues.z;
471            let dir = Vector::y();
472
473            if e == 0.0 {
474                (dir, 0.0)
475            } else {
476                (dir, 1.0 - vy / e)
477            }
478        } else {
479            let e = eigenvalues.x * eigenvalues.x + eigenvalues.y * eigenvalues.y;
480            let dir = Vector::z();
481
482            if e == 0.0 {
483                (dir, 0.0)
484            } else {
485                (dir, 1.0 - vz / e)
486            }
487        }
488    }
489
490    fn refine_axes_aligned_clipping_planes(
491        vset: &VoxelSet,
492        best_plane: &CutPlane,
493        downsampling: u32,
494        planes: &mut Vec<CutPlane>,
495    ) {
496        let min_v = vset.min_bb_voxels();
497        let max_v = vset.max_bb_voxels();
498
499        let best_id = best_plane.axis as usize;
500        let i0 = min_v[best_id].max(best_plane.index.saturating_sub(downsampling));
501        let i1 = max_v[best_id].min(best_plane.index + downsampling);
502
503        for i in i0..=i1 {
504            let plane = CutPlane {
505                abc: Vector::ith(best_id, 1.0),
506                axis: best_plane.axis,
507                d: -(vset.origin[best_id] + (i as Real + 0.5) * vset.scale),
508                index: i,
509            };
510            planes.push(plane);
511        }
512    }
513
514    // Returns the best plane, and the min concavity.
515    fn compute_best_clipping_plane(
516        &self,
517        input_voxels: &VoxelSet,
518        input_voxels_ch: &ConvexHull,
519        planes: &[CutPlane],
520        preferred_cutting_direction: &Vector<Real>,
521        w: Real,
522        alpha: Real,
523        beta: Real,
524        convex_hull_downsampling: u32,
525        params: &VHACDParameters,
526    ) -> (CutPlane, Real) {
527        let mut best_plane = planes[0];
528        let mut min_concavity = Real::MAX;
529        let mut i_best = -1;
530        let mut min_total = Real::MAX;
531
532        let mut left_ch;
533        let mut right_ch;
534        let mut left_ch_pts = Vec::new();
535        let mut right_ch_pts = Vec::new();
536        let mut left_voxels = VoxelSet::new();
537        let mut right_voxels = VoxelSet::new();
538        let mut on_surface_voxels = VoxelSet::new();
539
540        input_voxels.select_on_surface(&mut on_surface_voxels);
541
542        for (x, plane) in planes.iter().enumerate() {
543            // Compute convex hulls.
544            if params.convex_hull_approximation {
545                right_ch_pts.clear();
546                left_ch_pts.clear();
547
548                on_surface_voxels.intersect(
549                    plane,
550                    &mut right_ch_pts,
551                    &mut left_ch_pts,
552                    convex_hull_downsampling * 32,
553                );
554
555                clip_mesh(
556                    #[cfg(feature = "dim2")]
557                    input_voxels_ch,
558                    #[cfg(feature = "dim3")]
559                    &input_voxels_ch.0,
560                    plane,
561                    &mut right_ch_pts,
562                    &mut left_ch_pts,
563                );
564                right_ch = convex_hull(&right_ch_pts);
565                left_ch = convex_hull(&left_ch_pts);
566            } else {
567                on_surface_voxels.clip(plane, &mut right_voxels, &mut left_voxels);
568                right_ch = right_voxels.compute_convex_hull(convex_hull_downsampling);
569                left_ch = left_voxels.compute_convex_hull(convex_hull_downsampling);
570            }
571
572            let volume_left_ch = compute_volume(&left_ch);
573            let volume_right_ch = compute_volume(&right_ch);
574
575            // compute clipped volumes
576            let (volume_left, volume_right) = input_voxels.compute_clipped_volumes(plane);
577            let concavity_left = compute_concavity(volume_left, volume_left_ch, self.volume_ch0);
578            let concavity_right = compute_concavity(volume_right, volume_right_ch, self.volume_ch0);
579            let concavity = concavity_left + concavity_right;
580
581            // compute cost
582            let balance = alpha * (volume_left - volume_right).abs() / self.volume_ch0;
583            let d = w * plane.abc.dot(preferred_cutting_direction);
584            let symmetry = beta * d;
585            let total = concavity + balance + symmetry;
586
587            if total < min_total || (total == min_total && (x as i32) < i_best) {
588                min_concavity = concavity;
589                best_plane = *plane;
590                min_total = total;
591                i_best = x as i32;
592            }
593        }
594
595        (best_plane, min_concavity)
596    }
597
598    fn process_primitive_set(
599        &mut self,
600        params: &VHACDParameters,
601        first_iteration: bool,
602        parts: &mut Vec<VoxelSet>,
603        temp: &mut Vec<VoxelSet>,
604        mut voxels: VoxelSet,
605    ) {
606        let volume = voxels.compute_volume(); // Compute the volume for this primitive set
607        voxels.compute_bb(); // Compute the bounding box for this primitive set.
608        let voxels_convex_hull = voxels.compute_convex_hull(params.convex_hull_downsampling); // Generate the convex hull for this primitive set.
609
610        // Compute the volume of the convex hull
611        let volume_ch = compute_volume(&voxels_convex_hull);
612
613        // If this is the first iteration, store the volume of the base
614        if first_iteration {
615            self.volume_ch0 = volume_ch;
616        }
617
618        // Compute the concavity of this volume
619        let concavity = compute_concavity(volume, volume_ch, self.volume_ch0);
620
621        // Compute the volume error.
622        if concavity > params.concavity {
623            let eigenvalues = voxels.compute_principal_axes();
624            let (preferred_cutting_direction, w) =
625                Self::compute_preferred_cutting_direction(&eigenvalues);
626
627            let mut planes = Vec::new();
628            voxels.compute_axes_aligned_clipping_planes(params.plane_downsampling, &mut planes);
629
630            let (mut best_plane, mut min_concavity) = self.compute_best_clipping_plane(
631                &voxels,
632                &voxels_convex_hull,
633                &planes,
634                &preferred_cutting_direction,
635                w,
636                concavity * params.alpha,
637                concavity * params.beta,
638                params.convex_hull_downsampling,
639                params,
640            );
641
642            if params.plane_downsampling > 1 || params.convex_hull_downsampling > 1 {
643                let mut planes_ref = Vec::new();
644
645                Self::refine_axes_aligned_clipping_planes(
646                    &voxels,
647                    &best_plane,
648                    params.plane_downsampling,
649                    &mut planes_ref,
650                );
651
652                let best = self.compute_best_clipping_plane(
653                    &voxels,
654                    &voxels_convex_hull,
655                    &planes_ref,
656                    &preferred_cutting_direction,
657                    w,
658                    concavity * params.alpha,
659                    concavity * params.beta,
660                    1, // convex_hull_downsampling = 1
661                    params,
662                );
663
664                best_plane = best.0;
665                min_concavity = best.1;
666            }
667
668            if min_concavity > self.max_concavity {
669                self.max_concavity = min_concavity;
670            }
671
672            let mut best_left = VoxelSet::new();
673            let mut best_right = VoxelSet::new();
674
675            voxels.clip(&best_plane, &mut best_right, &mut best_left);
676
677            temp.push(best_left);
678            temp.push(best_right);
679        } else {
680            parts.push(voxels);
681        }
682    }
683
684    fn do_compute_acd(&mut self, params: &VHACDParameters, mut voxels: VoxelSet) {
685        let intersections = voxels.intersections.clone();
686        let mut input_parts = Vec::new();
687        let mut parts = Vec::new();
688        let mut temp = Vec::new();
689        input_parts.push(core::mem::take(&mut voxels));
690
691        let mut first_iteration = true;
692        self.volume_ch0 = 1.0;
693
694        // Compute the decomposition depth based on the number of convex hulls being requested.
695        let mut hull_count = 2;
696        let mut depth = 1;
697
698        while params.max_convex_hulls > hull_count {
699            depth += 1;
700            hull_count *= 2;
701        }
702
703        // We must always increment the decomposition depth one higher than the maximum number of hulls requested.
704        // The reason for this is as follows.
705        // Say, for example, the user requests 32 convex hulls exactly.  This would be a decomposition depth of 5.
706        // However, when we do that, we do *not* necessarily get 32 hulls as a result.  This is because, during
707        // the recursive descent of the binary tree, one or more of the leaf nodes may have no concavity and
708        // will not be split.  So, in this way, even with a decomposition depth of 5, you can produce fewer than
709        // 32 hulls.  So, in this case, we would set the decomposition depth to 6 (producing up to as high as 64 convex
710        // hulls). Then, the merge step which combines over-described hulls down to the user requested amount, we will end
711        // up getting exactly 32 convex hulls as a result. We could just allow the artist to directly control the
712        // decomposition depth directly, but this would be a bit too complex and the preference is simply to let them
713        // specify how many hulls they want and derive the solution from that.
714        depth += 1;
715
716        for _ in 0..depth {
717            if input_parts.is_empty() {
718                break;
719            }
720
721            for input_part in input_parts.drain(..) {
722                self.process_primitive_set(
723                    params,
724                    first_iteration,
725                    &mut parts,
726                    &mut temp,
727                    input_part,
728                );
729                first_iteration = false;
730            }
731
732            core::mem::swap(&mut input_parts, &mut temp);
733            // Note that temp is already clear because our previous for
734            // loop used `drain`. However we call `clear` here explicitly
735            // to make sure it still works if we remove the `drain` in the
736            // future.
737            temp.clear();
738        }
739
740        parts.append(&mut input_parts);
741        self.voxel_parts = parts;
742
743        for part in &mut self.voxel_parts {
744            part.intersections = intersections.clone();
745        }
746    }
747
748    // Returns a vector such that `result[i]` gives the index of the voxelized convex part that
749    // intersects it.
750    //
751    // If multiple convex parts intersect the same primitive, then `result[i]` is set to `u32::MAX`.
752    // This is used to avoid some useless triangle/segment cutting when computing the exact convex hull
753    // of a voxelized convex part.
754    fn classify_primitives(&self, num_primitives: usize) -> Vec<u32> {
755        if num_primitives == 0 {
756            return Vec::new();
757        }
758
759        const NO_CLASS: u32 = u32::MAX - 1;
760        const MULTICLASS: u32 = u32::MAX;
761
762        let mut primitive_classes = Vec::new();
763        primitive_classes.resize(num_primitives, NO_CLASS);
764
765        for (ipart, part) in self.voxel_parts.iter().enumerate() {
766            for voxel in &part.voxels {
767                let range = voxel.intersections_range.0..voxel.intersections_range.1;
768                for inter in &part.intersections[range] {
769                    let class = &mut primitive_classes[*inter as usize];
770                    if *class == NO_CLASS {
771                        *class = ipart as u32;
772                    } else if *class != ipart as u32 {
773                        *class = MULTICLASS;
774                    }
775                }
776            }
777        }
778
779        primitive_classes
780    }
781
782    /// Compute the intersections between the voxelized convex part of this decomposition,
783    /// and all the primitives from the original decomposed polyline/trimesh,
784    ///
785    /// This will panic if `keep_voxel_to_primitives_map` was set to `false` when initializing
786    /// `self`.
787    pub fn compute_primitive_intersections(
788        &self,
789        points: &[Point<Real>],
790        indices: &[[u32; DIM]],
791    ) -> Vec<Vec<Point<Real>>> {
792        self.voxel_parts
793            .iter()
794            .map(|part| part.compute_primitive_intersections(points, indices))
795            .collect()
796    }
797
798    /// Compute exact convex hulls using the original mesh geometry (2D version).
799    ///
800    /// This method generates convex hulls for each decomposed part by computing the convex
801    /// hull of the intersection between the voxelized part and the **original polyline primitives**.
802    /// This produces more accurate hulls than the voxel-based method, preserving the original
803    /// geometry's detail.
804    ///
805    /// # Requirements
806    ///
807    /// This method requires that `keep_voxel_to_primitives_map` was set to `true` when calling
808    /// [`VHACD::decompose`]. If it was `false`, this method will **panic**.
809    ///
810    /// # Parameters
811    ///
812    /// * `points` - The same vertex buffer used when calling [`VHACD::decompose`]
813    /// * `indices` - The same index buffer (segment indices `[u32; 2]`) used when calling
814    ///   [`VHACD::decompose`]
815    ///
816    /// # Returns
817    ///
818    /// A vector of convex polygons (one per decomposed part), where each polygon is
819    /// represented as a `Vec<Point<Real>>` containing the hull vertices in order.
820    ///
821    /// # Performance
822    ///
823    /// This is more expensive than [`compute_convex_hulls`](VHACD::compute_convex_hulls) because
824    /// it needs to compute intersections with the original primitives and then compute convex hulls.
825    /// However, it produces more accurate results.
826    ///
827    /// # Examples
828    ///
829    /// ```no_run
830    /// # #[cfg(all(feature = "dim2", feature = "f32"))] {
831    /// use parry2d::math::Point;
832    /// use parry2d::transformation::vhacd::{VHACD, VHACDParameters};
833    ///
834    /// // Define an L-shaped polyline
835    /// let vertices = vec![
836    ///     Point::new(0.0, 0.0), Point::new(2.0, 0.0),
837    ///     Point::new(2.0, 1.0), Point::new(1.0, 1.0),
838    ///     Point::new(1.0, 2.0), Point::new(0.0, 2.0),
839    /// ];
840    /// let indices = vec![
841    ///     [0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0],
842    /// ];
843    ///
844    /// // IMPORTANT: Set keep_voxel_to_primitives_map to true
845    /// let decomposition = VHACD::decompose(
846    ///     &VHACDParameters::default(),
847    ///     &vertices,
848    ///     &indices,
849    ///     true, // <-- Required for exact hulls
850    /// );
851    ///
852    /// // Compute exact convex hulls using original geometry
853    /// let exact_hulls = decomposition.compute_exact_convex_hulls(&vertices, &indices);
854    ///
855    /// for (i, hull) in exact_hulls.iter().enumerate() {
856    ///     println!("Hull {}: {} vertices", i, hull.len());
857    /// }
858    /// # }
859    /// ```
860    ///
861    /// # Panics
862    ///
863    /// Panics if `keep_voxel_to_primitives_map` was `false` during decomposition.
864    ///
865    /// # See Also
866    ///
867    /// - [`compute_convex_hulls`](VHACD::compute_convex_hulls): Faster voxel-based hulls
868    /// - [`compute_primitive_intersections`](VHACD::compute_primitive_intersections): Get raw intersection points
869    #[cfg(feature = "dim2")]
870    pub fn compute_exact_convex_hulls(
871        &self,
872        points: &[Point<Real>],
873        indices: &[[u32; DIM]],
874    ) -> Vec<Vec<Point<Real>>> {
875        self.voxel_parts
876            .iter()
877            .map(|part| part.compute_exact_convex_hull(points, indices))
878            .collect()
879    }
880
881    /// Compute exact convex hulls using the original mesh geometry (3D version).
882    ///
883    /// This method generates convex hulls for each decomposed part by computing the convex
884    /// hull of the intersection between the voxelized part and the **original triangle mesh
885    /// primitives**. This produces more accurate hulls than the voxel-based method, preserving
886    /// the original geometry's detail.
887    ///
888    /// # Requirements
889    ///
890    /// This method requires that `keep_voxel_to_primitives_map` was set to `true` when calling
891    /// [`VHACD::decompose`]. If it was `false`, this method will **panic**.
892    ///
893    /// # Parameters
894    ///
895    /// * `points` - The same vertex buffer used when calling [`VHACD::decompose`]
896    /// * `indices` - The same index buffer (triangle indices `[u32; 3]`) used when calling
897    ///   [`VHACD::decompose`]
898    ///
899    /// # Returns
900    ///
901    /// A vector of convex hulls (one per decomposed part), where each hull is represented as
902    /// a tuple `(vertices, indices)`:
903    /// - `vertices`: `Vec<Point<Real>>` - The hull vertices
904    /// - `indices`: `Vec<[u32; 3]>` - Triangle indices defining the hull surface
905    ///
906    /// # Performance
907    ///
908    /// This is more expensive than [`compute_convex_hulls`](VHACD::compute_convex_hulls) because
909    /// it needs to compute intersections with the original primitives and then compute convex hulls.
910    /// However, it produces more accurate results that better match the original mesh.
911    ///
912    /// # Examples
913    ///
914    /// ```no_run
915    /// # #[cfg(all(feature = "dim3", feature = "f32"))] {
916    /// use parry3d::math::Point;
917    /// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
918    ///
919    /// # let vertices = vec![
920    /// #     Point::new(0.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0),
921    /// #     Point::new(0.5, 1.0, 0.0), Point::new(0.5, 0.5, 1.0),
922    /// # ];
923    /// # let indices = vec![[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2]];
924    /// #
925    /// // IMPORTANT: Set keep_voxel_to_primitives_map to true
926    /// let decomposition = VHACD::decompose(
927    ///     &VHACDParameters::default(),
928    ///     &vertices,
929    ///     &indices,
930    ///     true, // <-- Required for exact hulls
931    /// );
932    ///
933    /// // Compute exact convex hulls using original geometry
934    /// let exact_hulls = decomposition.compute_exact_convex_hulls(&vertices, &indices);
935    ///
936    /// for (i, (verts, tris)) in exact_hulls.iter().enumerate() {
937    ///     println!("Hull {}: {} vertices, {} triangles", i, verts.len(), tris.len());
938    /// }
939    /// # }
940    /// ```
941    ///
942    /// # Panics
943    ///
944    /// Panics if `keep_voxel_to_primitives_map` was `false` during decomposition.
945    ///
946    /// # See Also
947    ///
948    /// - [`compute_convex_hulls`](VHACD::compute_convex_hulls): Faster voxel-based hulls
949    /// - [`compute_primitive_intersections`](VHACD::compute_primitive_intersections): Get raw intersection points
950    #[cfg(feature = "dim3")]
951    pub fn compute_exact_convex_hulls(
952        &self,
953        points: &[Point<Real>],
954        indices: &[[u32; DIM]],
955    ) -> Vec<(Vec<Point<Real>>, Vec<[u32; DIM]>)> {
956        self.voxel_parts
957            .iter()
958            .map(|part| part.compute_exact_convex_hull(points, indices))
959            .collect()
960    }
961
962    /// Compute convex hulls from the voxelized parts (2D version).
963    ///
964    /// This method generates convex polygons for each decomposed part by computing the convex
965    /// hull of the **voxel vertices**. This is faster than [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls)
966    /// but the hulls are based on the voxelized representation rather than the original geometry.
967    ///
968    /// # Parameters
969    ///
970    /// * `downsampling` - Controls how many voxels to skip when generating the hull.
971    ///   Higher values = fewer points = simpler hulls = faster computation.
972    ///   - `1`: Use all voxel vertices (highest quality, slowest)
973    ///   - `4`: Use every 4th voxel (good balance, recommended)
974    ///   - `8+`: Use fewer voxels (fastest, simpler hulls)
975    ///
976    ///   Values less than 1 are clamped to 1.
977    ///
978    /// # Returns
979    ///
980    /// A vector of convex polygons (one per decomposed part), where each polygon is
981    /// represented as a `Vec<Point<Real>>` containing the hull vertices in order.
982    ///
983    /// # Performance
984    ///
985    /// This method is faster than [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls)
986    /// because it works directly with voxel data without needing to intersect with original
987    /// mesh primitives. However, the resulting hulls are slightly less accurate.
988    ///
989    /// # When to Use
990    ///
991    /// Use this method when:
992    /// - You don't need the highest accuracy
993    /// - Performance is important
994    /// - You didn't enable `keep_voxel_to_primitives_map` during decomposition
995    /// - The voxel resolution is high enough for your needs
996    ///
997    /// Use [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls) when:
998    /// - You need the most accurate representation of the original geometry
999    /// - You have enabled `keep_voxel_to_primitives_map`
1000    ///
1001    /// # Examples
1002    ///
1003    /// ```no_run
1004    /// # #[cfg(all(feature = "dim2", feature = "f32"))] {
1005    /// use parry2d::math::Point;
1006    /// use parry2d::transformation::vhacd::{VHACD, VHACDParameters};
1007    ///
1008    /// // Define an L-shaped polyline
1009    /// let vertices = vec![
1010    ///     Point::new(0.0, 0.0), Point::new(2.0, 0.0),
1011    ///     Point::new(2.0, 1.0), Point::new(1.0, 1.0),
1012    ///     Point::new(1.0, 2.0), Point::new(0.0, 2.0),
1013    /// ];
1014    /// let indices = vec![
1015    ///     [0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0],
1016    /// ];
1017    ///
1018    /// let decomposition = VHACD::decompose(
1019    ///     &VHACDParameters::default(),
1020    ///     &vertices,
1021    ///     &indices,
1022    ///     false, // voxel-to-primitive mapping not needed
1023    /// );
1024    ///
1025    /// // Compute voxel-based convex hulls with moderate downsampling
1026    /// let hulls = decomposition.compute_convex_hulls(4);
1027    ///
1028    /// for (i, hull) in hulls.iter().enumerate() {
1029    ///     println!("Hull {}: {} vertices", i, hull.len());
1030    /// }
1031    /// # }
1032    /// ```
1033    ///
1034    /// # See Also
1035    ///
1036    /// - [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls): More accurate mesh-based hulls
1037    /// - [`voxel_parts`](VHACD::voxel_parts): Access the raw voxel data
1038    #[cfg(feature = "dim2")]
1039    pub fn compute_convex_hulls(&self, downsampling: u32) -> Vec<Vec<Point<Real>>> {
1040        let downsampling = downsampling.max(1);
1041        self.voxel_parts
1042            .iter()
1043            .map(|part| part.compute_convex_hull(downsampling))
1044            .collect()
1045    }
1046
1047    /// Compute convex hulls from the voxelized parts (3D version).
1048    ///
1049    /// This method generates convex meshes for each decomposed part by computing the convex
1050    /// hull of the **voxel vertices**. This is faster than [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls)
1051    /// but the hulls are based on the voxelized representation rather than the original geometry.
1052    ///
1053    /// # Parameters
1054    ///
1055    /// * `downsampling` - Controls how many voxels to skip when generating the hull.
1056    ///   Higher values = fewer points = simpler hulls = faster computation.
1057    ///   - `1`: Use all voxel vertices (highest quality, slowest)
1058    ///   - `4`: Use every 4th voxel (good balance, recommended)
1059    ///   - `8+`: Use fewer voxels (fastest, simpler hulls)
1060    ///
1061    ///   Values less than 1 are clamped to 1.
1062    ///
1063    /// # Returns
1064    ///
1065    /// A vector of convex hulls (one per decomposed part), where each hull is represented as
1066    /// a tuple `(vertices, indices)`:
1067    /// - `vertices`: `Vec<Point<Real>>` - The hull vertices
1068    /// - `indices`: `Vec<[u32; 3]>` - Triangle indices defining the hull surface
1069    ///
1070    /// # Performance
1071    ///
1072    /// This method is faster than [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls)
1073    /// because it works directly with voxel data without needing to intersect with original
1074    /// mesh primitives. The performance scales with the number of voxels and the downsampling factor.
1075    ///
1076    /// # When to Use
1077    ///
1078    /// Use this method when:
1079    /// - You don't need the highest accuracy
1080    /// - Performance is important
1081    /// - You didn't enable `keep_voxel_to_primitives_map` during decomposition
1082    /// - The voxel resolution is high enough for your needs
1083    /// - You're using this for real-time collision detection
1084    ///
1085    /// Use [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls) when:
1086    /// - You need the most accurate representation of the original geometry
1087    /// - You have enabled `keep_voxel_to_primitives_map`
1088    /// - Quality is more important than speed
1089    ///
1090    /// # Examples
1091    ///
1092    /// ```no_run
1093    /// # #[cfg(all(feature = "dim3", feature = "f32"))] {
1094    /// use parry3d::math::Point;
1095    /// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
1096    ///
1097    /// # let vertices = vec![
1098    /// #     Point::new(0.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0),
1099    /// #     Point::new(0.5, 1.0, 0.0), Point::new(0.5, 0.5, 1.0),
1100    /// # ];
1101    /// # let indices = vec![[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2]];
1102    /// #
1103    /// let decomposition = VHACD::decompose(
1104    ///     &VHACDParameters::default(),
1105    ///     &vertices,
1106    ///     &indices,
1107    ///     false, // voxel-to-primitive mapping not needed
1108    /// );
1109    ///
1110    /// // Compute voxel-based convex hulls with moderate downsampling
1111    /// let hulls = decomposition.compute_convex_hulls(4);
1112    ///
1113    /// for (i, (verts, tris)) in hulls.iter().enumerate() {
1114    ///     println!("Hull {}: {} vertices, {} triangles", i, verts.len(), tris.len());
1115    /// }
1116    /// # }
1117    /// ```
1118    ///
1119    /// ## Creating a Compound Shape for Collision Detection
1120    ///
1121    /// ```no_run
1122    /// # #[cfg(all(feature = "dim3", feature = "f32"))] {
1123    /// use parry3d::math::Point;
1124    /// use parry3d::shape::{SharedShape, Compound};
1125    /// use parry3d::transformation::vhacd::{VHACD, VHACDParameters};
1126    /// use parry3d::na::Isometry3;
1127    ///
1128    /// # let vertices = vec![
1129    /// #     Point::new(0.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0),
1130    /// #     Point::new(0.5, 1.0, 0.0), Point::new(0.5, 0.5, 1.0),
1131    /// # ];
1132    /// # let indices = vec![[0, 1, 2], [0, 2, 3], [0, 3, 1], [1, 3, 2]];
1133    /// #
1134    /// let decomposition = VHACD::decompose(
1135    ///     &VHACDParameters::default(),
1136    ///     &vertices,
1137    ///     &indices,
1138    ///     false,
1139    /// );
1140    ///
1141    /// let hulls = decomposition.compute_convex_hulls(4);
1142    /// # }
1143    /// ```
1144    ///
1145    /// # See Also
1146    ///
1147    /// - [`compute_exact_convex_hulls`](VHACD::compute_exact_convex_hulls): More accurate mesh-based hulls
1148    /// - [`voxel_parts`](VHACD::voxel_parts): Access the raw voxel data
1149    /// - [`SharedShape::convex_decomposition`](crate::shape::SharedShape::convex_decomposition): High-level API
1150    #[cfg(feature = "dim3")]
1151    pub fn compute_convex_hulls(
1152        &self,
1153        downsampling: u32,
1154    ) -> Vec<(Vec<Point<Real>>, Vec<[u32; DIM]>)> {
1155        let downsampling = downsampling.max(1);
1156        self.voxel_parts
1157            .iter()
1158            .map(|part| part.compute_convex_hull(downsampling))
1159            .collect()
1160    }
1161}
1162
1163fn compute_concavity(volume: Real, volume_ch: Real, volume0: Real) -> Real {
1164    (volume_ch - volume).abs() / volume0
1165}
1166
1167fn clip_mesh(
1168    points: &[Point<Real>],
1169    plane: &CutPlane,
1170    positive_part: &mut Vec<Point<Real>>,
1171    negative_part: &mut Vec<Point<Real>>,
1172) {
1173    for pt in points {
1174        let d = plane.abc.dot(&pt.coords) + plane.d;
1175
1176        if d > 0.0 {
1177            positive_part.push(*pt);
1178        } else if d < 0.0 {
1179            negative_part.push(*pt);
1180        } else {
1181            positive_part.push(*pt);
1182            negative_part.push(*pt);
1183        }
1184    }
1185}
1186
1187#[cfg(feature = "dim2")]
1188fn convex_hull(vertices: &[Point<Real>]) -> Vec<Point<Real>> {
1189    if vertices.len() > 1 {
1190        crate::transformation::convex_hull(vertices)
1191    } else {
1192        Vec::new()
1193    }
1194}
1195
1196#[cfg(feature = "dim3")]
1197fn convex_hull(vertices: &[Point<Real>]) -> (Vec<Point<Real>>, Vec<[u32; DIM]>) {
1198    if vertices.len() > 2 {
1199        crate::transformation::convex_hull(vertices)
1200    } else {
1201        (Vec::new(), Vec::new())
1202    }
1203}
1204
1205#[cfg(feature = "dim2")]
1206fn compute_volume(polygon: &[Point<Real>]) -> Real {
1207    if !polygon.is_empty() {
1208        crate::mass_properties::details::convex_polygon_area_and_center_of_mass(polygon).0
1209    } else {
1210        0.0
1211    }
1212}
1213
1214#[cfg(feature = "dim3")]
1215fn compute_volume(mesh: &(Vec<Point<Real>>, Vec<[u32; DIM]>)) -> Real {
1216    if !mesh.0.is_empty() {
1217        crate::mass_properties::details::trimesh_signed_volume_and_center_of_mass(&mesh.0, &mesh.1)
1218            .0
1219    } else {
1220        0.0
1221    }
1222}