src/px/pcond3.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 * Routines for computing and applying brightness mapping.
 */

#include "pcond.h"


#define CVRATIO         0.025           /* fraction of samples allowed > env. */

#define LN_10           2.30258509299404568402
#define exp10(x)        exp(LN_10*(x))

double  modhist[HISTRES];               /* modified histogram */
double  mhistot;                        /* modified histogram total */
double  cumf[HISTRES+1];                /* cumulative distribution function */


getfixations(fp)                /* load fixation history list */
FILE    *fp;
{
#define FIXHUNK         128
        RESOLU  fvres;
        int     pos[2];
        register int    px, py, i;
                                /* initialize our resolution struct */
        if ((fvres.rt=inpres.rt)&YMAJOR) {
                fvres.xr = fvxr;
                fvres.yr = fvyr;
        } else {
                fvres.xr = fvyr;
                fvres.yr = fvxr;
        }
                                /* read each picture position */
        while (fscanf(fp, "%d %d", &pos[0], &pos[1]) == 2) {
                                /* convert to closest index in foveal image */
                loc2pix(pos, &fvres,
                                (pos[0]+.5)/inpres.xr, (pos[1]+.5)/inpres.yr);
                                /* include nine neighborhood samples */
                for (px = pos[0]-1; px <= pos[0]+1; px++) {
                        if (px < 0 || px >= fvxr)
                                continue;
                        for (py = pos[1]-1; py <= pos[1]+1; py++) {
                                if (py < 0 || py >= fvyr)
                                        continue;
                                for (i = nfixations; i-- > 0; )
                                        if (fixlst[i][0] == px &&
                                                        fixlst[i][1] == py)
                                                break;
                                if (i >= 0)
                                        continue;       /* already there */
                                if (nfixations % FIXHUNK == 0) {
                                        if (nfixations)
                                                fixlst = (short (*)[2])
                                                        realloc((char *)fixlst,
                                                        (nfixations+FIXHUNK)*
                                                        2*sizeof(short));
                                        else
                                                fixlst = (short (*)[2])malloc(
                                                        FIXHUNK*2*sizeof(short)
                                                        );
                                        if (fixlst == NULL)
                                                syserror("malloc");
                                }
                                fixlst[nfixations][0] = px;
                                fixlst[nfixations][1] = py;
                                nfixations++;
                        }
                }
        }
        if (!feof(fp)) {
                fprintf(stderr, "%s: format error reading fixation data\n",
                                progname);
                exit(1);
        }
#undef  FIXHUNK
}


gethisto(fp)                    /* load precomputed luminance histogram */
FILE    *fp;
{
        double  histo[MAXPREHIST];
        double  histart, histep;
        double  l, b, lastb, w;
        int     n;
        register int    i;
                                        /* load data */
        for (i = 0; i < MAXPREHIST &&
                        fscanf(fp, "%lf %lf", &b, &histo[i]) == 2; i++) {
                if (i > 1 && fabs(b - lastb - histep) > .001) {
                        fprintf(stderr,
                                "%s: uneven step size in histogram data\n",
                                        progname);
                        exit(1);
                }
                if (i == 1)
                        if ((histep = b - (histart = lastb)) <= FTINY) {
                                fprintf(stderr,
                                        "%s: illegal step in histogram data\n",
                                                progname);
                                exit(1);
                        }
                lastb = b;
        }
        if (i < 2 || !feof(fp)) {
                fprintf(stderr,
                "%s: format/length error loading histogram (log10L %f at %d)\n",
                                progname, b, i);
                exit(1);
        }
        n = i;
        histart *= LN_10;
        histep *= LN_10;
                                        /* find extrema */
        for (i = 0; i < n && histo[i] <= FTINY; i++)
                ;
        bwmin = histart + (i-.001)*histep;
        for (i = n; i-- && histo[i] <= FTINY; )
                ;
        bwmax = histart + (i+1.001)*histep;
        if (bwmax > Bl(LMAX))
                bwmax = Bl(LMAX);
        if (bwmin < Bl(LMIN))
                bwmin = Bl(LMIN);
        else                            /* duplicate bottom bin */
                bwmin = bwmax - (bwmax-bwmin)*HISTRES/(HISTRES-1);
                                        /* convert histogram */
        bwavg = 0.; histot = 0.;
        for (i = 0; i < HISTRES; i++)
                bwhist[i] = 0.;
        for (i = 0, b = histart; i < n; i++, b += histep) {
                if (b < bwmin+FTINY) continue;
                if (b >= bwmax-FTINY) break;
                w = histo[i];
                bwavg += w*b;
                bwhist[bwhi(b)] += w;
                histot += w;
        }
        bwavg /= histot;
        if (bwmin > Bl(LMIN)+FTINY) {   /* add false samples at bottom */
                bwhist[1] *= 0.5;
                bwhist[0] += bwhist[1];
        }
}


double
centprob(x, y)                  /* center-weighting probability function */
int     x, y;
{
        double  xr, yr, p;
                                /* paraboloid, 0 at 90 degrees from center */
        xr = (x - .5*(fvxr-1))/90.;     /* 180 degree fisheye has fv?r == 90 */
        yr = (y - .5*(fvyr-1))/90.;
        p = 1. - xr*xr - yr*yr;
        return(p < 0. ? 0. : p);
}


comphist()                      /* create foveal sampling histogram */
{
        double  l, b, w, lwmin, lwmax;
        register int    x, y;
                                        /* check for precalculated histogram */
        if (what2do&DO_PREHIST)
                return;
        lwmin = 1e10;                   /* find extrema */
        lwmax = 0.;
        for (y = 0; y < fvyr; y++)
                for (x = 0; x < fvxr; x++) {
                        l = plum(fovscan(y)[x]);
                        if (l < lwmin) lwmin = l;
                        if (l > lwmax) lwmax = l;
                }
        lwmax *= 1.01;
        if (lwmax > LMAX)
                lwmax = LMAX;
        bwmax = Bl(lwmax);
        if (lwmin < LMIN) {
                lwmin = LMIN;
                bwmin = Bl(LMIN);
        } else {                        /* duplicate bottom bin */
                bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
                lwmin = Lb(bwmin);
        }
                                        /* (re)compute histogram */
        bwavg = 0.;
        histot = 0.;
        for (x = 0; x < HISTRES; x++)
                bwhist[x] = 0.;
                                        /* global average */
        if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
                for (y = 0; y < fvyr; y++)
                        for (x = 0; x < fvxr; x++) {
                                l = plum(fovscan(y)[x]);
                                if (l < lwmin+FTINY) continue;
                                if (l >= lwmax-FTINY) continue;
                                b = Bl(l);
                                w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
                                bwavg += w*b;
                                bwhist[bwhi(b)] += w;
                                histot += w;
                        }
                                        /* average fixation points */
        if (what2do&DO_FIXHIST && nfixations > 0) {
                if (histot > FTINY)
                        w = fixfrac/(1.-fixfrac)*histot/nfixations;
                else
                        w = 1.;
                for (x = 0; x < nfixations; x++) {
                        l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
                        if (l < lwmin+FTINY) continue;
                        if (l >= lwmax-FTINY) continue;
                        b = Bl(l);
                        bwavg += w*b;
                        bwhist[bwhi(b)] += w;
                        histot += w;
                }
        }
        bwavg /= histot;
        if (lwmin > LMIN+FTINY) {       /* add false samples at bottom */
                bwhist[1] *= 0.5;
                bwhist[0] += bwhist[1];
        }
}


mkcumf()                        /* make cumulative distribution function */
{
        register int    i;
        register double sum;

        mhistot = 0.;           /* compute modified total */
        for (i = 0; i < HISTRES; i++)
                mhistot += modhist[i];

        sum = 0.;               /* compute cumulative function */
        for (i = 0; i < HISTRES; i++) {
                cumf[i] = sum/mhistot;
                sum += modhist[i];
        }
        cumf[HISTRES] = 1.;
}


double
cf(b)                           /* return cumulative function at b */
double  b;
{
        double  x;
        register int    i;

        i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
        x -= (double)i;
        return(cumf[i]*(1.-x) + cumf[i+1]*x);
}


double
BLw(Lw)                         /* map world luminance to display brightness */
double  Lw;
{
        double  b;

        if (Lw <= LMIN || (b = Bl(Lw)) <= bwmin+FTINY)
                return(Bldmin);
        if (b >= bwmax-FTINY)
                return(Bldmax);
        return(Bldmin + cf(b)*(Bldmax-Bldmin));
}


double
htcontrs(La)            /* human threshold contrast sensitivity, dL(La) */
double  La;
{
        double  l10La, l10dL;
                                /* formula taken from Ferwerda et al. [SG96] */
        if (La < 1.148e-4)
                return(1.38e-3);
        l10La = log10(La);
        if (l10La < -1.44)              /* rod response regime */
                l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
        else if (l10La < -.0184)
                l10dL = l10La - .395;
        else if (l10La < 1.9)           /* cone response regime */
                l10dL = pow(.249*l10La + .65, 2.7) - .72;
        else
                l10dL = l10La - 1.255;

        return(exp10(l10dL));
}


double
clampf(Lw)                      /* histogram clamping function */
double  Lw;
{
        double  bLw, ratio;

        bLw = BLw(Lw);          /* apply current brightness map */
        ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
        return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
}

double
crfactor(Lw)                    /* contrast reduction factor */
double  Lw;
{
        int     i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
        double  bLw, ratio, Tdb;

        if (i <= 0)
                return(1.0);
        if (i >= HISTRES)
                return(1.0);
        bLw = BLw(Lw);
        ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
        Tdb = mhistot * (bwmax - bwmin) / HISTRES;
        return(modhist[i]*Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)/(Tdb*ratio));
}


#if ADJ_VEIL
mkcrfimage()                    /* compute contrast reduction factor image */
{
        int     i;
        float   *crfptr;
        COLOR   *fovptr;

        if (crfimg == NULL)
                crfimg = (float *)malloc(fvxr*fvyr*sizeof(float));
        if (crfimg == NULL)
                syserror("malloc");
        crfptr = crfimg;
        fovptr = fovimg;
        for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
                crfptr[0] = crfactor(plum(fovptr[0]));
}
#endif


int
mkbrmap()                       /* make dynamic range map */
{
        double  Tdb, b, s;
        double  ceiling, trimmings;
        register int    i;
                                        /* copy initial histogram */
        bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
        s = (bwmax - bwmin)/HISTRES;    /* s is delta b */
                                        /* loop until satisfactory */
        do {
                mkcumf();                       /* sync brightness mapping */
                if (mhistot <= histot*CVRATIO)
                        return(-1);             /* no compression needed! */
                Tdb = mhistot * s;
                trimmings = 0.;                 /* clip to envelope */
                for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
                        ceiling = Tdb*clampf(Lb(b));
                        if (modhist[i] > ceiling) {
                                trimmings += modhist[i] - ceiling;
                                modhist[i] = ceiling;
                        }
                }
        } while (trimmings > histot*CVRATIO);

#if ADJ_VEIL
        mkcrfimage();                   /* contrast reduction image */
#endif

        return(0);                      /* we got it */
}


scotscan(scan, xres)            /* apply scotopic color sensitivity loss */
COLOR   *scan;
int     xres;
{
        COLOR   ctmp;
        double  incolor, b, Lw;
        register int    i;

        for (i = 0; i < xres; i++) {
                Lw = plum(scan[i]);
                if (Lw >= TopMesopic)
                        incolor = 1.;
                else if (Lw <= BotMesopic)
                        incolor = 0.;
                else
                        incolor = (Lw - BotMesopic) /
                                        (TopMesopic - BotMesopic);
                if (incolor < 1.-FTINY) {
                        b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
                        if (lumf == rgblum) b /= WHTEFFICACY;
                        setcolor(ctmp, b, b, b);
                        if (incolor <= FTINY)
                                setcolor(scan[i], 0., 0., 0.);
                        else
                                scalecolor(scan[i], incolor);
                        addcolor(scan[i], ctmp);
                }
        }
}


mapscan(scan, xres)             /* apply tone mapping operator to scanline */
COLOR   *scan;
int     xres;
{
        double  mult, Lw, b;
        register int    x;

        for (x = 0; x < xres; x++) {
                Lw = plum(scan[x]);
                if (Lw < LMIN) {
                        setcolor(scan[x], 0., 0., 0.);
                        continue;
                }
                b = BLw(Lw);            /* apply brightness mapping */
                mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
                if (lumf == rgblum) mult *= WHTEFFICACY;
                scalecolor(scan[x], mult);
        }
}


putmapping(fp)                  /* put out mapping function */
FILE    *fp;
{
        double  b, s;
        register int    i;
        double  wlum, sf, dlum;

        sf = scalef*inpexp;
        if (lumf == cielum) sf *= WHTEFFICACY;
        s = (bwmax - bwmin)/HISTRES;
        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
                wlum = Lb(b);
                if (what2do&DO_LINEAR) {
                        dlum = sf*wlum;
                        if (dlum > ldmax) dlum = ldmax;
                        else if (dlum < ldmin) dlum = ldmin;
                        fprintf(fp, "%e %e\n", wlum, dlum);
                } else
                        fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
        }
}
Revision:	3.11
Committed:	Sat Feb 22 02:07:27 2003 UTC (22 years, 8 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R5
Changes since 3.10:	+73 -31 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	User	Rev	Content
1	greg	3.1	#ifndef lint
2	greg	3.11	static const char RCSid[] = "$Id$";
3	greg	3.1	#endif
4			/*
5			* Routines for computing and applying brightness mapping.
6			*/
7
8			#include "pcond.h"
9
10
11	greg	3.5	#define CVRATIO 0.025 /* fraction of samples allowed > env. */
12	greg	3.1
13	gwlarson	3.10	#define LN_10 2.30258509299404568402
14			#define exp10(x) exp(LN_10*(x))
15	greg	3.1
16	greg	3.11	double modhist[HISTRES]; /* modified histogram */
17	greg	3.5	double mhistot; /* modified histogram total */
18	greg	3.11	double cumf[HISTRES+1]; /* cumulative distribution function */
19	greg	3.1
20
21	greg	3.8	getfixations(fp) /* load fixation history list */
22			FILE *fp;
23			{
24			#define FIXHUNK 128
25			RESOLU fvres;
26			int pos[2];
27			register int px, py, i;
28			/* initialize our resolution struct */
29	greg	3.11	if ((fvres.rt=inpres.rt)&YMAJOR) {
30	greg	3.8	fvres.xr = fvxr;
31			fvres.yr = fvyr;
32			} else {
33			fvres.xr = fvyr;
34			fvres.yr = fvxr;
35			}
36			/* read each picture position */
37			while (fscanf(fp, "%d %d", &pos[0], &pos[1]) == 2) {
38			/* convert to closest index in foveal image */
39			loc2pix(pos, &fvres,
40			(pos[0]+.5)/inpres.xr, (pos[1]+.5)/inpres.yr);
41			/* include nine neighborhood samples */
42			for (px = pos[0]-1; px <= pos[0]+1; px++) {
43			if (px < 0 \|\| px >= fvxr)
44			continue;
45			for (py = pos[1]-1; py <= pos[1]+1; py++) {
46			if (py < 0 \|\| py >= fvyr)
47			continue;
48			for (i = nfixations; i-- > 0; )
49			if (fixlst[i][0] == px &&
50			fixlst[i][1] == py)
51			break;
52			if (i >= 0)
53			continue; /* already there */
54			if (nfixations % FIXHUNK == 0) {
55			if (nfixations)
56			fixlst = (short (*)[2])
57			realloc((char *)fixlst,
58			(nfixations+FIXHUNK)*
59			2*sizeof(short));
60			else
61			fixlst = (short (*)[2])malloc(
62			FIXHUNK2sizeof(short)
63			);
64			if (fixlst == NULL)
65			syserror("malloc");
66			}
67			fixlst[nfixations][0] = px;
68			fixlst[nfixations][1] = py;
69			nfixations++;
70			}
71			}
72			}
73			if (!feof(fp)) {
74			fprintf(stderr, "%s: format error reading fixation data\n",
75			progname);
76			exit(1);
77			}
78			#undef FIXHUNK
79			}
80
81
82	gwlarson	3.10	gethisto(fp) /* load precomputed luminance histogram */
83			FILE *fp;
84			{
85	greg	3.11	double histo[MAXPREHIST];
86	gwlarson	3.10	double histart, histep;
87			double l, b, lastb, w;
88			int n;
89			register int i;
90			/* load data */
91			for (i = 0; i < MAXPREHIST &&
92	greg	3.11	fscanf(fp, "%lf %lf", &b, &histo[i]) == 2; i++) {
93			if (i > 1 && fabs(b - lastb - histep) > .001) {
94	gwlarson	3.10	fprintf(stderr,
95			"%s: uneven step size in histogram data\n",
96			progname);
97			exit(1);
98			}
99			if (i == 1)
100	greg	3.11	if ((histep = b - (histart = lastb)) <= FTINY) {
101			fprintf(stderr,
102			"%s: illegal step in histogram data\n",
103			progname);
104			exit(1);
105			}
106	gwlarson	3.10	lastb = b;
107			}
108			if (i < 2 \|\| !feof(fp)) {
109			fprintf(stderr,
110			"%s: format/length error loading histogram (log10L %f at %d)\n",
111			progname, b, i);
112			exit(1);
113			}
114			n = i;
115			histart *= LN_10;
116			histep *= LN_10;
117			/* find extrema */
118			for (i = 0; i < n && histo[i] <= FTINY; i++)
119			;
120	greg	3.11	bwmin = histart + (i-.001)*histep;
121	gwlarson	3.10	for (i = n; i-- && histo[i] <= FTINY; )
122			;
123	greg	3.11	bwmax = histart + (i+1.001)*histep;
124	gwlarson	3.10	if (bwmax > Bl(LMAX))
125			bwmax = Bl(LMAX);
126			if (bwmin < Bl(LMIN))
127			bwmin = Bl(LMIN);
128			else /* duplicate bottom bin */
129			bwmin = bwmax - (bwmax-bwmin)*HISTRES/(HISTRES-1);
130			/* convert histogram */
131			bwavg = 0.; histot = 0.;
132			for (i = 0; i < HISTRES; i++)
133			bwhist[i] = 0.;
134			for (i = 0, b = histart; i < n; i++, b += histep) {
135	greg	3.11	if (b < bwmin+FTINY) continue;
136			if (b >= bwmax-FTINY) break;
137	gwlarson	3.10	w = histo[i];
138			bwavg += w*b;
139			bwhist[bwhi(b)] += w;
140			histot += w;
141			}
142			bwavg /= histot;
143			if (bwmin > Bl(LMIN)+FTINY) { /* add false samples at bottom */
144			bwhist[1] *= 0.5;
145			bwhist[0] += bwhist[1];
146			}
147			}
148
149
150	greg	3.8	double
151			centprob(x, y) /* center-weighting probability function */
152			int x, y;
153			{
154			double xr, yr, p;
155			/* paraboloid, 0 at 90 degrees from center */
156			xr = (x - .5(fvxr-1))/90.; / 180 degree fisheye has fv?r == 90 */
157			yr = (y - .5*(fvyr-1))/90.;
158			p = 1. - xrxr - yryr;
159			return(p < 0. ? 0. : p);
160			}
161
162
163			comphist() /* create foveal sampling histogram */
164			{
165			double l, b, w, lwmin, lwmax;
166			register int x, y;
167	gwlarson	3.10	/* check for precalculated histogram */
168			if (what2do&DO_PREHIST)
169			return;
170	greg	3.8	lwmin = 1e10; /* find extrema */
171			lwmax = 0.;
172			for (y = 0; y < fvyr; y++)
173			for (x = 0; x < fvxr; x++) {
174			l = plum(fovscan(y)[x]);
175			if (l < lwmin) lwmin = l;
176			if (l > lwmax) lwmax = l;
177			}
178	greg	3.9	lwmax *= 1.01;
179			if (lwmax > LMAX)
180			lwmax = LMAX;
181	greg	3.8	bwmax = Bl(lwmax);
182	greg	3.9	if (lwmin < LMIN) {
183			lwmin = LMIN;
184			bwmin = Bl(LMIN);
185			} else { /* duplicate bottom bin */
186			bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
187			lwmin = Lb(bwmin);
188			}
189	greg	3.8	/* (re)compute histogram */
190			bwavg = 0.;
191			histot = 0.;
192			for (x = 0; x < HISTRES; x++)
193			bwhist[x] = 0.;
194			/* global average */
195			if (!(what2do&DO_FIXHIST) \|\| fixfrac < 1.-FTINY)
196			for (y = 0; y < fvyr; y++)
197			for (x = 0; x < fvxr; x++) {
198			l = plum(fovscan(y)[x]);
199	greg	3.11	if (l < lwmin+FTINY) continue;
200			if (l >= lwmax-FTINY) continue;
201	greg	3.8	b = Bl(l);
202			w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
203	gwlarson	3.10	bwavg += w*b;
204	greg	3.8	bwhist[bwhi(b)] += w;
205			histot += w;
206			}
207			/* average fixation points */
208			if (what2do&DO_FIXHIST && nfixations > 0) {
209			if (histot > FTINY)
210			w = fixfrac/(1.-fixfrac)*histot/nfixations;
211			else
212			w = 1.;
213			for (x = 0; x < nfixations; x++) {
214			l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
215	greg	3.11	if (l < lwmin+FTINY) continue;
216			if (l >= lwmax-FTINY) continue;
217	greg	3.8	b = Bl(l);
218	gwlarson	3.10	bwavg += w*b;
219	greg	3.8	bwhist[bwhi(b)] += w;
220			histot += w;
221			}
222			}
223			bwavg /= histot;
224	greg	3.9	if (lwmin > LMIN+FTINY) { /* add false samples at bottom */
225			bwhist[1] *= 0.5;
226			bwhist[0] += bwhist[1];
227			}
228	greg	3.8	}
229
230
231	greg	3.1	mkcumf() /* make cumulative distribution function */
232			{
233			register int i;
234	greg	3.5	register double sum;
235	greg	3.1
236	greg	3.5	mhistot = 0.; /* compute modified total */
237			for (i = 0; i < HISTRES; i++)
238			mhistot += modhist[i];
239
240			sum = 0.; /* compute cumulative function */
241			for (i = 0; i < HISTRES; i++) {
242			cumf[i] = sum/mhistot;
243	greg	3.1	sum += modhist[i];
244			}
245			cumf[HISTRES] = 1.;
246			}
247
248
249			double
250			cf(b) /* return cumulative function at b */
251			double b;
252			{
253			double x;
254			register int i;
255
256			i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
257			x -= (double)i;
258			return(cumf[i](1.-x) + cumf[i+1]x);
259			}
260
261
262			double
263			BLw(Lw) /* map world luminance to display brightness */
264			double Lw;
265			{
266			double b;
267
268			if (Lw <= LMIN \|\| (b = Bl(Lw)) <= bwmin+FTINY)
269			return(Bldmin);
270			if (b >= bwmax-FTINY)
271			return(Bldmax);
272	greg	3.4	return(Bldmin + cf(b)*(Bldmax-Bldmin));
273	greg	3.1	}
274
275
276			double
277			htcontrs(La) /* human threshold contrast sensitivity, dL(La) */
278			double La;
279			{
280			double l10La, l10dL;
281			/* formula taken from Ferwerda et al. [SG96] */
282			if (La < 1.148e-4)
283			return(1.38e-3);
284			l10La = log10(La);
285			if (l10La < -1.44) /* rod response regime */
286			l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
287			else if (l10La < -.0184)
288			l10dL = l10La - .395;
289			else if (l10La < 1.9) /* cone response regime */
290			l10dL = pow(.249*l10La + .65, 2.7) - .72;
291			else
292			l10dL = l10La - 1.255;
293
294			return(exp10(l10dL));
295			}
296
297
298			double
299	greg	3.11	clampf(Lw) /* histogram clamping function */
300	greg	3.1	double Lw;
301			{
302			double bLw, ratio;
303
304			bLw = BLw(Lw); /* apply current brightness map */
305			ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
306			return(ratio/(Lb1(bLw)(Bldmax-Bldmin)Bl1(Lw)));
307			}
308
309	greg	3.11	double
310			crfactor(Lw) /* contrast reduction factor */
311			double Lw;
312			{
313			int i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
314			double bLw, ratio, Tdb;
315
316			if (i <= 0)
317			return(1.0);
318			if (i >= HISTRES)
319			return(1.0);
320			bLw = BLw(Lw);
321			ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
322			Tdb = mhistot * (bwmax - bwmin) / HISTRES;
323			return(modhist[i]Lb1(bLw)(Bldmax-Bldmin)Bl1(Lw)/(Tdbratio));
324			}
325
326
327			#if ADJ_VEIL
328			mkcrfimage() /* compute contrast reduction factor image */
329			{
330			int i;
331			float *crfptr;
332			COLOR *fovptr;
333
334			if (crfimg == NULL)
335			crfimg = (float )malloc(fvxrfvyr*sizeof(float));
336			if (crfimg == NULL)
337			syserror("malloc");
338			crfptr = crfimg;
339			fovptr = fovimg;
340			for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
341			crfptr[0] = crfactor(plum(fovptr[0]));
342			}
343			#endif
344
345	greg	3.1
346			int
347			mkbrmap() /* make dynamic range map */
348			{
349	greg	3.11	double Tdb, b, s;
350	greg	3.5	double ceiling, trimmings;
351	greg	3.1	register int i;
352			/* copy initial histogram */
353	greg	3.5	bcopy((char )bwhist, (char )modhist, sizeof(modhist));
354	greg	3.11	s = (bwmax - bwmin)/HISTRES; /* s is delta b */
355	greg	3.1	/* loop until satisfactory */
356	greg	3.5	do {
357			mkcumf(); /* sync brightness mapping */
358			if (mhistot <= histot*CVRATIO)
359			return(-1); /* no compression needed! */
360	greg	3.11	Tdb = mhistot * s;
361	greg	3.5	trimmings = 0.; /* clip to envelope */
362	greg	3.1	for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
363	greg	3.11	ceiling = Tdb*clampf(Lb(b));
364	greg	3.5	if (modhist[i] > ceiling) {
365			trimmings += modhist[i] - ceiling;
366			modhist[i] = ceiling;
367			}
368	greg	3.1	}
369	greg	3.5	} while (trimmings > histot*CVRATIO);
370
371	greg	3.11	#if ADJ_VEIL
372			mkcrfimage(); /* contrast reduction image */
373			#endif
374
375	greg	3.5	return(0); /* we got it */
376	greg	3.1	}
377
378
379			scotscan(scan, xres) /* apply scotopic color sensitivity loss */
380			COLOR *scan;
381			int xres;
382			{
383			COLOR ctmp;
384			double incolor, b, Lw;
385			register int i;
386
387			for (i = 0; i < xres; i++) {
388			Lw = plum(scan[i]);
389			if (Lw >= TopMesopic)
390			incolor = 1.;
391			else if (Lw <= BotMesopic)
392			incolor = 0.;
393			else
394			incolor = (Lw - BotMesopic) /
395			(TopMesopic - BotMesopic);
396			if (incolor < 1.-FTINY) {
397	greg	3.2	b = (1.-incolor)slum(scan[i])inpexp/SWNORM;
398			if (lumf == rgblum) b /= WHTEFFICACY;
399			setcolor(ctmp, b, b, b);
400	greg	3.1	if (incolor <= FTINY)
401			setcolor(scan[i], 0., 0., 0.);
402			else
403			scalecolor(scan[i], incolor);
404			addcolor(scan[i], ctmp);
405			}
406			}
407			}
408
409
410			mapscan(scan, xres) /* apply tone mapping operator to scanline */
411			COLOR *scan;
412			int xres;
413			{
414			double mult, Lw, b;
415	greg	3.11	register int x;
416	greg	3.1
417	greg	3.11	for (x = 0; x < xres; x++) {
418			Lw = plum(scan[x]);
419	greg	3.1	if (Lw < LMIN) {
420	greg	3.11	setcolor(scan[x], 0., 0., 0.);
421	greg	3.1	continue;
422			}
423	greg	3.11	b = BLw(Lw); /* apply brightness mapping */
424			mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
425	greg	3.1	if (lumf == rgblum) mult *= WHTEFFICACY;
426	greg	3.11	scalecolor(scan[x], mult);
427	greg	3.4	}
428			}
429
430
431			putmapping(fp) /* put out mapping function */
432			FILE *fp;
433			{
434			double b, s;
435			register int i;
436	greg	3.7	double wlum, sf, dlum;
437	greg	3.4
438			sf = scalef*inpexp;
439			if (lumf == cielum) sf *= WHTEFFICACY;
440			s = (bwmax - bwmin)/HISTRES;
441			for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
442			wlum = Lb(b);
443	greg	3.7	if (what2do&DO_LINEAR) {
444			dlum = sf*wlum;
445			if (dlum > ldmax) dlum = ldmax;
446			else if (dlum < ldmin) dlum = ldmin;
447			fprintf(fp, "%e %e\n", wlum, dlum);
448			} else
449	greg	3.4	fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
450	greg	3.1	}
451			}