src/px/pcond3.c

#ifndef lint
static const char       RCSid[] = "$Id: pcond3.c,v 3.19 2021/04/13 02:56:42 greg Exp $";
#endif
/*
 * Routines for computing and applying brightness mapping.
 */

#include <string.h>

#include "pcond.h"
#include "random.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 */

static double centprob(int x, int y);
static void mkcumf(void);
static double cf(double b);
static double BLw(double Lw);
static float *getlumsamp(int n);
#if ADJ_VEIL
static void mkcrfimage(void);
#endif


void
getfixations(           /* load fixation history list */
FILE    *fp
)
{
#define FIXHUNK         128
        RESOLU  fvres;
        int     pos[2];
        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((void *)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
}


void
gethisto(                       /* load precomputed luminance histogram */
        FILE    *fp
)
{
        double  histo[MAXPREHIST];
        double  histart, histep;
        double  b, lastb, w;
        int     n;
        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];
        }
}


static double
centprob(                       /* center-weighting probability function */
        int     x,
        int     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);
}


static float *
getlumsamp(int n)               /* get set of random point sample luminances */
{
        float   *ls = (float *)malloc(n*sizeof(float));
        COLR    *cscan = (COLR *)malloc(scanlen(&inpres)*sizeof(COLR));
        long    startpos = ftell(infp);
        long    npleft = (long)inpres.xr*inpres.yr;
        int     x;
        
        if ((ls == NULL) | (cscan == NULL))
                syserror("malloc");
        x = 0;                          /* read/convert samples */
        while (n > 0) {
                COLOR   col;
                int     sv = 0;
                double  rval, cumprob = 0;

                if (x <= 0 && fread2colrs(cscan, x=scanlen(&inpres), infp,
                                                NCSAMP, WLPART) < 0) {
                        fprintf(stderr, "%s: %s: scanline read error\n",
                                        progname, infn);
                        exit(1);
                }
                rval = frandom();       /* random distance to next sample */
                while ((cumprob += (1.-cumprob)*n/(npleft-sv)) < rval)
                        sv++;
                if (x < ++sv) {         /* out of pixels in this scanline */
                        npleft -= x;
                        x = 0;
                        continue;
                }
                x -= sv;
                colr_color(col, cscan[x]);
                ls[--n] = plum(col);
                npleft -= sv;
        }
        free(cscan);                    /* clean up and reset file pointer */
        if (fseek(infp, startpos, SEEK_SET) < 0)
                syserror("fseek");
        return(ls);
}


void
comphist(void)                  /* create foveal sampling histogram */
{
        double  l, b, w, lwmin, lwmax;
        float   *lumsamp;
        int     nlumsamp;
        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;
                }
                                        /* sample luminance pixels */
        nlumsamp = fvxr*fvyr*16;
        if (nlumsamp > inpres.xr*inpres.yr)
                nlumsamp = inpres.xr*inpres.yr;
        lumsamp = getlumsamp(nlumsamp);
        for (x = nlumsamp; x--; ) {
                l = lumsamp[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.;
        memset(bwhist, 0, sizeof(bwhist));
                                        /* 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;
                        }
                                        /* weight for point luminances */
                w = 1. * histot / nlumsamp;
                for (x = nlumsamp; x--; ) {
                        l = lumsamp[x];
                        if (l < lwmin+FTINY) continue;
                        if (l >= lwmax-FTINY) continue;
                        b = Bl(l);
                        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];
        }
        free(lumsamp);
}


static void
mkcumf(void)                    /* make cumulative distribution function */
{
        int     i;
        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.;
}


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

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


static double
BLw(                            /* 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(               /* 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(                 /* 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(                       /* 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
static void
mkcrfimage(void)                        /* 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(void)                   /* make dynamic range map */
{
        double  Tdb, b, s;
        double  ceiling, trimmings;
        int     i;
                                        /* copy initial histogram */
        memcpy((void *)modhist, (void *)bwhist, 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 > mhistot*CVRATIO);

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

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


void
scotscan(               /* apply scotopic color sensitivity loss */
        COLOR   *scan,
        int     xres
)
{
        COLOR   ctmp;
        double  incolor, b, Lw;
        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);
                }
        }
}


void
mapscan(                /* apply tone mapping operator to scanline */
        COLOR   *scan,
        int     xres
)
{
        double  mult, Lw, b;
        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);
        }
}


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