measures/inner/exact/
linear_map_3d.rs

1#[macro_export] // Don't add nor remove the first three lines and the last two lines.
2macro_rules! inner_define_linear_map_3d {
3    {} => {
4        /// Linear transformation of `Measure3d` objects.
5        pub struct LinearMap3d<Number = f64>
6        where
7            Number: ArithmeticOps,
8        {
9            pub c: [[Number; 3]; 3],
10        }
11
12        impl<Number> LinearMap3d<Number>
13        where
14            Number: ArithmeticOps,
15        {
16            /// Create a LinearMap3d from its 9 coefficients.
17            pub const fn new(coefficients: [[Number; 3]; 3]) -> Self {
18                Self { c: coefficients }
19            }
20
21            // Linear maps have no translations.
22
23            // Rotations
24
25            // Rotation by an angle measure around a unit vector.
26            // Precondition: unit_vector.squared_norm().value == 1
27            pub fn rotation<AngleUnit, AxisUnit>(
28                angle: Measure<AngleUnit, Number>,
29                unit_vector: Measure3d<AxisUnit, Number>,
30            ) -> Self
31            where
32                AngleUnit: AngleMeasurementUnit,
33                AxisUnit: MeasurementUnit<Property: VectorProperty>,
34            {
35                Self::rotation_by_radians_around_unit_vector(
36                    angle.convert::<Radian>().value,
37                    unit_vector.values[0],
38                    unit_vector.values[1],
39                    unit_vector.values[2],
40                )
41            }
42
43            // Projections
44
45            // Projection onto a line identified by a unit vector.
46            // Precondition: unit_vector.squared_norm().value == 1
47            pub fn projection_onto_line<Unit>(unit_vector: Measure3d<Unit, Number>) -> Self
48            where
49                Unit: MeasurementUnit<Property: VectorProperty>,
50            {
51                Self {
52                    c: [
53                        [
54                            unit_vector.values[0] * unit_vector.values[0],
55                            unit_vector.values[1] * unit_vector.values[0],
56                            unit_vector.values[2] * unit_vector.values[0],
57                        ],
58                        [
59                            unit_vector.values[0] * unit_vector.values[1],
60                            unit_vector.values[1] * unit_vector.values[1],
61                            unit_vector.values[2] * unit_vector.values[1],
62                        ],
63                        [
64                            unit_vector.values[0] * unit_vector.values[2],
65                            unit_vector.values[1] * unit_vector.values[2],
66                            unit_vector.values[2] * unit_vector.values[2],
67                        ],
68                    ],
69                }
70            }
71
72            // Projection onto a plane whose normal is identified by a unit vector.
73            // Precondition: unit_vector.squared_norm().value == 1
74            pub fn projection_onto_plane<Unit>(unit_vector: Measure3d<Unit, Number>) -> Self
75            where
76                Unit: MeasurementUnit<Property: VectorProperty>,
77            {
78                Self {
79                    c: [
80                        [
81                            Number::ONE - unit_vector.values[0] * unit_vector.values[0],
82                            -unit_vector.values[1] * unit_vector.values[0],
83                            -unit_vector.values[2] * unit_vector.values[0],
84                        ],
85                        [
86                            -unit_vector.values[0] * unit_vector.values[1],
87                            Number::ONE - unit_vector.values[1] * unit_vector.values[1],
88                            -unit_vector.values[2] * unit_vector.values[1],
89                        ],
90                        [
91                            -unit_vector.values[0] * unit_vector.values[2],
92                            -unit_vector.values[1] * unit_vector.values[2],
93                            Number::ONE - unit_vector.values[2] * unit_vector.values[2],
94                        ],
95                    ],
96                }
97            }
98
99            // Reflections
100
101            // Reflection over a line identified by a unit vector.
102            // Precondition: unit_vector.squared_norm().value == 1
103            pub fn reflection_over_line<AxisUnit>(unit_vector: Measure3d<AxisUnit, Number>) -> Self
104            where
105                AxisUnit: MeasurementUnit<Property: VectorProperty>,
106            {
107                let two = Number::ONE + Number::ONE;
108                Self {
109                    c: [
110                        [
111                            two * unit_vector.values[0] * unit_vector.values[0] - Number::ONE,
112                            two * unit_vector.values[1] * unit_vector.values[0],
113                            two * unit_vector.values[2] * unit_vector.values[0],
114                        ],
115                        [
116                            two * unit_vector.values[0] * unit_vector.values[1],
117                            two * unit_vector.values[1] * unit_vector.values[1] - Number::ONE,
118                            two * unit_vector.values[2] * unit_vector.values[1],
119                        ],
120                        [
121                            two * unit_vector.values[0] * unit_vector.values[2],
122                            two * unit_vector.values[1] * unit_vector.values[2],
123                            two * unit_vector.values[2] * unit_vector.values[2] - Number::ONE,
124                        ],
125                    ],
126                }
127            }
128
129            // Reflection over a plane whose normal is identified by a unit vector.
130            // Precondition: unit_vector.squared_norm().value == 1
131            pub fn reflection_over_plane<Unit>(unit_vector: Measure3d<Unit, Number>) -> Self
132            where
133                Unit: MeasurementUnit<Property: VectorProperty>,
134            {
135                let minus_two = -(Number::ONE + Number::ONE);
136                Self {
137                    c: [
138                        [
139                            minus_two * unit_vector.values[0] * unit_vector.values[0] + Number::ONE,
140                            minus_two * unit_vector.values[1] * unit_vector.values[0],
141                            minus_two * unit_vector.values[2] * unit_vector.values[0],
142                        ],
143                        [
144                            minus_two * unit_vector.values[0] * unit_vector.values[1],
145                            minus_two * unit_vector.values[1] * unit_vector.values[1] + Number::ONE,
146                            minus_two * unit_vector.values[2] * unit_vector.values[1],
147                        ],
148                        [
149                            minus_two * unit_vector.values[0] * unit_vector.values[2],
150                            minus_two * unit_vector.values[1] * unit_vector.values[2],
151                            minus_two * unit_vector.values[2] * unit_vector.values[2] + Number::ONE,
152                        ],
153                    ],
154                }
155            }
156
157            // Scaling by three factors.
158
159            pub fn scaling(factors: [Number; 3]) -> Self {
160                Self {
161                    c: [
162                        [factors[0], Number::ZERO, Number::ZERO],
163                        [Number::ZERO, factors[1], Number::ZERO],
164                        [Number::ZERO, Number::ZERO, factors[2]],
165                    ],
166                }
167            }
168
169            // Inversion
170
171            pub fn inverted(&self) -> Self {
172                let inv_determinant = Number::ONE
173                    / (self.c[0][0] * (self.c[1][1] * self.c[2][2] - self.c[1][2] * self.c[2][1])
174                        - self.c[0][1] * (self.c[1][0] * self.c[2][2] - self.c[1][2] * self.c[2][0])
175                        + self.c[0][2] * (self.c[1][0] * self.c[2][1] - self.c[1][1] * self.c[2][0]));
176                Self {
177                    c: [
178                        [
179                            (self.c[1][1] * self.c[2][2] - self.c[1][2] * self.c[2][1]) * inv_determinant,
180                            -(self.c[0][1] * self.c[2][2] - self.c[0][2] * self.c[2][1]) * inv_determinant,
181                            (self.c[0][1] * self.c[1][2] - self.c[0][2] * self.c[1][1]) * inv_determinant,
182                        ],
183                        [
184                            -(self.c[1][0] * self.c[2][2] - self.c[1][2] * self.c[2][0]) * inv_determinant,
185                            (self.c[0][0] * self.c[2][2] - self.c[0][2] * self.c[2][0]) * inv_determinant,
186                            -(self.c[0][0] * self.c[1][2] - self.c[0][2] * self.c[1][0]) * inv_determinant,
187                        ],
188                        [
189                            (self.c[1][0] * self.c[2][1] - self.c[1][1] * self.c[2][0]) * inv_determinant,
190                            -(self.c[0][0] * self.c[2][1] - self.c[0][1] * self.c[2][0]) * inv_determinant,
191                            (self.c[0][0] * self.c[1][1] - self.c[0][1] * self.c[1][0]) * inv_determinant,
192                        ],
193                    ],
194                }
195            }
196
197            // Composition of spacial linear transformations.
198            // Applying the resulting transformation is equivalent to apply first
199            // `other` and then `self`.
200            pub fn combined_with(&self, other: &LinearMap3d<Number>) -> Self {
201                Self {
202                    c: [
203                        [
204                            other.c[0][0] * self.c[0][0]
205                                + other.c[0][1] * self.c[1][0]
206                                + other.c[0][2] * self.c[2][0],
207                            other.c[0][0] * self.c[0][1]
208                                + other.c[0][1] * self.c[1][1]
209                                + other.c[0][2] * self.c[2][1],
210                            other.c[0][0] * self.c[0][2]
211                                + other.c[0][1] * self.c[1][2]
212                                + other.c[0][2] * self.c[2][2],
213                        ],
214                        [
215                            other.c[1][0] * self.c[0][0]
216                                + other.c[1][1] * self.c[1][0]
217                                + other.c[1][2] * self.c[2][0],
218                            other.c[1][0] * self.c[0][1]
219                                + other.c[1][1] * self.c[1][1]
220                                + other.c[1][2] * self.c[2][1],
221                            other.c[1][0] * self.c[0][2]
222                                + other.c[1][1] * self.c[1][2]
223                                + other.c[1][2] * self.c[2][2],
224                        ],
225                        [
226                            other.c[2][0] * self.c[0][0]
227                                + other.c[2][1] * self.c[1][0]
228                                + other.c[2][2] * self.c[2][0],
229                            other.c[2][0] * self.c[0][1]
230                                + other.c[2][1] * self.c[1][1]
231                                + other.c[2][2] * self.c[2][1],
232                            other.c[2][0] * self.c[0][2]
233                                + other.c[2][1] * self.c[1][2]
234                                + other.c[2][2] * self.c[2][2],
235                        ],
236                    ],
237                }
238            }
239
240            pub fn apply_to<Unit>(&self, m: Measure3d<Unit, Number>) -> Measure3d<Unit, Number>
241            where
242                Unit: MeasurementUnit<Property: VectorProperty>,
243            {
244                Measure3d::<Unit, Number>::new([
245                    self.c[0][0] * m.values[0] + self.c[0][1] * m.values[1] + self.c[0][2] * m.values[2],
246                    self.c[1][0] * m.values[0] + self.c[1][1] * m.values[1] + self.c[1][2] * m.values[2],
247                    self.c[2][0] * m.values[0] + self.c[2][1] * m.values[1] + self.c[2][2] * m.values[2],
248                ])
249            }
250
251            fn rotation_by_radians_around_unit_vector(
252                a: Number,
253                ux: Number,
254                uy: Number,
255                uz: Number,
256            ) -> Self {
257                let (sin_a, cos_a) = a.sin_cos();
258                let one_minus_cos_a = Number::ONE - cos_a;
259                Self {
260                    c: [
261                        [
262                            cos_a + ux * ux * one_minus_cos_a,
263                            ux * uy * one_minus_cos_a - uz * sin_a,
264                            ux * uz * one_minus_cos_a + uy * sin_a,
265                        ],
266                        [
267                            uy * ux * one_minus_cos_a + uz * sin_a,
268                            cos_a + uy * uy * one_minus_cos_a,
269                            uy * uz * one_minus_cos_a - ux * sin_a,
270                        ],
271                        [
272                            uz * ux * one_minus_cos_a - uy * sin_a,
273                            uz * uy * one_minus_cos_a + ux * sin_a,
274                            cos_a + uz * uz * one_minus_cos_a,
275                        ],
276                    ],
277                }
278            }
279        }
280
281        impl<Number> Default for LinearMap3d<Number>
282        where
283            Number: ArithmeticOps,
284        {
285            /// LinearMap3d::default() -> LinearMap3d
286            /// It returns the identity transformation.
287            fn default() -> Self {
288                Self::new([
289                    [Number::ONE, Number::ZERO, Number::ZERO],
290                    [Number::ZERO, Number::ONE, Number::ZERO],
291                    [Number::ZERO, Number::ZERO, Number::ONE],
292                ])
293            }
294        }
295
296        // LinearMap3d == LinearMap3d -> bool
297        impl<Number> PartialEq<LinearMap3d<Number>> for LinearMap3d<Number>
298        where
299            Number: ArithmeticOps,
300        {
301            fn eq(&self, other: &LinearMap3d<Number>) -> bool {
302                self.c == other.c
303            }
304        }
305
306        // LinearMap3d.clone() -> LinearMap3d
307        impl<Number> Clone for LinearMap3d<Number>
308        where
309            Number: ArithmeticOps,
310        {
311            fn clone(&self) -> Self {
312                Self { c: self.c.clone() }
313            }
314        }
315
316        impl From<LinearMap3d<f32>> for LinearMap3d<f64> {
317            fn from(m: LinearMap3d<f32>) -> Self {
318                Self::new([
319                    [m.c[0][0] as f64, m.c[0][1] as f64, m.c[0][2] as f64],
320                    [m.c[1][0] as f64, m.c[1][1] as f64, m.c[1][2] as f64],
321                    [m.c[2][0] as f64, m.c[2][1] as f64, m.c[2][2] as f64],
322                ])
323            }
324        }
325
326        /// format!("{}", LinearMap3d) -> String
327        /// LinearMap3d.to_string() -> String
328        impl<Number: ArithmeticOps> fmt::Display for LinearMap3d<Number> {
329            fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
330                write!(
331                    formatter,
332                    "{}",
333                    measures::matrix_utils::format_matrix::<3, 3, Number>(&self.c, "", 1)
334                )
335            }
336        }
337
338        // format!("{:?}", LinearMap3d)
339        impl<Number: ArithmeticOps> fmt::Debug for LinearMap3d<Number> {
340            fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
341                write!(
342                    formatter,
343                    "{}",
344                    measures::matrix_utils::format_matrix::<3, 3, Number>(&self.c, "", 1)
345                )
346            }
347        }
348    };
349}