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.12 by greg, Thu Apr 4 13:01:35 1996 UTC vs.
Revision 2.20 by greg, Thu Dec 7 21:15:54 2023 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1993 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  "standard.h"
10 > #include  "pfilt.h"
11  
12 < #include  "color.h"
16 <
12 > #define  RSCA           1.13    /* square-radius multiplier: sqrt(4/PI) */
13   #define  TEPS           0.2     /* threshold proximity goal */
14   #define  REPS           0.1     /* radius proximity goal */
15  
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 –
29 – extern int  xres, yres;         /* resolution of input */
30 –
31 – extern double  x_c, y_r;        /* conversion factors */
32 –
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 –
16   float  *gausstable;             /* gauss lookup table */
17  
18   float  *ringsum;                /* sum of ring values */
# Line 54 | Line 20 | short  *ringwt;                        /* weight (count) of ring values */
20   short  *ringndx;                /* ring index table */
21   float  *warr;                   /* array of pixel weights */
22  
23 < extern double  (*ourbright)();  /* brightness computation function */
23 > #define  lookgauss(x)           gausstable[(int)(20.*(x)+.5)]
24  
25 < double  pickfilt();
25 > static double pickfilt(double  p0);
26 > static void sumans(int  px, int  py, int  rcent, int  ccent, double  m);
27  
61 – #define  lookgauss(x)           gausstable[(int)(10.*(x)+.5)]
28  
29 <
30 < initmask()                      /* initialize gaussian lookup table */
29 > void
30 > initmask(void)                  /* initialize gaussian lookup table */
31   {
32          int  gtabsiz;
33          double  gaussN;
34          double  d;
35 <        register int  x;
35 >        int  x;
36  
37 <        gtabsiz = 150*CHECKRAD;
37 >        gtabsiz = 444*CHECKRAD*CHECKRAD;
38          gausstable = (float *)malloc(gtabsiz*sizeof(float));
39          if (gausstable == NULL)
40                  goto memerr;
41          d = x_c*y_r*0.25/(rad*rad);
42 <        gausstable[0] = exp(-d)/sqrt(d);
42 >        gausstable[0] = exp(-d);
43          for (x = 1; x < gtabsiz; x++)
44 <                if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
44 >                if (x*0.05 <= d)
45                          gausstable[x] = gausstable[0];
46 +                else
47 +                        gausstable[x] = exp(-x*0.05);
48          if (obarsize == 0)
49                  return;
50                                          /* compute integral of filter */
51 <        gaussN = PI*sqrt(d)*exp(-d);            /* plateau */
52 <        for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
53 <                gaussN += 0.1*2.0*PI*exp(-d*d);
51 >        gaussN = PI*d*exp(-d);                  /* plateau */
52 >        for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1)
53 >                gaussN += 0.1*2.0*PI*d*exp(-d*d);
54                                          /* normalize filter */
55          gaussN = x_c*y_r/(rad*rad*gaussN);
56          for (x = 0; x < gtabsiz; x++)
# Line 99 | Line 67 | initmask()                     /* initialize gaussian lookup table */
67          ringsum = (float *)malloc((orad+1)*sizeof(float));
68          ringwt = (short *)malloc((orad+1)*sizeof(short));
69          warr = (float *)malloc(obarsize*obarsize*sizeof(float));
70 <        if (ringsum == NULL | ringwt == 0 | warr == NULL)
70 >        if ((ringsum == NULL) | (ringwt == 0) | (warr == NULL))
71                  goto memerr;
72          return;
73   memerr:
# Line 108 | Line 76 | memerr:
76   }
77  
78  
79 < dobox(csum, xcent, ycent, c, r)                 /* simple box filter */
80 < COLOR  csum;
81 < int  xcent, ycent;
82 < int  c, r;
79 > void
80 > dobox(                  /* simple box filter */
81 >        COLOR  csum,
82 >        int  xcent,
83 >        int  ycent,
84 >        int  c,
85 >        int  r
86 > )
87   {
88          int  wsum;
89          double  d;
90          int  y;
91 <        register int  x, offs;
92 <        register COLOR  *scan;
91 >        int  x, offs;
92 >        COLOR   *scan;
93          
94          wsum = 0;
95          setcolor(csum, 0.0, 0.0, 0.0);
96          for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
97                  if (y < 0) continue;
98                  if (y >= yres) break;
99 <                d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
99 >                d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent);
100                  if (d < -0.5) continue;
101                  if (d >= 0.5) break;
102                  scan = scanin[y%barsize];
# Line 132 | Line 104 | int  c, r;
104                          offs = x < 0 ? xres : x >= xres ? -xres : 0;
105                          if (offs && !wrapfilt)
106                                  continue;
107 <                        d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
107 >                        d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent);
108                          if (d < -0.5) continue;
109                          if (d >= 0.5) break;
110                          wsum++;
111                          addcolor(csum, scan[x+offs]);
112                  }
113          }
114 <        if (wsum > 1)
115 <                scalecolor(csum, 1.0/wsum);
114 >        if (wsum > 1) {
115 >                d = 1.0/wsum;
116 >                scalecolor(csum, d);
117 >        }
118   }
119  
120  
121 < dogauss(csum, xcent, ycent, c, r)               /* gaussian filter */
122 < COLOR  csum;
123 < int  xcent, ycent;
124 < int  c, r;
121 > void
122 > dogauss(                /* gaussian filter */
123 >        COLOR  csum,
124 >        int  xcent,
125 >        int  ycent,
126 >        int  c,
127 >        int  r
128 > )
129   {
130          double  dy, dx, weight, wsum;
131          COLOR  ctmp;
132          int  y;
133 <        register int  x, offs;
134 <        register COLOR  *scan;
133 >        int  x, offs;
134 >        COLOR   *scan;
135  
136          wsum = FTINY;
137          setcolor(csum, 0.0, 0.0, 0.0);
# Line 174 | Line 152 | int  c, r;
152                          addcolor(csum, ctmp);
153                  }
154          }
155 <        scalecolor(csum, 1.0/wsum);
155 >        weight = 1.0/wsum;
156 >        scalecolor(csum, weight);
157   }
158  
159  
160 < dothresh(xcent, ycent, ccent, rcent)    /* gaussian threshold filter */
161 < int  xcent, ycent;
162 < int  ccent, rcent;
160 > void
161 > dothresh(       /* gaussian threshold filter */
162 >        int  xcent,
163 >        int  ycent,
164 >        int  ccent,
165 >        int  rcent
166 > )
167   {
168          double  d;
169          int  r, y, offs;
170 <        register int  c, x;
171 <        register float  *gscan;
170 >        int  c, x;
171 >        float  *gscan;
172                                          /* compute ring sums */
173 <        bzero((char *)ringsum, (orad+1)*sizeof(float));
174 <        bzero((char *)ringwt, (orad+1)*sizeof(short));
173 >        memset((char *)ringsum, '\0', (orad+1)*sizeof(float));
174 >        memset((char *)ringwt, '\0', (orad+1)*sizeof(short));
175          for (r = -orad; r <= orad; r++) {
176                  if (rcent+r < 0) continue;
177                  if (rcent+r >= nrows) break;
# Line 208 | Line 191 | int  ccent, rcent;
191          for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
192                  if (y < 0) continue;
193                  if (y >= yres) break;
194 <                d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent);
194 >                d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent);
195                  if (d < -0.5) continue;
196                  if (d >= 0.5) break;
197                  for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
198                          offs = x < 0 ? xres : x >= xres ? -xres : 0;
199                          if (offs && !wrapfilt)
200                                  continue;
201 <                        d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
201 >                        d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
202                          if (d < -0.5) continue;
203                          if (d >= 0.5) break;
204                          sumans(x, y, rcent, ccent,
# Line 225 | Line 208 | int  ccent, rcent;
208   }
209  
210  
211 < double
212 < pickfilt(p0)                    /* find filter multiplier for p0 */
213 < double  p0;
211 > static double
212 > pickfilt(                       /* find filter multiplier for p0 */
213 >        double  p0
214 > )
215   {
216          double  m = 1.0;
217          double  t, num, denom, avg, wsum;
218          double  mlimit[2];
219 <        int  ilimit = 4/TEPS;
220 <        register int  i;
219 >        int  ilimit = 4.0/TEPS;
220 >        int  i;
221                                  /* iterative search for m */
222          mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
223          do {
224                                          /* compute grey weighted average */
225 <                i = 1.25*CHECKRAD*rad*m + .5;
225 >                i = RSCA*CHECKRAD*rad*m + .5;
226                  if (i > orad) i = orad;
227                  avg = wsum = 0.0;
228                  while (i--) {
# Line 285 | Line 269 | double  p0;
269   }
270  
271  
272 < sumans(px, py, rcent, ccent, m)         /* sum input pixel to output */
273 < int  px, py;
274 < int  rcent, ccent;
275 < double  m;
272 > static void
273 > sumans(         /* sum input pixel to output */
274 >        int  px,
275 >        int  py,
276 >        int  rcent,
277 >        int  ccent,
278 >        double  m
279 > )
280   {
281 <        double  dy, dx;
281 >        double  dy2, dx;
282          COLOR  pval, ctmp;
283          int  ksiz, r, offs;
284          double  pc, pr, norm;
285 <        register int  i, c;
286 <        register COLOR  *scan;
285 >        int  i, c;
286 >        COLOR   *scan;
287          /*
288           * This normalization method fails at the picture borders because
289           * a different number of input pixels contribute there.
# Line 312 | Line 300 | double  m;
300          for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
301                  if (r < 0) continue;
302                  if (r >= nrows) break;
303 <                dy = (pr - (r+.5))/(m*rad);
303 >                dy2 = (pr - (r+.5))/(m*rad);
304 >                dy2 *= dy2;
305                  for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
306                          if (!wrapfilt) {
307                                  if (c < 0) continue;
308                                  if (c >= ncols) break;
309                          }
310                          dx = (pc - (c+.5))/(m*rad);
311 <                        norm += warr[i++] = lookgauss(dx*dx + dy*dy);
311 >                        norm += warr[i++] = lookgauss(dx*dx + dy2);
312                  }
313          }
314          norm = 1.0/norm;
# Line 333 | Line 322 | double  m;
322                  scan = scoutbar[r%obarsize];
323                  for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
324                          offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
325 <                        if (offs && !wrapfilt)
325 >                        if ((offs != 0) & !wrapfilt)
326                                  continue;
327                          copycolor(ctmp, pval);
328                          dx = norm*warr[i++];

Diff Legend

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