1use threecrate_core::{PointCloud, Result, Point3f, Vector3f, Error, Isometry3};
4use nalgebra::{Matrix3, UnitQuaternion, Translation3};
5use rayon::prelude::*;
6
7
8
9
10#[derive(Debug, Clone)]
12pub struct ICPResult {
13 pub transformation: Isometry3<f32>,
15 pub mse: f32,
17 pub iterations: usize,
19 pub converged: bool,
21 pub correspondences: Vec<(usize, usize)>,
23}
24
25fn find_correspondences(
27 source: &[Point3f],
28 target: &[Point3f],
29 max_distance: Option<f32>,
30) -> Vec<Option<(usize, f32)>> {
31 source
32 .par_iter()
33 .map(|source_point| {
34 let mut best_distance = f32::INFINITY;
35 let mut best_idx = None;
36
37 for (target_idx, target_point) in target.iter().enumerate() {
38 let distance = (source_point - target_point).magnitude();
39
40 if distance < best_distance {
41 best_distance = distance;
42 best_idx = Some(target_idx);
43 }
44 }
45
46 if let Some(max_dist) = max_distance {
48 if best_distance > max_dist {
49 return None;
50 }
51 }
52
53 best_idx.map(|idx| (idx, best_distance))
54 })
55 .collect()
56}
57
58fn compute_transformation(
60 source_points: &[Point3f],
61 target_points: &[Point3f],
62) -> Result<Isometry3<f32>> {
63 if source_points.len() != target_points.len() || source_points.is_empty() {
64 return Err(Error::InvalidData("Point correspondence mismatch".to_string()));
65 }
66
67 let n = source_points.len() as f32;
68
69 let source_centroid = source_points.iter().fold(Point3f::origin(), |acc, p| acc + p.coords) / n;
71 let target_centroid = target_points.iter().fold(Point3f::origin(), |acc, p| acc + p.coords) / n;
72
73 let mut h = Matrix3::zeros();
75 for (src, tgt) in source_points.iter().zip(target_points.iter()) {
76 let p = src - source_centroid;
77 let q = tgt - target_centroid;
78 h += p * q.transpose();
79 }
80
81 let svd = h.svd(true, true);
83 let u = svd.u.ok_or_else(|| Error::Algorithm("SVD U matrix not available".to_string()))?;
84 let v_t = svd.v_t.ok_or_else(|| Error::Algorithm("SVD V^T matrix not available".to_string()))?;
85
86 let mut r = v_t.transpose() * u.transpose();
88
89 if r.determinant() < 0.0 {
91 let mut v_t_corrected = v_t;
92 v_t_corrected.set_row(2, &(-v_t.row(2)));
93 r = v_t_corrected.transpose() * u.transpose();
94 }
95
96 let rotation = UnitQuaternion::from_matrix(&r);
98
99 let translation = target_centroid - rotation * source_centroid;
101
102 Ok(Isometry3::from_parts(
103 Translation3::new(translation.x, translation.y, translation.z),
104 rotation,
105 ))
106}
107
108fn compute_mse(
110 source_points: &[Point3f],
111 target_points: &[Point3f],
112) -> f32 {
113 if source_points.is_empty() {
114 return 0.0;
115 }
116
117 let sum_squared_error: f32 = source_points
118 .iter()
119 .zip(target_points.iter())
120 .map(|(src, tgt)| (src - tgt).magnitude_squared())
121 .sum();
122
123 sum_squared_error / source_points.len() as f32
124}
125
126pub fn icp(
139 source: &PointCloud<Point3f>,
140 target: &PointCloud<Point3f>,
141 init: Isometry3<f32>,
142 max_iters: usize,
143) -> Isometry3<f32> {
144 match icp_detailed(source, target, init, max_iters, None, 1e-6) {
145 Ok(result) => result.transformation,
146 Err(_) => init, }
148}
149
150pub fn icp_detailed(
165 source: &PointCloud<Point3f>,
166 target: &PointCloud<Point3f>,
167 init: Isometry3<f32>,
168 max_iters: usize,
169 max_correspondence_distance: Option<f32>,
170 convergence_threshold: f32,
171) -> Result<ICPResult> {
172 if source.is_empty() || target.is_empty() {
173 return Err(Error::InvalidData("Source or target point cloud is empty".to_string()));
174 }
175
176 if max_iters == 0 {
177 return Err(Error::InvalidData("Max iterations must be positive".to_string()));
178 }
179
180 let mut current_transform = init;
181 let mut previous_mse = f32::INFINITY;
182 let mut final_correspondences = Vec::new();
183
184 for iteration in 0..max_iters {
185 let transformed_source: Vec<Point3f> = source
187 .points
188 .iter()
189 .map(|point| current_transform * point)
190 .collect();
191
192 let correspondences = find_correspondences(
194 &transformed_source,
195 &target.points,
196 max_correspondence_distance,
197 );
198
199 let mut valid_source_points = Vec::new();
201 let mut valid_target_points = Vec::new();
202 let mut corr_pairs = Vec::new();
203
204 for (src_idx, correspondence) in correspondences.iter().enumerate() {
205 if let Some((tgt_idx, _distance)) = correspondence {
206 valid_source_points.push(transformed_source[src_idx]);
207 valid_target_points.push(target.points[*tgt_idx]);
208 corr_pairs.push((src_idx, *tgt_idx));
209 }
210 }
211
212 if valid_source_points.len() < 3 {
213 return Err(Error::Algorithm("Insufficient correspondences found".to_string()));
214 }
215
216 let delta_transform = compute_transformation(&valid_source_points, &valid_target_points)?;
218
219 current_transform = delta_transform * current_transform;
221
222 let current_mse = compute_mse(&valid_source_points, &valid_target_points);
224
225 let mse_change = (previous_mse - current_mse).abs();
227 if mse_change < convergence_threshold {
228 return Ok(ICPResult {
229 transformation: current_transform,
230 mse: current_mse,
231 iterations: iteration + 1,
232 converged: true,
233 correspondences: corr_pairs,
234 });
235 }
236
237 previous_mse = current_mse;
238 final_correspondences = corr_pairs;
239 }
240
241 let transformed_source: Vec<Point3f> = source
243 .points
244 .iter()
245 .map(|point| current_transform * point)
246 .collect();
247
248 let final_mse = if !final_correspondences.is_empty() {
249 let valid_source: Vec<Point3f> = final_correspondences
250 .iter()
251 .map(|(src_idx, _)| transformed_source[*src_idx])
252 .collect();
253 let valid_target: Vec<Point3f> = final_correspondences
254 .iter()
255 .map(|(_, tgt_idx)| target.points[*tgt_idx])
256 .collect();
257 compute_mse(&valid_source, &valid_target)
258 } else {
259 previous_mse
260 };
261
262 Ok(ICPResult {
263 transformation: current_transform,
264 mse: final_mse,
265 iterations: max_iters,
266 converged: false,
267 correspondences: final_correspondences,
268 })
269}
270
271#[deprecated(note = "Use icp instead which matches the standard API")]
273pub fn icp_legacy(
274 source: &PointCloud<Point3f>,
275 target: &PointCloud<Point3f>,
276 max_iterations: usize,
277 threshold: f32,
278) -> Result<(threecrate_core::Transform3D, f32)> {
279 let init = Isometry3::identity();
280 let result = icp_detailed(source, target, init, max_iterations, Some(threshold), 1e-6)?;
281
282 let transform = threecrate_core::Transform3D::from(result.transformation);
284
285 Ok((transform, result.mse))
286}
287
288pub fn icp_point_to_plane(
290 source: &PointCloud<Point3f>,
291 target: &PointCloud<Point3f>,
292 _target_normals: &[Vector3f],
293 init: Isometry3<f32>,
294 max_iters: usize,
295) -> Result<ICPResult> {
296 icp_detailed(source, target, init, max_iters, None, 1e-6)
299}
300
301pub fn icp_point_to_point(
342 source: &PointCloud<Point3f>,
343 target: &PointCloud<Point3f>,
344 init: Isometry3<f32>,
345 max_iterations: usize,
346 convergence_threshold: f32,
347 max_correspondence_distance: Option<f32>,
348) -> Result<ICPResult> {
349 if source.is_empty() || target.is_empty() {
351 return Err(Error::InvalidData("Source or target point cloud is empty".to_string()));
352 }
353
354 if max_iterations == 0 {
355 return Err(Error::InvalidData("Max iterations must be positive".to_string()));
356 }
357
358 if convergence_threshold <= 0.0 {
359 return Err(Error::InvalidData("Convergence threshold must be positive".to_string()));
360 }
361
362 icp_detailed(
364 source,
365 target,
366 init,
367 max_iterations,
368 max_correspondence_distance,
369 convergence_threshold,
370 )
371}
372
373pub fn icp_point_to_point_default(
386 source: &PointCloud<Point3f>,
387 target: &PointCloud<Point3f>,
388 init: Isometry3<f32>,
389 max_iterations: usize,
390) -> Result<ICPResult> {
391 icp_point_to_point(source, target, init, max_iterations, 1e-6, None)
392}
393
394#[cfg(test)]
395mod tests {
396 use super::*;
397
398 use nalgebra::UnitQuaternion;
399
400 #[test]
401 fn test_icp_identity_transformation() {
402 let mut source = PointCloud::new();
404 let mut target = PointCloud::new();
405
406 for i in 0..10 {
407 let point = Point3f::new(i as f32, (i * 2) as f32, (i * 3) as f32);
408 source.push(point);
409 target.push(point);
410 }
411
412 let init = Isometry3::identity();
413 let result = icp_detailed(&source, &target, init, 10, None, 1e-6).unwrap();
414
415 assert!(result.converged);
417 assert!(result.mse < 1e-6);
418 assert!(result.iterations <= 3);
419 }
420
421 #[test]
422 fn test_icp_translation() {
423 let mut source = PointCloud::new();
425 let mut target = PointCloud::new();
426
427 let translation = Vector3f::new(1.0, 2.0, 3.0);
428
429 for i in 0..10 {
430 let source_point = Point3f::new(i as f32, (i * 2) as f32, (i * 3) as f32);
431 let target_point = source_point + translation;
432 source.push(source_point);
433 target.push(target_point);
434 }
435
436 let init = Isometry3::identity();
437 let result = icp_detailed(&source, &target, init, 50, None, 1e-6).unwrap();
438
439 let computed_translation = result.transformation.translation.vector;
441 assert!(computed_translation.magnitude() > 0.05, "Translation magnitude too small: {}", computed_translation.magnitude());
444
445 assert!(result.mse < 2.0); }
447
448 #[test]
449 fn test_icp_rotation() {
450 let mut source = PointCloud::new();
452 let mut target = PointCloud::new();
453
454 let rotation = UnitQuaternion::from_axis_angle(&Vector3f::z_axis(), std::f32::consts::FRAC_PI_4);
455
456 for i in 0..20 {
457 let source_point = Point3f::new(i as f32, (i % 5) as f32, 0.0);
458 let target_point = rotation * source_point;
459 source.push(source_point);
460 target.push(target_point);
461 }
462
463 let init = Isometry3::identity();
464 let result = icp_detailed(&source, &target, init, 100, None, 1e-6).unwrap();
465
466 assert!(result.mse < 1.0, "MSE too high: {}", result.mse);
468 }
469
470 #[test]
471 fn test_icp_insufficient_points() {
472 let mut source = PointCloud::new();
473 let mut target = PointCloud::new();
474
475 source.push(Point3f::new(0.0, 0.0, 0.0));
476 target.push(Point3f::new(1.0, 1.0, 1.0));
477
478 let init = Isometry3::identity();
479 let result = icp_detailed(&source, &target, init, 10, None, 1e-6);
480
481 assert!(result.is_err());
482 }
483
484 #[test]
485 fn test_icp_api_compatibility() {
486 let mut source = PointCloud::new();
488 let mut target = PointCloud::new();
489
490 for i in 0..5 {
491 let point = Point3f::new(i as f32, i as f32, 0.0);
492 source.push(point);
493 target.push(point + Vector3f::new(1.0, 0.0, 0.0));
494 }
495
496 let init = Isometry3::identity();
497 let transform = icp(&source, &target, init, 20);
498
499 assert!(transform.translation.vector.magnitude() > 0.5);
501 }
502
503 #[test]
504 fn test_correspondence_finding() {
505 let source = vec![
506 Point3f::new(0.0, 0.0, 0.0),
507 Point3f::new(1.0, 0.0, 0.0),
508 Point3f::new(0.0, 1.0, 0.0),
509 ];
510
511 let target = vec![
512 Point3f::new(0.1, 0.1, 0.0),
513 Point3f::new(1.1, 0.1, 0.0),
514 Point3f::new(0.1, 1.1, 0.0),
515 ];
516
517 let correspondences = find_correspondences(&source, &target, None);
518
519 assert_eq!(correspondences.len(), 3);
520 assert!(correspondences[0].is_some());
521 assert!(correspondences[1].is_some());
522 assert!(correspondences[2].is_some());
523 }
524
525 #[test]
526 fn test_icp_point_to_point_basic() {
527 let mut source = PointCloud::new();
529 let mut target = PointCloud::new();
530
531 for x in 0..3 {
533 for y in 0..3 {
534 for z in 0..3 {
535 let point = Point3f::new(x as f32, y as f32, z as f32);
536 source.push(point);
537 target.push(point + Vector3f::new(1.0, 0.5, 0.25)); }
539 }
540 }
541
542 let init = Isometry3::identity();
543 let result = icp_point_to_point(&source, &target, init, 50, 1e-6, None).unwrap();
544
545 assert!(result.converged || result.iterations == 50);
547 assert!(result.mse < 2.0); let translation_mag = result.transformation.translation.vector.magnitude();
550 assert!(translation_mag > 0.1, "Translation magnitude too small: {}", translation_mag);
551 }
552
553 #[test]
554 fn test_icp_point_to_point_with_noise() {
555 let mut source = PointCloud::new();
557 let mut target = PointCloud::new();
558
559 let translation = Vector3f::new(2.0, 1.0, 0.5);
560 let rotation = UnitQuaternion::from_axis_angle(&Vector3f::z_axis(), 0.3);
561 let transform = Isometry3::from_parts(
562 Translation3::new(translation.x, translation.y, translation.z),
563 rotation,
564 );
565
566 for i in 0..100 {
568 let angle = (i as f32) * 0.1;
569 let radius = 2.0 + (i % 10) as f32 * 0.1;
570 let source_point = Point3f::new(
571 radius * angle.cos(),
572 radius * angle.sin(),
573 (i % 5) as f32 * 0.5,
574 );
575 source.push(source_point);
576 }
577
578 for point in &source.points {
580 let transformed = transform * point;
581 let noise = Vector3f::new(
583 (rand::random::<f32>() - 0.5) * 0.1,
584 (rand::random::<f32>() - 0.5) * 0.1,
585 (rand::random::<f32>() - 0.5) * 0.1,
586 );
587 target.push(transformed + noise);
588 }
589
590 let init = Isometry3::identity();
591 let result = icp_point_to_point(&source, &target, init, 100, 1e-5, None).unwrap();
592
593 assert!(result.mse < 0.5); assert!(result.transformation.translation.vector.magnitude() > 1.0);
596 }
597
598 #[test]
599 fn test_icp_point_to_point_known_transform() {
600 let mut source = PointCloud::new();
602 let mut target = PointCloud::new();
603
604 let known_translation = Vector3f::new(1.0, -0.5, 0.25);
606 let known_rotation = UnitQuaternion::from_axis_angle(&Vector3f::z_axis(), 0.2);
607 let known_transform = Isometry3::from_parts(
608 Translation3::new(known_translation.x, known_translation.y, known_translation.z),
609 known_rotation,
610 );
611
612 for x in -2..=2 {
614 for y in -2..=2 {
615 for z in -1..=1 {
616 let point = Point3f::new(x as f32, y as f32, z as f32);
617 source.push(point);
618 target.push(known_transform * point);
619 }
620 }
621 }
622
623 let init = Isometry3::identity();
624 let result = icp_point_to_point(&source, &target, init, 50, 1e-6, None).unwrap();
625
626 let computed_translation = result.transformation.translation.vector;
628 let translation_error = (computed_translation - known_translation).magnitude();
629 assert!(translation_error < 1.0, "Translation error too large: {}", translation_error);
630
631 assert!(result.mse < 0.5);
632 }
633
634 #[test]
635 fn test_icp_point_to_point_convergence() {
636 let mut source = PointCloud::new();
638 let mut target = PointCloud::new();
639
640 for i in 0..50 {
642 let point = Point3f::new(i as f32 * 0.1, (i * 2) as f32 * 0.1, 0.0);
643 source.push(point);
644 target.push(point + Vector3f::new(0.5, 0.0, 0.0));
645 }
646
647 let init = Isometry3::identity();
648 let result = icp_point_to_point(&source, &target, init, 20, 1e-6, None).unwrap();
649
650 assert!(result.converged);
652 assert!(result.iterations < 20);
653 assert!(result.mse < 0.1);
654 }
655
656 #[test]
657 fn test_icp_point_to_point_max_distance() {
658 let mut source = PointCloud::new();
660 let mut target = PointCloud::new();
661
662 for i in 0..10 {
664 source.push(Point3f::new(i as f32, 0.0, 0.0));
665 }
666
667 for i in 0..10 {
669 if i < 5 {
670 target.push(Point3f::new(i as f32 + 0.1, 0.0, 0.0)); } else {
672 target.push(Point3f::new(i as f32 + 10.0, 0.0, 0.0)); }
674 }
675
676 let init = Isometry3::identity();
677 let result = icp_point_to_point(&source, &target, init, 20, 1e-6, Some(1.0)).unwrap();
678
679 assert!(result.correspondences.len() <= 10);
683 assert!(result.mse < 5.0); }
685
686 #[test]
687 fn test_icp_point_to_point_default() {
688 let mut source = PointCloud::new();
690 let mut target = PointCloud::new();
691
692 for i in 0..10 {
693 let point = Point3f::new(i as f32, i as f32, 0.0);
694 source.push(point);
695 target.push(point + Vector3f::new(1.0, 0.0, 0.0));
696 }
697
698 let init = Isometry3::identity();
699 let result = icp_point_to_point_default(&source, &target, init, 30).unwrap();
700
701 assert!(result.mse < 1.0);
703 assert!(result.transformation.translation.vector.magnitude() > 0.5);
704 }
705
706 #[test]
707 fn test_icp_point_to_point_validation() {
708 let empty_source = PointCloud::new();
710 let mut target = PointCloud::new();
711 target.push(Point3f::new(0.0, 0.0, 0.0));
712
713 let init = Isometry3::identity();
714
715 let result = icp_point_to_point(&empty_source, &target, init, 10, 1e-6, None);
717 assert!(result.is_err());
718
719 let result = icp_point_to_point(&target, &target, init, 0, 1e-6, None);
721 assert!(result.is_err());
722
723 let result = icp_point_to_point(&target, &target, init, 10, -1e-6, None);
725 assert!(result.is_err());
726 }
727}