2 * Image Library -- Color Spaces
4 * (c) 2006 Pavel Charvat <pchar@ucw.cz>
6 * This software may be freely distributed and used according to the terms
7 * of the GNU Lesser General Public License.
12 #include "sherlock/sherlock.h"
14 #include "images/images.h"
15 #include "images/color.h"
17 struct color color_black = { .color_space = COLOR_SPACE_GRAYSCALE };
18 struct color color_white = { .c = { 255 }, .color_space = COLOR_SPACE_GRAYSCALE };
21 color_put_grayscale(byte *dest, struct color *color)
23 switch (color->color_space)
25 case COLOR_SPACE_GRAYSCALE:
26 dest[0] = color->c[0];
29 dest[0] = rgb_to_gray_func(color->c[0], color->c[1], color->c[2]);
37 color_put_rgb(byte *dest, struct color *color)
39 switch (color->color_space)
41 case COLOR_SPACE_GRAYSCALE:
42 dest[0] = dest[1] = dest[2] = color->c[0];
45 dest[0] = color->c[0];
46 dest[1] = color->c[1];
47 dest[2] = color->c[2];
55 color_put_color_space(byte *dest, struct color *color, enum color_space color_space)
59 case COLOR_SPACE_GRAYSCALE:
60 color_put_grayscale(dest, color);
63 color_put_rgb(dest, color);
70 /********************* EXACT CONVERSION ROUTINES **********************/
74 srgb_to_xyz_slow(double xyz[3], double srgb[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);
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];
89 xyz_to_srgb_slow(double srgb[3], double xyz[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++)
97 srgb[i] = 1.055 * pow(a[i], 1 / 2.4) - 0.055;
99 srgb[i] = 12.92 * a[i];
104 xyz_to_luv_slow(double luv[3], double xyz[3])
106 double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
108 luv[0] = luv[1] = luv[2] = 0;
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;
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] */
125 luv_to_xyz_slow(double xyz[3], double luv[3])
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)
134 var_y = (var_y - 16 / 116) / 7.787;
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);
141 /***************** OPTIMIZED SRGB -> LUV CONVERSION *********************/
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];
148 srgb_to_luv_init(void)
150 DBG("Initializing sRGB -> Luv table");
151 for (uns i = 0; i < 256; i++)
155 t = pow((t + 0.055) * (1 / 1.055), 2.4);
158 srgb_to_luv_tab1[i] = CLAMP(t * 0xfff + 0.5, 0, 0xfff);
160 for (uns i = 0; i < (9 << SRGB_TO_LUV_TAB2_SIZE); i++)
162 double t = i / (double)((9 << SRGB_TO_LUV_TAB2_SIZE) - 1);
164 t = 1.16 * pow(t, 1 / 3.) - 0.16;
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);
171 for (uns i = 0; i < (20 << SRGB_TO_LUV_TAB3_SIZE); i++)
173 srgb_to_luv_tab3[i] = i ? (13 << (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB3_SIZE)) / i : 0;
178 srgb_to_luv_pixels(byte *dest, byte *src, uns count)
182 srgb_to_luv_pixel(dest, src);
189 /************************ GRID INTERPOLATION ALGORITHM ************************/
191 struct color_grid_node *srgb_to_luv_grid;
192 struct color_interpolation_node *color_interpolation_table;
194 /* Returns volume of a given tetrahedron multiplied by 6 */
196 tetrahedron_volume(uns *v1, uns *v2, uns *v3, uns *v4)
198 int a[3], b[3], c[3];
199 for (uns i = 0; i < 3; i++)
201 a[i] = v2[i] - v1[i];
202 b[i] = v3[i] - v1[i];
203 c[i] = v4[i] - v1[i];
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;
213 interpolate_tetrahedron(struct color_interpolation_node *n, uns *p, const uns *c)
216 for (uns i = 0; i < 4; i++)
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;
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);
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;
232 for (j = 0; j < 4; j++)
235 for (uns i = 0; i < 4; i++)
237 n->ofs[i] = n->ofs[j];
241 interpolation_table_init(void)
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));
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]++)
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))
260 else if ((1 << COLOR_CONV_OFS) + p[0] <= p[1] + p[2])
262 else if ((1 << COLOR_CONV_OFS) + p[1] <= p[0] + p[2])
264 else if ((1 << COLOR_CONV_OFS) + p[2] <= p[0] + p[1])
268 interpolate_tetrahedron(n, p, tetrahedrons[index]);
273 typedef void color_conv_func(double dest[3], double src[3]);
276 conv_grid_init(struct color_grid_node **grid, color_conv_func func)
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++)
284 src[2] = k * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
285 for (uns j = 0; j < (1 << COLOR_CONV_SIZE); j++)
287 src[1] = j * (255/ (double)((1 << COLOR_CONV_SIZE) - 1));
288 for (uns i = 0; i < (1 << COLOR_CONV_SIZE); i++)
290 src[0] = i * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
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);
302 srgb_to_luv_func(double dest[3], double src[3])
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;
316 color_conv_init(void)
318 interpolation_table_init();
319 conv_grid_init(&srgb_to_luv_grid, srgb_to_luv_func);
323 color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid)
327 color_conv_pixel(dest, src, grid);
334 /**************************** TESTS *******************************/
340 conv_error(u32 color, struct color_grid_node *grid, color_conv_func func)
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++)
352 for (uns i = 0; i < 3; i++)
353 err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
357 typedef void test_fn(byte *dest, byte *src);
360 func_error(u32 color, test_fn test, color_conv_func func)
362 byte src[3], dest[3];
363 src[0] = color & 255;
364 src[1] = (color >> 8) & 255;
365 src[2] = (color >> 16) & 255;
367 double src2[3], dest2[3];
368 for (uns i = 0; i < 3; i++)
372 for (uns i = 0; i < 3; i++)
373 err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
378 test_grid(byte *name, struct color_grid_node *grid, color_conv_func func)
380 double max_err = 0, sum_err = 0;
382 for (uns i = 0; i < count; i++)
384 double err = conv_error(random_max(0x1000000), grid, func);
385 max_err = MAX(err, max_err);
388 DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
390 die("Too large error in %s conversion", name);
394 test_func(byte *name, test_fn test, color_conv_func func)
396 double max_err = 0, sum_err = 0;
398 for (uns i = 0; i < count; i++)
400 double err = func_error(random_max(0x1000000), test, func);
401 max_err = MAX(err, max_err);
404 DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
406 die("Too large error in %s conversion", name);
413 test_func("func sRGB -> Luv", srgb_to_luv_pixel, srgb_to_luv_func);
415 test_grid("grid sRGB -> Luv", srgb_to_luv_grid, srgb_to_luv_func);
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);
423 for (uns i = 0; i < TESTS; i++)
424 memcpy(b, a, CNT * 3);
425 DBG("memcpy time=%d", (uns)get_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());
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());