1 module photog.color;
2 
3 import mir.ndslice : Slice, sliced, SliceKind;
4 
5 import photog.utils;
6 
7 // dfmt off
8 /**
9 RGB working spaces.
10 */
11 enum WorkingSpace
12 {
13     sRgb = "sRgb"
14 }
15 
16 /**
17 Color space conversion matrices.
18 */
19 private static immutable
20 {
21     
22     auto sRgbRgb2Xyz = [
23         0.4124564, 0.3575761, 0.1804375,
24         0.2126729, 0.7151522, 0.0721750,
25         0.0193339, 0.1191920, 0.9503041
26     ];
27 
28     auto sRgbRgb2XyzReverse = [
29         0.1804375, 0.3575761, 0.4124564,
30         0.0721750, 0.7151522, 0.2126729,
31         0.9503041, 0.1191920, 0.0193339,
32     ];
33 	
34 	auto sRgbXyz2Rgb = [
35         3.2404542, -1.5371385, -0.4985314,
36 		-0.9692660, 1.8760108, 0.0415560,
37 		0.0556434, -0.2040259, 1.0572252
38     ];
39 }
40 
41 /**
42 Standard illuminants (XYZ).
43 */
44 enum Illuminant
45 {
46     d65 = [0.95047, 1.0, 1.08883]
47 }
48 
49 /**
50 Chromatic adaptation methods.
51 */
52 enum ChromAdaptMethod
53 {
54     bradford = "bradford"
55 }
56 
57 /**
58 Chromatic adaptation method matrices.
59 */
60 private static immutable
61 {
62     auto bradfordXyz2Lms = [
63         0.8951, 0.2664, -0.1614,
64         -0.7502, 1.7135, 0.0367,
65         0.0389, -0.0685, 1.0296
66     ];
67     
68     auto bradfordInvXyz2Lms = [
69         0.9869929, -0.1470543, 0.1599627,
70         0.4323053, 0.5183603, 0.0492912,
71         -0.0085287, 0.0400428, 0.9684867
72     ];
73 }
74 // dfmt on
75 
76 /**
77 Convert RGB input to XYZ.
78 */
79 Slice!(ReturnType*, 3) rgb2Xyz(WorkingSpace workingSpace = WorkingSpace.sRgb, ReturnType = double,
80         T)(Slice!(T, 3) input)
81 {
82     return rgbBgr2Xyz!(false, workingSpace, ReturnType)(input);
83 }
84 
85 ///
86 unittest
87 {
88     import std.math : approxEqual;
89 
90     // dfmt off
91     Slice!(double*, 3) rgb = [
92         1.0, 0.0, 0.0,
93         0.0, 1.0, 0.0,
94         0.0, 0.0, 1.0,
95         0.470588, 0.470588, 0.470588
96     ].sliced(4, 1, 3);
97 
98     Slice!(double*, 3) xyz = [
99         0.412456, 0.212673, 0.019334,
100         0.357576, 0.715152, 0.119192,
101         0.180438, 0.072175, 0.950304,
102         0.178518, 0.187821, 0.204505
103     ].sliced(4, 1, 3);
104     // dfmt on
105 
106     assert(approxEqual(rgb.rgb2Xyz, xyz));
107 }
108 
109 /**
110 Convert BGR input to XYZ.
111 */
112 Slice!(ReturnType*, 3) bgr2Xyz(WorkingSpace workingSpace = WorkingSpace.sRgb, ReturnType = double,
113         T)(Slice!(T, 3) input)
114 {
115     return rgbBgr2Xyz!(true, workingSpace, ReturnType)(input);
116 }
117 
118 ///
119 unittest
120 {
121     import std.math : approxEqual;
122 
123     // dfmt off
124     Slice!(double*, 3) bgr = [
125         0.0, 0.0, 1.0,
126         0.0, 1.0, 0.0,
127         1.0, 0.0, 0.0,
128         0.470588, 0.470588, 0.470588
129     ].sliced(4, 1, 3);
130 
131     Slice!(double*, 3) xyz = [
132         0.412456, 0.212673, 0.019334,
133         0.357576, 0.715152, 0.119192,
134         0.180438, 0.072175, 0.950304,
135         0.178518, 0.187821, 0.204505
136     ].sliced(4, 1, 3);
137     // dfmt on
138 
139     assert(approxEqual(bgr.bgr2Xyz, xyz));
140 }
141 
142 private Slice!(ReturnType*, 3) rgbBgr2Xyz(bool isBgr, WorkingSpace workingSpace, ReturnType, T)(
143         Slice!(T, 3) input)
144 in
145 {
146     import std.traits : isFloatingPoint;
147 
148     static assert(isFloatingPoint!ReturnType, "Return type must be floating point.");
149     assert(input.shape[2] == 3, "Input requires 3 channels.");
150 }
151 do
152 {
153     import mir.ndslice : each, uninitSlice, zip;
154 
155     auto output = uninitSlice!ReturnType(input.shape);
156     auto inputPack = pixelPack!3(input);
157     auto outputPack = pixelPack!3(output);
158     auto pixelZip = zip(inputPack, outputPack);
159 
160     static if (isBgr)
161         auto conversionMatrix = mixin(workingSpace ~ "Rgb2XyzReverse").sliced(3, 3);
162     else
163         auto conversionMatrix = mixin(workingSpace ~ "Rgb2Xyz").sliced(3, 3);
164 
165     pixelZip.each!((z) {
166         z.rgbBgr2XyzImpl!(isBgr, workingSpace)(conversionMatrix);
167     });
168 
169     return output;
170 }
171 
172 private void rgbBgr2XyzImpl(bool isBgr, WorkingSpace workingSpace, T, U)(T pixelZip,
173         Slice!(U, 2) conversionMatrix)
174 {
175     import kaleidic.lubeck : mtimes;
176     import mir.ndslice : each;
177 
178     // pixelZip[0] = input
179     // pixelZip[1] = output
180 
181     auto output = sliced(pixelZip[0].ptr, 3);
182 
183     // Linearize
184     output.each!((ref chnl) {
185         if (chnl <= 0.04045)
186             chnl = chnl / 12.92;
187         else
188             chnl = ((chnl + 0.055) / 1.055) ^^ 2.4;
189     });
190 
191     // Convert to XYZ
192     output = conversionMatrix.mtimes(output);
193 
194     pixelZip[1][] = output.field;
195 }
196 
197 /**
198 Convert XYZ input to RGB.
199 */
200 Slice!(ReturnType*, 3) xyz2Rgb(WorkingSpace workingSpace = WorkingSpace.sRgb, ReturnType = double,
201         T)(Slice!(T, 3) input)
202 {
203     return xyz2RgbBgr!(false, workingSpace, ReturnType)(input);
204 }
205 
206 ///
207 unittest
208 {
209     import std.math : approxEqual;
210 
211     // dfmt off
212     Slice!(double*, 3) rgb = [
213         1.0, 0.0, 0.0,
214         0.0, 1.0, 0.0,
215         0.0, 0.0, 1.0,
216         0.470588, 0.470588, 0.470588
217     ].sliced(4, 1, 3);
218 
219     Slice!(double*, 3) xyz = [
220         0.412456, 0.212673, 0.019334,
221         0.357576, 0.715152, 0.119192,
222         0.180438, 0.072175, 0.950304,
223         0.178518, 0.187821, 0.204505
224     ].sliced(4, 1, 3);
225     // dfmt on
226 
227     assert(approxEqual(xyz.xyz2Rgb, rgb, 0.01, 1e-04));
228 }
229 
230 /**
231 Convert XYZ input to BGR.
232 */
233 Slice!(ReturnType*, 3) xyz2Bgr(WorkingSpace workingSpace = WorkingSpace.sRgb, ReturnType = double,
234         T)(Slice!(T, 3) input)
235 {
236     return xyz2RgbBgr!(true, workingSpace, ReturnType)(input);
237 }
238 
239 ///
240 unittest
241 {
242     import std.math : approxEqual;
243 
244     // dfmt off
245     Slice!(double*, 3) bgr = [
246         0.0, 0.0, 1.0,
247         0.0, 1.0, 0.0,
248         1.0, 0.0, 0.0,
249         0.470588, 0.470588, 0.470588
250     ].sliced(4, 1, 3);
251 
252     Slice!(double*, 3) xyz = [
253         0.412456, 0.212673, 0.019334,
254         0.357576, 0.715152, 0.119192,
255         0.180438, 0.072175, 0.950304,
256         0.178518, 0.187821, 0.204505
257     ].sliced(4, 1, 3);
258     // dfmt on
259 
260     assert(approxEqual(xyz.xyz2Bgr, bgr, 0.01, 1e-04));
261 }
262 
263 private Slice!(ReturnType*, 3) xyz2RgbBgr(bool isBgr, WorkingSpace workingSpace, ReturnType, T)(
264         Slice!(T, 3) input)
265 in
266 {
267     import std.traits : isFloatingPoint;
268 
269     static assert(isFloatingPoint!ReturnType, "Return type must be floating point.");
270 }
271 do
272 {
273     import mir.ndslice : each, uninitSlice, zip;
274 
275     auto output = uninitSlice!ReturnType(input.shape);
276     auto inputPack = pixelPack!3(input);
277     auto outputPack = pixelPack!3(output);
278     auto pixelZip = zip(inputPack, outputPack);
279     auto conversionMatrix = mixin(workingSpace ~ "Xyz2Rgb").sliced(3, 3); // TODO: Known at compile-time. Send as template?
280     pixelZip.each!((z) {
281         z.xyz2RgbBgrImpl!(isBgr, workingSpace)(conversionMatrix);
282     });
283     // No need to unpack output. Underlying data has already been changed.
284     return output;
285 }
286 
287 private void xyz2RgbBgrImpl(bool isBgr, WorkingSpace workingSpace, T, U)(T pixelZip,
288         Slice!(U, 2) conversionMatrix)
289 {
290     import kaleidic.lubeck : mtimes;
291     import mir.ndslice : each;
292 
293     // pixelZip[0] = input
294     // pixelZip[1] = output
295 
296     // Convert to RGB
297     // TODO: Is manual matrix math faster here?
298     auto output = conversionMatrix.mtimes(sliced(pixelZip[0].ptr, 3));
299 
300     // Un-linearize
301     output.each!((ref chnl) {
302         if (chnl <= 0.0031308)
303             chnl = chnl * 12.92;
304         else
305             chnl = (1.055 * chnl ^^ (1 / 2.4)) - 0.055;
306     });
307 
308     static if (isBgr)
309     {
310         pixelZip[1][] = [output.field[2], output.field[1], output.field[0]];
311     }
312     else
313     {
314         pixelZip[1][] = output.field;
315     }
316 }
317 
318 /**
319 Chromatically adapt RGB input from the given source illuminant to the given destination illuminant.
320 
321 We employ the von Kries coefficient law for chromatic adaptation. A cone response under a source
322 illuminant is converted to one under a destination illuminant via diagonal scaling of the cone
323 response components. The output of the chromatic adaptation between a pair of illuminants can be 
324 tweaked by chromatic adaption methods which define the conversion to and from XYZ to cone responses
325 (LMS).
326 
327 References:
328     http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html
329     https://en.wikipedia.org/wiki/Von_Kries_coefficient_law#Chromatic_adaptation
330     https://en.wikipedia.org/wiki/CIE_1931_color_space#Color_matching_functions
331     https://en.wikipedia.org/wiki/LMS_color_space
332     http://scottburns.us/chromatic-adaptation-transform-by-reflectance-reconstruction/
333 
334 Params:
335     method = Method to use for chromatic adaptation. Bradford by default.
336     workingSpace = Working space for both the input and the eventual output. sRGB by default.
337     input = RBG image in floating point form.
338     srcIlluminant = XYZ illuminant under which input was recorded.
339     destIlluminant = XYZ illuminant that input should be chromatically adapted to.
340 
341 Returns:
342     Chromatically adapted floating point RGB image.
343 */
344 Slice!(Iterator, 3) chromAdapt(ChromAdaptMethod method = ChromAdaptMethod.bradford,
345         WorkingSpace workingSpace = WorkingSpace.sRgb, Iterator)(Slice!(Iterator,
346         3) input, double[] srcIlluminant, double[] destIlluminant)
347 in
348 {
349     import std.traits : isFloatingPoint;
350 
351     static assert(isFloatingPoint!(IteratorType!Iterator), "Input values must be floating point.");
352     assert(srcIlluminant.length == 3);
353     assert(destIlluminant.length == 3);
354 }
355 do
356 {
357     import kaleidic.lubeck : mtimes;
358     import mir.ndslice : diagonal, each, reshape, slice, uninitSlice, zip;
359 
360     auto xyz2Lms = mixin(method ~ "Xyz2Lms");
361     auto lms2Xyz = mixin(method ~ "InvXyz2Lms");
362 
363     auto lmsSrc = xyz2Lms.sliced(3, 3).mtimes(srcIlluminant.sliced(3, 1));
364     auto lmsDest = xyz2Lms.sliced(3, 3).mtimes(destIlluminant.sliced(3, 1));
365 
366     int err;
367     auto lmsGain = slice!double([3, 3], 0);
368     auto diag = lmsGain.diagonal;
369     diag[] = (lmsDest / lmsSrc).reshape([3], err);
370     assert(err == 0);
371 
372     auto transform = lms2Xyz.sliced(3, 3).mtimes(lmsGain).mtimes(xyz2Lms.sliced(3, 3));
373 
374     auto output = uninitSlice!(IteratorType!Iterator)(input.shape);
375     auto inputPack = pixelPack!3(input.rgb2Xyz!workingSpace);
376     auto outputPack = pixelPack!3(output);
377     auto pixelZip = zip(inputPack, outputPack);
378 
379     pixelZip.each!((z) { z.chromAdaptImpl(transform); });
380 
381     return output.xyz2Rgb!workingSpace;
382 }
383 
384 private void chromAdaptImpl(T, U)(T pixelZip, Slice!(U, 2) transform)
385 {
386     import kaleidic.lubeck : mtimes;
387 
388     auto output = sliced(pixelZip[0].ptr, 3);
389 
390     output = transform.mtimes(output);
391 
392     pixelZip[1][] = output.field;
393 }
394 
395 ///
396 unittest
397 {
398     import mir.ndslice : reshape, slice;
399 
400     // dfmt off
401     ubyte[] uImage = [
402         255, 0, 0,
403         0, 255, 0,
404         0, 0, 255,
405         120, 120, 120
406     ];
407 
408     auto image = uImage.sliced(4, 1, 3).toFloating;
409 
410     int err;
411     double[] srcIlluminant = image
412         .imageMean
413         .reshape([1, 1, 3], err) // TODO: Get rid of the need for reshaping.
414         .rgb2Xyz
415         .field;
416     assert(err == 0);
417     //dfmt on
418 
419     auto ca = chromAdapt(image, srcIlluminant, Illuminant.d65);
420     // TODO: Add assertion on chromAdapt output.
421 }