--- ray/src/common/tonemap.c	1997/04/16 20:28:03	3.3
+++ ray/src/common/tonemap.c	2003/09/18 20:18:19	3.14
@@ -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.14 2003/09/18 20:18:19 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"
@@ -18,20 +20,12 @@ static char SCCSid[] = "$SunId$ LBL";
 
 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
-tmErrorReturn(func, err)		/* error return (with message) */
-char	*func;
-int	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);
-}
+int	tmLastError;			/* last error incurred by library */
+char	*tmLastFunction;		/* error-generating function name */
 
 
 struct tmStruct *
@@ -40,7 +34,6 @@ int	flags;
 RGBPRIMP	monpri;
 double	gamval;
 {
-	static char	funcName[] = "tmInit";
 	COLORMAT	cmat;
 	register struct tmStruct	*tmnew;
 	register int	i;
@@ -57,37 +50,34 @@ double	gamval;
 		tmnew->clf[GRN] = rgb2xyzmat[1][1];
 		tmnew->clf[BLU] = rgb2xyzmat[1][2];
 	} else {
-		comprgb2xyzmat(cmat, tmnew->monpri=monpri);
+		comprgb2xyzWBmat(cmat, tmnew->monpri=monpri);
 		tmnew->clf[RED] = cmat[1][0];
 		tmnew->clf[GRN] = cmat[1][1];
 		tmnew->clf[BLU] = cmat[1][2];
 	}
-	tmnew->clfb[RED] = 256.*tmnew->clf[RED] + .5;
-	tmnew->clfb[GRN] = 256.*tmnew->clf[GRN] + .5;
-	tmnew->clfb[BLU] = 256.*tmnew->clf[BLU] + .5;
-	tmnew->clfb[EXP] = COLXS;
 						/* set gamma value */
 	if (gamval < MINGAM)
 		tmnew->mongam = DEFGAM;
 	else
 		tmnew->mongam = gamval;
-	for (i = TM_GAMTSZ; i--; )
-		tmnew->gamb[i] = 256.*pow((i+.5)/TM_GAMTSZ, 1./tmnew->mongam);
+						/* set color divisors */
+	for (i = 0; i < 3; i++)
+		tmnew->cdiv[i] = 256.*pow(tmnew->clf[i], 1./tmnew->mongam);
+
 						/* set input transform */
 	tmnew->inppri = tmnew->monpri;
 	tmnew->cmat[0][0] = tmnew->cmat[1][1] = tmnew->cmat[2][2] =
 			tmnew->inpsf = WHTEFFICACY;
-	tmnew->inpsfb = TM_BRTSCALE*log(tmnew->inpsf) + .5;
 	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->flags &= ~TM_F_NEEDMAT;
-	tmnew->brmin = tmnew->brmax = 0;
+	tmnew->hbrmin = 10; tmnew->hbrmax = -10;
 	tmnew->histo = NULL;
+	tmnew->mbrmin = 10; tmnew->mbrmax = -10;
 	tmnew->lumap = NULL;
-	tmnew->tmprev = NULL;
-
-	tmnew->flags |= TM_F_INITED;
-						/* make it current */
+						/* zero private data */
+	for (i = TM_MAXPKG; i--; )
+		tmnew->pd[i] = NULL;
+						/* make tmnew current */
 	tmnew->tmprev = tmTop;
 	return(tmTop = tmnew);
 }
@@ -110,17 +100,14 @@ double	sf;
 		returnOK;
 	tmTop->inppri = pri;			/* let's set it */
 	tmTop->inpsf = sf;
-	tmTop->inpsfb = TM_BRTSCALE*log(sf) + (sf>=1. ? .5 : -.5);
 
 	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.;
-			tmTop->clfb[CIEX] = tmTop->clfb[CIEZ] = 0;
-			tmTop->clfb[CIEY] = 255;
 		} else {
-			comprgb2xyzmat(tmTop->cmat, tmTop->monpri);
+			comprgb2xyzWBmat(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];
@@ -131,21 +118,28 @@ double	sf;
 		tmTop->cmat[1][2] = tmTop->cmat[2][0] = tmTop->cmat[2][1] = 0.;
 
 	} else if (tmTop->inppri == TM_XYZPRIM)	/* input is XYZ */
-		compxyz2rgbmat(tmTop->cmat, tmTop->monpri);
+		compxyz2rgbWBmat(tmTop->cmat, tmTop->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);
+		comprgb2rgbWBmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
 	}
 	for (i = 0; i < 3; i++)
 		for (j = 0; j < 3; j++)
 			tmTop->cmat[i][j] *= tmTop->inpsf;
-	if (tmTop->inppri != tmTop->monpri)
-		tmTop->flags |= TM_F_NEEDMAT;
-	else
-		tmTop->flags &= ~TM_F_NEEDMAT;
+						/* 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);
+		else
+			tmTop->cdiv[i] = 1;
+						/* notify packages */
+	for (i = tmNumPkgs; i--; )
+		if (tmTop->pd[i] != NULL && tmPkg[i]->NewSpace != NULL)
+			(*tmPkg[i]->NewSpace)(tmTop);
 	returnOK;
 }
 
@@ -155,7 +149,7 @@ tmClearHisto()			/* clear current histogram */
 {
 	if (tmTop == NULL || tmTop->histo == NULL)
 		return;
-	free((char *)tmTop->histo);
+	free((MEM_PTR)tmTop->histo);
 	tmTop->histo = NULL;
 }
 
@@ -168,6 +162,7 @@ COLOR	*scan;
 int	len;
 {
 	static char	funcName[] = "tmCvColors";
+	static COLOR	csmall = {.5*MINLUM, .5*MINLUM, .5*MINLUM};
 	COLOR	cmon;
 	double	lum, slum;
 	register double	d;
@@ -175,10 +170,10 @@ int	len;
 
 	if (tmTop == NULL)
 		returnErr(TM_E_TMINVAL);
-	if (ls == NULL | scan == NULL | len <= 0)
+	if ((ls == NULL) | (scan == NULL) | (len < 0))
 		returnErr(TM_E_ILLEGAL);
 	for (i = len; i--; ) {
-		if (tmTop->flags & TM_F_NEEDMAT)	/* get monitor RGB */
+		if (tmNeedMatrix(tmTop))		/* get monitor RGB */
 			colortrans(cmon, tmTop->cmat, scan[i]);
 		else {
 			cmon[RED] = tmTop->inpsf*scan[i][RED];
@@ -190,7 +185,7 @@ int	len;
 			tmTop->clf[GRN]*cmon[GRN] +
 			tmTop->clf[BLU]*cmon[BLU] ;
 							/* check range */
-		if (clipgamut(cmon, lum, CGAMUT_LOWER, cblack, cwhite))
+		if (clipgamut(cmon, lum, CGAMUT_LOWER, csmall, cwhite))
 			lum =	tmTop->clf[RED]*cmon[RED] +
 				tmTop->clf[GRN]*cmon[GRN] +
 				tmTop->clf[BLU]*cmon[BLU] ;
@@ -223,55 +218,82 @@ int	len;
 			cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
 		}
 		d = tmTop->clf[RED]*cmon[RED]/lum;
-		/* cs[3*i  ] = d>.999 ? 255 : 256.*pow(d, 1./tmTop->mongam); */
-		cs[3*i  ] = d>.999 ? 255 : tmTop->gamb[(int)(d*TM_GAMTSZ)];
+		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 : 256.*pow(d, 1./tmTop->mongam); */
-		cs[3*i+1] = d>.999 ? 255 : tmTop->gamb[(int)(d*TM_GAMTSZ)];
+		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 : 256.*pow(d, 1./tmTop->mongam); */
-		cs[3*i+2] = d>.999 ? 255 : tmTop->gamb[(int)(d*TM_GAMTSZ)];
+		cs[3*i+2] = d>=.999 ? 255 :
+				(int)(256.*pow(d, 1./tmTop->mongam));
 	}
 	returnOK;
 }
 
 
 int
+tmCvGrays(ls, scan, len)		/* convert float gray values */
+TMbright	*ls;
+float	*scan;
+int	len;
+{
+	static char	funcName[] = "tmCvGrays";
+	register double	d;
+	register int	i;
+
+	if (tmTop == 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(ls, len, wt)			/* add values to histogram */
 register TMbright	*ls;
 int	len;
 int	wt;
 {
 	static char	funcName[] = "tmAddHisto";
-	int	sum, oldorig, oldlen, horig, hlen;
+	int	oldorig=0, oldlen, horig, hlen;
 	register int	i, j;
 
-	if (len <= 0)
-		returnErr(TM_E_ILLEGAL);
 	if (tmTop == NULL)
 		returnErr(TM_E_TMINVAL);
+	if (len < 0)
+		returnErr(TM_E_ILLEGAL);
+	if (len == 0)
+		returnOK;
 						/* first, grow limits */
 	if (tmTop->histo == NULL) {
 		for (i = len; i-- && ls[i] < MINBRT; )
 			;
 		if (i < 0)
 			returnOK;
-		tmTop->brmin = tmTop->brmax = ls[i];
+		tmTop->hbrmin = tmTop->hbrmax = ls[i];
 		oldlen = 0;
 	} else {
-		oldorig = (tmTop->brmin-MINBRT)/HISTEP;
-		oldlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - oldorig;
+		oldorig = (tmTop->hbrmin-MINBRT)/HISTEP;
+		oldlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - oldorig;
 	}
 	for (i = len; i--; ) {
 		if ((j = ls[i]) < MINBRT)
 			continue;
-		if (j < tmTop->brmin)
-			tmTop->brmin = j;
-		else if (j > tmTop->brmax)
-			tmTop->brmax = j;
+		if (j < tmTop->hbrmin)
+			tmTop->hbrmin = j;
+		else if (j > tmTop->hbrmax)
+			tmTop->hbrmax = j;
 	}
-	horig = (tmTop->brmin-MINBRT)/HISTEP;
-	hlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - horig;
+	horig = (tmTop->hbrmin-MINBRT)/HISTEP;
+	hlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - horig;
 	if (hlen > oldlen) {			/* (re)allocate histogram */
 		register int	*newhist = (int *)calloc(hlen, sizeof(int));
 		if (newhist == NULL)
@@ -279,13 +301,9 @@ int	wt;
 		if (oldlen) {			/* copy and free old */
 			for (i = oldlen, j = i+oldorig-horig; i; )
 				newhist[--j] = tmTop->histo[--i];
-			free((char *)tmTop->histo);
+			free((MEM_PTR)tmTop->histo);
 		}
 		tmTop->histo = newhist;
-		if (tmTop->lumap != NULL) {	/* invalid tone map */
-			free((char *)tmTop->lumap);
-			tmTop->lumap = NULL;
-		}
 	}
 	if (wt == 0)
 		returnOK;
@@ -318,7 +336,50 @@ double	La;
 }
 
 
+static int
+tmNewMap()
+{
+	if (tmTop->lumap != NULL && (tmTop->mbrmax - tmTop->mbrmin) !=
+					(tmTop->hbrmax - tmTop->hbrmin)) {
+		free((MEM_PTR)tmTop->lumap);
+		tmTop->lumap = NULL;
+	}
+	tmTop->mbrmin = tmTop->hbrmin;
+	tmTop->mbrmax = tmTop->hbrmax;
+	if (tmTop->mbrmin > tmTop->mbrmax)
+		return 0;
+	if (tmTop->lumap == NULL)
+		tmTop->lumap = (unsigned short *)malloc(sizeof(unsigned short)*
+					(tmTop->mbrmax-tmTop->mbrmin+1));
+	return(tmTop->lumap != NULL);
+}
+
+
 int
+tmFixedMapping(expmult, gamval)
+double	expmult;
+double	gamval;
+{
+	static char	funcName[] = "tmFixedMapping";
+	double		d;
+	register int	i;
+	
+	if (!tmNewMap())
+		returnErr(TM_E_NOMEM);
+	if (expmult <= .0)
+		expmult = 1.;
+	if (gamval < MINGAM)
+		gamval = tmTop->mongam;
+	d = log(expmult/tmTop->inpsf);
+	for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
+		tmTop->lumap[i] = 256. * exp(
+			( d + (tmTop->mbrmin+i)*(1./TM_BRTSCALE) )
+			/ gamval );
+	returnOK;
+}
+
+
+int
 tmComputeMapping(gamval, Lddyn, Ldmax)
 double	gamval;
 double	Lddyn;
@@ -327,9 +388,10 @@ double	Ldmax;
 	static 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;
+	int32	histot;
+	double	sum;
 	register double	d;
 	register int	i, j;
 
@@ -343,9 +405,9 @@ double	Ldmax;
 	Ldmin = Ldmax/Lddyn;
 	logLddyn = log(Lddyn);
 	Ldavg = sqrt(Ldmax*Ldmin);
-	i = (tmTop->brmin-MINBRT)/HISTEP;
+	i = (tmTop->hbrmin-MINBRT)/HISTEP;
 	brt0 = MINBRT + HISTEP/2 + i*HISTEP;
-	histlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - i;
+	histlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - i;
 					/* histogram total and mean */
 	histot = 0; sum = 0;
 	j = brt0 + histlen*HISTEP;
@@ -353,77 +415,72 @@ double	Ldmax;
 		histot += tmTop->histo[i];
 		sum += (j -= HISTEP) * tmTop->histo[i];
 	}
-	threshold = histot*.025 + .5;
+	threshold = histot*0.005 + .5;
 	if (threshold < 4)
 		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];
+					/* allocate space for mapping */
+	if (!tmNewMap())
+		returnErr(TM_E_NOMEM);
+					/* use linear tone mapping? */
+	if (tmTop->flags & TM_F_LINEAR)
+		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] = 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];
+		}
+		cumf[histlen] = 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.;
 			}
-			cumf[i] = 1.;
-			Tr = histot * (double)(tmTop->brmax - tmTop->brmin) /
-				((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) {
-		tmTop->lumap = (unsigned short *)malloc(
-			(tmTop->brmax-tmTop->brmin+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->brmax-tmTop->brmin+1; i--; )
-			tmTop->lumap[i] = 256. * exp(
-				( d + (tmTop->brmin+i)/(double)TM_BRTSCALE )
-				/ gamval );
-	} else {
-						/* histogram adjustment */
-		for (i = tmTop->brmax-tmTop->brmin+1; i--; ) {
-			j = d = (double)i/(tmTop->brmax-tmTop->brmin)*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);
+					/* assign tone-mapping */
+	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);
 	}
-	if (!(tmTop->flags & TM_F_LINEAR)) {
-		free((char *)histo);
-		free((char *)cumf);
-	}
+	free((MEM_PTR)histo);		/* clean up and return */
+	free((MEM_PTR)cumf);
 	returnOK;
+linearmap:				/* linear tone-mapping */
+	if (tmTop->flags & TM_F_HCONTR)
+		d = htcontrs(Ldavg) / htcontrs(Lwavg);
+	else
+		d = Ldavg / Lwavg;
+	return(tmFixedMapping(tmTop->inpsf*d/Ldmax, gamval));
 }
 
 
@@ -435,30 +492,28 @@ register BYTE	*cs;
 int	len;
 {
 	static char	funcName[] = "tmMapPixels";
-	int	rdiv, gdiv, bdiv;
-	register int	li, pv;
+	register int32	li, pv;
 
 	if (tmTop == NULL || tmTop->lumap == NULL)
 		returnErr(TM_E_TMINVAL);
-	if (ps == NULL | ls == NULL | len <= 0)
+	if ((ps == NULL) | (ls == NULL) | (len < 0))
 		returnErr(TM_E_ILLEGAL);
-	rdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[RED])>>8];
-	gdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[GRN])>>8];
-	bdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[BLU])>>8];
 	while (len--) {
-		if ((li = *ls++) < tmTop->brmin)
-			li = tmTop->brmin;
-		else if (li > tmTop->brmax)
-			li = tmTop->brmax;
-		li = tmTop->lumap[li - tmTop->brmin];
+		if ((li = *ls++) < tmTop->mbrmin) {
+			li = 0;
+		} else {
+			if (li > tmTop->mbrmax)
+				li = tmTop->mbrmax;
+			li = tmTop->lumap[li - tmTop->mbrmin];
+		}
 		if (cs == TM_NOCHROM)
 			*ps++ = li>255 ? 255 : li;
 		else {
-			pv = *cs++ * li / rdiv;
+			pv = *cs++ * li / tmTop->cdiv[RED];
 			*ps++ = pv>255 ? 255 : pv;
-			pv = *cs++ * li / gdiv;
+			pv = *cs++ * li / tmTop->cdiv[GRN];
 			*ps++ = pv>255 ? 255 : pv;
-			pv = *cs++ * li / bdiv;
+			pv = *cs++ * li / tmTop->cdiv[BLU];
 			*ps++ = pv>255 ? 255 : pv;
 		}
 	}
@@ -483,7 +538,7 @@ register struct tmStruct	*tms;
 {
 	register struct tmStruct	*tms2;
 					/* special cases first */
-	if (tms == NULL | tmTop == NULL)
+	if ((tms == NULL) | (tmTop == NULL))
 		return(0);
 	if (tms == tmTop) {
 		tmTop = tms->tmprev;
@@ -514,21 +569,24 @@ tmDup()				/* duplicate top tone mapping */
 		return(NULL);
 	*tmnew = *tmTop;		/* copy everything */
 	if (tmnew->histo != NULL) {	/* duplicate histogram */
-		len = (tmnew->brmax-MINBRT)/HISTEP + 1 -
-				(tmnew->brmin-MINBRT)/HISTEP;
+		len = (tmnew->hbrmax-MINBRT)/HISTEP + 1 -
+				(tmnew->hbrmin-MINBRT)/HISTEP;
 		tmnew->histo = (int *)malloc(len*sizeof(int));
 		if (tmnew->histo != NULL)
 			for (i = len; i--; )
 				tmnew->histo[i] = tmTop->histo[i];
 	}
 	if (tmnew->lumap != NULL) {	/* duplicate luminance mapping */
-		len = tmnew->brmax-tmnew->brmin+1;
+		len = tmnew->mbrmax-tmnew->mbrmin+1;
 		tmnew->lumap = (unsigned short *)malloc(
 						len*sizeof(unsigned short) );
 		if (tmnew->lumap != NULL)
 			for (i = len; i--; )
 				tmnew->lumap[i] = tmTop->lumap[i];
 	}
+					/* clear package data */
+	for (i = tmNumPkgs; i--; )
+		tmnew->pd[i] = NULL;
 	tmnew->tmprev = tmTop;		/* make copy current */
 	return(tmTop = tmnew);
 }
@@ -540,7 +598,7 @@ register struct tmStruct	*tms;
 {
 	static char	funcName[] = "tmPush";
 					/* check validity */
-	if (tms == NULL || !(tms->flags & TM_F_INITED))
+	if (tms == NULL)
 		returnErr(TM_E_ILLEGAL);
 	if (tms == tmTop)		/* check necessity */
 		returnOK;
@@ -557,6 +615,7 @@ void
 tmDone(tms)			/* done with tone mapping -- destroy it */
 register struct tmStruct	*tms;
 {
+	register int	i;
 					/* NULL arg. is equiv. to tmTop */
 	if (tms == NULL && (tms = tmTop) == NULL)
 		return;
@@ -564,9 +623,47 @@ register struct tmStruct	*tms;
 	(void)tmPull(tms);
 					/* free tables */
 	if (tms->histo != NULL)
-		free((char *)tms->histo);
+		free((MEM_PTR)tms->histo);
 	if (tms->lumap != NULL)
-		free((char *)tms->lumap);
-	tms->flags = 0;
-	free((char *)tms);		/* free basic structure */
+		free((MEM_PTR)tms->lumap);
+					/* free private data */
+	for (i = tmNumPkgs; i--; )
+		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 */
+{
+	register int	i;
+
+	if (tmMesofact[BMESUPPER-BMESLOWER-1])
+		return;
+
+	for (i = BMESLOWER; i < BMESUPPER; i++)
+		tmMesofact[i-BMESLOWER] = 256. *
+				(tmLuminance(i) - LMESLOWER) /
+				(LMESUPPER - LMESLOWER);
+}
+
+
+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);
 }