1 //!This module provides simple matrix operations on 3x3 matrix to aid in chromatic adaptation and
2 //!conversion calculations.
3
4 use float::Float;
5
6 use core::marker::PhantomData;
7
8 use {Component, Xyz};
9 use white_point::WhitePoint;
10 use rgb::{Primaries, Rgb, RgbSpace};
11 use encoding::Linear;
12 use convert::IntoColor;
13
14 ///A 9 element array representing a 3x3 matrix
15 pub type Mat3<T> = [T; 9];
16
17 ///Multiply the 3x3 matrix with the XYZ color
multiply_xyz<Swp: WhitePoint, Dwp: WhitePoint, T: Component + Float>( c: &Mat3<T>, f: &Xyz<Swp, T>, ) -> Xyz<Dwp, T>18 pub fn multiply_xyz<Swp: WhitePoint, Dwp: WhitePoint, T: Component + Float>(
19 c: &Mat3<T>,
20 f: &Xyz<Swp, T>,
21 ) -> Xyz<Dwp, T> {
22 Xyz {
23 x: (c[0] * f.x) + (c[1] * f.y) + (c[2] * f.z),
24 y: (c[3] * f.x) + (c[4] * f.y) + (c[5] * f.z),
25 z: (c[6] * f.x) + (c[7] * f.y) + (c[8] * f.z),
26 white_point: PhantomData,
27 }
28 }
29 ///Multiply the 3x3 matrix with the XYZ color into RGB color
multiply_xyz_to_rgb<S: RgbSpace, T: Component + Float>( c: &Mat3<T>, f: &Xyz<S::WhitePoint, T>, ) -> Rgb<Linear<S>, T>30 pub fn multiply_xyz_to_rgb<S: RgbSpace, T: Component + Float>(
31 c: &Mat3<T>,
32 f: &Xyz<S::WhitePoint, T>,
33 ) -> Rgb<Linear<S>, T> {
34 Rgb {
35 red: (c[0] * f.x) + (c[1] * f.y) + (c[2] * f.z),
36 green: (c[3] * f.x) + (c[4] * f.y) + (c[5] * f.z),
37 blue: (c[6] * f.x) + (c[7] * f.y) + (c[8] * f.z),
38 standard: PhantomData,
39 }
40 }
41 ///Multiply the 3x3 matrix with the RGB into XYZ color
multiply_rgb_to_xyz<S: RgbSpace, T: Component + Float>( c: &Mat3<T>, f: &Rgb<Linear<S>, T>, ) -> Xyz<S::WhitePoint, T>42 pub fn multiply_rgb_to_xyz<S: RgbSpace, T: Component + Float>(
43 c: &Mat3<T>,
44 f: &Rgb<Linear<S>, T>,
45 ) -> Xyz<S::WhitePoint, T> {
46 Xyz {
47 x: (c[0] * f.red) + (c[1] * f.green) + (c[2] * f.blue),
48 y: (c[3] * f.red) + (c[4] * f.green) + (c[5] * f.blue),
49 z: (c[6] * f.red) + (c[7] * f.green) + (c[8] * f.blue),
50 white_point: PhantomData,
51 }
52 }
53
54 ///Multiply a 3x3 matrix with another 3x3 matrix
multiply_3x3<T: Float>(c: &Mat3<T>, f: &Mat3<T>) -> Mat3<T>55 pub fn multiply_3x3<T: Float>(c: &Mat3<T>, f: &Mat3<T>) -> Mat3<T> {
56 let mut out = [T::zero(); 9];
57 out[0] = c[0] * f[0] + c[1] * f[3] + c[2] * f[6];
58 out[1] = c[0] * f[1] + c[1] * f[4] + c[2] * f[7];
59 out[2] = c[0] * f[2] + c[1] * f[5] + c[2] * f[8];
60
61 out[3] = c[3] * f[0] + c[4] * f[3] + c[5] * f[6];
62 out[4] = c[3] * f[1] + c[4] * f[4] + c[5] * f[7];
63 out[5] = c[3] * f[2] + c[4] * f[5] + c[5] * f[8];
64
65 out[6] = c[6] * f[0] + c[7] * f[3] + c[8] * f[6];
66 out[7] = c[6] * f[1] + c[7] * f[4] + c[8] * f[7];
67 out[8] = c[6] * f[2] + c[7] * f[5] + c[8] * f[8];
68
69 out
70 }
71
72 ///Invert a 3x3 matrix and panic if matrix is not invertable.
matrix_inverse<T: Float>(a: &Mat3<T>) -> Mat3<T>73 pub fn matrix_inverse<T: Float>(a: &Mat3<T>) -> Mat3<T> {
74 let d0 = a[4] * a[8] - a[5] * a[7];
75 let d1 = a[3] * a[8] - a[5] * a[6];
76 let d2 = a[3] * a[7] - a[4] * a[6];
77 let det = a[0] * d0 - a[1] * d1 + a[2] * d2;
78 if !det.is_normal() {
79 panic!("The given matrix is not invertible")
80 }
81 let d3 = a[1] * a[8] - a[2] * a[7];
82 let d4 = a[0] * a[8] - a[2] * a[6];
83 let d5 = a[0] * a[7] - a[1] * a[6];
84 let d6 = a[1] * a[5] - a[2] * a[4];
85 let d7 = a[0] * a[5] - a[2] * a[3];
86 let d8 = a[0] * a[4] - a[1] * a[3];
87
88 [
89 d0 / det,
90 -d3 / det,
91 d6 / det,
92 -d1 / det,
93 d4 / det,
94 -d7 / det,
95 d2 / det,
96 -d5 / det,
97 d8 / det,
98 ]
99 }
100
101 ///Geneartes to Srgb to Xyz transformation matrix for the given white point
rgb_to_xyz_matrix<S: RgbSpace, T: Component + Float>() -> Mat3<T>102 pub fn rgb_to_xyz_matrix<S: RgbSpace, T: Component + Float>() -> Mat3<T> {
103 let r: Xyz<S::WhitePoint, T> = S::Primaries::red().into_xyz();
104 let g: Xyz<S::WhitePoint, T> = S::Primaries::green().into_xyz();
105 let b: Xyz<S::WhitePoint, T> = S::Primaries::blue().into_xyz();
106
107 let mut transform_matrix = mat3_from_primaries(r, g, b);
108
109 let s_matrix: Rgb<Linear<S>, T> = multiply_xyz_to_rgb(
110 &matrix_inverse(&transform_matrix),
111 &S::WhitePoint::get_xyz(),
112 );
113 transform_matrix[0] = transform_matrix[0] * s_matrix.red;
114 transform_matrix[1] = transform_matrix[1] * s_matrix.green;
115 transform_matrix[2] = transform_matrix[2] * s_matrix.blue;
116 transform_matrix[3] = transform_matrix[3] * s_matrix.red;
117 transform_matrix[4] = transform_matrix[4] * s_matrix.green;
118 transform_matrix[5] = transform_matrix[5] * s_matrix.blue;
119 transform_matrix[6] = transform_matrix[6] * s_matrix.red;
120 transform_matrix[7] = transform_matrix[7] * s_matrix.green;
121 transform_matrix[8] = transform_matrix[8] * s_matrix.blue;
122
123 transform_matrix
124 }
125
126 #[cfg_attr(rustfmt, rustfmt_skip)]
mat3_from_primaries<T: Component + Float, Wp: WhitePoint>(r: Xyz<Wp, T>, g: Xyz<Wp, T>, b: Xyz<Wp, T>) -> Mat3<T>127 fn mat3_from_primaries<T: Component + Float, Wp: WhitePoint>(r: Xyz<Wp, T>, g: Xyz<Wp, T>, b: Xyz<Wp, T>) -> Mat3<T> {
128 [
129 r.x, g.x, b.x,
130 r.y, g.y, b.y,
131 r.z, g.z, b.z,
132 ]
133 }
134
135 #[cfg(test)]
136 mod test {
137 use Xyz;
138 use rgb::Rgb;
139 use encoding::{Linear, Srgb};
140 use chromatic_adaptation::AdaptInto;
141 use white_point::D50;
142 use super::{matrix_inverse, multiply_xyz, rgb_to_xyz_matrix, multiply_3x3};
143
144 #[test]
matrix_multiply_3x3()145 fn matrix_multiply_3x3() {
146 let inp1 = [1.0, 2.0, 3.0, 3.0, 2.0, 1.0, 2.0, 1.0, 3.0];
147 let inp2 = [4.0, 5.0, 6.0, 6.0, 5.0, 4.0, 4.0, 6.0, 5.0];
148 let expected = [28.0, 33.0, 29.0, 28.0, 31.0, 31.0, 26.0, 33.0, 31.0];
149
150 let computed = multiply_3x3(&inp1, &inp2);
151 assert_eq!(expected, computed)
152 }
153
154 #[test]
matrix_multiply_xyz()155 fn matrix_multiply_xyz() {
156 let inp1 = [0.1, 0.2, 0.3, 0.3, 0.2, 0.1, 0.2, 0.1, 0.3];
157 let inp2 = Xyz::new(0.4, 0.6, 0.8);
158
159 let expected = Xyz::new(0.4, 0.32, 0.38);
160
161 let computed = multiply_xyz(&inp1, &inp2);
162 assert_relative_eq!(expected, computed)
163 }
164
165 #[test]
matrix_inverse_check_1()166 fn matrix_inverse_check_1() {
167 let input: [f64; 9] = [3.0, 0.0, 2.0, 2.0, 0.0, -2.0, 0.0, 1.0, 1.0];
168
169 let expected: [f64; 9] = [0.2, 0.2, 0.0, -0.2, 0.3, 1.0, 0.2, -0.3, 0.0];
170 let computed = matrix_inverse(&input);
171 assert_eq!(expected, computed);
172 }
173 #[test]
matrix_inverse_check_2()174 fn matrix_inverse_check_2() {
175 let input: [f64; 9] = [1.0, 0.0, 1.0, 0.0, 2.0, 1.0, 1.0, 1.0, 1.0];
176
177 let expected: [f64; 9] = [-1.0, -1.0, 2.0, -1.0, 0.0, 1.0, 2.0, 1.0, -2.0];
178 let computed = matrix_inverse(&input);
179 assert_eq!(expected, computed);
180 }
181 #[test]
182 #[should_panic]
matrix_inverse_panic()183 fn matrix_inverse_panic() {
184 let input: [f64; 9] = [1.0, 0.0, 0.0, 2.0, 0.0, 0.0, -4.0, 6.0, 1.0];
185 matrix_inverse(&input);
186 }
187
188 #[cfg_attr(rustfmt, rustfmt_skip)]
189 #[test]
d65_rgb_conversion_matrix()190 fn d65_rgb_conversion_matrix() {
191 let expected = [
192 0.4124564, 0.3575761, 0.1804375,
193 0.2126729, 0.7151522, 0.0721750,
194 0.0193339, 0.1191920, 0.9503041
195 ];
196 let computed = rgb_to_xyz_matrix::<Srgb, f64>();
197 for (e, c) in expected.iter().zip(computed.iter()) {
198 assert_relative_eq!(e, c, epsilon = 0.000001)
199 }
200 }
201
202 #[test]
d65_to_d50()203 fn d65_to_d50() {
204 let input: Rgb<Linear<Srgb>> = Rgb::new(1.0, 1.0, 1.0);
205 let expected: Rgb<Linear<(Srgb, D50)>> = Rgb::new(1.0, 1.0, 1.0);
206
207 let computed: Rgb<Linear<(Srgb, D50)>> = input.adapt_into();
208 assert_relative_eq!(expected, computed, epsilon = 0.000001);
209 }
210 }
211