src/px/pcond3.c

/* Copyright (c) 1996 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 */


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
shiftdir(bw)            /* compute shift direction for histogram */
double  bw;
{
        if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
                return(1);
        return(-1);
}


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;

        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)
                        fprintf(fp, "%e %e\n", wlum, sf*wlum);
                else
                        fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
        }
}
Revision:	3.5
Committed:	Thu Jan 9 13:56:30 1997 UTC (27 years, 3 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.4:	+27 -84 lines
Log Message:	changed to histogram truncation made veiling reflection modify FOV sample image
#	User	Rev	Content
1	greg	3.1	/* Copyright (c) 1996 Regents of the University of California */
2
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			mkcumf() /* make cumulative distribution function */
24			{
25			register int i;
26	greg	3.5	register double sum;
27	greg	3.1
28	greg	3.5	mhistot = 0.; /* compute modified total */
29			for (i = 0; i < HISTRES; i++)
30			mhistot += modhist[i];
31
32			sum = 0.; /* compute cumulative function */
33			for (i = 0; i < HISTRES; i++) {
34			cumf[i] = sum/mhistot;
35	greg	3.1	sum += modhist[i];
36			}
37			cumf[HISTRES] = 1.;
38			}
39
40
41			double
42			cf(b) /* return cumulative function at b */
43			double b;
44			{
45			double x;
46			register int i;
47
48			i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
49			x -= (double)i;
50			return(cumf[i](1.-x) + cumf[i+1]x);
51			}
52
53
54			double
55			BLw(Lw) /* map world luminance to display brightness */
56			double Lw;
57			{
58			double b;
59
60			if (Lw <= LMIN \|\| (b = Bl(Lw)) <= bwmin+FTINY)
61			return(Bldmin);
62			if (b >= bwmax-FTINY)
63			return(Bldmax);
64	greg	3.4	return(Bldmin + cf(b)*(Bldmax-Bldmin));
65	greg	3.1	}
66
67
68			double
69			htcontrs(La) /* human threshold contrast sensitivity, dL(La) */
70			double La;
71			{
72			double l10La, l10dL;
73			/* formula taken from Ferwerda et al. [SG96] */
74			if (La < 1.148e-4)
75			return(1.38e-3);
76			l10La = log10(La);
77			if (l10La < -1.44) /* rod response regime */
78			l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
79			else if (l10La < -.0184)
80			l10dL = l10La - .395;
81			else if (l10La < 1.9) /* cone response regime */
82			l10dL = pow(.249*l10La + .65, 2.7) - .72;
83			else
84			l10dL = l10La - 1.255;
85
86			return(exp10(l10dL));
87			}
88
89
90			double
91			clampf(Lw) /* derivative clamping function */
92			double Lw;
93			{
94			double bLw, ratio;
95
96			bLw = BLw(Lw); /* apply current brightness map */
97			ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
98			return(ratio/(Lb1(bLw)(Bldmax-Bldmin)Bl1(Lw)));
99			}
100
101
102			int
103			shiftdir(bw) /* compute shift direction for histogram */
104			double bw;
105			{
106			if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
107			return(1);
108			return(-1);
109			}
110
111
112			int
113			mkbrmap() /* make dynamic range map */
114			{
115			double T, b, s;
116	greg	3.5	double ceiling, trimmings;
117	greg	3.1	register int i;
118			/* copy initial histogram */
119	greg	3.5	bcopy((char )bwhist, (char )modhist, sizeof(modhist));
120	greg	3.1	s = (bwmax - bwmin)/HISTRES;
121			/* loop until satisfactory */
122	greg	3.5	do {
123			mkcumf(); /* sync brightness mapping */
124			if (mhistot <= histot*CVRATIO)
125			return(-1); /* no compression needed! */
126			T = mhistot * (bwmax - bwmin) / HISTRES;
127			trimmings = 0.; /* clip to envelope */
128	greg	3.1	for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
129	greg	3.5	ceiling = T*clampf(Lb(b));
130			if (modhist[i] > ceiling) {
131			trimmings += modhist[i] - ceiling;
132			modhist[i] = ceiling;
133			}
134	greg	3.1	}
135	greg	3.5	} while (trimmings > histot*CVRATIO);
136
137			return(0); /* we got it */
138	greg	3.1	}
139
140
141			scotscan(scan, xres) /* apply scotopic color sensitivity loss */
142			COLOR *scan;
143			int xres;
144			{
145			COLOR ctmp;
146			double incolor, b, Lw;
147			register int i;
148
149			for (i = 0; i < xres; i++) {
150			Lw = plum(scan[i]);
151			if (Lw >= TopMesopic)
152			incolor = 1.;
153			else if (Lw <= BotMesopic)
154			incolor = 0.;
155			else
156			incolor = (Lw - BotMesopic) /
157			(TopMesopic - BotMesopic);
158			if (incolor < 1.-FTINY) {
159	greg	3.2	b = (1.-incolor)slum(scan[i])inpexp/SWNORM;
160			if (lumf == rgblum) b /= WHTEFFICACY;
161			setcolor(ctmp, b, b, b);
162	greg	3.1	if (incolor <= FTINY)
163			setcolor(scan[i], 0., 0., 0.);
164			else
165			scalecolor(scan[i], incolor);
166			addcolor(scan[i], ctmp);
167			}
168			}
169			}
170
171
172			mapscan(scan, xres) /* apply tone mapping operator to scanline */
173			COLOR *scan;
174			int xres;
175			{
176			double mult, Lw, b;
177			register int i;
178
179			for (i = 0; i < xres; i++) {
180			Lw = plum(scan[i]);
181			if (Lw < LMIN) {
182			setcolor(scan[i], 0., 0., 0.);
183			continue;
184			}
185			b = BLw(Lw);
186			mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
187			if (lumf == rgblum) mult *= WHTEFFICACY;
188			scalecolor(scan[i], mult);
189	greg	3.4	}
190			}
191
192
193			putmapping(fp) /* put out mapping function */
194			FILE *fp;
195			{
196			double b, s;
197			register int i;
198			double wlum, sf;
199
200			sf = scalef*inpexp;
201			if (lumf == cielum) sf *= WHTEFFICACY;
202			s = (bwmax - bwmin)/HISTRES;
203			for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
204			wlum = Lb(b);
205			if (what2do&DO_LINEAR)
206			fprintf(fp, "%e %e\n", wlum, sf*wlum);
207			else
208			fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
209	greg	3.1	}
210			}