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

Comparing ray/src/px/pf3.c (file contents):
Revision 2.9 by greg, Fri Jul 9 10:27:52 1993 UTC vs.
Revision 2.10 by greg, Thu Jul 15 18:32:06 1993 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 1992 Regents of the University of California */
1 > /* Copyright (c) 1993 Regents of the University of California */
2  
3   #ifndef lint
4   static char SCCSid[] = "$SunId$ LBL";
# Line 15 | Line 15 | static char SCCSid[] = "$SunId$ LBL";
15   #include  "color.h"
16  
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 222 | Line 223 | pickfilt(p0)                   /* find filter multiplier for p0 */
223   double  p0;
224   {
225          double  m = 1.0;
226 <        double  t, avg, wsum;
226 >        double  t, num, denom, avg, wsum;
227 >        double  maxm = orad/rad/CHECKRAD;
228 >        double  mlimit[2];
229          int  ilimit = 4/TEPS;
230          register int  i;
231                                  /* iterative search for m */
232 +        mlimit[0] = 1.0; mlimit[1] = maxm;
233          do {
234                                          /* compute grey weighted average */
235                  i = 1.25*CHECKRAD*rad*m + .5;
# Line 237 | Line 241 | double  p0;
241                          avg += t*ringsum[i];
242                          wsum += t*ringwt[i];
243                  }
244 +                if (avg < 1e-20)                /* zero inclusive average */
245 +                        return(1.0);
246                  avg /= wsum;
247                                          /* check threshold */
248 <                t = p0 - avg;
249 <                if (t < 0.0) t = -t;
250 <                t *= gausstable[0]/(m*m*avg);
251 <                if (t <= thresh && (m <= 1.0+FTINY ||
252 <                                (thresh-t)/thresh <= TEPS))
253 <                        break;
254 <                if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY ||
255 <                                (t-thresh)/thresh <= TEPS))
256 <                        break;
248 >                denom = m*m/gausstable[0] - p0/avg;
249 >                if (denom <= FTINY)             /* zero exclusive average */
250 >                        return(maxm);
251 >                num = p0/avg - 1.0;
252 >                if (num < 0.0) num = -num;
253 >                t = num/denom;
254 >                if (t <= thresh) {
255 >                        if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS)
256 >                                break;
257 >                } else {
258 >                        if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS)
259 >                                break;
260 >                }
261 >                t = m;                  /* remember current m */
262                                          /* next guesstimate */
263 <                m *= sqrt(t/thresh);
264 <                if (m < 1.0) m = 1.0;
265 <                else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD;
263 >                m = sqrt(gausstable[0]*(num/thresh + p0/avg));
264 >                if (m < t) {            /* bound it */
265 >                        if (m <= mlimit[0]+FTINY)
266 >                                m = 0.5*(mlimit[0] + t);
267 >                        mlimit[1] = t;
268 >                } else {
269 >                        if (m >= mlimit[1]-FTINY)
270 >                                m = 0.5*(mlimit[1] + t);
271 >                        mlimit[0] = t;
272 >                }
273          } while (--ilimit > 0);
274          return(m);
275   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines (old)
> Changed lines (new)