images/color.c

   1 /*
   2  *      Image Library -- Color Spaces
   3  *
   4  *      (c) 2006 Pavel Charvat <pchar@ucw.cz>
   5  *
   6  *      This software may be freely distributed and used according to the terms
   7  *      of the GNU Lesser General Public License.
   8  */
   9
  10 #undef LOCAL_DEBUG
  11
  12 #include "lib/lib.h"
  13 #include "images/images.h"
  14 #include "images/color.h"
  15 #include "images/error.h"
  16 #include "images/math.h"
  17
  18 #include <string.h>
  19 #include <math.h>
  20
  21 uns color_space_channels[COLOR_SPACE_MAX] = {
  22   [COLOR_SPACE_UNKNOWN] = 0,
  23   [COLOR_SPACE_UNKNOWN_1] = 1,
  24   [COLOR_SPACE_UNKNOWN_2] = 2,
  25   [COLOR_SPACE_UNKNOWN_3] = 3,
  26   [COLOR_SPACE_UNKNOWN_4] = 4,
  27   [COLOR_SPACE_GRAYSCALE] = 1,
  28   [COLOR_SPACE_RGB] = 3,
  29   [COLOR_SPACE_XYZ] = 3,
  30   [COLOR_SPACE_LAB] = 3,
  31   [COLOR_SPACE_YCBCR] = 3,
  32   [COLOR_SPACE_CMYK] = 4,
  33   [COLOR_SPACE_YCCK] = 4,
  34 };
  35
  36 byte *color_space_name[COLOR_SPACE_MAX] = {
  37   [COLOR_SPACE_UNKNOWN] = "Unknown",
  38   [COLOR_SPACE_UNKNOWN_1] = "1-channel",
  39   [COLOR_SPACE_UNKNOWN_2] = "2-channels",
  40   [COLOR_SPACE_UNKNOWN_3] = "3-channels",
  41   [COLOR_SPACE_UNKNOWN_4] = "4-channels",
  42   [COLOR_SPACE_GRAYSCALE] = "Grayscale",
  43   [COLOR_SPACE_RGB] = "RGB",
  44   [COLOR_SPACE_XYZ] = "XYZ",
  45   [COLOR_SPACE_LAB] = "LAB",
  46   [COLOR_SPACE_YCBCR] = "YCbCr",
  47   [COLOR_SPACE_CMYK] = "CMYK",
  48   [COLOR_SPACE_YCCK] = "YCCK",
  49 };
  50
  51 byte *
  52 color_space_id_to_name(uns id)
  53 {
  54   ASSERT(id < COLOR_SPACE_MAX);
  55   return color_space_name[id];
  56 }
  57
  58 uns
  59 color_space_name_to_id(byte *name)
  60 {
  61   for (uns i = 1; i < COLOR_SPACE_MAX; i++)
  62     if (color_space_name[i] && !strcasecmp(name, color_space_name[i]))
  63       return i;
  64   return 0;
  65 }
  66
  67 struct color color_black = { .color_space = COLOR_SPACE_GRAYSCALE };
  68 struct color color_white = { .c = { 255 }, .color_space = COLOR_SPACE_GRAYSCALE };
  69
  70 int
  71 color_get(struct color *color, byte *src, uns src_space)
  72 {
  73   color->color_space = src_space;
  74   memcpy(color->c, src, color_space_channels[src_space]);
  75   return 1;
  76 }
  77
  78 int
  79 color_put(struct image_context *ctx, struct color *color, byte *dest, uns dest_space)
  80 {
  81   switch (dest_space)
  82     {
  83       case COLOR_SPACE_GRAYSCALE:
  84         switch (color->color_space)
  85           {
  86             case COLOR_SPACE_GRAYSCALE:
  87               dest[0] = color->c[0];
  88               return 1;
  89             case COLOR_SPACE_RGB:
  90               dest[0] = rgb_to_gray_func(color->c[0], color->c[1], color->c[2]);
  91               return 1;
  92           }
  93         break;
  94       case COLOR_SPACE_RGB:
  95         switch (color->color_space)
  96           {
  97             case COLOR_SPACE_GRAYSCALE:
  98               dest[0] = dest[1] = dest[2] = color->c[0];
  99               return 1;
 100             case COLOR_SPACE_RGB:
 101               dest[0] = color->c[0];
 102               dest[1] = color->c[1];
 103               dest[2] = color->c[2];
 104               return 1;
 105             case COLOR_SPACE_CMYK:
 106               {
 107                 double rgb[3], cmyk[4];
 108                 for (uns i = 0; i < 4; i++)
 109                   cmyk[i] = color->c[i] * (1.0 / 255);
 110                 cmyk_to_rgb_exact(rgb, cmyk);
 111                 for (uns i = 0; i < 3; i++)
 112                   dest[i] = CLAMP(rgb[i] * 255, 0, 255);
 113               }
 114               return 1;
 115           }
 116         break;
 117       case COLOR_SPACE_CMYK:
 118         switch (color->color_space)
 119           {
 120             case COLOR_SPACE_GRAYSCALE:
 121               dest[0] = dest[1] = dest[2] = 0;
 122               dest[3] = 255 - color->c[0];
 123               return 1;
 124             case COLOR_SPACE_RGB:
 125               {
 126                 double rgb[3], cmyk[4];
 127                 for (uns i = 0; i < 3; i++)
 128                   rgb[i] = color->c[i] * (1.0 / 255);
 129                 rgb_to_cmyk_exact(cmyk, rgb);
 130                 for (uns i = 0; i < 4; i++)
 131                   dest[i] = CLAMP(cmyk[i] * 255, 0, 255);
 132               }
 133               return 1;
 134           }
 135         break;
 136     }
 137   if (dest_space != COLOR_SPACE_RGB )
 138     {
 139       /* Try to convert the color via RGB */
 140       struct color rgb;
 141       if (!color_put(ctx, color, rgb.c, COLOR_SPACE_RGB))
 142         return 0;
 143       rgb.color_space = COLOR_SPACE_RGB;
 144       return color_put(ctx, &rgb, dest, dest_space);
 145     }
 146   IMAGE_ERROR(ctx, IMAGE_ERROR_INVALID_PIXEL_FORMAT, "Conversion from %s to %s is not supported",
 147       color_space_id_to_name(color->color_space), color_space_id_to_name(color->color_space));
 148   return 0;
 149 }
 150
 151
 152 /********************* IMAGE CONVERSION ROUTINES **********************/
 153
 154 struct image_conv_options image_conv_defaults = {
 155   .flags = IMAGE_CONV_COPY_ALPHA | IMAGE_CONV_FILL_ALPHA | IMAGE_CONV_APPLY_ALPHA,
 156   .background = { .color_space = COLOR_SPACE_GRAYSCALE } };
 157
 158 /* Grayscale <-> RGB */
 159
 160 #define IMAGE_WALK_PREFIX(x) walk_##x
 161 #define IMAGE_WALK_FUNC_NAME image_conv_gray_1_to_rgb_n
 162 #define IMAGE_WALK_DOUBLE
 163 #define IMAGE_WALK_SEC_COL_STEP 1
 164 #define IMAGE_WALK_UNROLL 4
 165 #define IMAGE_WALK_DO_STEP do{ walk_pos[0] = walk_pos[1] = walk_pos[2] = walk_sec_pos[0]; }while(0)
 166 #include "images/image-walk.h"
 167
 168 #define IMAGE_WALK_PREFIX(x) walk_##x
 169 #define IMAGE_WALK_FUNC_NAME image_conv_rgb_n_to_gray_1
 170 #define IMAGE_WALK_DOUBLE
 171 #define IMAGE_WALK_COL_STEP 1
 172 #define IMAGE_WALK_UNROLL 2
 173 #define IMAGE_WALK_DO_STEP do{ walk_pos[0] = rgb_to_gray_func(walk_sec_pos[0], walk_sec_pos[1], walk_sec_pos[2]); }while(0)
 174 #include "images/image-walk.h"
 175
 176 /* YCbCr <-> RGB */
 177
 178 static inline void
 179 pixel_conv_ycbcr_to_rgb(byte *dest, byte *src)
 180 {
 181   /* R = Y                + 1.40200 * Cr
 182    * G = Y - 0.34414 * Cb - 0.71414 * Cr
 183    * B = Y + 1.77200 * Cb */
 184   int y = src[0], cb = src[1] - 128, cr = src[2] - 128;
 185   dest[0] = CLAMP(y + (91881 * cr) / 0x10000, 0, 255);
 186   dest[1] = CLAMP(y - (22553 * cb + 46801 * cr) / 0x10000, 0, 255);
 187   dest[2] = CLAMP(y + (116129 * cb) / 0x10000, 0, 255);
 188 }
 189
 190 #define IMAGE_WALK_PREFIX(x) walk_##x
 191 #define IMAGE_WALK_FUNC_NAME image_conv_ycbcr_n_to_rgb_n
 192 #define IMAGE_WALK_DOUBLE
 193 #define IMAGE_WALK_DO_STEP do{ pixel_conv_ycbcr_to_rgb(walk_pos, walk_sec_pos); }while(0)
 194 #include "images/image-walk.h"
 195
 196 static inline void
 197 pixel_conv_rgb_to_ycbcr(byte *dest, byte *src)
 198 {
 199   /* Y  =  0.29900 * R + 0.58700 * G + 0.11400 * B
 200    * Cb = -0.16874 * R - 0.33126 * G + 0.50000 * B + CENTER
 201    * Cr =  0.50000 * R - 0.41869 * G - 0.08131 * B + CENTER */
 202   uns r = src[0], g = src[1], b = src[2];
 203   dest[0] = (19595 * r + 38470 * g + 7471 * b) / 0x10000;
 204   dest[1] = (0x800000 + 0x8000 * b - 11058 * r - 21710 * g) / 0x10000;
 205   dest[2] = (0x800000 + 0x8000 * r - 27439 * g - 5329 * b) / 0x10000;
 206 }
 207
 208 #define IMAGE_WALK_PREFIX(x) walk_##x
 209 #define IMAGE_WALK_FUNC_NAME image_conv_rgb_n_to_ycbcr_n
 210 #define IMAGE_WALK_DOUBLE
 211 #define IMAGE_WALK_DO_STEP do{ pixel_conv_rgb_to_ycbcr(walk_pos, walk_sec_pos); }while(0)
 212 #include "images/image-walk.h"
 213
 214 /* CMYK <-> RGB */
 215
 216 static inline void
 217 pixel_conv_cmyk_to_rgb(byte *dest, byte *src)
 218 {
 219   uns d = (255 - src[3]) * (0xffffffffU / 255 /255);
 220   dest[0] = d * (255 - src[0]) >> 24;
 221   dest[1] = d * (255 - src[1]) >> 24;
 222   dest[2] = d * (255 - src[2]) >> 24;
 223 }
 224
 225 #define IMAGE_WALK_PREFIX(x) walk_##x
 226 #define IMAGE_WALK_FUNC_NAME image_conv_cmyk_4_to_rgb_n
 227 #define IMAGE_WALK_DOUBLE
 228 #define IMAGE_WALK_SEC_COL_STEP 4
 229 #define IMAGE_WALK_DO_STEP do{ pixel_conv_cmyk_to_rgb(walk_pos, walk_sec_pos); }while(0)
 230 #include "images/image-walk.h"
 231
 232 static inline void
 233 pixel_conv_rgb_to_cmyk(byte *dest, byte *src)
 234 {
 235   uns k = MAX(src[0], src[1]);
 236   k = MAX(k, src[2]);
 237   uns d = fast_div_u32_u8(0x7fffffffU, k); /* == 0 for zero K */
 238   dest[0] = (d * (k - src[0])) >> 23;
 239   dest[1] = (d * (k - src[1])) >> 23;
 240   dest[2] = (d * (k - src[2])) >> 23;
 241   dest[3] = 255 - k;
 242 }
 243
 244 #define IMAGE_WALK_PREFIX(x) walk_##x
 245 #define IMAGE_WALK_FUNC_NAME image_conv_rgb_n_to_cmyk_4
 246 #define IMAGE_WALK_DOUBLE
 247 #define IMAGE_WALK_COL_STEP 4
 248 #define IMAGE_WALK_DO_STEP do{ pixel_conv_rgb_to_cmyk(walk_pos, walk_sec_pos); }while(0)
 249 #include "images/image-walk.h"
 250
 251 /* YCCK <-> RGB */
 252
 253 static inline void
 254 pixel_conv_ycck_to_rgb(byte *dest, byte *src)
 255 {
 256   int y = src[0], cb = src[1] - 128, cr = src[2] - 128;
 257   uns d = (255 - src[3]) * (0xffffffffU / 255 /255);
 258   dest[0] = (d * CLAMP(y + (91881 * cr) / 0x10000, 0, 255) >> 24);
 259   dest[1] = (d * CLAMP(y - (22553 * cb + 46801 * cr) / 0x10000, 0, 255) >> 24);
 260   dest[2] = (d * CLAMP(y + (116129 * cb) / 0x10000, 0, 255) >> 24);
 261 }
 262
 263 #define IMAGE_WALK_PREFIX(x) walk_##x
 264 #define IMAGE_WALK_FUNC_NAME image_conv_ycck_4_to_rgb_n
 265 #define IMAGE_WALK_DOUBLE
 266 #define IMAGE_WALK_SEC_COL_STEP 4
 267 #define IMAGE_WALK_DO_STEP do{ pixel_conv_ycck_to_rgb(walk_pos, walk_sec_pos); }while(0)
 268 #include "images/image-walk.h"
 269
 270 static inline void
 271 pixel_conv_rgb_to_ycck(byte *dest, byte *src)
 272 {
 273   uns k = MAX(src[0], src[1]);
 274   k = MAX(k, src[2]);
 275   uns d = fast_div_u32_u8(0x7fffffffU, k); /* == 0 for zero K */
 276   uns r = 255 - ((d * (k - src[0])) >> 23);
 277   uns g = 255 - ((d * (k - src[1])) >> 23);
 278   uns b = 255 - ((d * (k - src[2])) >> 23);
 279   dest[0] = (19595 * r + 38470 * g + 7471 * b) / 0x10000;
 280   dest[1] = (0x800000 + 0x8000 * b - 11058 * r - 21710 * g) / 0x10000;
 281   dest[2] = (0x800000 + 0x8000 * r - 27439 * g - 5329 * b) / 0x10000;
 282   dest[3] = 255 - k;
 283 }
 284
 285 #define IMAGE_WALK_PREFIX(x) walk_##x
 286 #define IMAGE_WALK_FUNC_NAME image_conv_rgb_n_to_ycck_4
 287 #define IMAGE_WALK_DOUBLE
 288 #define IMAGE_WALK_COL_STEP 4
 289 #define IMAGE_WALK_DO_STEP do{ pixel_conv_rgb_to_ycck(walk_pos, walk_sec_pos); }while(0)
 290 #include "images/image-walk.h"
 291
 292 /* Main functions */
 293
 294 static int
 295 image_conv_color_space(struct image_context *ctx UNUSED, struct image *dest, struct image *src, struct image_conv_options *opt UNUSED)
 296 {
 297   switch (dest->flags & IMAGE_COLOR_SPACE)
 298     {
 299       case COLOR_SPACE_GRAYSCALE:
 300         switch (src->flags & IMAGE_COLOR_SPACE)
 301           {
 302             case COLOR_SPACE_RGB:
 303               if (dest->pixel_size == 1)
 304                 {
 305                   image_conv_rgb_n_to_gray_1(dest, src);
 306                   return 1;
 307                 }
 308               break;
 309           }
 310         break;
 311       case COLOR_SPACE_RGB:
 312         switch (src->flags & IMAGE_CHANNELS_FORMAT)
 313           {
 314             case COLOR_SPACE_GRAYSCALE:
 315               if (src->pixel_size == 1)
 316                 {
 317                   image_conv_gray_1_to_rgb_n(dest, src);
 318                   return 1;
 319                 }
 320               break;
 321             case COLOR_SPACE_YCBCR:
 322               image_conv_ycbcr_n_to_rgb_n(dest, src);
 323               return 1;
 324             case COLOR_SPACE_CMYK:
 325               if (src->pixel_size == 4)
 326                 {
 327                   image_conv_cmyk_4_to_rgb_n(dest, src);
 328                   return 1;
 329                 }
 330               break;
 331             case COLOR_SPACE_YCCK:
 332               if (src->pixel_size == 4)
 333                 {
 334                   image_conv_ycck_4_to_rgb_n(dest, src);
 335                   return 1;
 336                 }
 337               break;
 338         }
 339         break;
 340       case COLOR_SPACE_YCBCR:
 341         switch (src->flags & IMAGE_CHANNELS_FORMAT)
 342           {
 343             case COLOR_SPACE_RGB:
 344               image_conv_rgb_n_to_ycbcr_n(dest, src);
 345               return 1;
 346           }
 347         break;
 348       case COLOR_SPACE_CMYK:
 349         switch (src->flags & IMAGE_CHANNELS_FORMAT)
 350           {
 351             case COLOR_SPACE_RGB:
 352               if (dest->pixel_size == 4)
 353                 {
 354                   image_conv_rgb_n_to_cmyk_4(dest, src);
 355                   return 1;
 356                 }
 357               break;
 358           }
 359         break;
 360       case COLOR_SPACE_YCCK:
 361         switch (src->flags & IMAGE_CHANNELS_FORMAT)
 362           {
 363             case COLOR_SPACE_RGB:
 364               if (dest->pixel_size == 4)
 365                 {
 366                   image_conv_rgb_n_to_ycck_4(dest, src);
 367                   return 1;
 368                 }
 369               break;
 370           }
 371         break;
 372     }
 373   return 0;
 374 }
 375
 376 static void
 377 image_conv_copy(struct image *dest, struct image *src)
 378 {
 379   if (dest->pixels == src->pixels)
 380     return;
 381   else if (dest->pixel_size != src->pixel_size)
 382     {
 383       uns channels = MIN(dest->channels, src->channels);
 384       switch (channels)
 385         {
 386           case 1:
 387             {
 388 #             define IMAGE_WALK_PREFIX(x) walk_##x
 389 #             define IMAGE_WALK_INLINE
 390 #             define IMAGE_WALK_DOUBLE
 391 #             define IMAGE_WALK_IMAGE dest
 392 #             define IMAGE_WALK_SEC_IMAGE src
 393 #             define IMAGE_WALK_UNROLL 4
 394 #             define IMAGE_WALK_DO_STEP do{ walk_pos[0] = walk_sec_pos[0]; }while(0)
 395 #             include "images/image-walk.h"
 396             }
 397             return;
 398           case 2:
 399 #             define IMAGE_WALK_PREFIX(x) walk_##x
 400 #             define IMAGE_WALK_INLINE
 401 #             define IMAGE_WALK_DOUBLE
 402 #             define IMAGE_WALK_IMAGE dest
 403 #             define IMAGE_WALK_SEC_IMAGE src
 404 #             define IMAGE_WALK_UNROLL 4
 405 #             define IMAGE_WALK_DO_STEP do{ walk_pos[0] = walk_sec_pos[0]; walk_pos[1] = walk_sec_pos[1]; }while(0)
 406 #             include "images/image-walk.h"
 407             return;
 408           case 3:
 409 #             define IMAGE_WALK_PREFIX(x) walk_##x
 410 #             define IMAGE_WALK_INLINE
 411 #             define IMAGE_WALK_DOUBLE
 412 #             define IMAGE_WALK_IMAGE dest
 413 #             define IMAGE_WALK_SEC_IMAGE src
 414 #             define IMAGE_WALK_UNROLL 2
 415 #             define IMAGE_WALK_DO_STEP do{ walk_pos[0] = walk_sec_pos[0]; walk_pos[1] = walk_sec_pos[1]; walk_pos[2] = walk_sec_pos[2]; }while(0)
 416 #             include "images/image-walk.h"
 417             return;
 418           case 4:
 419 #             define IMAGE_WALK_PREFIX(x) walk_##x
 420 #             define IMAGE_WALK_INLINE
 421 #             define IMAGE_WALK_DOUBLE
 422 #             define IMAGE_WALK_IMAGE dest
 423 #             define IMAGE_WALK_SEC_IMAGE src
 424 #             define IMAGE_WALK_UNROLL 2
 425 #             define IMAGE_WALK_DO_STEP do{ walk_pos[0] = walk_sec_pos[0]; walk_pos[1] = walk_sec_pos[1]; walk_pos[2] = walk_sec_pos[2]; walk_pos[3] = walk_sec_pos[3]; }while(0)
 426 #             include "images/image-walk.h"
 427             return;
 428           default:
 429 #             define IMAGE_WALK_PREFIX(x) walk_##x
 430 #             define IMAGE_WALK_INLINE
 431 #             define IMAGE_WALK_DOUBLE
 432 #             define IMAGE_WALK_IMAGE dest
 433 #             define IMAGE_WALK_SEC_IMAGE src
 434 #             define IMAGE_WALK_DO_STEP do{ for (uns i = 0; i < channels; i++) walk_pos[i] = walk_sec_pos[i]; }while(0)
 435 #             include "images/image-walk.h"
 436             return;
 437         }
 438     }
 439   else if (dest->row_size != src->row_size || ((dest->flags | src->flags) & IMAGE_GAPS_PROTECTED))
 440     {
 441       byte *s = src->pixels;
 442       byte *d = dest->pixels;
 443       for (uns row = src->rows; row--; )
 444         {
 445           memcpy(d, s, src->row_pixels_size);
 446           d += dest->row_size;
 447           s += src->row_size;
 448         }
 449     }
 450   else if (dest->pixels != src->pixels)
 451     memcpy(dest->pixels, src->pixels, src->image_size);
 452 }
 453
 454 static void
 455 image_conv_fill_alpha(struct image *dest)
 456 {
 457   switch (dest->channels)
 458     {
 459       case 2:
 460         if (dest->pixel_size == 2)
 461           {
 462 #           define IMAGE_WALK_PREFIX(x) walk_##x
 463 #           define IMAGE_WALK_INLINE
 464 #           define IMAGE_WALK_IMAGE dest
 465 #           define IMAGE_WALK_COL_STEP 2
 466 #           define IMAGE_WALK_UNROLL 4
 467 #           define IMAGE_WALK_DO_STEP do{ walk_pos[1] = 255; }while(0)
 468 #           include "images/image-walk.h"
 469             return;
 470           }
 471         break;
 472       case 4:
 473         if (dest->pixel_size == 4)
 474           {
 475 #           define IMAGE_WALK_PREFIX(x) walk_##x
 476 #           define IMAGE_WALK_INLINE
 477 #           define IMAGE_WALK_IMAGE dest
 478 #           define IMAGE_WALK_COL_STEP 4
 479 #           define IMAGE_WALK_UNROLL 4
 480 #           define IMAGE_WALK_DO_STEP do{ walk_pos[3] = 255; }while(0)
 481 #           include "images/image-walk.h"
 482             return;
 483           }
 484         break;
 485     }
 486   {
 487 #   define IMAGE_WALK_PREFIX(x) walk_##x
 488 #   define IMAGE_WALK_INLINE
 489 #   define IMAGE_WALK_IMAGE dest
 490 #   define IMAGE_WALK_UNROLL 4
 491 #   define IMAGE_WALK_DO_STEP do{ walk_pos[dest->channels - 1] = 255; }while(0)
 492 #   include "images/image-walk.h"
 493   }
 494 }
 495
 496 static void
 497 image_conv_copy_alpha(struct image *dest, struct image *src)
 498 {
 499   if (dest->pixels != src->pixels || dest->channels != src->channels)
 500     {
 501 #     define IMAGE_WALK_PREFIX(x) walk_##x
 502 #     define IMAGE_WALK_INLINE
 503 #     define IMAGE_WALK_DOUBLE
 504 #     define IMAGE_WALK_IMAGE dest
 505 #     define IMAGE_WALK_SEC_IMAGE src
 506 #     define IMAGE_WALK_UNROLL 4
 507 #     define IMAGE_WALK_DO_STEP do{ walk_pos[dest->channels - 1] = walk_sec_pos[src->channels - 1]; }while(0)
 508 #     include "images/image-walk.h"
 509     }
 510 }
 511
 512 static inline uns
 513 image_conv_alpha_func(uns value, uns alpha, uns acoef, uns bcoef)
 514 {
 515   return ((uns)(acoef + (int)alpha * (int)(value - bcoef)) * (0xffffffffU / 255 / 255)) >> 24;
 516 }
 517
 518 static int
 519 image_conv_apply_alpha_from(struct image_context *ctx, struct image *dest, struct image *src, struct image_conv_options *opt)
 520 {
 521   if (!opt->background.color_space)
 522     return 1;
 523   byte background[IMAGE_MAX_CHANNELS];
 524   if (unlikely(!color_put(ctx, &opt->background, background, dest->flags & IMAGE_COLOR_SPACE)))
 525     return 0;
 526   uns a[IMAGE_MAX_CHANNELS], b[IMAGE_MAX_CHANNELS];
 527   for (uns i = 0; i < dest->channels; i++)
 528     a[i] = 255 * (b[i] = background[i]);
 529   switch (dest->channels)
 530     {
 531       case 1:
 532         {
 533 #         define IMAGE_WALK_PREFIX(x) walk_##x
 534 #         define IMAGE_WALK_INLINE
 535 #         define IMAGE_WALK_IMAGE dest
 536 #         define IMAGE_WALK_SEC_IMAGE src
 537 #         define IMAGE_WALK_DOUBLE
 538 #         define IMAGE_WALK_UNROLL 2
 539 #         define IMAGE_WALK_DO_STEP do{ \
 540               walk_pos[0] = image_conv_alpha_func(walk_pos[0], walk_sec_pos[src->channels - 1], a[0], b[0]); }while(0)
 541 #         include "images/image-walk.h"
 542         }
 543         return 1;
 544       case 3:
 545         {
 546 #         define IMAGE_WALK_PREFIX(x) walk_##x
 547 #         define IMAGE_WALK_INLINE
 548 #         define IMAGE_WALK_IMAGE dest
 549 #         define IMAGE_WALK_SEC_IMAGE src
 550 #         define IMAGE_WALK_DOUBLE
 551 #         define IMAGE_WALK_DO_STEP do{ \
 552               walk_pos[0] = image_conv_alpha_func(walk_pos[0], walk_sec_pos[src->channels - 1], a[0], b[0]); \
 553               walk_pos[1] = image_conv_alpha_func(walk_pos[1], walk_sec_pos[src->channels - 1], a[1], b[1]); \
 554               walk_pos[2] = image_conv_alpha_func(walk_pos[2], walk_sec_pos[src->channels - 1], a[2], b[2]); }while(0)
 555 #         include "images/image-walk.h"
 556         }
 557         return 1;
 558     }
 559   {
 560 #   define IMAGE_WALK_PREFIX(x) walk_##x
 561 #   define IMAGE_WALK_INLINE
 562 #   define IMAGE_WALK_IMAGE dest
 563 #   define IMAGE_WALK_SEC_IMAGE src
 564 #   define IMAGE_WALK_DOUBLE
 565 #   define IMAGE_WALK_DO_STEP do{ for (uns i = 0; i < dest->channels; i++) \
 566         walk_pos[i] = image_conv_alpha_func(walk_pos[i], walk_sec_pos[src->channels - 1], a[i], b[i]); }while(0)
 567 #   include "images/image-walk.h"
 568   }
 569   return 1;
 570 }
 571
 572 static int
 573 image_conv_apply_alpha_to(struct image_context *ctx, struct image *dest, struct image *src, struct image_conv_options *opt)
 574 {
 575   if (!opt->background.color_space)
 576     {
 577       image_conv_copy(dest, src);
 578       return 1;
 579     }
 580   byte background[IMAGE_MAX_CHANNELS];
 581   if (unlikely(!color_put(ctx, &opt->background, background, dest->flags & IMAGE_COLOR_SPACE)))
 582     return 0;
 583   uns a[IMAGE_MAX_CHANNELS], b[IMAGE_MAX_CHANNELS];
 584   for (uns i = 0; i < dest->channels; i++)
 585     a[i] = 255 * (b[i] = background[i]);
 586   switch (dest->channels)
 587     {
 588       case 1:
 589         {
 590 #         define IMAGE_WALK_PREFIX(x) walk_##x
 591 #         define IMAGE_WALK_INLINE
 592 #         define IMAGE_WALK_IMAGE dest
 593 #         define IMAGE_WALK_SEC_IMAGE src
 594 #         define IMAGE_WALK_DOUBLE
 595 #         define IMAGE_WALK_UNROLL 2
 596 #         define IMAGE_WALK_DO_STEP do{ \
 597               walk_pos[0] = image_conv_alpha_func(walk_sec_pos[0], walk_sec_pos[src->channels - 1], a[0], b[0]); }while(0)
 598 #         include "images/image-walk.h"
 599         }
 600         return 1;
 601       case 3:
 602         {
 603 #         define IMAGE_WALK_PREFIX(x) walk_##x
 604 #         define IMAGE_WALK_INLINE
 605 #         define IMAGE_WALK_IMAGE dest
 606 #         define IMAGE_WALK_SEC_IMAGE src
 607 #         define IMAGE_WALK_DOUBLE
 608 #         define IMAGE_WALK_DO_STEP do{ \
 609               walk_pos[0] = image_conv_alpha_func(walk_sec_pos[0], walk_sec_pos[src->channels - 1], a[0], b[0]); \
 610               walk_pos[1] = image_conv_alpha_func(walk_sec_pos[1], walk_sec_pos[src->channels - 1], a[1], b[1]); \
 611               walk_pos[2] = image_conv_alpha_func(walk_sec_pos[2], walk_sec_pos[src->channels - 1], a[2], b[2]); }while(0)
 612 #         include "images/image-walk.h"
 613         }
 614         return 1;
 615     }
 616   {
 617 #   define IMAGE_WALK_PREFIX(x) walk_##x
 618 #   define IMAGE_WALK_INLINE
 619 #   define IMAGE_WALK_IMAGE dest
 620 #   define IMAGE_WALK_SEC_IMAGE src
 621 #   define IMAGE_WALK_DOUBLE
 622 #   define IMAGE_WALK_DO_STEP do{ for (uns i = 0; i < dest->channels; i++) \
 623         walk_pos[i] = image_conv_alpha_func(walk_sec_pos[i], walk_sec_pos[src->channels - 1], a[i], b[i]); }while(0)
 624 #   include "images/image-walk.h"
 625   }
 626   return 1;
 627 }
 628
 629 int
 630 image_conv(struct image_context *ctx, struct image *dest, struct image *src, struct image_conv_options *opt)
 631 {
 632   ASSERT(dest->cols == src->cols && dest->rows == src->rows);
 633   if (!((dest->flags ^ src->flags) & IMAGE_COLOR_SPACE))
 634     {
 635       if (!(src->flags & IMAGE_ALPHA) || (dest->flags & IMAGE_ALPHA))
 636         image_conv_copy(dest, src);
 637       else if (unlikely(!image_conv_apply_alpha_to(ctx, dest, src, opt)))
 638         return 0;
 639     }
 640   else
 641     {
 642       if (!(src->flags & IMAGE_ALPHA))
 643         {
 644           if (unlikely(!image_conv_color_space(ctx, dest, src, opt)))
 645             goto error;
 646           if ((dest->flags & IMAGE_ALPHA) && (opt->flags & IMAGE_CONV_FILL_ALPHA))
 647             image_conv_fill_alpha(dest);
 648         }
 649       else
 650         {
 651           if (dest->flags & IMAGE_ALPHA)
 652             {
 653               if (dest->channels <= src->channels)
 654                 {
 655                   if (unlikely(!image_conv_color_space(ctx, dest, src, opt)))
 656                     goto error;
 657                   if (opt->flags & IMAGE_CONV_COPY_ALPHA)
 658                     image_conv_copy_alpha(dest, src);
 659                   else if (opt->flags & IMAGE_CONV_FILL_ALPHA)
 660                     image_conv_fill_alpha(dest);
 661                 }
 662               else
 663                 {
 664                   if (opt->flags & IMAGE_CONV_COPY_ALPHA)
 665                     image_conv_copy_alpha(dest, src);
 666                   else
 667                     image_conv_fill_alpha(dest);
 668                   if (unlikely(!image_conv_color_space(ctx, dest, src, opt)))
 669                     goto error;
 670                 }
 671             }
 672           else
 673             {
 674               if (dest->channels <= src->channels)
 675                 {
 676                   if (unlikely(!image_conv_color_space(ctx, dest, src, opt)))
 677                     goto error;
 678                   if (unlikely(!image_conv_apply_alpha_from(ctx, dest, src, opt)))
 679                     return 0;
 680                 }
 681               else
 682                 {
 683                   if (unlikely(!image_conv_apply_alpha_to(ctx, dest, src, opt)))
 684                     return 0;
 685                   if (unlikely(!image_conv_color_space(ctx, dest, dest, opt)))
 686                     goto error;
 687                 }
 688             }
 689         }
 690     }
 691   return 1;
 692 error:
 693   IMAGE_ERROR(ctx, IMAGE_ERROR_INVALID_PIXEL_FORMAT, "Image conversion not supported for such pixel formats");
 694   return 0;
 695 }
 696
 697 /********************* EXACT CONVERSION ROUTINES **********************/
 698
 699 /* Reference whites */
 700 #define COLOR_ILLUMINANT_A      0.44757, 0.40744
 701 #define COLOR_ILLUMINANT_B      0.34840, 0.35160
 702 #define COLOR_ILLUMINANT_C      0.31006, 0.31615
 703 #define COLOR_ILLUMINANT_D50    0.34567, 0.35850
 704 #define COLOR_ILLUMINANT_D55    0.33242, 0.34743
 705 #define COLOR_ILLUMINANT_D65    0.31273, 0.32902
 706 #define COLOR_ILLUMINANT_D75    0.29902, 0.31485
 707 #define COLOR_ILLUMINANT_9300K  0.28480, 0.29320
 708 #define COLOR_ILLUMINANT_E      (1./3.), (1./3.)
 709 #define COLOR_ILLUMINANT_F2     0.37207, 0.37512
 710 #define COLOR_ILLUMINANT_F7     0.31285, 0.32918
 711 #define COLOR_ILLUMINANT_F11    0.38054, 0.37691
 712
 713 const double
 714   color_illuminant_d50[2] = {COLOR_ILLUMINANT_D50},
 715   color_illuminant_d65[2] = {COLOR_ILLUMINANT_D65},
 716   color_illuminant_e[2] = {COLOR_ILLUMINANT_E};
 717
 718 /* RGB profiles (many missing) */
 719 const struct color_space_info
 720   color_adobe_rgb_info = {"Adobe RGB", {{0.6400, 0.3300}, {0.2100, 0.7100}, {0.1500, 0.0600}, {COLOR_ILLUMINANT_D65}}, {0.45, 0.45, 0, 0, 0}},
 721   color_apple_rgb_info = {"Apple RGB", {{0.6250, 0.3400}, {0.2800, 0.5950}, {0.1550, 0.0700}, {COLOR_ILLUMINANT_D65}}, {0.56, 0.56, 0, 0, 0}},
 722   color_cie_rgb_info = {"CIE RGB", {{0.7350, 0.2650}, {0.2740, 0.7170}, {0.1670, 0.0090}, {COLOR_ILLUMINANT_E}}, {0.45, 0.45, 0, 0, 0}},
 723   color_color_match_rgb_info = {"ColorMatch RGB", {{0.6300, 0.3400}, {0.2950, 0.6050}, {0.1500, 0.0750}, {COLOR_ILLUMINANT_D50}}, {0.56, 0.56, 0, 0, 0}},
 724   color_srgb_info = {"sRGB", {{0.6400, 0.3300}, {0.3000, 0.6000}, {0.1500, 0.0600}, {COLOR_ILLUMINANT_D65}}, {0.45, 0.42, 0.055, 0.003, 12.92}};
 725
 726 #define CLIP(x, min, max) (((x) < (min)) ? (min) : ((x) > (max)) ? (max) : (x))
 727
 728 static inline void
 729 clip(double a[3])
 730 {
 731   a[0] = CLIP(a[0], 0, 1);
 732   a[1] = CLIP(a[1], 0, 1);
 733   a[2] = CLIP(a[2], 0, 1);
 734 }
 735
 736 static inline void
 737 correct_gamma_simple(double dest[3], double src[3], const struct color_space_gamma_info *info)
 738 {
 739   dest[0] = pow(src[0], info->simple_gamma);
 740   dest[1] = pow(src[1], info->simple_gamma);
 741   dest[2] = pow(src[2], info->simple_gamma);
 742 }
 743
 744 static inline void
 745 invert_gamma_simple(double dest[3], double src[3], const struct color_space_gamma_info *info)
 746 {
 747   dest[0] = pow(src[0], 1 / info->simple_gamma);
 748   dest[1] = pow(src[1], 1 / info->simple_gamma);
 749   dest[2] = pow(src[2], 1 / info->simple_gamma);
 750 }
 751
 752 static inline void
 753 correct_gamma_detailed(double dest[3], double src[3], const struct color_space_gamma_info *info)
 754 {
 755   for (uns i = 0; i < 3; i++)
 756     if (src[i] > info->transition)
 757       dest[i] = (1 + info->offset) * pow(src[i], info->detailed_gamma) - info->offset;
 758     else
 759       dest[i] = info->slope * src[i];
 760 }
 761
 762 static inline void
 763 invert_gamma_detailed(double dest[3], double src[3], const struct color_space_gamma_info *info)
 764 {
 765   for (uns i = 0; i < 3; i++)
 766     if (src[i] > info->transition * info->slope)
 767       dest[i] = pow((src[i] + info->offset) / (1 + info->offset), 1 / info->detailed_gamma);
 768     else
 769       dest[i] = src[i] / info->slope;
 770 }
 771
 772 static inline void
 773 apply_matrix(double dest[3], double src[3], double matrix[9])
 774 {
 775   dest[0] = src[0] * matrix[0] + src[1] * matrix[1] + src[2] * matrix[2];
 776   dest[1] = src[0] * matrix[3] + src[1] * matrix[4] + src[2] * matrix[5];
 777   dest[2] = src[0] * matrix[6] + src[1] * matrix[7] + src[2] * matrix[8];
 778 }
 779
 780 void
 781 color_invert_matrix(double dest[9], double matrix[9])
 782 {
 783   double *i = dest, *m = matrix;
 784   double a0 = m[4] * m[8] - m[5] * m[7];
 785   double a1 = m[3] * m[8] - m[5] * m[6];
 786   double a2 = m[3] * m[7] - m[4] * m[6];
 787   double d = 1 / (m[0] * a0 - m[1] * a1 + m[2] * a2);
 788   i[0] = d * a0;
 789   i[3] = -d * a1;
 790   i[6] = d * a2;
 791   i[1] = -d * (m[1] * m[8] - m[2] * m[7]);
 792   i[4] = d * (m[0] * m[8] - m[2] * m[6]);
 793   i[7] = -d * (m[0] * m[7] - m[1] * m[6]);
 794   i[2] = d * (m[1] * m[5] - m[2] * m[4]);
 795   i[5] = -d * (m[0] * m[5] - m[2] * m[3]);
 796   i[8] = d * (m[0] * m[4] - m[1] * m[3]);
 797 }
 798
 799 static void
 800 mul_matrices(double r[9], double a[9], double b[9])
 801 {
 802   r[0] = a[0] * b[0] + a[1] * b[3] + a[2] * b[6];
 803   r[1] = a[0] * b[1] + a[1] * b[4] + a[2] * b[7];
 804   r[2] = a[0] * b[2] + a[1] * b[5] + a[2] * b[8];
 805   r[3] = a[3] * b[0] + a[4] * b[3] + a[5] * b[6];
 806   r[4] = a[3] * b[1] + a[4] * b[4] + a[5] * b[7];
 807   r[5] = a[3] * b[2] + a[4] * b[5] + a[5] * b[8];
 808   r[6] = a[6] * b[0] + a[7] * b[3] + a[8] * b[6];
 809   r[7] = a[6] * b[1] + a[7] * b[4] + a[8] * b[7];
 810   r[8] = a[6] * b[2] + a[7] * b[5] + a[8] * b[8];
 811 }
 812
 813 /* computes conversion matrix from a given color space to CIE XYZ */
 814 void
 815 color_compute_color_space_to_xyz_matrix(double matrix[9], const struct color_space_chromacity_info *space)
 816 {
 817   double wX = space->white[0] / space->white[1];
 818   double wZ = (1 - space->white[0] - space->white[1]) / space->white[1];
 819   double a[9], b[9];
 820   a[0] = space->prim1[0]; a[3] = space->prim1[1]; a[6] = 1 - a[0] - a[3];
 821   a[1] = space->prim2[0]; a[4] = space->prim2[1]; a[7] = 1 - a[1] - a[4];
 822   a[2] = space->prim3[0]; a[5] = space->prim3[1]; a[8] = 1 - a[2] - a[5];
 823   color_invert_matrix(b, a);
 824   double ra = wX * b[0] + b[1] + wZ * b[2];
 825   double rb = wX * b[3] + b[4] + wZ * b[5];
 826   double rc = wX * b[6] + b[7] + wZ * b[8];
 827   matrix[0] = a[0] * ra;
 828   matrix[1] = a[1] * rb;
 829   matrix[2] = a[2] * rc;
 830   matrix[3] = a[3] * ra;
 831   matrix[4] = a[4] * rb;
 832   matrix[5] = a[5] * rc;
 833   matrix[6] = a[6] * ra;
 834   matrix[7] = a[7] * rb;
 835   matrix[8] = a[8] * rc;
 836 }
 837
 838 /* computes matrix to join transformations with different reference whites */
 839 void
 840 color_compute_bradford_matrix(double matrix[9], const double source[2], const double dest[2])
 841 {
 842   /* cone response matrix and its inversion */
 843   static double r[9] = {
 844     0.8951, 0.2664, -0.1614,
 845     -0.7502, 1.7135, 0.0367,
 846     0.0389, -0.0685, 1.0296};
 847   //static double i[9] = {0.9870, -0.1471, 0.1600, 0.4323, 0.5184, 0.0493, -0.0085, 0.0400, 0.9685};
 848   double i[9];
 849   color_invert_matrix(i, r);
 850   double aX = source[0] / source[1];
 851   double aZ = (1 - source[0] - source[1]) / source[1];
 852   double bX = dest[0] / dest[1];
 853   double bZ = (1 - dest[0] - dest[1]) / dest[1];
 854   double x = (r[0] * bX + r[1] + r[2] * bZ) / (r[0] * aX + r[1] + r[2] * aZ);
 855   double y = (r[3] * bX + r[4] + r[5] * bZ) / (r[3] * aX + r[4] + r[5] * aZ);
 856   double z = (r[6] * bX + r[7] + r[8] * bZ) / (r[6] * aX + r[7] + r[8] * aZ);
 857   double m[9];
 858   m[0] = i[0] * x; m[1] = i[1] * y; m[2] = i[2] * z;
 859   m[3] = i[3] * x; m[4] = i[4] * y; m[5] = i[5] * z;
 860   m[6] = i[6] * x; m[7] = i[7] * y; m[8] = i[8] * z;
 861   mul_matrices(matrix, m, r);
 862 }
 863
 864 void
 865 color_compute_color_spaces_conversion_matrix(double matrix[9], const struct color_space_chromacity_info *src, const struct color_space_chromacity_info *dest)
 866 {
 867   double a_to_xyz[9], b_to_xyz[9], xyz_to_b[9], bradford[9], m[9];
 868   color_compute_color_space_to_xyz_matrix(a_to_xyz, src);
 869   color_compute_color_space_to_xyz_matrix(b_to_xyz, dest);
 870   color_invert_matrix(xyz_to_b, b_to_xyz);
 871   if (src->white[0] == dest->white[0] && src->white[1] == dest->white[1])
 872     mul_matrices(matrix, a_to_xyz, xyz_to_b);
 873   else
 874     {
 875       color_compute_bradford_matrix(bradford, src->white, dest->white);
 876       mul_matrices(m, a_to_xyz, bradford);
 877       mul_matrices(matrix, m, xyz_to_b);
 878     }
 879 }
 880
 881 /* sRGB to XYZ */
 882 void
 883 srgb_to_xyz_exact(double xyz[3], double srgb[3])
 884 {
 885   static double matrix[9] = {
 886     0.41248031, 0.35756952, 0.18043951,
 887     0.21268516, 0.71513904, 0.07217580,
 888     0.01933501, 0.11918984, 0.95031473};
 889   double srgb_lin[3];
 890   invert_gamma_detailed(srgb_lin, srgb, &color_srgb_info.gamma);
 891   apply_matrix(xyz, srgb_lin, matrix);
 892   xyz_to_srgb_exact(srgb_lin, xyz);
 893 }
 894
 895 /* XYZ to sRGB */
 896 void
 897 xyz_to_srgb_exact(double srgb[3], double xyz[3])
 898 {
 899   static double matrix[9] = {
 900      3.24026666, -1.53704957, -0.49850256,
 901     -0.96928381,  1.87604525,  0.04155678,
 902      0.05564281, -0.20402363,  1.05721334};
 903   double srgb_lin[3];
 904   apply_matrix(srgb_lin, xyz, matrix);
 905   clip(srgb_lin);
 906   correct_gamma_detailed(srgb, srgb_lin, &color_srgb_info.gamma);
 907 }
 908
 909 /* XYZ to CIE-Luv */
 910 void
 911 xyz_to_luv_exact(double luv[3], double xyz[3])
 912 {
 913   double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
 914   if (sum < 0.000001)
 915     luv[0] = luv[1] = luv[2] = 0;
 916   else
 917     {
 918       double var_u = 4 * xyz[0] / sum;
 919       double var_v = 9 * xyz[1] / sum;
 920       if (xyz[1] > 0.008856)
 921         luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
 922       else
 923         luv[0] = (116 * 7.787) * xyz[1];
 924      luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
 925      luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
 926      /* intervals [0..100], [-134..220], [-140..122] */
 927    }
 928 }
 929
 930 /* CIE-Luv to XYZ */
 931 void
 932 luv_to_xyz_exact(double xyz[3], double luv[3])
 933 {
 934   double var_u = luv[1] / (13 * luv[0]) + (4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z));
 935   double var_v = luv[2] / (13 * luv[0]) + (9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z));
 936   double var_y = (luv[0] + 16) / 116;
 937   double pow_y = var_y * var_y * var_y;
 938   if (pow_y > 0.008856)
 939     var_y = pow_y;
 940   else
 941     var_y = (var_y - 16 / 116) / 7.787;
 942   xyz[1] = var_y;
 943   xyz[0] = -(9 * xyz[1] * var_u) / ((var_u - 4) * var_v - var_u * var_v);
 944   xyz[2] = (9 * xyz[1] - 15 * var_v * xyz[1] - var_v * xyz[0]) / (3 * var_v);
 945 }
 946
 947 /* RGB to CMYK - a very simple version, not too accureate  */
 948 void
 949 rgb_to_cmyk_exact(double cmyk[4], double rgb[3])
 950 {
 951   cmyk[0] = 1 - rgb[0];
 952   cmyk[1] = 1 - rgb[1];
 953   cmyk[2] = 1 - rgb[2];
 954   cmyk[3] = MIN(cmyk[0], cmyk[1]);
 955   cmyk[3] = MIN(cmyk[3], cmyk[2]);
 956   if (cmyk[3] > 0.9999)
 957     {
 958       cmyk[3] = 1;
 959       cmyk[0] = cmyk[1] = cmyk[2] = 0;
 960     }
 961   else
 962     {
 963       double d = 1 / (1 - cmyk[3]);
 964       for (uns i = 0; i < 3; i++)
 965         cmyk[i] = d * (cmyk[i] - cmyk[3]);
 966     }
 967 }
 968
 969 /* CMYK to RGB */
 970 void
 971 cmyk_to_rgb_exact(double rgb[3], double cmyk[4])
 972 {
 973   double d = 1 - cmyk[1];
 974   for (uns i = 0; i < 3; i++)
 975     rgb[i] = d * (1 - cmyk[i]);
 976 }
 977
 978 /***************** OPTIMIZED SRGB -> LUV CONVERSION *********************/
 979
 980 u16 srgb_to_luv_tab1[256];
 981 u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE];
 982 u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE];
 983
 984 void
 985 srgb_to_luv_init(void)
 986 {
 987   DBG("Initializing sRGB -> Luv table");
 988   for (uns i = 0; i < 256; i++)
 989     {
 990       double t = i / 255.;
 991       if (t > 0.04045)
 992         t = pow((t + 0.055) * (1 / 1.055), 2.4);
 993       else
 994         t = t * (1 / 12.92);
 995       srgb_to_luv_tab1[i] = CLAMP(t * 0xfff + 0.5, 0, 0xfff);
 996     }
 997   for (uns i = 0; i < (9 << SRGB_TO_LUV_TAB2_SIZE); i++)
 998     {
 999       double t = i / (double)((9 << SRGB_TO_LUV_TAB2_SIZE) - 1);
1000       if (t > 0.008856)
1001         t = 1.16 * pow(t, 1 / 3.) - 0.16;
1002       else
1003         t = (1.16 * 7.787) * t;
1004       srgb_to_luv_tab2[i] =
1005         CLAMP(t * ((1 << SRGB_TO_LUV_TAB2_SCALE) - 1) + 0.5,
1006           0, (1 << SRGB_TO_LUV_TAB2_SCALE) - 1);
1007     }
1008   for (uns i = 0; i < (20 << SRGB_TO_LUV_TAB3_SIZE); i++)
1009     {
1010       srgb_to_luv_tab3[i] = i ? (13 << (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB3_SIZE)) / i : 0;
1011     }
1012 }
1013
1014 void
1015 srgb_to_luv_pixels(byte *dest, byte *src, uns count)
1016 {
1017   while (count--)
1018     {
1019       srgb_to_luv_pixel(dest, src);
1020       dest += 3;
1021       src += 3;
1022     }
1023 }
1024
1025
1026 /************************ GRID INTERPOLATION ALGORITHM ************************/
1027
1028 struct color_grid_node *srgb_to_luv_grid;
1029 struct color_interpolation_node *color_interpolation_table;
1030
1031 /* Returns volume of a given tetrahedron multiplied by 6 */
1032 static inline uns
1033 tetrahedron_volume(uns *v1, uns *v2, uns *v3, uns *v4)
1034 {
1035   int a[3], b[3], c[3];
1036   for (uns i = 0; i < 3; i++)
1037     {
1038       a[i] = v2[i] - v1[i];
1039       b[i] = v3[i] - v1[i];
1040       c[i] = v4[i] - v1[i];
1041     }
1042   int result =
1043     a[0] * (b[1] * c[2] - b[2] * c[1]) -
1044     a[1] * (b[0] * c[2] - b[2] * c[0]) +
1045     a[2] * (b[0] * c[1] - b[1] * c[0]);
1046   return (result > 0) ? result : -result;
1047 }
1048
1049 static void
1050 interpolate_tetrahedron(struct color_interpolation_node *n, uns *p, const uns *c)
1051 {
1052   uns v[4][3];
1053   for (uns i = 0; i < 4; i++)
1054     {
1055       v[i][0] = (c[i] & 0001) ? (1 << COLOR_CONV_OFS) : 0;
1056       v[i][1] = (c[i] & 0010) ? (1 << COLOR_CONV_OFS) : 0;
1057       v[i][2] = (c[i] & 0100) ? (1 << COLOR_CONV_OFS) : 0;
1058       n->ofs[i] =
1059         ((c[i] & 0001) ? 1 : 0) +
1060         ((c[i] & 0010) ? (1 << COLOR_CONV_SIZE) : 0) +
1061         ((c[i] & 0100) ? (1 << (COLOR_CONV_SIZE * 2)) : 0);
1062     }
1063   uns vol = tetrahedron_volume(v[0], v[1], v[2], v[3]);
1064   n->mul[0] = ((tetrahedron_volume(p, v[1], v[2], v[3]) << 8) + (vol >> 1)) / vol;
1065   n->mul[1] = ((tetrahedron_volume(v[0], p, v[2], v[3]) << 8) + (vol >> 1)) / vol;
1066   n->mul[2] = ((tetrahedron_volume(v[0], v[1], p, v[3]) << 8) + (vol >> 1)) / vol;
1067   n->mul[3] = ((tetrahedron_volume(v[0], v[1], v[2], p) << 8) + (vol >> 1)) / vol;
1068   uns j;
1069   for (j = 0; j < 4; j++)
1070     if (n->mul[j])
1071       break;
1072   for (uns i = 0; i < 4; i++)
1073     if (n->mul[i] == 0)
1074       n->ofs[i] = n->ofs[j];
1075 }
1076
1077 static void
1078 interpolation_table_init(void)
1079 {
1080   DBG("Initializing color interpolation table");
1081   struct color_interpolation_node *n = color_interpolation_table =
1082     xmalloc(sizeof(struct color_interpolation_node) << (COLOR_CONV_OFS * 3));
1083   uns p[3];
1084   for (p[2] = 0; p[2] < (1 << COLOR_CONV_OFS); p[2]++)
1085     for (p[1] = 0; p[1] < (1 << COLOR_CONV_OFS); p[1]++)
1086       for (p[0] = 0; p[0] < (1 << COLOR_CONV_OFS); p[0]++)
1087         {
1088           uns index;
1089           static const uns tetrahedra[5][4] = {
1090             {0000, 0001, 0010, 0100},
1091             {0110, 0111, 0100, 0010},
1092             {0101, 0100, 0111, 0001},
1093             {0011, 0010, 0001, 0111},
1094             {0111, 0001, 0010, 0100}};
1095           if (p[0] + p[1] + p[2] <= (1 << COLOR_CONV_OFS))
1096             index = 0;
1097           else if ((1 << COLOR_CONV_OFS) + p[0] <= p[1] + p[2])
1098             index = 1;
1099           else if ((1 << COLOR_CONV_OFS) + p[1] <= p[0] + p[2])
1100             index = 2;
1101           else if ((1 << COLOR_CONV_OFS) + p[2] <= p[0] + p[1])
1102             index = 3;
1103           else
1104             index = 4;
1105           interpolate_tetrahedron(n, p, tetrahedra[index]);
1106           n++;
1107         }
1108 }
1109
1110 typedef void color_conv_func(double dest[3], double src[3]);
1111
1112 static void
1113 conv_grid_init(struct color_grid_node **grid, color_conv_func func)
1114 {
1115   if (*grid)
1116     return;
1117   struct color_grid_node *g = *grid = xmalloc((sizeof(struct color_grid_node)) << (COLOR_CONV_SIZE * 3));
1118   double src[3], dest[3];
1119   for (uns k = 0; k < (1 << COLOR_CONV_SIZE); k++)
1120     {
1121       src[2] = k * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
1122       for (uns j = 0; j < (1 << COLOR_CONV_SIZE); j++)
1123         {
1124           src[1] = j * (255/ (double)((1 << COLOR_CONV_SIZE) - 1));
1125           for (uns i = 0; i < (1 << COLOR_CONV_SIZE); i++)
1126             {
1127               src[0] = i * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
1128               func(dest, src);
1129               g->val[0] = CLAMP(dest[0] + 0.5, 0, 255);
1130               g->val[1] = CLAMP(dest[1] + 0.5, 0, 255);
1131               g->val[2] = CLAMP(dest[2] + 0.5, 0, 255);
1132               g++;
1133             }
1134         }
1135     }
1136 }
1137
1138 static void
1139 srgb_to_luv_func(double dest[3], double src[3])
1140 {
1141   double srgb[3], xyz[3], luv[3];
1142   srgb[0] = src[0] / 255.;
1143   srgb[1] = src[1] / 255.;
1144   srgb[2] = src[2] / 255.;
1145   srgb_to_xyz_exact(xyz, srgb);
1146   xyz_to_luv_exact(luv, xyz);
1147   dest[0] = luv[0] * 2.55;
1148   dest[1] = luv[1] * (2.55 / 4) + 128;
1149   dest[2] = luv[2] * (2.55 / 4) + 128;
1150 }
1151
1152 void
1153 color_conv_init(void)
1154 {
1155   interpolation_table_init();
1156   conv_grid_init(&srgb_to_luv_grid, srgb_to_luv_func);
1157 }
1158
1159 void
1160 color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid)
1161 {
1162   while (count--)
1163     {
1164       color_conv_pixel(dest, src, grid);
1165       dest += 3;
1166       src += 3;
1167     }
1168 }
1169
1170
1171 /**************************** TESTS *******************************/
1172
1173 #ifdef TEST
1174 #include <string.h>
1175
1176 static double
1177 conv_error(u32 color, struct color_grid_node *grid, color_conv_func func)
1178 {
1179   byte src[3], dest[3];
1180   src[0] = color & 255;
1181   src[1] = (color >> 8) & 255;
1182   src[2] = (color >> 16) & 255;
1183   color_conv_pixel(dest, src, grid);
1184   double src2[3], dest2[3];
1185   for (uns i = 0; i < 3; i++)
1186     src2[i] = src[i];
1187   func(dest2, src2);
1188   double err = 0;
1189   for (uns i = 0; i < 3; i++)
1190     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
1191   return err;
1192 }
1193
1194 typedef void test_fn(byte *dest, byte *src);
1195
1196 static double
1197 func_error(u32 color, test_fn test, color_conv_func func)
1198 {
1199   byte src[3], dest[3];
1200   src[0] = color & 255;
1201   src[1] = (color >> 8) & 255;
1202   src[2] = (color >> 16) & 255;
1203   test(dest, src);
1204   double src2[3], dest2[3];
1205   for (uns i = 0; i < 3; i++)
1206     src2[i] = src[i];
1207   func(dest2, src2);
1208   double err = 0;
1209   for (uns i = 0; i < 3; i++)
1210     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
1211   return err;
1212 }
1213
1214 static void
1215 test_grid(byte *name, struct color_grid_node *grid, color_conv_func func)
1216 {
1217   double max_err = 0, sum_err = 0;
1218   uns count = 100000;
1219   for (uns i = 0; i < count; i++)
1220     {
1221       double err = conv_error(random_max(0x1000000), grid, func);
1222       max_err = MAX(err, max_err);
1223       sum_err += err;
1224     }
1225   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
1226   if (max_err > 12)
1227     die("Too large error in %s conversion", name);
1228 }
1229
1230 static void
1231 test_func(byte *name, test_fn test, color_conv_func func)
1232 {
1233   double max_err = 0, sum_err = 0;
1234   uns count = 100000;
1235   for (uns i = 0; i < count; i++)
1236     {
1237       double err = func_error(random_max(0x1000000), test, func);
1238       max_err = MAX(err, max_err);
1239       sum_err += err;
1240     }
1241   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
1242   if (max_err > 12)
1243     die("Too large error in %s conversion", name);
1244 }
1245
1246 int
1247 main(void)
1248 {
1249   srgb_to_luv_init();
1250   test_func("func sRGB -> Luv", srgb_to_luv_pixel, srgb_to_luv_func);
1251   color_conv_init();
1252   test_grid("grid sRGB -> Luv", srgb_to_luv_grid, srgb_to_luv_func);
1253 #ifdef LOCAL_DEBUG
1254 #define CNT 1000000
1255 #define TESTS 10
1256   byte *a = xmalloc(3 * CNT), *b = xmalloc(3 * CNT);
1257   for (uns i = 0; i < 3 * CNT; i++)
1258     a[i] = random_max(256);
1259   timestamp_t timer;
1260   init_timer(&timer);
1261   for (uns i = 0; i < TESTS; i++)
1262     memcpy(b, a, CNT * 3);
1263   DBG("memcpy time=%d", get_timer(&timer));
1264   init_timer(&timer);
1265   for (uns i = 0; i < TESTS; i++)
1266     srgb_to_luv_pixels(b, a, CNT);
1267   DBG("direct time=%d", get_timer(&timer));
1268   init_timer(&timer);
1269   for (uns i = 0; i < TESTS; i++)
1270     color_conv_pixels(b, a, CNT, srgb_to_luv_grid);
1271   DBG("grid time=%d", get_timer(&timer));
1272 #endif
1273   return 0;
1274 }
1275 #endif
1276