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

Diff Legend

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