faster and more accurate sRGB <-> XYZ

just use one table, since all colour channels are the same now
have more points in the float -> int direction, fewer in int -> float
faster out of range detection

ChangeLog

@@ -7,6 +7,7 @@
   im_dECMC_fromLab(), im_dE00_from_Lab()
   as classes
 - added vips_colour_convert(), replaces all derived conversions
+- faster and more accurate sRGB <-> XYZ conversion
 - dzsave can write zoomify and google maps layout as well
 - tilecache supports threaded access
 - openslide2vips gets underlying tile size from openslide

TODO

@@ -1,11 +1,10 @@
-- see vips_col_sRGB2XYZ()
+- add something to dzsave to control pyramid depth ... can then use it for
-
+  straight tiling
-	t_r2Yr[] is a 1501 element array, but we just index with a uchar
-
-	trim this down!
 
 - add mono as a colourspace? also rad?
 
+- rename INTERPRETATION_UCS to _CMC
+
 - something to test if an image is in a supported colourspace?
 
   at the moment we only look at interpretation, but for things like labq 

libvips/colour/LabQ2sRGB.c

@@ -10,6 +10,8 @@
  * 21/9/12
  * 	- redone as a class
  * 	- sRGB only, support for other RGBs is now via lcms
+ * 1/11/12
+ * 	- faster and more accurate sRGB <-> XYZ conversion
  */
 
 /*
@@ -51,6 +53,8 @@
 
 #include "colour.h"
 
+#define TABLE_SIZE (20000)
+
 typedef VipsColourCode VipsLabQ2sRGB;
 typedef VipsColourCodeClass VipsLabQ2sRGBClass;
 
@@ -88,16 +92,11 @@ struct im_col_display {
  */
 struct im_col_tab_disp {
 	/*< private >*/
-	float	t_Yr2r[1501];		/* Conversion of Yr to r */
+	float	t_Y2v[TABLE_SIZE];	/* Conversion of Y to v */
-	float	t_Yg2g[1501];		/* Conversion of Yg to g */
+	float	t_v2Y[256];		/* Conversion of v to Y */
-	float	t_Yb2b[1501];		/* Conversion of Yb to b */
-	float	t_r2Yr[1501];		/* Conversion of r to Yr */
-	float	t_g2Yg[1501];		/* Conversion of g to Yg */
-	float	t_b2Yb[1501];		/* Conversion of b to Yb */
 	float	mat_XYZ2lum[3][3];	/* XYZ to Yr, Yg, Yb matrix */
 	float	mat_lum2XYZ[3][3];	/* Yr, Yg, Yb to XYZ matrix */
-	float rstep, gstep, bstep;
+	float rstep;			/* Scale Y by this to fit TABLE_SIZE */
-	float ristep, gistep, bistep;
 };
 
 /* Do our own indexing of the arrays below to make sure we get efficient mults.
@@ -140,7 +139,6 @@ calcul_tables( void *client )
 
 	int i, j;
 	float a, ga_i, ga, c, f, yo, p;
-	float maxr, maxg, maxb;
 	double **temp;
 
 	c = (d->d_B - 100.0) / 500.0;
@@ -153,54 +151,13 @@ calcul_tables( void *client )
 	p = d->d_P / 100.0;
 	f = d->d_Vrwr / p;
 
-	maxr = (float) d->d_Vrwr;
+	table->rstep = a / (TABLE_SIZE - 1);
-	table->ristep = maxr / 1500.0;
-	table->rstep = a / 1500.0;
 
-	for( i = 0; i < 1501; i++ )
+	for( i = 0; i < TABLE_SIZE; i++ )
-		table->t_Yr2r[i] = f * (pow( i * table->rstep / a, ga_i ) - c);
+		table->t_Y2v[i] = f * (pow( i * table->rstep / a, ga_i ) - c);
 
-	for( i = 0; i < 1501; i++ )
+	for( i = 0; i < 256; i++ )
-		table->t_r2Yr[i] = yo + 
+		table->t_v2Y[i] = yo + a * pow( i / f + c, ga );
-			a * pow( i * table->ristep / f + c, ga );
-
-	/**** Green ****/
-	yo = d->d_Y0G;
-	a = d->d_YCG - yo;
-	ga = d->d_gammaG;
-	ga_i = 1.0 / ga;
-	p = d->d_P / 100.0;
-	f = d->d_Vrwg / p;
-
-	maxg = (float)d->d_Vrwg;
-	table->gistep = maxg / 1500.0;
-	table->gstep = a / 1500.0;
-
-	for( i = 0; i < 1501; i++ )
-		table->t_Yg2g[i] = f * (pow( i * table->gstep / a, ga_i ) - c);
-
-	for( i = 0; i < 1501; i++ )
-		table->t_g2Yg[i] = yo + 
-			a * pow( i * table->gistep / f + c, ga );
-
-	/**** Blue ****/
-	yo = d->d_Y0B;
-	a = d->d_YCB - yo;
-	ga = d->d_gammaB;
-	ga_i = 1.0 / ga;
-	p = d->d_P / 100.0;
-	f = d->d_Vrwb / p;
-
-	maxb = (float)d->d_Vrwb;
-	table->bistep = maxb / 1500.0;
-	table->bstep = a / 1500.0;
-
-	for( i = 0; i < 1501; i++ )
-		table->t_Yb2b[i] = f * (pow( i * table->bstep / a, ga_i ) - c);
-
-	for( i = 0; i < 1501; i++ )
-		table->t_b2Yb[i] = yo + 
-			a * pow( i * table->bistep / f + c, ga );
 
 	if( !(temp = im_dmat_alloc( 0, 2, 0, 2 )) )
 		return( NULL );
@@ -256,18 +213,14 @@ vips_col_sRGB2XYZ( int r, int g, int b, float *X, float *Y, float *Z )
 	float Yr, Yg, Yb;
 	int i;
 
-  	r = VIPS_CLIP( 0, r, 255 );
+  	i = VIPS_CLIP( 0, r, 255 );
-  	g = VIPS_CLIP( 0, g, 255 );
+	Yr = table->t_v2Y[i];
-  	b = VIPS_CLIP( 0, b, 255 );
 
-	i = r / table->ristep;
+  	i = VIPS_CLIP( 0, g, 255 );
-	Yr = table->t_r2Yr[i];
+	Yg = table->t_v2Y[i];
 
-	i = g / table->gistep;
+  	i = VIPS_CLIP( 0, b, 255 );
-	Yg = table->t_g2Yg[i];
+	Yb = table->t_v2Y[i];
-
-	i = b / table->bistep;
-	Yb = table->t_b2Yb[i];
 
 	*X = mat[0] * Yr + mat[1] * Yg + mat[2] * Yb;
 	*Y = mat[3] * Yr + mat[4] * Yg + mat[5] * Yb;
@@ -300,35 +253,33 @@ vips_col_XYZ2sRGB( float X, float Y, float Z,
 	Yg = mat[3] * X + mat[4] * Y + mat[5] * Z;
 	Yb = mat[6] * X + mat[7] * Y + mat[8] * Z;
 
-	/* Any negatives? If yes, set the out-of-range flag and bump up.
+	/* Clip range, set the out-of-range flag.
 	 */
-	if( Yr < d->d_Y0R ) { 
+#define CLIP( L, V, H ) { \
-		or = 1; 
+	if( (V) < (L) ) { \
-		Yr = d->d_Y0R; 
+		(V) = (L); \
-	}
+		or = 1; \
-	if( Yg < d->d_Y0G ) { 
+	} \
-		or = 1; 
+	if( (V) > (H) ) { \
-		Yg = d->d_Y0G; 
+		(V) = (H); \
-	}
+		or = 1; \
-	if( Yb < d->d_Y0B ) { 
+	} \
-		or = 1; 
+}
-		Yb = d->d_Y0B; 
-	}
 
 	/* Work out colour value (0-Vrw) to feed the tube to get that
 	 * luminosity. 
 	 */
 	Yint = (Yr - d->d_Y0R) / table->rstep;
-	Yint = VIPS_CLIP( 0, Yint, 1500 ); 
+	CLIP( 0, Yint, TABLE_SIZE - 1);
-	r = IM_RINT( table->t_Yr2r[Yint] );
+	r = VIPS_RINT( table->t_Y2v[Yint] );
 
-	Yint = (Yg - d->d_Y0G) / table->gstep;
+	Yint = (Yg - d->d_Y0G) / table->rstep;
-	Yint = VIPS_CLIP( 0, Yint, 1500 ); 
+	CLIP( 0, Yint, TABLE_SIZE - 1);
-	g = IM_RINT( table->t_Yg2g[Yint] );
+	g = VIPS_RINT( table->t_Y2v[Yint] );
 
-	Yint = (Yb - d->d_Y0B) / table->bstep;
+	Yint = (Yb - d->d_Y0B) / table->rstep;
-	Yint = VIPS_CLIP( 0, Yint, 1500 ); 
+	CLIP( 0, Yint, TABLE_SIZE - 1);
-	b = IM_RINT( table->t_Yb2b[Yint] );
+	b = VIPS_RINT( table->t_Y2v[Yint] );
 
 	*r_ret = r;
 	*g_ret = g;
@@ -468,9 +419,9 @@ vips_LabQ2sRGB_init( VipsLabQ2sRGB *LabQ2sRGB )
  * @in: input image
  * @out: output image
  *
- * Unpack a LabQ (#IM_CODING_LABQ) image to a three-band short image.
+ * Unpack a LabQ (#VIPS_CODING_LABQ) image to a three-band short image.
  *
- * See also: im_LabS2LabQ(), im_LabQ2sRGB(), im_rad2float().
+ * See also: vips_LabS2LabQ(), vips_LabQ2sRGB(), vips_rad2float().
  *
  * Returns: 0 on success, -1 on error.
  */