1 /**
2 * @file colorspace.c
3 * @author Pascal Getreuer 2005-2010 <getreuer@gmail.com>
4 *
5 * == Summary ==
6 * This file implements routines for color transformations between the spaces
7 * sRGB, Y'UV, Y'CbCr, Y'PbPr, Y'DbDr, Y'IQ, HSV, HSL, HSI, CIEXYZ, CIELAB,
8 * CIELUV, CIELCH, and CIECAT02 LMS.
9 *
10 * == Usage ==
11 * First call GetColorTransform, specifying the source and destination color
12 * spaces as "dest<-src" or "src->dest". Then call ApplyColorTransform to
13 * perform the transform:
14 @code
15 double S[3] = {173, 0.8, 0.5};
16 double D[3];
17 colortransform Trans;
18
19 if(!(GetColorTransform(&Trans, "HSI -> Lab")))
20 {
21 printf("Invalid syntax or unknown color space\n");
22 return;
23 }
24
25 ApplyColorTransform(Trans, &D[0], &D[1], &D[2], S[0], S[1], S[2]);
26 @endcode
27 * "num" is a typedef defined at the beginning of colorspace.h that may be set
28 * to either double or float, depending on the application.
29 *
30 * Specific transformation routines can also be called directly. The following
31 * converts an sRGB color to CIELAB and then back to sRGB:
32 @code
33 double R = 0.85, G = 0.32, B = 0.5;
34 double L, a, b;
35 Rgb2Lab(&L, &a, &b, R, G, B);
36 Lab2Rgb(&R, &G, &B, L, a, b);
37 @endcode
38 * Generally, the calling syntax is
39 @code
40 Foo2Bar(&B0, &B1, &B2, F0, F1, F2);
41 @endcode
42 * where (F0,F1,F2) are the coordinates of a color in space "Foo" and
43 * (B0,B1,B2) are the transformed coordinates in space "Bar." For any
44 * transformation routine, its inverse has the sytax
45 @code
46 Bar2Foo(&F0, &F1, &F2, B0, B1, B2);
47 @endcode
48 *
49 * The conversion routines are consistently named with the first letter of a
50 * color space capitalized with following letters in lower case and omitting
51 * prime symbols. For example, "Rgb2Ydbdr" converts sRGB to Y'DbDr. For
52 * any transformation routine Foo2Bar, its inverse is Bar2Foo.
53 *
54 * All transformations assume a two degree observer angle and a D65 illuminant.
55 * The white point can be changed by modifying the WHITEPOINT_X, WHITEPOINT_Y,
56 * WHITEPOINT_Z definitions at the beginning of colorspace.h.
57 *
58 * == List of transformation routines ==
59 * - Rgb2Yuv(double *Y, double *U, double *V, double R, double G, double B)
60 * - Rgb2Ycbcr(double *Y, double *Cb, double *Cr, double R, double G, double B)
61 * - Rgb2Jpegycbcr(double *Y, double *Cb, double *Cr, double R, double G, double B)
62 * - Rgb2Ypbpr(double *Y, double *Pb, double *Pr, double R, double G, double B)
63 * - Rgb2Ydbdr(double *Y, double *Db, double *Dr, double R, double G, double B)
64 * - Rgb2Yiq(double *Y, double *I, double *Q, double R, double G, double B)
65 * - Rgb2Hsv(double *H, double *S, double *V, double R, double G, double B)
66 * - Rgb2Hsl(double *H, double *S, double *L, double R, double G, double B)
67 * - Rgb2Hsi(double *H, double *S, double *I, double R, double G, double B)
68 * - Rgb2Xyz(double *X, double *Y, double *Z, double R, double G, double B)
69 * - Xyz2Lab(double *L, double *a, double *b, double X, double Y, double Z)
70 * - Xyz2Luv(double *L, double *u, double *v, double X, double Y, double Z)
71 * - Xyz2Lch(double *L, double *C, double *h, double X, double Y, double Z)
72 * - Xyz2Cat02lms(double *L, double *M, double *S, double X, double Y, double Z)
73 * - Rgb2Lab(double *L, double *a, double *b, double R, double G, double B)
74 * - Rgb2Luv(double *L, double *u, double *v, double R, double G, double B)
75 * - Rgb2Lch(double *L, double *C, double *h, double R, double G, double B)
76 * - Rgb2Cat02lms(double *L, double *M, double *S, double R, double G, double B)
77 * (Similarly for the inverse transformations.)
78 *
79 * It is possible to transform between two arbitrary color spaces by first
80 * transforming from the source space to sRGB and then transforming from
81 * sRGB to the desired destination space. For transformations between CIE
82 * color spaces, it is convenient to use XYZ as the intermediate space. This
83 * is the strategy used by GetColorTransform and ApplyColorTransform.
84 *
85 * == References ==
86 * The definitions of these spaces and the many of the transformation formulas
87 * can be found in
88 *
89 * Poynton, "Frequently Asked Questions About Gamma"
90 * http://www.poynton.com/notes/colour_and_gamma/GammaFAQ.html
91 *
92 * Poynton, "Frequently Asked Questions About Color"
93 * http://www.poynton.com/notes/colour_and_gamma/ColorFAQ.html
94 *
95 * and Wikipedia articles
96 * http://en.wikipedia.org/wiki/SRGB
97 * http://en.wikipedia.org/wiki/YUV
98 * http://en.wikipedia.org/wiki/YCbCr
99 * http://en.wikipedia.org/wiki/YPbPr
100 * http://en.wikipedia.org/wiki/YDbDr
101 * http://en.wikipedia.org/wiki/YIQ
102 * http://en.wikipedia.org/wiki/HSL_and_HSV
103 * http://en.wikipedia.org/wiki/CIE_1931_color_space
104 * http://en.wikipedia.org/wiki/Lab_color_space
105 * http://en.wikipedia.org/wiki/CIELUV_color_space
106 * http://en.wikipedia.org/wiki/LMS_color_space
107 *
108 * == License (BSD) ==
109 * Copyright (c) 2005-2010, Pascal Getreuer
110 * All rights reserved.
111 *
112 * Redistribution and use in source and binary forms, with or without
113 * modification, are permitted provided that the following conditions are met:
114 *
115 * - Redistributions of source code must retain the above copyright
116 * notice, this list of conditions and the following disclaimer.
117 * - Redistributions in binary form must reproduce the above copyright
118 * notice, this list of conditions and the following disclaimer in
119 * the documentation and/or other materials provided with the distribution.
120 *
121 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
122 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
123 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
124 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
125 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
126 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
127 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
128 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
129 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
130 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
131 * POSSIBILITY OF SUCH DAMAGE.
132 */
133 #include <math.h>
134 #include <stdio.h>
135 #include <string.h>
136 #include <ctype.h>
137
138 #include "gtkmm2ext/colorspace.h"
139
140 namespace Gtkmm2ext {
141
142 /** @brief Min of A and B */
143 #define MIN(A,B) (((A) <= (B)) ? (A) : (B))
144
145 /** @brief Max of A and B */
146 #define MAX(A,B) (((A) >= (B)) ? (A) : (B))
147
148 /** @brief Min of A, B, and C */
149 #define MIN3(A,B,C) (((A) <= (B)) ? MIN(A,C) : MIN(B,C))
150
151 /** @brief Max of A, B, and C */
152 #define MAX3(A,B,C) (((A) >= (B)) ? MAX(A,C) : MAX(B,C))
153
154 #ifndef M_PI
155 /** @brief The constant pi */
156 #define M_PI 3.14159265358979323846264338327950288
157 #endif
158
159 /**
160 * @brief sRGB gamma correction, transforms R to R'
161 * http://en.wikipedia.org/wiki/SRGB
162 */
163 #define GAMMACORRECTION(t) \
164 (((t) <= 0.0031306684425005883) ? \
165 (12.92*(t)) : (1.055*pow((t), 0.416666666666666667) - 0.055))
166
167 /**
168 * @brief Inverse sRGB gamma correction, transforms R' to R
169 */
170 #define INVGAMMACORRECTION(t) \
171 (((t) <= 0.0404482362771076) ? \
172 ((t)/12.92) : pow(((t) + 0.055)/1.055, 2.4))
173
174 /**
175 * @brief CIE L*a*b* f function (used to convert XYZ to L*a*b*)
176 * http://en.wikipedia.org/wiki/Lab_color_space
177 */
178 #define LABF(t) \
179 ((t >= 8.85645167903563082e-3) ? \
180 pow(t,0.333333333333333) : (841.0/108.0)*(t) + (4.0/29.0))
181
182 /**
183 * @brief CIE L*a*b* inverse f function
184 * http://en.wikipedia.org/wiki/Lab_color_space
185 */
186 #define LABINVF(t) \
187 ((t >= 0.206896551724137931) ? \
188 ((t)*(t)*(t)) : (108.0/841.0)*((t) - (4.0/29.0)))
189
190 /** @brief u'v' coordinates of the white point for CIE Lu*v* */
191 #define WHITEPOINT_U ((4*WHITEPOINT_X) \
192 /(WHITEPOINT_X + 15*WHITEPOINT_Y + 3*WHITEPOINT_Z))
193 #define WHITEPOINT_V ((9*WHITEPOINT_Y) \
194 /(WHITEPOINT_X + 15*WHITEPOINT_Y + 3*WHITEPOINT_Z))
195
196 /** @brief Enumeration of the supported color spaces */
197 #define UNKNOWN_SPACE 0
198 #define RGB_SPACE 1
199 #define YUV_SPACE 2
200 #define YCBCR_SPACE 3
201 #define JPEGYCBCR_SPACE 4
202 #define YPBPR_SPACE 5
203 #define YDBDR_SPACE 6
204 #define YIQ_SPACE 7
205 #define HSV_SPACE 8
206 #define HSL_SPACE 9
207 #define HSI_SPACE 10
208 #define XYZ_SPACE 11
209 #define LAB_SPACE 12
210 #define LUV_SPACE 13
211 #define LCH_SPACE 14
212 #define CAT02LMS_SPACE 15
213
214 #define NUM_TRANSFORM_PAIRS 18
215
216
217
218 /*
219 * == Linear color transformations ==
220 *
221 * The following routines implement transformations between sRGB and
222 * the linearly-related color spaces Y'UV, Y'PbPr, Y'DbDr, and Y'IQ.
223 */
224
225
226 /**
227 * @brief Convert sRGB to NTSC/PAL Y'UV Luma + Chroma
228 *
229 * @param Y, U, V pointers to hold the result
230 * @param R, G, B the input sRGB values
231 *
232 * Wikipedia: http://en.wikipedia.org/wiki/YUV
233 */
Rgb2Yuv(double * Y,double * U,double * V,double R,double G,double B)234 void Rgb2Yuv(double *Y, double *U, double *V, double R, double G, double B)
235 {
236 *Y = (double)( 0.299*R + 0.587*G + 0.114*B);
237 *U = (double)(-0.147*R - 0.289*G + 0.436*B);
238 *V = (double)( 0.615*R - 0.515*G - 0.100*B);
239 }
240
241
242 /**
243 * @brief Convert NTSC/PAL Y'UV to sRGB
244 *
245 * @param R, G, B pointers to hold the result
246 * @param Y, U, V the input YUV values
247 */
Yuv2Rgb(double * R,double * G,double * B,double Y,double U,double V)248 void Yuv2Rgb(double *R, double *G, double *B, double Y, double U, double V)
249 {
250 *R = (double)(Y - 3.945707070708279e-05*U + 1.1398279671717170825*V);
251 *G = (double)(Y - 0.3946101641414141437*U - 0.5805003156565656797*V);
252 *B = (double)(Y + 2.0319996843434342537*U - 4.813762626262513e-04*V);
253 }
254
255
256 /** @brief sRGB to Y'CbCr Luma + Chroma */
Rgb2Ycbcr(double * Y,double * Cb,double * Cr,double R,double G,double B)257 void Rgb2Ycbcr(double *Y, double *Cb, double *Cr, double R, double G, double B)
258 {
259 *Y = (double)( 65.481*R + 128.553*G + 24.966*B + 16);
260 *Cb = (double)(-37.797*R - 74.203*G + 112.0 *B + 128);
261 *Cr = (double)(112.0 *R - 93.786*G - 18.214*B + 128);
262 }
263
264
265 /** @brief Y'CbCr to sRGB */
Ycbcr2Rgb(double * R,double * G,double * B,double Y,double Cr,double Cb)266 void Ycbcr2Rgb(double *R, double *G, double *B, double Y, double Cr, double Cb)
267 {
268 Y -= 16;
269 Cb -= 128;
270 Cr -= 128;
271 *R = (double)(0.00456621004566210107*Y + 1.1808799897946415e-09*Cr + 0.00625892896994393634*Cb);
272 *G = (double)(0.00456621004566210107*Y - 0.00153632368604490212*Cr - 0.00318811094965570701*Cb);
273 *B = (double)(0.00456621004566210107*Y + 0.00791071623355474145*Cr + 1.1977497040190077e-08*Cb);
274 }
275
276
277 /** @brief sRGB to JPEG-Y'CbCr Luma + Chroma */
Rgb2Jpegycbcr(double * Y,double * Cb,double * Cr,double R,double G,double B)278 void Rgb2Jpegycbcr(double *Y, double *Cb, double *Cr, double R, double G, double B)
279 {
280 Rgb2Ypbpr(Y, Cb, Cr, R, G, B);
281 *Cb += (double)0.5;
282 *Cr += (double)0.5;
283 }
284
285 /** @brief JPEG-Y'CbCr to sRGB */
Jpegycbcr2Rgb(double * R,double * G,double * B,double Y,double Cb,double Cr)286 void Jpegycbcr2Rgb(double *R, double *G, double *B, double Y, double Cb, double Cr)
287 {
288 Cb -= (double)0.5;
289 Cr -= (double)0.5;
290 Ypbpr2Rgb(R, G, B, Y, Cb, Cr);
291 }
292
293
294 /** @brief sRGB to Y'PbPr Luma (ITU-R BT.601) + Chroma */
Rgb2Ypbpr(double * Y,double * Pb,double * Pr,double R,double G,double B)295 void Rgb2Ypbpr(double *Y, double *Pb, double *Pr, double R, double G, double B)
296 {
297 *Y = (double)( 0.299 *R + 0.587 *G + 0.114 *B);
298 *Pb = (double)(-0.1687367*R - 0.331264*G + 0.5 *B);
299 *Pr = (double)( 0.5 *R - 0.418688*G - 0.081312*B);
300 }
301
302
303 /** @brief Y'PbPr to sRGB */
Ypbpr2Rgb(double * R,double * G,double * B,double Y,double Pb,double Pr)304 void Ypbpr2Rgb(double *R, double *G, double *B, double Y, double Pb, double Pr)
305 {
306 *R = (double)(0.99999999999914679361*Y - 1.2188941887145875e-06*Pb + 1.4019995886561440468*Pr);
307 *G = (double)(0.99999975910502514331*Y - 0.34413567816504303521*Pb - 0.71413649331646789076*Pr);
308 *B = (double)(1.00000124040004623180*Y + 1.77200006607230409200*Pb + 2.1453384174593273e-06*Pr);
309 }
310
311
312 /** @brief sRGB to SECAM Y'DbDr Luma + Chroma */
Rgb2Ydbdr(double * Y,double * Db,double * Dr,double R,double G,double B)313 void Rgb2Ydbdr(double *Y, double *Db, double *Dr, double R, double G, double B)
314 {
315 *Y = (double)( 0.299*R + 0.587*G + 0.114*B);
316 *Db = (double)(-0.450*R - 0.883*G + 1.333*B);
317 *Dr = (double)(-1.333*R + 1.116*G + 0.217*B);
318 }
319
320
321 /** @brief SECAM Y'DbDr to sRGB */
Ydbdr2Rgb(double * R,double * G,double * B,double Y,double Db,double Dr)322 void Ydbdr2Rgb(double *R, double *G, double *B, double Y, double Db, double Dr)
323 {
324 *R = (double)(Y + 9.2303716147657e-05*Db - 0.52591263066186533*Dr);
325 *G = (double)(Y - 0.12913289889050927*Db + 0.26789932820759876*Dr);
326 *B = (double)(Y + 0.66467905997895482*Db - 7.9202543533108e-05*Dr);
327 }
328
329
330 /** @brief sRGB to NTSC YIQ */
Rgb2Yiq(double * Y,double * I,double * Q,double R,double G,double B)331 void Rgb2Yiq(double *Y, double *I, double *Q, double R, double G, double B)
332 {
333 *Y = (double)(0.299 *R + 0.587 *G + 0.114 *B);
334 *I = (double)(0.595716*R - 0.274453*G - 0.321263*B);
335 *Q = (double)(0.211456*R - 0.522591*G + 0.311135*B);
336 }
337
338
339 /** @brief Convert NTSC YIQ to sRGB */
Yiq2Rgb(double * R,double * G,double * B,double Y,double I,double Q)340 void Yiq2Rgb(double *R, double *G, double *B, double Y, double I, double Q)
341 {
342 *R = (double)(Y + 0.9562957197589482261*I + 0.6210244164652610754*Q);
343 *G = (double)(Y - 0.2721220993185104464*I - 0.6473805968256950427*Q);
344 *B = (double)(Y - 1.1069890167364901945*I + 1.7046149983646481374*Q);
345 }
346
347
348
349 /*
350 * == Hue Saturation Value/Lightness/Intensity color transformations ==
351 *
352 * The following routines implement transformations between sRGB and
353 * color spaces HSV, HSL, and HSI.
354 */
355
356
357 /**
358 * @brief Convert an sRGB color to Hue-Saturation-Value (HSV)
359 *
360 * @param H, S, V pointers to hold the result
361 * @param R, G, B the input sRGB values scaled in [0,1]
362 *
363 * This routine transforms from sRGB to the hexcone HSV color space. The
364 * sRGB values are assumed to be between 0 and 1. The output values are
365 * H = hexagonal hue angle (0 <= H < 360),
366 * S = C/V (0 <= S <= 1),
367 * V = max(R',G',B') (0 <= V <= 1),
368 * where C = max(R',G',B') - min(R',G',B'). The inverse color transformation
369 * is given by Hsv2Rgb.
370 *
371 * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
372 */
Rgb2Hsv(double * H,double * S,double * V,double R,double G,double B)373 void Rgb2Hsv(double *H, double *S, double *V, double R, double G, double B)
374 {
375 double Max = MAX3(R, G, B);
376 double Min = MIN3(R, G, B);
377 double C = Max - Min;
378
379
380 *V = Max;
381
382 if(C > 0)
383 {
384 if(Max == R)
385 {
386 *H = (G - B) / C;
387
388 if(G < B)
389 *H += 6;
390 }
391 else if(Max == G)
392 *H = 2 + (B - R) / C;
393 else
394 *H = 4 + (R - G) / C;
395
396 *H *= 60;
397 *S = C / Max;
398 }
399 else
400 *H = *S = 0;
401 }
402
403
404 /**
405 * @brief Convert a Hue-Saturation-Value (HSV) color to sRGB
406 *
407 * @param R, G, B pointers to hold the result
408 * @param H, S, V the input HSV values
409 *
410 * The input values are assumed to be scaled as
411 * 0 <= H < 360,
412 * 0 <= S <= 1,
413 * 0 <= V <= 1.
414 * The output sRGB values are scaled between 0 and 1. This is the inverse
415 * transformation of Rgb2Hsv.
416 *
417 * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
418 */
Hsv2Rgb(double * R,double * G,double * B,double H,double S,double V)419 void Hsv2Rgb(double *R, double *G, double *B, double H, double S, double V)
420 {
421 double C = S * V;
422 double Min = V - C;
423 double X;
424
425
426 H -= 360*floor(H/360);
427 H /= 60;
428 X = C*(1 - fabs(H - 2*floor(H/2) - 1));
429
430 switch((int)H)
431 {
432 case 0:
433 *R = Min + C;
434 *G = Min + X;
435 *B = Min;
436 break;
437 case 1:
438 *R = Min + X;
439 *G = Min + C;
440 *B = Min;
441 break;
442 case 2:
443 *R = Min;
444 *G = Min + C;
445 *B = Min + X;
446 break;
447 case 3:
448 *R = Min;
449 *G = Min + X;
450 *B = Min + C;
451 break;
452 case 4:
453 *R = Min + X;
454 *G = Min;
455 *B = Min + C;
456 break;
457 case 5:
458 *R = Min + C;
459 *G = Min;
460 *B = Min + X;
461 break;
462 default:
463 *R = *G = *B = 0;
464 }
465 }
466
467
468 /**
469 * @brief Convert an sRGB color to Hue-Saturation-Lightness (HSL)
470 *
471 * @param H, S, L pointers to hold the result
472 * @param R, G, B the input sRGB values scaled in [0,1]
473 *
474 * This routine transforms from sRGB to the double hexcone HSL color space
475 * The sRGB values are assumed to be between 0 and 1. The outputs are
476 * H = hexagonal hue angle (0 <= H < 360),
477 * S = { C/(2L) if L <= 1/2 (0 <= S <= 1),
478 * { C/(2 - 2L) if L > 1/2
479 * L = (max(R',G',B') + min(R',G',B'))/2 (0 <= L <= 1),
480 * where C = max(R',G',B') - min(R',G',B'). The inverse color transformation
481 * is given by Hsl2Rgb.
482 *
483 * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
484 */
Rgb2Hsl(double * H,double * S,double * L,double R,double G,double B)485 void Rgb2Hsl(double *H, double *S, double *L, double R, double G, double B)
486 {
487 double Max = MAX3(R, G, B);
488 double Min = MIN3(R, G, B);
489 double C = Max - Min;
490
491
492 *L = (Max + Min)/2;
493
494 if(C > 0)
495 {
496 if(Max == R)
497 {
498 *H = (G - B) / C;
499
500 if(G < B)
501 *H += 6;
502 }
503 else if(Max == G)
504 *H = 2 + (B - R) / C;
505 else
506 *H = 4 + (R - G) / C;
507
508 *H *= 60;
509 *S = (*L <= 0.5) ? (C/(2*(*L))) : (C/(2 - 2*(*L)));
510 }
511 else
512 *H = *S = 0;
513 }
514
515
516 /**
517 * @brief Convert a Hue-Saturation-Lightness (HSL) color to sRGB
518 *
519 * @param R, G, B pointers to hold the result
520 * @param H, S, L the input HSL values
521 *
522 * The input values are assumed to be scaled as
523 * 0 <= H < 360,
524 * 0 <= S <= 1,
525 * 0 <= L <= 1.
526 * The output sRGB values are scaled between 0 and 1. This is the inverse
527 * transformation of Rgb2Hsl.
528 *
529 * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
530 */
Hsl2Rgb(double * R,double * G,double * B,double H,double S,double L)531 void Hsl2Rgb(double *R, double *G, double *B, double H, double S, double L)
532 {
533 double C = (L <= 0.5) ? (2*L*S) : ((2 - 2*L)*S);
534 double Min = L - 0.5*C;
535 double X;
536
537
538 H -= 360*floor(H/360);
539 H /= 60;
540 X = C*(1 - fabs(H - 2*floor(H/2) - 1));
541
542 switch((int)H)
543 {
544 case 0:
545 *R = Min + C;
546 *G = Min + X;
547 *B = Min;
548 break;
549 case 1:
550 *R = Min + X;
551 *G = Min + C;
552 *B = Min;
553 break;
554 case 2:
555 *R = Min;
556 *G = Min + C;
557 *B = Min + X;
558 break;
559 case 3:
560 *R = Min;
561 *G = Min + X;
562 *B = Min + C;
563 break;
564 case 4:
565 *R = Min + X;
566 *G = Min;
567 *B = Min + C;
568 break;
569 case 5:
570 *R = Min + C;
571 *G = Min;
572 *B = Min + X;
573 break;
574 default:
575 *R = *G = *B = 0;
576 }
577 }
578
579
580 /**
581 * @brief Convert an sRGB color to Hue-Saturation-Intensity (HSI)
582 *
583 * @param H, S, I pointers to hold the result
584 * @param R, G, B the input sRGB values scaled in [0,1]
585 *
586 * This routine transforms from sRGB to the cylindrical HSI color space. The
587 * sRGB values are assumed to be between 0 and 1. The output values are
588 * H = polar hue angle (0 <= H < 360),
589 * S = 1 - min(R',G',B')/I (0 <= S <= 1),
590 * I = (R'+G'+B')/3 (0 <= I <= 1).
591 * The inverse color transformation is given by Hsi2Rgb.
592 *
593 * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
594 */
Rgb2Hsi(double * H,double * S,double * I,double R,double G,double B)595 void Rgb2Hsi(double *H, double *S, double *I, double R, double G, double B)
596 {
597 double alpha = 0.5*(2*R - G - B);
598 double beta = 0.866025403784439*(G - B);
599
600
601 *I = (R + G + B)/3;
602
603 if(*I > 0)
604 {
605 *S = 1 - MIN3(R,G,B) / *I;
606 *H = atan2(beta, alpha)*(180/M_PI);
607
608 if(*H < 0)
609 *H += 360;
610 }
611 else
612 *H = *S = 0;
613 }
614
615
616 /**
617 * @brief Convert a Hue-Saturation-Intesity (HSI) color to sRGB
618 *
619 * @param R, G, B pointers to hold the result
620 * @param H, S, I the input HSI values
621 *
622 * The input values are assumed to be scaled as
623 * 0 <= H < 360,
624 * 0 <= S <= 1,
625 * 0 <= I <= 1.
626 * The output sRGB values are scaled between 0 and 1. This is the inverse
627 * transformation of Rgb2Hsi.
628 *
629 * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
630 */
Hsi2Rgb(double * R,double * G,double * B,double H,double S,double I)631 void Hsi2Rgb(double *R, double *G, double *B, double H, double S, double I)
632 {
633 H -= 360*floor(H/360);
634
635 if(H < 120)
636 {
637 *B = I*(1 - S);
638 *R = I*(1 + S*cos(H*(M_PI/180))/cos((60 - H)*(M_PI/180)));
639 *G = 3*I - *R - *B;
640 }
641 else if(H < 240)
642 {
643 H -= 120;
644 *R = I*(1 - S);
645 *G = I*(1 + S*cos(H*(M_PI/180))/cos((60 - H)*(M_PI/180)));
646 *B = 3*I - *R - *G;
647 }
648 else
649 {
650 H -= 240;
651 *G = I*(1 - S);
652 *B = I*(1 + S*cos(H*(M_PI/180))/cos((60 - H)*(M_PI/180)));
653 *R = 3*I - *G - *B;
654 }
655 }
656
657
658 /*
659 * == CIE color transformations ==
660 *
661 * The following routines implement transformations between sRGB and
662 * the CIE color spaces XYZ, L*a*b, L*u*v*, and L*C*H*. These
663 * transforms assume a 2 degree observer angle and a D65 illuminant.
664 */
665
666
667 /**
668 * @brief Transform sRGB to CIE XYZ with the D65 white point
669 *
670 * @param X, Y, Z pointers to hold the result
671 * @param R, G, B the input sRGB values
672 *
673 * Poynton, "Frequently Asked Questions About Color," page 10
674 * Wikipedia: http://en.wikipedia.org/wiki/SRGB
675 * Wikipedia: http://en.wikipedia.org/wiki/CIE_1931_color_space
676 */
Rgb2Xyz(double * X,double * Y,double * Z,double R,double G,double B)677 void Rgb2Xyz(double *X, double *Y, double *Z, double R, double G, double B)
678 {
679 R = INVGAMMACORRECTION(R);
680 G = INVGAMMACORRECTION(G);
681 B = INVGAMMACORRECTION(B);
682 *X = (double)(0.4123955889674142161*R + 0.3575834307637148171*G + 0.1804926473817015735*B);
683 *Y = (double)(0.2125862307855955516*R + 0.7151703037034108499*G + 0.07220049864333622685*B);
684 *Z = (double)(0.01929721549174694484*R + 0.1191838645808485318*G + 0.9504971251315797660*B);
685 }
686
687
688 /**
689 * @brief Transform CIE XYZ to sRGB with the D65 white point
690 *
691 * @param R, G, B pointers to hold the result
692 * @param X, Y, Z the input XYZ values
693 *
694 * Official sRGB specification (IEC 61966-2-1:1999)
695 * Poynton, "Frequently Asked Questions About Color," page 10
696 * Wikipedia: http://en.wikipedia.org/wiki/SRGB
697 * Wikipedia: http://en.wikipedia.org/wiki/CIE_1931_color_space
698 */
Xyz2Rgb(double * R,double * G,double * B,double X,double Y,double Z)699 void Xyz2Rgb(double *R, double *G, double *B, double X, double Y, double Z)
700 {
701 double R1, B1, G1, Min;
702
703
704 R1 = (double)( 3.2406*X - 1.5372*Y - 0.4986*Z);
705 G1 = (double)(-0.9689*X + 1.8758*Y + 0.0415*Z);
706 B1 = (double)( 0.0557*X - 0.2040*Y + 1.0570*Z);
707
708 Min = MIN3(R1, G1, B1);
709
710 /* Force nonnegative values so that gamma correction is well-defined. */
711 if(Min < 0)
712 {
713 R1 -= Min;
714 G1 -= Min;
715 B1 -= Min;
716 }
717
718 /* Transform from RGB to R'G'B' */
719 *R = GAMMACORRECTION(R1);
720 *G = GAMMACORRECTION(G1);
721 *B = GAMMACORRECTION(B1);
722 }
723
724
725 /**
726 * Convert CIE XYZ to CIE L*a*b* (CIELAB) with the D65 white point
727 *
728 * @param L, a, b pointers to hold the result
729 * @param X, Y, Z the input XYZ values
730 *
731 * Wikipedia: http://en.wikipedia.org/wiki/Lab_color_space
732 */
Xyz2Lab(double * L,double * a,double * b,double X,double Y,double Z)733 void Xyz2Lab(double *L, double *a, double *b, double X, double Y, double Z)
734 {
735 X /= WHITEPOINT_X;
736 Y /= WHITEPOINT_Y;
737 Z /= WHITEPOINT_Z;
738 X = LABF(X);
739 Y = LABF(Y);
740 Z = LABF(Z);
741 *L = 116*Y - 16;
742 *a = 500*(X - Y);
743 *b = 200*(Y - Z);
744 }
745
746
747 /**
748 * Convert CIE L*a*b* (CIELAB) to CIE XYZ with the D65 white point
749 *
750 * @param X, Y, Z pointers to hold the result
751 * @param L, a, b the input L*a*b* values
752 *
753 * Wikipedia: http://en.wikipedia.org/wiki/Lab_color_space
754 */
Lab2Xyz(double * X,double * Y,double * Z,double L,double a,double b)755 void Lab2Xyz(double *X, double *Y, double *Z, double L, double a, double b)
756 {
757 L = (L + 16)/116;
758 a = L + a/500;
759 b = L - b/200;
760 *X = WHITEPOINT_X*LABINVF(a);
761 *Y = WHITEPOINT_Y*LABINVF(L);
762 *Z = WHITEPOINT_Z*LABINVF(b);
763 }
764
765
766 /**
767 * Convert CIE XYZ to CIE L*u*v* (CIELUV) with the D65 white point
768 *
769 * @param L, u, v pointers to hold the result
770 * @param X, Y, Z the input XYZ values
771 *
772 * Wikipedia: http://en.wikipedia.org/wiki/CIELUV_color_space
773 */
Xyz2Luv(double * L,double * u,double * v,double X,double Y,double Z)774 void Xyz2Luv(double *L, double *u, double *v, double X, double Y, double Z)
775 {
776 double u1, v1, Denom;
777
778
779 if((Denom = X + 15*Y + 3*Z) > 0)
780 {
781 u1 = (4*X) / Denom;
782 v1 = (9*Y) / Denom;
783 }
784 else
785 u1 = v1 = 0;
786
787 Y /= WHITEPOINT_Y;
788 Y = LABF(Y);
789 *L = 116*Y - 16;
790 *u = 13*(*L)*(u1 - WHITEPOINT_U);
791 *v = 13*(*L)*(v1 - WHITEPOINT_V);
792 }
793
794
795 /**
796 * Convert CIE L*u*v* (CIELUV) to CIE XYZ with the D65 white point
797 *
798 * @param X, Y, Z pointers to hold the result
799 * @param L, u, v the input L*u*v* values
800 *
801 * Wikipedia: http://en.wikipedia.org/wiki/CIELUV_color_space
802 */
Luv2Xyz(double * X,double * Y,double * Z,double L,double u,double v)803 void Luv2Xyz(double *X, double *Y, double *Z, double L, double u, double v)
804 {
805 *Y = (L + 16)/116;
806 *Y = WHITEPOINT_Y*LABINVF(*Y);
807
808 if(L != 0)
809 {
810 u /= L;
811 v /= L;
812 }
813
814 u = u/13 + WHITEPOINT_U;
815 v = v/13 + WHITEPOINT_V;
816 *X = (*Y) * ((9*u)/(4*v));
817 *Z = (*Y) * ((3 - 0.75*u)/v - 5);
818 }
819
820
821 /**
822 * Convert CIE XYZ to CIE L*C*H* with the D65 white point
823 *
824 * @param L, C, H pointers to hold the result
825 * @param X, Y, Z the input XYZ values
826 *
827 * CIE L*C*H* is related to CIE L*a*b* by
828 * a* = C* cos(H* pi/180),
829 * b* = C* sin(H* pi/180).
830 */
Xyz2Lch(double * L,double * C,double * H,double X,double Y,double Z)831 void Xyz2Lch(double *L, double *C, double *H, double X, double Y, double Z)
832 {
833 double a, b;
834
835
836 Xyz2Lab(L, &a, &b, X, Y, Z);
837 *C = sqrt(a*a + b*b);
838 *H = atan2(b, a)*180.0/M_PI;
839
840 if(*H < 0)
841 *H += 360;
842 }
843
844 /**
845 * Convert CIE L*C*H* to CIE XYZ with the D65 white point
846 *
847 * @param X, Y, Z pointers to hold the result
848 * @param L, C, H the input L*C*H* values
849 */
Lch2Xyz(double * X,double * Y,double * Z,double L,double C,double H)850 void Lch2Xyz(double *X, double *Y, double *Z, double L, double C, double H)
851 {
852 double a = C * cos(H*(M_PI/180.0));
853 double b = C * sin(H*(M_PI/180.0));
854
855
856 Lab2Xyz(X, Y, Z, L, a, b);
857 }
858
859
860 /** @brief XYZ to CAT02 LMS */
Xyz2Cat02lms(double * L,double * M,double * S,double X,double Y,double Z)861 void Xyz2Cat02lms(double *L, double *M, double *S, double X, double Y, double Z)
862 {
863 *L = (double)( 0.7328*X + 0.4296*Y - 0.1624*Z);
864 *M = (double)(-0.7036*X + 1.6975*Y + 0.0061*Z);
865 *S = (double)( 0.0030*X + 0.0136*Y + 0.9834*Z);
866 }
867
868
869 /** @brief CAT02 LMS to XYZ */
Cat02lms2Xyz(double * X,double * Y,double * Z,double L,double M,double S)870 void Cat02lms2Xyz(double *X, double *Y, double *Z, double L, double M, double S)
871 {
872 *X = (double)( 1.096123820835514*L - 0.278869000218287*M + 0.182745179382773*S);
873 *Y = (double)( 0.454369041975359*L + 0.473533154307412*M + 0.072097803717229*S);
874 *Z = (double)(-0.009627608738429*L - 0.005698031216113*M + 1.015325639954543*S);
875 }
876
877
878 /*
879 * == Glue functions for multi-stage transforms ==
880 */
881
Rgb2Lab(double * L,double * a,double * b,double R,double G,double B)882 void Rgb2Lab(double *L, double *a, double *b, double R, double G, double B)
883 {
884 double X, Y, Z;
885 Rgb2Xyz(&X, &Y, &Z, R, G, B);
886 Xyz2Lab(L, a, b, X, Y, Z);
887 }
888
889
Lab2Rgb(double * R,double * G,double * B,double L,double a,double b)890 void Lab2Rgb(double *R, double *G, double *B, double L, double a, double b)
891 {
892 double X, Y, Z;
893 Lab2Xyz(&X, &Y, &Z, L, a, b);
894 Xyz2Rgb(R, G, B, X, Y, Z);
895 }
896
897
Rgb2Luv(double * L,double * u,double * v,double R,double G,double B)898 void Rgb2Luv(double *L, double *u, double *v, double R, double G, double B)
899 {
900 double X, Y, Z;
901 Rgb2Xyz(&X, &Y, &Z, R, G, B);
902 Xyz2Luv(L, u, v, X, Y, Z);
903 }
904
905
Luv2Rgb(double * R,double * G,double * B,double L,double u,double v)906 void Luv2Rgb(double *R, double *G, double *B, double L, double u, double v)
907 {
908 double X, Y, Z;
909 Luv2Xyz(&X, &Y, &Z, L, u, v);
910 Xyz2Rgb(R, G, B, X, Y, Z);
911 }
912
Rgb2Lch(double * L,double * C,double * H,double R,double G,double B)913 void Rgb2Lch(double *L, double *C, double *H, double R, double G, double B)
914 {
915 double X, Y, Z;
916 Rgb2Xyz(&X, &Y, &Z, R, G, B);
917 Xyz2Lch(L, C, H, X, Y, Z);
918 }
919
920
Lch2Rgb(double * R,double * G,double * B,double L,double C,double H)921 void Lch2Rgb(double *R, double *G, double *B, double L, double C, double H)
922 {
923 double X, Y, Z;
924 Lch2Xyz(&X, &Y, &Z, L, C, H);
925 Xyz2Rgb(R, G, B, X, Y, Z);
926 }
927
928
Rgb2Cat02lms(double * L,double * M,double * S,double R,double G,double B)929 void Rgb2Cat02lms(double *L, double *M, double *S, double R, double G, double B)
930 {
931 double X, Y, Z;
932 Rgb2Xyz(&X, &Y, &Z, R, G, B);
933 Xyz2Cat02lms(L, M, S, X, Y, Z);
934 }
935
936
Cat02lms2Rgb(double * R,double * G,double * B,double L,double M,double S)937 void Cat02lms2Rgb(double *R, double *G, double *B, double L, double M, double S)
938 {
939 double X, Y, Z;
940 Cat02lms2Xyz(&X, &Y, &Z, L, M, S);
941 Xyz2Rgb(R, G, B, X, Y, Z);
942 }
943
944 } /* namespace */
945