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.8 by greg, Wed Jan 29 13:22:13 1997 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 1996 Regents of the University of California */
1 > /* Copyright (c) 1997 Regents of the University of California */
2  
3   #ifndef lint
4   static char SCCSid[] = "$SunId$ LBL";
# Line 11 | Line 11 | static char SCCSid[] = "$SunId$ LBL";
11   #include "pcond.h"
12  
13  
14 < #define CVRATIO         0.025           /* fraction of pixels allowed > env. */
14 > #define CVRATIO         0.025           /* fraction of samples allowed > env. */
15  
16 #define BotMesopic      5.62e-3         /* top of scotopic range */
17 #define TopMesopic      5.62            /* bottom of photopic range */
18
16   #define exp10(x)        exp(2.302585093*(x))
17  
18 < int     modhist[HISTRES];               /* modified histogram */
18 > float   modhist[HISTRES];               /* modified histogram */
19 > double  mhistot;                        /* modified histogram total */
20   float   cumf[HISTRES+1];                /* cumulative distribution function */
21  
22  
23 + 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 +        if ((fvres.or=inpres.or)&YMAJOR) {
32 +                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 +                                                        realloc((char *)fixlst,
60 +                                                        (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 + double
85 + centprob(x, y)                  /* center-weighting probability function */
86 + int     x, y;
87 + {
88 +        double  xr, yr, p;
89 +                                /* paraboloid, 0 at 90 degrees from center */
90 +        xr = (x - .5*(fvxr-1))/90.;     /* 180 degree fisheye has fv?r == 90 */
91 +        yr = (y - .5*(fvyr-1))/90.;
92 +        p = 1. - xr*xr - yr*yr;
93 +        return(p < 0. ? 0. : p);
94 + }
95 +
96 +
97 + comphist()                      /* create foveal sampling histogram */
98 + {
99 +        double  l, b, w, lwmin, lwmax;
100 +        register int    x, y;
101 +
102 +        lwmin = 1e10;                   /* find extrema */
103 +        lwmax = 0.;
104 +        for (y = 0; y < fvyr; y++)
105 +                for (x = 0; x < fvxr; x++) {
106 +                        l = plum(fovscan(y)[x]);
107 +                        if (l < lwmin) lwmin = l;
108 +                        if (l > lwmax) lwmax = l;
109 +                }
110 +        lwmin -= FTINY;
111 +        lwmax += FTINY;
112 +        if (lwmin < LMIN) lwmin = LMIN;
113 +        if (lwmax > LMAX) lwmax = LMAX;
114 +        bwmin = Bl(lwmin);
115 +        bwmax = Bl(lwmax);
116 +                                        /* (re)compute histogram */
117 +        bwavg = 0.;
118 +        histot = 0.;
119 +        for (x = 0; x < HISTRES; x++)
120 +                bwhist[x] = 0.;
121 +                                        /* global average */
122 +        if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
123 +                for (y = 0; y < fvyr; y++)
124 +                        for (x = 0; x < fvxr; x++) {
125 +                                l = plum(fovscan(y)[x]);
126 +                                if (l < lwmin) continue;
127 +                                if (l > lwmax) continue;
128 +                                b = Bl(l);
129 +                                bwavg += b;
130 +                                w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
131 +                                bwhist[bwhi(b)] += w;
132 +                                histot += w;
133 +                        }
134 +                                        /* average fixation points */
135 +        if (what2do&DO_FIXHIST && nfixations > 0) {
136 +                if (histot > FTINY)
137 +                        w = fixfrac/(1.-fixfrac)*histot/nfixations;
138 +                else
139 +                        w = 1.;
140 +                for (x = 0; x < nfixations; x++) {
141 +                        l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
142 +                        if (l < lwmin) continue;
143 +                        if (l > lwmax) continue;
144 +                        b = Bl(l);
145 +                        bwavg += b;
146 +                        bwhist[bwhi(b)] += w;
147 +                        histot += w;
148 +                }
149 +        }
150 +        bwavg /= histot;
151 + }
152 +
153 +
154   mkcumf()                        /* make cumulative distribution function */
155   {
156          register int    i;
157 <        register long   sum;
157 >        register double sum;
158  
159 <        cumf[0] = 0.;
160 <        sum = modhist[0];
161 <        for (i = 1; i < HISTRES; i++) {
162 <                cumf[i] = (double)sum/histot;
159 >        mhistot = 0.;           /* compute modified total */
160 >        for (i = 0; i < HISTRES; i++)
161 >                mhistot += modhist[i];
162 >
163 >        sum = 0.;               /* compute cumulative function */
164 >        for (i = 0; i < HISTRES; i++) {
165 >                cumf[i] = sum/mhistot;
166                  sum += modhist[i];
167          }
168          cumf[HISTRES] = 1.;
# Line 60 | Line 192 | double Lw;
192                  return(Bldmin);
193          if (b >= bwmax-FTINY)
194                  return(Bldmax);
195 <        return(Bldmin + cf(Bl(Lw))*(Bldmax-Bldmin));
195 >        return(Bldmin + cf(b)*(Bldmax-Bldmin));
196   }
197  
198  
# Line 99 | Line 231 | double Lw;
231  
232  
233   int
102 shiftdir(bw)            /* compute shift direction for histogram */
103 double  bw;
104 {
105        if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
106                return(1);
107        return(-1);
108 }
109
110
111 int
234   mkbrmap()                       /* make dynamic range map */
235   {
114        int     hdiffs[HISTRES], above, below;
236          double  T, b, s;
237 <        int     maxd, maxi, sd;
237 >        double  ceiling, trimmings;
238          register int    i;
239                                          /* copy initial histogram */
240 <        for (i = 0; i < HISTRES; i++)
120 <                modhist[i] = bwhist[i];
121 <        T = histot * (bwmax - bwmin) / HISTRES;
240 >        bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
241          s = (bwmax - bwmin)/HISTRES;
242                                          /* loop until satisfactory */
243 <        for ( ; ; ) {
244 <                mkcumf();               /* sync brightness mapping */
245 <                above = below = 0;      /* compute visibility overflow */
243 >        do {
244 >                mkcumf();                       /* sync brightness mapping */
245 >                if (mhistot <= histot*CVRATIO)
246 >                        return(-1);             /* no compression needed! */
247 >                T = mhistot * (bwmax - bwmin) / HISTRES;
248 >                trimmings = 0.;                 /* clip to envelope */
249                  for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
250 <                        hdiffs[i] = modhist[i] - (int)(T*clampf(Lb(b)) + .5);
251 <                        if (hdiffs[i] > 0) above += hdiffs[i];
252 <                        else below -= hdiffs[i];
250 >                        ceiling = T*clampf(Lb(b));
251 >                        if (modhist[i] > ceiling) {
252 >                                trimmings += modhist[i] - ceiling;
253 >                                modhist[i] = ceiling;
254 >                        }
255                  }
256 <                if (above <= histot*CVRATIO)
257 <                        break;          /* close enough */
258 <                if (above-below >= 0)
135 <                        return(-1);     /* Houston, we have a problem.... */
136 <                /* original looped here as well (BEGIN_L2) */
137 <                maxd = 0;               /* find largest overvis */
138 <                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);
256 >        } while (trimmings > histot*CVRATIO);
257 >
258 >        return(0);                      /* we got it */
259   }
260  
261  
# Line 180 | Line 277 | int    xres;
277                          incolor = (Lw - BotMesopic) /
278                                          (TopMesopic - BotMesopic);
279                  if (incolor < 1.-FTINY) {
280 +                        b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
281 +                        if (lumf == rgblum) b /= WHTEFFICACY;
282 +                        setcolor(ctmp, b, b, b);
283                          if (incolor <= FTINY)
284                                  setcolor(scan[i], 0., 0., 0.);
285                          else
286                                  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);
287                          addcolor(scan[i], ctmp);
288                  }
289          }
# Line 214 | Line 311 | int    xres;
311   }
312  
313  
314 < #ifdef DEBUG
315 < doplots()                       /* generate debugging plots */
314 > putmapping(fp)                  /* put out mapping function */
315 > FILE    *fp;
316   {
317 <        double  T, b, s;
221 <        FILE    *fp;
222 <        char    fname[128];
317 >        double  b, s;
318          register int    i;
319 +        double  wlum, sf, dlum;
320  
321 <        T = histot * (bwmax - bwmin) / HISTRES;
321 >        sf = scalef*inpexp;
322 >        if (lumf == cielum) sf *= WHTEFFICACY;
323          s = (bwmax - bwmin)/HISTRES;
324 <
325 <        sprintf(fname, "%s_hist.plt", infn);
326 <        if ((fp = fopen(fname, "w")) == NULL)
327 <                syserror(fname);
328 <        fputs("include=curve.plt\n", fp);
329 <        fputs("title=\"Brightness Frequency Distribution\"\n", fp);
330 <        fprintf(fp, "subtitle=%s\n", infn);
331 <        fputs("ymin=0\n", fp);
332 <        fputs("xlabel=\"Perceptual Brightness B(Lw)\"\n", fp);
333 <        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);
324 >        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
325 >                wlum = Lb(b);
326 >                if (what2do&DO_LINEAR) {
327 >                        dlum = sf*wlum;
328 >                        if (dlum > ldmax) dlum = ldmax;
329 >                        else if (dlum < ldmin) dlum = ldmin;
330 >                        fprintf(fp, "%e %e\n", wlum, dlum);
331 >                } else
332 >                        fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
333 >        }
334   }
250 #endif

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines