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.1 by greg, Thu Oct 3 16:52:50 1996 UTC vs.
Revision 3.11 by greg, Sat Feb 22 02:07:27 2003 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   */
# Line 11 | Line 8 | static char SCCSid[] = "$SunId$ LBL";
8   #include "pcond.h"
9  
10  
11 < #define CVRATIO         0.025           /* fraction of pixels allowed > env. */
11 > #define CVRATIO         0.025           /* fraction of samples allowed > env. */
12  
13 < #define BotMesopic      5.62e-3         /* top of scotopic range */
14 < #define TopMesopic      5.62            /* bottom of photopic range */
13 > #define LN_10           2.30258509299404568402
14 > #define exp10(x)        exp(LN_10*(x))
15  
16 < #define exp10(x)        exp(2.302585093*(x))
16 > double  modhist[HISTRES];               /* modified histogram */
17 > double  mhistot;                        /* modified histogram total */
18 > double  cumf[HISTRES+1];                /* cumulative distribution function */
19  
21 – int     modhist[HISTRES];               /* modified histogram */
22 – float   cumf[HISTRES+1];                /* cumulative distribution function */
20  
21 + getfixations(fp)                /* load fixation history list */
22 + FILE    *fp;
23 + {
24 + #define FIXHUNK         128
25 +        RESOLU  fvres;
26 +        int     pos[2];
27 +        register int    px, py, i;
28 +                                /* initialize our resolution struct */
29 +        if ((fvres.rt=inpres.rt)&YMAJOR) {
30 +                fvres.xr = fvxr;
31 +                fvres.yr = fvyr;
32 +        } else {
33 +                fvres.xr = fvyr;
34 +                fvres.yr = fvxr;
35 +        }
36 +                                /* read each picture position */
37 +        while (fscanf(fp, "%d %d", &pos[0], &pos[1]) == 2) {
38 +                                /* convert to closest index in foveal image */
39 +                loc2pix(pos, &fvres,
40 +                                (pos[0]+.5)/inpres.xr, (pos[1]+.5)/inpres.yr);
41 +                                /* include nine neighborhood samples */
42 +                for (px = pos[0]-1; px <= pos[0]+1; px++) {
43 +                        if (px < 0 || px >= fvxr)
44 +                                continue;
45 +                        for (py = pos[1]-1; py <= pos[1]+1; py++) {
46 +                                if (py < 0 || py >= fvyr)
47 +                                        continue;
48 +                                for (i = nfixations; i-- > 0; )
49 +                                        if (fixlst[i][0] == px &&
50 +                                                        fixlst[i][1] == py)
51 +                                                break;
52 +                                if (i >= 0)
53 +                                        continue;       /* already there */
54 +                                if (nfixations % FIXHUNK == 0) {
55 +                                        if (nfixations)
56 +                                                fixlst = (short (*)[2])
57 +                                                        realloc((char *)fixlst,
58 +                                                        (nfixations+FIXHUNK)*
59 +                                                        2*sizeof(short));
60 +                                        else
61 +                                                fixlst = (short (*)[2])malloc(
62 +                                                        FIXHUNK*2*sizeof(short)
63 +                                                        );
64 +                                        if (fixlst == NULL)
65 +                                                syserror("malloc");
66 +                                }
67 +                                fixlst[nfixations][0] = px;
68 +                                fixlst[nfixations][1] = py;
69 +                                nfixations++;
70 +                        }
71 +                }
72 +        }
73 +        if (!feof(fp)) {
74 +                fprintf(stderr, "%s: format error reading fixation data\n",
75 +                                progname);
76 +                exit(1);
77 +        }
78 + #undef  FIXHUNK
79 + }
80  
81 +
82 + gethisto(fp)                    /* load precomputed luminance histogram */
83 + FILE    *fp;
84 + {
85 +        double  histo[MAXPREHIST];
86 +        double  histart, histep;
87 +        double  l, b, lastb, w;
88 +        int     n;
89 +        register int    i;
90 +                                        /* load data */
91 +        for (i = 0; i < MAXPREHIST &&
92 +                        fscanf(fp, "%lf %lf", &b, &histo[i]) == 2; i++) {
93 +                if (i > 1 && fabs(b - lastb - histep) > .001) {
94 +                        fprintf(stderr,
95 +                                "%s: uneven step size in histogram data\n",
96 +                                        progname);
97 +                        exit(1);
98 +                }
99 +                if (i == 1)
100 +                        if ((histep = b - (histart = lastb)) <= FTINY) {
101 +                                fprintf(stderr,
102 +                                        "%s: illegal step in histogram data\n",
103 +                                                progname);
104 +                                exit(1);
105 +                        }
106 +                lastb = b;
107 +        }
108 +        if (i < 2 || !feof(fp)) {
109 +                fprintf(stderr,
110 +                "%s: format/length error loading histogram (log10L %f at %d)\n",
111 +                                progname, b, i);
112 +                exit(1);
113 +        }
114 +        n = i;
115 +        histart *= LN_10;
116 +        histep *= LN_10;
117 +                                        /* find extrema */
118 +        for (i = 0; i < n && histo[i] <= FTINY; i++)
119 +                ;
120 +        bwmin = histart + (i-.001)*histep;
121 +        for (i = n; i-- && histo[i] <= FTINY; )
122 +                ;
123 +        bwmax = histart + (i+1.001)*histep;
124 +        if (bwmax > Bl(LMAX))
125 +                bwmax = Bl(LMAX);
126 +        if (bwmin < Bl(LMIN))
127 +                bwmin = Bl(LMIN);
128 +        else                            /* duplicate bottom bin */
129 +                bwmin = bwmax - (bwmax-bwmin)*HISTRES/(HISTRES-1);
130 +                                        /* convert histogram */
131 +        bwavg = 0.; histot = 0.;
132 +        for (i = 0; i < HISTRES; i++)
133 +                bwhist[i] = 0.;
134 +        for (i = 0, b = histart; i < n; i++, b += histep) {
135 +                if (b < bwmin+FTINY) continue;
136 +                if (b >= bwmax-FTINY) break;
137 +                w = histo[i];
138 +                bwavg += w*b;
139 +                bwhist[bwhi(b)] += w;
140 +                histot += w;
141 +        }
142 +        bwavg /= histot;
143 +        if (bwmin > Bl(LMIN)+FTINY) {   /* add false samples at bottom */
144 +                bwhist[1] *= 0.5;
145 +                bwhist[0] += bwhist[1];
146 +        }
147 + }
148 +
149 +
150 + double
151 + centprob(x, y)                  /* center-weighting probability function */
152 + int     x, y;
153 + {
154 +        double  xr, yr, p;
155 +                                /* paraboloid, 0 at 90 degrees from center */
156 +        xr = (x - .5*(fvxr-1))/90.;     /* 180 degree fisheye has fv?r == 90 */
157 +        yr = (y - .5*(fvyr-1))/90.;
158 +        p = 1. - xr*xr - yr*yr;
159 +        return(p < 0. ? 0. : p);
160 + }
161 +
162 +
163 + comphist()                      /* create foveal sampling histogram */
164 + {
165 +        double  l, b, w, lwmin, lwmax;
166 +        register int    x, y;
167 +                                        /* check for precalculated histogram */
168 +        if (what2do&DO_PREHIST)
169 +                return;
170 +        lwmin = 1e10;                   /* find extrema */
171 +        lwmax = 0.;
172 +        for (y = 0; y < fvyr; y++)
173 +                for (x = 0; x < fvxr; x++) {
174 +                        l = plum(fovscan(y)[x]);
175 +                        if (l < lwmin) lwmin = l;
176 +                        if (l > lwmax) lwmax = l;
177 +                }
178 +        lwmax *= 1.01;
179 +        if (lwmax > LMAX)
180 +                lwmax = LMAX;
181 +        bwmax = Bl(lwmax);
182 +        if (lwmin < LMIN) {
183 +                lwmin = LMIN;
184 +                bwmin = Bl(LMIN);
185 +        } else {                        /* duplicate bottom bin */
186 +                bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
187 +                lwmin = Lb(bwmin);
188 +        }
189 +                                        /* (re)compute histogram */
190 +        bwavg = 0.;
191 +        histot = 0.;
192 +        for (x = 0; x < HISTRES; x++)
193 +                bwhist[x] = 0.;
194 +                                        /* global average */
195 +        if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
196 +                for (y = 0; y < fvyr; y++)
197 +                        for (x = 0; x < fvxr; x++) {
198 +                                l = plum(fovscan(y)[x]);
199 +                                if (l < lwmin+FTINY) continue;
200 +                                if (l >= lwmax-FTINY) continue;
201 +                                b = Bl(l);
202 +                                w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
203 +                                bwavg += w*b;
204 +                                bwhist[bwhi(b)] += w;
205 +                                histot += w;
206 +                        }
207 +                                        /* average fixation points */
208 +        if (what2do&DO_FIXHIST && nfixations > 0) {
209 +                if (histot > FTINY)
210 +                        w = fixfrac/(1.-fixfrac)*histot/nfixations;
211 +                else
212 +                        w = 1.;
213 +                for (x = 0; x < nfixations; x++) {
214 +                        l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
215 +                        if (l < lwmin+FTINY) continue;
216 +                        if (l >= lwmax-FTINY) continue;
217 +                        b = Bl(l);
218 +                        bwavg += w*b;
219 +                        bwhist[bwhi(b)] += w;
220 +                        histot += w;
221 +                }
222 +        }
223 +        bwavg /= histot;
224 +        if (lwmin > LMIN+FTINY) {       /* add false samples at bottom */
225 +                bwhist[1] *= 0.5;
226 +                bwhist[0] += bwhist[1];
227 +        }
228 + }
229 +
230 +
231   mkcumf()                        /* make cumulative distribution function */
232   {
233          register int    i;
234 <        register long   sum;
234 >        register double sum;
235  
236 <        cumf[0] = 0.;
237 <        sum = modhist[0];
238 <        for (i = 1; i < HISTRES; i++) {
239 <                cumf[i] = (double)sum/histot;
236 >        mhistot = 0.;           /* compute modified total */
237 >        for (i = 0; i < HISTRES; i++)
238 >                mhistot += modhist[i];
239 >
240 >        sum = 0.;               /* compute cumulative function */
241 >        for (i = 0; i < HISTRES; i++) {
242 >                cumf[i] = sum/mhistot;
243                  sum += modhist[i];
244          }
245          cumf[HISTRES] = 1.;
# Line 60 | Line 269 | double Lw;
269                  return(Bldmin);
270          if (b >= bwmax-FTINY)
271                  return(Bldmax);
272 <        return(Bldmin + cf(Bl(Lw))*(Bldmax-Bldmin));
272 >        return(Bldmin + cf(b)*(Bldmax-Bldmin));
273   }
274  
275  
# Line 87 | Line 296 | double La;
296  
297  
298   double
299 < clampf(Lw)              /* derivative clamping function */
299 > clampf(Lw)                      /* histogram clamping function */
300   double  Lw;
301   {
302          double  bLw, ratio;
# Line 97 | Line 306 | double Lw;
306          return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
307   }
308  
309 + double
310 + crfactor(Lw)                    /* contrast reduction factor */
311 + double  Lw;
312 + {
313 +        int     i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
314 +        double  bLw, ratio, Tdb;
315  
316 < int
317 < shiftdir(bw)            /* compute shift direction for histogram */
318 < double  bw;
316 >        if (i <= 0)
317 >                return(1.0);
318 >        if (i >= HISTRES)
319 >                return(1.0);
320 >        bLw = BLw(Lw);
321 >        ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
322 >        Tdb = mhistot * (bwmax - bwmin) / HISTRES;
323 >        return(modhist[i]*Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)/(Tdb*ratio));
324 > }
325 >
326 >
327 > #if ADJ_VEIL
328 > mkcrfimage()                    /* compute contrast reduction factor image */
329   {
330 <        if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
331 <                return(1);
332 <        return(-1);
330 >        int     i;
331 >        float   *crfptr;
332 >        COLOR   *fovptr;
333 >
334 >        if (crfimg == NULL)
335 >                crfimg = (float *)malloc(fvxr*fvyr*sizeof(float));
336 >        if (crfimg == NULL)
337 >                syserror("malloc");
338 >        crfptr = crfimg;
339 >        fovptr = fovimg;
340 >        for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
341 >                crfptr[0] = crfactor(plum(fovptr[0]));
342   }
343 + #endif
344  
345  
346   int
347   mkbrmap()                       /* make dynamic range map */
348   {
349 <        int     hdiffs[HISTRES], above, below;
350 <        double  T, b, s;
116 <        int     maxd, maxi, sd;
349 >        double  Tdb, b, s;
350 >        double  ceiling, trimmings;
351          register int    i;
352                                          /* copy initial histogram */
353 <        for (i = 0; i < HISTRES; i++)
354 <                modhist[i] = bwhist[i];
121 <        T = histot * (bwmax - bwmin) / HISTRES;
122 <        s = (bwmax - bwmin)/HISTRES;
353 >        bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
354 >        s = (bwmax - bwmin)/HISTRES;    /* s is delta b */
355                                          /* loop until satisfactory */
356 <        for ( ; ; ) {
357 <                mkcumf();               /* sync brightness mapping */
358 <                above = below = 0;      /* compute visibility overflow */
356 >        do {
357 >                mkcumf();                       /* sync brightness mapping */
358 >                if (mhistot <= histot*CVRATIO)
359 >                        return(-1);             /* no compression needed! */
360 >                Tdb = mhistot * s;
361 >                trimmings = 0.;                 /* clip to envelope */
362                  for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
363 <                        hdiffs[i] = modhist[i] - (int)(T*clampf(Lb(b)) + .5);
364 <                        if (hdiffs[i] > 0) above += hdiffs[i];
365 <                        else below -= hdiffs[i];
363 >                        ceiling = Tdb*clampf(Lb(b));
364 >                        if (modhist[i] > ceiling) {
365 >                                trimmings += modhist[i] - ceiling;
366 >                                modhist[i] = ceiling;
367 >                        }
368                  }
369 <                if (above <= histot*CVRATIO)
370 <                        break;          /* close enough */
371 <                if (above-below >= 0)
372 <                        return(-1);     /* Houston, we have a problem.... */
373 <                /* original looped here as well (BEGIN_L2) */
374 <                maxd = 0;               /* find largest overvis */
375 <                for (i = 0; i < HISTRES; i++)
139 <                        if (hdiffs[i] > maxd)
140 <                                maxd = hdiffs[maxi=i];
141 <                /* broke loop here when (maxd == 0) (BREAK_L2) */
142 <                for (sd = shiftdir((maxi+.5)/HISTRES*(bwmax-bwmin)+bwmin);
143 <                                hdiffs[maxi] == maxd; sd = -sd)
144 <                        for (i = maxi+sd; i >= 0 & i < HISTRES; i += sd)
145 <                                if (hdiffs[i] < 0) {
146 <                                        if (hdiffs[i] <= -maxd) {
147 <                                                modhist[i] += maxd;
148 <                                                modhist[maxi] -= maxd;
149 <                                                hdiffs[i] += maxd;
150 <                                                hdiffs[maxi] = 0;
151 <                                        } else {
152 <                                                modhist[maxi] += hdiffs[i];
153 <                                                modhist[i] -= hdiffs[i];
154 <                                                hdiffs[maxi] += hdiffs[i];
155 <                                                hdiffs[i] = 0;
156 <                                        }
157 <                                        break;
158 <                                }
159 <                /* (END_L2) */
160 <        }
161 <        return(0);
369 >        } while (trimmings > histot*CVRATIO);
370 >
371 > #if ADJ_VEIL
372 >        mkcrfimage();                   /* contrast reduction image */
373 > #endif
374 >
375 >        return(0);                      /* we got it */
376   }
377  
378  
# Line 180 | Line 394 | int    xres;
394                          incolor = (Lw - BotMesopic) /
395                                          (TopMesopic - BotMesopic);
396                  if (incolor < 1.-FTINY) {
397 +                        b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
398 +                        if (lumf == rgblum) b /= WHTEFFICACY;
399 +                        setcolor(ctmp, b, b, b);
400                          if (incolor <= FTINY)
401                                  setcolor(scan[i], 0., 0., 0.);
402                          else
403                                  scalecolor(scan[i], incolor);
187 –                        b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
188 –                        if (lumf == rgblum) b /= WHTEFFICACY;
189 –                        setcolor(ctmp, b, b, b);
404                          addcolor(scan[i], ctmp);
405                  }
406          }
# Line 198 | Line 412 | COLOR  *scan;
412   int     xres;
413   {
414          double  mult, Lw, b;
415 <        register int    i;
415 >        register int    x;
416  
417 <        for (i = 0; i < xres; i++) {
418 <                Lw = plum(scan[i]);
417 >        for (x = 0; x < xres; x++) {
418 >                Lw = plum(scan[x]);
419                  if (Lw < LMIN) {
420 <                        setcolor(scan[i], 0., 0., 0.);
420 >                        setcolor(scan[x], 0., 0., 0.);
421                          continue;
422                  }
423 <                b = BLw(Lw);
424 <                mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
423 >                b = BLw(Lw);            /* apply brightness mapping */
424 >                mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
425                  if (lumf == rgblum) mult *= WHTEFFICACY;
426 <                scalecolor(scan[i], mult);
426 >                scalecolor(scan[x], mult);
427          }
428   }
429  
430  
431 < #ifdef DEBUG
432 < doplots()                       /* generate debugging plots */
431 > putmapping(fp)                  /* put out mapping function */
432 > FILE    *fp;
433   {
434 <        double  T, b, s;
221 <        FILE    *fp;
222 <        char    fname[128];
434 >        double  b, s;
435          register int    i;
436 +        double  wlum, sf, dlum;
437  
438 <        T = histot * (bwmax - bwmin) / HISTRES;
438 >        sf = scalef*inpexp;
439 >        if (lumf == cielum) sf *= WHTEFFICACY;
440          s = (bwmax - bwmin)/HISTRES;
441 <
442 <        sprintf(fname, "%s_hist.plt", infn);
443 <        if ((fp = fopen(fname, "w")) == NULL)
444 <                syserror(fname);
445 <        fputs("include=curve.plt\n", fp);
446 <        fputs("title=\"Brightness Frequency Distribution\"\n", fp);
447 <        fprintf(fp, "subtitle=%s\n", infn);
448 <        fputs("ymin=0\n", fp);
449 <        fputs("xlabel=\"Perceptual Brightness B(Lw)\"\n", fp);
450 <        fputs("ylabel=\"Frequency Count\"\n", fp);
237 <        fputs("Alabel=\"Histogram\"\n", fp);
238 <        fputs("Alintype=0\n", fp);
239 <        fputs("Blabel=\"Envelope\"\n", fp);
240 <        fputs("Bsymsize=0\n", fp);
241 <        fputs("Adata=\n", fp);
242 <        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
243 <                fprintf(fp, "\t%f %d\n", b, modhist[i]);
244 <        fputs(";\nBdata=\n", fp);
245 <        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
246 <                fprintf(fp, "\t%f %f\n", b, T*clampf(Lb(b)));
247 <        fputs(";\n", fp);
248 <        fclose(fp);
441 >        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
442 >                wlum = Lb(b);
443 >                if (what2do&DO_LINEAR) {
444 >                        dlum = sf*wlum;
445 >                        if (dlum > ldmax) dlum = ldmax;
446 >                        else if (dlum < ldmin) dlum = ldmin;
447 >                        fprintf(fp, "%e %e\n", wlum, dlum);
448 >                } else
449 >                        fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
450 >        }
451   }
250 – #endif

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines