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 1.1 by greg, Thu Feb 2 10:49:26 1989 UTC vs.
Revision 2.13 by greg, Thu Dec 12 14:24:00 1996 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 1986 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 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  "color.h"
16  
17 < #define  FTINY          1e-7
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 +
22   extern double  rad;             /* output pixel radius for filtering */
23  
24 + extern double  thresh;          /* maximum contribution for subpixel */
25 +
26   extern int  nrows;              /* number of rows for output */
27   extern int  ncols;              /* number of columns for output */
28  
24 extern int  boxfilt;            /* do box filtering? */
25
29   extern int  xres, yres;         /* resolution of input */
30  
31   extern double  x_c, y_r;        /* conversion factors */
32  
33 < extern int  xrad;               /* x window size */
34 < extern int  yrad;               /* y window size */
33 > extern int  xrad;               /* x search radius */
34 > extern int  yrad;               /* y search radius */
35 > extern int  xbrad;              /* x box size */
36 > extern int  ybrad;              /* y box size */
37  
38   extern int  barsize;            /* size of input scan bar */
39   extern COLOR  **scanin;         /* input scan bar */
40   extern COLOR  *scanout;         /* output scan line */
41 + extern COLOR  **scoutbar;       /* output scan bar (if thresh > 0) */
42 + 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  *exptable;               /* exponent table */
50 > float  *gausstable;             /* gauss lookup table */
51  
52 < #define  lookexp(x)             exptable[(int)(-10.*(x)+.5)]
52 > float  *ringsum;                /* sum of ring values */
53 > 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)]
62 +
63 +
64   initmask()                      /* initialize gaussian lookup table */
65   {
66 <        extern char  *malloc();
67 <        extern double  exp();
66 >        int  gtabsiz;
67 >        double  gaussN;
68 >        double  d;
69          register int  x;
70  
71 <        exptable = (float *)malloc(100*sizeof(float));
72 <        if (exptable == NULL) {
73 <                fprintf(stderr, "%s: out of memory in initmask\n", progname);
74 <                quit(1);
75 <        }
76 <        for (x = 0; x < 100; x++)
77 <                exptable[x] = exp(-x*0.1);
71 >        gtabsiz = 150*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);
77 >        for (x = 1; x < gtabsiz; x++)
78 >                if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
79 >                        gausstable[x] = gausstable[0];
80 >        if (obarsize == 0)
81 >                return;
82 >                                        /* compute integral of filter */
83 >        gaussN = PI*sqrt(d)*exp(-d);            /* plateau */
84 >        for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
85 >                gaussN += 0.1*2.0*PI*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 62 | 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;
120 <        register COLOR  *scan;
121 <
116 >        int  wsum;
117 >        double  d;
118 >        int  y;
119 >        register int  x, offs;
120 >        register COLOR  *scan;
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)
# Line 97 | Line 149 | COLOR  csum;
149   int  xcent, ycent;
150   int  c, r;
151   {
152 <        static double  dy, dx, weight, wsum;
153 <        static COLOR  ctmp;
154 <        static int  y;
155 <        register int  x;
156 <        register COLOR  *scan;
152 >        double  dy, dx, weight, wsum;
153 >        COLOR  ctmp;
154 >        int  y;
155 >        register int  x, offs;
156 >        register COLOR  *scan;
157  
158          wsum = FTINY;
159          setcolor(csum, 0.0, 0.0, 0.0);
160 <        for (y = ycent+1-yrad; y <= ycent+yrad; y++) {
161 <                if (y < 0 || y >= yres)
162 <                        continue;
163 <                dy = (y_r*y - r)/rad;
160 >        for (y = ycent-yrad; y <= ycent+yrad; y++) {
161 >                if (y < 0) continue;
162 >                if (y >= yres) break;
163 >                dy = (y_r*(y+.5) - (r+.5))/rad;
164                  scan = scanin[y%barsize];
165 <                for (x = xcent+1-xrad; x <= xcent+xrad; x++) {
166 <                        if (x < 0 || x >= xres)
165 >                for (x = xcent-xrad; x <= xcent+xrad; x++) {
166 >                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
167 >                        if (offs && !wrapfilt)
168                                  continue;
169 <                        dx = (x_c*x - c)/rad;
170 <                        weight = lookexp(-(dx*dx + dy*dy));
169 >                        dx = (x_c*(x+.5) - (c+.5))/rad;
170 >                        weight = lookgauss(dx*dx + dy*dy);
171                          wsum += weight;
172 <                        copycolor(ctmp, scan[x]);
172 >                        copycolor(ctmp, scan[x+offs]);
173                          scalecolor(ctmp, weight);
174                          addcolor(csum, ctmp);
175                  }
176          }
177          scalecolor(csum, 1.0/wsum);
178 + }
179 +
180 +
181 + dothresh(xcent, ycent, ccent, rcent)    /* gaussian threshold filter */
182 + int  xcent, ycent;
183 + int  ccent, rcent;
184 + {
185 +        double  d;
186 +        int  r, y, offs;
187 +        register int  c, x;
188 +        register float  *gscan;
189 +                                        /* compute ring sums */
190 +        bzero((char *)ringsum, (orad+1)*sizeof(float));
191 +        bzero((char *)ringwt, (orad+1)*sizeof(short));
192 +        for (r = -orad; r <= orad; r++) {
193 +                if (rcent+r < 0) continue;
194 +                if (rcent+r >= nrows) break;
195 +                gscan = greybar[(rcent+r)%obarsize];
196 +                for (c = -orad; c <= orad; c++) {
197 +                        offs = ccent+c < 0 ? ncols :
198 +                                        ccent+c >= ncols ? -ncols : 0;
199 +                        if (offs && !wrapfilt)
200 +                                continue;
201 +                        x = ringndx[c*c + r*r];
202 +                        if (x < 0) continue;
203 +                        ringsum[x] += gscan[ccent+c+offs];
204 +                        ringwt[x]++;
205 +                }
206 +        }
207 +                                        /* filter each subpixel */
208 +        for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
209 +                if (y < 0) continue;
210 +                if (y >= yres) break;
211 +                d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent);
212 +                if (d < -0.5) continue;
213 +                if (d >= 0.5) break;
214 +                for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
215 +                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
216 +                        if (offs && !wrapfilt)
217 +                                continue;
218 +                        d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
219 +                        if (d < -0.5) continue;
220 +                        if (d >= 0.5) break;
221 +                        sumans(x, y, rcent, ccent,
222 +                        pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
223 +                }
224 +        }
225 + }
226 +
227 +
228 + double
229 + pickfilt(p0)                    /* find filter multiplier for p0 */
230 + double  p0;
231 + {
232 +        double  m = 1.0;
233 +        double  t, num, denom, avg, wsum;
234 +        double  mlimit[2];
235 +        int  ilimit = 4/TEPS;
236 +        register int  i;
237 +                                /* iterative search for m */
238 +        mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
239 +        do {
240 +                                        /* compute grey weighted average */
241 +                i = 1.25*CHECKRAD*rad*m + .5;
242 +                if (i > orad) i = orad;
243 +                avg = wsum = 0.0;
244 +                while (i--) {
245 +                        t = (i+.5)/(m*rad);
246 +                        t = lookgauss(t*t);
247 +                        avg += t*ringsum[i];
248 +                        wsum += t*ringwt[i];
249 +                }
250 +                if (avg < 1e-20)                /* zero inclusive average */
251 +                        return(1.0);
252 +                avg /= wsum;
253 +                                        /* check threshold */
254 +                denom = m*m/gausstable[0] - p0/avg;
255 +                if (denom <= FTINY) {           /* zero exclusive average */
256 +                        if (m >= mlimit[1]-REPS)
257 +                                break;
258 +                        m = mlimit[1];
259 +                        continue;
260 +                }
261 +                num = p0/avg - 1.0;
262 +                if (num < 0.0) num = -num;
263 +                t = num/denom;
264 +                if (t <= thresh) {
265 +                        if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS)
266 +                                break;
267 +                } else {
268 +                        if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS)
269 +                                break;
270 +                }
271 +                t = m;                  /* remember current m */
272 +                                        /* next guesstimate */
273 +                m = sqrt(gausstable[0]*(num/thresh + p0/avg));
274 +                if (m < t) {            /* bound it */
275 +                        if (m <= mlimit[0]+FTINY)
276 +                                m = 0.5*(mlimit[0] + t);
277 +                        mlimit[1] = t;
278 +                } else {
279 +                        if (m >= mlimit[1]-FTINY)
280 +                                m = 0.5*(mlimit[1] + t);
281 +                        mlimit[0] = t;
282 +                }
283 +        } while (--ilimit > 0);
284 +        return(m);
285 + }
286 +
287 +
288 + sumans(px, py, rcent, ccent, m)         /* sum input pixel to output */
289 + int  px, py;
290 + int  rcent, ccent;
291 + double  m;
292 + {
293 +        double  dy, dx;
294 +        COLOR  pval, ctmp;
295 +        int  ksiz, r, offs;
296 +        double  pc, pr, norm;
297 +        register int  i, c;
298 +        register COLOR  *scan;
299 +        /*
300 +         * This normalization method fails at the picture borders because
301 +         * a different number of input pixels contribute there.
302 +         */
303 +        scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0);
304 +        copycolor(pval, scan[px]);
305 +        pc = x_c*(px+.5);
306 +        pr = y_r*(py+.5);
307 +        ksiz = CHECKRAD*m*rad + 1;
308 +        if (ksiz > orad) ksiz = orad;
309 +                                                /* compute normalization */
310 +        norm = 0.0;
311 +        i = 0;
312 +        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
313 +                if (r < 0) continue;
314 +                if (r >= nrows) break;
315 +                dy = (pr - (r+.5))/(m*rad);
316 +                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
317 +                        if (!wrapfilt) {
318 +                                if (c < 0) continue;
319 +                                if (c >= ncols) break;
320 +                        }
321 +                        dx = (pc - (c+.5))/(m*rad);
322 +                        norm += warr[i++] = lookgauss(dx*dx + dy*dy);
323 +                }
324 +        }
325 +        norm = 1.0/norm;
326 +        if (x_c < 1.0) norm *= x_c;
327 +        if (y_r < 1.0) norm *= y_r;
328 +                                                /* sum pixels */
329 +        i = 0;
330 +        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
331 +                if (r < 0) continue;
332 +                if (r >= nrows) break;
333 +                scan = scoutbar[r%obarsize];
334 +                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
335 +                        offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
336 +                        if (offs && !wrapfilt)
337 +                                continue;
338 +                        copycolor(ctmp, pval);
339 +                        dx = norm*warr[i++];
340 +                        scalecolor(ctmp, dx);
341 +                        addcolor(scan[c+offs], ctmp);
342 +                }
343 +        }
344   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines