--- ray/src/common/tonemap.c	1997/11/17 14:02:12	3.6
+++ ray/src/common/tonemap.c	2009/02/09 20:23:51	3.33
@@ -1,14 +1,16 @@
-/* Copyright (c) 1997 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: tonemap.c,v 3.33 2009/02/09 20:23:51 greg Exp $";
 #endif
-
 /*
  * Tone mapping functions.
  * See tonemap.h for detailed function descriptions.
+ * Added von Kries white-balance calculations 10/01 (GW).
+ *
+ * Externals declared in tonemap.h
  */
 
+#include "copyright.h"
+
 #include	<stdio.h>
 #include	<math.h>
 #include	"tmprivat.h"
@@ -16,49 +18,34 @@ static char SCCSid[] = "$SunId$ LBL";
 
 #define	exp10(x)	exp(M_LN10*(x))
 
-struct tmStruct	*tmTop = NULL;		/* current tone mapping stack */
-
 					/* our list of conversion packages */
 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]
 
-int
-tmErrorReturn(func, err)		/* error return (with message) */
-char	*func;
-int	err;
-{
-	tmLastFunction = func;
-	tmLastError = err;
-	if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
-		return(err);
-	fputs(func, stderr);
-	fputs(": ", stderr);
-	fputs(tmErrorMessage[err], stderr);
-	fputs("!\n", stderr);
-	return(err);
-}
 
-
-struct tmStruct *
-tmInit(flags, monpri, gamval)		/* initialize new tone mapping */
-int	flags;
-RGBPRIMP	monpri;
-double	gamval;
+TMstruct *
+tmInit(					/* initialize new tone mapping */
+int	flags,
+RGBPRIMP	monpri,
+double	gamval
+)
 {
-	static char	funcName[] = "tmInit";
 	COLORMAT	cmat;
-	register struct tmStruct	*tmnew;
-	register int	i;
+	TMstruct	*tmnew;
+	int	i;
 						/* allocate structure */
-	tmnew = (struct tmStruct *)malloc(sizeof(struct tmStruct));
+	tmnew = (TMstruct *)malloc(sizeof(TMstruct));
 	if (tmnew == NULL)
 		return(NULL);
 
 	tmnew->flags = flags & ~TM_F_UNIMPL;
+	if (tmnew->flags & TM_F_BW)
+		tmnew->flags &= ~TM_F_MESOPIC;
 						/* set monitor transform */
 	if (monpri == NULL || monpri == stdprims || tmnew->flags & TM_F_BW) {
 		tmnew->monpri = stdprims;
@@ -86,226 +73,333 @@ 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->hbrmin = tmnew->hbrmax = 0;
+	tmnew->inpdat = NULL;
+	tmnew->hbrmin = 10; tmnew->hbrmax = -10;
 	tmnew->histo = NULL;
-	tmnew->mbrmin = tmnew->mbrmax = 0;
+	tmnew->mbrmin = 10; tmnew->mbrmax = -10;
 	tmnew->lumap = NULL;
 						/* zero private data */
 	for (i = TM_MAXPKG; i--; )
 		tmnew->pd[i] = NULL;
-						/* make tmnew current */
-	tmnew->tmprev = tmTop;
-	return(tmTop = tmnew);
+	tmnew->lastError = TM_E_OK;
+	tmnew->lastFunc = "NoErr";
+						/* return new TMstruct */
+	return(tmnew);
 }
 
 
 int
-tmSetSpace(pri, sf)		/* set input color space for conversions */
-RGBPRIMP	pri;
-double	sf;
+tmSetSpace(			/* set input color space for conversions */
+TMstruct	*tms,
+RGBPRIMP	pri,
+double	sf,
+MEM_PTR	dat
+)
 {
-	static char	funcName[] = "tmSetSpace";
-	register int	i, j;
+	static const char funcName[] = "tmSetSpace";
+	int	i, j;
 						/* error check */
-	if (tmTop == NULL)
+	if (tms == NULL)
 		returnErr(TM_E_TMINVAL);
 	if (sf <= 1e-12)
 		returnErr(TM_E_ILLEGAL);
 						/* check if no change */
-	if (pri == tmTop->inppri && FEQ(sf, tmTop->inpsf))
+	if (pri == tms->inppri && FEQ(sf, tms->inpsf) && dat == tms->inpdat)
 		returnOK;
-	tmTop->inppri = pri;			/* let's set it */
-	tmTop->inpsf = sf;
+	tms->inppri = pri;			/* let's set it */
+	tms->inpsf = sf;
+	tms->inpdat = dat;
 
-	if (tmTop->flags & TM_F_BW) {		/* color doesn't matter */
-		tmTop->monpri = tmTop->inppri;		/* eliminate xform */
-		if (tmTop->inppri == TM_XYZPRIM) {
-			tmTop->clf[CIEX] = tmTop->clf[CIEZ] = 0.;
-			tmTop->clf[CIEY] = 1.;
+	if (tms->flags & TM_F_BW) {		/* color doesn't matter */
+		tms->monpri = tms->inppri;		/* eliminate xform */
+		if (tms->inppri == TM_XYZPRIM) {
+			tms->clf[CIEX] = tms->clf[CIEZ] = 0.;
+			tms->clf[CIEY] = 1.;
 		} else {
-			comprgb2xyzmat(tmTop->cmat, tmTop->monpri);
-			tmTop->clf[RED] = tmTop->cmat[1][0];
-			tmTop->clf[GRN] = tmTop->cmat[1][1];
-			tmTop->clf[BLU] = tmTop->cmat[1][2];
+			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];
 		}
-		tmTop->cmat[0][0] = tmTop->cmat[1][1] = tmTop->cmat[2][2] =
-				tmTop->inpsf;
-		tmTop->cmat[0][1] = tmTop->cmat[0][2] = tmTop->cmat[1][0] =
-		tmTop->cmat[1][2] = tmTop->cmat[2][0] = tmTop->cmat[2][1] = 0.;
+		tms->cmat[0][0] = tms->cmat[1][1] = tms->cmat[2][2] =
+				tms->inpsf;
+		tms->cmat[0][1] = tms->cmat[0][2] = tms->cmat[1][0] =
+		tms->cmat[1][2] = tms->cmat[2][0] = tms->cmat[2][1] = 0.;
 
-	} else if (tmTop->inppri == TM_XYZPRIM)	/* input is XYZ */
-		compxyz2rgbmat(tmTop->cmat, tmTop->monpri);
+	} else if (tms->inppri == TM_XYZPRIM)	/* input is XYZ */
+		compxyz2rgbWBmat(tms->cmat, tms->monpri);
 
 	else {					/* input is RGB */
-		if (tmTop->inppri != tmTop->monpri &&
-				PRIMEQ(tmTop->inppri, tmTop->monpri))
-			tmTop->inppri = tmTop->monpri;	/* no xform */
-		comprgb2rgbmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
+		if (tms->inppri != tms->monpri &&
+				PRIMEQ(tms->inppri, tms->monpri))
+			tms->inppri = tms->monpri;	/* no xform */
+		comprgb2rgbWBmat(tms->cmat, tms->inppri, tms->monpri);
 	}
 	for (i = 0; i < 3; i++)
 		for (j = 0; j < 3; j++)
-			tmTop->cmat[i][j] *= tmTop->inpsf;
+			tms->cmat[i][j] *= tms->inpsf;
 						/* set color divisors */
 	for (i = 0; i < 3; i++)
-		if (tmTop->clf[i] > .001)
-			tmTop->cdiv[i] =
-				256.*pow(tmTop->clf[i], 1./tmTop->mongam);
+		if (tms->clf[i] > .001)
+			tms->cdiv[i] =
+				256.*pow(tms->clf[i], 1./tms->mongam);
 		else
-			tmTop->cdiv[i] = 1;
+			tms->cdiv[i] = 1;
 						/* notify packages */
 	for (i = tmNumPkgs; i--; )
-		if (tmTop->pd[i] != NULL && tmPkg[i]->NewSpace != NULL)
-			(*tmPkg[i]->NewSpace)(tmTop);
+		if (tms->pd[i] != NULL && tmPkg[i]->NewSpace != NULL)
+			(*tmPkg[i]->NewSpace)(tms);
 	returnOK;
 }
 
 
 void
-tmClearHisto()			/* clear current histogram */
+tmClearHisto(				/* clear current histogram */
+TMstruct	*tms
+)
 {
-	if (tmTop == NULL || tmTop->histo == NULL)
+	if (tms == NULL || tms->histo == NULL)
 		return;
-	free((MEM_PTR)tmTop->histo);
-	tmTop->histo = NULL;
+	free((MEM_PTR)tms->histo);
+	tms->histo = NULL;
 }
 
 
+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
-tmCvColors(ls, cs, scan, len)		/* convert float colors */
-TMbright	*ls;
-BYTE	*cs;
-COLOR	*scan;
-int	len;
+tmCvLums(				/* convert luminances using lookup */
+TMbright	*ls,
+float		*scan,
+int		len
+)
 {
-	static char	funcName[] = "tmCvColors";
-	static COLOR	csmall = {1e-6, 1e-6, 1e-6};
+	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,
+BYTE	*cs,
+COLOR	*scan,
+int	len
+)
+{
+	static const char funcName[] = "tmCvColors";
+	static BYTE	gamtab[1024];
+	static double	curgam = .0;
 	COLOR	cmon;
-	double	lum, slum;
-	register double	d;
-	register int	i;
+	float	lum, slum, d;
+	int	i;
 
-	if (tmTop == NULL)
+	if (tms == NULL)
 		returnErr(TM_E_TMINVAL);
-	if (ls == NULL | scan == NULL | len <= 0)
+	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(tmTop))		/* get monitor RGB */
-			colortrans(cmon, tmTop->cmat, scan[i]);
-		else {
-			cmon[RED] = tmTop->inpsf*scan[i][RED];
-			cmon[GRN] = tmTop->inpsf*scan[i][GRN];
-			cmon[BLU] = tmTop->inpsf*scan[i][BLU];
-		}
-							/* world luminance */
-		lum =	tmTop->clf[RED]*cmon[RED] +
-			tmTop->clf[GRN]*cmon[GRN] +
-			tmTop->clf[BLU]*cmon[BLU] ;
-							/* check range */
-		if (clipgamut(cmon, lum, CGAMUT_LOWER, csmall, cwhite))
-			lum =	tmTop->clf[RED]*cmon[RED] +
-				tmTop->clf[GRN]*cmon[GRN] +
-				tmTop->clf[BLU]*cmon[BLU] ;
-		if (lum < MINLUM) {
-			ls[i] = MINBRT-1;		/* bogus value */
-			lum = MINLUM;
+		if (tmNeedMatrix(tms)) {		/* get monitor RGB */
+			colortrans(cmon, tms->cmat, scan[i]);
 		} else {
-			d = TM_BRTSCALE*log(lum);	/* encode it */
-			ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
+			cmon[RED] = tms->inpsf*scan[i][RED];
+			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] ;
+		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 (tmTop->flags & TM_F_MESOPIC && lum < LMESUPPER) {
+		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 (tmTop->flags & TM_F_BW)
+				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;
 			}
-		} else if (tmTop->flags & TM_F_BW) {
+		} else if (tms->flags & TM_F_BW) {
 			cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
 		}
-		d = tmTop->clf[RED]*cmon[RED]/lum;
-		cs[3*i  ] = d>=.999 ? 255 :
-				(int)(256.*pow(d, 1./tmTop->mongam));
-		d = tmTop->clf[GRN]*cmon[GRN]/lum;
-		cs[3*i+1] = d>=.999 ? 255 :
-				(int)(256.*pow(d, 1./tmTop->mongam));
-		d = tmTop->clf[BLU]*cmon[BLU]/lum;
-		cs[3*i+2] = d>=.999 ? 255 :
-				(int)(256.*pow(d, 1./tmTop->mongam));
+		d = tms->clf[RED]*cmon[RED]/lum;
+		cs[3*i  ] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
+		d = tms->clf[GRN]*cmon[GRN]/lum;
+		cs[3*i+1] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
+		d = tms->clf[BLU]*cmon[BLU]/lum;
+		cs[3*i+2] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
 	}
 	returnOK;
 }
 
 
 int
-tmAddHisto(ls, len, wt)			/* add values to histogram */
-register TMbright	*ls;
-int	len;
-int	wt;
+tmAddHisto(				/* add values to histogram */
+TMstruct	*tms,
+TMbright	*ls,
+int	len,
+int	wt
+)
 {
-	static char	funcName[] = "tmAddHisto";
-	int	sum, oldorig, oldlen, horig, hlen;
-	register int	i, j;
+	static const char funcName[] = "tmAddHisto";
+	int	oldorig=0, oldlen, horig, hlen;
+	int	i, j;
 
-	if (len <= 0)
-		returnErr(TM_E_ILLEGAL);
-	if (tmTop == NULL)
+	if (tms == NULL)
 		returnErr(TM_E_TMINVAL);
+	if (len < 0)
+		returnErr(TM_E_ILLEGAL);
+	if (len == 0)
+		returnOK;
 						/* first, grow limits */
-	if (tmTop->histo == NULL) {
+	if (tms->histo == NULL) {
 		for (i = len; i-- && ls[i] < MINBRT; )
 			;
 		if (i < 0)
 			returnOK;
-		tmTop->hbrmin = tmTop->hbrmax = ls[i];
+		tms->hbrmin = tms->hbrmax = ls[i];
 		oldlen = 0;
 	} else {
-		oldorig = (tmTop->hbrmin-MINBRT)/HISTEP;
-		oldlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - oldorig;
+		oldorig = HISTI(tms->hbrmin);
+		oldlen = HISTI(tms->hbrmax) + 1 - oldorig;
 	}
 	for (i = len; i--; ) {
 		if ((j = ls[i]) < MINBRT)
 			continue;
-		if (j < tmTop->hbrmin)
-			tmTop->hbrmin = j;
-		else if (j > tmTop->hbrmax)
-			tmTop->hbrmax = j;
+		if (j < tms->hbrmin)
+			tms->hbrmin = j;
+		else if (j > tms->hbrmax)
+			tms->hbrmax = j;
 	}
-	horig = (tmTop->hbrmin-MINBRT)/HISTEP;
-	hlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - horig;
+	horig = HISTI(tms->hbrmin);
+	hlen = HISTI(tms->hbrmax) + 1 - horig;
 	if (hlen > oldlen) {			/* (re)allocate histogram */
-		register int	*newhist = (int *)calloc(hlen, sizeof(int));
+		int	*newhist = (int *)calloc(hlen, sizeof(int));
 		if (newhist == NULL)
 			returnErr(TM_E_NOMEM);
 		if (oldlen) {			/* copy and free old */
 			for (i = oldlen, j = i+oldorig-horig; i; )
-				newhist[--j] = tmTop->histo[--i];
-			free((MEM_PTR)tmTop->histo);
+				newhist[--j] = tms->histo[--i];
+			free((MEM_PTR)tms->histo);
 		}
-		tmTop->histo = newhist;
+		tms->histo = newhist;
 	}
 	if (wt == 0)
 		returnOK;
 	for (i = len; i--; )			/* add in new counts */
 		if (ls[i] >= MINBRT)
-			tmTop->histo[ (ls[i]-MINBRT)/HISTEP - horig ] += wt;
+			tms->histo[ HISTI(ls[i]) - horig ] += wt;
 	returnOK;
 }
 
 
 static double
-htcontrs(La)		/* human threshold contrast sensitivity, dL(La) */
-double	La;
+htcontrs(		/* human threshold contrast sensitivity, dL(La) */
+double	La
+)
 {
 	double	l10La, l10dL;
 				/* formula taken from Ferwerda et al. [SG96] */
@@ -326,144 +420,173 @@ double	La;
 
 
 int
-tmComputeMapping(gamval, Lddyn, Ldmax)
-double	gamval;
-double	Lddyn;
-double	Ldmax;
+tmFixedMapping(			/* compute fixed, linear tone-mapping */
+TMstruct	*tms,
+double	expmult,
+double	gamval
+)
 {
-	static char	funcName[] = "tmComputeMapping";
+	static const char funcName[] = "tmFixedMapping";
+	double		d;
+	int	i;
+	
+	if (!tmNewMap(tms))
+		returnErr(TM_E_NOMEM);
+	if (expmult <= .0)
+		expmult = 1.;
+	if (gamval < MINGAM)
+		gamval = tms->mongam;
+	d = log(expmult/tms->inpsf);
+	for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
+		double	val = 256. * exp(
+			( d + (tms->mbrmin+i)*(1./TM_BRTSCALE) )
+			/ gamval);
+		tms->lumap[i] = val >= (double)0xffff ? 0xffff : (int)val;
+	}
+	returnOK;
+}
+
+
+int
+tmComputeMapping(			/* compute histogram tone-mapping */
+TMstruct	*tms,
+double	gamval,
+double	Lddyn,
+double	Ldmax
+)
+{
+	static const char funcName[] = "tmComputeMapping";
 	int	*histo;
 	float	*cumf;
-	int	brt0, histlen, histot, threshold, ceiling, trimmings;
+	int	brt0, histlen, threshold, ceiling, trimmings;
 	double	logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
-	int4	sum;
-	register double	d;
-	register int	i, j;
+	int32	histot;
+	double	sum;
+	double	d;
+	int	i, j;
 
-	if (tmTop == NULL || tmTop->histo == NULL)
+	if (tms == NULL || tms->histo == NULL)
 		returnErr(TM_E_TMINVAL);
 					/* check arguments */
 	if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
 	if (Ldmax < MINLDMAX) Ldmax = DEFLDMAX;
-	if (gamval < MINGAM) gamval = tmTop->mongam;
+	if (gamval < MINGAM) gamval = tms->mongam;
 					/* compute handy values */
 	Ldmin = Ldmax/Lddyn;
 	logLddyn = log(Lddyn);
 	Ldavg = sqrt(Ldmax*Ldmin);
-	i = (tmTop->hbrmin-MINBRT)/HISTEP;
-	brt0 = MINBRT + HISTEP/2 + i*HISTEP;
-	histlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - i;
+	i = HISTI(tms->hbrmin);
+	brt0 = HISTV(i);
+	histlen = HISTI(tms->hbrmax) + 1 - i;
 					/* histogram total and mean */
 	histot = 0; sum = 0;
 	j = brt0 + histlen*HISTEP;
 	for (i = histlen; i--; ) {
-		histot += tmTop->histo[i];
-		sum += (j -= HISTEP) * tmTop->histo[i];
+		histot += tms->histo[i];
+		sum += (double)(j -= HISTEP) * tms->histo[i];
 	}
-	threshold = histot*.025 + .5;
-	if (threshold < 4)
+	threshold = histot*0.005 + .5;
+	if (!histot)
 		returnErr(TM_E_TMFAIL);
 	Lwavg = tmLuminance( (double)sum / histot );
-	if (!(tmTop->flags & TM_F_LINEAR)) {	/* clamp histogram */
-		histo = (int *)malloc(histlen*sizeof(int));
-		cumf = (float *)malloc((histlen+1)*sizeof(float));
-		if (histo == NULL | cumf == NULL)
-			returnErr(TM_E_NOMEM);
-		for (i = histlen; i--; )	/* make malleable copy */
-			histo[i] = tmTop->histo[i];
-		do {				/* iterate to solution */
-			sum = 0;		/* cumulative probability */
-			for (i = 0; i < histlen; i++) {
-				cumf[i] = (double)sum/histot;
-				sum += histo[i];
+					/* use linear tone mapping? */
+	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));
+	cumf = (float *)malloc((histlen+2)*sizeof(float));
+	if ((histo == NULL) | (cumf == NULL))
+		returnErr(TM_E_NOMEM);
+	cumf[histlen+1] = 1.;		/* guard for assignment code */
+	for (i = histlen; i--; )	/* make malleable copy */
+		histo[i] = tms->histo[i];
+	do {				/* iterate to solution */
+		sum = 0;		/* cumulative probability */
+		for (i = 0; i < histlen; i++) {
+			cumf[i] = (double)sum/histot;
+			sum += histo[i];
+		}
+		cumf[histlen] = 1.;
+		Tr = histot * (double)(tms->hbrmax - tms->hbrmin) /
+			((double)histlen*TM_BRTSCALE) / logLddyn;
+		ceiling = Tr + 1.;
+		trimmings = 0;		/* clip to envelope */
+		for (i = histlen; i--; ) {
+			if (tms->flags & TM_F_HCONTR) {
+				Lw = tmLuminance(brt0 + i*HISTEP);
+				Ld = Ldmin * exp( logLddyn *
+					.5*(cumf[i]+cumf[i+1]) );
+				ceiling = Tr * (htcontrs(Ld) * Lw) /
+					(htcontrs(Lw) * Ld) + 1.;
 			}
-			cumf[i] = 1.;
-			Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
-				((double)histlen*TM_BRTSCALE) / logLddyn;
-			ceiling = Tr + 1.;
-			trimmings = 0;			/* clip to envelope */
-			for (i = histlen; i--; ) {
-				if (tmTop->flags & TM_F_HCONTR) {
-					Lw = tmLuminance(brt0 + i*HISTEP);
-					Ld = Ldmin * exp( logLddyn *
-						.5*(cumf[i]+cumf[i+1]) );
-					ceiling = Tr * (htcontrs(Ld) * Lw) /
-						(htcontrs(Lw) * Ld) + 1.;
-				}
-				if (histo[i] > ceiling) {
-					trimmings += histo[i] - ceiling;
-					histo[i] = ceiling;
-				}
+			if (histo[i] > ceiling) {
+				trimmings += histo[i] - ceiling;
+				histo[i] = ceiling;
 			}
-		} while ((histot -= trimmings) > threshold &&
-						trimmings > threshold);
-	}
-						/* allocate luminance map */
-	if (tmTop->lumap != NULL)
-		free((MEM_PTR)tmTop->lumap);
-	tmTop->mbrmin = tmTop->hbrmin;
-	tmTop->mbrmax = tmTop->hbrmax;
-	tmTop->lumap = (unsigned short *)malloc(
-		(tmTop->mbrmax-tmTop->mbrmin+1)*sizeof(unsigned short) );
-	if (tmTop->lumap == NULL)
-		returnErr(TM_E_NOMEM);
-	if (tmTop->flags & TM_F_LINEAR || histot <= threshold) {
-						/* linear tone mapping */
-		if (tmTop->flags & TM_F_HCONTR)
-			d = htcontrs(Ldavg) / htcontrs(Lwavg);
-		else
-			d = Ldavg / Lwavg;
-		d = log(d/Ldmax);
-		for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
-			tmTop->lumap[i] = 256. * exp(
-				( d + (tmTop->mbrmin+i)/(double)TM_BRTSCALE )
-				/ gamval );
-	} else {
-						/* histogram adjustment */
-		for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; ) {
-			j = d = (double)i/(tmTop->mbrmax-tmTop->mbrmin)*histlen;
-			d -= (double)j;
-			Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
-			d = (Ld - Ldmin)/(Ldmax - Ldmin);
-			tmTop->lumap[i] = 256.*pow(d, 1./gamval);
 		}
+					/* check if we're out of data */
+		if ((histot -= trimmings) <= threshold) {
+			free((MEM_PTR)histo);
+			free((MEM_PTR)cumf);
+			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;
+		d -= (double)j;
+		Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
+		d = (Ld - Ldmin)/(Ldmax - Ldmin);
+		tms->lumap[i] = 256.*pow(d, 1./gamval);
 	}
-	if (!(tmTop->flags & TM_F_LINEAR)) {
-		free((MEM_PTR)histo);
-		free((MEM_PTR)cumf);
-	}
+	free((MEM_PTR)histo);		/* clean up and return */
+	free((MEM_PTR)cumf);
 	returnOK;
+linearmap:				/* linear tone-mapping */
+	if (tms->flags & TM_F_HCONTR)
+		d = htcontrs(Ldavg) / htcontrs(Lwavg);
+	else
+		d = Ldavg / Lwavg;
+	return(tmFixedMapping(tms, tms->inpsf*d/Ldmax, gamval));
 }
 
 
 int
-tmMapPixels(ps, ls, cs, len)
-register BYTE	*ps;
-TMbright	*ls;
-register BYTE	*cs;
-int	len;
+tmMapPixels(			/* apply tone-mapping to pixel(s) */
+TMstruct	*tms,
+BYTE	*ps,
+TMbright	*ls,
+BYTE	*cs,
+int	len
+)
 {
-	static char	funcName[] = "tmMapPixels";
-	register int4	li, pv;
+	static const char funcName[] = "tmMapPixels";
+	int32	li, pv;
 
-	if (tmTop == NULL || tmTop->lumap == NULL)
+	if (tms == NULL || tms->lumap == NULL)
 		returnErr(TM_E_TMINVAL);
-	if (ps == NULL | ls == NULL | len <= 0)
+	if ((ps == NULL) | (ls == NULL) | (len < 0))
 		returnErr(TM_E_ILLEGAL);
 	while (len--) {
-		if ((li = *ls++) < tmTop->mbrmin)
-			li = tmTop->mbrmin;
-		else if (li > tmTop->mbrmax)
-			li = tmTop->mbrmax;
-		li = tmTop->lumap[li - tmTop->mbrmin];
+		if ((li = *ls++) < tms->mbrmin) {
+			li = 0;
+		} else {
+			if (li > tms->mbrmax)
+				li = tms->mbrmax;
+			li = tms->lumap[li - tms->mbrmin];
+		}
 		if (cs == TM_NOCHROM)
 			*ps++ = li>255 ? 255 : li;
 		else {
-			pv = *cs++ * li / tmTop->cdiv[RED];
+			pv = *cs++ * li / tms->cdiv[RED];
 			*ps++ = pv>255 ? 255 : pv;
-			pv = *cs++ * li / tmTop->cdiv[GRN];
+			pv = *cs++ * li / tms->cdiv[GRN];
 			*ps++ = pv>255 ? 255 : pv;
-			pv = *cs++ * li / tmTop->cdiv[BLU];
+			pv = *cs++ * li / tms->cdiv[BLU];
 			*ps++ = pv>255 ? 255 : pv;
 		}
 	}
@@ -471,60 +594,27 @@ int	len;
 }
 
 
-struct tmStruct *
-tmPop()				/* pop top tone mapping off stack */
+TMstruct *
+tmDup(				/* duplicate top tone mapping */
+TMstruct	*tms
+)
 {
-	register struct tmStruct	*tms;
-
-	if ((tms = tmTop) != NULL)
-		tmTop = tms->tmprev;
-	return(tms);
-}
-
-
-int
-tmPull(tms)			/* pull a tone mapping from stack */
-register struct tmStruct	*tms;
-{
-	register struct tmStruct	*tms2;
-					/* special cases first */
-	if (tms == NULL | tmTop == NULL)
-		return(0);
-	if (tms == tmTop) {
-		tmTop = tms->tmprev;
-		tms->tmprev = NULL;
-		return(1);
-	}
-	for (tms2 = tmTop; tms2->tmprev != NULL; tms2 = tms2->tmprev)
-		if (tms == tms2->tmprev) {	/* remove it */
-			tms2->tmprev = tms->tmprev;
-			tms->tmprev = NULL;
-			return(1);
-		}
-	return(0);			/* not found on stack */
-}
-
-
-struct tmStruct *
-tmDup()				/* duplicate top tone mapping */
-{
 	int	len;
-	register int	i;
-	register struct tmStruct	*tmnew;
+	int	i;
+	TMstruct	*tmnew;
 
-	if (tmTop == NULL)		/* anything to duplicate? */
+	if (tms == NULL)		/* anything to duplicate? */
 		return(NULL);
-	tmnew = (struct tmStruct *)malloc(sizeof(struct tmStruct));
+	tmnew = (TMstruct *)malloc(sizeof(TMstruct));
 	if (tmnew == NULL)
 		return(NULL);
-	*tmnew = *tmTop;		/* copy everything */
+	*tmnew = *tms;		/* copy everything */
 	if (tmnew->histo != NULL) {	/* duplicate histogram */
-		len = (tmnew->hbrmax-MINBRT)/HISTEP + 1 -
-				(tmnew->hbrmin-MINBRT)/HISTEP;
+		len = HISTI(tmnew->hbrmax) + 1 - HISTI(tmnew->hbrmin);
 		tmnew->histo = (int *)malloc(len*sizeof(int));
 		if (tmnew->histo != NULL)
 			for (i = len; i--; )
-				tmnew->histo[i] = tmTop->histo[i];
+				tmnew->histo[i] = tms->histo[i];
 	}
 	if (tmnew->lumap != NULL) {	/* duplicate luminance mapping */
 		len = tmnew->mbrmax-tmnew->mbrmin+1;
@@ -532,45 +622,24 @@ tmDup()				/* duplicate top tone mapping */
 						len*sizeof(unsigned short) );
 		if (tmnew->lumap != NULL)
 			for (i = len; i--; )
-				tmnew->lumap[i] = tmTop->lumap[i];
+				tmnew->lumap[i] = tms->lumap[i];
 	}
 					/* clear package data */
 	for (i = tmNumPkgs; i--; )
 		tmnew->pd[i] = NULL;
-	tmnew->tmprev = tmTop;		/* make copy current */
-	return(tmTop = tmnew);
+					/* return copy */
+	return(tmnew);
 }
 
 
-int
-tmPush(tms)			/* push tone mapping on top of stack */
-register struct tmStruct	*tms;
-{
-	static char	funcName[] = "tmPush";
-					/* check validity */
-	if (tms == NULL)
-		returnErr(TM_E_ILLEGAL);
-	if (tms == tmTop)		/* check necessity */
-		returnOK;
-					/* pull if already in stack */
-	(void)tmPull(tms);
-					/* push it on top */
-	tms->tmprev = tmTop;
-	tmTop = tms;
-	returnOK;
-}
-
-
 void
 tmDone(tms)			/* done with tone mapping -- destroy it */
-register struct tmStruct	*tms;
+TMstruct	*tms;
 {
-	register int	i;
-					/* NULL arg. is equiv. to tmTop */
-	if (tms == NULL && (tms = tmTop) == NULL)
+	int	i;
+					/* NULL arg. is equiv. to tms */
+	if (tms == NULL)
 		return;
-					/* take out of stack if present */
-	(void)tmPull(tms);
 					/* free tables */
 	if (tms->histo != NULL)
 		free((MEM_PTR)tms->histo);
@@ -581,4 +650,64 @@ register struct tmStruct	*tms;
 		if (tms->pd[i] != NULL)
 			(*tmPkg[i]->Free)(tms->pd[i]);
 	free((MEM_PTR)tms);		/* free basic structure */
+}
+
+/******************** Shared but Private library routines *********************/
+
+BYTE	tmMesofact[BMESUPPER-BMESLOWER];
+
+void
+tmMkMesofact()				/* build mesopic lookup factor table */
+{
+	int	i;
+
+	if (tmMesofact[BMESUPPER-BMESLOWER-1])
+		return;
+
+	for (i = BMESLOWER; i < BMESUPPER; i++)
+		tmMesofact[i-BMESLOWER] = 256. *
+				(tmLuminance(i) - LMESLOWER) /
+				(LMESUPPER - LMESLOWER);
+}
+
+
+int
+tmNewMap(			/* allocate new tone-mapping array */
+TMstruct	*tms
+)
+{
+	if (tms->lumap != NULL && (tms->mbrmax - tms->mbrmin) !=
+					(tms->hbrmax - tms->hbrmin)) {
+		free((MEM_PTR)tms->lumap);
+		tms->lumap = NULL;
+	}
+	tms->mbrmin = tms->hbrmin;
+	tms->mbrmax = tms->hbrmax;
+	if (tms->mbrmin > tms->mbrmax)
+		return 0;
+	if (tms->lumap == NULL)
+		tms->lumap = (unsigned short *)malloc(sizeof(unsigned short)*
+					(tms->mbrmax-tms->mbrmin+1));
+	return(tms->lumap != NULL);
+}
+
+
+int
+tmErrorReturn(				/* error return (with message) */
+const char	*func,
+TMstruct	*tms,
+int	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);
+	fputs("!\n", stderr);
+	return(err);
 }