libs/canvas/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 "canvas/colorspace.h"
 139
 140 /** @brief Min of A and B */
 141 #define MIN(A,B)        (((A) <= (B)) ? (A) : (B))
 142
 143 /** @brief Max of A and B */
 144 #define MAX(A,B)        (((A) >= (B)) ? (A) : (B))
 145
 146 /** @brief Min of A, B, and C */
 147 #define MIN3(A,B,C)     (((A) <= (B)) ? MIN(A,C) : MIN(B,C))
 148
 149 /** @brief Max of A, B, and C */
 150 #define MAX3(A,B,C)     (((A) >= (B)) ? MAX(A,C) : MAX(B,C))
 151
 152 #ifndef M_PI
 153 /** @brief The constant pi */
 154 #define M_PI    3.14159265358979323846264338327950288
 155 #endif
 156
 157 /**
 158  * @brief sRGB gamma correction, transforms R to R'
 159  * http://en.wikipedia.org/wiki/SRGB
 160  */
 161 #define GAMMACORRECTION(t)      \
 162         (((t) <= 0.0031306684425005883) ? \
 163         (12.92*(t)) : (1.055*pow((t), 0.416666666666666667) - 0.055))
 164
 165 /**
 166  * @brief Inverse sRGB gamma correction, transforms R' to R
 167  */
 168 #define INVGAMMACORRECTION(t)   \
 169         (((t) <= 0.0404482362771076) ? \
 170         ((t)/12.92) : pow(((t) + 0.055)/1.055, 2.4))
 171
 172 /**
 173  * @brief CIE L*a*b* f function (used to convert XYZ to L*a*b*)
 174  * http://en.wikipedia.org/wiki/Lab_color_space
 175  */
 176 #define LABF(t) \
 177         ((t >= 8.85645167903563082e-3) ? \
 178         pow(t,0.333333333333333) : (841.0/108.0)*(t) + (4.0/29.0))
 179
 180 /**
 181  * @brief CIE L*a*b* inverse f function
 182  * http://en.wikipedia.org/wiki/Lab_color_space
 183  */
 184 #define LABINVF(t)      \
 185         ((t >= 0.206896551724137931) ? \
 186         ((t)*(t)*(t)) : (108.0/841.0)*((t) - (4.0/29.0)))
 187
 188 /** @brief u'v' coordinates of the white point for CIE Lu*v* */
 189 #define WHITEPOINT_U    ((4*WHITEPOINT_X) \
 190         /(WHITEPOINT_X + 15*WHITEPOINT_Y + 3*WHITEPOINT_Z))
 191 #define WHITEPOINT_V    ((9*WHITEPOINT_Y) \
 192         /(WHITEPOINT_X + 15*WHITEPOINT_Y + 3*WHITEPOINT_Z))
 193
 194 /** @brief Enumeration of the supported color spaces */
 195 #define UNKNOWN_SPACE   0
 196 #define RGB_SPACE               1
 197 #define YUV_SPACE               2
 198 #define YCBCR_SPACE             3
 199 #define JPEGYCBCR_SPACE 4
 200 #define YPBPR_SPACE             5
 201 #define YDBDR_SPACE             6
 202 #define YIQ_SPACE               7
 203 #define HSV_SPACE               8
 204 #define HSL_SPACE               9
 205 #define HSI_SPACE               10
 206 #define XYZ_SPACE               11
 207 #define LAB_SPACE               12
 208 #define LUV_SPACE               13
 209 #define LCH_SPACE               14
 210 #define CAT02LMS_SPACE  15
 211
 212 #define NUM_TRANSFORM_PAIRS             18
 213
 214
 215
 216 /*
 217  * == Linear color transformations ==
 218  *
 219  * The following routines implement transformations between sRGB and
 220  * the linearly-related color spaces Y'UV, Y'PbPr, Y'DbDr, and Y'IQ.
 221  */
 222
 223
 224 /**
 225  * @brief Convert sRGB to NTSC/PAL Y'UV Luma + Chroma
 226  *
 227  * @param Y, U, V pointers to hold the result
 228  * @param R, G, B the input sRGB values
 229  *
 230  * Wikipedia: http://en.wikipedia.org/wiki/YUV
 231  */
 232 void Rgb2Yuv(double *Y, double *U, double *V, double R, double G, double B)
 233 {
 234         *Y = (double)( 0.299*R + 0.587*G + 0.114*B);
 235         *U = (double)(-0.147*R - 0.289*G + 0.436*B);
 236         *V = (double)( 0.615*R - 0.515*G - 0.100*B);
 237 }
 238
 239
 240 /**
 241  * @brief Convert NTSC/PAL Y'UV to sRGB
 242  *
 243  * @param R, G, B pointers to hold the result
 244  * @param Y, U, V the input YUV values
 245  */
 246 void Yuv2Rgb(double *R, double *G, double *B, double Y, double U, double V)
 247 {
 248         *R = (double)(Y - 3.945707070708279e-05*U + 1.1398279671717170825*V);
 249         *G = (double)(Y - 0.3946101641414141437*U - 0.5805003156565656797*V);
 250         *B = (double)(Y + 2.0319996843434342537*U - 4.813762626262513e-04*V);
 251 }
 252
 253
 254 /** @brief sRGB to Y'CbCr Luma + Chroma */
 255 void Rgb2Ycbcr(double *Y, double *Cb, double *Cr, double R, double G, double B)
 256 {
 257         *Y  = (double)( 65.481*R + 128.553*G +  24.966*B +  16);
 258         *Cb = (double)(-37.797*R -  74.203*G + 112.0  *B + 128);
 259         *Cr = (double)(112.0  *R -  93.786*G -  18.214*B + 128);
 260 }
 261
 262
 263 /** @brief Y'CbCr to sRGB */
 264 void Ycbcr2Rgb(double *R, double *G, double *B, double Y, double Cr, double Cb)
 265 {
 266         Y -= 16;
 267         Cb -= 128;
 268         Cr -= 128;
 269         *R = (double)(0.00456621004566210107*Y + 1.1808799897946415e-09*Cr + 0.00625892896994393634*Cb);
 270         *G = (double)(0.00456621004566210107*Y - 0.00153632368604490212*Cr - 0.00318811094965570701*Cb);
 271         *B = (double)(0.00456621004566210107*Y + 0.00791071623355474145*Cr + 1.1977497040190077e-08*Cb);
 272 }
 273
 274
 275 /** @brief sRGB to JPEG-Y'CbCr Luma + Chroma */
 276 void Rgb2Jpegycbcr(double *Y, double *Cb, double *Cr, double R, double G, double B)
 277 {
 278         Rgb2Ypbpr(Y, Cb, Cr, R, G, B);
 279         *Cb += (double)0.5;
 280         *Cr += (double)0.5;
 281 }
 282
 283 /** @brief JPEG-Y'CbCr to sRGB */
 284 void Jpegycbcr2Rgb(double *R, double *G, double *B, double Y, double Cb, double Cr)
 285 {
 286         Cb -= (double)0.5;
 287         Cr -= (double)0.5;
 288         Ypbpr2Rgb(R, G, B, Y, Cb, Cr);
 289 }
 290
 291
 292 /** @brief sRGB to Y'PbPr Luma (ITU-R BT.601) + Chroma */
 293 void Rgb2Ypbpr(double *Y, double *Pb, double *Pr, double R, double G, double B)
 294 {
 295         *Y  = (double)( 0.299    *R + 0.587   *G + 0.114   *B);
 296         *Pb = (double)(-0.1687367*R - 0.331264*G + 0.5     *B);
 297         *Pr = (double)( 0.5      *R - 0.418688*G - 0.081312*B);
 298 }
 299
 300
 301 /** @brief Y'PbPr to sRGB */
 302 void Ypbpr2Rgb(double *R, double *G, double *B, double Y, double Pb, double Pr)
 303 {
 304         *R = (double)(0.99999999999914679361*Y - 1.2188941887145875e-06*Pb + 1.4019995886561440468*Pr);
 305         *G = (double)(0.99999975910502514331*Y - 0.34413567816504303521*Pb - 0.71413649331646789076*Pr);
 306         *B = (double)(1.00000124040004623180*Y + 1.77200006607230409200*Pb + 2.1453384174593273e-06*Pr);
 307 }
 308
 309
 310 /** @brief sRGB to SECAM Y'DbDr Luma + Chroma */
 311 void Rgb2Ydbdr(double *Y, double *Db, double *Dr, double R, double G, double B)
 312 {
 313         *Y  = (double)( 0.299*R + 0.587*G + 0.114*B);
 314         *Db = (double)(-0.450*R - 0.883*G + 1.333*B);
 315         *Dr = (double)(-1.333*R + 1.116*G + 0.217*B);
 316 }
 317
 318
 319 /** @brief SECAM Y'DbDr to sRGB */
 320 void Ydbdr2Rgb(double *R, double *G, double *B, double Y, double Db, double Dr)
 321 {
 322         *R = (double)(Y + 9.2303716147657e-05*Db - 0.52591263066186533*Dr);
 323         *G = (double)(Y - 0.12913289889050927*Db + 0.26789932820759876*Dr);
 324         *B = (double)(Y + 0.66467905997895482*Db - 7.9202543533108e-05*Dr);
 325 }
 326
 327
 328 /** @brief sRGB to NTSC YIQ */
 329 void Rgb2Yiq(double *Y, double *I, double *Q, double R, double G, double B)
 330 {
 331         *Y = (double)(0.299   *R + 0.587   *G + 0.114   *B);
 332         *I = (double)(0.595716*R - 0.274453*G - 0.321263*B);
 333         *Q = (double)(0.211456*R - 0.522591*G + 0.311135*B);
 334 }
 335
 336
 337 /** @brief Convert NTSC YIQ to sRGB */
 338 void Yiq2Rgb(double *R, double *G, double *B, double Y, double I, double Q)
 339 {
 340         *R = (double)(Y + 0.9562957197589482261*I + 0.6210244164652610754*Q);
 341         *G = (double)(Y - 0.2721220993185104464*I - 0.6473805968256950427*Q);
 342         *B = (double)(Y - 1.1069890167364901945*I + 1.7046149983646481374*Q);
 343 }
 344
 345
 346
 347 /*
 348  * == Hue Saturation Value/Lightness/Intensity color transformations ==
 349  *
 350  * The following routines implement transformations between sRGB and
 351  * color spaces HSV, HSL, and HSI.
 352  */
 353
 354
 355 /**
 356  * @brief Convert an sRGB color to Hue-Saturation-Value (HSV)
 357  *
 358  * @param H, S, V pointers to hold the result
 359  * @param R, G, B the input sRGB values scaled in [0,1]
 360  *
 361  * This routine transforms from sRGB to the hexcone HSV color space.  The
 362  * sRGB values are assumed to be between 0 and 1.  The output values are
 363  *   H = hexagonal hue angle   (0 <= H < 360),
 364  *   S = C/V                   (0 <= S <= 1),
 365  *   V = max(R',G',B')         (0 <= V <= 1),
 366  * where C = max(R',G',B') - min(R',G',B').  The inverse color transformation
 367  * is given by Hsv2Rgb.
 368  *
 369  * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
 370  */
 371 void Rgb2Hsv(double *H, double *S, double *V, double R, double G, double B)
 372 {
 373         double Max = MAX3(R, G, B);
 374         double Min = MIN3(R, G, B);
 375         double C = Max - Min;
 376
 377
 378         *V = Max;
 379
 380         if(C > 0)
 381         {
 382                 if(Max == R)
 383                 {
 384                         *H = (G - B) / C;
 385
 386                         if(G < B)
 387                                 *H += 6;
 388                 }
 389                 else if(Max == G)
 390                         *H = 2 + (B - R) / C;
 391                 else
 392                         *H = 4 + (R - G) / C;
 393
 394                 *H *= 60;
 395                 *S = C / Max;
 396         }
 397         else
 398                 *H = *S = 0;
 399 }
 400
 401
 402 /**
 403  * @brief Convert a Hue-Saturation-Value (HSV) color to sRGB
 404  *
 405  * @param R, G, B pointers to hold the result
 406  * @param H, S, V the input HSV values
 407  *
 408  * The input values are assumed to be scaled as
 409  *    0 <= H < 360,
 410  *    0 <= S <= 1,
 411  *    0 <= V <= 1.
 412  * The output sRGB values are scaled between 0 and 1.  This is the inverse
 413  * transformation of Rgb2Hsv.
 414  *
 415  * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
 416  */
 417 void Hsv2Rgb(double *R, double *G, double *B, double H, double S, double V)
 418 {
 419         double C = S * V;
 420         double Min = V - C;
 421         double X;
 422
 423
 424         H -= 360*floor(H/360);
 425         H /= 60;
 426         X = C*(1 - fabs(H - 2*floor(H/2) - 1));
 427
 428         switch((int)H)
 429         {
 430         case 0:
 431                 *R = Min + C;
 432                 *G = Min + X;
 433                 *B = Min;
 434                 break;
 435         case 1:
 436                 *R = Min + X;
 437                 *G = Min + C;
 438                 *B = Min;
 439                 break;
 440         case 2:
 441                 *R = Min;
 442                 *G = Min + C;
 443                 *B = Min + X;
 444                 break;
 445         case 3:
 446                 *R = Min;
 447                 *G = Min + X;
 448                 *B = Min + C;
 449                 break;
 450         case 4:
 451                 *R = Min + X;
 452                 *G = Min;
 453                 *B = Min + C;
 454                 break;
 455         case 5:
 456                 *R = Min + C;
 457                 *G = Min;
 458                 *B = Min + X;
 459                 break;
 460         default:
 461                 *R = *G = *B = 0;
 462         }
 463 }
 464
 465
 466 /**
 467  * @brief Convert an sRGB color to Hue-Saturation-Lightness (HSL)
 468  *
 469  * @param H, S, L pointers to hold the result
 470  * @param R, G, B the input sRGB values scaled in [0,1]
 471  *
 472  * This routine transforms from sRGB to the double hexcone HSL color space
 473  * The sRGB values are assumed to be between 0 and 1.  The outputs are
 474  *   H = hexagonal hue angle                (0 <= H < 360),
 475  *   S = { C/(2L)     if L <= 1/2           (0 <= S <= 1),
 476  *       { C/(2 - 2L) if L >  1/2
 477  *   L = (max(R',G',B') + min(R',G',B'))/2  (0 <= L <= 1),
 478  * where C = max(R',G',B') - min(R',G',B').  The inverse color transformation
 479  * is given by Hsl2Rgb.
 480  *
 481  * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
 482  */
 483 void Rgb2Hsl(double *H, double *S, double *L, double R, double G, double B)
 484 {
 485         double Max = MAX3(R, G, B);
 486         double Min = MIN3(R, G, B);
 487         double C = Max - Min;
 488
 489
 490         *L = (Max + Min)/2;
 491
 492         if(C > 0)
 493         {
 494                 if(Max == R)
 495                 {
 496                         *H = (G - B) / C;
 497
 498                         if(G < B)
 499                                 *H += 6;
 500                 }
 501                 else if(Max == G)
 502                         *H = 2 + (B - R) / C;
 503                 else
 504                         *H = 4 + (R - G) / C;
 505
 506                 *H *= 60;
 507                 *S = (*L <= 0.5) ? (C/(2*(*L))) : (C/(2 - 2*(*L)));
 508         }
 509         else
 510                 *H = *S = 0;
 511 }
 512
 513
 514 /**
 515  * @brief Convert a Hue-Saturation-Lightness (HSL) color to sRGB
 516  *
 517  * @param R, G, B pointers to hold the result
 518  * @param H, S, L the input HSL values
 519  *
 520  * The input values are assumed to be scaled as
 521  *    0 <= H < 360,
 522  *    0 <= S <= 1,
 523  *    0 <= L <= 1.
 524  * The output sRGB values are scaled between 0 and 1.  This is the inverse
 525  * transformation of Rgb2Hsl.
 526  *
 527  * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
 528  */
 529 void Hsl2Rgb(double *R, double *G, double *B, double H, double S, double L)
 530 {
 531         double C = (L <= 0.5) ? (2*L*S) : ((2 - 2*L)*S);
 532         double Min = L - 0.5*C;
 533         double X;
 534
 535
 536         H -= 360*floor(H/360);
 537         H /= 60;
 538         X = C*(1 - fabs(H - 2*floor(H/2) - 1));
 539
 540         switch((int)H)
 541         {
 542         case 0:
 543                 *R = Min + C;
 544                 *G = Min + X;
 545                 *B = Min;
 546                 break;
 547         case 1:
 548                 *R = Min + X;
 549                 *G = Min + C;
 550                 *B = Min;
 551                 break;
 552         case 2:
 553                 *R = Min;
 554                 *G = Min + C;
 555                 *B = Min + X;
 556                 break;
 557         case 3:
 558                 *R = Min;
 559                 *G = Min + X;
 560                 *B = Min + C;
 561                 break;
 562         case 4:
 563                 *R = Min + X;
 564                 *G = Min;
 565                 *B = Min + C;
 566                 break;
 567         case 5:
 568                 *R = Min + C;
 569                 *G = Min;
 570                 *B = Min + X;
 571                 break;
 572         default:
 573                 *R = *G = *B = 0;
 574         }
 575 }
 576
 577
 578 /**
 579  * @brief Convert an sRGB color to Hue-Saturation-Intensity (HSI)
 580  *
 581  * @param H, S, I pointers to hold the result
 582  * @param R, G, B the input sRGB values scaled in [0,1]
 583  *
 584  * This routine transforms from sRGB to the cylindrical HSI color space.  The
 585  * sRGB values are assumed to be between 0 and 1.  The output values are
 586  *   H = polar hue angle         (0 <= H < 360),
 587  *   S = 1 - min(R',G',B')/I     (0 <= S <= 1),
 588  *   I = (R'+G'+B')/3            (0 <= I <= 1).
 589  * The inverse color transformation is given by Hsi2Rgb.
 590  *
 591  * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
 592  */
 593 void Rgb2Hsi(double *H, double *S, double *I, double R, double G, double B)
 594 {
 595         double alpha = 0.5*(2*R - G - B);
 596         double beta = 0.866025403784439*(G - B);
 597
 598
 599         *I = (R + G + B)/3;
 600
 601         if(*I > 0)
 602         {
 603                 *S = 1 - MIN3(R,G,B) / *I;
 604                 *H = atan2(beta, alpha)*(180/M_PI);
 605
 606                 if(*H < 0)
 607                         *H += 360;
 608         }
 609         else
 610                 *H = *S = 0;
 611 }
 612
 613
 614 /**
 615  * @brief Convert a Hue-Saturation-Intesity (HSI) color to sRGB
 616  *
 617  * @param R, G, B pointers to hold the result
 618  * @param H, S, I the input HSI values
 619  *
 620  * The input values are assumed to be scaled as
 621  *    0 <= H < 360,
 622  *    0 <= S <= 1,
 623  *    0 <= I <= 1.
 624  * The output sRGB values are scaled between 0 and 1.  This is the inverse
 625  * transformation of Rgb2Hsi.
 626  *
 627  * Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV
 628  */
 629 void Hsi2Rgb(double *R, double *G, double *B, double H, double S, double I)
 630 {
 631         H -= 360*floor(H/360);
 632
 633         if(H < 120)
 634         {
 635                 *B = I*(1 - S);
 636                 *R = I*(1 + S*cos(H*(M_PI/180))/cos((60 - H)*(M_PI/180)));
 637                 *G = 3*I - *R - *B;
 638         }
 639         else if(H < 240)
 640         {
 641                 H -= 120;
 642                 *R = I*(1 - S);
 643                 *G = I*(1 + S*cos(H*(M_PI/180))/cos((60 - H)*(M_PI/180)));
 644                 *B = 3*I - *R - *G;
 645         }
 646         else
 647         {
 648                 H -= 240;
 649                 *G = I*(1 - S);
 650                 *B = I*(1 + S*cos(H*(M_PI/180))/cos((60 - H)*(M_PI/180)));
 651                 *R = 3*I - *G - *B;
 652         }
 653 }
 654
 655
 656 /*
 657  * == CIE color transformations ==
 658  *
 659  * The following routines implement transformations between sRGB and
 660  * the CIE color spaces XYZ, L*a*b, L*u*v*, and L*C*H*.  These
 661  * transforms assume a 2 degree observer angle and a D65 illuminant.
 662  */
 663
 664
 665 /**
 666  * @brief Transform sRGB to CIE XYZ with the D65 white point
 667  *
 668  * @param X, Y, Z pointers to hold the result
 669  * @param R, G, B the input sRGB values
 670  *
 671  * Poynton, "Frequently Asked Questions About Color," page 10
 672  * Wikipedia: http://en.wikipedia.org/wiki/SRGB
 673  * Wikipedia: http://en.wikipedia.org/wiki/CIE_1931_color_space
 674  */
 675 void Rgb2Xyz(double *X, double *Y, double *Z, double R, double G, double B)
 676 {
 677         R = INVGAMMACORRECTION(R);
 678         G = INVGAMMACORRECTION(G);
 679         B = INVGAMMACORRECTION(B);
 680         *X = (double)(0.4123955889674142161*R + 0.3575834307637148171*G + 0.1804926473817015735*B);
 681         *Y = (double)(0.2125862307855955516*R + 0.7151703037034108499*G + 0.07220049864333622685*B);
 682         *Z = (double)(0.01929721549174694484*R + 0.1191838645808485318*G + 0.9504971251315797660*B);
 683 }
 684
 685
 686 /**
 687  * @brief Transform CIE XYZ to sRGB with the D65 white point
 688  *
 689  * @param R, G, B pointers to hold the result
 690  * @param X, Y, Z the input XYZ values
 691  *
 692  * Official sRGB specification (IEC 61966-2-1:1999)
 693  * Poynton, "Frequently Asked Questions About Color," page 10
 694  * Wikipedia: http://en.wikipedia.org/wiki/SRGB
 695  * Wikipedia: http://en.wikipedia.org/wiki/CIE_1931_color_space
 696  */
 697 void Xyz2Rgb(double *R, double *G, double *B, double X, double Y, double Z)
 698 {
 699         double R1, B1, G1, Min;
 700
 701
 702         R1 = (double)( 3.2406*X - 1.5372*Y - 0.4986*Z);
 703         G1 = (double)(-0.9689*X + 1.8758*Y + 0.0415*Z);
 704         B1 = (double)( 0.0557*X - 0.2040*Y + 1.0570*Z);
 705
 706         Min = MIN3(R1, G1, B1);
 707
 708         /* Force nonnegative values so that gamma correction is well-defined. */
 709         if(Min < 0)
 710         {
 711                 R1 -= Min;
 712                 G1 -= Min;
 713                 B1 -= Min;
 714         }
 715
 716         /* Transform from RGB to R'G'B' */
 717         *R = GAMMACORRECTION(R1);
 718         *G = GAMMACORRECTION(G1);
 719         *B = GAMMACORRECTION(B1);
 720 }
 721
 722
 723 /**
 724  * Convert CIE XYZ to CIE L*a*b* (CIELAB) with the D65 white point
 725  *
 726  * @param L, a, b pointers to hold the result
 727  * @param X, Y, Z the input XYZ values
 728  *
 729  * Wikipedia: http://en.wikipedia.org/wiki/Lab_color_space
 730  */
 731 void Xyz2Lab(double *L, double *a, double *b, double X, double Y, double Z)
 732 {
 733         X /= WHITEPOINT_X;
 734         Y /= WHITEPOINT_Y;
 735         Z /= WHITEPOINT_Z;
 736         X = LABF(X);
 737         Y = LABF(Y);
 738         Z = LABF(Z);
 739         *L = 116*Y - 16;
 740         *a = 500*(X - Y);
 741         *b = 200*(Y - Z);
 742 }
 743
 744
 745 /**
 746  * Convert CIE L*a*b* (CIELAB) to CIE XYZ with the D65 white point
 747  *
 748  * @param X, Y, Z pointers to hold the result
 749  * @param L, a, b the input L*a*b* values
 750  *
 751  * Wikipedia: http://en.wikipedia.org/wiki/Lab_color_space
 752  */
 753 void Lab2Xyz(double *X, double *Y, double *Z, double L, double a, double b)
 754 {
 755         L = (L + 16)/116;
 756         a = L + a/500;
 757         b = L - b/200;
 758         *X = WHITEPOINT_X*LABINVF(a);
 759         *Y = WHITEPOINT_Y*LABINVF(L);
 760         *Z = WHITEPOINT_Z*LABINVF(b);
 761 }
 762
 763
 764 /**
 765  * Convert CIE XYZ to CIE L*u*v* (CIELUV) with the D65 white point
 766  *
 767  * @param L, u, v pointers to hold the result
 768  * @param X, Y, Z the input XYZ values
 769  *
 770  * Wikipedia: http://en.wikipedia.org/wiki/CIELUV_color_space
 771  */
 772 void Xyz2Luv(double *L, double *u, double *v, double X, double Y, double Z)
 773 {
 774         double u1, v1, Denom;
 775
 776
 777         if((Denom = X + 15*Y + 3*Z) > 0)
 778         {
 779                 u1 = (4*X) / Denom;
 780                 v1 = (9*Y) / Denom;
 781         }
 782         else
 783                 u1 = v1 = 0;
 784
 785         Y /= WHITEPOINT_Y;
 786         Y = LABF(Y);
 787         *L = 116*Y - 16;
 788         *u = 13*(*L)*(u1 - WHITEPOINT_U);
 789         *v = 13*(*L)*(v1 - WHITEPOINT_V);
 790 }
 791
 792
 793 /**
 794  * Convert CIE L*u*v* (CIELUV) to CIE XYZ with the D65 white point
 795  *
 796  * @param X, Y, Z pointers to hold the result
 797  * @param L, u, v the input L*u*v* values
 798  *
 799  * Wikipedia: http://en.wikipedia.org/wiki/CIELUV_color_space
 800  */
 801 void Luv2Xyz(double *X, double *Y, double *Z, double L, double u, double v)
 802 {
 803         *Y = (L + 16)/116;
 804         *Y = WHITEPOINT_Y*LABINVF(*Y);
 805
 806         if(L != 0)
 807         {
 808                 u /= L;
 809                 v /= L;
 810         }
 811
 812         u = u/13 + WHITEPOINT_U;
 813         v = v/13 + WHITEPOINT_V;
 814         *X = (*Y) * ((9*u)/(4*v));
 815         *Z = (*Y) * ((3 - 0.75*u)/v - 5);
 816 }
 817
 818
 819 /**
 820  * Convert CIE XYZ to CIE L*C*H* with the D65 white point
 821  *
 822  * @param L, C, H pointers to hold the result
 823  * @param X, Y, Z the input XYZ values
 824  *
 825  * CIE L*C*H* is related to CIE L*a*b* by
 826  *    a* = C* cos(H* pi/180),
 827  *    b* = C* sin(H* pi/180).
 828  */
 829 void Xyz2Lch(double *L, double *C, double *H, double X, double Y, double Z)
 830 {
 831         double a, b;
 832
 833
 834         Xyz2Lab(L, &a, &b, X, Y, Z);
 835         *C = sqrt(a*a + b*b);
 836         *H = atan2(b, a)*180.0/M_PI;
 837
 838         if(*H < 0)
 839                 *H += 360;
 840 }
 841
 842 /**
 843  * Convert CIE L*C*H* to CIE XYZ with the D65 white point
 844  *
 845  * @param X, Y, Z pointers to hold the result
 846  * @param L, C, H the input L*C*H* values
 847  */
 848 void Lch2Xyz(double *X, double *Y, double *Z, double L, double C, double H)
 849 {
 850         double a = C * cos(H*(M_PI/180.0));
 851         double b = C * sin(H*(M_PI/180.0));
 852
 853
 854         Lab2Xyz(X, Y, Z, L, a, b);
 855 }
 856
 857
 858 /** @brief XYZ to CAT02 LMS */
 859 void Xyz2Cat02lms(double *L, double *M, double *S, double X, double Y, double Z)
 860 {
 861         *L = (double)( 0.7328*X + 0.4296*Y - 0.1624*Z);
 862         *M = (double)(-0.7036*X + 1.6975*Y + 0.0061*Z);
 863         *S = (double)( 0.0030*X + 0.0136*Y + 0.9834*Z);
 864 }
 865
 866
 867 /** @brief CAT02 LMS to XYZ */
 868 void Cat02lms2Xyz(double *X, double *Y, double *Z, double L, double M, double S)
 869 {
 870         *X = (double)( 1.096123820835514*L - 0.278869000218287*M + 0.182745179382773*S);
 871         *Y = (double)( 0.454369041975359*L + 0.473533154307412*M + 0.072097803717229*S);
 872         *Z = (double)(-0.009627608738429*L - 0.005698031216113*M + 1.015325639954543*S);
 873 }
 874
 875
 876 /*
 877  * == Glue functions for multi-stage transforms ==
 878  */
 879
 880 void Rgb2Lab(double *L, double *a, double *b, double R, double G, double B)
 881 {
 882         double X, Y, Z;
 883         Rgb2Xyz(&X, &Y, &Z, R, G, B);
 884         Xyz2Lab(L, a, b, X, Y, Z);
 885 }
 886
 887
 888 void Lab2Rgb(double *R, double *G, double *B, double L, double a, double b)
 889 {
 890         double X, Y, Z;
 891         Lab2Xyz(&X, &Y, &Z, L, a, b);
 892         Xyz2Rgb(R, G, B, X, Y, Z);
 893 }
 894
 895
 896 void Rgb2Luv(double *L, double *u, double *v, double R, double G, double B)
 897 {
 898         double X, Y, Z;
 899         Rgb2Xyz(&X, &Y, &Z, R, G, B);
 900         Xyz2Luv(L, u, v, X, Y, Z);
 901 }
 902
 903
 904 void Luv2Rgb(double *R, double *G, double *B, double L, double u, double v)
 905 {
 906         double X, Y, Z;
 907         Luv2Xyz(&X, &Y, &Z, L, u, v);
 908         Xyz2Rgb(R, G, B, X, Y, Z);
 909 }
 910
 911 void Rgb2Lch(double *L, double *C, double *H, double R, double G, double B)
 912 {
 913         double X, Y, Z;
 914         Rgb2Xyz(&X, &Y, &Z, R, G, B);
 915         Xyz2Lch(L, C, H, X, Y, Z);
 916 }
 917
 918
 919 void Lch2Rgb(double *R, double *G, double *B, double L, double C, double H)
 920 {
 921         double X, Y, Z;
 922         Lch2Xyz(&X, &Y, &Z, L, C, H);
 923         Xyz2Rgb(R, G, B, X, Y, Z);
 924 }
 925
 926
 927 void Rgb2Cat02lms(double *L, double *M, double *S, double R, double G, double B)
 928 {
 929         double X, Y, Z;
 930         Rgb2Xyz(&X, &Y, &Z, R, G, B);
 931         Xyz2Cat02lms(L, M, S, X, Y, Z);
 932 }
 933
 934
 935 void Cat02lms2Rgb(double *R, double *G, double *B, double L, double M, double S)
 936 {
 937         double X, Y, Z;
 938         Cat02lms2Xyz(&X, &Y, &Z, L, M, S);
 939         Xyz2Rgb(R, G, B, X, Y, Z);
 940 }