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.3 by greg, Fri Oct 11 10:47:00 1996 UTC vs.
Revision 3.10 by gwlarson, Thu Mar 12 15:47:34 1998 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 exp10(x)        exp(2.302585093*(x))
16 > #define LN_10           2.30258509299404568402
17 > #define exp10(x)        exp(LN_10*(x))
18  
19 < int     modhist[HISTRES];               /* modified histogram */
19 > float   modhist[HISTRES];               /* modified histogram */
20 > double  mhistot;                        /* modified histogram total */
21   float   cumf[HISTRES+1];                /* cumulative distribution function */
22  
23  
24 + getfixations(fp)                /* load fixation history list */
25 + FILE    *fp;
26 + {
27 + #define FIXHUNK         128
28 +        RESOLU  fvres;
29 +        int     pos[2];
30 +        register int    px, py, i;
31 +                                /* initialize our resolution struct */
32 +        if ((fvres.or=inpres.or)&YMAJOR) {
33 +                fvres.xr = fvxr;
34 +                fvres.yr = fvyr;
35 +        } else {
36 +                fvres.xr = fvyr;
37 +                fvres.yr = fvxr;
38 +        }
39 +                                /* read each picture position */
40 +        while (fscanf(fp, "%d %d", &pos[0], &pos[1]) == 2) {
41 +                                /* convert to closest index in foveal image */
42 +                loc2pix(pos, &fvres,
43 +                                (pos[0]+.5)/inpres.xr, (pos[1]+.5)/inpres.yr);
44 +                                /* include nine neighborhood samples */
45 +                for (px = pos[0]-1; px <= pos[0]+1; px++) {
46 +                        if (px < 0 || px >= fvxr)
47 +                                continue;
48 +                        for (py = pos[1]-1; py <= pos[1]+1; py++) {
49 +                                if (py < 0 || py >= fvyr)
50 +                                        continue;
51 +                                for (i = nfixations; i-- > 0; )
52 +                                        if (fixlst[i][0] == px &&
53 +                                                        fixlst[i][1] == py)
54 +                                                break;
55 +                                if (i >= 0)
56 +                                        continue;       /* already there */
57 +                                if (nfixations % FIXHUNK == 0) {
58 +                                        if (nfixations)
59 +                                                fixlst = (short (*)[2])
60 +                                                        realloc((char *)fixlst,
61 +                                                        (nfixations+FIXHUNK)*
62 +                                                        2*sizeof(short));
63 +                                        else
64 +                                                fixlst = (short (*)[2])malloc(
65 +                                                        FIXHUNK*2*sizeof(short)
66 +                                                        );
67 +                                        if (fixlst == NULL)
68 +                                                syserror("malloc");
69 +                                }
70 +                                fixlst[nfixations][0] = px;
71 +                                fixlst[nfixations][1] = py;
72 +                                nfixations++;
73 +                        }
74 +                }
75 +        }
76 +        if (!feof(fp)) {
77 +                fprintf(stderr, "%s: format error reading fixation data\n",
78 +                                progname);
79 +                exit(1);
80 +        }
81 + #undef  FIXHUNK
82 + }
83 +
84 +
85 + gethisto(fp)                    /* load precomputed luminance histogram */
86 + FILE    *fp;
87 + {
88 +        float   histo[MAXPREHIST];
89 +        double  histart, histep;
90 +        double  l, b, lastb, w;
91 +        int     n;
92 +        register int    i;
93 +                                        /* load data */
94 +        for (i = 0; i < MAXPREHIST &&
95 +                        fscanf(fp, "%lf %f", &b, &histo[i]) == 2; i++) {
96 +                if (i > 1 && fabs(b - lastb - histep) > 1e-3) {
97 +                        fprintf(stderr,
98 +                                "%s: uneven step size in histogram data\n",
99 +                                        progname);
100 +                        exit(1);
101 +                }
102 +                if (i == 1)
103 +                        histep = b - (histart = lastb);
104 +                lastb = b;
105 +        }
106 +        if (i < 2 || !feof(fp)) {
107 +                fprintf(stderr,
108 +                "%s: format/length error loading histogram (log10L %f at %d)\n",
109 +                                progname, b, i);
110 +                exit(1);
111 +        }
112 +        n = i;
113 +        histart *= LN_10;
114 +        histep *= LN_10;
115 +                                        /* find extrema */
116 +        for (i = 0; i < n && histo[i] <= FTINY; i++)
117 +                ;
118 +        bwmin = histart + i*histep;
119 +        for (i = n; i-- && histo[i] <= FTINY; )
120 +                ;
121 +        bwmax = histart + i*histep;
122 +        if (bwmax > Bl(LMAX))
123 +                bwmax = Bl(LMAX);
124 +        if (bwmin < Bl(LMIN))
125 +                bwmin = Bl(LMIN);
126 +        else                            /* duplicate bottom bin */
127 +                bwmin = bwmax - (bwmax-bwmin)*HISTRES/(HISTRES-1);
128 +                                        /* convert histogram */
129 +        bwavg = 0.; histot = 0.;
130 +        for (i = 0; i < HISTRES; i++)
131 +                bwhist[i] = 0.;
132 +        for (i = 0, b = histart; i < n; i++, b += histep) {
133 +                if (b < bwmin) continue;
134 +                if (b > bwmax) break;
135 +                w = histo[i];
136 +                bwavg += w*b;
137 +                bwhist[bwhi(b)] += w;
138 +                histot += w;
139 +        }
140 +        bwavg /= histot;
141 +        if (bwmin > Bl(LMIN)+FTINY) {   /* add false samples at bottom */
142 +                bwhist[1] *= 0.5;
143 +                bwhist[0] += bwhist[1];
144 +        }
145 + }
146 +
147 +
148 + double
149 + centprob(x, y)                  /* center-weighting probability function */
150 + int     x, y;
151 + {
152 +        double  xr, yr, p;
153 +                                /* paraboloid, 0 at 90 degrees from center */
154 +        xr = (x - .5*(fvxr-1))/90.;     /* 180 degree fisheye has fv?r == 90 */
155 +        yr = (y - .5*(fvyr-1))/90.;
156 +        p = 1. - xr*xr - yr*yr;
157 +        return(p < 0. ? 0. : p);
158 + }
159 +
160 +
161 + comphist()                      /* create foveal sampling histogram */
162 + {
163 +        double  l, b, w, lwmin, lwmax;
164 +        register int    x, y;
165 +                                        /* check for precalculated histogram */
166 +        if (what2do&DO_PREHIST)
167 +                return;
168 +        lwmin = 1e10;                   /* find extrema */
169 +        lwmax = 0.;
170 +        for (y = 0; y < fvyr; y++)
171 +                for (x = 0; x < fvxr; x++) {
172 +                        l = plum(fovscan(y)[x]);
173 +                        if (l < lwmin) lwmin = l;
174 +                        if (l > lwmax) lwmax = l;
175 +                }
176 +        lwmax *= 1.01;
177 +        if (lwmax > LMAX)
178 +                lwmax = LMAX;
179 +        bwmax = Bl(lwmax);
180 +        if (lwmin < LMIN) {
181 +                lwmin = LMIN;
182 +                bwmin = Bl(LMIN);
183 +        } else {                        /* duplicate bottom bin */
184 +                bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
185 +                lwmin = Lb(bwmin);
186 +        }
187 +                                        /* (re)compute histogram */
188 +        bwavg = 0.;
189 +        histot = 0.;
190 +        for (x = 0; x < HISTRES; x++)
191 +                bwhist[x] = 0.;
192 +                                        /* global average */
193 +        if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
194 +                for (y = 0; y < fvyr; y++)
195 +                        for (x = 0; x < fvxr; x++) {
196 +                                l = plum(fovscan(y)[x]);
197 +                                if (l < lwmin) continue;
198 +                                if (l > lwmax) continue;
199 +                                b = Bl(l);
200 +                                w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
201 +                                bwavg += w*b;
202 +                                bwhist[bwhi(b)] += w;
203 +                                histot += w;
204 +                        }
205 +                                        /* average fixation points */
206 +        if (what2do&DO_FIXHIST && nfixations > 0) {
207 +                if (histot > FTINY)
208 +                        w = fixfrac/(1.-fixfrac)*histot/nfixations;
209 +                else
210 +                        w = 1.;
211 +                for (x = 0; x < nfixations; x++) {
212 +                        l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
213 +                        if (l < lwmin) continue;
214 +                        if (l > lwmax) continue;
215 +                        b = Bl(l);
216 +                        bwavg += w*b;
217 +                        bwhist[bwhi(b)] += w;
218 +                        histot += w;
219 +                }
220 +        }
221 +        bwavg /= histot;
222 +        if (lwmin > LMIN+FTINY) {       /* add false samples at bottom */
223 +                bwhist[1] *= 0.5;
224 +                bwhist[0] += bwhist[1];
225 +        }
226 + }
227 +
228 +
229   mkcumf()                        /* make cumulative distribution function */
230   {
231          register int    i;
232 <        register long   sum;
232 >        register double sum;
233  
234 <        cumf[0] = 0.;
235 <        sum = modhist[0];
236 <        for (i = 1; i < HISTRES; i++) {
237 <                cumf[i] = (double)sum/histot;
234 >        mhistot = 0.;           /* compute modified total */
235 >        for (i = 0; i < HISTRES; i++)
236 >                mhistot += modhist[i];
237 >
238 >        sum = 0.;               /* compute cumulative function */
239 >        for (i = 0; i < HISTRES; i++) {
240 >                cumf[i] = sum/mhistot;
241                  sum += modhist[i];
242          }
243          cumf[HISTRES] = 1.;
# Line 57 | Line 267 | double Lw;
267                  return(Bldmin);
268          if (b >= bwmax-FTINY)
269                  return(Bldmax);
270 <        return(Bldmin + cf(Bl(Lw))*(Bldmax-Bldmin));
270 >        return(Bldmin + cf(b)*(Bldmax-Bldmin));
271   }
272  
273  
# Line 96 | Line 306 | double Lw;
306  
307  
308   int
99 – shiftdir(bw)            /* compute shift direction for histogram */
100 – double  bw;
101 – {
102 –        if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
103 –                return(1);
104 –        return(-1);
105 – }
106 –
107 –
108 – int
309   mkbrmap()                       /* make dynamic range map */
310   {
111 –        int     hdiffs[HISTRES], above, below;
311          double  T, b, s;
312 <        int     maxd, maxi, sd;
312 >        double  ceiling, trimmings;
313          register int    i;
314                                          /* copy initial histogram */
315 <        for (i = 0; i < HISTRES; i++)
117 <                modhist[i] = bwhist[i];
118 <        T = histot * (bwmax - bwmin) / HISTRES;
315 >        bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
316          s = (bwmax - bwmin)/HISTRES;
317                                          /* loop until satisfactory */
318 <        for ( ; ; ) {
319 <                mkcumf();               /* sync brightness mapping */
320 <                above = below = 0;      /* compute visibility overflow */
318 >        do {
319 >                mkcumf();                       /* sync brightness mapping */
320 >                if (mhistot <= histot*CVRATIO)
321 >                        return(-1);             /* no compression needed! */
322 >                T = mhistot * (bwmax - bwmin) / HISTRES;
323 >                trimmings = 0.;                 /* clip to envelope */
324                  for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
325 <                        hdiffs[i] = modhist[i] - (int)(T*clampf(Lb(b)) + .5);
326 <                        if (hdiffs[i] > 0) above += hdiffs[i];
327 <                        else below -= hdiffs[i];
325 >                        ceiling = T*clampf(Lb(b));
326 >                        if (modhist[i] > ceiling) {
327 >                                trimmings += modhist[i] - ceiling;
328 >                                modhist[i] = ceiling;
329 >                        }
330                  }
331 <                if (above <= histot*CVRATIO)
332 <                        break;          /* close enough */
333 <                if (above-below >= 0)
132 <                        return(-1);     /* Houston, we have a problem.... */
133 <                /* original looped here as well (BEGIN_L2) */
134 <                maxd = 0;               /* find largest overvis */
135 <                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);
331 >        } while (trimmings > histot*CVRATIO);
332 >
333 >        return(0);                      /* we got it */
334   }
335  
336  
# Line 211 | Line 386 | int    xres;
386   }
387  
388  
389 < #ifdef DEBUG
390 < doplots()                       /* generate debugging plots */
389 > putmapping(fp)                  /* put out mapping function */
390 > FILE    *fp;
391   {
392 <        double  T, b, s;
218 <        FILE    *fp;
219 <        char    fname[128];
392 >        double  b, s;
393          register int    i;
394 +        double  wlum, sf, dlum;
395  
396 <        T = histot * (bwmax - bwmin) / HISTRES;
396 >        sf = scalef*inpexp;
397 >        if (lumf == cielum) sf *= WHTEFFICACY;
398          s = (bwmax - bwmin)/HISTRES;
399 <
400 <        sprintf(fname, "%s_hist.plt", infn);
401 <        if ((fp = fopen(fname, "w")) == NULL)
402 <                syserror(fname);
403 <        fputs("include=curve.plt\n", fp);
404 <        fputs("title=\"Brightness Frequency Distribution\"\n", fp);
405 <        fprintf(fp, "subtitle=%s\n", infn);
406 <        fputs("ymin=0\n", fp);
407 <        fputs("xlabel=\"Perceptual Brightness B(Lw)\"\n", fp);
408 <        fputs("ylabel=\"Frequency Count\"\n", fp);
234 <        fputs("Alabel=\"Histogram\"\n", fp);
235 <        fputs("Alintype=0\n", fp);
236 <        fputs("Blabel=\"Envelope\"\n", fp);
237 <        fputs("Bsymsize=0\n", fp);
238 <        fputs("Adata=\n", fp);
239 <        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
240 <                fprintf(fp, "\t%f %d\n", b, modhist[i]);
241 <        fputs(";\nBdata=\n", fp);
242 <        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
243 <                fprintf(fp, "\t%f %f\n", b, T*clampf(Lb(b)));
244 <        fputs(";\n", fp);
245 <        fclose(fp);
399 >        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
400 >                wlum = Lb(b);
401 >                if (what2do&DO_LINEAR) {
402 >                        dlum = sf*wlum;
403 >                        if (dlum > ldmax) dlum = ldmax;
404 >                        else if (dlum < ldmin) dlum = ldmin;
405 >                        fprintf(fp, "%e %e\n", wlum, dlum);
406 >                } else
407 >                        fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
408 >        }
409   }
247 – #endif

Diff Legend

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