ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pcond3.c
(Generate patch)

Comparing ray/src/px/pcond3.c (file contents):
Revision 3.4 by greg, Wed Jan 8 17:57:01 1997 UTC vs.
Revision 3.15 by schorsch, Sun Mar 28 20:33:14 2004 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1996 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  
12  
13 < #define CVRATIO         0.025           /* fraction of pixels allowed > env. */
13 > #define CVRATIO         0.025           /* fraction of samples allowed > env. */
14  
15 < #define exp10(x)        exp(2.302585093*(x))
15 > #define LN_10           2.30258509299404568402
16 > #define exp10(x)        exp(LN_10*(x))
17  
18 < int     modhist[HISTRES];               /* modified histogram */
19 < float   cumf[HISTRES+1];                /* cumulative distribution function */
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 < mkcumf()                        /* make cumulative distribution function */
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 <        register long   sum;
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 <        cumf[0] = 0.;
165 <        sum = modhist[0];
166 <        for (i = 1; i < HISTRES; i++) {
167 <                cumf[i] = (double)sum/histot;
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 < double
269 < cf(b)                           /* return cumulative function at b */
270 < double  b;
268 > static double
269 > cf(                             /* return cumulative function at b */
270 >        double  b
271 > )
272   {
273          double  x;
274          register int    i;
# Line 47 | Line 279 | double b;
279   }
280  
281  
282 < double
283 < BLw(Lw)                         /* map world luminance to display brightness */
284 < double  Lw;
282 > static double
283 > BLw(                            /* map world luminance to display brightness */
284 >        double  Lw
285 > )
286   {
287          double  b;
288  
# Line 61 | Line 294 | double Lw;
294   }
295  
296  
297 < double
298 < htcontrs(La)            /* human threshold contrast sensitivity, dL(La) */
299 < double  La;
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] */
# Line 83 | Line 317 | double La;
317   }
318  
319  
320 < double
321 < clampf(Lw)              /* derivative clamping function */
322 < double  Lw;
320 > extern double
321 > clampf(                 /* histogram clamping function */
322 >        double  Lw
323 > )
324   {
325          double  bLw, ratio;
326  
# Line 94 | Line 329 | double 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 < int
341 < shiftdir(bw)            /* compute shift direction for histogram */
342 < double  bw;
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 <        if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
356 <                return(1);
357 <        return(-1);
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 < int
372 < mkbrmap()                       /* make dynamic range map */
371 > extern int
372 > mkbrmap(void)                   /* make dynamic range map */
373   {
374 <        int     hdiffs[HISTRES], above, below;
375 <        double  T, b, s;
113 <        int     maxd, maxi, sd;
374 >        double  Tdb, b, s;
375 >        double  ceiling, trimmings;
376          register int    i;
377                                          /* copy initial histogram */
378 <        for (i = 0; i < HISTRES; i++)
379 <                modhist[i] = bwhist[i];
118 <        T = histot * (bwmax - bwmin) / HISTRES;
119 <        s = (bwmax - bwmin)/HISTRES;
378 >        memcpy((void *)modhist, (void *)bwhist, sizeof(modhist));
379 >        s = (bwmax - bwmin)/HISTRES;    /* s is delta b */
380                                          /* loop until satisfactory */
381 <        for ( ; ; ) {
382 <                mkcumf();               /* sync brightness mapping */
383 <                above = below = 0;      /* compute visibility overflow */
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 <                        hdiffs[i] = modhist[i] - (int)(T*clampf(Lb(b)) + .5);
389 <                        if (hdiffs[i] > 0) above += hdiffs[i];
390 <                        else below -= hdiffs[i];
388 >                        ceiling = Tdb*clampf(Lb(b));
389 >                        if (modhist[i] > ceiling) {
390 >                                trimmings += modhist[i] - ceiling;
391 >                                modhist[i] = ceiling;
392 >                        }
393                  }
394 <                if (above <= histot*CVRATIO)
395 <                        break;          /* close enough */
396 <                if (above-below >= 0)
397 <                        return(-1);     /* Houston, we have a problem.... */
398 <                /* original looped here as well (BEGIN_L2) */
399 <                maxd = 0;               /* find largest overvis */
400 <                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);
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 < scotscan(scan, xres)            /* apply scotopic color sensitivity loss */
405 < COLOR   *scan;
406 < int     xres;
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;
# Line 190 | Line 434 | int    xres;
434   }
435  
436  
437 < mapscan(scan, xres)             /* apply tone mapping operator to scanline */
438 < COLOR   *scan;
439 < int     xres;
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    i;
444 >        register int    x;
445  
446 <        for (i = 0; i < xres; i++) {
447 <                Lw = plum(scan[i]);
446 >        for (x = 0; x < xres; x++) {
447 >                Lw = plum(scan[x]);
448                  if (Lw < LMIN) {
449 <                        setcolor(scan[i], 0., 0., 0.);
449 >                        setcolor(scan[x], 0., 0., 0.);
450                          continue;
451                  }
452 <                b = BLw(Lw);
453 <                mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
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[i], mult);
455 >                scalecolor(scan[x], mult);
456          }
457   }
458  
459  
460 < putmapping(fp)                  /* put out mapping function */
461 < FILE    *fp;
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;
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 <                        fprintf(fp, "%e %e\n", wlum, sf*wlum);
476 <                else
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   }
232
233
234 #ifdef DEBUG
235 doplots()                       /* generate debugging plots */
236 {
237        double  T, b, s;
238        FILE    *fp;
239        char    fname[128];
240        register int    i;
241
242        T = histot * (bwmax - bwmin) / HISTRES;
243        s = (bwmax - bwmin)/HISTRES;
244
245        sprintf(fname, "%s_hist.plt", infn);
246        if ((fp = fopen(fname, "w")) == NULL)
247                syserror(fname);
248        fputs("include=curve.plt\n", fp);
249        fputs("title=\"Brightness Frequency Distribution\"\n", fp);
250        fprintf(fp, "subtitle=%s\n", infn);
251        fputs("ymin=0\n", fp);
252        fputs("xlabel=\"Perceptual Brightness B(Lw)\"\n", fp);
253        fputs("ylabel=\"Frequency Count\"\n", fp);
254        fputs("Alabel=\"Histogram\"\n", fp);
255        fputs("Alintype=0\n", fp);
256        fputs("Blabel=\"Envelope\"\n", fp);
257        fputs("Bsymsize=0\n", fp);
258        fputs("Adata=\n", fp);
259        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
260                fprintf(fp, "\t%f %d\n", b, modhist[i]);
261        fputs(";\nBdata=\n", fp);
262        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
263                fprintf(fp, "\t%f %f\n", b, T*clampf(Lb(b)));
264        fputs(";\n", fp);
265        fclose(fp);
266 }
267 #endif

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines