src/px/pcond3.c

/* Copyright (c) 1997 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 * Routines for computing and applying brightness mapping.
 */

#include "pcond.h"


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

#define exp10(x)        exp(2.302585093*(x))

float   modhist[HISTRES];               /* modified histogram */
double  mhistot;                        /* modified histogram total */
float   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.or=inpres.or)&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
}


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;

        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) continue;
                                if (l > lwmax) continue;
                                b = Bl(l);
                                bwavg += b;
                                w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
                                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) continue;
                        if (l > lwmax) continue;
                        b = Bl(l);
                        bwavg += 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)              /* derivative 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)));
}


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

        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    i;

        for (i = 0; i < xres; i++) {
                Lw = plum(scan[i]);
                if (Lw < LMIN) {
                        setcolor(scan[i], 0., 0., 0.);
                        continue;
                }
                b = BLw(Lw);
                mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
                if (lumf == rgblum) mult *= WHTEFFICACY;
                scalecolor(scan[i], 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.9
Committed:	Fri Apr 11 18:34:39 1997 UTC (28 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.8:	+14 -5 lines
Log Message:	added an extra histogram bin at lower limit to avoid clipping
#	User	Rev	Content
1	greg	3.7	/* Copyright (c) 1997 Regents of the University of California */
2	greg	3.1
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* Routines for computing and applying brightness mapping.
9			*/
10
11			#include "pcond.h"
12
13
14	greg	3.5	#define CVRATIO 0.025 /* fraction of samples allowed > env. */
15	greg	3.1
16			#define exp10(x) exp(2.302585093*(x))
17
18	greg	3.5	float modhist[HISTRES]; /* modified histogram */
19			double mhistot; /* modified histogram total */
20	greg	3.1	float cumf[HISTRES+1]; /* cumulative distribution function */
21
22
23	greg	3.8	getfixations(fp) /* load fixation history list */
24			FILE *fp;
25			{
26			#define FIXHUNK 128
27			RESOLU fvres;
28			int pos[2];
29			register int px, py, i;
30			/* initialize our resolution struct */
31			if ((fvres.or=inpres.or)&YMAJOR) {
32			fvres.xr = fvxr;
33			fvres.yr = fvyr;
34			} else {
35			fvres.xr = fvyr;
36			fvres.yr = fvxr;
37			}
38			/* read each picture position */
39			while (fscanf(fp, "%d %d", &pos[0], &pos[1]) == 2) {
40			/* convert to closest index in foveal image */
41			loc2pix(pos, &fvres,
42			(pos[0]+.5)/inpres.xr, (pos[1]+.5)/inpres.yr);
43			/* include nine neighborhood samples */
44			for (px = pos[0]-1; px <= pos[0]+1; px++) {
45			if (px < 0 \|\| px >= fvxr)
46			continue;
47			for (py = pos[1]-1; py <= pos[1]+1; py++) {
48			if (py < 0 \|\| py >= fvyr)
49			continue;
50			for (i = nfixations; i-- > 0; )
51			if (fixlst[i][0] == px &&
52			fixlst[i][1] == py)
53			break;
54			if (i >= 0)
55			continue; /* already there */
56			if (nfixations % FIXHUNK == 0) {
57			if (nfixations)
58			fixlst = (short (*)[2])
59			realloc((char *)fixlst,
60			(nfixations+FIXHUNK)*
61			2*sizeof(short));
62			else
63			fixlst = (short (*)[2])malloc(
64			FIXHUNK2sizeof(short)
65			);
66			if (fixlst == NULL)
67			syserror("malloc");
68			}
69			fixlst[nfixations][0] = px;
70			fixlst[nfixations][1] = py;
71			nfixations++;
72			}
73			}
74			}
75			if (!feof(fp)) {
76			fprintf(stderr, "%s: format error reading fixation data\n",
77			progname);
78			exit(1);
79			}
80			#undef FIXHUNK
81			}
82
83
84			double
85			centprob(x, y) /* center-weighting probability function */
86			int x, y;
87			{
88			double xr, yr, p;
89			/* paraboloid, 0 at 90 degrees from center */
90			xr = (x - .5(fvxr-1))/90.; / 180 degree fisheye has fv?r == 90 */
91			yr = (y - .5*(fvyr-1))/90.;
92			p = 1. - xrxr - yryr;
93			return(p < 0. ? 0. : p);
94			}
95
96
97			comphist() /* create foveal sampling histogram */
98			{
99			double l, b, w, lwmin, lwmax;
100			register int x, y;
101
102			lwmin = 1e10; /* find extrema */
103			lwmax = 0.;
104			for (y = 0; y < fvyr; y++)
105			for (x = 0; x < fvxr; x++) {
106			l = plum(fovscan(y)[x]);
107			if (l < lwmin) lwmin = l;
108			if (l > lwmax) lwmax = l;
109			}
110	greg	3.9	lwmax *= 1.01;
111			if (lwmax > LMAX)
112			lwmax = LMAX;
113	greg	3.8	bwmax = Bl(lwmax);
114	greg	3.9	if (lwmin < LMIN) {
115			lwmin = LMIN;
116			bwmin = Bl(LMIN);
117			} else { /* duplicate bottom bin */
118			bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
119			lwmin = Lb(bwmin);
120			}
121	greg	3.8	/* (re)compute histogram */
122			bwavg = 0.;
123			histot = 0.;
124			for (x = 0; x < HISTRES; x++)
125			bwhist[x] = 0.;
126			/* global average */
127			if (!(what2do&DO_FIXHIST) \|\| fixfrac < 1.-FTINY)
128			for (y = 0; y < fvyr; y++)
129			for (x = 0; x < fvxr; x++) {
130			l = plum(fovscan(y)[x]);
131			if (l < lwmin) continue;
132			if (l > lwmax) continue;
133			b = Bl(l);
134			bwavg += b;
135			w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
136			bwhist[bwhi(b)] += w;
137			histot += w;
138			}
139			/* average fixation points */
140			if (what2do&DO_FIXHIST && nfixations > 0) {
141			if (histot > FTINY)
142			w = fixfrac/(1.-fixfrac)*histot/nfixations;
143			else
144			w = 1.;
145			for (x = 0; x < nfixations; x++) {
146			l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
147			if (l < lwmin) continue;
148			if (l > lwmax) continue;
149			b = Bl(l);
150			bwavg += b;
151			bwhist[bwhi(b)] += w;
152			histot += w;
153			}
154			}
155			bwavg /= histot;
156	greg	3.9	if (lwmin > LMIN+FTINY) { /* add false samples at bottom */
157			bwhist[1] *= 0.5;
158			bwhist[0] += bwhist[1];
159			}
160	greg	3.8	}
161
162
163	greg	3.1	mkcumf() /* make cumulative distribution function */
164			{
165			register int i;
166	greg	3.5	register double sum;
167	greg	3.1
168	greg	3.5	mhistot = 0.; /* compute modified total */
169			for (i = 0; i < HISTRES; i++)
170			mhistot += modhist[i];
171
172			sum = 0.; /* compute cumulative function */
173			for (i = 0; i < HISTRES; i++) {
174			cumf[i] = sum/mhistot;
175	greg	3.1	sum += modhist[i];
176			}
177			cumf[HISTRES] = 1.;
178			}
179
180
181			double
182			cf(b) /* return cumulative function at b */
183			double b;
184			{
185			double x;
186			register int i;
187
188			i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
189			x -= (double)i;
190			return(cumf[i](1.-x) + cumf[i+1]x);
191			}
192
193
194			double
195			BLw(Lw) /* map world luminance to display brightness */
196			double Lw;
197			{
198			double b;
199
200			if (Lw <= LMIN \|\| (b = Bl(Lw)) <= bwmin+FTINY)
201			return(Bldmin);
202			if (b >= bwmax-FTINY)
203			return(Bldmax);
204	greg	3.4	return(Bldmin + cf(b)*(Bldmax-Bldmin));
205	greg	3.1	}
206
207
208			double
209			htcontrs(La) /* human threshold contrast sensitivity, dL(La) */
210			double La;
211			{
212			double l10La, l10dL;
213			/* formula taken from Ferwerda et al. [SG96] */
214			if (La < 1.148e-4)
215			return(1.38e-3);
216			l10La = log10(La);
217			if (l10La < -1.44) /* rod response regime */
218			l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
219			else if (l10La < -.0184)
220			l10dL = l10La - .395;
221			else if (l10La < 1.9) /* cone response regime */
222			l10dL = pow(.249*l10La + .65, 2.7) - .72;
223			else
224			l10dL = l10La - 1.255;
225
226			return(exp10(l10dL));
227			}
228
229
230			double
231			clampf(Lw) /* derivative clamping function */
232			double Lw;
233			{
234			double bLw, ratio;
235
236			bLw = BLw(Lw); /* apply current brightness map */
237			ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
238			return(ratio/(Lb1(bLw)(Bldmax-Bldmin)Bl1(Lw)));
239			}
240
241
242			int
243			mkbrmap() /* make dynamic range map */
244			{
245			double T, b, s;
246	greg	3.5	double ceiling, trimmings;
247	greg	3.1	register int i;
248			/* copy initial histogram */
249	greg	3.5	bcopy((char )bwhist, (char )modhist, sizeof(modhist));
250	greg	3.1	s = (bwmax - bwmin)/HISTRES;
251			/* loop until satisfactory */
252	greg	3.5	do {
253			mkcumf(); /* sync brightness mapping */
254			if (mhistot <= histot*CVRATIO)
255			return(-1); /* no compression needed! */
256			T = mhistot * (bwmax - bwmin) / HISTRES;
257			trimmings = 0.; /* clip to envelope */
258	greg	3.1	for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
259	greg	3.5	ceiling = T*clampf(Lb(b));
260			if (modhist[i] > ceiling) {
261			trimmings += modhist[i] - ceiling;
262			modhist[i] = ceiling;
263			}
264	greg	3.1	}
265	greg	3.5	} while (trimmings > histot*CVRATIO);
266
267			return(0); /* we got it */
268	greg	3.1	}
269
270
271			scotscan(scan, xres) /* apply scotopic color sensitivity loss */
272			COLOR *scan;
273			int xres;
274			{
275			COLOR ctmp;
276			double incolor, b, Lw;
277			register int i;
278
279			for (i = 0; i < xres; i++) {
280			Lw = plum(scan[i]);
281			if (Lw >= TopMesopic)
282			incolor = 1.;
283			else if (Lw <= BotMesopic)
284			incolor = 0.;
285			else
286			incolor = (Lw - BotMesopic) /
287			(TopMesopic - BotMesopic);
288			if (incolor < 1.-FTINY) {
289	greg	3.2	b = (1.-incolor)slum(scan[i])inpexp/SWNORM;
290			if (lumf == rgblum) b /= WHTEFFICACY;
291			setcolor(ctmp, b, b, b);
292	greg	3.1	if (incolor <= FTINY)
293			setcolor(scan[i], 0., 0., 0.);
294			else
295			scalecolor(scan[i], incolor);
296			addcolor(scan[i], ctmp);
297			}
298			}
299			}
300
301
302			mapscan(scan, xres) /* apply tone mapping operator to scanline */
303			COLOR *scan;
304			int xres;
305			{
306			double mult, Lw, b;
307			register int i;
308
309			for (i = 0; i < xres; i++) {
310			Lw = plum(scan[i]);
311			if (Lw < LMIN) {
312			setcolor(scan[i], 0., 0., 0.);
313			continue;
314			}
315			b = BLw(Lw);
316			mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
317			if (lumf == rgblum) mult *= WHTEFFICACY;
318			scalecolor(scan[i], mult);
319	greg	3.4	}
320			}
321
322
323			putmapping(fp) /* put out mapping function */
324			FILE *fp;
325			{
326			double b, s;
327			register int i;
328	greg	3.7	double wlum, sf, dlum;
329	greg	3.4
330			sf = scalef*inpexp;
331			if (lumf == cielum) sf *= WHTEFFICACY;
332			s = (bwmax - bwmin)/HISTRES;
333			for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
334			wlum = Lb(b);
335	greg	3.7	if (what2do&DO_LINEAR) {
336			dlum = sf*wlum;
337			if (dlum > ldmax) dlum = ldmax;
338			else if (dlum < ldmin) dlum = ldmin;
339			fprintf(fp, "%e %e\n", wlum, dlum);
340			} else
341	greg	3.4	fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
342	greg	3.1	}
343			}