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

Comparing ray/src/px/pf3.c (file contents):
Revision 2.8 by greg, Fri Jun 25 17:08:40 1993 UTC vs.
Revision 2.15 by greg, Tue Feb 25 00:26:05 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1992 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   *  pf3.c - routines for gaussian and box filtering
6   *
# Line 14 | Line 11 | static char SCCSid[] = "$SunId$ LBL";
11  
12   #include  "color.h"
13  
14 < #define  TEPS           0.25    /* threshold proximity goal */
14 > #define  RSCA           1.13    /* square-radius multiplier: sqrt(4/PI) */
15 > #define  TEPS           0.2     /* threshold proximity goal */
16 > #define  REPS           0.1     /* radius proximity goal */
17  
18   extern double  CHECKRAD;        /* radius over which gaussian is summed */
19  
# Line 42 | Line 41 | extern float  **greybar;       /* grey-averaged input values
41   extern int  obarsize;           /* size of output scan bar */
42   extern int  orad;               /* output window radius */
43  
44 + extern int  wrapfilt;           /* wrap filter horizontally? */
45 +
46   extern char  *progname;
47  
48   float  *gausstable;             /* gauss lookup table */
# Line 51 | Line 52 | short  *ringwt;                        /* weight (count) of ring values */
52   short  *ringndx;                /* ring index table */
53   float  *warr;                   /* array of pixel weights */
54  
55 + extern double  (*ourbright)();  /* brightness computation function */
56 +
57   double  pickfilt();
58  
59   #define  lookgauss(x)           gausstable[(int)(10.*(x)+.5)]
# Line 63 | Line 66 | initmask()                     /* initialize gaussian lookup table */
66          double  d;
67          register int  x;
68  
69 <        gtabsiz = 150*CHECKRAD;
69 >        gtabsiz = 111*CHECKRAD*CHECKRAD;
70          gausstable = (float *)malloc(gtabsiz*sizeof(float));
71          if (gausstable == NULL)
72                  goto memerr;
73          d = x_c*y_r*0.25/(rad*rad);
74 <        gausstable[0] = exp(-d)/sqrt(d);
74 >        gausstable[0] = exp(-d);
75          for (x = 1; x < gtabsiz; x++)
76 <                if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
76 >                if (x*0.1 <= d)
77                          gausstable[x] = gausstable[0];
78 +                else
79 +                        gausstable[x] = exp(-x*0.1);
80          if (obarsize == 0)
81                  return;
82                                          /* compute integral of filter */
83 <        gaussN = PI*sqrt(d)*exp(-d);            /* plateau */
84 <        for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
85 <                gaussN += 0.1*2.0*PI*exp(-d*d);
83 >        gaussN = PI*d*exp(-d);                  /* plateau */
84 >        for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1)
85 >                gaussN += 0.1*2.0*PI*d*exp(-d*d);
86                                          /* normalize filter */
87          gaussN = x_c*y_r/(rad*rad*gaussN);
88          for (x = 0; x < gtabsiz; x++)
# Line 111 | Line 116 | int  c, r;
116          int  wsum;
117          double  d;
118          int  y;
119 <        register int  x;
119 >        register int  x, offs;
120          register COLOR  *scan;
121          
122          wsum = 0;
# Line 119 | Line 124 | int  c, r;
124          for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
125                  if (y < 0) continue;
126                  if (y >= yres) break;
127 <                d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
127 >                d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent);
128                  if (d < -0.5) continue;
129                  if (d >= 0.5) break;
130                  scan = scanin[y%barsize];
131                  for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
132 <                        if (x < 0) continue;
133 <                        if (x >= xres) break;
134 <                        d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
132 >                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
133 >                        if (offs && !wrapfilt)
134 >                                continue;
135 >                        d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent);
136                          if (d < -0.5) continue;
137                          if (d >= 0.5) break;
138                          wsum++;
139 <                        addcolor(csum, scan[x]);
139 >                        addcolor(csum, scan[x+offs]);
140                  }
141          }
142 <        if (wsum > 1)
143 <                scalecolor(csum, 1.0/wsum);
142 >        if (wsum > 1) {
143 >                d = 1.0/wsum;
144 >                scalecolor(csum, d);
145 >        }
146   }
147  
148  
# Line 146 | Line 154 | int  c, r;
154          double  dy, dx, weight, wsum;
155          COLOR  ctmp;
156          int  y;
157 <        register int  x;
157 >        register int  x, offs;
158          register COLOR  *scan;
159  
160          wsum = FTINY;
# Line 157 | Line 165 | int  c, r;
165                  dy = (y_r*(y+.5) - (r+.5))/rad;
166                  scan = scanin[y%barsize];
167                  for (x = xcent-xrad; x <= xcent+xrad; x++) {
168 <                        if (x < 0) continue;
169 <                        if (x >= xres) break;
168 >                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
169 >                        if (offs && !wrapfilt)
170 >                                continue;
171                          dx = (x_c*(x+.5) - (c+.5))/rad;
172                          weight = lookgauss(dx*dx + dy*dy);
173                          wsum += weight;
174 <                        copycolor(ctmp, scan[x]);
174 >                        copycolor(ctmp, scan[x+offs]);
175                          scalecolor(ctmp, weight);
176                          addcolor(csum, ctmp);
177                  }
178          }
179 <        scalecolor(csum, 1.0/wsum);
179 >        weight = 1.0/wsum;
180 >        scalecolor(csum, weight);
181   }
182  
183  
# Line 176 | Line 186 | int  xcent, ycent;
186   int  ccent, rcent;
187   {
188          double  d;
189 <        int  r, y;
189 >        int  r, y, offs;
190          register int  c, x;
191          register float  *gscan;
182 – #define pval gscan
192                                          /* compute ring sums */
193          bzero((char *)ringsum, (orad+1)*sizeof(float));
194          bzero((char *)ringwt, (orad+1)*sizeof(short));
# Line 188 | Line 197 | int  ccent, rcent;
197                  if (rcent+r >= nrows) break;
198                  gscan = greybar[(rcent+r)%obarsize];
199                  for (c = -orad; c <= orad; c++) {
200 <                        if (ccent+c < 0) continue;
201 <                        if (ccent+c >= ncols) break;
200 >                        offs = ccent+c < 0 ? ncols :
201 >                                        ccent+c >= ncols ? -ncols : 0;
202 >                        if (offs && !wrapfilt)
203 >                                continue;
204                          x = ringndx[c*c + r*r];
205                          if (x < 0) continue;
206 <                        ringsum[x] += gscan[ccent+c];
206 >                        ringsum[x] += gscan[ccent+c+offs];
207                          ringwt[x]++;
208                  }
209          }
# Line 200 | Line 211 | int  ccent, rcent;
211          for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
212                  if (y < 0) continue;
213                  if (y >= yres) break;
214 <                d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent);
214 >                d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent);
215                  if (d < -0.5) continue;
216                  if (d >= 0.5) break;
217                  for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
218 <                        if (x < 0) continue;
219 <                        if (x >= xres) break;
220 <                        d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
218 >                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
219 >                        if (offs && !wrapfilt)
220 >                                continue;
221 >                        d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
222                          if (d < -0.5) continue;
223                          if (d >= 0.5) break;
224 <                        pval = scanin[y%barsize][x];
225 <                        sumans(x, y, rcent, ccent, pickfilt(bright(pval)));
224 >                        sumans(x, y, rcent, ccent,
225 >                        pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
226                  }
227          }
216 – #undef pval
228   }
229  
230  
# Line 222 | Line 233 | pickfilt(p0)                   /* find filter multiplier for p0 */
233   double  p0;
234   {
235          double  m = 1.0;
236 <        double  t, avg, wsum;
237 <        int  ilimit = 5/TEPS;
236 >        double  t, num, denom, avg, wsum;
237 >        double  mlimit[2];
238 >        int  ilimit = 4.0/TEPS;
239          register int  i;
240                                  /* iterative search for m */
241 +        mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
242          do {
243                                          /* compute grey weighted average */
244 <                i = 1.25*CHECKRAD*rad*m + .5;
244 >                i = RSCA*CHECKRAD*rad*m + .5;
245                  if (i > orad) i = orad;
246                  avg = wsum = 0.0;
247                  while (i--) {
# Line 237 | Line 250 | double  p0;
250                          avg += t*ringsum[i];
251                          wsum += t*ringwt[i];
252                  }
253 +                if (avg < 1e-20)                /* zero inclusive average */
254 +                        return(1.0);
255                  avg /= wsum;
256                                          /* check threshold */
257 <                t = p0 - avg;
258 <                if (t < 0.0) t = -t;
259 <                t *= gausstable[0]/(m*m*avg);
260 <                if (t <= thresh && (m <= 1.0+FTINY ||
261 <                                (thresh-t)/thresh <= TEPS))
262 <                        break;
263 <                if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY ||
264 <                                (t-thresh)/thresh <= TEPS))
265 <                        break;
257 >                denom = m*m/gausstable[0] - p0/avg;
258 >                if (denom <= FTINY) {           /* zero exclusive average */
259 >                        if (m >= mlimit[1]-REPS)
260 >                                break;
261 >                        m = mlimit[1];
262 >                        continue;
263 >                }
264 >                num = p0/avg - 1.0;
265 >                if (num < 0.0) num = -num;
266 >                t = num/denom;
267 >                if (t <= thresh) {
268 >                        if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS)
269 >                                break;
270 >                } else {
271 >                        if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS)
272 >                                break;
273 >                }
274 >                t = m;                  /* remember current m */
275                                          /* next guesstimate */
276 <                m *= sqrt(t/thresh);
277 <                if (m < 1.0) m = 1.0;
278 <                else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD;
276 >                m = sqrt(gausstable[0]*(num/thresh + p0/avg));
277 >                if (m < t) {            /* bound it */
278 >                        if (m <= mlimit[0]+FTINY)
279 >                                m = 0.5*(mlimit[0] + t);
280 >                        mlimit[1] = t;
281 >                } else {
282 >                        if (m >= mlimit[1]-FTINY)
283 >                                m = 0.5*(mlimit[1] + t);
284 >                        mlimit[0] = t;
285 >                }
286          } while (--ilimit > 0);
287          return(m);
288   }
# Line 262 | Line 293 | int  px, py;
293   int  rcent, ccent;
294   double  m;
295   {
296 <        double  dy, dx;
296 >        double  dy2, dx;
297          COLOR  pval, ctmp;
298 <        int  ksiz, r;
298 >        int  ksiz, r, offs;
299          double  pc, pr, norm;
300          register int  i, c;
301          register COLOR  *scan;
302 <
303 <        copycolor(pval, scanin[py%barsize][px]);
302 >        /*
303 >         * This normalization method fails at the picture borders because
304 >         * a different number of input pixels contribute there.
305 >         */
306 >        scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0);
307 >        copycolor(pval, scan[px]);
308          pc = x_c*(px+.5);
309          pr = y_r*(py+.5);
310          ksiz = CHECKRAD*m*rad + 1;
# Line 280 | Line 315 | double  m;
315          for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
316                  if (r < 0) continue;
317                  if (r >= nrows) break;
318 <                dy = (pr - (r+.5))/(m*rad);
318 >                dy2 = (pr - (r+.5))/(m*rad);
319 >                dy2 *= dy2;
320                  for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
321 <                        if (c < 0) continue;
322 <                        if (c >= ncols) break;
321 >                        if (!wrapfilt) {
322 >                                if (c < 0) continue;
323 >                                if (c >= ncols) break;
324 >                        }
325                          dx = (pc - (c+.5))/(m*rad);
326 <                        norm += warr[i++] = lookgauss(dx*dx + dy*dy);
326 >                        norm += warr[i++] = lookgauss(dx*dx + dy2);
327                  }
328          }
329          norm = 1.0/norm;
# Line 298 | Line 336 | double  m;
336                  if (r >= nrows) break;
337                  scan = scoutbar[r%obarsize];
338                  for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
339 <                        if (c < 0) continue;
340 <                        if (c >= ncols) break;
339 >                        offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
340 >                        if (offs && !wrapfilt)
341 >                                continue;
342                          copycolor(ctmp, pval);
343                          dx = norm*warr[i++];
344                          scalecolor(ctmp, dx);
345 <                        addcolor(scan[c], ctmp);
345 >                        addcolor(scan[c+offs], ctmp);
346                  }
347          }
348   }

Diff Legend

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