--- ray/src/common/tonemap.c 1997/11/17 14:02:12 3.6
+++ ray/src/common/tonemap.c 2003/02/22 02:07:22 3.9
@@ -1,14 +1,71 @@
-/* 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.9 2003/02/22 02:07:22 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
*/
+/* ====================================================================
+ * The Radiance Software License, Version 1.0
+ *
+ * Copyright (c) 1990 - 2002 The Regents of the University of California,
+ * through Lawrence Berkeley National Laboratory. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in
+ * the documentation and/or other materials provided with the
+ * distribution.
+ *
+ * 3. The end-user documentation included with the redistribution,
+ * if any, must include the following acknowledgment:
+ * "This product includes Radiance software
+ * (http://radsite.lbl.gov/)
+ * developed by the Lawrence Berkeley National Laboratory
+ * (http://www.lbl.gov/)."
+ * Alternately, this acknowledgment may appear in the software itself,
+ * if and wherever such third-party acknowledgments normally appear.
+ *
+ * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
+ * and "The Regents of the University of California" must
+ * not be used to endorse or promote products derived from this
+ * software without prior written permission. For written
+ * permission, please contact radiance@radsite.lbl.gov.
+ *
+ * 5. Products derived from this software may not be called "Radiance",
+ * nor may "Radiance" appear in their name, without prior written
+ * permission of Lawrence Berkeley National Laboratory.
+ *
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
+ * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+ * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+ * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+ * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ * ====================================================================
+ *
+ * This software consists of voluntary contributions made by many
+ * individuals on behalf of Lawrence Berkeley National Laboratory. For more
+ * information on Lawrence Berkeley National Laboratory, please see
+ * .
+ */
+
#include
#include
#include "tmprivat.h"
@@ -26,30 +83,12 @@ int tmLastError; /* last error incurred by library *
char *tmLastFunction; /* error-generating function name */
-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;
{
- static char funcName[] = "tmInit";
COLORMAT cmat;
register struct tmStruct *tmnew;
register int i;
@@ -66,7 +105,7 @@ 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];
@@ -86,9 +125,9 @@ 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->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--; )
@@ -123,7 +162,7 @@ double sf;
tmTop->clf[CIEX] = tmTop->clf[CIEZ] = 0.;
tmTop->clf[CIEY] = 1.;
} 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];
@@ -134,13 +173,13 @@ 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++)
@@ -186,7 +225,7 @@ 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 (tmNeedMatrix(tmTop)) /* get monitor RGB */
@@ -248,19 +287,46 @@ int len;
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; )
@@ -325,7 +391,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;
@@ -334,9 +443,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;
+ int4 histot;
+ double sum;
register double d;
register int i, j;
@@ -364,75 +474,68 @@ double Ldmax;
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->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);
+ /* 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((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 (tmTop->flags & TM_F_HCONTR)
+ d = htcontrs(Ldavg) / htcontrs(Lwavg);
+ else
+ d = Ldavg / Lwavg;
+ return(tmFixedMapping(tmTop->inpsf*d/Ldmax, gamval));
}
@@ -448,14 +551,16 @@ int len;
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);
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++) < 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 {
@@ -581,4 +686,39 @@ 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 */
+{
+ 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);
}