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

dssim_core/
tolab.rs
