moxcms/
jzazbz.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 3/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::Xyz;
30use crate::jzczhz::Jzczhz;
31use crate::mlaf::mlaf;
32use num_traits::Pow;
33use pxfm::{dirty_powf, f_cbrtf, f_powf};
34use std::ops::{
35    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
36};
37
38#[inline]
39fn perceptual_quantizer(x: f32) -> f32 {
40    if x <= 0. {
41        return 0.;
42    }
43    let xx = dirty_powf(x * 1e-4, 0.1593017578125);
44    let rs = dirty_powf(
45        mlaf(0.8359375, 18.8515625, xx) / mlaf(1., 18.6875, xx),
46        134.034375,
47    );
48    if rs.is_nan() {
49        return 0.;
50    }
51    rs
52}
53
54#[inline]
55fn perceptual_quantizer_inverse(x: f32) -> f32 {
56    if x <= 0. {
57        return 0.;
58    }
59    let xx = dirty_powf(x, 7.460772656268214e-03);
60    let rs = 1e4
61        * dirty_powf(
62            (0.8359375 - xx) / mlaf(-18.8515625, 18.6875, xx),
63            6.277394636015326,
64        );
65    if rs.is_nan() {
66        return 0.;
67    }
68    rs
69}
70
71#[repr(C)]
72#[derive(Debug, Copy, Clone, PartialOrd, PartialEq, Default)]
73/// Represents Jzazbz
74pub struct Jzazbz {
75    /// Jz(lightness) generally expects to be between `0.0..1.0`.
76    pub jz: f32,
77    /// Az generally expects to be between `-0.5..0.5`.
78    pub az: f32,
79    /// Bz generally expects to be between `-0.5..0.5`.
80    pub bz: f32,
81}
82
83impl Jzazbz {
84    /// Constructs new instance
85    #[inline]
86    pub fn new(jz: f32, az: f32, bz: f32) -> Jzazbz {
87        Jzazbz { jz, az, bz }
88    }
89
90    /// Creates new [Jzazbz] from CIE [Xyz].
91    ///
92    /// JzAzBz is defined in D65 white point, adapt XYZ if needed first.
93    #[inline]
94    pub fn from_xyz(xyz: Xyz) -> Jzazbz {
95        Self::from_xyz_with_display_luminance(xyz, 200.)
96    }
97
98    /// Creates new [Jzazbz] from CIE [Xyz].
99    ///
100    /// JzAzBz is defined in D65 white point, adapt XYZ if needed first.
101    #[inline]
102    pub fn from_xyz_with_display_luminance(xyz: Xyz, display_luminance: f32) -> Jzazbz {
103        let abs_xyz = xyz * display_luminance;
104        let lp = perceptual_quantizer(mlaf(
105            mlaf(0.674207838 * abs_xyz.x, 0.382799340, abs_xyz.y),
106            -0.047570458,
107            abs_xyz.z,
108        ));
109        let mp = perceptual_quantizer(mlaf(
110            mlaf(0.149284160 * abs_xyz.x, 0.739628340, abs_xyz.y),
111            0.083327300,
112            abs_xyz.z,
113        ));
114        let sp = perceptual_quantizer(mlaf(
115            mlaf(0.070941080 * abs_xyz.x, 0.174768000, abs_xyz.y),
116            0.670970020,
117            abs_xyz.z,
118        ));
119        let iz = 0.5 * (lp + mp);
120        let az = mlaf(mlaf(3.524000 * lp, -4.066708, mp), 0.542708, sp);
121        let bz = mlaf(mlaf(0.199076 * lp, 1.096799, mp), -1.295875, sp);
122        let jz = (0.44 * iz) / mlaf(1., -0.56, iz) - 1.6295499532821566e-11;
123        Jzazbz::new(jz, az, bz)
124    }
125
126    /// Converts [Jzazbz] to [Xyz] D65
127    #[inline]
128    pub fn to_xyz(&self, display_luminance: f32) -> Xyz {
129        let jz = self.jz + 1.6295499532821566e-11;
130
131        let iz = jz / mlaf(0.44f32, 0.56, jz);
132        let l = perceptual_quantizer_inverse(mlaf(
133            mlaf(iz, 1.386050432715393e-1, self.az),
134            5.804731615611869e-2,
135            self.bz,
136        ));
137        let m = perceptual_quantizer_inverse(mlaf(
138            mlaf(iz, -1.386050432715393e-1, self.az),
139            -5.804731615611891e-2,
140            self.bz,
141        ));
142        let s = perceptual_quantizer_inverse(mlaf(
143            mlaf(iz, -9.601924202631895e-2, self.az),
144            -8.118918960560390e-1,
145            self.bz,
146        ));
147        let x = mlaf(
148            mlaf(1.661373055774069e+00 * l, -9.145230923250668e-01, m),
149            2.313620767186147e-01,
150            s,
151        );
152        let y = mlaf(
153            mlaf(-3.250758740427037e-01 * l, 1.571847038366936e+00, m),
154            -2.182538318672940e-01,
155            s,
156        );
157        let z = mlaf(
158            mlaf(-9.098281098284756e-02 * l, -3.127282905230740e-01, m),
159            1.522766561305260e+00,
160            s,
161        );
162        let rel_luminance = 1f32 / display_luminance;
163        Xyz::new(x, y, z) * rel_luminance
164    }
165
166    /// Converts into *Jzczhz*
167    #[inline]
168    pub fn to_jzczhz(&self) -> Jzczhz {
169        Jzczhz::from_jzazbz(*self)
170    }
171
172    #[inline]
173    pub fn euclidean_distance(&self, other: Self) -> f32 {
174        let djz = self.jz - other.jz;
175        let daz = self.az - other.az;
176        let dbz = self.bz - other.bz;
177        (djz * djz + daz * daz + dbz * dbz).sqrt()
178    }
179
180    #[inline]
181    pub fn taxicab_distance(&self, other: Self) -> f32 {
182        let djz = self.jz - other.jz;
183        let daz = self.az - other.az;
184        let dbz = self.bz - other.bz;
185        djz.abs() + daz.abs() + dbz.abs()
186    }
187}
188
189impl Index<usize> for Jzazbz {
190    type Output = f32;
191
192    #[inline]
193    fn index(&self, index: usize) -> &f32 {
194        match index {
195            0 => &self.jz,
196            1 => &self.az,
197            2 => &self.bz,
198            _ => panic!("Index out of bounds for Jzazbz"),
199        }
200    }
201}
202
203impl IndexMut<usize> for Jzazbz {
204    #[inline]
205    fn index_mut(&mut self, index: usize) -> &mut f32 {
206        match index {
207            0 => &mut self.jz,
208            1 => &mut self.az,
209            2 => &mut self.bz,
210            _ => panic!("Index out of bounds for Jzazbz"),
211        }
212    }
213}
214
215impl Add<f32> for Jzazbz {
216    type Output = Jzazbz;
217
218    #[inline]
219    fn add(self, rhs: f32) -> Self::Output {
220        Jzazbz::new(self.jz + rhs, self.az + rhs, self.bz + rhs)
221    }
222}
223
224impl Sub<f32> for Jzazbz {
225    type Output = Jzazbz;
226
227    #[inline]
228    fn sub(self, rhs: f32) -> Self::Output {
229        Jzazbz::new(self.jz - rhs, self.az - rhs, self.bz - rhs)
230    }
231}
232
233impl Mul<f32> for Jzazbz {
234    type Output = Jzazbz;
235
236    #[inline]
237    fn mul(self, rhs: f32) -> Self::Output {
238        Jzazbz::new(self.jz * rhs, self.az * rhs, self.bz * rhs)
239    }
240}
241
242impl Div<f32> for Jzazbz {
243    type Output = Jzazbz;
244
245    #[inline]
246    fn div(self, rhs: f32) -> Self::Output {
247        Jzazbz::new(self.jz / rhs, self.az / rhs, self.bz / rhs)
248    }
249}
250
251impl Add<Jzazbz> for Jzazbz {
252    type Output = Jzazbz;
253
254    #[inline]
255    fn add(self, rhs: Jzazbz) -> Self::Output {
256        Jzazbz::new(self.jz + rhs.jz, self.az + rhs.az, self.bz + rhs.bz)
257    }
258}
259
260impl Sub<Jzazbz> for Jzazbz {
261    type Output = Jzazbz;
262
263    #[inline]
264    fn sub(self, rhs: Jzazbz) -> Self::Output {
265        Jzazbz::new(self.jz - rhs.jz, self.az - rhs.az, self.bz - rhs.bz)
266    }
267}
268
269impl Mul<Jzazbz> for Jzazbz {
270    type Output = Jzazbz;
271
272    #[inline]
273    fn mul(self, rhs: Jzazbz) -> Self::Output {
274        Jzazbz::new(self.jz * rhs.jz, self.az * rhs.az, self.bz * rhs.bz)
275    }
276}
277
278impl Div<Jzazbz> for Jzazbz {
279    type Output = Jzazbz;
280
281    #[inline]
282    fn div(self, rhs: Jzazbz) -> Self::Output {
283        Jzazbz::new(self.jz / rhs.jz, self.az / rhs.az, self.bz / rhs.bz)
284    }
285}
286
287impl AddAssign<Jzazbz> for Jzazbz {
288    #[inline]
289    fn add_assign(&mut self, rhs: Jzazbz) {
290        self.jz += rhs.jz;
291        self.az += rhs.az;
292        self.bz += rhs.bz;
293    }
294}
295
296impl SubAssign<Jzazbz> for Jzazbz {
297    #[inline]
298    fn sub_assign(&mut self, rhs: Jzazbz) {
299        self.jz -= rhs.jz;
300        self.az -= rhs.az;
301        self.bz -= rhs.bz;
302    }
303}
304
305impl MulAssign<Jzazbz> for Jzazbz {
306    #[inline]
307    fn mul_assign(&mut self, rhs: Jzazbz) {
308        self.jz *= rhs.jz;
309        self.az *= rhs.az;
310        self.bz *= rhs.bz;
311    }
312}
313
314impl DivAssign<Jzazbz> for Jzazbz {
315    #[inline]
316    fn div_assign(&mut self, rhs: Jzazbz) {
317        self.jz /= rhs.jz;
318        self.az /= rhs.az;
319        self.bz /= rhs.bz;
320    }
321}
322
323impl AddAssign<f32> for Jzazbz {
324    #[inline]
325    fn add_assign(&mut self, rhs: f32) {
326        self.jz += rhs;
327        self.az += rhs;
328        self.bz += rhs;
329    }
330}
331
332impl SubAssign<f32> for Jzazbz {
333    #[inline]
334    fn sub_assign(&mut self, rhs: f32) {
335        self.jz -= rhs;
336        self.az -= rhs;
337        self.bz -= rhs;
338    }
339}
340
341impl MulAssign<f32> for Jzazbz {
342    #[inline]
343    fn mul_assign(&mut self, rhs: f32) {
344        self.jz *= rhs;
345        self.az *= rhs;
346        self.bz *= rhs;
347    }
348}
349
350impl DivAssign<f32> for Jzazbz {
351    #[inline]
352    fn div_assign(&mut self, rhs: f32) {
353        self.jz /= rhs;
354        self.az /= rhs;
355        self.bz /= rhs;
356    }
357}
358
359impl Neg for Jzazbz {
360    type Output = Jzazbz;
361
362    #[inline]
363    fn neg(self) -> Self::Output {
364        Jzazbz::new(-self.jz, -self.az, -self.bz)
365    }
366}
367
368impl Jzazbz {
369    #[inline]
370    pub fn sqrt(&self) -> Jzazbz {
371        Jzazbz::new(self.jz.sqrt(), self.az.sqrt(), self.bz.sqrt())
372    }
373
374    #[inline]
375    pub fn cbrt(&self) -> Jzazbz {
376        Jzazbz::new(f_cbrtf(self.jz), f_cbrtf(self.az), f_cbrtf(self.bz))
377    }
378}
379
380impl Pow<f32> for Jzazbz {
381    type Output = Jzazbz;
382
383    #[inline]
384    fn pow(self, rhs: f32) -> Self::Output {
385        Jzazbz::new(
386            f_powf(self.jz, rhs),
387            f_powf(self.az, rhs),
388            f_powf(self.bz, rhs),
389        )
390    }
391}
392
393impl Pow<Jzazbz> for Jzazbz {
394    type Output = Jzazbz;
395
396    #[inline]
397    fn pow(self, rhs: Jzazbz) -> Self::Output {
398        Jzazbz::new(
399            f_powf(self.jz, rhs.jz),
400            f_powf(self.az, self.az),
401            f_powf(self.bz, self.bz),
402        )
403    }
404}
405
406#[cfg(test)]
407mod tests {
408    use super::*;
409
410    #[test]
411    fn jzazbz_round() {
412        let xyz = Xyz::new(0.5, 0.4, 0.3);
413        let jzazbz = Jzazbz::from_xyz_with_display_luminance(xyz, 253f32);
414        let old_xyz = jzazbz.to_xyz(253f32);
415        assert!(
416            (xyz.x - old_xyz.x).abs() <= 1e-3,
417            "{:?} != {:?}",
418            xyz,
419            old_xyz
420        );
421        assert!(
422            (xyz.y - old_xyz.y).abs() <= 1e-3,
423            "{:?} != {:?}",
424            xyz,
425            old_xyz
426        );
427        assert!(
428            (xyz.z - old_xyz.z).abs() <= 1e-3,
429            "{:?} != {:?}",
430            xyz,
431            old_xyz
432        );
433    }
434}
moxcms/jzazbz.rs

moxcms/
jzazbz.rs
