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