ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pcond3.c
Revision: 3.15
Committed: Sun Mar 28 20:33:14 2004 UTC (20 years ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad3R7P2, rad3R7P1, rad4R2, rad4R1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9, rad4R2P1
Changes since 3.14: +66 -37 lines
Log Message:
Continued ANSIfication, and other fixes and clarifications.

File Contents

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