Skip to content

Commit

Permalink
Use F64 compute scale for const matrix inversion
Browse files Browse the repository at this point in the history
The tree jump test is now it's original margin
  • Loading branch information
Beinsezii committed Oct 25, 2024
1 parent a3adc2d commit 1a99f1e
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 11 deletions.
32 changes: 22 additions & 10 deletions src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -255,35 +255,47 @@ const JZAZBZ_P: f32 = 1.7 * PQEOTF_M2;

// ### MATRICES ### {{{

/// Const matrix upcast.
/// What I wouldn't give for .map() to work
const fn m32_to_m64(m: [[f32; 3]; 3]) -> [[f64; 3]; 3] {
[
[m[0][0] as f64, m[0][1] as f64, m[0][2] as f64],
[m[1][0] as f64, m[1][1] as f64, m[1][2] as f64],
[m[2][0] as f64, m[2][1] as f64, m[2][2] as f64],
]
}

/// Matrix determinant using only constant math
const fn det(m: [[f32; 3]; 3]) -> f32 {
const fn det(m: [[f64; 3]; 3]) -> f64 {
m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
+ m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])
}

/// Matrix inversion using only constant math
/// Panics if determinant is zero
/// Compute scale of f64
const fn inv(m: [[f32; 3]; 3]) -> [[f32; 3]; 3] {
let m = m32_to_m64(m);
let d = det(m);
if d == 0.0 {
panic!("Matrix is singular and has no inverse")
}

[
[
(m[1][1] * m[2][2] - m[1][2] * m[2][1]) / d,
(m[0][2] * m[2][1] - m[0][1] * m[2][2]) / d,
(m[0][1] * m[1][2] - m[0][2] * m[1][1]) / d,
((m[1][1] * m[2][2] - m[1][2] * m[2][1]) / d) as f32,
((m[0][2] * m[2][1] - m[0][1] * m[2][2]) / d) as f32,
((m[0][1] * m[1][2] - m[0][2] * m[1][1]) / d) as f32,
],
[
(m[1][2] * m[2][0] - m[1][0] * m[2][2]) / d,
(m[0][0] * m[2][2] - m[0][2] * m[2][0]) / d,
(m[0][2] * m[1][0] - m[0][0] * m[1][2]) / d,
((m[1][2] * m[2][0] - m[1][0] * m[2][2]) / d) as f32,
((m[0][0] * m[2][2] - m[0][2] * m[2][0]) / d) as f32,
((m[0][2] * m[1][0] - m[0][0] * m[1][2]) / d) as f32,
],
[
(m[1][0] * m[2][1] - m[1][1] * m[2][0]) / d,
(m[0][1] * m[2][0] - m[0][0] * m[2][1]) / d,
(m[0][0] * m[1][1] - m[0][1] * m[1][0]) / d,
((m[1][0] * m[2][1] - m[1][1] * m[2][0]) / d) as f32,
((m[0][1] * m[2][0] - m[0][0] * m[2][1]) / d) as f32,
((m[0][0] * m[1][1] - m[0][1] * m[1][0]) / d) as f32,
],
]
}
Expand Down
2 changes: 1 addition & 1 deletion src/tests.rs
Original file line number Diff line number Diff line change
Expand Up @@ -298,7 +298,7 @@ fn tree_jump() {
convert_space_chunked::<f64, 3>($from_space, $to_space, &mut input);
// strange this is 1e-3 while indiv is 1e-2
// also skip places where hue can wrap
pix_cmp(&input, $to_data, 5e-3, &[0, 1, 7])
pix_cmp(&input, $to_data, 1e-3, &[0, 1, 7])
};
}

Expand Down

0 comments on commit 1a99f1e

Please sign in to comment.