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 "sherlock/sherlock.h"
  13 #include "lib/math.h"
  14 #include "images/images.h"
  15 #include "images/color.h"
  16 #include "images/error.h"
  17 #include "images/math.h"
  18
  19 #include <string.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 /* Main functions */
 252
 253 static int
 254 image_conv_color_space(struct image_context *ctx UNUSED, struct image *dest, struct image *src, struct image_conv_options *opt UNUSED)
 255 {
 256   switch (dest->flags & IMAGE_COLOR_SPACE)
 257     {
 258       case COLOR_SPACE_GRAYSCALE:
 259         switch (src->flags & IMAGE_COLOR_SPACE)
 260           {
 261             case COLOR_SPACE_RGB:
 262               if (dest->pixel_size == 1)
 263                 {
 264                   image_conv_rgb_n_to_gray_1(dest, src);
 265                   return 1;
 266                 }
 267               break;
 268           }
 269         break;
 270       case COLOR_SPACE_RGB:
 271         switch (src->flags & IMAGE_CHANNELS_FORMAT)
 272           {
 273             case COLOR_SPACE_GRAYSCALE:
 274               if (src->pixel_size == 1)
 275                 {
 276                   image_conv_gray_1_to_rgb_n(dest, src);
 277                   return 1;
 278                 }
 279               break;
 280             case COLOR_SPACE_YCBCR:
 281               image_conv_ycbcr_n_to_rgb_n(dest, src);
 282               return 1;
 283             case COLOR_SPACE_CMYK:
 284               if (src->pixel_size == 4)
 285                 {
 286                   image_conv_cmyk_4_to_rgb_n(dest, src);
 287                   return 1;
 288                 }
 289               break;
 290         }
 291         break;
 292       case COLOR_SPACE_YCBCR:
 293         switch (src->flags & IMAGE_CHANNELS_FORMAT)
 294           {
 295             case COLOR_SPACE_RGB:
 296               image_conv_rgb_n_to_ycbcr_n(dest, src);
 297               return 1;
 298           }
 299         break;
 300       case COLOR_SPACE_CMYK:
 301         switch (src->flags & IMAGE_CHANNELS_FORMAT)
 302           {
 303             case COLOR_SPACE_RGB:
 304               if (dest->pixel_size == 4)
 305                 {
 306                   image_conv_rgb_n_to_cmyk_4(dest, src);
 307                   return 1;
 308                 }
 309               break;
 310           }
 311         break;
 312     }
 313   return 0;
 314 }
 315
 316 static void
 317 image_conv_copy(struct image *dest, struct image *src)
 318 {
 319   if (dest->pixels == src->pixels)
 320     return;
 321   else if (dest->pixel_size != src->pixel_size)
 322     {
 323       uns channels = MIN(dest->channels, src->channels);
 324       switch (channels)
 325         {
 326           case 1:
 327             {
 328 #             define IMAGE_WALK_PREFIX(x) walk_##x
 329 #             define IMAGE_WALK_INLINE
 330 #             define IMAGE_WALK_DOUBLE
 331 #             define IMAGE_WALK_IMAGE dest
 332 #             define IMAGE_WALK_SEC_IMAGE src
 333 #             define IMAGE_WALK_UNROLL 4
 334 #             define IMAGE_WALK_DO_STEP do{ walk_pos[0] = walk_sec_pos[0]; }while(0)
 335 #             include "images/image-walk.h"
 336             }
 337             return;
 338           case 2:
 339 #             define IMAGE_WALK_PREFIX(x) walk_##x
 340 #             define IMAGE_WALK_INLINE
 341 #             define IMAGE_WALK_DOUBLE
 342 #             define IMAGE_WALK_IMAGE dest
 343 #             define IMAGE_WALK_SEC_IMAGE src
 344 #             define IMAGE_WALK_UNROLL 4
 345 #             define IMAGE_WALK_DO_STEP do{ walk_pos[0] = walk_sec_pos[0]; walk_pos[1] = walk_sec_pos[1]; }while(0)
 346 #             include "images/image-walk.h"
 347             return;
 348           case 3:
 349 #             define IMAGE_WALK_PREFIX(x) walk_##x
 350 #             define IMAGE_WALK_INLINE
 351 #             define IMAGE_WALK_DOUBLE
 352 #             define IMAGE_WALK_IMAGE dest
 353 #             define IMAGE_WALK_SEC_IMAGE src
 354 #             define IMAGE_WALK_UNROLL 2
 355 #             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)
 356 #             include "images/image-walk.h"
 357             return;
 358           case 4:
 359 #             define IMAGE_WALK_PREFIX(x) walk_##x
 360 #             define IMAGE_WALK_INLINE
 361 #             define IMAGE_WALK_DOUBLE
 362 #             define IMAGE_WALK_IMAGE dest
 363 #             define IMAGE_WALK_SEC_IMAGE src
 364 #             define IMAGE_WALK_UNROLL 2
 365 #             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)
 366 #             include "images/image-walk.h"
 367             return;
 368           default:
 369 #             define IMAGE_WALK_PREFIX(x) walk_##x
 370 #             define IMAGE_WALK_INLINE
 371 #             define IMAGE_WALK_DOUBLE
 372 #             define IMAGE_WALK_IMAGE dest
 373 #             define IMAGE_WALK_SEC_IMAGE src
 374 #             define IMAGE_WALK_DO_STEP do{ for (uns i = 0; i < channels; i++) walk_pos[i] = walk_sec_pos[i]; }while(0)
 375 #             include "images/image-walk.h"
 376             return;
 377         }
 378     }
 379   else if (dest->row_size != src->row_size || ((dest->flags | src->flags) & IMAGE_GAPS_PROTECTED))
 380     {
 381       byte *s = src->pixels;
 382       byte *d = dest->pixels;
 383       for (uns row = src->rows; row--; )
 384         {
 385           memcpy(d, s, src->row_pixels_size);
 386           d += dest->row_size;
 387           s += src->row_size;
 388         }
 389     }
 390   else if (dest->pixels != src->pixels)
 391     memcpy(dest->pixels, src->pixels, src->image_size);
 392 }
 393
 394 static void
 395 image_conv_fill_alpha(struct image *dest)
 396 {
 397   switch (dest->channels)
 398     {
 399       case 2:
 400         if (dest->pixel_size == 2)
 401           {
 402 #           define IMAGE_WALK_PREFIX(x) walk_##x
 403 #           define IMAGE_WALK_INLINE
 404 #           define IMAGE_WALK_IMAGE dest
 405 #           define IMAGE_WALK_COL_STEP 2
 406 #           define IMAGE_WALK_UNROLL 4
 407 #           define IMAGE_WALK_DO_STEP do{ walk_pos[1] = 255; }while(0)
 408 #           include "images/image-walk.h"
 409             return;
 410           }
 411         break;
 412       case 4:
 413         if (dest->pixel_size == 4)
 414           {
 415 #           define IMAGE_WALK_PREFIX(x) walk_##x
 416 #           define IMAGE_WALK_INLINE
 417 #           define IMAGE_WALK_IMAGE dest
 418 #           define IMAGE_WALK_COL_STEP 4
 419 #           define IMAGE_WALK_UNROLL 4
 420 #           define IMAGE_WALK_DO_STEP do{ walk_pos[3] = 255; }while(0)
 421 #           include "images/image-walk.h"
 422             return;
 423           }
 424         break;
 425     }
 426   {
 427 #   define IMAGE_WALK_PREFIX(x) walk_##x
 428 #   define IMAGE_WALK_INLINE
 429 #   define IMAGE_WALK_IMAGE dest
 430 #   define IMAGE_WALK_UNROLL 4
 431 #   define IMAGE_WALK_DO_STEP do{ walk_pos[dest->channels - 1] = 255; }while(0)
 432 #   include "images/image-walk.h"
 433   }
 434 }
 435
 436 static void
 437 image_conv_copy_alpha(struct image *dest, struct image *src)
 438 {
 439   if (dest->pixels != src->pixels || dest->channels != src->channels)
 440     {
 441 #     define IMAGE_WALK_PREFIX(x) walk_##x
 442 #     define IMAGE_WALK_INLINE
 443 #     define IMAGE_WALK_DOUBLE
 444 #     define IMAGE_WALK_IMAGE dest
 445 #     define IMAGE_WALK_SEC_IMAGE src
 446 #     define IMAGE_WALK_UNROLL 4
 447 #     define IMAGE_WALK_DO_STEP do{ walk_pos[dest->channels - 1] = walk_sec_pos[src->channels - 1]; }while(0)
 448 #     include "images/image-walk.h"
 449     }
 450 }
 451
 452 static inline uns
 453 image_conv_alpha_func(uns value, uns alpha, uns acoef, uns bcoef)
 454 {
 455   return ((uns)(acoef + (int)alpha * (int)(value - bcoef)) * (0xffffffffU / 255 / 255)) >> 24;
 456 }
 457
 458 static int
 459 image_conv_apply_alpha_from(struct image_context *ctx, struct image *dest, struct image *src, struct image_conv_options *opt)
 460 {
 461   if (!opt->background.color_space)
 462     return 1;
 463   byte background[IMAGE_MAX_CHANNELS];
 464   if (unlikely(!color_put(ctx, &opt->background, background, dest->flags & IMAGE_COLOR_SPACE)))
 465     return 0;
 466   uns a[IMAGE_MAX_CHANNELS], b[IMAGE_MAX_CHANNELS];
 467   for (uns i = 0; i < dest->channels; i++)
 468     a[i] = 255 * (b[i] = background[i]);
 469   switch (dest->channels)
 470     {
 471       case 1:
 472         {
 473 #         define IMAGE_WALK_PREFIX(x) walk_##x
 474 #         define IMAGE_WALK_INLINE
 475 #         define IMAGE_WALK_IMAGE dest
 476 #         define IMAGE_WALK_SEC_IMAGE src
 477 #         define IMAGE_WALK_DOUBLE
 478 #         define IMAGE_WALK_UNROLL 2
 479 #         define IMAGE_WALK_DO_STEP do{ \
 480               walk_pos[0] = image_conv_alpha_func(walk_pos[0], walk_sec_pos[src->channels - 1], a[0], b[0]); }while(0)
 481 #         include "images/image-walk.h"
 482         }
 483         return 1;
 484       case 3:
 485         {
 486 #         define IMAGE_WALK_PREFIX(x) walk_##x
 487 #         define IMAGE_WALK_INLINE
 488 #         define IMAGE_WALK_IMAGE dest
 489 #         define IMAGE_WALK_SEC_IMAGE src
 490 #         define IMAGE_WALK_DOUBLE
 491 #         define IMAGE_WALK_DO_STEP do{ \
 492               walk_pos[0] = image_conv_alpha_func(walk_pos[0], walk_sec_pos[src->channels - 1], a[0], b[0]); \
 493               walk_pos[1] = image_conv_alpha_func(walk_pos[1], walk_sec_pos[src->channels - 1], a[1], b[1]); \
 494               walk_pos[2] = image_conv_alpha_func(walk_pos[2], walk_sec_pos[src->channels - 1], a[2], b[2]); }while(0)
 495 #         include "images/image-walk.h"
 496         }
 497         return 1;
 498     }
 499   {
 500 #   define IMAGE_WALK_PREFIX(x) walk_##x
 501 #   define IMAGE_WALK_INLINE
 502 #   define IMAGE_WALK_IMAGE dest
 503 #   define IMAGE_WALK_SEC_IMAGE src
 504 #   define IMAGE_WALK_DOUBLE
 505 #   define IMAGE_WALK_DO_STEP do{ for (uns i = 0; i < dest->channels; i++) \
 506         walk_pos[i] = image_conv_alpha_func(walk_pos[i], walk_sec_pos[src->channels - 1], a[i], b[i]); }while(0)
 507 #   include "images/image-walk.h"
 508   }
 509   return 1;
 510 }
 511
 512 static int
 513 image_conv_apply_alpha_to(struct image_context *ctx, struct image *dest, struct image *src, struct image_conv_options *opt)
 514 {
 515   if (!opt->background.color_space)
 516     {
 517       image_conv_copy(dest, src);
 518       return 1;
 519     }
 520   byte background[IMAGE_MAX_CHANNELS];
 521   if (unlikely(!color_put(ctx, &opt->background, background, dest->flags & IMAGE_COLOR_SPACE)))
 522     return 0;
 523   uns a[IMAGE_MAX_CHANNELS], b[IMAGE_MAX_CHANNELS];
 524   for (uns i = 0; i < dest->channels; i++)
 525     a[i] = 255 * (b[i] = background[i]);
 526   switch (dest->channels)
 527     {
 528       case 1:
 529         {
 530 #         define IMAGE_WALK_PREFIX(x) walk_##x
 531 #         define IMAGE_WALK_INLINE
 532 #         define IMAGE_WALK_IMAGE dest
 533 #         define IMAGE_WALK_SEC_IMAGE src
 534 #         define IMAGE_WALK_DOUBLE
 535 #         define IMAGE_WALK_UNROLL 2
 536 #         define IMAGE_WALK_DO_STEP do{ \
 537               walk_pos[0] = image_conv_alpha_func(walk_sec_pos[0], walk_sec_pos[src->channels - 1], a[0], b[0]); }while(0)
 538 #         include "images/image-walk.h"
 539         }
 540         return 1;
 541       case 3:
 542         {
 543 #         define IMAGE_WALK_PREFIX(x) walk_##x
 544 #         define IMAGE_WALK_INLINE
 545 #         define IMAGE_WALK_IMAGE dest
 546 #         define IMAGE_WALK_SEC_IMAGE src
 547 #         define IMAGE_WALK_DOUBLE
 548 #         define IMAGE_WALK_DO_STEP do{ \
 549               walk_pos[0] = image_conv_alpha_func(walk_sec_pos[0], walk_sec_pos[src->channels - 1], a[0], b[0]); \
 550               walk_pos[1] = image_conv_alpha_func(walk_sec_pos[1], walk_sec_pos[src->channels - 1], a[1], b[1]); \
 551               walk_pos[2] = image_conv_alpha_func(walk_sec_pos[2], walk_sec_pos[src->channels - 1], a[2], b[2]); }while(0)
 552 #         include "images/image-walk.h"
 553         }
 554         return 1;
 555     }
 556   {
 557 #   define IMAGE_WALK_PREFIX(x) walk_##x
 558 #   define IMAGE_WALK_INLINE
 559 #   define IMAGE_WALK_IMAGE dest
 560 #   define IMAGE_WALK_SEC_IMAGE src
 561 #   define IMAGE_WALK_DOUBLE
 562 #   define IMAGE_WALK_DO_STEP do{ for (uns i = 0; i < dest->channels; i++) \
 563         walk_pos[i] = image_conv_alpha_func(walk_sec_pos[i], walk_sec_pos[src->channels - 1], a[i], b[i]); }while(0)
 564 #   include "images/image-walk.h"
 565   }
 566   return 1;
 567 }
 568
 569 int
 570 image_conv(struct image_context *ctx, struct image *dest, struct image *src, struct image_conv_options *opt)
 571 {
 572   ASSERT(dest->cols == src->cols && dest->rows == src->rows);
 573   if (!((dest->flags ^ src->flags) & IMAGE_COLOR_SPACE))
 574     {
 575       if (!(src->flags & IMAGE_ALPHA) || (dest->flags & IMAGE_ALPHA))
 576         image_conv_copy(dest, src);
 577       else if (unlikely(!image_conv_apply_alpha_to(ctx, dest, src, opt)))
 578         return 0;
 579     }
 580   else
 581     {
 582       if (!(src->flags & IMAGE_ALPHA))
 583         {
 584           if (unlikely(!image_conv_color_space(ctx, dest, src, opt)))
 585             goto error;
 586           if ((dest->flags & IMAGE_ALPHA) && (opt->flags & IMAGE_CONV_FILL_ALPHA))
 587             image_conv_fill_alpha(dest);
 588         }
 589       else
 590         {
 591           if (dest->flags & IMAGE_ALPHA)
 592             {
 593               if (dest->channels <= src->channels)
 594                 {
 595                   if (unlikely(!image_conv_color_space(ctx, dest, src, opt)))
 596                     goto error;
 597                   if (opt->flags & IMAGE_CONV_COPY_ALPHA)
 598                     image_conv_copy_alpha(dest, src);
 599                   else if (opt->flags & IMAGE_CONV_FILL_ALPHA)
 600                     image_conv_fill_alpha(dest);
 601                 }
 602               else
 603                 {
 604                   if (opt->flags & IMAGE_CONV_COPY_ALPHA)
 605                     image_conv_copy_alpha(dest, src);
 606                   else
 607                     image_conv_fill_alpha(dest);
 608                   if (unlikely(!image_conv_color_space(ctx, dest, src, opt)))
 609                     goto error;
 610                 }
 611             }
 612           else
 613             {
 614               if (dest->channels <= src->channels)
 615                 {
 616                   if (unlikely(!image_conv_color_space(ctx, dest, src, opt)))
 617                     goto error;
 618                   if (unlikely(!image_conv_apply_alpha_from(ctx, dest, src, opt)))
 619                     return 0;
 620                 }
 621               else
 622                 {
 623                   if (unlikely(!image_conv_apply_alpha_to(ctx, dest, src, opt)))
 624                     return 0;
 625                   if (unlikely(!image_conv_color_space(ctx, dest, dest, opt)))
 626                     goto error;
 627                 }
 628             }
 629         }
 630     }
 631   return 1;
 632 error:
 633   IMAGE_ERROR(ctx, IMAGE_ERROR_INVALID_PIXEL_FORMAT, "Image conversion not supported for such pixel formats");
 634   return 0;
 635 }
 636
 637 /********************* EXACT CONVERSION ROUTINES **********************/
 638
 639 /* Reference whites */
 640 #define COLOR_ILLUMINANT_A      0.44757, 0.40744
 641 #define COLOR_ILLUMINANT_B      0.34840, 0.35160
 642 #define COLOR_ILLUMINANT_C      0.31006, 0.31615
 643 #define COLOR_ILLUMINANT_D50    0.34567, 0.35850
 644 #define COLOR_ILLUMINANT_D55    0.33242, 0.34743
 645 #define COLOR_ILLUMINANT_D65    0.31273, 0.32902
 646 #define COLOR_ILLUMINANT_D75    0.29902, 0.31485
 647 #define COLOR_ILLUMINANT_9300K  0.28480, 0.29320
 648 #define COLOR_ILLUMINANT_E      (1./3.), (1./3.)
 649 #define COLOR_ILLUMINANT_F2     0.37207, 0.37512
 650 #define COLOR_ILLUMINANT_F7     0.31285, 0.32918
 651 #define COLOR_ILLUMINANT_F11    0.38054, 0.37691
 652
 653 const double
 654   color_illuminant_d50[2] = {COLOR_ILLUMINANT_D50},
 655   color_illuminant_d65[2] = {COLOR_ILLUMINANT_D65},
 656   color_illuminant_e[2] = {COLOR_ILLUMINANT_E};
 657
 658 /* RGB profiles (many missing) */
 659 const struct color_space_info
 660   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}},
 661   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}},
 662   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}},
 663   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}},
 664   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}};
 665
 666 #define CLIP(x, min, max) (((x) < (min)) ? (min) : ((x) > (max)) ? (max) : (x))
 667
 668 static inline void
 669 clip(double a[3])
 670 {
 671   a[0] = CLIP(a[0], 0, 1);
 672   a[1] = CLIP(a[1], 0, 1);
 673   a[2] = CLIP(a[2], 0, 1);
 674 }
 675
 676 static inline void
 677 correct_gamma_simple(double dest[3], double src[3], const struct color_space_gamma_info *info)
 678 {
 679   dest[0] = pow(src[0], info->simple_gamma);
 680   dest[1] = pow(src[1], info->simple_gamma);
 681   dest[2] = pow(src[2], info->simple_gamma);
 682 }
 683
 684 static inline void
 685 invert_gamma_simple(double dest[3], double src[3], const struct color_space_gamma_info *info)
 686 {
 687   dest[0] = pow(src[0], 1 / info->simple_gamma);
 688   dest[1] = pow(src[1], 1 / info->simple_gamma);
 689   dest[2] = pow(src[2], 1 / info->simple_gamma);
 690 }
 691
 692 static inline void
 693 correct_gamma_detailed(double dest[3], double src[3], const struct color_space_gamma_info *info)
 694 {
 695   for (uns i = 0; i < 3; i++)
 696     if (src[i] > info->transition)
 697       dest[i] = (1 + info->offset) * pow(src[i], info->detailed_gamma) - info->offset;
 698     else
 699       dest[i] = info->slope * src[i];
 700 }
 701
 702 static inline void
 703 invert_gamma_detailed(double dest[3], double src[3], const struct color_space_gamma_info *info)
 704 {
 705   for (uns i = 0; i < 3; i++)
 706     if (src[i] > info->transition * info->slope)
 707       dest[i] = pow((src[i] + info->offset) / (1 + info->offset), 1 / info->detailed_gamma);
 708     else
 709       dest[i] = src[i] / info->slope;
 710 }
 711
 712 static inline void
 713 apply_matrix(double dest[3], double src[3], double matrix[9])
 714 {
 715   dest[0] = src[0] * matrix[0] + src[1] * matrix[1] + src[2] * matrix[2];
 716   dest[1] = src[0] * matrix[3] + src[1] * matrix[4] + src[2] * matrix[5];
 717   dest[2] = src[0] * matrix[6] + src[1] * matrix[7] + src[2] * matrix[8];
 718 }
 719
 720 void
 721 color_invert_matrix(double dest[9], double matrix[9])
 722 {
 723   double *i = dest, *m = matrix;
 724   double a0 = m[4] * m[8] - m[5] * m[7];
 725   double a1 = m[3] * m[8] - m[5] * m[6];
 726   double a2 = m[3] * m[7] - m[4] * m[6];
 727   double d = 1 / (m[0] * a0 - m[1] * a1 + m[2] * a2);
 728   i[0] = d * a0;
 729   i[3] = -d * a1;
 730   i[6] = d * a2;
 731   i[1] = -d * (m[1] * m[8] - m[2] * m[7]);
 732   i[4] = d * (m[0] * m[8] - m[2] * m[6]);
 733   i[7] = -d * (m[0] * m[7] - m[1] * m[6]);
 734   i[2] = d * (m[1] * m[5] - m[2] * m[4]);
 735   i[5] = -d * (m[0] * m[5] - m[2] * m[3]);
 736   i[8] = d * (m[0] * m[4] - m[1] * m[3]);
 737 }
 738
 739 static void
 740 mul_matrices(double r[9], double a[9], double b[9])
 741 {
 742   r[0] = a[0] * b[0] + a[1] * b[3] + a[2] * b[6];
 743   r[1] = a[0] * b[1] + a[1] * b[4] + a[2] * b[7];
 744   r[2] = a[0] * b[2] + a[1] * b[5] + a[2] * b[8];
 745   r[3] = a[3] * b[0] + a[4] * b[3] + a[5] * b[6];
 746   r[4] = a[3] * b[1] + a[4] * b[4] + a[5] * b[7];
 747   r[5] = a[3] * b[2] + a[4] * b[5] + a[5] * b[8];
 748   r[6] = a[6] * b[0] + a[7] * b[3] + a[8] * b[6];
 749   r[7] = a[6] * b[1] + a[7] * b[4] + a[8] * b[7];
 750   r[8] = a[6] * b[2] + a[7] * b[5] + a[8] * b[8];
 751 }
 752
 753 /* computes conversion matrix from a given color space to CIE XYZ */
 754 void
 755 color_compute_color_space_to_xyz_matrix(double matrix[9], const struct color_space_chromacity_info *space)
 756 {
 757   double wX = space->white[0] / space->white[1];
 758   double wZ = (1 - space->white[0] - space->white[1]) / space->white[1];
 759   double a[9], b[9];
 760   a[0] = space->prim1[0]; a[3] = space->prim1[1]; a[6] = 1 - a[0] - a[3];
 761   a[1] = space->prim2[0]; a[4] = space->prim2[1]; a[7] = 1 - a[1] - a[4];
 762   a[2] = space->prim3[0]; a[5] = space->prim3[1]; a[8] = 1 - a[2] - a[5];
 763   color_invert_matrix(b, a);
 764   double ra = wX * b[0] + b[1] + wZ * b[2];
 765   double rb = wX * b[3] + b[4] + wZ * b[5];
 766   double rc = wX * b[6] + b[7] + wZ * b[8];
 767   matrix[0] = a[0] * ra;
 768   matrix[1] = a[1] * rb;
 769   matrix[2] = a[2] * rc;
 770   matrix[3] = a[3] * ra;
 771   matrix[4] = a[4] * rb;
 772   matrix[5] = a[5] * rc;
 773   matrix[6] = a[6] * ra;
 774   matrix[7] = a[7] * rb;
 775   matrix[8] = a[8] * rc;
 776 }
 777
 778 /* computes matrix to join transformations with different reference whites */
 779 void
 780 color_compute_bradford_matrix(double matrix[9], const double source[2], const double dest[2])
 781 {
 782   /* cone response matrix and its inversion */
 783   static double r[9] = {
 784     0.8951, 0.2664, -0.1614,
 785     -0.7502, 1.7135, 0.0367,
 786     0.0389, -0.0685, 1.0296};
 787   //static double i[9] = {0.9870, -0.1471, 0.1600, 0.4323, 0.5184, 0.0493, -0.0085, 0.0400, 0.9685};
 788   double i[9];
 789   color_invert_matrix(i, r);
 790   double aX = source[0] / source[1];
 791   double aZ = (1 - source[0] - source[1]) / source[1];
 792   double bX = dest[0] / dest[1];
 793   double bZ = (1 - dest[0] - dest[1]) / dest[1];
 794   double x = (r[0] * bX + r[1] + r[2] * bZ) / (r[0] * aX + r[1] + r[2] * aZ);
 795   double y = (r[3] * bX + r[4] + r[5] * bZ) / (r[3] * aX + r[4] + r[5] * aZ);
 796   double z = (r[6] * bX + r[7] + r[8] * bZ) / (r[6] * aX + r[7] + r[8] * aZ);
 797   double m[9];
 798   m[0] = i[0] * x; m[1] = i[1] * y; m[2] = i[2] * z;
 799   m[3] = i[3] * x; m[4] = i[4] * y; m[5] = i[5] * z;
 800   m[6] = i[6] * x; m[7] = i[7] * y; m[8] = i[8] * z;
 801   mul_matrices(matrix, m, r);
 802 }
 803
 804 void
 805 color_compute_color_spaces_conversion_matrix(double matrix[9], const struct color_space_chromacity_info *src, const struct color_space_chromacity_info *dest)
 806 {
 807   double a_to_xyz[9], b_to_xyz[9], xyz_to_b[9], bradford[9], m[9];
 808   color_compute_color_space_to_xyz_matrix(a_to_xyz, src);
 809   color_compute_color_space_to_xyz_matrix(b_to_xyz, dest);
 810   color_invert_matrix(xyz_to_b, b_to_xyz);
 811   if (src->white[0] == dest->white[0] && src->white[1] == dest->white[1])
 812     mul_matrices(matrix, a_to_xyz, xyz_to_b);
 813   else
 814     {
 815       color_compute_bradford_matrix(bradford, src->white, dest->white);
 816       mul_matrices(m, a_to_xyz, bradford);
 817       mul_matrices(matrix, m, xyz_to_b);
 818     }
 819 }
 820
 821 /* sRGB to XYZ */
 822 void
 823 srgb_to_xyz_exact(double xyz[3], double srgb[3])
 824 {
 825   static double matrix[9] = {
 826     0.41248031, 0.35756952, 0.18043951,
 827     0.21268516, 0.71513904, 0.07217580,
 828     0.01933501, 0.11918984, 0.95031473};
 829   double srgb_lin[3];
 830   invert_gamma_detailed(srgb_lin, srgb, &color_srgb_info.gamma);
 831   apply_matrix(xyz, srgb_lin, matrix);
 832   xyz_to_srgb_exact(srgb_lin, xyz);
 833 }
 834
 835 /* XYZ to sRGB */
 836 void
 837 xyz_to_srgb_exact(double srgb[3], double xyz[3])
 838 {
 839   static double matrix[9] = {
 840      3.24026666, -1.53704957, -0.49850256,
 841     -0.96928381,  1.87604525,  0.04155678,
 842      0.05564281, -0.20402363,  1.05721334};
 843   double srgb_lin[3];
 844   apply_matrix(srgb_lin, xyz, matrix);
 845   clip(srgb_lin);
 846   correct_gamma_detailed(srgb, srgb_lin, &color_srgb_info.gamma);
 847 }
 848
 849 /* XYZ to CIE-Luv */
 850 void
 851 xyz_to_luv_exact(double luv[3], double xyz[3])
 852 {
 853   double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
 854   if (sum < 0.000001)
 855     luv[0] = luv[1] = luv[2] = 0;
 856   else
 857     {
 858       double var_u = 4 * xyz[0] / sum;
 859       double var_v = 9 * xyz[1] / sum;
 860       if (xyz[1] > 0.008856)
 861         luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
 862       else
 863         luv[0] = (116 * 7.787) * xyz[1];
 864      luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
 865      luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
 866      /* intervals [0..100], [-134..220], [-140..122] */
 867    }
 868 }
 869
 870 /* CIE-Luv to XYZ */
 871 void
 872 luv_to_xyz_exact(double xyz[3], double luv[3])
 873 {
 874   double var_u = luv[1] / (13 * luv[0]) + (4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z));
 875   double var_v = luv[2] / (13 * luv[0]) + (9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z));
 876   double var_y = (luv[0] + 16) / 116;
 877   double pow_y = var_y * var_y * var_y;
 878   if (pow_y > 0.008856)
 879     var_y = pow_y;
 880   else
 881     var_y = (var_y - 16 / 116) / 7.787;
 882   xyz[1] = var_y;
 883   xyz[0] = -(9 * xyz[1] * var_u) / ((var_u - 4) * var_v - var_u * var_v);
 884   xyz[2] = (9 * xyz[1] - 15 * var_v * xyz[1] - var_v * xyz[0]) / (3 * var_v);
 885 }
 886
 887 /* RGB to CMYK - a very simple version, not too accureate  */
 888 void
 889 rgb_to_cmyk_exact(double cmyk[4], double rgb[3])
 890 {
 891   cmyk[0] = 1 - rgb[0];
 892   cmyk[1] = 1 - rgb[1];
 893   cmyk[2] = 1 - rgb[2];
 894   cmyk[3] = MIN(cmyk[0], cmyk[1]);
 895   cmyk[3] = MIN(cmyk[3], cmyk[2]);
 896   if (cmyk[3] > 0.9999)
 897     {
 898       cmyk[3] = 1;
 899       cmyk[0] = cmyk[1] = cmyk[2] = 0;
 900     }
 901   else
 902     {
 903       double d = 1 / (1 - cmyk[3]);
 904       for (uns i = 0; i < 3; i++)
 905         cmyk[i] = d * (cmyk[i] - cmyk[3]);
 906     }
 907 }
 908
 909 /* CMYK to RGB */
 910 void
 911 cmyk_to_rgb_exact(double rgb[3], double cmyk[4])
 912 {
 913   double d = 1 - cmyk[1];
 914   for (uns i = 0; i < 3; i++)
 915     rgb[i] = d * (1 - cmyk[i]);
 916 }
 917
 918 /***************** OPTIMIZED SRGB -> LUV CONVERSION *********************/
 919
 920 u16 srgb_to_luv_tab1[256];
 921 u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE];
 922 u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE];
 923
 924 void
 925 srgb_to_luv_init(void)
 926 {
 927   DBG("Initializing sRGB -> Luv table");
 928   for (uns i = 0; i < 256; i++)
 929     {
 930       double t = i / 255.;
 931       if (t > 0.04045)
 932         t = pow((t + 0.055) * (1 / 1.055), 2.4);
 933       else
 934         t = t * (1 / 12.92);
 935       srgb_to_luv_tab1[i] = CLAMP(t * 0xfff + 0.5, 0, 0xfff);
 936     }
 937   for (uns i = 0; i < (9 << SRGB_TO_LUV_TAB2_SIZE); i++)
 938     {
 939       double t = i / (double)((9 << SRGB_TO_LUV_TAB2_SIZE) - 1);
 940       if (t > 0.008856)
 941         t = 1.16 * pow(t, 1 / 3.) - 0.16;
 942       else
 943         t = (1.16 * 7.787) * t;
 944       srgb_to_luv_tab2[i] =
 945         CLAMP(t * ((1 << SRGB_TO_LUV_TAB2_SCALE) - 1) + 0.5,
 946           0, (1 << SRGB_TO_LUV_TAB2_SCALE) - 1);
 947     }
 948   for (uns i = 0; i < (20 << SRGB_TO_LUV_TAB3_SIZE); i++)
 949     {
 950       srgb_to_luv_tab3[i] = i ? (13 << (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB3_SIZE)) / i : 0;
 951     }
 952 }
 953
 954 void
 955 srgb_to_luv_pixels(byte *dest, byte *src, uns count)
 956 {
 957   while (count--)
 958     {
 959       srgb_to_luv_pixel(dest, src);
 960       dest += 3;
 961       src += 3;
 962     }
 963 }
 964
 965
 966 /************************ GRID INTERPOLATION ALGORITHM ************************/
 967
 968 struct color_grid_node *srgb_to_luv_grid;
 969 struct color_interpolation_node *color_interpolation_table;
 970
 971 /* Returns volume of a given tetrahedron multiplied by 6 */
 972 static inline uns
 973 tetrahedron_volume(uns *v1, uns *v2, uns *v3, uns *v4)
 974 {
 975   int a[3], b[3], c[3];
 976   for (uns i = 0; i < 3; i++)
 977     {
 978       a[i] = v2[i] - v1[i];
 979       b[i] = v3[i] - v1[i];
 980       c[i] = v4[i] - v1[i];
 981     }
 982   int result =
 983     a[0] * (b[1] * c[2] - b[2] * c[1]) -
 984     a[1] * (b[0] * c[2] - b[2] * c[0]) +
 985     a[2] * (b[0] * c[1] - b[1] * c[0]);
 986   return (result > 0) ? result : -result;
 987 }
 988
 989 static void
 990 interpolate_tetrahedron(struct color_interpolation_node *n, uns *p, const uns *c)
 991 {
 992   uns v[4][3];
 993   for (uns i = 0; i < 4; i++)
 994     {
 995       v[i][0] = (c[i] & 0001) ? (1 << COLOR_CONV_OFS) : 0;
 996       v[i][1] = (c[i] & 0010) ? (1 << COLOR_CONV_OFS) : 0;
 997       v[i][2] = (c[i] & 0100) ? (1 << COLOR_CONV_OFS) : 0;
 998       n->ofs[i] =
 999         ((c[i] & 0001) ? 1 : 0) +
1000         ((c[i] & 0010) ? (1 << COLOR_CONV_SIZE) : 0) +
1001         ((c[i] & 0100) ? (1 << (COLOR_CONV_SIZE * 2)) : 0);
1002     }
1003   uns vol = tetrahedron_volume(v[0], v[1], v[2], v[3]);
1004   n->mul[0] = ((tetrahedron_volume(p, v[1], v[2], v[3]) << 8) + (vol >> 1)) / vol;
1005   n->mul[1] = ((tetrahedron_volume(v[0], p, v[2], v[3]) << 8) + (vol >> 1)) / vol;
1006   n->mul[2] = ((tetrahedron_volume(v[0], v[1], p, v[3]) << 8) + (vol >> 1)) / vol;
1007   n->mul[3] = ((tetrahedron_volume(v[0], v[1], v[2], p) << 8) + (vol >> 1)) / vol;
1008   uns j;
1009   for (j = 0; j < 4; j++)
1010     if (n->mul[j])
1011       break;
1012   for (uns i = 0; i < 4; i++)
1013     if (n->mul[i] == 0)
1014       n->ofs[i] = n->ofs[j];
1015 }
1016
1017 static void
1018 interpolation_table_init(void)
1019 {
1020   DBG("Initializing color interpolation table");
1021   struct color_interpolation_node *n = color_interpolation_table =
1022     xmalloc(sizeof(struct color_interpolation_node) << (COLOR_CONV_OFS * 3));
1023   uns p[3];
1024   for (p[2] = 0; p[2] < (1 << COLOR_CONV_OFS); p[2]++)
1025     for (p[1] = 0; p[1] < (1 << COLOR_CONV_OFS); p[1]++)
1026       for (p[0] = 0; p[0] < (1 << COLOR_CONV_OFS); p[0]++)
1027         {
1028           uns index;
1029           static const uns tetrahedra[5][4] = {
1030             {0000, 0001, 0010, 0100},
1031             {0110, 0111, 0100, 0010},
1032             {0101, 0100, 0111, 0001},
1033             {0011, 0010, 0001, 0111},
1034             {0111, 0001, 0010, 0100}};
1035           if (p[0] + p[1] + p[2] <= (1 << COLOR_CONV_OFS))
1036             index = 0;
1037           else if ((1 << COLOR_CONV_OFS) + p[0] <= p[1] + p[2])
1038             index = 1;
1039           else if ((1 << COLOR_CONV_OFS) + p[1] <= p[0] + p[2])
1040             index = 2;
1041           else if ((1 << COLOR_CONV_OFS) + p[2] <= p[0] + p[1])
1042             index = 3;
1043           else
1044             index = 4;
1045           interpolate_tetrahedron(n, p, tetrahedra[index]);
1046           n++;
1047         }
1048 }
1049
1050 typedef void color_conv_func(double dest[3], double src[3]);
1051
1052 static void
1053 conv_grid_init(struct color_grid_node **grid, color_conv_func func)
1054 {
1055   if (*grid)
1056     return;
1057   struct color_grid_node *g = *grid = xmalloc((sizeof(struct color_grid_node)) << (COLOR_CONV_SIZE * 3));
1058   double src[3], dest[3];
1059   for (uns k = 0; k < (1 << COLOR_CONV_SIZE); k++)
1060     {
1061       src[2] = k * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
1062       for (uns j = 0; j < (1 << COLOR_CONV_SIZE); j++)
1063         {
1064           src[1] = j * (255/ (double)((1 << COLOR_CONV_SIZE) - 1));
1065           for (uns i = 0; i < (1 << COLOR_CONV_SIZE); i++)
1066             {
1067               src[0] = i * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
1068               func(dest, src);
1069               g->val[0] = CLAMP(dest[0] + 0.5, 0, 255);
1070               g->val[1] = CLAMP(dest[1] + 0.5, 0, 255);
1071               g->val[2] = CLAMP(dest[2] + 0.5, 0, 255);
1072               g++;
1073             }
1074         }
1075     }
1076 }
1077
1078 static void
1079 srgb_to_luv_func(double dest[3], double src[3])
1080 {
1081   double srgb[3], xyz[3], luv[3];
1082   srgb[0] = src[0] / 255.;
1083   srgb[1] = src[1] / 255.;
1084   srgb[2] = src[2] / 255.;
1085   srgb_to_xyz_exact(xyz, srgb);
1086   xyz_to_luv_exact(luv, xyz);
1087   dest[0] = luv[0] * 2.55;
1088   dest[1] = luv[1] * (2.55 / 4) + 128;
1089   dest[2] = luv[2] * (2.55 / 4) + 128;
1090 }
1091
1092 void
1093 color_conv_init(void)
1094 {
1095   interpolation_table_init();
1096   conv_grid_init(&srgb_to_luv_grid, srgb_to_luv_func);
1097 }
1098
1099 void
1100 color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid)
1101 {
1102   while (count--)
1103     {
1104       color_conv_pixel(dest, src, grid);
1105       dest += 3;
1106       src += 3;
1107     }
1108 }
1109
1110
1111 /**************************** TESTS *******************************/
1112
1113 #ifdef TEST
1114 #include <string.h>
1115
1116 static double
1117 conv_error(u32 color, struct color_grid_node *grid, color_conv_func func)
1118 {
1119   byte src[3], dest[3];
1120   src[0] = color & 255;
1121   src[1] = (color >> 8) & 255;
1122   src[2] = (color >> 16) & 255;
1123   color_conv_pixel(dest, src, grid);
1124   double src2[3], dest2[3];
1125   for (uns i = 0; i < 3; i++)
1126     src2[i] = src[i];
1127   func(dest2, src2);
1128   double err = 0;
1129   for (uns i = 0; i < 3; i++)
1130     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
1131   return err;
1132 }
1133
1134 typedef void test_fn(byte *dest, byte *src);
1135
1136 static double
1137 func_error(u32 color, test_fn test, color_conv_func func)
1138 {
1139   byte src[3], dest[3];
1140   src[0] = color & 255;
1141   src[1] = (color >> 8) & 255;
1142   src[2] = (color >> 16) & 255;
1143   test(dest, src);
1144   double src2[3], dest2[3];
1145   for (uns i = 0; i < 3; i++)
1146     src2[i] = src[i];
1147   func(dest2, src2);
1148   double err = 0;
1149   for (uns i = 0; i < 3; i++)
1150     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
1151   return err;
1152 }
1153
1154 static void
1155 test_grid(byte *name, struct color_grid_node *grid, color_conv_func func)
1156 {
1157   double max_err = 0, sum_err = 0;
1158   uns count = 100000;
1159   for (uns i = 0; i < count; i++)
1160     {
1161       double err = conv_error(random_max(0x1000000), grid, func);
1162       max_err = MAX(err, max_err);
1163       sum_err += err;
1164     }
1165   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
1166   if (max_err > 12)
1167     die("Too large error in %s conversion", name);
1168 }
1169
1170 static void
1171 test_func(byte *name, test_fn test, color_conv_func func)
1172 {
1173   double max_err = 0, sum_err = 0;
1174   uns count = 100000;
1175   for (uns i = 0; i < count; i++)
1176     {
1177       double err = func_error(random_max(0x1000000), test, func);
1178       max_err = MAX(err, max_err);
1179       sum_err += err;
1180     }
1181   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
1182   if (max_err > 12)
1183     die("Too large error in %s conversion", name);
1184 }
1185
1186 int
1187 main(void)
1188 {
1189   srgb_to_luv_init();
1190   test_func("func sRGB -> Luv", srgb_to_luv_pixel, srgb_to_luv_func);
1191   color_conv_init();
1192   test_grid("grid sRGB -> Luv", srgb_to_luv_grid, srgb_to_luv_func);
1193 #ifdef LOCAL_DEBUG
1194 #define CNT 1000000
1195 #define TESTS 10
1196   byte *a = xmalloc(3 * CNT), *b = xmalloc(3 * CNT);
1197   for (uns i = 0; i < 3 * CNT; i++)
1198     a[i] = random_max(256);
1199   init_timer();
1200   for (uns i = 0; i < TESTS; i++)
1201     memcpy(b, a, CNT * 3);
1202   DBG("memcpy time=%d", (uns)get_timer());
1203   init_timer();
1204   for (uns i = 0; i < TESTS; i++)
1205     srgb_to_luv_pixels(b, a, CNT);
1206   DBG("direct time=%d", (uns)get_timer());
1207   init_timer();
1208   for (uns i = 0; i < TESTS; i++)
1209     color_conv_pixels(b, a, CNT, srgb_to_luv_grid);
1210   DBG("grid time=%d", (uns)get_timer());
1211 #endif
1212   return 0;
1213 }
1214 #endif
1215