ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/px/pcond3.c

Comparing ray/src/px/pcond3.c (file contents):
Revision 3.10 by gwlarson, Thu Mar 12 15:47:34 1998 UTC vs.
Revision 3.18 by greg, Sat Jun 27 03:37:24 2020 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1997 Regents of the University of California */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ LBL";
2	>	static const char       RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		* Routines for computing and applying brightness mapping.
6		*/
7
8	+	#include <string.h>
9	+
10		#include "pcond.h"
11	+	#include "random.h"
12
13
14		#define CVRATIO         0.025           /* fraction of samples allowed > env. */
#	Line 16 | Line 16 | static char SCCSid[] = "$SunId$ LBL";
16		#define LN_10           2.30258509299404568402
17		#define exp10(x)        exp(LN_10*(x))
18
19	<	float   modhist[HISTRES];               /* modified histogram */
19	>	double  modhist[HISTRES];               /* modified histogram */
20		double  mhistot;                        /* modified histogram total */
21	<	float   cumf[HISTRES+1];                /* cumulative distribution function */
21	>	double  cumf[HISTRES+1];                /* cumulative distribution function */
22
23	+	static double centprob(int x, int y);
24	+	static void mkcumf(void);
25	+	static double cf(double b);
26	+	static double BLw(double Lw);
27	+	static float *getlumsamp(int n);
28	+	#if ADJ_VEIL
29	+	static void mkcrfimage(void);
30	+	#endif
31
32	<	getfixations(fp)                /* load fixation history list */
33	<	FILE    *fp;
32	>
33	>
34	>	void
35	>	getfixations(           /* load fixation history list */
36	>	FILE    *fp
37	>	)
38		{
39		#define FIXHUNK         128
40		RESOLU  fvres;
41		int     pos[2];
42	<	register int    px, py, i;
42	>	int     px, py, i;
43		/* initialize our resolution struct */
44	<	if ((fvres.or=inpres.or)&YMAJOR) {
44	>	if ((fvres.rt=inpres.rt)&YMAJOR) {
45		fvres.xr = fvxr;
46		fvres.yr = fvyr;
47		} else {
#	Line 57 | Line 69 | FILE   *fp;
69		if (nfixations % FIXHUNK == 0) {
70		if (nfixations)
71		fixlst = (short (*)[2])
72	<	realloc((char *)fixlst,
72	>	realloc((void *)fixlst,
73		(nfixations+FIXHUNK)*
74		2*sizeof(short));
75		else
#	Line 82 | Line 94 | FILE   *fp;
94		}
95
96
97	<	gethisto(fp)                    /* load precomputed luminance histogram */
98	<	FILE    *fp;
97	>	void
98	>	gethisto(                       /* load precomputed luminance histogram */
99	>	FILE    *fp
100	>	)
101		{
102	<	float   histo[MAXPREHIST];
102	>	double  histo[MAXPREHIST];
103		double  histart, histep;
104	<	double  l, b, lastb, w;
104	>	double  b, lastb, w;
105		int     n;
106	<	register int    i;
106	>	int     i;
107		/* load data */
108		for (i = 0; i < MAXPREHIST &&
109	<	fscanf(fp, "%lf %f", &b, &histo[i]) == 2; i++) {
110	<	if (i > 1 && fabs(b - lastb - histep) > 1e-3) {
109	>	fscanf(fp, "%lf %lf", &b, &histo[i]) == 2; i++) {
110	>	if (i > 1 && fabs(b - lastb - histep) > .001) {
111		fprintf(stderr,
112		"%s: uneven step size in histogram data\n",
113		progname);
114		exit(1);
115		}
116		if (i == 1)
117	<	histep = b - (histart = lastb);
117	>	if ((histep = b - (histart = lastb)) <= FTINY) {
118	>	fprintf(stderr,
119	>	"%s: illegal step in histogram data\n",
120	>	progname);
121	>	exit(1);
122	>	}
123		lastb = b;
124		}
125		if (i < 2 || !feof(fp)) {
#	Line 115 | Line 134 | FILE   *fp;
134		/* find extrema */
135		for (i = 0; i < n && histo[i] <= FTINY; i++)
136		;
137	<	bwmin = histart + i*histep;
137	>	bwmin = histart + (i-.001)*histep;
138		for (i = n; i-- && histo[i] <= FTINY; )
139		;
140	<	bwmax = histart + i*histep;
140	>	bwmax = histart + (i+1.001)*histep;
141		if (bwmax > Bl(LMAX))
142		bwmax = Bl(LMAX);
143		if (bwmin < Bl(LMIN))
#	Line 130 | Line 149 | FILE   *fp;
149		for (i = 0; i < HISTRES; i++)
150		bwhist[i] = 0.;
151		for (i = 0, b = histart; i < n; i++, b += histep) {
152	<	if (b < bwmin) continue;
153	<	if (b > bwmax) break;
152	>	if (b < bwmin+FTINY) continue;
153	>	if (b >= bwmax-FTINY) break;
154		w = histo[i];
155		bwavg += w*b;
156		bwhist[bwhi(b)] += w;
#	Line 145 | Line 164 | FILE   *fp;
164		}
165
166
167	<	double
168	<	centprob(x, y)                  /* center-weighting probability function */
169	<	int     x, y;
167	>	static double
168	>	centprob(                       /* center-weighting probability function */
169	>	int     x,
170	>	int     y
171	>	)
172		{
173		double  xr, yr, p;
174		/* paraboloid, 0 at 90 degrees from center */
#	Line 158 | Line 179 | int    x, y;
179		}
180
181
182	<	comphist()                      /* create foveal sampling histogram */
182	>	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	+	if (x <= 0 && freadcolrs(cscan, x=scanlen(&inpres), infp) < 0) {
200	+	fprintf(stderr, "%s: %s: scanline read error\n",
201	+	progname, infn);
202	+	exit(1);
203	+	}
204	+	rval = frandom();       /* random distance to next sample */
205	+	while ((cumprob += (1.-cumprob)*n/(npleft-sv)) < rval)
206	+	sv++;
207	+	if (x < ++sv) {         /* out of pixels in this scanline */
208	+	npleft -= x;
209	+	x = 0;
210	+	continue;
211	+	}
212	+	x -= sv;
213	+	colr_color(col, cscan[x]);
214	+	ls[--n] = plum(col);
215	+	npleft -= sv;
216	+	}
217	+	free(cscan);                    /* clean up and reset file pointer */
218	+	if (fseek(infp, startpos, SEEK_SET) < 0)
219	+	syserror("fseek");
220	+	return(ls);
221	+	}
222	+
223	+
224	+	void
225	+	comphist(void)                  /* create foveal sampling histogram */
226	+	{
227		double  l, b, w, lwmin, lwmax;
228	<	register int    x, y;
228	>	float   *lumsamp;
229	>	int     nlumsamp;
230	>	int     x, y;
231		/* check for precalculated histogram */
232		if (what2do&DO_PREHIST)
233		return;
#	Line 173 | Line 239 | comphist()                     /* create foveal sampling histogram */
239		if (l < lwmin) lwmin = l;
240		if (l > lwmax) lwmax = l;
241		}
242	+	/* sample luminance pixels */
243	+	nlumsamp = fvxr*fvyr*16;
244	+	if (nlumsamp > inpres.xr*inpres.yr)
245	+	nlumsamp = inpres.xr*inpres.yr;
246	+	lumsamp = getlumsamp(nlumsamp);
247	+	for (x = nlumsamp; x--; ) {
248	+	l = lumsamp[x];
249	+	if (l < lwmin) lwmin = l;
250	+	if (l > lwmax) lwmax = l;
251	+	}
252		lwmax *= 1.01;
253		if (lwmax > LMAX)
254		lwmax = LMAX;
#	Line 187 | Line 263 | comphist()                     /* create foveal sampling histogram */
263		/* (re)compute histogram */
264		bwavg = 0.;
265		histot = 0.;
266	<	for (x = 0; x < HISTRES; x++)
191	<	bwhist[x] = 0.;
266	>	memset(bwhist, 0, sizeof(bwhist));
267		/* global average */
268	<	if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
268	>	if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY) {
269		for (y = 0; y < fvyr; y++)
270		for (x = 0; x < fvxr; x++) {
271		l = plum(fovscan(y)[x]);
272	<	if (l < lwmin) continue;
273	<	if (l > lwmax) continue;
272	>	if (l < lwmin+FTINY) continue;
273	>	if (l >= lwmax-FTINY) continue;
274		b = Bl(l);
275		w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
276		bwavg += w*b;
277		bwhist[bwhi(b)] += w;
278		histot += w;
279		}
280	+	/* weight for point luminances */
281	+	w = 1. * histot / nlumsamp;
282	+	for (x = nlumsamp; x--; ) {
283	+	l = lumsamp[x];
284	+	if (l < lwmin+FTINY) continue;
285	+	if (l >= lwmax-FTINY) continue;
286	+	b = Bl(l);
287	+	bwavg += w*b;
288	+	bwhist[bwhi(b)] += w;
289	+	histot += w;
290	+	}
291	+	}
292		/* average fixation points */
293		if (what2do&DO_FIXHIST && nfixations > 0) {
294		if (histot > FTINY)
#	Line 210 | Line 297 | comphist()                     /* create foveal sampling histogram */
297		w = 1.;
298		for (x = 0; x < nfixations; x++) {
299		l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
300	<	if (l < lwmin) continue;
301	<	if (l > lwmax) continue;
300	>	if (l < lwmin+FTINY) continue;
301	>	if (l >= lwmax-FTINY) continue;
302		b = Bl(l);
303		bwavg += w*b;
304		bwhist[bwhi(b)] += w;
#	Line 223 | Line 310 | comphist()                     /* create foveal sampling histogram */
310		bwhist[1] *= 0.5;
311		bwhist[0] += bwhist[1];
312		}
313	+	free(lumsamp);
314		}
315
316
317	<	mkcumf()                        /* make cumulative distribution function */
317	>	static void
318	>	mkcumf(void)                    /* make cumulative distribution function */
319		{
320	<	register int    i;
321	<	register double sum;
320	>	int     i;
321	>	double  sum;
322
323		mhistot = 0.;           /* compute modified total */
324		for (i = 0; i < HISTRES; i++)
#	Line 244 | Line 333 | mkcumf()                       /* make cumulative distribution function */
333		}
334
335
336	<	double
337	<	cf(b)                           /* return cumulative function at b */
338	<	double  b;
336	>	static double
337	>	cf(                             /* return cumulative function at b */
338	>	double  b
339	>	)
340		{
341		double  x;
342	<	register int    i;
342	>	int     i;
343
344		i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
345		x -= (double)i;
#	Line 257 | Line 347 | double b;
347		}
348
349
350	<	double
351	<	BLw(Lw)                         /* map world luminance to display brightness */
352	<	double  Lw;
350	>	static double
351	>	BLw(                            /* map world luminance to display brightness */
352	>	double  Lw
353	>	)
354		{
355		double  b;
356
#	Line 272 | Line 363 | double Lw;
363
364
365		double
366	<	htcontrs(La)            /* human threshold contrast sensitivity, dL(La) */
367	<	double  La;
366	>	htcontrs(               /* human threshold contrast sensitivity, dL(La) */
367	>	double  La
368	>	)
369		{
370		double  l10La, l10dL;
371		/* formula taken from Ferwerda et al. [SG96] */
#	Line 294 | Line 386 | double La;
386
387
388		double
389	<	clampf(Lw)              /* derivative clamping function */
390	<	double  Lw;
389	>	clampf(                 /* histogram clamping function */
390	>	double  Lw
391	>	)
392		{
393		double  bLw, ratio;
394
#	Line 304 | Line 397 | double Lw;
397		return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
398		}
399
400	+	double
401	+	crfactor(                       /* contrast reduction factor */
402	+	double  Lw
403	+	)
404	+	{
405	+	int     i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
406	+	double  bLw, ratio, Tdb;
407
408	+	if (i <= 0)
409	+	return(1.0);
410	+	if (i >= HISTRES)
411	+	return(1.0);
412	+	bLw = BLw(Lw);
413	+	ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
414	+	Tdb = mhistot * (bwmax - bwmin) / HISTRES;
415	+	return(modhist[i]*Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)/(Tdb*ratio));
416	+	}
417	+
418	+
419	+	#if ADJ_VEIL
420	+	static void
421	+	mkcrfimage(void)                        /* compute contrast reduction factor image */
422	+	{
423	+	int     i;
424	+	float   *crfptr;
425	+	COLOR   *fovptr;
426	+
427	+	if (crfimg == NULL)
428	+	crfimg = (float *)malloc(fvxr*fvyr*sizeof(float));
429	+	if (crfimg == NULL)
430	+	syserror("malloc");
431	+	crfptr = crfimg;
432	+	fovptr = fovimg;
433	+	for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
434	+	crfptr[0] = crfactor(plum(fovptr[0]));
435	+	}
436	+	#endif
437	+
438	+
439		int
440	<	mkbrmap()                       /* make dynamic range map */
440	>	mkbrmap(void)                   /* make dynamic range map */
441		{
442	<	double  T, b, s;
442	>	double  Tdb, b, s;
443		double  ceiling, trimmings;
444	<	register int    i;
444	>	int     i;
445		/* copy initial histogram */
446	<	bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
447	<	s = (bwmax - bwmin)/HISTRES;
446	>	memcpy((void *)modhist, (void *)bwhist, sizeof(modhist));
447	>	s = (bwmax - bwmin)/HISTRES;    /* s is delta b */
448		/* loop until satisfactory */
449		do {
450		mkcumf();                       /* sync brightness mapping */
451		if (mhistot <= histot*CVRATIO)
452		return(-1);             /* no compression needed! */
453	<	T = mhistot * (bwmax - bwmin) / HISTRES;
453	>	Tdb = mhistot * s;
454		trimmings = 0.;                 /* clip to envelope */
455		for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
456	<	ceiling = T*clampf(Lb(b));
456	>	ceiling = Tdb*clampf(Lb(b));
457		if (modhist[i] > ceiling) {
458		trimmings += modhist[i] - ceiling;
459		modhist[i] = ceiling;
#	Line 330 | Line 461 | mkbrmap()                      /* make dynamic range map */
461		}
462		} while (trimmings > histot*CVRATIO);
463
464	+	#if ADJ_VEIL
465	+	mkcrfimage();                   /* contrast reduction image */
466	+	#endif
467	+
468		return(0);                      /* we got it */
469		}
470
471
472	<	scotscan(scan, xres)            /* apply scotopic color sensitivity loss */
473	<	COLOR   *scan;
474	<	int     xres;
472	>	void
473	>	scotscan(               /* apply scotopic color sensitivity loss */
474	>	COLOR   *scan,
475	>	int     xres
476	>	)
477		{
478		COLOR   ctmp;
479		double  incolor, b, Lw;
480	<	register int    i;
480	>	int     i;
481
482		for (i = 0; i < xres; i++) {
483		Lw = plum(scan[i]);
#	Line 365 | Line 502 | int    xres;
502		}
503
504
505	<	mapscan(scan, xres)             /* apply tone mapping operator to scanline */
506	<	COLOR   *scan;
507	<	int     xres;
505	>	void
506	>	mapscan(                /* apply tone mapping operator to scanline */
507	>	COLOR   *scan,
508	>	int     xres
509	>	)
510		{
511		double  mult, Lw, b;
512	<	register int    i;
512	>	int     x;
513
514	<	for (i = 0; i < xres; i++) {
515	<	Lw = plum(scan[i]);
514	>	for (x = 0; x < xres; x++) {
515	>	Lw = plum(scan[x]);
516		if (Lw < LMIN) {
517	<	setcolor(scan[i], 0., 0., 0.);
517	>	setcolor(scan[x], 0., 0., 0.);
518		continue;
519		}
520	<	b = BLw(Lw);
521	<	mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
520	>	b = BLw(Lw);            /* apply brightness mapping */
521	>	mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
522		if (lumf == rgblum) mult *= WHTEFFICACY;
523	<	scalecolor(scan[i], mult);
523	>	scalecolor(scan[x], mult);
524		}
525		}
526
527
528	<	putmapping(fp)                  /* put out mapping function */
529	<	FILE    *fp;
528	>	void
529	>	putmapping(                     /* put out mapping function */
530	>	FILE    *fp
531	>	)
532		{
533		double  b, s;
534	<	register int    i;
534	>	int     i;
535		double  wlum, sf, dlum;
536
537		sf = scalef*inpexp;

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)