checkerboard_calibrate/calibrate/
extrinsics.rs1use nalgebra::{Matrix3, Rotation3, Vector3};
14
15#[derive(Clone, Copy, Debug)]
17pub struct Extrinsics {
18 pub rotation: Rotation3<f64>,
19 pub translation: Vector3<f64>,
20}
21
22pub fn init_extrinsics(k: &Matrix3<f64>, h: &Matrix3<f64>) -> Option<Extrinsics> {
26 let kinv = k.try_inverse()?;
27
28 let kh1 = kinv * h.column(0);
29 let kh2 = kinv * h.column(1);
30 let kh3 = kinv * h.column(2);
31
32 let n1 = kh1.norm();
33 if n1 < f64::EPSILON {
34 return None;
35 }
36 let lambda = 2.0 / (n1 + kh2.norm());
39
40 let mut r1 = lambda * kh1;
41 let mut r2 = lambda * kh2;
42 let mut t = lambda * kh3;
43
44 if t.z < 0.0 {
46 r1 = -r1;
47 r2 = -r2;
48 t = -t;
49 }
50 let r3 = r1.cross(&r2);
51
52 let mut q = Matrix3::zeros();
54 q.set_column(0, &r1);
55 q.set_column(1, &r2);
56 q.set_column(2, &r3);
57
58 let svd = q.svd(true, true);
59 let u = svd.u?;
60 let v_t = svd.v_t?;
61 let mut rot = u * v_t;
62 if rot.determinant() < 0.0 {
63 let mut u2 = u;
65 let last = u2.ncols() - 1;
66 let col = -u2.column(last);
67 u2.set_column(last, &col);
68 rot = u2 * v_t;
69 }
70
71 let rotation = Rotation3::from_matrix_unchecked(rot);
72 Some(Extrinsics {
73 rotation,
74 translation: t,
75 })
76}
77
78#[cfg(test)]
79mod tests {
80 use super::*;
81
82 fn homography_for_view(k: &Matrix3<f64>, r: &Rotation3<f64>, t: &Vector3<f64>) -> Matrix3<f64> {
83 let rm = r.matrix();
84 let mut m = Matrix3::zeros();
85 m.set_column(0, &rm.column(0));
86 m.set_column(1, &rm.column(1));
87 m.set_column(2, t);
88 k * m
89 }
90
91 #[test]
92 fn recovers_pose() {
93 let k = Matrix3::new(520.0, 0.0, 319.5, 0.0, 510.0, 239.5, 0.0, 0.0, 1.0);
94 let r_true = Rotation3::from_euler_angles(0.12, -0.22, 0.07);
95 let t_true = Vector3::new(-1.3, 0.6, 8.0);
96
97 let h = homography_for_view(&k, &r_true, &t_true);
98 let ext = init_extrinsics(&k, &h).expect("extrinsics");
99
100 approx::assert_abs_diff_eq!(ext.translation, t_true, epsilon = 1e-9);
101 approx::assert_abs_diff_eq!(ext.rotation.matrix(), r_true.matrix(), epsilon = 1e-9);
102 }
103}