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