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.12 by greg, Thu Apr 4 13:01:35 1996 UTC vs.
Revision 2.16 by schorsch, Mon Jun 30 14:59:12 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1993 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  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  
# Line 68 | 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 124 | 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];
# Line 132 | Line 134 | int  c, r;
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 : (double)(x - xcent);
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+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 174 | Line 178 | int  c, r;
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 187 | Line 192 | int  ccent, rcent;
192          register int  c, x;
193          register float  *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 >= nrows) break;
# Line 208 | 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                          offs = x < 0 ? xres : x >= xres ? -xres : 0;
221                          if (offs && !wrapfilt)
222                                  continue;
223 <                        d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
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                          sumans(x, y, rcent, ccent,
# Line 232 | Line 237 | double  p0;
237          double  m = 1.0;
238          double  t, num, denom, avg, wsum;
239          double  mlimit[2];
240 <        int  ilimit = 4/TEPS;
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;
246 >                i = RSCA*CHECKRAD*rad*m + .5;
247                  if (i > orad) i = orad;
248                  avg = wsum = 0.0;
249                  while (i--) {
# Line 290 | 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, offs;
301          double  pc, pr, norm;
# Line 312 | 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 (!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;

Diff Legend

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