pbrt_r3/core/transform/
decompose.rs

1use super::matrix4x4::*;
2use crate::core::base::*;
3use crate::core::quaternion::*;
4
5fn invertible(m: &Matrix4x4) -> bool {
6    m.inverse().is_some()
7}
8
9fn length(x: Float, y: Float, z: Float) -> Float {
10    (x * x + y * y + z * z).sqrt()
11}
12
13fn suppress_for_scale(m: Matrix4x4) -> Matrix4x4 {
14    let mut mm = Matrix4x4::identity();
15    for i in 0..3 {
16        mm.m[4 * i + i] = m.m[4 * i + i];
17    }
18    return mm;
19}
20
21// Decompose a 4x4 matrix into translation, rotation, and scale.
22pub fn decompose(
23    m: &Matrix4x4,
24    epsilon: Float,
25    max_count: i32,
26) -> Option<(Vector3f, Quaternion, Vector3f)> {
27    // Extract translation _T_ from transformation matrix
28    let t = Vector3f::new(m.m[4 * 0 + 3], m.m[4 * 1 + 3], m.m[4 * 2 + 3]);
29
30    // Compute new transformation matrix _M_ without translation
31    let mut mm = *m;
32    for i in 0..3 {
33        mm.m[4 * i + 3] = 0.0;
34    }
35    mm.m[15] = 1.0;
36
37    // Extract rotation _R_ from transformation matrix
38    let mut r = mm;
39    // pbrt-r3
40    let sx = length(r.m[4 * 0 + 0], r.m[4 * 1 + 0], r.m[4 * 2 + 0]);
41    let sy = length(r.m[4 * 0 + 1], r.m[4 * 1 + 1], r.m[4 * 2 + 1]);
42    let sz = length(r.m[4 * 0 + 2], r.m[4 * 1 + 2], r.m[4 * 2 + 2]);
43    if sx != 0.0 {
44        r.m[4 * 0 + 0] /= sx;
45        r.m[4 * 0 + 1] /= sx;
46        r.m[4 * 0 + 2] /= sx;
47    }
48    if sy != 0.0 {
49        r.m[4 * 1 + 0] /= sy;
50        r.m[4 * 1 + 1] /= sy;
51        r.m[4 * 1 + 2] /= sy;
52    }
53    if sz != 0.0 {
54        r.m[4 * 2 + 0] /= sz;
55        r.m[4 * 2 + 1] /= sz;
56        r.m[4 * 2 + 2] /= sz;
57    }
58    // pbrt-r3
59    let mut count = 0;
60    let mut norm: Float = 0.0;
61    loop {
62        // Compute inverse of _R_ and check for singularity
63        assert!(invertible(&r));
64        if let Some(r_it) = r.transpose().inverse() {
65            // Compute next matrix _Rnext_ in series
66            let mut r_next = r;
67            for i in 0..4 {
68                for j in 0..4 {
69                    r_next.m[4 * i + j] = lerp(0.5, r.m[4 * i + j], r_it.m[4 * i + j]);
70                }
71            }
72
73            // pbrt-r3
74            assert!(invertible(&r_next));
75            if let Some(ir_next) = r_next.inverse() {
76                let s = suppress_for_scale(ir_next * mm);
77                if let Some(is) = s.inverse() {
78                    r_next = mm * is;
79                }
80            }
81            let q = Quaternion::from(r_next);
82            let q = q.normalize();
83            r_next = q.to_matrix();
84            // pbrt-r3
85
86            // Compute norm of difference between _R_ and _Rnext_
87            for i in 0..3 {
88                let n = Float::abs(r.m[4 * i + 0] - r_next.m[4 * i + 0])
89                    + Float::abs(r.m[4 * i + 1] - r_next.m[4 * i + 1])
90                    + Float::abs(r.m[4 * i + 2] - r_next.m[4 * i + 2]);
91                norm = Float::max(norm, n);
92            }
93            r = r_next;
94        } else {
95            return None;
96        }
97
98        count += 1;
99        if !(count < max_count && norm > epsilon) {
100            break;
101        }
102    }
103
104    if let Some(ir) = r.inverse() {
105        let s = suppress_for_scale(ir * mm);
106        if let Some(is) = s.inverse() {
107            r = mm * is;
108        }
109        let q = Quaternion::from(r);
110        let q = q.normalize(); //pbrt-r3
111        let s = Vector3f::new(s.m[0], s.m[5], s.m[10]);
112        return Some((t, q, s));
113    }
114    return None;
115}
pbrt_r3/core/transform/decompose.rs

pbrt_r3/core/transform/
decompose.rs
