parry3d/transformation/vhacd/
vhacd.rs1use 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
38pub struct VHACD {
40 voxel_parts: Vec<VoxelSet>,
42 volume_ch0: Real,
43 max_concavity: Real,
44}
45
46impl VHACD {
47 pub fn decompose(
59 params: &VHACDParameters,
60 points: &[Point<Real>],
61 indices: &[[u32; DIM]],
62 keep_voxel_to_primitives_map: bool,
63 ) -> Self {
64 let voxelized = VoxelizedVolume::voxelize(
70 points,
71 indices,
72 params.resolution,
73 params.fill_mode,
74 keep_voxel_to_primitives_map,
75 );
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 pub fn from_voxels(params: &VHACDParameters, voxels: VoxelSet) -> Self {
90 let mut result = Self {
91 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 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 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 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 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 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(); voxels.compute_bb(); let voxels_convex_hull = voxels.compute_convex_hull(params.convex_hull_downsampling); let volume_ch = compute_volume(&voxels_convex_hull);
290
291 if first_iteration {
293 self.volume_ch0 = volume_ch;
294 }
295
296 let concavity = compute_concavity(volume, volume_ch, self.volume_ch0);
298
299 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, 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 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 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 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 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 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 #[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 #[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 #[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 #[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}