images/color.c

   1 /*
   2  *      Image Library -- Color Spaces
   3  *
   4  *      (c) 2006 Pavel Charvat <pchar@ucw.cz>
   5  *
   6  *      This software may be freely distributed and used according to the terms
   7  *      of the GNU Lesser General Public License.
   8  */
   9
  10 #undef LOCAL_DEBUG
  11
  12 #include "sherlock/sherlock.h"
  13 #include "lib/math.h"
  14 #include "images/images.h"
  15 #include "images/color.h"
  16
  17 #include <string.h>
  18
  19 uns color_space_channels[COLOR_SPACE_MAX] = {
  20   [COLOR_SPACE_UNKNOWN] = 0,
  21   [COLOR_SPACE_UNKNOWN_1] = 1,
  22   [COLOR_SPACE_UNKNOWN_2] = 2,
  23   [COLOR_SPACE_UNKNOWN_3] = 3,
  24   [COLOR_SPACE_UNKNOWN_4] = 4,
  25   [COLOR_SPACE_GRAYSCALE] = 1,
  26   [COLOR_SPACE_RGB] = 3,
  27   [COLOR_SPACE_XYZ] = 3,
  28   [COLOR_SPACE_LAB] = 3,
  29   [COLOR_SPACE_YCBCR] = 3,
  30   [COLOR_SPACE_CMYK] = 4,
  31   [COLOR_SPACE_YCCK] = 4,
  32 };
  33
  34 byte *color_space_name[COLOR_SPACE_MAX] = {
  35   [COLOR_SPACE_UNKNOWN] = "Unknown",
  36   [COLOR_SPACE_UNKNOWN_1] = "1-channel",
  37   [COLOR_SPACE_UNKNOWN_2] = "2-channels",
  38   [COLOR_SPACE_UNKNOWN_3] = "3-channels",
  39   [COLOR_SPACE_UNKNOWN_4] = "4-channels",
  40   [COLOR_SPACE_GRAYSCALE] = "Grayscale",
  41   [COLOR_SPACE_RGB] = "RGB",
  42   [COLOR_SPACE_XYZ] = "XYZ",
  43   [COLOR_SPACE_LAB] = "LAB",
  44   [COLOR_SPACE_YCBCR] = "YCbCr",
  45   [COLOR_SPACE_CMYK] = "CMYK",
  46   [COLOR_SPACE_YCCK] = "YCCK",
  47 };
  48
  49 byte *
  50 color_space_id_to_name(uns id)
  51 {
  52   ASSERT(id < COLOR_SPACE_MAX);
  53   return color_space_name[id];
  54 }
  55
  56 uns
  57 color_space_name_to_id(byte *name)
  58 {
  59   for (uns i = 1; i < COLOR_SPACE_MAX; i++)
  60     if (!strcasecmp(name, color_space_name[i]))
  61       return i;
  62   return 0;
  63 }
  64
  65 struct color color_black = { .color_space = COLOR_SPACE_GRAYSCALE };
  66 struct color color_white = { .c = { 255 }, .color_space = COLOR_SPACE_GRAYSCALE };
  67
  68 inline void
  69 color_put_grayscale(byte *dest, struct color *color)
  70 {
  71   switch (color->color_space)
  72     {
  73       case COLOR_SPACE_GRAYSCALE:
  74         dest[0] = color->c[0];
  75         break;
  76       case COLOR_SPACE_RGB:
  77         dest[0] = rgb_to_gray_func(color->c[0], color->c[1], color->c[2]);
  78         break;
  79       default:
  80         ASSERT(0);
  81     }
  82 }
  83
  84 inline void
  85 color_put_rgb(byte *dest, struct color *color)
  86 {
  87   switch (color->color_space)
  88     {
  89       case COLOR_SPACE_GRAYSCALE:
  90         dest[0] = dest[1] = dest[2] = color->c[0];
  91         break;
  92       case COLOR_SPACE_RGB:
  93         dest[0] = color->c[0];
  94         dest[1] = color->c[1];
  95         dest[2] = color->c[2];
  96         break;
  97       default:
  98         ASSERT(0);
  99     }
 100 }
 101
 102 void
 103 color_put_color_space(byte *dest, struct color *color, uns color_space)
 104 {
 105   switch (color_space)
 106     {
 107       case COLOR_SPACE_GRAYSCALE:
 108         color_put_grayscale(dest, color);
 109         break;
 110       case COLOR_SPACE_RGB:
 111         color_put_rgb(dest, color);
 112         break;
 113       default:
 114         ASSERT(0);
 115     }
 116 }
 117
 118 /********************* EXACT CONVERSION ROUTINES **********************/
 119
 120 /* Reference whites */
 121 #define COLOR_ILLUMINANT_A      0.44757, 0.40744
 122 #define COLOR_ILLUMINANT_B      0.34840, 0.35160
 123 #define COLOR_ILLUMINANT_C      0.31006, 0.31615
 124 #define COLOR_ILLUMINANT_D50    0.34567, 0.35850
 125 #define COLOR_ILLUMINANT_D55    0.33242, 0.34743
 126 #define COLOR_ILLUMINANT_D65    0.31273, 0.32902
 127 #define COLOR_ILLUMINANT_D75    0.29902, 0.31485
 128 #define COLOR_ILLUMINANT_9300K  0.28480, 0.29320
 129 #define COLOR_ILLUMINANT_E      (1./3.), (1./3.)
 130 #define COLOR_ILLUMINANT_F2     0.37207, 0.37512
 131 #define COLOR_ILLUMINANT_F7     0.31285, 0.32918
 132 #define COLOR_ILLUMINANT_F11    0.38054, 0.37691
 133
 134 const double
 135   color_illuminant_d50[2] = {COLOR_ILLUMINANT_D50},
 136   color_illuminant_d65[2] = {COLOR_ILLUMINANT_D65},
 137   color_illuminant_e[2] = {COLOR_ILLUMINANT_E};
 138
 139 /* RGB profiles (many missing) */
 140 const struct color_space_info
 141   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}},
 142   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}},
 143   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}},
 144   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}},
 145   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}};
 146
 147 #define CLIP(x, min, max) (((x) < (min)) ? (min) : ((x) > (max)) ? (max) : (x))
 148
 149 static inline void
 150 clip(double a[3])
 151 {
 152   a[0] = CLIP(a[0], 0, 1);
 153   a[1] = CLIP(a[1], 0, 1);
 154   a[2] = CLIP(a[2], 0, 1);
 155 }
 156
 157 static inline void
 158 correct_gamma_simple(double dest[3], double src[3], const struct color_space_gamma_info *info)
 159 {
 160   dest[0] = pow(src[0], info->simple_gamma);
 161   dest[1] = pow(src[1], info->simple_gamma);
 162   dest[2] = pow(src[2], info->simple_gamma);
 163 }
 164
 165 static inline void
 166 invert_gamma_simple(double dest[3], double src[3], const struct color_space_gamma_info *info)
 167 {
 168   dest[0] = pow(src[0], 1 / info->simple_gamma);
 169   dest[1] = pow(src[1], 1 / info->simple_gamma);
 170   dest[2] = pow(src[2], 1 / info->simple_gamma);
 171 }
 172
 173 static inline void
 174 correct_gamma_detailed(double dest[3], double src[3], const struct color_space_gamma_info *info)
 175 {
 176   for (uns i = 0; i < 3; i++)
 177     if (src[i] > info->transition)
 178       dest[i] = (1 + info->offset) * pow(src[i], info->detailed_gamma) - info->offset;
 179     else
 180       dest[i] = info->slope * src[i];
 181 }
 182
 183 static inline void
 184 invert_gamma_detailed(double dest[3], double src[3], const struct color_space_gamma_info *info)
 185 {
 186   for (uns i = 0; i < 3; i++)
 187     if (src[i] > info->transition * info->slope)
 188       dest[i] = pow((src[i] + info->offset) / (1 + info->offset), 1 / info->detailed_gamma);
 189     else
 190       dest[i] = src[i] / info->slope;
 191 }
 192
 193 static inline void
 194 apply_matrix(double dest[3], double src[3], double matrix[9])
 195 {
 196   dest[0] = src[0] * matrix[0] + src[1] * matrix[1] + src[2] * matrix[2];
 197   dest[1] = src[0] * matrix[3] + src[1] * matrix[4] + src[2] * matrix[5];
 198   dest[2] = src[0] * matrix[6] + src[1] * matrix[7] + src[2] * matrix[8];
 199 }
 200
 201 void
 202 color_invert_matrix(double dest[9], double matrix[9])
 203 {
 204   double *i = dest, *m = matrix;
 205   double a0 = m[4] * m[8] - m[5] * m[7];
 206   double a1 = m[3] * m[8] - m[5] * m[6];
 207   double a2 = m[3] * m[7] - m[4] * m[6];
 208   double d = 1 / (m[0] * a0 - m[1] * a1 + m[2] * a2);
 209   i[0] = d * a0;
 210   i[3] = -d * a1;
 211   i[6] = d * a2;
 212   i[1] = -d * (m[1] * m[8] - m[2] * m[7]);
 213   i[4] = d * (m[0] * m[8] - m[2] * m[6]);
 214   i[7] = -d * (m[0] * m[7] - m[1] * m[6]);
 215   i[2] = d * (m[1] * m[5] - m[2] * m[4]);
 216   i[5] = -d * (m[0] * m[5] - m[2] * m[3]);
 217   i[8] = d * (m[0] * m[4] - m[1] * m[3]);
 218 }
 219
 220 static void
 221 mul_matrices(double r[9], double a[9], double b[9])
 222 {
 223   r[0] = a[0] * b[0] + a[1] * b[3] + a[2] * b[6];
 224   r[1] = a[0] * b[1] + a[1] * b[4] + a[2] * b[7];
 225   r[2] = a[0] * b[2] + a[1] * b[5] + a[2] * b[8];
 226   r[3] = a[3] * b[0] + a[4] * b[3] + a[5] * b[6];
 227   r[4] = a[3] * b[1] + a[4] * b[4] + a[5] * b[7];
 228   r[5] = a[3] * b[2] + a[4] * b[5] + a[5] * b[8];
 229   r[6] = a[6] * b[0] + a[7] * b[3] + a[8] * b[6];
 230   r[7] = a[6] * b[1] + a[7] * b[4] + a[8] * b[7];
 231   r[8] = a[6] * b[2] + a[7] * b[5] + a[8] * b[8];
 232 }
 233
 234 /* computes conversion matrix from a given color space to CIE XYZ */
 235 void
 236 color_compute_color_space_to_xyz_matrix(double matrix[9], const struct color_space_chromacity_info *space)
 237 {
 238   double wX = space->white[0] / space->white[1];
 239   double wZ = (1 - space->white[0] - space->white[1]) / space->white[1];
 240   double a[9], b[9];
 241   a[0] = space->prim1[0]; a[3] = space->prim1[1]; a[6] = 1 - a[0] - a[3];
 242   a[1] = space->prim2[0]; a[4] = space->prim2[1]; a[7] = 1 - a[1] - a[4];
 243   a[2] = space->prim3[0]; a[5] = space->prim3[1]; a[8] = 1 - a[2] - a[5];
 244   color_invert_matrix(b, a);
 245   double ra = wX * b[0] + b[1] + wZ * b[2];
 246   double rb = wX * b[3] + b[4] + wZ * b[5];
 247   double rc = wX * b[6] + b[7] + wZ * b[8];
 248   matrix[0] = a[0] * ra;
 249   matrix[1] = a[1] * rb;
 250   matrix[2] = a[2] * rc;
 251   matrix[3] = a[3] * ra;
 252   matrix[4] = a[4] * rb;
 253   matrix[5] = a[5] * rc;
 254   matrix[6] = a[6] * ra;
 255   matrix[7] = a[7] * rb;
 256   matrix[8] = a[8] * rc;
 257 }
 258
 259 /* computes matrix to join transformations with different reference whites */
 260 void
 261 color_compute_bradford_matrix(double matrix[9], const double source[2], const double dest[2])
 262 {
 263   /* cone response matrix and its inversion */
 264   static double r[9] = {
 265     0.8951, 0.2664, -0.1614,
 266     -0.7502, 1.7135, 0.0367,
 267     0.0389, -0.0685, 1.0296};
 268   //static double i[9] = {0.9870, -0.1471, 0.1600, 0.4323, 0.5184, 0.0493, -0.0085, 0.0400, 0.9685};
 269   double i[9];
 270   color_invert_matrix(i, r);
 271   double aX = source[0] / source[1];
 272   double aZ = (1 - source[0] - source[1]) / source[1];
 273   double bX = dest[0] / dest[1];
 274   double bZ = (1 - dest[0] - dest[1]) / dest[1];
 275   double x = (r[0] * bX + r[1] + r[2] * bZ) / (r[0] * aX + r[1] + r[2] * aZ);
 276   double y = (r[3] * bX + r[4] + r[5] * bZ) / (r[3] * aX + r[4] + r[5] * aZ);
 277   double z = (r[6] * bX + r[7] + r[8] * bZ) / (r[6] * aX + r[7] + r[8] * aZ);
 278   double m[9];
 279   m[0] = i[0] * x; m[1] = i[1] * y; m[2] = i[2] * z;
 280   m[3] = i[3] * x; m[4] = i[4] * y; m[5] = i[5] * z;
 281   m[6] = i[6] * x; m[7] = i[7] * y; m[8] = i[8] * z;
 282   mul_matrices(matrix, m, r);
 283 }
 284
 285 void
 286 color_compute_color_spaces_conversion_matrix(double matrix[9], const struct color_space_chromacity_info *src, const struct color_space_chromacity_info *dest)
 287 {
 288   double a_to_xyz[9], b_to_xyz[9], xyz_to_b[9], bradford[9], m[9];
 289   color_compute_color_space_to_xyz_matrix(a_to_xyz, src);
 290   color_compute_color_space_to_xyz_matrix(b_to_xyz, dest);
 291   color_invert_matrix(xyz_to_b, b_to_xyz);
 292   if (src->white[0] == dest->white[0] && src->white[1] == dest->white[1])
 293     mul_matrices(matrix, a_to_xyz, xyz_to_b);
 294   else
 295     {
 296       color_compute_bradford_matrix(bradford, src->white, dest->white);
 297       mul_matrices(m, a_to_xyz, bradford);
 298       mul_matrices(matrix, m, xyz_to_b);
 299     }
 300 }
 301
 302 /* sRGB to XYZ */
 303 void
 304 srgb_to_xyz_exact(double xyz[3], double srgb[3])
 305 {
 306   static double matrix[9] = {
 307     0.41248031, 0.35756952, 0.18043951,
 308     0.21268516, 0.71513904, 0.07217580,
 309     0.01933501, 0.11918984, 0.95031473};
 310   double srgb_lin[3];
 311   invert_gamma_detailed(srgb_lin, srgb, &color_srgb_info.gamma);
 312   apply_matrix(xyz, srgb_lin, matrix);
 313   xyz_to_srgb_exact(srgb_lin, xyz);
 314 }
 315
 316 /* XYZ to sRGB */
 317 void
 318 xyz_to_srgb_exact(double srgb[3], double xyz[3])
 319 {
 320   static double matrix[9] = {
 321      3.24026666, -1.53704957, -0.49850256,
 322     -0.96928381,  1.87604525,  0.04155678,
 323      0.05564281, -0.20402363,  1.05721334};
 324   double srgb_lin[3];
 325   apply_matrix(srgb_lin, xyz, matrix);
 326   clip(srgb_lin);
 327   correct_gamma_detailed(srgb, srgb_lin, &color_srgb_info.gamma);
 328 }
 329
 330 /* XYZ to CIE-Luv */
 331 void
 332 xyz_to_luv_exact(double luv[3], double xyz[3])
 333 {
 334   double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
 335   if (sum < 0.000001)
 336     luv[0] = luv[1] = luv[2] = 0;
 337   else
 338     {
 339       double var_u = 4 * xyz[0] / sum;
 340       double var_v = 9 * xyz[1] / sum;
 341       if (xyz[1] > 0.008856)
 342         luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
 343       else
 344         luv[0] = (116 * 7.787) * xyz[1];
 345      luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
 346      luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
 347      /* intervals [0..100], [-134..220], [-140..122] */
 348    }
 349 }
 350
 351 /* CIE-Luv to XYZ */
 352 void
 353 luv_to_xyz_exact(double xyz[3], double luv[3])
 354 {
 355   double var_u = luv[1] / (13 * luv[0]) + (4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z));
 356   double var_v = luv[2] / (13 * luv[0]) + (9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z));
 357   double var_y = (luv[0] + 16) / 116;
 358   double pow_y = var_y * var_y * var_y;
 359   if (pow_y > 0.008856)
 360     var_y = pow_y;
 361   else
 362     var_y = (var_y - 16 / 116) / 7.787;
 363   xyz[1] = var_y;
 364   xyz[0] = -(9 * xyz[1] * var_u) / ((var_u - 4) * var_v - var_u * var_v);
 365   xyz[2] = (9 * xyz[1] - 15 * var_v * xyz[1] - var_v * xyz[0]) / (3 * var_v);
 366 }
 367
 368
 369 /***************** OPTIMIZED SRGB -> LUV CONVERSION *********************/
 370
 371 u16 srgb_to_luv_tab1[256];
 372 u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE];
 373 u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE];
 374
 375 void
 376 srgb_to_luv_init(void)
 377 {
 378   DBG("Initializing sRGB -> Luv table");
 379   for (uns i = 0; i < 256; i++)
 380     {
 381       double t = i / 255.;
 382       if (t > 0.04045)
 383         t = pow((t + 0.055) * (1 / 1.055), 2.4);
 384       else
 385         t = t * (1 / 12.92);
 386       srgb_to_luv_tab1[i] = CLAMP(t * 0xfff + 0.5, 0, 0xfff);
 387     }
 388   for (uns i = 0; i < (9 << SRGB_TO_LUV_TAB2_SIZE); i++)
 389     {
 390       double t = i / (double)((9 << SRGB_TO_LUV_TAB2_SIZE) - 1);
 391       if (t > 0.008856)
 392         t = 1.16 * pow(t, 1 / 3.) - 0.16;
 393       else
 394         t = (1.16 * 7.787) * t;
 395       srgb_to_luv_tab2[i] =
 396         CLAMP(t * ((1 << SRGB_TO_LUV_TAB2_SCALE) - 1) + 0.5,
 397           0, (1 << SRGB_TO_LUV_TAB2_SCALE) - 1);
 398     }
 399   for (uns i = 0; i < (20 << SRGB_TO_LUV_TAB3_SIZE); i++)
 400     {
 401       srgb_to_luv_tab3[i] = i ? (13 << (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB3_SIZE)) / i : 0;
 402     }
 403 }
 404
 405 void
 406 srgb_to_luv_pixels(byte *dest, byte *src, uns count)
 407 {
 408   while (count--)
 409     {
 410       srgb_to_luv_pixel(dest, src);
 411       dest += 3;
 412       src += 3;
 413     }
 414 }
 415
 416
 417 /************************ GRID INTERPOLATION ALGORITHM ************************/
 418
 419 struct color_grid_node *srgb_to_luv_grid;
 420 struct color_interpolation_node *color_interpolation_table;
 421
 422 /* Returns volume of a given tetrahedron multiplied by 6 */
 423 static inline uns
 424 tetrahedron_volume(uns *v1, uns *v2, uns *v3, uns *v4)
 425 {
 426   int a[3], b[3], c[3];
 427   for (uns i = 0; i < 3; i++)
 428     {
 429       a[i] = v2[i] - v1[i];
 430       b[i] = v3[i] - v1[i];
 431       c[i] = v4[i] - v1[i];
 432     }
 433   int result =
 434     a[0] * (b[1] * c[2] - b[2] * c[1]) -
 435     a[1] * (b[0] * c[2] - b[2] * c[0]) +
 436     a[2] * (b[0] * c[1] - b[1] * c[0]);
 437   return (result > 0) ? result : -result;
 438 }
 439
 440 static void
 441 interpolate_tetrahedron(struct color_interpolation_node *n, uns *p, const uns *c)
 442 {
 443   uns v[4][3];
 444   for (uns i = 0; i < 4; i++)
 445     {
 446       v[i][0] = (c[i] & 0001) ? (1 << COLOR_CONV_OFS) : 0;
 447       v[i][1] = (c[i] & 0010) ? (1 << COLOR_CONV_OFS) : 0;
 448       v[i][2] = (c[i] & 0100) ? (1 << COLOR_CONV_OFS) : 0;
 449       n->ofs[i] =
 450         ((c[i] & 0001) ? 1 : 0) +
 451         ((c[i] & 0010) ? (1 << COLOR_CONV_SIZE) : 0) +
 452         ((c[i] & 0100) ? (1 << (COLOR_CONV_SIZE * 2)) : 0);
 453     }
 454   uns vol = tetrahedron_volume(v[0], v[1], v[2], v[3]);
 455   n->mul[0] = ((tetrahedron_volume(p, v[1], v[2], v[3]) << 8) + (vol >> 1)) / vol;
 456   n->mul[1] = ((tetrahedron_volume(v[0], p, v[2], v[3]) << 8) + (vol >> 1)) / vol;
 457   n->mul[2] = ((tetrahedron_volume(v[0], v[1], p, v[3]) << 8) + (vol >> 1)) / vol;
 458   n->mul[3] = ((tetrahedron_volume(v[0], v[1], v[2], p) << 8) + (vol >> 1)) / vol;
 459   uns j;
 460   for (j = 0; j < 4; j++)
 461     if (n->mul[j])
 462       break;
 463   for (uns i = 0; i < 4; i++)
 464     if (n->mul[i] == 0)
 465       n->ofs[i] = n->ofs[j];
 466 }
 467
 468 static void
 469 interpolation_table_init(void)
 470 {
 471   DBG("Initializing color interpolation table");
 472   struct color_interpolation_node *n = color_interpolation_table =
 473     xmalloc(sizeof(struct color_interpolation_node) << (COLOR_CONV_OFS * 3));
 474   uns p[3];
 475   for (p[2] = 0; p[2] < (1 << COLOR_CONV_OFS); p[2]++)
 476     for (p[1] = 0; p[1] < (1 << COLOR_CONV_OFS); p[1]++)
 477       for (p[0] = 0; p[0] < (1 << COLOR_CONV_OFS); p[0]++)
 478         {
 479           uns index;
 480           static const uns tetrahedra[5][4] = {
 481             {0000, 0001, 0010, 0100},
 482             {0110, 0111, 0100, 0010},
 483             {0101, 0100, 0111, 0001},
 484             {0011, 0010, 0001, 0111},
 485             {0111, 0001, 0010, 0100}};
 486           if (p[0] + p[1] + p[2] <= (1 << COLOR_CONV_OFS))
 487             index = 0;
 488           else if ((1 << COLOR_CONV_OFS) + p[0] <= p[1] + p[2])
 489             index = 1;
 490           else if ((1 << COLOR_CONV_OFS) + p[1] <= p[0] + p[2])
 491             index = 2;
 492           else if ((1 << COLOR_CONV_OFS) + p[2] <= p[0] + p[1])
 493             index = 3;
 494           else
 495             index = 4;
 496           interpolate_tetrahedron(n, p, tetrahedra[index]);
 497           n++;
 498         }
 499 }
 500
 501 typedef void color_conv_func(double dest[3], double src[3]);
 502
 503 static void
 504 conv_grid_init(struct color_grid_node **grid, color_conv_func func)
 505 {
 506   if (*grid)
 507     return;
 508   struct color_grid_node *g = *grid = xmalloc((sizeof(struct color_grid_node)) << (COLOR_CONV_SIZE * 3));
 509   double src[3], dest[3];
 510   for (uns k = 0; k < (1 << COLOR_CONV_SIZE); k++)
 511     {
 512       src[2] = k * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
 513       for (uns j = 0; j < (1 << COLOR_CONV_SIZE); j++)
 514         {
 515           src[1] = j * (255/ (double)((1 << COLOR_CONV_SIZE) - 1));
 516           for (uns i = 0; i < (1 << COLOR_CONV_SIZE); i++)
 517             {
 518               src[0] = i * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
 519               func(dest, src);
 520               g->val[0] = CLAMP(dest[0] + 0.5, 0, 255);
 521               g->val[1] = CLAMP(dest[1] + 0.5, 0, 255);
 522               g->val[2] = CLAMP(dest[2] + 0.5, 0, 255);
 523               g++;
 524             }
 525         }
 526     }
 527 }
 528
 529 static void
 530 srgb_to_luv_func(double dest[3], double src[3])
 531 {
 532   double srgb[3], xyz[3], luv[3];
 533   srgb[0] = src[0] / 255.;
 534   srgb[1] = src[1] / 255.;
 535   srgb[2] = src[2] / 255.;
 536   srgb_to_xyz_exact(xyz, srgb);
 537   xyz_to_luv_exact(luv, xyz);
 538   dest[0] = luv[0] * 2.55;
 539   dest[1] = luv[1] * (2.55 / 4) + 128;
 540   dest[2] = luv[2] * (2.55 / 4) + 128;
 541 }
 542
 543 void
 544 color_conv_init(void)
 545 {
 546   interpolation_table_init();
 547   conv_grid_init(&srgb_to_luv_grid, srgb_to_luv_func);
 548 }
 549
 550 void
 551 color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid)
 552 {
 553   while (count--)
 554     {
 555       color_conv_pixel(dest, src, grid);
 556       dest += 3;
 557       src += 3;
 558     }
 559 }
 560
 561
 562 /**************************** TESTS *******************************/
 563
 564 #ifdef TEST
 565 #include <string.h>
 566
 567 static double
 568 conv_error(u32 color, struct color_grid_node *grid, color_conv_func func)
 569 {
 570   byte src[3], dest[3];
 571   src[0] = color & 255;
 572   src[1] = (color >> 8) & 255;
 573   src[2] = (color >> 16) & 255;
 574   color_conv_pixel(dest, src, grid);
 575   double src2[3], dest2[3];
 576   for (uns i = 0; i < 3; i++)
 577     src2[i] = src[i];
 578   func(dest2, src2);
 579   double err = 0;
 580   for (uns i = 0; i < 3; i++)
 581     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
 582   return err;
 583 }
 584
 585 typedef void test_fn(byte *dest, byte *src);
 586
 587 static double
 588 func_error(u32 color, test_fn test, color_conv_func func)
 589 {
 590   byte src[3], dest[3];
 591   src[0] = color & 255;
 592   src[1] = (color >> 8) & 255;
 593   src[2] = (color >> 16) & 255;
 594   test(dest, src);
 595   double src2[3], dest2[3];
 596   for (uns i = 0; i < 3; i++)
 597     src2[i] = src[i];
 598   func(dest2, src2);
 599   double err = 0;
 600   for (uns i = 0; i < 3; i++)
 601     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
 602   return err;
 603 }
 604
 605 static void
 606 test_grid(byte *name, struct color_grid_node *grid, color_conv_func func)
 607 {
 608   double max_err = 0, sum_err = 0;
 609   uns count = 100000;
 610   for (uns i = 0; i < count; i++)
 611     {
 612       double err = conv_error(random_max(0x1000000), grid, func);
 613       max_err = MAX(err, max_err);
 614       sum_err += err;
 615     }
 616   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
 617   if (max_err > 12)
 618     die("Too large error in %s conversion", name);
 619 }
 620
 621 static void
 622 test_func(byte *name, test_fn test, color_conv_func func)
 623 {
 624   double max_err = 0, sum_err = 0;
 625   uns count = 100000;
 626   for (uns i = 0; i < count; i++)
 627     {
 628       double err = func_error(random_max(0x1000000), test, func);
 629       max_err = MAX(err, max_err);
 630       sum_err += err;
 631     }
 632   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
 633   if (max_err > 12)
 634     die("Too large error in %s conversion", name);
 635 }
 636
 637 int
 638 main(void)
 639 {
 640   srgb_to_luv_init();
 641   test_func("func sRGB -> Luv", srgb_to_luv_pixel, srgb_to_luv_func);
 642   color_conv_init();
 643   test_grid("grid sRGB -> Luv", srgb_to_luv_grid, srgb_to_luv_func);
 644 #ifdef LOCAL_DEBUG
 645 #define CNT 1000000
 646 #define TESTS 10
 647   byte *a = xmalloc(3 * CNT), *b = xmalloc(3 * CNT);
 648   for (uns i = 0; i < 3 * CNT; i++)
 649     a[i] = random_max(256);
 650   init_timer();
 651   for (uns i = 0; i < TESTS; i++)
 652     memcpy(b, a, CNT * 3);
 653   DBG("memcpy time=%d", (uns)get_timer());
 654   init_timer();
 655   for (uns i = 0; i < TESTS; i++)
 656     srgb_to_luv_pixels(b, a, CNT);
 657   DBG("direct time=%d", (uns)get_timer());
 658   init_timer();
 659   for (uns i = 0; i < TESTS; i++)
 660     color_conv_pixels(b, a, CNT, srgb_to_luv_grid);
 661   DBG("grid time=%d", (uns)get_timer());
 662 #endif
 663   return 0;
 664 }
 665 #endif
 666