From: Pavel Charvat Date: Mon, 1 May 2006 08:56:40 +0000 (+0200) Subject: Faster sRGB -> Luv colorspace conversion. X-Git-Tag: holmes-import~648 X-Git-Url: http://mj.ucw.cz/gitweb/?a=commitdiff_plain;h=d1fd4b0760bb5d06d7eed80f759de005d2ad44f9;p=libucw.git Faster sRGB -> Luv colorspace conversion. I tried to implement two different algorithms. The first one directly evaluates the conversion formula while the second one is general interpolation of continous 3-dimensional function f: [0,255]^3 -> [0,255]^3 given by a grid of sample values. --- diff --git a/images/color.c b/images/color.c new file mode 100644 index 00000000..457e2e2c --- /dev/null +++ b/images/color.c @@ -0,0 +1,343 @@ +/* + * Image Library -- Color Spaces + * + * (c) 2006 Pavel Charvat + * + * This software may be freely distributed and used according to the terms + * of the GNU Lesser General Public License. + * + * Reference: + * - http://www.tecgraf.puc-rio.br/~mgattass/color/ColorIndex.html + */ + +#undef LOCAL_DEBUG + +#include "sherlock/sherlock.h" +#include "lib/math.h" +#include "images/color.h" + +u16 srgb_to_luv_tab1[256]; +u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE]; +u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE]; + +struct color_grid_node *srgb_to_luv_grid; +struct color_interpolation_node *color_interpolation_table; + +/* sRGB to XYZ */ +void +srgb_to_xyz_slow(double xyz[3], double srgb[3]) +{ + double a[3]; + for (uns i = 0; i < 3; i++) + if (srgb[i] > 0.04045) + a[i] = pow((srgb[i] + 0.055) * (1 / 1.055), 2.4); + else + a[i] = srgb[i] * (1 / 12.92); + xyz[0] = SRGB_XYZ_XR * a[0] + SRGB_XYZ_XG * a[1] + SRGB_XYZ_XB * a[2]; + xyz[1] = SRGB_XYZ_YR * a[0] + SRGB_XYZ_YG * a[1] + SRGB_XYZ_YB * a[2]; + xyz[2] = SRGB_XYZ_ZR * a[0] + SRGB_XYZ_ZG * a[1] + SRGB_XYZ_ZB * a[2]; +} + +/* XYZ to CIE-Luv */ +void +xyz_to_luv_slow(double luv[3], double xyz[3]) +{ + double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2]; + if (sum < 0.000001) + luv[0] = luv[1] = luv[2] = 0; + else + { + double var_u = 4 * xyz[0] / sum; + double var_v = 9 * xyz[1] / sum; + if (xyz[1] > 0.008856) + luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16; + else + luv[0] = (116 * 7.787) * xyz[1]; + luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z))); + luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z))); + /* intervals [0..100], [-134..220], [-140..122] */ + } +} + +void +srgb_to_luv_init(void) +{ + DBG("Initializing sRGB -> Luv table"); + for (uns i = 0; i < 256; i++) + { + double t = i / 255.; + if (t > 0.04045) + t = pow((t + 0.055) * (1 / 1.055), 2.4); + else + t = t * (1 / 12.92); + srgb_to_luv_tab1[i] = CLAMP(t * 0xfff + 0.5, 0, 0xfff); + } + for (uns i = 0; i < (9 << SRGB_TO_LUV_TAB2_SIZE); i++) + { + double t = i / (double)((9 << SRGB_TO_LUV_TAB2_SIZE) - 1); + if (t > 0.008856) + t = 1.16 * pow(t, 1 / 3.) - 0.16; + else + t = (1.16 * 7.787) * t; + srgb_to_luv_tab2[i] = + CLAMP(t * ((1 << SRGB_TO_LUV_TAB2_SCALE) - 1) + 0.5, + 0, (1 << SRGB_TO_LUV_TAB2_SCALE) - 1); + } + for (uns i = 0; i < (20 << SRGB_TO_LUV_TAB3_SIZE); i++) + { + srgb_to_luv_tab3[i] = i ? (13 << (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB3_SIZE)) / i : 0; + } +} + +void +srgb_to_luv_pixels(byte *dest, byte *src, uns count) +{ + while (count--) + { + srgb_to_luv_pixel(dest, src); + dest += 3; + src += 3; + } +} + +/* Returns volume of a given tetrahedron multiplied by 6 */ +static inline uns +tetrahedron_volume(uns *v1, uns *v2, uns *v3, uns *v4) +{ + int a[3], b[3], c[3]; + for (uns i = 0; i < 3; i++) + { + a[i] = v2[i] - v1[i]; + b[i] = v3[i] - v1[i]; + c[i] = v4[i] - v1[i]; + } + int result = + a[0] * (b[1] * c[2] - b[2] * c[1]) - + a[1] * (b[0] * c[2] - b[2] * c[0]) + + a[2] * (b[0] * c[1] - b[1] * c[0]); + return (result > 0) ? result : -result; +} + +static void +interpolate_tetrahedron(struct color_interpolation_node *n, uns *p, const uns *c) +{ + uns v[4][3]; + for (uns i = 0; i < 4; i++) + { + v[i][0] = (c[i] & 0001) ? (1 << COLOR_CONV_OFS) : 0; + v[i][1] = (c[i] & 0010) ? (1 << COLOR_CONV_OFS) : 0; + v[i][2] = (c[i] & 0100) ? (1 << COLOR_CONV_OFS) : 0; + n->ofs[i] = + ((c[i] & 0001) ? 1 : 0) + + ((c[i] & 0010) ? (1 << COLOR_CONV_SIZE) : 0) + + ((c[i] & 0100) ? (1 << (COLOR_CONV_SIZE * 2)) : 0); + } + uns vol = tetrahedron_volume(v[0], v[1], v[2], v[3]); + n->mul[0] = ((tetrahedron_volume(p, v[1], v[2], v[3]) << 8) + (vol >> 1)) / vol; + n->mul[1] = ((tetrahedron_volume(v[0], p, v[2], v[3]) << 8) + (vol >> 1)) / vol; + n->mul[2] = ((tetrahedron_volume(v[0], v[1], p, v[3]) << 8) + (vol >> 1)) / vol; + n->mul[3] = ((tetrahedron_volume(v[0], v[1], v[2], p) << 8) + (vol >> 1)) / vol; + uns j; + for (j = 0; j < 4; j++) + if (n->mul[j]) + break; + for (uns i = 0; i < 4; i++) + if (n->mul[i] == 0) + n->ofs[i] = n->ofs[j]; +} + +static void +interpolation_table_init(void) +{ + DBG("Initializing color interpolation table"); + struct color_interpolation_node *n = color_interpolation_table = + xmalloc(sizeof(struct color_interpolation_node) << (COLOR_CONV_OFS * 3)); + uns p[3]; + for (p[2] = 0; p[2] < (1 << COLOR_CONV_OFS); p[2]++) + for (p[1] = 0; p[1] < (1 << COLOR_CONV_OFS); p[1]++) + for (p[0] = 0; p[0] < (1 << COLOR_CONV_OFS); p[0]++) + { + uns index; + static const uns tetrahedrons[5][4] = { + {0000, 0001, 0010, 0100}, + {0110, 0111, 0100, 0010}, + {0101, 0100, 0111, 0001}, + {0011, 0010, 0001, 0111}, + {0111, 0001, 0010, 0100}}; + if (p[0] + p[1] + p[2] <= (1 << COLOR_CONV_OFS)) + index = 0; + else if ((1 << COLOR_CONV_OFS) + p[0] <= p[1] + p[2]) + index = 1; + else if ((1 << COLOR_CONV_OFS) + p[1] <= p[0] + p[2]) + index = 2; + else if ((1 << COLOR_CONV_OFS) + p[2] <= p[0] + p[1]) + index = 3; + else + index = 4; + interpolate_tetrahedron(n, p, tetrahedrons[index]); + n++; + } +} + +typedef void color_conv_func(double dest[3], double src[3]); + +static void +conv_grid_init(struct color_grid_node **grid, color_conv_func func) +{ + if (*grid) + return; + struct color_grid_node *g = *grid = xmalloc((sizeof(struct color_grid_node)) << (COLOR_CONV_SIZE * 3)); + double src[3], dest[3]; + for (uns k = 0; k < (1 << COLOR_CONV_SIZE); k++) + { + src[2] = k * (255 / (double)((1 << COLOR_CONV_SIZE) - 1)); + for (uns j = 0; j < (1 << COLOR_CONV_SIZE); j++) + { + src[1] = j * (255/ (double)((1 << COLOR_CONV_SIZE) - 1)); + for (uns i = 0; i < (1 << COLOR_CONV_SIZE); i++) + { + src[0] = i * (255 / (double)((1 << COLOR_CONV_SIZE) - 1)); + func(dest, src); + g->val[0] = CLAMP(dest[0] + 0.5, 0, 255); + g->val[1] = CLAMP(dest[1] + 0.5, 0, 255); + g->val[2] = CLAMP(dest[2] + 0.5, 0, 255); + g++; + } + } + } +} + +static void +srgb_to_luv_func(double dest[3], double src[3]) +{ + double srgb[3], xyz[3], luv[3]; + srgb[0] = src[0] / 255.; + srgb[1] = src[1] / 255.; + srgb[2] = src[2] / 255.; + srgb_to_xyz_slow(xyz, srgb); + xyz_to_luv_slow(luv, xyz); + dest[0] = luv[0] * 2.55; + dest[1] = luv[1] * (2.55 / 4) + 128; + dest[2] = luv[2] * (2.55 / 4) + 128; +} + +void +color_conv_init(void) +{ + interpolation_table_init(); + conv_grid_init(&srgb_to_luv_grid, srgb_to_luv_func); +} + +void +color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid) +{ + while (count--) + { + color_conv_pixel(dest, src, grid); + dest += 3; + src += 3; + } +} + +#ifdef TEST +#include + +static double +conv_error(u32 color, struct color_grid_node *grid, color_conv_func func) +{ + byte src[3], dest[3]; + src[0] = color & 255; + src[1] = (color >> 8) & 255; + src[2] = (color >> 16) & 255; + color_conv_pixel(dest, src, grid); + double src2[3], dest2[3]; + for (uns i = 0; i < 3; i++) + src2[i] = src[i]; + func(dest2, src2); + double err = 0; + for (uns i = 0; i < 3; i++) + err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]); + return err; +} + +typedef void test_fn(byte *dest, byte *src); + +static double +func_error(u32 color, test_fn test, color_conv_func func) +{ + byte src[3], dest[3]; + src[0] = color & 255; + src[1] = (color >> 8) & 255; + src[2] = (color >> 16) & 255; + test(dest, src); + double src2[3], dest2[3]; + for (uns i = 0; i < 3; i++) + src2[i] = src[i]; + func(dest2, src2); + double err = 0; + for (uns i = 0; i < 3; i++) + err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]); + return err; +} + +static void +test_grid(byte *name, struct color_grid_node *grid, color_conv_func func) +{ + double max_err = 0, sum_err = 0; + uns count = 100000; + for (uns i = 0; i < count; i++) + { + double err = conv_error(random_max(0x1000000), grid, func); + max_err = MAX(err, max_err); + sum_err += err; + } + DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count); + if (max_err > 12) + die("Too large error in %s conversion", name); +} + +static void +test_func(byte *name, test_fn test, color_conv_func func) +{ + double max_err = 0, sum_err = 0; + uns count = 100000; + for (uns i = 0; i < count; i++) + { + double err = func_error(random_max(0x1000000), test, func); + max_err = MAX(err, max_err); + sum_err += err; + } + DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count); + if (max_err > 12) + die("Too large error in %s conversion", name); +} + +int +main(void) +{ + srgb_to_luv_init(); + test_func("func sRGB -> Luv", srgb_to_luv_pixel, srgb_to_luv_func); + color_conv_init(); + test_grid("grid sRGB -> Luv", srgb_to_luv_grid, srgb_to_luv_func); +#ifdef LOCAL_DEBUG +#define CNT 1000000 + byte *a = xmalloc(3 * CNT), *b = xmalloc(3 * CNT); + init_timer(); + for (uns i = 0; i < 20; i++) + memcpy(b, a, CNT * 3); + DBG("memcpy time=%d", (uns)get_timer()); + for (uns i = 0; i < 3 * CNT; i++) + a[i] = random_max(256); + init_timer(); + for (uns i = 0; i < 20; i++) + srgb_to_luv_pixels(b, a, CNT); + DBG("direct time=%d", (uns)get_timer()); + init_timer(); + for (uns i = 0; i < 20; i++) + color_conv_pixels(b, a, CNT, srgb_to_luv_grid); + DBG("grid time=%d", (uns)get_timer()); +#endif + return 0; +} +#endif + diff --git a/images/color.h b/images/color.h new file mode 100644 index 00000000..cb64a9c2 --- /dev/null +++ b/images/color.h @@ -0,0 +1,135 @@ +/* + * Image Library -- Color Spaces + * + * (c) 2006 Pavel Charvat + * + * This software may be freely distributed and used according to the terms + * of the GNU Lesser General Public License. + * + * FIXME: + * - fix theoretical problems with rounding errors in srgb_to_luv_pixel() + * - SIMD should help to speed up conversion of large arrays + * - maybe try to generate a long switch in color_conv_pixel() + * with optimized entries instead of access to interpolation table + */ + +#ifndef _IMAGES_COLOR_H +#define _IMAGES_COLOR_H + +/* Exact slow conversion routines */ +void srgb_to_xyz_slow(double dest[3], double src[3]); +void xyz_to_luv_slow(double dest[3], double src[3]); + +/* Reference white */ +#define REF_WHITE_X 0.96422 +#define REF_WHITE_Y 1. +#define REF_WHITE_Z 0.82521 + +/* sRGB -> XYZ matrix */ +#define SRGB_XYZ_XR 0.412424 +#define SRGB_XYZ_XG 0.357579 +#define SRGB_XYZ_XB 0.180464 +#define SRGB_XYZ_YR 0.212656 +#define SRGB_XYZ_YG 0.715158 +#define SRGB_XYZ_YB 0.072186 +#define SRGB_XYZ_ZR 0.019332 +#define SRGB_XYZ_ZG 0.119193 +#define SRGB_XYZ_ZB 0.950444 + + +/*********************** OPTIMIZED CONVERSION ROUTINES **********************/ + +/* sRGB -> Luv parameters */ +#define SRGB_TO_LUV_TAB2_SIZE 9 +#define SRGB_TO_LUV_TAB2_SCALE 11 +#define SRGB_TO_LUV_TAB3_SIZE 8 +#define SRGB_TO_LUV_TAB3_SCALE (39 - SRGB_TO_LUV_TAB2_SCALE - SRGB_TO_LUV_TAB3_SIZE) + +extern u16 srgb_to_luv_tab1[256]; +extern u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE]; +extern u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE]; + +void srgb_to_luv_init(void); +void srgb_to_luv_pixels(byte *dest, byte *src, uns count); + +static inline void +srgb_to_luv_pixel(byte *dest, byte *src) +{ + uns r = srgb_to_luv_tab1[src[0]]; + uns g = srgb_to_luv_tab1[src[1]]; + uns b = srgb_to_luv_tab1[src[2]]; + uns x = + (uns)(4 * SRGB_XYZ_XR * 0xffff) * r + + (uns)(4 * SRGB_XYZ_XG * 0xffff) * g + + (uns)(4 * SRGB_XYZ_XB * 0xffff) * b; + uns y = + (uns)(9 * SRGB_XYZ_YR * 0xffff) * r + + (uns)(9 * SRGB_XYZ_YG * 0xffff) * g + + (uns)(9 * SRGB_XYZ_YB * 0xffff) * b; + uns l = srgb_to_luv_tab2[y >> (28 - SRGB_TO_LUV_TAB2_SIZE)]; + dest[0] = l >> (SRGB_TO_LUV_TAB2_SCALE - 8); + uns sum = + (uns)((SRGB_XYZ_XR + 15 * SRGB_XYZ_YR + 3 * SRGB_XYZ_ZR) * 0x7fff) * r + + (uns)((SRGB_XYZ_XG + 15 * SRGB_XYZ_YG + 3 * SRGB_XYZ_ZG) * 0x7fff) * g + + (uns)((SRGB_XYZ_XB + 15 * SRGB_XYZ_YB + 3 * SRGB_XYZ_ZB) * 0x7fff) * b; + uns s = srgb_to_luv_tab3[sum >> (27 - SRGB_TO_LUV_TAB3_SIZE)]; + int xs = ((u64)x * s) >> 32; + int ys = ((u64)y * s) >> 32; + int xw = ((4 * 13) << (SRGB_TO_LUV_TAB3_SCALE - 4)) * + REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z); + int yw = ((9 * 13) << (SRGB_TO_LUV_TAB3_SCALE - 4)) * + REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z); + int u = (int)(l) * (xs - xw); + int v = (int)(l) * (ys - yw); + dest[1] = 128 + (u >> (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB2_SCALE - 10)); + dest[2] = 128 + (v >> (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB2_SCALE - 10)); +} + + +/****************** GENERAL INTERPOLATION IN 3D GRID ********************/ + +#define COLOR_CONV_SIZE 5 /* 128K conversion grid size */ +#define COLOR_CONV_OFS 3 /* 8K interpolation table size */ + +struct color_grid_node { + byte val[4]; +}; + +struct color_interpolation_node { + u16 ofs[4]; + u16 mul[4]; +}; + +extern struct color_grid_node *srgb_to_luv_grid; +extern struct color_interpolation_node *color_interpolation_table; + +void color_conv_init(void); +void color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid); + +#define COLOR_CONV_SCALE_CONST (((((1 << COLOR_CONV_SIZE) - 1) << 16) + (1 << (16 - COLOR_CONV_OFS))) / 255) + +static inline void +color_conv_pixel(byte *dest, byte *src, struct color_grid_node *grid) +{ + uns s0 = src[0] * COLOR_CONV_SCALE_CONST; + uns s1 = src[1] * COLOR_CONV_SCALE_CONST; + uns s2 = src[2] * COLOR_CONV_SCALE_CONST; + struct color_grid_node *g0, *g1, *g2, *g3, *g = grid + + ((s0 >> 16) + ((s1 >> 16) << COLOR_CONV_SIZE) + ((s2 >> 16) << (2 * COLOR_CONV_SIZE))); + struct color_interpolation_node *n = color_interpolation_table + + (((s0 & (0x10000 - (0x10000 >> COLOR_CONV_OFS))) >> (16 - COLOR_CONV_OFS)) + + ((s1 & (0x10000 - (0x10000 >> COLOR_CONV_OFS))) >> (16 - 2 * COLOR_CONV_OFS)) + + ((s2 & (0x10000 - (0x10000 >> COLOR_CONV_OFS))) >> (16 - 3 * COLOR_CONV_OFS))); + g0 = g + n->ofs[0]; + g1 = g + n->ofs[1]; + g2 = g + n->ofs[2]; + g3 = g + n->ofs[3]; + dest[0] = (g0->val[0] * n->mul[0] + g1->val[0] * n->mul[1] + + g2->val[0] * n->mul[2] + g3->val[0] * n->mul[3] + 128) >> 8; + dest[1] = (g0->val[1] * n->mul[0] + g1->val[1] * n->mul[1] + + g2->val[1] * n->mul[2] + g3->val[1] * n->mul[3] + 128) >> 8; + dest[2] = (g0->val[2] * n->mul[0] + g1->val[2] * n->mul[1] + + g2->val[2] * n->mul[2] + g3->val[2] * n->mul[3] + 128) >> 8; +} + +#endif diff --git a/images/color.t b/images/color.t new file mode 100644 index 00000000..42596089 --- /dev/null +++ b/images/color.t @@ -0,0 +1,4 @@ +# Tests for color conversion module + +Run: obj/images/color-t +