theon 0.2.0

Abstraction of Euclidean spaces.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

//! LAPACK and non-trivial linear algebra.

#![cfg(all(feature = "lapack", target_os = "linux"))]

use ndarray::{Array, Ix2};
use ndarray_linalg::convert;
use ndarray_linalg::layout::MatrixLayout;
use ndarray_linalg::svd::SVDInto;
use typenum::type_operators::Cmp;
use typenum::{Greater, Unsigned, U2};

use crate::adjunct::{FromItems, IntoItems};
use crate::query::{Plane, Unit};
use crate::space::{EuclideanSpace, FiniteDimensional, Scalar, Vector};

/// Scalar types that can be used with LAPACK.
pub trait Lapack: ndarray_linalg::types::Lapack + ndarray_linalg::types::Scalar {}

impl<T> Lapack for T where T: ndarray_linalg::types::Lapack + ndarray_linalg::types::Scalar {}

impl<S> Plane<S>
where
    S: EuclideanSpace + FiniteDimensional,
    <S as FiniteDimensional>::N: Cmp<U2, Output = Greater>,
{
    pub fn from_points<I>(points: I) -> Option<Self>
    where
        Scalar<S>: Lapack,
        Vector<S>: FromItems + IntoItems,
        I: AsRef<[S]>,
    {
        svd_ev_plane(points)
    }
}

/// Maps columnar data into a two-dimensional array.
///
/// Produces a two-dimensional array that forms a matrix from each input column.
fn map_into_array<I, T, U, F>(columns: I, f: F) -> Option<Array<U::Item, Ix2>>
where
    I: AsRef<[T]>,
    U: FiniteDimensional + IntoItems,
    F: Fn(&T) -> U,
{
    let columns = columns.as_ref();
    let n = columns.len();
    convert::into_matrix(
        MatrixLayout::F {
            col: n as i32,
            lda: <U as FiniteDimensional>::N::USIZE as i32,
        },
        columns
            .iter()
            .map(f)
            .flat_map(|column| column.into_items())
            .collect(),
    )
    .ok()
}

// TODO: Handle edge cases and improve error handling.
/// Computes a best-fit plane from a set of points.
///
/// The plane is fit using least squares via a singular value decomposition.
fn svd_ev_plane<S, I>(points: I) -> Option<Plane<S>>
where
    S: EuclideanSpace + FiniteDimensional,
    <S as FiniteDimensional>::N: Cmp<U2, Output = Greater>,
    Scalar<S>: Lapack,
    Vector<S>: FromItems + IntoItems,
    I: AsRef<[S]>,
{
    let points = points.as_ref();
    let centroid = EuclideanSpace::centroid(points.iter().cloned())?;
    let m = map_into_array(points, |point| *point - centroid)?;
    // TODO: Fails at runtime if `V^T` is not requested.
    if let Ok((Some(u), sigma, _)) = m.svd_into(true, true) {
        let i = sigma
            .iter()
            .enumerate()
            .min_by(|(_, v1), (_, v2)| v1.partial_cmp(v2).unwrap())?
            .0;
        if i < u.ncols() {
            let normal = Vector::<S>::from_items(u.column(i).into_iter().cloned())?;
            Some(Plane {
                origin: centroid,
                normal: Unit::try_from_inner(normal)?,
            })
        }
        else {
            None
        }
    }
    else {
        None
    }
}

#[cfg(test)]
mod tests {
    use approx::assert_abs_diff_eq;
    use nalgebra::Point3;

    use crate::query::Plane;
    use crate::space::{EuclideanSpace, Vector};

    type E3 = Point3<f64>;

    #[test]
    fn determined_svd_ev_plane_e3() {
        let plane = Plane::<E3>::from_points(vec![
            EuclideanSpace::from_xyz(1.0, 0.0, 0.0),
            EuclideanSpace::from_xyz(0.0, 1.0, 0.0),
            EuclideanSpace::from_xyz(0.0, 0.0, 0.0),
        ])
        .unwrap();
        assert_abs_diff_eq!(
            Vector::<E3>::z(),
            // The axial direction of the normal is ambiguous.
            plane.normal.get().map(f64::abs),
        );
    }

    #[test]
    fn overdetermined_svd_ev_plane_e3() {
        let plane = Plane::<E3>::from_points(vec![
            EuclideanSpace::from_xyz(1.0, 1.0, 0.0),
            EuclideanSpace::from_xyz(2.0, 1.0, 0.0),
            EuclideanSpace::from_xyz(3.0, 1.0, 0.0),
            EuclideanSpace::from_xyz(2.0, 1.0, 0.0),
            EuclideanSpace::from_xyz(2.0, 2.0, 0.0),
            EuclideanSpace::from_xyz(2.0, 3.0, 0.0),
        ])
        .unwrap();
        assert_abs_diff_eq!(
            Vector::<E3>::z(),
            // The axial direction of the normal is ambiguous.
            plane.normal.get().map(f64::abs),
        );
    }
}