pasture_algorithms/
normal_estimation.rs

1// The normal estimation algorithm is inspired by the PCL library (https://pointclouds.org/)
2use core::panic;
3use kd_tree::{self, KdPoint, KdTree};
4use num_traits::{self};
5use pasture_core::containers::{
6    BorrowedBuffer, BorrowedBufferExt, BorrowedMutBufferExt, HashMapBuffer, OwningBuffer,
7};
8use pasture_core::layout::{attributes::POSITION_3D, PointType};
9use pasture_core::nalgebra::{DMatrix, Vector3};
10use std::result::Result;
11
12/// Normal Estimation Algorithm
13/// returns a vector of quintuplets where each quintuplet has the following values: (current point, normal_vector, curvature).
14/// It iterates over all points in the buffer and constructs new point buffers of size k_nn.
15/// Make sure that knn >= 3 holds as it is not possible to construct a plane over less than three point
16///
17/// # Panics
18///
19/// Panics if the number of k nearest neighbors is less than 3 or the point cloud has less than 3 elements.
20///
21/// # Examples
22///
23/// ```
24/// # use kd_tree::{self, KdPoint};
25/// # use pasture_core::nalgebra::Vector3;
26/// # use pasture_core::{containers::*, layout::PointType};
27/// # use pasture_algorithms::normal_estimation::compute_normals;
28/// # use pasture_derive::PointType;
29/// # use typenum;
30///
31/// #[repr(C, packed)]
32/// #[derive(PointType, Debug, Clone, Copy, bytemuck::AnyBitPattern, bytemuck::NoUninit)]
33/// struct SimplePoint {
34///     #[pasture(BUILTIN_POSITION_3D)]
35///     pub position: Vector3<f64>,
36///     #[pasture(BUILTIN_INTENSITY)]
37///     pub intensity: u16,
38/// }
39/// impl KdPoint for SimplePoint {
40///     type Scalar = f64;
41///     type Dim = typenum::U3;
42///     fn at(&self, k: usize) -> f64 {
43///         let position = self.position;
44///         position[k]
45///     }
46/// }
47/// fn main() {
48///     let points = vec![
49///         SimplePoint {
50///             position: Vector3::new(11.0, 22.0, 154.0),
51///             intensity: 42,
52///         },
53///         SimplePoint {
54///             position: Vector3::new(12.0, 23.0, 0.0),
55///             intensity: 84,
56///         },
57///         SimplePoint {
58///             position: Vector3::new(103.0, 84.0, 2.0),
59///             intensity: 84,
60///         },
61///         SimplePoint {
62///             position: Vector3::new(101.0, 0.0, 1.0),
63///             intensity: 84,
64///         },
65///     ];
66///
67///     let interleaved = points.into_iter().collect::<VectorBuffer>();
68///
69///     let solution_vec = compute_normals::<VectorBuffer, SimplePoint>(&interleaved, 4);
70///     for solution in solution_vec {
71///    println!(
72///        "Point: {:?}, n_x: {}, n_y: {}, n_z: {}, curvature: {}",
73///        solution.0, solution.0[0], solution.0[1], solution.0[2], solution.1
74///    );
75/// }
76/// }
77/// ```
78
79pub fn compute_normals<'a, T: BorrowedBuffer<'a>, P: PointType + KdPoint + Copy>(
80    point_cloud: &'a T,
81    k_nn: usize,
82) -> Vec<(Vector3<f64>, f64)>
83where
84    P::Scalar: num_traits::Float,
85{
86    if point_cloud.len() < 3 {
87        panic!("The point cloud is too small. Please use a point cloud that has 3 or more points!");
88    }
89    if k_nn < 3 {
90        panic!("The k nearest neigbors attribute is too small!");
91    }
92
93    // this is the solution that will be returned
94    let mut points_with_normals_curvature = vec![];
95
96    // transform point cloud in vector of points
97    let mut points: Vec<[f64; 3]> = vec![];
98    for point in point_cloud.view_attribute::<Vector3<f64>>(&POSITION_3D) {
99        points.push(*point.as_ref());
100    }
101
102    // construct kd tree over the vector of points.
103    let cloud_as_kd_tree = KdTree::build_by_ordered_float(points);
104
105    // iterate over all points in the point cloud and and calculate the k nearest neighbors with the constructed kd tree
106    for point in point_cloud.view_attribute::<Vector3<f64>>(&POSITION_3D) {
107        let r: &[f64; 3] = point.as_ref();
108        let nearest_points = cloud_as_kd_tree.nearests(r, k_nn);
109
110        // stores the k nearest neighbors in a PointStorage
111        let mut k_nn_buffer =
112            HashMapBuffer::with_capacity(nearest_points.len(), point_cloud.point_layout().clone());
113        k_nn_buffer.resize(nearest_points.len());
114
115        for (index, point) in nearest_points.iter().enumerate() {
116            k_nn_buffer.view_attribute_mut(&POSITION_3D).set_at(
117                index,
118                Vector3::new(point.item[0], point.item[1], point.item[2]),
119            );
120        }
121
122        // coordinates of the surface normal and curvature
123        let (surface_normal, curvature) = normal_estimation(&k_nn_buffer);
124
125        // solution vector
126        points_with_normals_curvature.push((surface_normal, curvature));
127    }
128
129    points_with_normals_curvature
130}
131
132/// checks whether a given point cloud has points with coordinates that are Not a Number
133fn is_dense<'a, T: BorrowedBuffer<'a>>(point_cloud: &'a T) -> bool {
134    for point in point_cloud.view_attribute::<Vector3<f64>>(&POSITION_3D) {
135        if point.x.is_nan() || point.y.is_nan() || point.z.is_nan() {
136            return false;
137        }
138    }
139    true
140}
141
142/// checks whether a given point has finite coordinates
143fn is_finite(point: &Vector3<f64>) -> bool {
144    if point.x.is_finite() && point.y.is_finite() && point.z.is_finite() {
145        return true;
146    }
147    false
148}
149
150/// Computes the centroid for a given point cloud.
151/// The centroid is the point that has the same distance to all other points in the point cloud.
152/// Iterates over all points in the 'point_cloud'.
153///
154/// # Panics
155///
156/// If the point_cloud is empty
157///
158/// # Examples
159///
160/// ```
161/// # use pasture_core::{containers::*, layout::PointType};
162/// # use pasture_core::nalgebra::Vector3;
163/// # use pasture_derive::PointType;
164/// # use pasture_algorithms::normal_estimation::compute_centroid;
165///
166/// #[repr(C, packed)]
167/// #[derive(PointType, Debug, Clone, Copy, bytemuck::AnyBitPattern, bytemuck::NoUninit)]
168/// struct SimplePoint {
169///     #[pasture(BUILTIN_POSITION_3D)]
170///     pub position: Vector3<f64>,
171///     #[pasture(BUILTIN_INTENSITY)]
172///     pub intensity: u16,
173/// }
174/// let points: Vec<SimplePoint> = vec![
175///     SimplePoint {
176///         position: Vector3::new(1.0, 0.0, 0.0),
177///         intensity: 42,
178///     },
179///     SimplePoint {
180///         position: Vector3::new(0.0, 1.0, 0.0),
181///         intensity: 84,
182///     },
183///     SimplePoint {
184///         position: Vector3::new(1.0, 1.0, 0.0),
185///         intensity: 84,
186///     },
187///     SimplePoint {
188///         position: Vector3::new(-1.0, 0.0, 0.0),
189///         intensity: 84,
190///     },
191/// ];
192
193/// let interleaved = points.into_iter().collect::<VectorBuffer>();
194
195/// let centroid = compute_centroid(&interleaved);
196///
197/// ```
198pub fn compute_centroid<'a, T: BorrowedBuffer<'a>>(point_cloud: &'a T) -> Vector3<f64> {
199    if point_cloud.is_empty() {
200        panic!("The point cloud is empty!");
201    }
202
203    let mut centroid = Vector3::<f64>::zeros();
204    let mut temp_centroid = Vector3::<f64>::zeros();
205
206    if is_dense(point_cloud) {
207        // add all points up
208        for point in point_cloud.view_attribute::<Vector3<f64>>(&POSITION_3D) {
209            temp_centroid[0] += point.x;
210            temp_centroid[1] += point.y;
211            temp_centroid[2] += point.z;
212        }
213
214        //normalize over all points
215        centroid[0] = temp_centroid[0] / point_cloud.len() as f64;
216        centroid[1] = temp_centroid[1] / point_cloud.len() as f64;
217        centroid[2] = temp_centroid[2] / point_cloud.len() as f64;
218    } else {
219        let mut points_in_cloud = 0;
220        for point in point_cloud.view_attribute::<Vector3<f64>>(&POSITION_3D) {
221            if is_finite(&point) {
222                // add all points up
223                temp_centroid[0] += point.x;
224                temp_centroid[1] += point.y;
225                temp_centroid[2] += point.z;
226                points_in_cloud += 1;
227            }
228        }
229
230        // normalize over all points
231        centroid[0] = temp_centroid[0] / points_in_cloud as f64;
232        centroid[1] = temp_centroid[1] / points_in_cloud as f64;
233        centroid[2] = temp_centroid[2] / points_in_cloud as f64;
234    }
235
236    centroid
237}
238
239/// compute the covariance matrix for a given point cloud which is a measure of spread out the points are
240fn compute_covariance_matrix<'a, T: BorrowedBuffer<'a>>(
241    point_cloud: &'a T,
242) -> Result<DMatrix<f64>, &'static str> {
243    let mut covariance_matrix = DMatrix::<f64>::zeros(3, 3);
244    let mut point_count = 0;
245
246    // compute the centroid of the point cloud
247    let centroid = compute_centroid(point_cloud);
248
249    if is_dense(point_cloud) {
250        point_count = point_cloud.len();
251        let mut diff_mean = Vector3::<f64>::zeros();
252        for point in point_cloud.view_attribute::<Vector3<f64>>(&POSITION_3D) {
253            // calculate difference from the centroid for each point
254            diff_mean[0] = point.x - centroid[0];
255            diff_mean[1] = point.y - centroid[1];
256            diff_mean[2] = point.z - centroid[2];
257
258            covariance_matrix[(1, 1)] += diff_mean[1] * diff_mean[1];
259            covariance_matrix[(1, 2)] += diff_mean[1] * diff_mean[2];
260            covariance_matrix[(2, 2)] += diff_mean[2] * diff_mean[2];
261
262            let diff_x = diff_mean[0];
263            diff_mean.iter_mut().for_each(|x| *x *= diff_x);
264
265            covariance_matrix[(0, 0)] += diff_mean[0];
266            covariance_matrix[(0, 1)] += diff_mean[1];
267            covariance_matrix[(0, 2)] += diff_mean[2];
268        }
269    } else {
270        // in this case we dont know the number of points in the point cloud that are finite
271        for point in point_cloud.view_attribute::<Vector3<f64>>(&POSITION_3D) {
272            if !is_finite(&point) {
273                continue;
274            }
275            // only compute the covariance matrix for finite points
276            let mut diff_mean = Vector3::<f64>::zeros();
277            // calculate difference from the centroid for each point
278            diff_mean[0] = point.x - centroid[0];
279            diff_mean[1] = point.y - centroid[1];
280            diff_mean[2] = point.z - centroid[2];
281
282            covariance_matrix[(1, 1)] += diff_mean[1] * diff_mean[1];
283            covariance_matrix[(1, 2)] += diff_mean[1] * diff_mean[2];
284            covariance_matrix[(2, 2)] += diff_mean[2] * diff_mean[2];
285
286            let diff_x = diff_mean[0];
287            diff_mean.iter_mut().for_each(|x| *x *= diff_x);
288
289            covariance_matrix[(0, 0)] += diff_mean[0];
290            covariance_matrix[(0, 1)] += diff_mean[1];
291            covariance_matrix[(0, 2)] += diff_mean[2];
292            point_count += 1;
293        }
294    }
295
296    if point_count < 3 {
297        return Err("The number of valid (finite and non-NaN values) points in a k nearest neighborhood is not enough to span a plane!");
298    }
299
300    covariance_matrix[(1, 0)] = covariance_matrix[(0, 1)];
301    covariance_matrix[(2, 0)] = covariance_matrix[(0, 2)];
302    covariance_matrix[(2, 1)] = covariance_matrix[(1, 2)];
303
304    Ok(covariance_matrix)
305}
306
307/// find the eigen value solution if the highest degree of the polynomial is 2
308fn solve_polynomial_quadratic(coefficient_2: f64, coefficient_1: f64) -> Vector3<f64> {
309    let mut eigen_values = Vector3::<f64>::zeros();
310
311    eigen_values[0] = 0.0;
312
313    let mut delta = coefficient_2 * coefficient_2 - 4.0 * coefficient_1;
314
315    if delta < 0.0 {
316        delta = 0.0;
317    }
318
319    let sqrt_delta = f64::sqrt(delta);
320
321    eigen_values[2] = 0.5 * (coefficient_2 + sqrt_delta);
322    eigen_values[1] = 0.5 * (coefficient_2 - sqrt_delta);
323
324    eigen_values
325}
326
327/// solve the polynomial to find the eigen values for a given covariance matrix
328fn solve_polynomial(covariance_matrix: &DMatrix<f64>) -> Vector3<f64> {
329    let coefficient_0 = covariance_matrix[(0, 0)]
330        * covariance_matrix[(1, 1)]
331        * covariance_matrix[(2, 2)]
332        + 2.0 * covariance_matrix[(0, 1)] * covariance_matrix[(0, 2)] * covariance_matrix[(1, 2)]
333        - covariance_matrix[(0, 0)] * covariance_matrix[(1, 2)] * covariance_matrix[(1, 2)]
334        - covariance_matrix[(1, 1)] * covariance_matrix[(0, 2)] * covariance_matrix[(0, 2)]
335        - covariance_matrix[(2, 2)] * covariance_matrix[(0, 1)] * covariance_matrix[(0, 1)];
336    let coefficient_1 = covariance_matrix[(0, 0)] * covariance_matrix[(1, 1)]
337        - covariance_matrix[(0, 1)] * covariance_matrix[(0, 1)]
338        + covariance_matrix[(0, 0)] * covariance_matrix[(2, 2)]
339        - covariance_matrix[(0, 2)] * covariance_matrix[(0, 2)]
340        + covariance_matrix[(1, 1)] * covariance_matrix[(2, 2)]
341        - covariance_matrix[(1, 2)] * covariance_matrix[(1, 2)];
342    let coefficient_2 =
343        covariance_matrix[(0, 0)] + covariance_matrix[(1, 1)] + covariance_matrix[(2, 2)];
344
345    // check if one eigen value solution is zero
346    if coefficient_0.abs() < std::f64::EPSILON {
347        solve_polynomial_quadratic(coefficient_2, coefficient_1)
348    } else {
349        let mut eigen_values = Vector3::<f64>::zeros();
350
351        let one_third = 1.0 / 3.0;
352        let sqrt_3 = f64::sqrt(3.0);
353
354        let coefficient_2_third = coefficient_2 * one_third;
355        let mut alpha_third = (coefficient_1 - coefficient_2 * coefficient_2_third) * one_third;
356        if alpha_third > 0.0 {
357            alpha_third = 0.0;
358        }
359
360        let half_beta = 0.5
361            * (coefficient_0
362                + coefficient_2_third
363                    * (2.0 * coefficient_2_third * coefficient_2_third - coefficient_1));
364
365        let mut q = half_beta * half_beta + alpha_third * alpha_third * alpha_third;
366        if q > 0.0 {
367            q = 0.0;
368        }
369
370        // calculate eigen values
371        let rho = f64::sqrt(-alpha_third);
372        let theta = f64::atan2(f64::sqrt(-q), half_beta) * one_third;
373        let cosine_of_theta = f64::cos(theta);
374        let sine_of_theta = f64::sin(theta);
375
376        eigen_values[0] = coefficient_2_third + 2.0 * rho * cosine_of_theta;
377        eigen_values[1] = coefficient_2_third - rho * (cosine_of_theta + sqrt_3 * sine_of_theta);
378        eigen_values[2] = coefficient_2_third - rho * (cosine_of_theta - sqrt_3 * sine_of_theta);
379
380        // sort increasing so that eigen_values[0] is the smallest eigen value
381        eigen_values
382            .as_mut_slice()
383            .sort_by(|a, b| a.partial_cmp(b).unwrap());
384
385        // if the smallest eigen value is zero or less the solution is a quadratic
386        if eigen_values[0] <= 0.0 {
387            solve_polynomial_quadratic(coefficient_2, coefficient_1)
388        } else {
389            eigen_values
390        }
391    }
392}
393
394/// calculates the largest eigen vector for a given matrix
395fn get_largest_eigen_vector(scaled_matrix: &DMatrix<f64>) -> Vector3<f64> {
396    let rows = vec![
397        scaled_matrix.row(0).cross(&scaled_matrix.row(1)),
398        scaled_matrix.row(0).cross(&scaled_matrix.row(2)),
399        scaled_matrix.row(1).cross(&scaled_matrix.row(2)),
400    ];
401
402    let mut cross_product = DMatrix::<f64>::zeros(3, 3);
403
404    // write rows of cross product
405    for it in cross_product.row_iter_mut().zip(rows) {
406        let (mut cross_row, row) = it;
407        // row from rows vector is written to the row of the cross product
408        cross_row.copy_from(&row);
409    }
410
411    // find largest eigen vector
412    let mut largest_eigen_vec = cross_product.row(0);
413    for row in cross_product.row_iter() {
414        if row.norm() > largest_eigen_vec.norm() {
415            largest_eigen_vec = row;
416        }
417    }
418
419    // set eigen vector to largest vector
420    let mut eigen_vector = Vector3::<f64>::zeros();
421    for i in 0..eigen_vector.len() {
422        eigen_vector[i] = largest_eigen_vec[i];
423    }
424
425    eigen_vector
426}
427
428/// for a given 3x3 matrix the functions calculates the eigen vector of the smallest eigen value that can be found
429fn eigen_3x3(covariance_matrix: &DMatrix<f64>) -> (f64, Vector3<f64>) {
430    let scale = covariance_matrix.abs().max();
431    let mut covariance_matrix_scaled = DMatrix::<f64>::zeros(3, 3);
432    for i in 0..covariance_matrix.len() {
433        covariance_matrix_scaled[i] = covariance_matrix[i] / scale;
434    }
435
436    // scale the matrix down
437    for (index, value) in covariance_matrix.iter().enumerate() {
438        covariance_matrix_scaled[index] = value / scale;
439    }
440
441    let eigen_values = solve_polynomial(covariance_matrix);
442    // undo scale for smallest eigen vector
443    let eigen_value = eigen_values[0] * scale;
444
445    // subtract the smallest eigen value from the diagonal of the matrix
446    covariance_matrix_scaled
447        .diagonal()
448        .iter_mut()
449        .for_each(|x| *x -= eigen_values[0]);
450
451    let eigen_vector = get_largest_eigen_vector(&covariance_matrix_scaled);
452    (eigen_value, eigen_vector)
453}
454
455/// calculates the orientation of the surface as a normal vector with components n_x, n_y and n_z as well as the curvature of the surface for a given covariance matrix
456fn solve_plane_parameter(covariance_matrix: &DMatrix<f64>) -> (Vector3<f64>, f64) {
457    let (eigen_value, eigen_vector) = eigen_3x3(covariance_matrix);
458
459    let eigen_sum = covariance_matrix[0] + covariance_matrix[4] + covariance_matrix[8];
460    let curvature = if eigen_sum != 0.0 {
461        (eigen_value / eigen_sum).abs()
462    } else {
463        0.0
464    };
465
466    (eigen_vector, curvature)
467}
468
469/// calculates the normal vectors and the curvature of the surface for the given point cloud
470fn normal_estimation<'a, T: BorrowedBuffer<'a>>(point_cloud: &'a T) -> (Vector3<f64>, f64) {
471    let covariance_matrix = compute_covariance_matrix(point_cloud).unwrap();
472
473    let (eigen_vector, curvature) = solve_plane_parameter(&covariance_matrix);
474
475    (eigen_vector, curvature)
476}
477
478#[cfg(test)]
479mod tests {
480
481    use pasture_core::{containers::VectorBuffer, nalgebra::Matrix3, nalgebra::Vector3};
482    use pasture_derive::PointType;
483
484    use super::*;
485
486    #[repr(C, packed)]
487    #[derive(PointType, Debug, Clone, Copy, bytemuck::AnyBitPattern, bytemuck::NoUninit)]
488    pub struct SimplePoint {
489        #[pasture(BUILTIN_POSITION_3D)]
490        pub position: Vector3<f64>,
491        #[pasture(BUILTIN_INTENSITY)]
492        pub intensity: u16,
493    }
494    impl KdPoint for SimplePoint {
495        type Scalar = f64;
496        type Dim = typenum::U3;
497        fn at(&self, k: usize) -> f64 {
498            let position = self.position;
499            position[k]
500        }
501    }
502
503    #[test]
504    fn test_compute_normal_sub() {
505        let points: Vec<SimplePoint> = vec![
506            SimplePoint {
507                position: Vector3::new(1.0, 0.0, 0.0),
508                intensity: 42,
509            },
510            SimplePoint {
511                position: Vector3::new(0.0, 1.0, 0.0),
512                intensity: 84,
513            },
514            SimplePoint {
515                position: Vector3::new(1.0, 1.0, 0.0),
516                intensity: 84,
517            },
518            SimplePoint {
519                position: Vector3::new(-1.0, 0.0, 0.0),
520                intensity: 84,
521            },
522        ];
523
524        let interleaved = points.into_iter().collect::<VectorBuffer>();
525
526        let centroid = compute_centroid(&interleaved);
527        let result_centroid = Vector3::<f64>::new(0.25, 0.5, 0.0);
528        assert_eq!(centroid, result_centroid);
529
530        let covariance_matrix = compute_covariance_matrix(&interleaved).unwrap();
531        let result = Matrix3::new(
532            0.6875 * 4.0,
533            0.125 * 4.0,
534            0.0,
535            0.125 * 4.0,
536            0.25 * 4.0,
537            0.0,
538            0.0,
539            0.0,
540            0.0,
541        );
542        assert_eq!(covariance_matrix, result);
543
544        let (normal_vector, curvature) = solve_plane_parameter(&covariance_matrix);
545
546        assert_eq!(normal_vector[0], 0.0);
547        assert_eq!(normal_vector[1], 0.0);
548        assert_ne!(normal_vector[2], 0.0);
549        assert_eq!(curvature, 0.0);
550    }
551
552    #[test]
553    fn test_covariance_error() {
554        let points: Vec<SimplePoint> = vec![
555            SimplePoint {
556                position: Vector3::new(f64::NAN, 0.0, 0.0),
557                intensity: 42,
558            },
559            SimplePoint {
560                position: Vector3::new(0.0, 1.0, f64::NAN),
561                intensity: 84,
562            },
563            SimplePoint {
564                position: Vector3::new(1.0, 1.0, f64::NAN),
565                intensity: 84,
566            },
567            SimplePoint {
568                position: Vector3::new(-1.0, f64::NAN, 0.0),
569                intensity: 84,
570            },
571        ];
572
573        let interleaved = points.into_iter().collect::<VectorBuffer>();
574
575        let result = compute_covariance_matrix(&interleaved);
576        let expected_result = Err("The number of valid (finite and non-NaN values) points in a k nearest neighborhood is not enough to span a plane!");
577        assert_eq!(result, expected_result);
578    }
579
580    #[test]
581    fn test_compute_normal() {
582        let points: Vec<SimplePoint> = vec![
583            SimplePoint {
584                position: Vector3::new(1.0, 0.0, 0.0),
585                intensity: 42,
586            },
587            SimplePoint {
588                position: Vector3::new(0.0, 1.0, 0.0),
589                intensity: 84,
590            },
591            SimplePoint {
592                position: Vector3::new(1.0, 1.0, 0.0),
593                intensity: 84,
594            },
595            SimplePoint {
596                position: Vector3::new(-1.0, 0.0, 0.0),
597                intensity: 84,
598            },
599        ];
600
601        let interleaved = points.into_iter().collect::<VectorBuffer>();
602
603        let solution_vec = compute_normals::<VectorBuffer, SimplePoint>(&interleaved, 3);
604        for solution in solution_vec {
605            assert_eq!(solution.0[0], 0.0);
606            assert_eq!(solution.0[1], 0.0);
607            assert_ne!(solution.0[2], 0.0);
608            assert_eq!(solution.1, 0.0);
609        }
610    }
611
612    #[test]
613    #[should_panic(
614        expected = "The point cloud is too small. Please use a point cloud that has 3 or more points!"
615    )]
616    fn test_compute_normal_not_enough_points_1() {
617        let points: Vec<SimplePoint> = vec![SimplePoint {
618            position: Vector3::new(1.0, 0.0, 0.0),
619            intensity: 42,
620        }];
621
622        let interleaved = points.into_iter().collect::<VectorBuffer>();
623
624        let _solution_vec = compute_normals::<VectorBuffer, SimplePoint>(&interleaved, 3);
625    }
626    #[test]
627    #[should_panic(
628        expected = "The point cloud is too small. Please use a point cloud that has 3 or more points!"
629    )]
630    fn test_compute_normal_not_enough_points_2() {
631        let points: Vec<SimplePoint> = vec![
632            SimplePoint {
633                position: Vector3::new(1.0, 0.0, 0.0),
634                intensity: 42,
635            },
636            SimplePoint {
637                position: Vector3::new(0.0, 1.0, 0.0),
638                intensity: 84,
639            },
640        ];
641
642        let interleaved = points.into_iter().collect::<VectorBuffer>();
643
644        let _solution_vec = compute_normals::<VectorBuffer, SimplePoint>(&interleaved, 3);
645    }
646
647    #[test]
648    #[should_panic(expected = "The k nearest neigbors attribute is too small!")]
649    fn test_compute_normal_knn_too_small_1() {
650        let points: Vec<SimplePoint> = vec![
651            SimplePoint {
652                position: Vector3::new(1.0, 0.0, 0.0),
653                intensity: 42,
654            },
655            SimplePoint {
656                position: Vector3::new(0.0, 1.0, 0.0),
657                intensity: 84,
658            },
659            SimplePoint {
660                position: Vector3::new(1.0, 1.0, 0.0),
661                intensity: 84,
662            },
663            SimplePoint {
664                position: Vector3::new(-1.0, 0.0, 0.0),
665                intensity: 84,
666            },
667        ];
668
669        let interleaved = points.into_iter().collect::<VectorBuffer>();
670
671        let _solution_vec = compute_normals::<VectorBuffer, SimplePoint>(&interleaved, 1);
672    }
673    #[test]
674    #[should_panic(expected = "The k nearest neigbors attribute is too small!")]
675    fn test_compute_normal_knn_too_small_2() {
676        let points: Vec<SimplePoint> = vec![
677            SimplePoint {
678                position: Vector3::new(1.0, 0.0, 0.0),
679                intensity: 42,
680            },
681            SimplePoint {
682                position: Vector3::new(0.0, 1.0, 0.0),
683                intensity: 84,
684            },
685            SimplePoint {
686                position: Vector3::new(1.0, 1.0, 0.0),
687                intensity: 84,
688            },
689            SimplePoint {
690                position: Vector3::new(-1.0, 0.0, 0.0),
691                intensity: 84,
692            },
693        ];
694
695        let interleaved = points.into_iter().collect::<VectorBuffer>();
696
697        let _solution_vec = compute_normals::<VectorBuffer, SimplePoint>(&interleaved, 2);
698    }
699}
pasture_algorithms/normal_estimation.rs

pasture_algorithms/
normal_estimation.rs
