--- ray/src/common/tonemap.c	2005/01/07 20:33:02	3.16
+++ ray/src/common/tonemap.c	2007/11/27 02:42:54	3.30
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: tonemap.c,v 3.16 2005/01/07 20:33:02 greg Exp $";
+static const char	RCSid[] = "$Id: tonemap.c,v 3.30 2007/11/27 02:42:54 greg Exp $";
 #endif
 /*
  * Tone mapping functions.
@@ -22,10 +22,12 @@ static const char	RCSid[] = "$Id: tonemap.c,v 3.16 200
 struct tmPackage	*tmPkg[TM_MAXPKG];
 int	tmNumPkgs = 0;			/* number of registered packages */
 
-int	tmLastError;			/* last error incurred by library */
-char	*tmLastFunction;		/* error-generating function name */
+					/* luminance->brightness lookup */
+static TMbright		*tmFloat2BrtLUT = NULL;
 
+#define tmCvLumLUfp(pf)	tmFloat2BrtLUT[*(int32 *)(pf) >> 15]
 
+
 TMstruct *
 tmInit(					/* initialize new tone mapping */
 int	flags,
@@ -51,7 +53,7 @@ double	gamval
 		tmnew->clf[GRN] = rgb2xyzmat[1][1];
 		tmnew->clf[BLU] = rgb2xyzmat[1][2];
 	} else {
-		comprgb2xyzWBmat(cmat, tmnew->monpri=monpri);
+		comprgb2xyzmat(cmat, tmnew->monpri=monpri);
 		tmnew->clf[RED] = cmat[1][0];
 		tmnew->clf[GRN] = cmat[1][1];
 		tmnew->clf[BLU] = cmat[1][2];
@@ -71,6 +73,7 @@ double	gamval
 			tmnew->inpsf = WHTEFFICACY;
 	tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
 	tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
+	tmnew->inpdat = NULL;
 	tmnew->hbrmin = 10; tmnew->hbrmax = -10;
 	tmnew->histo = NULL;
 	tmnew->mbrmin = 10; tmnew->mbrmax = -10;
@@ -78,6 +81,8 @@ double	gamval
 						/* zero private data */
 	for (i = TM_MAXPKG; i--; )
 		tmnew->pd[i] = NULL;
+	tmnew->lastError = TM_E_OK;
+	tmnew->lastFunc = "NoErr";
 						/* return new TMstruct */
 	return(tmnew);
 }
@@ -87,10 +92,11 @@ int
 tmSetSpace(			/* set input color space for conversions */
 TMstruct	*tms,
 RGBPRIMP	pri,
-double	sf
+double	sf,
+MEM_PTR	dat
 )
 {
-	static char	funcName[] = "tmSetSpace";
+	static const char funcName[] = "tmSetSpace";
 	int	i, j;
 						/* error check */
 	if (tms == NULL)
@@ -98,10 +104,11 @@ double	sf
 	if (sf <= 1e-12)
 		returnErr(TM_E_ILLEGAL);
 						/* check if no change */
-	if (pri == tms->inppri && FEQ(sf, tms->inpsf))
+	if (pri == tms->inppri && FEQ(sf, tms->inpsf) && dat == tms->inpdat)
 		returnOK;
 	tms->inppri = pri;			/* let's set it */
 	tms->inpsf = sf;
+	tms->inpdat = dat;
 
 	if (tms->flags & TM_F_BW) {		/* color doesn't matter */
 		tms->monpri = tms->inppri;		/* eliminate xform */
@@ -109,7 +116,7 @@ double	sf
 			tms->clf[CIEX] = tms->clf[CIEZ] = 0.;
 			tms->clf[CIEY] = 1.;
 		} else {
-			comprgb2xyzWBmat(tms->cmat, tms->monpri);
+			comprgb2xyzmat(tms->cmat, tms->monpri);
 			tms->clf[RED] = tms->cmat[1][0];
 			tms->clf[GRN] = tms->cmat[1][1];
 			tms->clf[BLU] = tms->cmat[1][2];
@@ -147,7 +154,7 @@ double	sf
 
 
 void
-tmClearHisto(			/* clear current histogram */
+tmClearHisto(				/* clear current histogram */
 TMstruct	*tms
 )
 {
@@ -158,7 +165,93 @@ TMstruct	*tms
 }
 
 
+TMbright
+tmCvLuminance(				/* convert a single luminance */
+double	lum
+)
+{
+	double	d;
+
+#ifdef isfinite
+	if (!isfinite(lum) || lum <= TM_NOLUM)
+#else
+	if (lum <= TM_NOLUM)
+#endif
+		return(TM_NOBRT);
+	d = TM_BRTSCALE*log(lum);
+	if (d > 0.)
+		return((TMbright)(d+.5));
+	return((TMbright)(d-.5));
+}
+
+
 int
+tmCvLums(				/* convert luminances using lookup */
+TMbright	*ls,
+float		*scan,
+int		len
+)
+{
+	if (tmFloat2BrtLUT == NULL) {	/* initialize lookup table */
+		int32	i;
+		tmFloat2BrtLUT = (TMbright *)malloc(sizeof(TMbright)*0x10000);
+		if (tmFloat2BrtLUT == NULL)
+			return(TM_E_NOMEM);
+		for (i = 0; i < 0x10000; i++) {
+			int32	l = (i<<1 | 1) << 14;
+#ifndef isfinite
+			if ((l & 0x7f800000) == 0x7f800000)
+				tmFloat2BrtLUT[i] = TM_NOBRT;
+			else
+#endif
+			tmFloat2BrtLUT[i] = tmCvLuminance(*(float *)&l);
+		}
+	}
+	if (len <= 0)
+		return(TM_E_OK);
+	if ((ls == NULL) | (scan == NULL))
+		return(TM_E_ILLEGAL);
+	while (len--) {
+		if (*scan <= TM_NOLUM) {
+			*ls++ = TM_NOBRT;
+			++scan;
+			continue;
+		}
+		*ls++ = tmCvLumLUfp(scan++);
+	}
+	return(TM_E_OK);
+}
+
+
+int
+tmCvGrays(				/* convert float gray values */
+TMstruct	*tms,
+TMbright	*ls,
+float		*scan,
+int		len
+)
+{
+	static const char funcName[] = "tmCvGrays";
+	int	i;
+
+	if (tms == NULL)
+		returnErr(TM_E_TMINVAL);
+	if ((ls == NULL) | (scan == NULL) | (len < 0))
+		returnErr(TM_E_ILLEGAL);
+	if (tmFloat2BrtLUT == NULL)			/* initialize */
+		tmCvLums(NULL, NULL, 0);
+	for (i = len; i--; ) {
+		float	lum = tms->inpsf * scan[i];
+		if (lum <= TM_NOLUM)
+			ls[i] = TM_NOBRT;
+		else
+			ls[i] = tmCvLumLUfp(&lum);
+	}
+	returnOK;
+}
+
+
+int
 tmCvColors(				/* convert float colors */
 TMstruct	*tms,
 TMbright	*ls,
@@ -167,17 +260,24 @@ COLOR	*scan,
 int	len
 )
 {
-	static char	funcName[] = "tmCvColors";
-	static COLOR	csmall = {.5*MINLUM, .5*MINLUM, .5*MINLUM};
+	static const char funcName[] = "tmCvColors";
+	static BYTE	gamtab[1024];
+	static double	curgam = .0;
 	COLOR	cmon;
-	double	lum, slum;
-	double	d;
+	float	lum, slum, d;
 	int	i;
 
 	if (tms == NULL)
 		returnErr(TM_E_TMINVAL);
 	if ((ls == NULL) | (scan == NULL) | (len < 0))
 		returnErr(TM_E_ILLEGAL);
+	if (tmFloat2BrtLUT == NULL)			/* initialize */
+		tmCvLums(NULL, NULL, 0);
+	if (cs != TM_NOCHROM && fabs(tms->mongam - curgam) > .02) {
+		curgam = tms->mongam;			/* (re)build table */
+		for (i = 1024; i--; )
+			gamtab[i] = (int)(256.*pow((i+.5)/1024., 1./curgam));
+	}
 	for (i = len; i--; ) {
 		if (tmNeedMatrix(tms)) {		/* get monitor RGB */
 			colortrans(cmon, tms->cmat, scan[i]);
@@ -186,36 +286,38 @@ int	len
 			cmon[GRN] = tms->inpsf*scan[i][GRN];
 			cmon[BLU] = tms->inpsf*scan[i][BLU];
 		}
+#ifdef isfinite
+		if (!isfinite(cmon[RED]) || cmon[RED] < .0f) cmon[RED] = .0f;
+		if (!isfinite(cmon[GRN]) || cmon[GRN] < .0f) cmon[GRN] = .0f;
+		if (!isfinite(cmon[BLU]) || cmon[BLU] < .0f) cmon[BLU] = .0f;
+#else
+		if (cmon[RED] < .0f) cmon[RED] = .0f;
+		if (cmon[GRN] < .0f) cmon[GRN] = .0f;
+		if (cmon[BLU] < .0f) cmon[BLU] = .0f;
+#endif
 							/* world luminance */
 		lum =	tms->clf[RED]*cmon[RED] +
 			tms->clf[GRN]*cmon[GRN] +
 			tms->clf[BLU]*cmon[BLU] ;
-							/* check range */
-		if (clipgamut(cmon, lum, CGAMUT_LOWER, csmall, cwhite))
-			lum =	tms->clf[RED]*cmon[RED] +
-				tms->clf[GRN]*cmon[GRN] +
-				tms->clf[BLU]*cmon[BLU] ;
-		if (lum < MINLUM) {
-			ls[i] = MINBRT-1;		/* bogus value */
-			lum = MINLUM;
-		} else {
-			d = TM_BRTSCALE*log(lum);	/* encode it */
-			ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
-		}
+		if (lum <= TM_NOLUM) {			/* convert brightness */
+			lum = cmon[RED] = cmon[GRN] = cmon[BLU] = TM_NOLUM;
+			ls[i] = TM_NOBRT;
+		} else
+			ls[i] = tmCvLumLUfp(&lum);
 		if (cs == TM_NOCHROM)			/* no color? */
 			continue;
 		if (tms->flags & TM_F_MESOPIC && lum < LMESUPPER) {
 			slum = scotlum(cmon);		/* mesopic adj. */
-			if (lum < LMESLOWER)
+			if (lum < LMESLOWER) {
 				cmon[RED] = cmon[GRN] = cmon[BLU] = slum;
-			else {
+			} else {
 				d = (lum - LMESLOWER)/(LMESUPPER - LMESLOWER);
 				if (tms->flags & TM_F_BW)
 					cmon[RED] = cmon[GRN] =
 							cmon[BLU] = d*lum;
 				else
 					scalecolor(cmon, d);
-				d = (1.-d)*slum;
+				d = (1.f-d)*slum;
 				cmon[RED] += d;
 				cmon[GRN] += d;
 				cmon[BLU] += d;
@@ -224,47 +326,17 @@ int	len
 			cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
 		}
 		d = tms->clf[RED]*cmon[RED]/lum;
-		cs[3*i  ] = d>=.999 ? 255 :
-				(int)(256.*pow(d, 1./tms->mongam));
+		cs[3*i  ] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
 		d = tms->clf[GRN]*cmon[GRN]/lum;
-		cs[3*i+1] = d>=.999 ? 255 :
-				(int)(256.*pow(d, 1./tms->mongam));
+		cs[3*i+1] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
 		d = tms->clf[BLU]*cmon[BLU]/lum;
-		cs[3*i+2] = d>=.999 ? 255 :
-				(int)(256.*pow(d, 1./tms->mongam));
+		cs[3*i+2] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
 	}
 	returnOK;
 }
 
 
 int
-tmCvGrays(				/* convert float gray values */
-TMstruct	*tms,
-TMbright	*ls,
-float	*scan,
-int	len
-)
-{
-	static char	funcName[] = "tmCvGrays";
-	double	d;
-	int	i;
-
-	if (tms == NULL)
-		returnErr(TM_E_TMINVAL);
-	if ((ls == NULL) | (scan == NULL) | (len < 0))
-		returnErr(TM_E_ILLEGAL);
-	for (i = len; i--; )
-		if (scan[i] <= TM_NOLUM) {
-			ls[i] = TM_NOBRT;		/* bogus value */
-		} else {
-			d = TM_BRTSCALE*log(scan[i]);	/* encode it */
-			ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
-		}
-	returnOK;
-}
-
-
-int
 tmAddHisto(				/* add values to histogram */
 TMstruct	*tms,
 TMbright	*ls,
@@ -272,7 +344,7 @@ int	len,
 int	wt
 )
 {
-	static char	funcName[] = "tmAddHisto";
+	static const char funcName[] = "tmAddHisto";
 	int	oldorig=0, oldlen, horig, hlen;
 	int	i, j;
 
@@ -375,7 +447,7 @@ double	expmult,
 double	gamval
 )
 {
-	static char	funcName[] = "tmFixedMapping";
+	static const char funcName[] = "tmFixedMapping";
 	double		d;
 	int	i;
 	
@@ -386,10 +458,12 @@ double	gamval
 	if (gamval < MINGAM)
 		gamval = tms->mongam;
 	d = log(expmult/tms->inpsf);
-	for (i = tms->mbrmax-tms->mbrmin+1; i--; )
-		tms->lumap[i] = 256. * exp(
+	for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
+		double	val = 256. * exp(
 			( d + (tms->mbrmin+i)*(1./TM_BRTSCALE) )
-			/ gamval );
+			/ gamval);
+		tms->lumap[i] = val >= (double)0xffff ? 0xffff : (int)val;
+	}
 	returnOK;
 }
 
@@ -402,7 +476,7 @@ double	Lddyn,
 double	Ldmax
 )
 {
-	static char	funcName[] = "tmComputeMapping";
+	static const char funcName[] = "tmComputeMapping";
 	int	*histo;
 	float	*cumf;
 	int	brt0, histlen, threshold, ceiling, trimmings;
@@ -430,17 +504,15 @@ double	Ldmax
 	j = brt0 + histlen*HISTEP;
 	for (i = histlen; i--; ) {
 		histot += tms->histo[i];
-		sum += (j -= HISTEP) * tms->histo[i];
+		sum += (double)(j -= HISTEP) * tms->histo[i];
 	}
 	threshold = histot*0.005 + .5;
-	if (threshold < 4)
+	if (!histot)
 		returnErr(TM_E_TMFAIL);
 	Lwavg = tmLuminance( (double)sum / histot );
-					/* allocate space for mapping */
-	if (!tmNewMap(tms))
-		returnErr(TM_E_NOMEM);
 					/* use linear tone mapping? */
-	if (tms->flags & TM_F_LINEAR)
+	if (tms->flags & TM_F_LINEAR || threshold < 4 ||
+			tms->hbrmax - tms->hbrmin < TM_BRTSCALE*logLddyn)
 		goto linearmap;
 					/* clamp histogram */
 	histo = (int *)malloc(histlen*sizeof(int));
@@ -481,6 +553,9 @@ double	Ldmax
 			goto linearmap;
 		}
 	} while (trimmings > threshold);
+					/* allocate space for mapping */
+	if (!tmNewMap(tms))
+		returnErr(TM_E_NOMEM);
 					/* assign tone-mapping */
 	for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
 		j = d = (double)i/(tms->mbrmax-tms->mbrmin)*histlen;
@@ -510,7 +585,7 @@ BYTE	*cs,
 int	len
 )
 {
-	static char	funcName[] = "tmMapPixels";
+	static const char funcName[] = "tmMapPixels";
 	int32	li, pv;
 
 	if (tms == NULL || tms->lumap == NULL)
@@ -622,15 +697,17 @@ tmMkMesofact()				/* build mesopic lookup factor table
 
 int
 tmErrorReturn(				/* error return (with message) */
-char	*func,
+const char	*func,
 TMstruct	*tms,
 int	err
 )
 {
-	tmLastFunction = func;
-	tmLastError = err;
-	if (tms != NULL && tms->flags & TM_F_NOSTDERR)
-		return(err);
+	if (tms != NULL) {
+		tms->lastFunc = func;
+		tms->lastError = err;
+		if (tms->flags & TM_F_NOSTDERR)
+			return(err);
+	}
 	fputs(func, stderr);
 	fputs(": ", stderr);
 	fputs(tmErrorMessage[err], stderr);