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 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, enum color_space 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 /* sRGB to XYZ */
  73 void
  74 srgb_to_xyz_slow(double xyz[3], double srgb[3])
  75 {
  76   double a[3];
  77   for (uns i = 0; i < 3; i++)
  78     if (srgb[i] > 0.04045)
  79       a[i] = pow((srgb[i] + 0.055) * (1 / 1.055), 2.4);
  80     else
  81       a[i] = srgb[i] * (1 / 12.92);
  82   xyz[0] = SRGB_XYZ_XR * a[0] + SRGB_XYZ_XG * a[1] + SRGB_XYZ_XB * a[2];
  83   xyz[1] = SRGB_XYZ_YR * a[0] + SRGB_XYZ_YG * a[1] + SRGB_XYZ_YB * a[2];
  84   xyz[2] = SRGB_XYZ_ZR * a[0] + SRGB_XYZ_ZG * a[1] + SRGB_XYZ_ZB * a[2];
  85 }
  86
  87 /* XYZ to sRGB */
  88 void
  89 xyz_to_srgb_slow(double srgb[3], double xyz[3])
  90 {
  91   double a[3];
  92   a[0] =  3.2406 * xyz[0] + -1.5372 * xyz[1] + -0.4986 * xyz[2];
  93   a[1] = -0.9689 * xyz[0] +  1.8758 * xyz[1] +  0.0415 * xyz[2];
  94   a[2] =  0.0557 * xyz[0] + -0.2040 * xyz[1] +  1.0570 * xyz[2];
  95   for (uns i = 0; i < 3; i++)
  96     if (a[i] > 0.0031308)
  97       srgb[i] = 1.055 * pow(a[i], 1 / 2.4) - 0.055;
  98     else
  99       srgb[i] = 12.92 * a[i];
 100 }
 101
 102 /* XYZ to CIE-Luv */
 103 void
 104 xyz_to_luv_slow(double luv[3], double xyz[3])
 105 {
 106   double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
 107   if (sum < 0.000001)
 108     luv[0] = luv[1] = luv[2] = 0;
 109   else
 110     {
 111       double var_u = 4 * xyz[0] / sum;
 112       double var_v = 9 * xyz[1] / sum;
 113       if (xyz[1] > 0.008856)
 114         luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
 115       else
 116         luv[0] = (116 * 7.787) * xyz[1];
 117      luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
 118      luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
 119      /* intervals [0..100], [-134..220], [-140..122] */
 120    }
 121 }
 122
 123 /* CIE-Luv to XYZ */
 124 void
 125 luv_to_xyz_slow(double xyz[3], double luv[3])
 126 {
 127   double var_u = luv[1] / (13 * luv[0]) + (4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z));
 128   double var_v = luv[2] / (13 * luv[0]) + (9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z));
 129   double var_y = (luv[0] + 16) / 116;
 130   double pow_y = var_y * var_y * var_y;
 131   if (pow_y > 0.008856)
 132     var_y = pow_y;
 133   else
 134     var_y = (var_y - 16 / 116) / 7.787;
 135   xyz[1] = var_y;
 136   xyz[0] = -(9 * xyz[1] * var_u) / ((var_u - 4) * var_v - var_u * var_v);
 137   xyz[2] = (9 * xyz[1] - 15 * var_v * xyz[1] - var_v * xyz[0]) / (3 * var_v);
 138 }
 139
 140
 141 /***************** OPTIMIZED SRGB -> LUV CONVERSION *********************/
 142
 143 u16 srgb_to_luv_tab1[256];
 144 u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE];
 145 u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE];
 146
 147 void
 148 srgb_to_luv_init(void)
 149 {
 150   DBG("Initializing sRGB -> Luv table");
 151   for (uns i = 0; i < 256; i++)
 152     {
 153       double t = i / 255.;
 154       if (t > 0.04045)
 155         t = pow((t + 0.055) * (1 / 1.055), 2.4);
 156       else
 157         t = t * (1 / 12.92);
 158       srgb_to_luv_tab1[i] = CLAMP(t * 0xfff + 0.5, 0, 0xfff);
 159     }
 160   for (uns i = 0; i < (9 << SRGB_TO_LUV_TAB2_SIZE); i++)
 161     {
 162       double t = i / (double)((9 << SRGB_TO_LUV_TAB2_SIZE) - 1);
 163       if (t > 0.008856)
 164         t = 1.16 * pow(t, 1 / 3.) - 0.16;
 165       else
 166         t = (1.16 * 7.787) * t;
 167       srgb_to_luv_tab2[i] =
 168         CLAMP(t * ((1 << SRGB_TO_LUV_TAB2_SCALE) - 1) + 0.5,
 169           0, (1 << SRGB_TO_LUV_TAB2_SCALE) - 1);
 170     }
 171   for (uns i = 0; i < (20 << SRGB_TO_LUV_TAB3_SIZE); i++)
 172     {
 173       srgb_to_luv_tab3[i] = i ? (13 << (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB3_SIZE)) / i : 0;
 174     }
 175 }
 176
 177 void
 178 srgb_to_luv_pixels(byte *dest, byte *src, uns count)
 179 {
 180   while (count--)
 181     {
 182       srgb_to_luv_pixel(dest, src);
 183       dest += 3;
 184       src += 3;
 185     }
 186 }
 187
 188
 189 /************************ GRID INTERPOLATION ALGORITHM ************************/
 190
 191 struct color_grid_node *srgb_to_luv_grid;
 192 struct color_interpolation_node *color_interpolation_table;
 193
 194 /* Returns volume of a given tetrahedron multiplied by 6 */
 195 static inline uns
 196 tetrahedron_volume(uns *v1, uns *v2, uns *v3, uns *v4)
 197 {
 198   int a[3], b[3], c[3];
 199   for (uns i = 0; i < 3; i++)
 200     {
 201       a[i] = v2[i] - v1[i];
 202       b[i] = v3[i] - v1[i];
 203       c[i] = v4[i] - v1[i];
 204     }
 205   int result =
 206     a[0] * (b[1] * c[2] - b[2] * c[1]) -
 207     a[1] * (b[0] * c[2] - b[2] * c[0]) +
 208     a[2] * (b[0] * c[1] - b[1] * c[0]);
 209   return (result > 0) ? result : -result;
 210 }
 211
 212 static void
 213 interpolate_tetrahedron(struct color_interpolation_node *n, uns *p, const uns *c)
 214 {
 215   uns v[4][3];
 216   for (uns i = 0; i < 4; i++)
 217     {
 218       v[i][0] = (c[i] & 0001) ? (1 << COLOR_CONV_OFS) : 0;
 219       v[i][1] = (c[i] & 0010) ? (1 << COLOR_CONV_OFS) : 0;
 220       v[i][2] = (c[i] & 0100) ? (1 << COLOR_CONV_OFS) : 0;
 221       n->ofs[i] =
 222         ((c[i] & 0001) ? 1 : 0) +
 223         ((c[i] & 0010) ? (1 << COLOR_CONV_SIZE) : 0) +
 224         ((c[i] & 0100) ? (1 << (COLOR_CONV_SIZE * 2)) : 0);
 225     }
 226   uns vol = tetrahedron_volume(v[0], v[1], v[2], v[3]);
 227   n->mul[0] = ((tetrahedron_volume(p, v[1], v[2], v[3]) << 8) + (vol >> 1)) / vol;
 228   n->mul[1] = ((tetrahedron_volume(v[0], p, v[2], v[3]) << 8) + (vol >> 1)) / vol;
 229   n->mul[2] = ((tetrahedron_volume(v[0], v[1], p, v[3]) << 8) + (vol >> 1)) / vol;
 230   n->mul[3] = ((tetrahedron_volume(v[0], v[1], v[2], p) << 8) + (vol >> 1)) / vol;
 231   uns j;
 232   for (j = 0; j < 4; j++)
 233     if (n->mul[j])
 234       break;
 235   for (uns i = 0; i < 4; i++)
 236     if (n->mul[i] == 0)
 237       n->ofs[i] = n->ofs[j];
 238 }
 239
 240 static void
 241 interpolation_table_init(void)
 242 {
 243   DBG("Initializing color interpolation table");
 244   struct color_interpolation_node *n = color_interpolation_table =
 245     xmalloc(sizeof(struct color_interpolation_node) << (COLOR_CONV_OFS * 3));
 246   uns p[3];
 247   for (p[2] = 0; p[2] < (1 << COLOR_CONV_OFS); p[2]++)
 248     for (p[1] = 0; p[1] < (1 << COLOR_CONV_OFS); p[1]++)
 249       for (p[0] = 0; p[0] < (1 << COLOR_CONV_OFS); p[0]++)
 250         {
 251           uns index;
 252           static const uns tetrahedrons[5][4] = {
 253             {0000, 0001, 0010, 0100},
 254             {0110, 0111, 0100, 0010},
 255             {0101, 0100, 0111, 0001},
 256             {0011, 0010, 0001, 0111},
 257             {0111, 0001, 0010, 0100}};
 258           if (p[0] + p[1] + p[2] <= (1 << COLOR_CONV_OFS))
 259             index = 0;
 260           else if ((1 << COLOR_CONV_OFS) + p[0] <= p[1] + p[2])
 261             index = 1;
 262           else if ((1 << COLOR_CONV_OFS) + p[1] <= p[0] + p[2])
 263             index = 2;
 264           else if ((1 << COLOR_CONV_OFS) + p[2] <= p[0] + p[1])
 265             index = 3;
 266           else
 267             index = 4;
 268           interpolate_tetrahedron(n, p, tetrahedrons[index]);
 269           n++;
 270         }
 271 }
 272
 273 typedef void color_conv_func(double dest[3], double src[3]);
 274
 275 static void
 276 conv_grid_init(struct color_grid_node **grid, color_conv_func func)
 277 {
 278   if (*grid)
 279     return;
 280   struct color_grid_node *g = *grid = xmalloc((sizeof(struct color_grid_node)) << (COLOR_CONV_SIZE * 3));
 281   double src[3], dest[3];
 282   for (uns k = 0; k < (1 << COLOR_CONV_SIZE); k++)
 283     {
 284       src[2] = k * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
 285       for (uns j = 0; j < (1 << COLOR_CONV_SIZE); j++)
 286         {
 287           src[1] = j * (255/ (double)((1 << COLOR_CONV_SIZE) - 1));
 288           for (uns i = 0; i < (1 << COLOR_CONV_SIZE); i++)
 289             {
 290               src[0] = i * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
 291               func(dest, src);
 292               g->val[0] = CLAMP(dest[0] + 0.5, 0, 255);
 293               g->val[1] = CLAMP(dest[1] + 0.5, 0, 255);
 294               g->val[2] = CLAMP(dest[2] + 0.5, 0, 255);
 295               g++;
 296             }
 297         }
 298     }
 299 }
 300
 301 static void
 302 srgb_to_luv_func(double dest[3], double src[3])
 303 {
 304   double srgb[3], xyz[3], luv[3];
 305   srgb[0] = src[0] / 255.;
 306   srgb[1] = src[1] / 255.;
 307   srgb[2] = src[2] / 255.;
 308   srgb_to_xyz_slow(xyz, srgb);
 309   xyz_to_luv_slow(luv, xyz);
 310   dest[0] = luv[0] * 2.55;
 311   dest[1] = luv[1] * (2.55 / 4) + 128;
 312   dest[2] = luv[2] * (2.55 / 4) + 128;
 313 }
 314
 315 void
 316 color_conv_init(void)
 317 {
 318   interpolation_table_init();
 319   conv_grid_init(&srgb_to_luv_grid, srgb_to_luv_func);
 320 }
 321
 322 void
 323 color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid)
 324 {
 325   while (count--)
 326     {
 327       color_conv_pixel(dest, src, grid);
 328       dest += 3;
 329       src += 3;
 330     }
 331 }
 332
 333
 334 /**************************** TESTS *******************************/
 335
 336 #ifdef TEST
 337 #include <string.h>
 338
 339 static double
 340 conv_error(u32 color, struct color_grid_node *grid, color_conv_func func)
 341 {
 342   byte src[3], dest[3];
 343   src[0] = color & 255;
 344   src[1] = (color >> 8) & 255;
 345   src[2] = (color >> 16) & 255;
 346   color_conv_pixel(dest, src, grid);
 347   double src2[3], dest2[3];
 348   for (uns i = 0; i < 3; i++)
 349     src2[i] = src[i];
 350   func(dest2, src2);
 351   double err = 0;
 352   for (uns i = 0; i < 3; i++)
 353     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
 354   return err;
 355 }
 356
 357 typedef void test_fn(byte *dest, byte *src);
 358
 359 static double
 360 func_error(u32 color, test_fn test, color_conv_func func)
 361 {
 362   byte src[3], dest[3];
 363   src[0] = color & 255;
 364   src[1] = (color >> 8) & 255;
 365   src[2] = (color >> 16) & 255;
 366   test(dest, src);
 367   double src2[3], dest2[3];
 368   for (uns i = 0; i < 3; i++)
 369     src2[i] = src[i];
 370   func(dest2, src2);
 371   double err = 0;
 372   for (uns i = 0; i < 3; i++)
 373     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
 374   return err;
 375 }
 376
 377 static void
 378 test_grid(byte *name, struct color_grid_node *grid, color_conv_func func)
 379 {
 380   double max_err = 0, sum_err = 0;
 381   uns count = 100000;
 382   for (uns i = 0; i < count; i++)
 383     {
 384       double err = conv_error(random_max(0x1000000), grid, func);
 385       max_err = MAX(err, max_err);
 386       sum_err += err;
 387     }
 388   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
 389   if (max_err > 12)
 390     die("Too large error in %s conversion", name);
 391 }
 392
 393 static void
 394 test_func(byte *name, test_fn test, color_conv_func func)
 395 {
 396   double max_err = 0, sum_err = 0;
 397   uns count = 100000;
 398   for (uns i = 0; i < count; i++)
 399     {
 400       double err = func_error(random_max(0x1000000), test, func);
 401       max_err = MAX(err, max_err);
 402       sum_err += err;
 403     }
 404   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
 405   if (max_err > 12)
 406     die("Too large error in %s conversion", name);
 407 }
 408
 409 int
 410 main(void)
 411 {
 412   srgb_to_luv_init();
 413   test_func("func sRGB -> Luv", srgb_to_luv_pixel, srgb_to_luv_func);
 414   color_conv_init();
 415   test_grid("grid sRGB -> Luv", srgb_to_luv_grid, srgb_to_luv_func);
 416 #ifdef LOCAL_DEBUG
 417 #define CNT 1000000
 418 #define TESTS 10
 419   byte *a = xmalloc(3 * CNT), *b = xmalloc(3 * CNT);
 420   for (uns i = 0; i < 3 * CNT; i++)
 421     a[i] = random_max(256);
 422   init_timer();
 423   for (uns i = 0; i < TESTS; i++)
 424     memcpy(b, a, CNT * 3);
 425   DBG("memcpy time=%d", (uns)get_timer());
 426   init_timer();
 427   for (uns i = 0; i < TESTS; i++)
 428     srgb_to_luv_pixels(b, a, CNT);
 429   DBG("direct time=%d", (uns)get_timer());
 430   init_timer();
 431   for (uns i = 0; i < TESTS; i++)
 432     color_conv_pixels(b, a, CNT, srgb_to_luv_grid);
 433   DBG("grid time=%d", (uns)get_timer());
 434 #endif
 435   return 0;
 436 }
 437 #endif
 438