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.3 by greg, Fri Oct 2 16:23:41 1992 UTC vs.
Revision 2.14 by greg, Sat Feb 22 02:07:27 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   *
7   *     8/13/86
8   */
9  
10 < #include  <stdio.h>
10 > #include  "standard.h"
11  
15 #include  <math.h>
16
12   #include  "color.h"
13  
14 < #define  FTINY          1e-7
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 +
20   extern double  rad;             /* output pixel radius for filtering */
21  
22 + extern double  thresh;          /* maximum contribution for subpixel */
23 +
24   extern int  nrows;              /* number of rows for output */
25   extern int  ncols;              /* number of columns for output */
26  
26 extern int  boxfilt;            /* do box filtering? */
27
27   extern int  xres, yres;         /* resolution of input */
28  
29   extern double  x_c, y_r;        /* conversion factors */
30  
31 < extern int  xrad;               /* x window size */
32 < extern int  yrad;               /* y window size */
31 > extern int  xrad;               /* x search radius */
32 > extern int  yrad;               /* y search radius */
33 > extern int  xbrad;              /* x box size */
34 > extern int  ybrad;              /* y box size */
35  
36   extern int  barsize;            /* size of input scan bar */
37   extern COLOR  **scanin;         /* input scan bar */
38   extern COLOR  *scanout;         /* output scan line */
39 + extern COLOR  **scoutbar;       /* output scan bar (if thresh > 0) */
40 + 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  *exptable;               /* exponent table */
48 > float  *gausstable;             /* gauss lookup table */
49  
50 < #define  lookexp(x)             exptable[(int)(-10.*(x)+.5)]
50 > float  *ringsum;                /* sum of ring values */
51 > 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)]
60 +
61 +
62   initmask()                      /* initialize gaussian lookup table */
63   {
64 <        extern char  *malloc();
64 >        int  gtabsiz;
65 >        double  gaussN;
66 >        double  d;
67          register int  x;
68  
69 <        exptable = (float *)malloc(100*sizeof(float));
70 <        if (exptable == NULL) {
71 <                fprintf(stderr, "%s: out of memory in initmask\n", progname);
72 <                quit(1);
73 <        }
74 <        for (x = 0; x < 100; x++)
75 <                exptable[x] = exp(-x*0.1);
69 >        gtabsiz = 111*CHECKRAD*CHECKRAD/(rad*rad);
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);
75 >        for (x = 1; x < gtabsiz; x++)
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*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++)
89 >                gausstable[x] *= gaussN;
90 >                                        /* create ring averages table */
91 >        ringndx = (short *)malloc((2*orad*orad+1)*sizeof(short));
92 >        if (ringndx == NULL)
93 >                goto memerr;
94 >        for (x = 2*orad*orad+1; --x > orad*orad; )
95 >                ringndx[x] = -1;
96 >        do
97 >                ringndx[x] = sqrt((double)x);
98 >        while (x--);
99 >        ringsum = (float *)malloc((orad+1)*sizeof(float));
100 >        ringwt = (short *)malloc((orad+1)*sizeof(short));
101 >        warr = (float *)malloc(obarsize*obarsize*sizeof(float));
102 >        if (ringsum == NULL | ringwt == 0 | warr == NULL)
103 >                goto memerr;
104 >        return;
105 > memerr:
106 >        fprintf(stderr, "%s: out of memory in initmask\n", progname);
107 >        quit(1);
108   }
109  
110  
# Line 63 | Line 113 | COLOR  csum;
113   int  xcent, ycent;
114   int  c, r;
115   {
116 <        static int  wsum;
117 <        static double  d;
118 <        static int  y;
119 <        register int  x;
116 >        int  wsum;
117 >        double  d;
118 >        int  y;
119 >        register int  x, offs;
120          register COLOR  *scan;
121 <
121 >        
122          wsum = 0;
123          setcolor(csum, 0.0, 0.0, 0.0);
124 <        for (y = ycent+1-yrad; y <= ycent+yrad; y++) {
125 <                if (y < 0 || y >= yres)
126 <                        continue;
127 <                d = y_r < 1.0 ? y_r*y - r : y - ycent;
128 <                if (d > 0.5+FTINY || d < -0.5-FTINY)
129 <                        continue;
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+.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-xrad; x <= xcent+xrad; x++) {
132 <                        if (x < 0 || x >= xres)
131 >                for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
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 : x - xcent;
136 <                        if (d > 0.5+FTINY || d < -0.5-FTINY)
137 <                                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 98 | Line 151 | COLOR  csum;
151   int  xcent, ycent;
152   int  c, r;
153   {
154 <        static double  dy, dx, weight, wsum;
155 <        static COLOR  ctmp;
156 <        static int  y;
157 <        register int  x;
154 >        double  dy, dx, weight, wsum;
155 >        COLOR  ctmp;
156 >        int  y;
157 >        register int  x, offs;
158          register COLOR  *scan;
159  
160          wsum = FTINY;
161          setcolor(csum, 0.0, 0.0, 0.0);
162          for (y = ycent-yrad; y <= ycent+yrad; y++) {
163 <                if (y < 0 || y >= yres)
164 <                        continue;
163 >                if (y < 0) continue;
164 >                if (y >= yres) break;
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 || x >= xres)
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 = lookexp(-(dx*dx + dy*dy));
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 >
184 > dothresh(xcent, ycent, ccent, rcent)    /* gaussian threshold filter */
185 > int  xcent, ycent;
186 > int  ccent, rcent;
187 > {
188 >        double  d;
189 >        int  r, y, offs;
190 >        register int  c, x;
191 >        register float  *gscan;
192 >                                        /* compute ring sums */
193 >        bzero((char *)ringsum, (orad+1)*sizeof(float));
194 >        bzero((char *)ringwt, (orad+1)*sizeof(short));
195 >        for (r = -orad; r <= orad; r++) {
196 >                if (rcent+r < 0) continue;
197 >                if (rcent+r >= nrows) break;
198 >                gscan = greybar[(rcent+r)%obarsize];
199 >                for (c = -orad; c <= orad; c++) {
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+offs];
207 >                        ringwt[x]++;
208 >                }
209 >        }
210 >                                        /* filter each subpixel */
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+.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 >                        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 >                        sumans(x, y, rcent, ccent,
225 >                        pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
226 >                }
227 >        }
228 > }
229 >
230 >
231 > double
232 > pickfilt(p0)                    /* find filter multiplier for p0 */
233 > double  p0;
234 > {
235 >        double  m = 1.0;
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 = RSCA*CHECKRAD*rad*m + .5;
245 >                if (i > orad) i = orad;
246 >                avg = wsum = 0.0;
247 >                while (i--) {
248 >                        t = (i+.5)/(m*rad);
249 >                        t = lookgauss(t*t);
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 >                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(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 > }
289 >
290 >
291 > sumans(px, py, rcent, ccent, m)         /* sum input pixel to output */
292 > int  px, py;
293 > int  rcent, ccent;
294 > double  m;
295 > {
296 >        double  dy2, dx;
297 >        COLOR  pval, ctmp;
298 >        int  ksiz, r, offs;
299 >        double  pc, pr, norm;
300 >        register int  i, c;
301 >        register COLOR  *scan;
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;
311 >        if (ksiz > orad) ksiz = orad;
312 >                                                /* compute normalization */
313 >        norm = 0.0;
314 >        i = 0;
315 >        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
316 >                if (r < 0) continue;
317 >                if (r >= nrows) break;
318 >                dy2 = (pr - (r+.5))/(m*rad);
319 >                dy2 *= dy2;
320 >                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
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 + dy2);
327 >                }
328 >        }
329 >        norm = 1.0/norm;
330 >        if (x_c < 1.0) norm *= x_c;
331 >        if (y_r < 1.0) norm *= y_r;
332 >                                                /* sum pixels */
333 >        i = 0;
334 >        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
335 >                if (r < 0) continue;
336 >                if (r >= nrows) break;
337 >                scan = scoutbar[r%obarsize];
338 >                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
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+offs], ctmp);
346 >                }
347 >        }
348   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines