dssim_core/
tolab.rs

1#![allow(non_snake_case)]
2#![allow(non_upper_case_globals)]
3
4use crate::image::ToRGB;
5use crate::image::RGBAPLU;
6use crate::image::RGBLU;
7use imgref::*;
8#[cfg(not(feature = "threads"))]
9use crate::lieon as rayon;
10use rayon::prelude::*;
11
12const D65x: f32 = 0.9505;
13const D65y: f32 = 1.0;
14const D65z: f32 = 1.089;
15
16pub type GBitmap = ImgVec<f32>;
17pub(crate) trait ToLAB {
18    fn to_lab(&self) -> (f32, f32, f32);
19}
20
21#[inline(always)]
22fn fma_matrix(r: f32, rx: f32, g: f32, gx: f32, b: f32, bx: f32) -> f32 {
23    b.mul_add(bx, g.mul_add(gx, r * rx))
24}
25
26const EPSILON: f32 = 216. / 24389.;
27const K: f32 = 24389. / (27. * 116.); // http://www.brucelindbloom.com/LContinuity.html
28
29impl ToLAB for RGBLU {
30    fn to_lab(&self) -> (f32, f32, f32) {
31        let fx = fma_matrix(self.r, 0.4124 / D65x, self.g, 0.3576 / D65x, self.b, 0.1805 / D65x);
32        let fy = fma_matrix(self.r, 0.2126 / D65y, self.g, 0.7152 / D65y, self.b, 0.0722 / D65y);
33        let fz = fma_matrix(self.r, 0.0193 / D65z, self.g, 0.1192 / D65z, self.b, 0.9505 / D65z);
34
35        let X = if fx > EPSILON { cbrt_poly(fx) - 16. / 116. } else { K * fx };
36        let Y = if fy > EPSILON { cbrt_poly(fy) - 16. / 116. } else { K * fy };
37        let Z = if fz > EPSILON { cbrt_poly(fz) - 16. / 116. } else { K * fz };
38
39        let lab = (
40            (Y * 1.05f32), // 1.05 instead of 1.16 to boost color importance without pushing colors outside of 1.0 range
41            (500.0 / 220.0f32).mul_add(X - Y, 86.2 / 220.0f32), /* 86 is a fudge to make the value positive */
42            (200.0 / 220.0f32).mul_add(Y - Z, 107.9 / 220.0f32), /* 107 is a fudge to make the value positive */
43        );
44        debug_assert!(lab.0 <= 1.0 && lab.1 <= 1.0 && lab.2 <= 1.0);
45        lab
46    }
47}
48
49#[inline]
50fn cbrt_poly(x: f32) -> f32 {
51    // Polynomial approximation
52    let poly = [0.2f32, 1.51, -0.5];
53    let y = poly[2].mul_add(x, poly[1]).mul_add(x, poly[0]);
54
55    // 2x Halley's Method
56    let y3 = y * y * y;
57    let y = y * 2.0f32.mul_add(x, y3) / 2.0f32.mul_add(y3, x);
58    let y3 = y * y * y;
59    let y = y * 2.0f32.mul_add(x, y3) / 2.0f32.mul_add(y3, x);
60    debug_assert!(y < 1.001);
61    debug_assert!(x < 216. / 24389. || y >= 16. / 116.);
62    y
63}
64
65/// Convert image to L\*a\*b\* planar
66///
67/// It should return 1 (gray) or 3 (color) planes.
68pub trait ToLABBitmap {
69    fn to_lab(&self) -> Vec<GBitmap>;
70}
71
72impl ToLABBitmap for ImgVec<RGBAPLU> {
73    #[inline(always)]
74    fn to_lab(&self) -> Vec<GBitmap> {
75        self.as_ref().to_lab()
76    }
77}
78
79impl ToLABBitmap for ImgVec<RGBLU> {
80    #[inline(always)]
81    fn to_lab(&self) -> Vec<GBitmap> {
82        self.as_ref().to_lab()
83    }
84}
85impl ToLABBitmap for GBitmap {
86    fn to_lab(&self) -> Vec<GBitmap> {
87        debug_assert!(self.width() > 0);
88        let f = |fy| {
89            if fy > EPSILON { (cbrt_poly(fy) - 16. / 116.) * 1.16 } else { (K * 1.16) * fy }
90        };
91
92        #[cfg(feature = "threads")]
93        let out = (0..self.height()).into_par_iter().flat_map_iter(|y| {
94            self[y].iter().map(|&fy| f(fy))
95        }).collect();
96
97        #[cfg(not(feature = "threads"))]
98        let out = self.pixels().map(f).collect();
99
100        vec![Self::new(out, self.width(), self.height())]
101    }
102}
103
104#[inline(never)]
105fn rgb_to_lab<T: Copy + Sync + Send + 'static, F>(img: ImgRef<'_, T>, cb: F) -> Vec<GBitmap>
106    where F: Fn(T, usize) -> (f32, f32, f32) + Sync + Send + 'static
107{
108    let width = img.width();
109    assert!(width > 0);
110    let height = img.height();
111    let area = width * height;
112
113    let mut out_l = Vec::with_capacity(area);
114    let mut out_a = Vec::with_capacity(area);
115    let mut out_b = Vec::with_capacity(area);
116
117    // For output width == stride
118    out_l.spare_capacity_mut().par_chunks_exact_mut(width).take(height).zip(
119        out_a.spare_capacity_mut().par_chunks_exact_mut(width).take(height).zip(
120            out_b.spare_capacity_mut().par_chunks_exact_mut(width).take(height))
121    ).enumerate()
122    .for_each(|(y, (l_row, (a_row, b_row)))| {
123        let in_row = &img.rows().nth(y).unwrap()[0..width];
124        let l_row = &mut l_row[0..width];
125        let a_row = &mut a_row[0..width];
126        let b_row = &mut b_row[0..width];
127        for x in 0..width {
128            let n = (x+11) ^ (y+11);
129            let (l,a,b) = cb(in_row[x], n);
130            l_row[x].write(l);
131            a_row[x].write(a);
132            b_row[x].write(b);
133        }
134    });
135
136    unsafe { out_l.set_len(area) };
137    unsafe { out_a.set_len(area) };
138    unsafe { out_b.set_len(area) };
139
140    vec![
141        Img::new(out_l, width, height),
142        Img::new(out_a, width, height),
143        Img::new(out_b, width, height),
144    ]
145}
146
147impl ToLABBitmap for ImgRef<'_, RGBAPLU> {
148    #[inline]
149    fn to_lab(&self) -> Vec<GBitmap> {
150        rgb_to_lab(*self, |px, n|{
151            px.to_rgb(n).to_lab()
152        })
153    }
154}
155
156impl ToLABBitmap for ImgRef<'_, RGBLU> {
157    #[inline]
158    fn to_lab(&self) -> Vec<GBitmap> {
159        rgb_to_lab(*self, |px, _n|{
160            px.to_lab()
161        })
162    }
163}
164
165#[test]
166fn cbrts1() {
167    let mut totaldiff = 0.;
168    let mut maxdiff: f64 = 0.;
169    for i in (0..=10001).rev() {
170        let x = (f64::from(i) / 10001.) as f32;
171        let a = cbrt_poly(x);
172        let actual = a * a * a;
173        let expected = x;
174        let absdiff = (f64::from(expected) - f64::from(actual)).abs();
175        assert!(absdiff < 0.0002, "{expected} - {actual} = {} @ {x}", expected - actual);
176        if i % 400 == 0 {
177            println!("{:+0.3}", (expected - actual) * 255.);
178        }
179        totaldiff += absdiff;
180        maxdiff = maxdiff.max(absdiff);
181    }
182    println!("1={totaldiff:0.6}; {maxdiff:0.8}");
183    assert!(totaldiff < 0.0025, "{totaldiff}");
184}
185
186#[test]
187fn cbrts2() {
188    let mut totaldiff = 0.;
189    let mut maxdiff: f64 = 0.;
190    for i in (2000..=10001).rev() {
191        let x = f64::from(i) / 10001.;
192        let actual = f64::from(cbrt_poly(x as f32));
193        let expected = x.cbrt();
194        let absdiff = (expected - actual).abs();
195        totaldiff += absdiff;
196        maxdiff = maxdiff.max(absdiff);
197        assert!(absdiff < 0.0000005, "{expected} - {actual} = {} @ {x}", expected - actual);
198    }
199    println!("2={totaldiff:0.6}; {maxdiff:0.8}");
200    assert!(totaldiff < 0.0025, "{totaldiff}");
201}
dssim_core/tolab.rs

dssim_core/
tolab.rs
