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 }