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.9 by greg, Fri Jul 9 10:27:52 1993 UTC vs.
Revision 2.18 by schorsch, Sun Mar 28 20:33:14 2004 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   *
# 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 + #include  "pfilt.h"
16  
17 + #define  RSCA           1.13    /* square-radius multiplier: sqrt(4/PI) */
18   #define  TEPS           0.2     /* threshold proximity goal */
19 + #define  REPS           0.1     /* radius proximity goal */
20  
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 –
28 – extern int  xres, yres;         /* resolution of input */
29 –
30 – extern double  x_c, y_r;        /* conversion factors */
31 –
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 –
21   float  *gausstable;             /* gauss lookup table */
22  
23   float  *ringsum;                /* sum of ring values */
# Line 51 | Line 25 | short  *ringwt;                        /* weight (count) of ring values */
25   short  *ringndx;                /* ring index table */
26   float  *warr;                   /* array of pixel weights */
27  
54 – double  pickfilt();
55 –
28   #define  lookgauss(x)           gausstable[(int)(10.*(x)+.5)]
29  
30 + static double pickfilt(double  p0);
31 + static void sumans(int  px, int  py, int  rcent, int  ccent, double  m);
32  
33 < initmask()                      /* initialize gaussian lookup table */
33 >
34 > extern void
35 > initmask(void)                  /* initialize gaussian lookup table */
36   {
37          int  gtabsiz;
38          double  gaussN;
39          double  d;
40          register int  x;
41  
42 <        gtabsiz = 150*CHECKRAD;
42 >        gtabsiz = 111*CHECKRAD*CHECKRAD;
43          gausstable = (float *)malloc(gtabsiz*sizeof(float));
44          if (gausstable == NULL)
45                  goto memerr;
46          d = x_c*y_r*0.25/(rad*rad);
47 <        gausstable[0] = exp(-d)/sqrt(d);
47 >        gausstable[0] = exp(-d);
48          for (x = 1; x < gtabsiz; x++)
49 <                if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
49 >                if (x*0.1 <= d)
50                          gausstable[x] = gausstable[0];
51 +                else
52 +                        gausstable[x] = exp(-x*0.1);
53          if (obarsize == 0)
54                  return;
55                                          /* compute integral of filter */
56 <        gaussN = PI*sqrt(d)*exp(-d);            /* plateau */
57 <        for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
58 <                gaussN += 0.1*2.0*PI*exp(-d*d);
56 >        gaussN = PI*d*exp(-d);                  /* plateau */
57 >        for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1)
58 >                gaussN += 0.1*2.0*PI*d*exp(-d*d);
59                                          /* normalize filter */
60          gaussN = x_c*y_r/(rad*rad*gaussN);
61          for (x = 0; x < gtabsiz; x++)
# Line 94 | Line 72 | initmask()                     /* initialize gaussian lookup table */
72          ringsum = (float *)malloc((orad+1)*sizeof(float));
73          ringwt = (short *)malloc((orad+1)*sizeof(short));
74          warr = (float *)malloc(obarsize*obarsize*sizeof(float));
75 <        if (ringsum == NULL | ringwt == 0 | warr == NULL)
75 >        if ((ringsum == NULL) | (ringwt == 0) | (warr == NULL))
76                  goto memerr;
77          return;
78   memerr:
# Line 103 | Line 81 | memerr:
81   }
82  
83  
84 < dobox(csum, xcent, ycent, c, r)                 /* simple box filter */
85 < COLOR  csum;
86 < int  xcent, ycent;
87 < int  c, r;
84 > extern void
85 > dobox(                  /* simple box filter */
86 >        COLOR  csum,
87 >        int  xcent,
88 >        int  ycent,
89 >        int  c,
90 >        int  r
91 > )
92   {
93          int  wsum;
94          double  d;
95          int  y;
96 <        register int  x;
96 >        register int  x, offs;
97          register COLOR  *scan;
98          
99          wsum = 0;
# Line 119 | Line 101 | int  c, r;
101          for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
102                  if (y < 0) continue;
103                  if (y >= yres) break;
104 <                d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
104 >                d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent);
105                  if (d < -0.5) continue;
106                  if (d >= 0.5) break;
107                  scan = scanin[y%barsize];
108                  for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
109 <                        if (x < 0) continue;
110 <                        if (x >= xres) break;
111 <                        d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
109 >                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
110 >                        if (offs && !wrapfilt)
111 >                                continue;
112 >                        d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent);
113                          if (d < -0.5) continue;
114                          if (d >= 0.5) break;
115                          wsum++;
116 <                        addcolor(csum, scan[x]);
116 >                        addcolor(csum, scan[x+offs]);
117                  }
118          }
119 <        if (wsum > 1)
120 <                scalecolor(csum, 1.0/wsum);
119 >        if (wsum > 1) {
120 >                d = 1.0/wsum;
121 >                scalecolor(csum, d);
122 >        }
123   }
124  
125  
126 < dogauss(csum, xcent, ycent, c, r)               /* gaussian filter */
127 < COLOR  csum;
128 < int  xcent, ycent;
129 < int  c, r;
126 > extern void
127 > dogauss(                /* gaussian filter */
128 >        COLOR  csum,
129 >        int  xcent,
130 >        int  ycent,
131 >        int  c,
132 >        int  r
133 > )
134   {
135          double  dy, dx, weight, wsum;
136          COLOR  ctmp;
137          int  y;
138 <        register int  x;
138 >        register int  x, offs;
139          register COLOR  *scan;
140  
141          wsum = FTINY;
# Line 157 | Line 146 | int  c, r;
146                  dy = (y_r*(y+.5) - (r+.5))/rad;
147                  scan = scanin[y%barsize];
148                  for (x = xcent-xrad; x <= xcent+xrad; x++) {
149 <                        if (x < 0) continue;
150 <                        if (x >= xres) break;
149 >                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
150 >                        if (offs && !wrapfilt)
151 >                                continue;
152                          dx = (x_c*(x+.5) - (c+.5))/rad;
153                          weight = lookgauss(dx*dx + dy*dy);
154                          wsum += weight;
155 <                        copycolor(ctmp, scan[x]);
155 >                        copycolor(ctmp, scan[x+offs]);
156                          scalecolor(ctmp, weight);
157                          addcolor(csum, ctmp);
158                  }
159          }
160 <        scalecolor(csum, 1.0/wsum);
160 >        weight = 1.0/wsum;
161 >        scalecolor(csum, weight);
162   }
163  
164  
165 < dothresh(xcent, ycent, ccent, rcent)    /* gaussian threshold filter */
166 < int  xcent, ycent;
167 < int  ccent, rcent;
165 > extern void
166 > dothresh(       /* gaussian threshold filter */
167 >        int  xcent,
168 >        int  ycent,
169 >        int  ccent,
170 >        int  rcent
171 > )
172   {
173          double  d;
174 <        int  r, y;
174 >        int  r, y, offs;
175          register int  c, x;
176          register float  *gscan;
182 – #define pval gscan
177                                          /* compute ring sums */
178 <        bzero((char *)ringsum, (orad+1)*sizeof(float));
179 <        bzero((char *)ringwt, (orad+1)*sizeof(short));
178 >        memset((char *)ringsum, '\0', (orad+1)*sizeof(float));
179 >        memset((char *)ringwt, '\0', (orad+1)*sizeof(short));
180          for (r = -orad; r <= orad; r++) {
181                  if (rcent+r < 0) continue;
182                  if (rcent+r >= nrows) break;
183                  gscan = greybar[(rcent+r)%obarsize];
184                  for (c = -orad; c <= orad; c++) {
185 <                        if (ccent+c < 0) continue;
186 <                        if (ccent+c >= ncols) break;
185 >                        offs = ccent+c < 0 ? ncols :
186 >                                        ccent+c >= ncols ? -ncols : 0;
187 >                        if (offs && !wrapfilt)
188 >                                continue;
189                          x = ringndx[c*c + r*r];
190                          if (x < 0) continue;
191 <                        ringsum[x] += gscan[ccent+c];
191 >                        ringsum[x] += gscan[ccent+c+offs];
192                          ringwt[x]++;
193                  }
194          }
# Line 200 | Line 196 | int  ccent, rcent;
196          for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
197                  if (y < 0) continue;
198                  if (y >= yres) break;
199 <                d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent);
199 >                d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent);
200                  if (d < -0.5) continue;
201                  if (d >= 0.5) break;
202                  for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
203 <                        if (x < 0) continue;
204 <                        if (x >= xres) break;
205 <                        d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
203 >                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
204 >                        if (offs && !wrapfilt)
205 >                                continue;
206 >                        d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
207                          if (d < -0.5) continue;
208                          if (d >= 0.5) break;
209 <                        pval = scanin[y%barsize][x];
210 <                        sumans(x, y, rcent, ccent, pickfilt(bright(pval)));
209 >                        sumans(x, y, rcent, ccent,
210 >                        pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
211                  }
212          }
216 – #undef pval
213   }
214  
215  
216 < double
217 < pickfilt(p0)                    /* find filter multiplier for p0 */
218 < double  p0;
216 > static double
217 > pickfilt(                       /* find filter multiplier for p0 */
218 >        double  p0
219 > )
220   {
221          double  m = 1.0;
222 <        double  t, avg, wsum;
223 <        int  ilimit = 4/TEPS;
222 >        double  t, num, denom, avg, wsum;
223 >        double  mlimit[2];
224 >        int  ilimit = 4.0/TEPS;
225          register int  i;
226                                  /* iterative search for m */
227 +        mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
228          do {
229                                          /* compute grey weighted average */
230 <                i = 1.25*CHECKRAD*rad*m + .5;
230 >                i = RSCA*CHECKRAD*rad*m + .5;
231                  if (i > orad) i = orad;
232                  avg = wsum = 0.0;
233                  while (i--) {
# Line 237 | Line 236 | double  p0;
236                          avg += t*ringsum[i];
237                          wsum += t*ringwt[i];
238                  }
239 +                if (avg < 1e-20)                /* zero inclusive average */
240 +                        return(1.0);
241                  avg /= wsum;
242                                          /* check threshold */
243 <                t = p0 - avg;
244 <                if (t < 0.0) t = -t;
245 <                t *= gausstable[0]/(m*m*avg);
246 <                if (t <= thresh && (m <= 1.0+FTINY ||
247 <                                (thresh-t)/thresh <= TEPS))
248 <                        break;
249 <                if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY ||
250 <                                (t-thresh)/thresh <= TEPS))
251 <                        break;
243 >                denom = m*m/gausstable[0] - p0/avg;
244 >                if (denom <= FTINY) {           /* zero exclusive average */
245 >                        if (m >= mlimit[1]-REPS)
246 >                                break;
247 >                        m = mlimit[1];
248 >                        continue;
249 >                }
250 >                num = p0/avg - 1.0;
251 >                if (num < 0.0) num = -num;
252 >                t = num/denom;
253 >                if (t <= thresh) {
254 >                        if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS)
255 >                                break;
256 >                } else {
257 >                        if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS)
258 >                                break;
259 >                }
260 >                t = m;                  /* remember current m */
261                                          /* next guesstimate */
262 <                m *= sqrt(t/thresh);
263 <                if (m < 1.0) m = 1.0;
264 <                else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD;
262 >                m = sqrt(gausstable[0]*(num/thresh + p0/avg));
263 >                if (m < t) {            /* bound it */
264 >                        if (m <= mlimit[0]+FTINY)
265 >                                m = 0.5*(mlimit[0] + t);
266 >                        mlimit[1] = t;
267 >                } else {
268 >                        if (m >= mlimit[1]-FTINY)
269 >                                m = 0.5*(mlimit[1] + t);
270 >                        mlimit[0] = t;
271 >                }
272          } while (--ilimit > 0);
273          return(m);
274   }
275  
276  
277 < sumans(px, py, rcent, ccent, m)         /* sum input pixel to output */
278 < int  px, py;
279 < int  rcent, ccent;
280 < double  m;
277 > static void
278 > sumans(         /* sum input pixel to output */
279 >        int  px,
280 >        int  py,
281 >        int  rcent,
282 >        int  ccent,
283 >        double  m
284 > )
285   {
286 <        double  dy, dx;
286 >        double  dy2, dx;
287          COLOR  pval, ctmp;
288 <        int  ksiz, r;
288 >        int  ksiz, r, offs;
289          double  pc, pr, norm;
290          register int  i, c;
291          register COLOR  *scan;
292 <
293 <        copycolor(pval, scanin[py%barsize][px]);
292 >        /*
293 >         * This normalization method fails at the picture borders because
294 >         * a different number of input pixels contribute there.
295 >         */
296 >        scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0);
297 >        copycolor(pval, scan[px]);
298          pc = x_c*(px+.5);
299          pr = y_r*(py+.5);
300          ksiz = CHECKRAD*m*rad + 1;
# Line 280 | Line 305 | double  m;
305          for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
306                  if (r < 0) continue;
307                  if (r >= nrows) break;
308 <                dy = (pr - (r+.5))/(m*rad);
308 >                dy2 = (pr - (r+.5))/(m*rad);
309 >                dy2 *= dy2;
310                  for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
311 <                        if (c < 0) continue;
312 <                        if (c >= ncols) break;
311 >                        if (!wrapfilt) {
312 >                                if (c < 0) continue;
313 >                                if (c >= ncols) break;
314 >                        }
315                          dx = (pc - (c+.5))/(m*rad);
316 <                        norm += warr[i++] = lookgauss(dx*dx + dy*dy);
316 >                        norm += warr[i++] = lookgauss(dx*dx + dy2);
317                  }
318          }
319          norm = 1.0/norm;
# Line 298 | Line 326 | double  m;
326                  if (r >= nrows) break;
327                  scan = scoutbar[r%obarsize];
328                  for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
329 <                        if (c < 0) continue;
330 <                        if (c >= ncols) break;
329 >                        offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
330 >                        if (offs && !wrapfilt)
331 >                                continue;
332                          copycolor(ctmp, pval);
333                          dx = norm*warr[i++];
334                          scalecolor(ctmp, dx);
335 <                        addcolor(scan[c], ctmp);
335 >                        addcolor(scan[c+offs], ctmp);
336                  }
337          }
338   }

Diff Legend

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