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 pixels allowed > env. */

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

int     modhist[HISTRES];               /* modified histogram */
float   cumf[HISTRES+1];                /* cumulative distribution function */


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

        cumf[0] = 0.;
        sum = modhist[0];
        for (i = 1; i < HISTRES; i++) {
                cumf[i] = (double)sum/histot;
                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(Bl(Lw))*(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 */
{
        int     hdiffs[HISTRES], above, below;
        double  T, b, s;
        int     maxd, maxi, sd;
        register int    i;
                                        /* copy initial histogram */
        for (i = 0; i < HISTRES; i++)
                modhist[i] = bwhist[i];
        T = histot * (bwmax - bwmin) / HISTRES;
        s = (bwmax - bwmin)/HISTRES;
                                        /* loop until satisfactory */
        for ( ; ; ) {
                mkcumf();               /* sync brightness mapping */
                above = below = 0;      /* compute visibility overflow */
                for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
                        hdiffs[i] = modhist[i] - (int)(T*clampf(Lb(b)) + .5);
                        if (hdiffs[i] > 0) above += hdiffs[i];
                        else below -= hdiffs[i];
                }
                if (above <= histot*CVRATIO)
                        break;          /* close enough */
                if (above-below >= 0)
                        return(-1);     /* Houston, we have a problem.... */
                /* original looped here as well (BEGIN_L2) */
                maxd = 0;               /* find largest overvis */
                for (i = 0; i < HISTRES; i++)
                        if (hdiffs[i] > maxd)
                                maxd = hdiffs[maxi=i];
                /* broke loop here when (maxd == 0) (BREAK_L2) */
                for (sd = shiftdir((maxi+.5)/HISTRES*(bwmax-bwmin)+bwmin);
                                hdiffs[maxi] == maxd; sd = -sd)
                        for (i = maxi+sd; i >= 0 & i < HISTRES; i += sd)
                                if (hdiffs[i] < 0) {
                                        if (hdiffs[i] <= -maxd) {
                                                modhist[i] += maxd;
                                                modhist[maxi] -= maxd;
                                                hdiffs[i] += maxd;
                                                hdiffs[maxi] = 0;
                                        } else {
                                                modhist[maxi] += hdiffs[i];
                                                modhist[i] -= hdiffs[i];
                                                hdiffs[maxi] += hdiffs[i];
                                                hdiffs[i] = 0;
                                        }
                                        break;
                                }
                /* (END_L2) */
        }
        return(0);
}


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);
        }
}


#ifdef DEBUG
doplots()                       /* generate debugging plots */
{
        double  T, b, s;
        FILE    *fp;
        char    fname[128];
        register int    i;

        T = histot * (bwmax - bwmin) / HISTRES;
        s = (bwmax - bwmin)/HISTRES;

        sprintf(fname, "%s_hist.plt", infn);
        if ((fp = fopen(fname, "w")) == NULL)
                syserror(fname);
        fputs("include=curve.plt\n", fp);
        fputs("title=\"Brightness Frequency Distribution\"\n", fp);
        fprintf(fp, "subtitle=%s\n", infn);
        fputs("ymin=0\n", fp);
        fputs("xlabel=\"Perceptual Brightness B(Lw)\"\n", fp);
        fputs("ylabel=\"Frequency Count\"\n", fp);
        fputs("Alabel=\"Histogram\"\n", fp);
        fputs("Alintype=0\n", fp);
        fputs("Blabel=\"Envelope\"\n", fp);
        fputs("Bsymsize=0\n", fp);
        fputs("Adata=\n", fp);
        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
                fprintf(fp, "\t%f %d\n", b, modhist[i]);
        fputs(";\nBdata=\n", fp);
        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
                fprintf(fp, "\t%f %f\n", b, T*clampf(Lb(b)));
        fputs(";\n", fp);
        fclose(fp);
}
#endif
Revision:	3.3
Committed:	Fri Oct 11 10:47:00 1996 UTC (27 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.2:	+0 -3 lines
Log Message:	corrected tests in check2do()
#	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			#define CVRATIO 0.025 /* fraction of pixels allowed > env. */
15
16			#define exp10(x) exp(2.302585093*(x))
17
18			int modhist[HISTRES]; /* modified histogram */
19			float cumf[HISTRES+1]; /* cumulative distribution function */
20
21
22			mkcumf() /* make cumulative distribution function */
23			{
24			register int i;
25			register long sum;
26
27			cumf[0] = 0.;
28			sum = modhist[0];
29			for (i = 1; i < HISTRES; i++) {
30			cumf[i] = (double)sum/histot;
31			sum += modhist[i];
32			}
33			cumf[HISTRES] = 1.;
34			}
35
36
37			double
38			cf(b) /* return cumulative function at b */
39			double b;
40			{
41			double x;
42			register int i;
43
44			i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
45			x -= (double)i;
46			return(cumf[i](1.-x) + cumf[i+1]x);
47			}
48
49
50			double
51			BLw(Lw) /* map world luminance to display brightness */
52			double Lw;
53			{
54			double b;
55
56			if (Lw <= LMIN \|\| (b = Bl(Lw)) <= bwmin+FTINY)
57			return(Bldmin);
58			if (b >= bwmax-FTINY)
59			return(Bldmax);
60			return(Bldmin + cf(Bl(Lw))*(Bldmax-Bldmin));
61			}
62
63
64			double
65			htcontrs(La) /* human threshold contrast sensitivity, dL(La) */
66			double La;
67			{
68			double l10La, l10dL;
69			/* formula taken from Ferwerda et al. [SG96] */
70			if (La < 1.148e-4)
71			return(1.38e-3);
72			l10La = log10(La);
73			if (l10La < -1.44) /* rod response regime */
74			l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
75			else if (l10La < -.0184)
76			l10dL = l10La - .395;
77			else if (l10La < 1.9) /* cone response regime */
78			l10dL = pow(.249*l10La + .65, 2.7) - .72;
79			else
80			l10dL = l10La - 1.255;
81
82			return(exp10(l10dL));
83			}
84
85
86			double
87			clampf(Lw) /* derivative clamping function */
88			double Lw;
89			{
90			double bLw, ratio;
91
92			bLw = BLw(Lw); /* apply current brightness map */
93			ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
94			return(ratio/(Lb1(bLw)(Bldmax-Bldmin)Bl1(Lw)));
95			}
96
97
98			int
99			shiftdir(bw) /* compute shift direction for histogram */
100			double bw;
101			{
102			if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
103			return(1);
104			return(-1);
105			}
106
107
108			int
109			mkbrmap() /* make dynamic range map */
110			{
111			int hdiffs[HISTRES], above, below;
112			double T, b, s;
113			int maxd, maxi, sd;
114			register int i;
115			/* copy initial histogram */
116			for (i = 0; i < HISTRES; i++)
117			modhist[i] = bwhist[i];
118			T = histot * (bwmax - bwmin) / HISTRES;
119			s = (bwmax - bwmin)/HISTRES;
120			/* loop until satisfactory */
121			for ( ; ; ) {
122			mkcumf(); /* sync brightness mapping */
123			above = below = 0; /* compute visibility overflow */
124			for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
125			hdiffs[i] = modhist[i] - (int)(T*clampf(Lb(b)) + .5);
126			if (hdiffs[i] > 0) above += hdiffs[i];
127			else below -= hdiffs[i];
128			}
129			if (above <= histot*CVRATIO)
130			break; /* close enough */
131			if (above-below >= 0)
132			return(-1); /* Houston, we have a problem.... */
133			/* original looped here as well (BEGIN_L2) */
134			maxd = 0; /* find largest overvis */
135			for (i = 0; i < HISTRES; i++)
136			if (hdiffs[i] > maxd)
137			maxd = hdiffs[maxi=i];
138			/* broke loop here when (maxd == 0) (BREAK_L2) */
139			for (sd = shiftdir((maxi+.5)/HISTRES*(bwmax-bwmin)+bwmin);
140			hdiffs[maxi] == maxd; sd = -sd)
141			for (i = maxi+sd; i >= 0 & i < HISTRES; i += sd)
142			if (hdiffs[i] < 0) {
143			if (hdiffs[i] <= -maxd) {
144			modhist[i] += maxd;
145			modhist[maxi] -= maxd;
146			hdiffs[i] += maxd;
147			hdiffs[maxi] = 0;
148			} else {
149			modhist[maxi] += hdiffs[i];
150			modhist[i] -= hdiffs[i];
151			hdiffs[maxi] += hdiffs[i];
152			hdiffs[i] = 0;
153			}
154			break;
155			}
156			/* (END_L2) */
157			}
158			return(0);
159			}
160
161
162			scotscan(scan, xres) /* apply scotopic color sensitivity loss */
163			COLOR *scan;
164			int xres;
165			{
166			COLOR ctmp;
167			double incolor, b, Lw;
168			register int i;
169
170			for (i = 0; i < xres; i++) {
171			Lw = plum(scan[i]);
172			if (Lw >= TopMesopic)
173			incolor = 1.;
174			else if (Lw <= BotMesopic)
175			incolor = 0.;
176			else
177			incolor = (Lw - BotMesopic) /
178			(TopMesopic - BotMesopic);
179			if (incolor < 1.-FTINY) {
180	greg	3.2	b = (1.-incolor)slum(scan[i])inpexp/SWNORM;
181			if (lumf == rgblum) b /= WHTEFFICACY;
182			setcolor(ctmp, b, b, b);
183	greg	3.1	if (incolor <= FTINY)
184			setcolor(scan[i], 0., 0., 0.);
185			else
186			scalecolor(scan[i], incolor);
187			addcolor(scan[i], ctmp);
188			}
189			}
190			}
191
192
193			mapscan(scan, xres) /* apply tone mapping operator to scanline */
194			COLOR *scan;
195			int xres;
196			{
197			double mult, Lw, b;
198			register int i;
199
200			for (i = 0; i < xres; i++) {
201			Lw = plum(scan[i]);
202			if (Lw < LMIN) {
203			setcolor(scan[i], 0., 0., 0.);
204			continue;
205			}
206			b = BLw(Lw);
207			mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
208			if (lumf == rgblum) mult *= WHTEFFICACY;
209			scalecolor(scan[i], mult);
210			}
211			}
212
213
214			#ifdef DEBUG
215			doplots() /* generate debugging plots */
216			{
217			double T, b, s;
218			FILE *fp;
219			char fname[128];
220			register int i;
221
222			T = histot * (bwmax - bwmin) / HISTRES;
223			s = (bwmax - bwmin)/HISTRES;
224
225			sprintf(fname, "%s_hist.plt", infn);
226			if ((fp = fopen(fname, "w")) == NULL)
227			syserror(fname);
228			fputs("include=curve.plt\n", fp);
229			fputs("title=\"Brightness Frequency Distribution\"\n", fp);
230			fprintf(fp, "subtitle=%s\n", infn);
231			fputs("ymin=0\n", fp);
232			fputs("xlabel=\"Perceptual Brightness B(Lw)\"\n", fp);
233			fputs("ylabel=\"Frequency Count\"\n", fp);
234			fputs("Alabel=\"Histogram\"\n", fp);
235			fputs("Alintype=0\n", fp);
236			fputs("Blabel=\"Envelope\"\n", fp);
237			fputs("Bsymsize=0\n", fp);
238			fputs("Adata=\n", fp);
239			for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
240			fprintf(fp, "\t%f %d\n", b, modhist[i]);
241			fputs(";\nBdata=\n", fp);
242			for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
243			fprintf(fp, "\t%f %f\n", b, T*clampf(Lb(b)));
244			fputs(";\n", fp);
245			fclose(fp);
246			}
247			#endif