ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pcond3.c
Revision: 3.14
Committed: Mon Jun 30 14:59:12 2003 UTC (20 years, 10 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 3.13: +4 -2 lines
Log Message:
Replaced most outdated BSD function calls with their posix equivalents, and cleaned up a few other platform dependencies.

File Contents

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