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