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.6 by greg, Wed Jun 23 12:17:09 1993 UTC vs.
Revision 2.7 by greg, Fri Jun 25 15:01:20 1993 UTC

# Line 10 | Line 10 | static char SCCSid[] = "$SunId$ LBL";
10   *     8/13/86
11   */
12  
13 < #include  <stdio.h>
13 > #include  "standard.h"
14  
15 – #include  <math.h>
16 –
15   #include  "color.h"
16  
17 < #define  FTINY          1e-7
17 > #define  TEPS           0.25    /* threshold proximity goal */
18  
19 + extern double  CHECKRAD;        /* radius over which gaussian is summed */
20 +
21   extern double  rad;             /* output pixel radius for filtering */
22  
23 + extern double  thresh;          /* maximum contribution for subpixel */
24 +
25   extern int  nrows;              /* number of rows for output */
26   extern int  ncols;              /* number of columns for output */
27  
26 – extern int  boxfilt;            /* do box filtering? */
27 –
28   extern int  xres, yres;         /* resolution of input */
29  
30   extern double  x_c, y_r;        /* conversion factors */
31  
32 < extern int  xrad;               /* x window size */
33 < extern int  yrad;               /* y window size */
32 > extern int  xrad;               /* x search radius */
33 > extern int  yrad;               /* y search radius */
34 > extern int  xbrad;              /* x box size */
35 > extern int  ybrad;              /* y box size */
36  
37   extern int  barsize;            /* size of input scan bar */
38   extern COLOR  **scanin;         /* input scan bar */
39   extern COLOR  *scanout;         /* output scan line */
40 + extern COLOR  **scoutbar;       /* output scan bar (if thresh > 0) */
41 + extern float  **greybar;        /* grey-averaged input values */
42 + extern int  obarsize;           /* size of output scan bar */
43 + extern int  orad;               /* output window radius */
44  
45   extern char  *progname;
46  
47   float  *gausstable;             /* gauss lookup table */
48  
49 + float  *ringsum;                /* sum of ring values */
50 + short  *ringwt;                 /* weight (count) of ring values */
51 + short  *ringndx;                /* ring index table */
52 + float  *warr;                   /* array of pixel weights */
53 +
54 + double  pickfilt();
55 +
56   #define  lookgauss(x)           gausstable[(int)(10.*(x)+.5)]
57  
58  
59   initmask()                      /* initialize gaussian lookup table */
60   {
61 <        extern char  *malloc();
61 >        int  gtabsiz;
62 >        double  gaussN;
63          double  d;
64          register int  x;
65  
66 <        gausstable = (float *)malloc(100*sizeof(float));
67 <        if (gausstable == NULL) {
68 <                fprintf(stderr, "%s: out of memory in initmask\n", progname);
69 <                quit(1);
70 <        }
57 <        d = x_c*y_r*0.25/rad/rad;
66 >        gtabsiz = 150*CHECKRAD;
67 >        gausstable = (float *)malloc(gtabsiz*sizeof(float));
68 >        if (gausstable == NULL)
69 >                goto memerr;
70 >        d = x_c*y_r*0.25/(rad*rad);
71          gausstable[0] = exp(-d)/sqrt(d);
72 <        for (x = 1; x < 100; x++)
72 >        for (x = 1; x < gtabsiz; x++)
73                  if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
74                          gausstable[x] = gausstable[0];
75 +        if (obarsize == 0)
76 +                return;
77 +                                        /* compute integral of filter */
78 +        gaussN = PI*sqrt(d)*exp(-d);            /* plateau */
79 +        for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
80 +                gaussN += 0.1*2.0*PI*exp(-d*d);
81 +                                        /* normalize filter */
82 +        gaussN = x_c*y_r/(rad*rad*gaussN);
83 +        for (x = 0; x < gtabsiz; x++)
84 +                gausstable[x] *= gaussN;
85 +                                        /* create ring averages table */
86 +        ringndx = (short *)malloc((2*orad*orad+1)*sizeof(short));
87 +        if (ringndx == NULL)
88 +                goto memerr;
89 +        for (x = 2*orad*orad+1; --x > orad*orad; )
90 +                ringndx[x] = -1;
91 +        do
92 +                ringndx[x] = sqrt((double)x) + 0.5;
93 +        while (x--);
94 +        ringsum = (float *)malloc((orad+1)*sizeof(float));
95 +        ringwt = (short *)malloc((orad+1)*sizeof(short));
96 +        warr = (float *)malloc(obarsize*obarsize*sizeof(float));
97 +        if (ringsum == NULL | ringwt == 0 | warr == NULL)
98 +                goto memerr;
99 +        return;
100 + memerr:
101 +        fprintf(stderr, "%s: out of memory in initmask\n", progname);
102 +        quit(1);
103   }
104  
105  
# Line 72 | Line 113 | int  c, r;
113          int  y;
114          register int  x;
115          register COLOR  *scan;
116 <
116 >        
117          wsum = 0;
118          setcolor(csum, 0.0, 0.0, 0.0);
119 <        for (y = ycent+1-yrad; y <= ycent+yrad; y++) {
119 >        for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
120                  if (y < 0) continue;
121                  if (y >= yres) break;
122                  d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
123 <                if (d > 0.5+FTINY || d < -0.5-FTINY)
124 <                        continue;
123 >                if (d < -0.5) continue;
124 >                if (d >= 0.5) break;
125                  scan = scanin[y%barsize];
126 <                for (x = xcent+1-xrad; x <= xcent+xrad; x++) {
126 >                for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
127                          if (x < 0) continue;
128                          if (x >= xres) break;
129                          d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
130 <                        if (d > 0.5+FTINY || d < -0.5-FTINY)
131 <                                continue;
130 >                        if (d < -0.5) continue;
131 >                        if (d >= 0.5) break;
132                          wsum++;
133                          addcolor(csum, scan[x]);
134                  }
# Line 127 | Line 168 | int  c, r;
168                  }
169          }
170          scalecolor(csum, 1.0/wsum);
171 + }
172 +
173 +
174 + dothresh(xcent, ycent, ccent, rcent)    /* gaussian threshold filter */
175 + int  xcent, ycent;
176 + int  ccent, rcent;
177 + {
178 +        double  d;
179 +        int  r, y;
180 +        register int  c, x;
181 +        register float  *gscan;
182 + #define pval gscan
183 +                                        /* compute ring sums */
184 +        bzero((char *)ringsum, (orad+1)*sizeof(float));
185 +        bzero((char *)ringwt, (orad+1)*sizeof(short));
186 +        for (r = -orad; r <= orad; r++) {
187 +                if (rcent+r < 0) continue;
188 +                if (rcent+r >= ncols) break;
189 +                gscan = greybar[(rcent+r)%obarsize];
190 +                for (c = -orad; c <= orad; c++) {
191 +                        if (ccent+c < 0) continue;
192 +                        if (ccent+c >= ncols) break;
193 +                        x = ringndx[c*c + r*r];
194 +                        if (x < 0) continue;
195 +                        ringsum[x] += gscan[ccent+c];
196 +                        ringwt[x]++;
197 +                }
198 +        }
199 +                                        /* filter each subpixel */
200 +        for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
201 +                if (y < 0) continue;
202 +                if (y >= yres) break;
203 +                d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent);
204 +                if (d < -0.5) continue;
205 +                if (d >= 0.5) break;
206 +                for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
207 +                        if (x < 0) continue;
208 +                        if (x >= xres) break;
209 +                        d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
210 +                        if (d < -0.5) continue;
211 +                        if (d >= 0.5) break;
212 +                        pval = scanin[y%barsize][x];
213 +                        sumans(x, y, rcent, ccent, pickfilt(bright(pval)));
214 +                }
215 +        }
216 + #undef pval
217 + }
218 +
219 +
220 + double
221 + pickfilt(p0)                    /* find filter multiplier for p0 */
222 + double  p0;
223 + {
224 +        double  m = 1.0;
225 +        double  t, avg, wsum;
226 +        int  ilimit = 5/TEPS;
227 +        register int  i;
228 +                                /* iterative search for m */
229 +        do {
230 +                                        /* compute grey weighted average */
231 +                i = 1.25*CHECKRAD*rad*m + .5;
232 +                if (i > orad) i = orad+1;
233 +                avg = wsum = 0.0;
234 +                while (i--) {
235 +                        t = i/(m*rad);
236 +                        t = lookgauss(t*t);
237 +                        avg += t*ringsum[i];
238 +                        wsum += t*ringwt[i];
239 +                }
240 +                avg /= wsum;
241 +                                        /* check threshold */
242 +                t = p0 - avg;
243 +                if (t < 0.0) t = -t;
244 +                t *= gausstable[0]/(m*m*avg);
245 +                if (t <= thresh && (m <= 1.0+FTINY ||
246 +                                (thresh-t)/thresh <= TEPS))
247 +                        break;
248 +                if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY ||
249 +                                (t-thresh)/thresh <= TEPS))
250 +                        break;
251 +                                        /* next guesstimate */
252 +                m *= sqrt(t/thresh);
253 +                if (m < 1.0) m = 1.0;
254 +                else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD;
255 +        } while (--ilimit > 0);
256 +        return(m);
257 + }
258 +
259 +
260 + sumans(px, py, rcent, ccent, m)         /* sum input pixel to output */
261 + int  px, py;
262 + int  rcent, ccent;
263 + double  m;
264 + {
265 +        double  dy, dx;
266 +        COLOR  pval, ctmp;
267 +        int  ksiz, r;
268 +        double  pc, pr, norm;
269 +        register int  i, c;
270 +        register COLOR  *scan;
271 +
272 +        copycolor(pval, scanin[py%barsize][px]);
273 +        pc = x_c*(px+.5);
274 +        pr = y_r*(py+.5);
275 +        ksiz = CHECKRAD*m*rad + 1;
276 +        if (ksiz > orad) ksiz = orad;
277 +                                                /* compute normalization */
278 +        norm = 0.0;
279 +        i = 0;
280 +        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
281 +                if (r < 0) continue;
282 +                if (r >= nrows) break;
283 +                dy = (pr - (r+.5))/(m*rad);
284 +                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
285 +                        if (c < 0) continue;
286 +                        if (c >= ncols) break;
287 +                        dx = (pc - (c+.5))/(m*rad);
288 +                        norm += warr[i++] = lookgauss(dx*dx + dy*dy);
289 +                }
290 +        }
291 +        norm = 1.0/norm;
292 +        if (x_c < 1.0) norm *= x_c;
293 +        if (y_r < 1.0) norm *= y_r;
294 +                                                /* sum pixels */
295 +        i = 0;
296 +        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
297 +                if (r < 0) continue;
298 +                if (r >= nrows) break;
299 +                scan = scoutbar[r%obarsize];
300 +                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
301 +                        if (c < 0) continue;
302 +                        if (c >= ncols) break;
303 +                        copycolor(ctmp, pval);
304 +                        dx = norm*warr[i++];
305 +                        scalecolor(ctmp, dx);
306 +                        addcolor(scan[c], ctmp);
307 +                }
308 +        }
309   }

Diff Legend

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