vecmat/transform/
linear.rs

1#[cfg(feature = "rand")]
2use crate::distr::{Invertible, Normal};
3use crate::{
4    traits::{Dot, Normalize},
5    transform::{Reorder, Shift, Directional},
6    Matrix, Transform, Vector,
7};
8#[cfg(feature = "approx")]
9use approx::{abs_diff_eq, AbsDiffEq};
10use core::ops::Neg;
11use num_traits::{Float, Num, NumCast, One, Inv};
12#[cfg(feature = "rand")]
13use rand_::{distributions::Distribution, Rng};
14
15/// Linear transformation.
16#[repr(transparent)]
17#[derive(Clone, Copy, PartialEq, Debug)]
18pub struct Linear<T, const N: usize> {
19    lin: Matrix<T, N, N>,
20}
21
22pub type Linear2<T> = Linear<T, 2>;
23pub type Linear3<T> = Linear<T, 3>;
24pub type Linear4<T> = Linear<T, 4>;
25
26impl<T, const N: usize> Linear<T, N> {
27    pub fn from_matrix(lin: Matrix<T, N, N>) -> Self {
28        Self { lin }
29    }
30    pub fn into_matrix(self) -> Matrix<T, N, N> {
31        self.lin
32    }
33}
34impl<T, const N: usize> From<Matrix<T, N, N>> for Linear<T, N> {
35    fn from(lin: Matrix<T, N, N>) -> Self {
36        Self::from_matrix(lin)
37    }
38}
39impl<T, const N: usize> From<Linear<T, N>> for Matrix<T, N, N> {
40    fn from(lin: Linear<T, N>) -> Matrix<T, N, N> {
41        lin.into_matrix()
42    }
43}
44
45impl<T, const N: usize> Transform<Vector<T, N>> for Linear<T, N>
46where
47    T: Neg<Output = T> + Num + Copy,
48{
49    fn identity() -> Self {
50        Self { lin: Matrix::one() }
51    }
52    fn inv(self) -> Self {
53        Self {
54            lin: self.lin.inv(),
55        }
56    }
57    fn apply(&self, pos: Vector<T, N>) -> Vector<T, N> {
58        self.lin.dot(pos)
59    }
60    fn deriv(&self, _pos: Vector<T, N>, dir: Vector<T, N>) -> Vector<T, N> {
61        self.apply(dir)
62    }
63    fn chain(self, other: Self) -> Self {
64        Self {
65            lin: self.lin.dot(other.lin),
66        }
67    }
68}
69
70impl<T, const N: usize> Linear<T, N>
71where
72    T: Neg<Output = T> + Num + Copy,
73    Matrix<T, N, N>: Inv<Output=Matrix<T, N, N>>,
74{
75    pub fn normal_transform(self) -> Self {
76        Self { lin: self.lin.inv().transpose() }
77    }
78}
79
80impl<T, const N: usize> Directional<Vector<T, N>> for Linear<T, N>
81where
82    Self: Transform<Vector<T, N>>,
83    T: Neg<Output = T> + Num + Copy,
84    Vector<T, N>: Normalize,
85    Matrix<T, N, N>: Inv<Output=Matrix<T, N, N>>,
86{
87    fn apply_dir(&self, pos: Vector<T, N>, dir: Vector<T, N>) -> Vector<T, N> {
88        self.deriv(pos, dir).normalize()
89    }
90    fn apply_normal(&self, _: Vector<T, N>, normal: Vector<T, N>) -> Vector<T, N> {
91        self.normal_transform().apply(normal).normalize()
92    }
93}
94
95#[cfg(feature = "rand")]
96impl<T, const N: usize> Distribution<Linear<T, N>> for Normal
97where
98    Normal: Distribution<Matrix<T, N, N>>,
99{
100    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Linear<T, N> {
101        Linear::from_matrix(self.sample(rng))
102    }
103}
104#[cfg(feature = "rand")]
105impl<T, const N: usize> Distribution<Linear<T, N>> for Invertible
106where
107    Invertible: Distribution<Matrix<T, N, N>>,
108{
109    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Linear<T, N> {
110        Linear::from_matrix(self.sample(rng))
111    }
112}
113
114#[cfg(feature = "approx")]
115impl<T, const N: usize> AbsDiffEq for Linear<T, N>
116where
117    T: AbsDiffEq<Epsilon = T> + Copy,
118{
119    type Epsilon = T;
120    fn default_epsilon() -> Self::Epsilon {
121        T::default_epsilon()
122    }
123    fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
124        abs_diff_eq!(self.lin, other.lin, epsilon = epsilon)
125    }
126}
127
128impl<T> Linear<T, 3>
129where
130    T: Float,
131{
132    /// Returns the transformation that rotates `-z`-axis to `dir`
133    /// and `y`-axis to `up`.
134    pub fn look_at(dir: Vector<T, 3>, up: Vector<T, 3>) -> Self {
135        let right = dir.cross(up).normalize();
136        let strict_up = right.cross(dir);
137        Self::from(Matrix::from([right, strict_up, -dir]).transpose())
138    }
139}
140impl<T> Linear<T, 3>
141where
142    T: Float + NumCast,
143{
144    /// Returns any of transformations that rotate `dir` to `-z`-axis.
145    pub fn look_at_any(dir: Vector<T, 3>) -> Self {
146        if dir.z().abs() < T::from(0.5).unwrap() {
147            Self::look_at(dir, Vector::from([T::zero(), T::zero(), T::one()]))
148        } else {
149            Self::look_at(dir, Vector::from([T::zero(), T::one(), T::zero()]))
150        }
151    }
152}
153
154impl<T, const N: usize> Reorder<Linear<T, N>, Vector<T, N>> for Shift<T, N>
155where
156    Linear<T, N>: Transform<Vector<T, N>> + Copy,
157    Self: Transform<Vector<T, N>>,
158{
159    fn reorder(self, other: Linear<T, N>) -> (Linear<T, N>, Shift<T, N>) {
160        (other, other.inv().apply(self.into_vector()).into())
161    }
162}
163impl<T, const N: usize> Reorder<Shift<T, N>, Vector<T, N>> for Linear<T, N>
164where
165    Self: Transform<Vector<T, N>>,
166    Shift<T, N>: Transform<Vector<T, N>>,
167{
168    fn reorder(self, other: Shift<T, N>) -> (Shift<T, N>, Linear<T, N>) {
169        (self.apply(other.into_vector()).into(), self)
170    }
171}
172
173#[cfg(all(test, feature = "rand", feature = "approx"))]
174mod tests {
175    use super::*;
176    use crate::distr::{Normal, Unit};
177    use approx::assert_abs_diff_eq;
178    use rand_::prelude::*;
179    use rand_xorshift::XorShiftRng;
180
181    const SAMPLE_ATTEMPTS: usize = 256;
182
183    #[test]
184    fn linearity() {
185        const EPS: f64 = 1e-14;
186        let mut rng = XorShiftRng::seed_from_u64(0xBEE);
187
188        for _ in 0..SAMPLE_ATTEMPTS {
189            let m: Matrix<f64, 3, 3> = rng.sample(&Normal);
190            let x: Vector<f64, 3> = rng.sample(&Normal);
191            let a: f64 = rng.sample(&Normal);
192
193            assert_abs_diff_eq!(
194                Linear::from(m * a).apply(x),
195                Linear::from(m).apply(x * a),
196                epsilon = EPS
197            );
198            assert_abs_diff_eq!(
199                Linear::from(m * a).apply(x),
200                Linear::from(m).apply(x) * a,
201                epsilon = EPS
202            );
203        }
204    }
205
206    #[test]
207    fn chaining() {
208        const EPS: f64 = 1e-14;
209        let mut rng = XorShiftRng::seed_from_u64(0xBEA);
210
211        for _ in 0..SAMPLE_ATTEMPTS {
212            let a: Linear<f64, 3> = rng.sample(&Normal);
213            let b: Linear<f64, 3> = rng.sample(&Normal);
214            let c: Vector<f64, 3> = rng.sample(&Normal);
215
216            assert_abs_diff_eq!(a.chain(Linear::identity()), a, epsilon = EPS);
217            assert_abs_diff_eq!(Linear::identity().chain(b), b, epsilon = EPS);
218            assert_abs_diff_eq!(a.chain(b).apply(c), a.apply(b.apply(c)), epsilon = EPS);
219        }
220    }
221
222    #[test]
223    fn inversion() {
224        const EPS: f64 = 1e-12;
225        let mut rng = XorShiftRng::seed_from_u64(0xBEB);
226
227        for _ in 0..SAMPLE_ATTEMPTS {
228            let a: Linear<f64, 3> = rng.sample(&Invertible);
229            let x: Vector<f64, 3> = rng.sample(&Normal);
230
231            assert_abs_diff_eq!(a.chain(a.inv()), Linear::identity(), epsilon = EPS);
232            assert_abs_diff_eq!(a.inv().chain(a), Linear::identity(), epsilon = EPS);
233            assert_abs_diff_eq!(a.inv().apply(a.apply(x)), x, epsilon = EPS);
234            assert_abs_diff_eq!(a.apply(a.inv().apply(x)), x, epsilon = EPS);
235        }
236    }
237
238    #[test]
239    fn look_to_the_direction() {
240        const EPS: f64 = 1e-14;
241        let mut rng = XorShiftRng::seed_from_u64(0xBEC);
242
243        for _ in 0..SAMPLE_ATTEMPTS {
244            let d: Vector<f64, 3> = rng.sample(&Unit);
245            let m = Linear::look_at_any(d);
246
247            assert_abs_diff_eq!(m.apply(Vector::from([0.0, 0.0, -1.0])), d, epsilon = EPS);
248        }
249    }
250}
vecmat/transform/linear.rs

vecmat/transform/
linear.rs
