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/color.h"
  15
  16
  17 /********************* EXACT CONVERSION ROUTINES **********************/
  18
  19 /* sRGB to XYZ */
  20 void
  21 srgb_to_xyz_slow(double xyz[3], double srgb[3])
  22 {
  23   double a[3];
  24   for (uns i = 0; i < 3; i++)
  25     if (srgb[i] > 0.04045)
  26       a[i] = pow((srgb[i] + 0.055) * (1 / 1.055), 2.4);
  27     else
  28       a[i] = srgb[i] * (1 / 12.92);
  29   xyz[0] = SRGB_XYZ_XR * a[0] + SRGB_XYZ_XG * a[1] + SRGB_XYZ_XB * a[2];
  30   xyz[1] = SRGB_XYZ_YR * a[0] + SRGB_XYZ_YG * a[1] + SRGB_XYZ_YB * a[2];
  31   xyz[2] = SRGB_XYZ_ZR * a[0] + SRGB_XYZ_ZG * a[1] + SRGB_XYZ_ZB * a[2];
  32 }
  33
  34 /* XYZ to CIE-Luv */
  35 void
  36 xyz_to_luv_slow(double luv[3], double xyz[3])
  37 {
  38   double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
  39   if (sum < 0.000001)
  40     luv[0] = luv[1] = luv[2] = 0;
  41   else
  42     {
  43       double var_u = 4 * xyz[0] / sum;
  44       double var_v = 9 * xyz[1] / sum;
  45       if (xyz[1] > 0.008856)
  46         luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
  47       else
  48         luv[0] = (116 * 7.787) * xyz[1];
  49      luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
  50      luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
  51      /* intervals [0..100], [-134..220], [-140..122] */
  52    }
  53 }
  54
  55
  56 /***************** OPTIMIZED SRGB -> LUV CONVERSION *********************/
  57
  58 u16 srgb_to_luv_tab1[256];
  59 u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE];
  60 u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE];
  61
  62 void
  63 srgb_to_luv_init(void)
  64 {
  65   DBG("Initializing sRGB -> Luv table");
  66   for (uns i = 0; i < 256; i++)
  67     {
  68       double t = i / 255.;
  69       if (t > 0.04045)
  70         t = pow((t + 0.055) * (1 / 1.055), 2.4);
  71       else
  72         t = t * (1 / 12.92);
  73       srgb_to_luv_tab1[i] = CLAMP(t * 0xfff + 0.5, 0, 0xfff);
  74     }
  75   for (uns i = 0; i < (9 << SRGB_TO_LUV_TAB2_SIZE); i++)
  76     {
  77       double t = i / (double)((9 << SRGB_TO_LUV_TAB2_SIZE) - 1);
  78       if (t > 0.008856)
  79         t = 1.16 * pow(t, 1 / 3.) - 0.16;
  80       else
  81         t = (1.16 * 7.787) * t;
  82       srgb_to_luv_tab2[i] =
  83         CLAMP(t * ((1 << SRGB_TO_LUV_TAB2_SCALE) - 1) + 0.5,
  84           0, (1 << SRGB_TO_LUV_TAB2_SCALE) - 1);
  85     }
  86   for (uns i = 0; i < (20 << SRGB_TO_LUV_TAB3_SIZE); i++)
  87     {
  88       srgb_to_luv_tab3[i] = i ? (13 << (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB3_SIZE)) / i : 0;
  89     }
  90 }
  91
  92 void
  93 srgb_to_luv_pixels(byte *dest, byte *src, uns count)
  94 {
  95   while (count--)
  96     {
  97       srgb_to_luv_pixel(dest, src);
  98       dest += 3;
  99       src += 3;
 100     }
 101 }
 102
 103
 104 /************************ GRID INTERPOLATION ALGORITHM ************************/
 105
 106 struct color_grid_node *srgb_to_luv_grid;
 107 struct color_interpolation_node *color_interpolation_table;
 108
 109 /* Returns volume of a given tetrahedron multiplied by 6 */
 110 static inline uns
 111 tetrahedron_volume(uns *v1, uns *v2, uns *v3, uns *v4)
 112 {
 113   int a[3], b[3], c[3];
 114   for (uns i = 0; i < 3; i++)
 115     {
 116       a[i] = v2[i] - v1[i];
 117       b[i] = v3[i] - v1[i];
 118       c[i] = v4[i] - v1[i];
 119     }
 120   int result =
 121     a[0] * (b[1] * c[2] - b[2] * c[1]) -
 122     a[1] * (b[0] * c[2] - b[2] * c[0]) +
 123     a[2] * (b[0] * c[1] - b[1] * c[0]);
 124   return (result > 0) ? result : -result;
 125 }
 126
 127 static void
 128 interpolate_tetrahedron(struct color_interpolation_node *n, uns *p, const uns *c)
 129 {
 130   uns v[4][3];
 131   for (uns i = 0; i < 4; i++)
 132     {
 133       v[i][0] = (c[i] & 0001) ? (1 << COLOR_CONV_OFS) : 0;
 134       v[i][1] = (c[i] & 0010) ? (1 << COLOR_CONV_OFS) : 0;
 135       v[i][2] = (c[i] & 0100) ? (1 << COLOR_CONV_OFS) : 0;
 136       n->ofs[i] =
 137         ((c[i] & 0001) ? 1 : 0) +
 138         ((c[i] & 0010) ? (1 << COLOR_CONV_SIZE) : 0) +
 139         ((c[i] & 0100) ? (1 << (COLOR_CONV_SIZE * 2)) : 0);
 140     }
 141   uns vol = tetrahedron_volume(v[0], v[1], v[2], v[3]);
 142   n->mul[0] = ((tetrahedron_volume(p, v[1], v[2], v[3]) << 8) + (vol >> 1)) / vol;
 143   n->mul[1] = ((tetrahedron_volume(v[0], p, v[2], v[3]) << 8) + (vol >> 1)) / vol;
 144   n->mul[2] = ((tetrahedron_volume(v[0], v[1], p, v[3]) << 8) + (vol >> 1)) / vol;
 145   n->mul[3] = ((tetrahedron_volume(v[0], v[1], v[2], p) << 8) + (vol >> 1)) / vol;
 146   uns j;
 147   for (j = 0; j < 4; j++)
 148     if (n->mul[j])
 149       break;
 150   for (uns i = 0; i < 4; i++)
 151     if (n->mul[i] == 0)
 152       n->ofs[i] = n->ofs[j];
 153 }
 154
 155 static void
 156 interpolation_table_init(void)
 157 {
 158   DBG("Initializing color interpolation table");
 159   struct color_interpolation_node *n = color_interpolation_table =
 160     xmalloc(sizeof(struct color_interpolation_node) << (COLOR_CONV_OFS * 3));
 161   uns p[3];
 162   for (p[2] = 0; p[2] < (1 << COLOR_CONV_OFS); p[2]++)
 163     for (p[1] = 0; p[1] < (1 << COLOR_CONV_OFS); p[1]++)
 164       for (p[0] = 0; p[0] < (1 << COLOR_CONV_OFS); p[0]++)
 165         {
 166           uns index;
 167           static const uns tetrahedrons[5][4] = {
 168             {0000, 0001, 0010, 0100},
 169             {0110, 0111, 0100, 0010},
 170             {0101, 0100, 0111, 0001},
 171             {0011, 0010, 0001, 0111},
 172             {0111, 0001, 0010, 0100}};
 173           if (p[0] + p[1] + p[2] <= (1 << COLOR_CONV_OFS))
 174             index = 0;
 175           else if ((1 << COLOR_CONV_OFS) + p[0] <= p[1] + p[2])
 176             index = 1;
 177           else if ((1 << COLOR_CONV_OFS) + p[1] <= p[0] + p[2])
 178             index = 2;
 179           else if ((1 << COLOR_CONV_OFS) + p[2] <= p[0] + p[1])
 180             index = 3;
 181           else
 182             index = 4;
 183           interpolate_tetrahedron(n, p, tetrahedrons[index]);
 184           n++;
 185         }
 186 }
 187
 188 typedef void color_conv_func(double dest[3], double src[3]);
 189
 190 static void
 191 conv_grid_init(struct color_grid_node **grid, color_conv_func func)
 192 {
 193   if (*grid)
 194     return;
 195   struct color_grid_node *g = *grid = xmalloc((sizeof(struct color_grid_node)) << (COLOR_CONV_SIZE * 3));
 196   double src[3], dest[3];
 197   for (uns k = 0; k < (1 << COLOR_CONV_SIZE); k++)
 198     {
 199       src[2] = k * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
 200       for (uns j = 0; j < (1 << COLOR_CONV_SIZE); j++)
 201         {
 202           src[1] = j * (255/ (double)((1 << COLOR_CONV_SIZE) - 1));
 203           for (uns i = 0; i < (1 << COLOR_CONV_SIZE); i++)
 204             {
 205               src[0] = i * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
 206               func(dest, src);
 207               g->val[0] = CLAMP(dest[0] + 0.5, 0, 255);
 208               g->val[1] = CLAMP(dest[1] + 0.5, 0, 255);
 209               g->val[2] = CLAMP(dest[2] + 0.5, 0, 255);
 210               g++;
 211             }
 212         }
 213     }
 214 }
 215
 216 static void
 217 srgb_to_luv_func(double dest[3], double src[3])
 218 {
 219   double srgb[3], xyz[3], luv[3];
 220   srgb[0] = src[0] / 255.;
 221   srgb[1] = src[1] / 255.;
 222   srgb[2] = src[2] / 255.;
 223   srgb_to_xyz_slow(xyz, srgb);
 224   xyz_to_luv_slow(luv, xyz);
 225   dest[0] = luv[0] * 2.55;
 226   dest[1] = luv[1] * (2.55 / 4) + 128;
 227   dest[2] = luv[2] * (2.55 / 4) + 128;
 228 }
 229
 230 void
 231 color_conv_init(void)
 232 {
 233   interpolation_table_init();
 234   conv_grid_init(&srgb_to_luv_grid, srgb_to_luv_func);
 235 }
 236
 237 void
 238 color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid)
 239 {
 240   while (count--)
 241     {
 242       color_conv_pixel(dest, src, grid);
 243       dest += 3;
 244       src += 3;
 245     }
 246 }
 247
 248
 249 /**************************** TESTS *******************************/
 250
 251 #ifdef TEST
 252 #include <string.h>
 253
 254 static double
 255 conv_error(u32 color, struct color_grid_node *grid, color_conv_func func)
 256 {
 257   byte src[3], dest[3];
 258   src[0] = color & 255;
 259   src[1] = (color >> 8) & 255;
 260   src[2] = (color >> 16) & 255;
 261   color_conv_pixel(dest, src, grid);
 262   double src2[3], dest2[3];
 263   for (uns i = 0; i < 3; i++)
 264     src2[i] = src[i];
 265   func(dest2, src2);
 266   double err = 0;
 267   for (uns i = 0; i < 3; i++)
 268     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
 269   return err;
 270 }
 271
 272 typedef void test_fn(byte *dest, byte *src);
 273
 274 static double
 275 func_error(u32 color, test_fn test, color_conv_func func)
 276 {
 277   byte src[3], dest[3];
 278   src[0] = color & 255;
 279   src[1] = (color >> 8) & 255;
 280   src[2] = (color >> 16) & 255;
 281   test(dest, src);
 282   double src2[3], dest2[3];
 283   for (uns i = 0; i < 3; i++)
 284     src2[i] = src[i];
 285   func(dest2, src2);
 286   double err = 0;
 287   for (uns i = 0; i < 3; i++)
 288     err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
 289   return err;
 290 }
 291
 292 static void
 293 test_grid(byte *name, struct color_grid_node *grid, color_conv_func func)
 294 {
 295   double max_err = 0, sum_err = 0;
 296   uns count = 100000;
 297   for (uns i = 0; i < count; i++)
 298     {
 299       double err = conv_error(random_max(0x1000000), grid, func);
 300       max_err = MAX(err, max_err);
 301       sum_err += err;
 302     }
 303   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
 304   if (max_err > 12)
 305     die("Too large error in %s conversion", name);
 306 }
 307
 308 static void
 309 test_func(byte *name, test_fn test, color_conv_func func)
 310 {
 311   double max_err = 0, sum_err = 0;
 312   uns count = 100000;
 313   for (uns i = 0; i < count; i++)
 314     {
 315       double err = func_error(random_max(0x1000000), test, func);
 316       max_err = MAX(err, max_err);
 317       sum_err += err;
 318     }
 319   DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
 320   if (max_err > 12)
 321     die("Too large error in %s conversion", name);
 322 }
 323
 324 int
 325 main(void)
 326 {
 327   srgb_to_luv_init();
 328   test_func("func sRGB -> Luv", srgb_to_luv_pixel, srgb_to_luv_func);
 329   color_conv_init();
 330   test_grid("grid sRGB -> Luv", srgb_to_luv_grid, srgb_to_luv_func);
 331 #ifdef LOCAL_DEBUG
 332 #define CNT 1000000
 333 #define TESTS 10
 334   byte *a = xmalloc(3 * CNT), *b = xmalloc(3 * CNT);
 335   for (uns i = 0; i < 3 * CNT; i++)
 336     a[i] = random_max(256);
 337   init_timer();
 338   for (uns i = 0; i < TESTS; i++)
 339     memcpy(b, a, CNT * 3);
 340   DBG("memcpy time=%d", (uns)get_timer());
 341   init_timer();
 342   for (uns i = 0; i < TESTS; i++)
 343     srgb_to_luv_pixels(b, a, CNT);
 344   DBG("direct time=%d", (uns)get_timer());
 345   init_timer();
 346   for (uns i = 0; i < TESTS; i++)
 347     color_conv_pixels(b, a, CNT, srgb_to_luv_grid);
 348   DBG("grid time=%d", (uns)get_timer());
 349 #endif
 350   return 0;
 351 }
 352 #endif
 353