libs/gtkmm2ext/colorspace.cc

   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  */
 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  */
 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 */
 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 */
 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 */
 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 */
 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 */
 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 */
 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 */
 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 */
 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 */
 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 */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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 */
 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 */
 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
 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
 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
 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
 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
 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
 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
 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
 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 */