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 */


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.7
Committed:	Sat Jan 11 10:39:25 1997 UTC (28 years, 10 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.6:	+8 -5 lines
Log Message:	fixed writing of linear ramp with -m option
#	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			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			mkbrmap() /* make dynamic range map */
104			{
105			double T, b, s;
106	greg	3.5	double ceiling, trimmings;
107	greg	3.1	register int i;
108			/* copy initial histogram */
109	greg	3.5	bcopy((char )bwhist, (char )modhist, sizeof(modhist));
110	greg	3.1	s = (bwmax - bwmin)/HISTRES;
111			/* loop until satisfactory */
112	greg	3.5	do {
113			mkcumf(); /* sync brightness mapping */
114			if (mhistot <= histot*CVRATIO)
115			return(-1); /* no compression needed! */
116			T = mhistot * (bwmax - bwmin) / HISTRES;
117			trimmings = 0.; /* clip to envelope */
118	greg	3.1	for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
119	greg	3.5	ceiling = T*clampf(Lb(b));
120			if (modhist[i] > ceiling) {
121			trimmings += modhist[i] - ceiling;
122			modhist[i] = ceiling;
123			}
124	greg	3.1	}
125	greg	3.5	} while (trimmings > histot*CVRATIO);
126
127			return(0); /* we got it */
128	greg	3.1	}
129
130
131			scotscan(scan, xres) /* apply scotopic color sensitivity loss */
132			COLOR *scan;
133			int xres;
134			{
135			COLOR ctmp;
136			double incolor, b, Lw;
137			register int i;
138
139			for (i = 0; i < xres; i++) {
140			Lw = plum(scan[i]);
141			if (Lw >= TopMesopic)
142			incolor = 1.;
143			else if (Lw <= BotMesopic)
144			incolor = 0.;
145			else
146			incolor = (Lw - BotMesopic) /
147			(TopMesopic - BotMesopic);
148			if (incolor < 1.-FTINY) {
149	greg	3.2	b = (1.-incolor)slum(scan[i])inpexp/SWNORM;
150			if (lumf == rgblum) b /= WHTEFFICACY;
151			setcolor(ctmp, b, b, b);
152	greg	3.1	if (incolor <= FTINY)
153			setcolor(scan[i], 0., 0., 0.);
154			else
155			scalecolor(scan[i], incolor);
156			addcolor(scan[i], ctmp);
157			}
158			}
159			}
160
161
162			mapscan(scan, xres) /* apply tone mapping operator to scanline */
163			COLOR *scan;
164			int xres;
165			{
166			double mult, Lw, b;
167			register int i;
168
169			for (i = 0; i < xres; i++) {
170			Lw = plum(scan[i]);
171			if (Lw < LMIN) {
172			setcolor(scan[i], 0., 0., 0.);
173			continue;
174			}
175			b = BLw(Lw);
176			mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
177			if (lumf == rgblum) mult *= WHTEFFICACY;
178			scalecolor(scan[i], mult);
179	greg	3.4	}
180			}
181
182
183			putmapping(fp) /* put out mapping function */
184			FILE *fp;
185			{
186			double b, s;
187			register int i;
188	greg	3.7	double wlum, sf, dlum;
189	greg	3.4
190			sf = scalef*inpexp;
191			if (lumf == cielum) sf *= WHTEFFICACY;
192			s = (bwmax - bwmin)/HISTRES;
193			for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
194			wlum = Lb(b);
195	greg	3.7	if (what2do&DO_LINEAR) {
196			dlum = sf*wlum;
197			if (dlum > ldmax) dlum = ldmax;
198			else if (dlum < ldmin) dlum = ldmin;
199			fprintf(fp, "%e %e\n", wlum, dlum);
200			} else
201	greg	3.4	fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
202	greg	3.1	}
203			}