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