ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pf3.c
Revision: 2.18
Committed: Sun Mar 28 20:33:14 2004 UTC (20 years ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad5R1, rad3R7P2, rad3R7P1, rad4R2, rad4R1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9, rad4R2P1
Changes since 2.17: +42 -54 lines
Log Message:
Continued ANSIfication, and other fixes and clarifications.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: pf3.c,v 2.17 2003/07/27 22:12:03 schorsch Exp $";
3 #endif
4 /*
5 * pf3.c - routines for gaussian and box filtering
6 *
7 * 8/13/86
8 */
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
21 float *gausstable; /* gauss lookup table */
22
23 float *ringsum; /* sum of ring values */
24 short *ringwt; /* weight (count) of ring values */
25 short *ringndx; /* ring index table */
26 float *warr; /* array of pixel weights */
27
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
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 = 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);
48 for (x = 1; x < gtabsiz; x++)
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*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++)
62 gausstable[x] *= gaussN;
63 /* create ring averages table */
64 ringndx = (short *)malloc((2*orad*orad+1)*sizeof(short));
65 if (ringndx == NULL)
66 goto memerr;
67 for (x = 2*orad*orad+1; --x > orad*orad; )
68 ringndx[x] = -1;
69 do
70 ringndx[x] = sqrt((double)x);
71 while (x--);
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))
76 goto memerr;
77 return;
78 memerr:
79 fprintf(stderr, "%s: out of memory in initmask\n", progname);
80 quit(1);
81 }
82
83
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, offs;
97 register COLOR *scan;
98
99 wsum = 0;
100 setcolor(csum, 0.0, 0.0, 0.0);
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+.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 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+offs]);
117 }
118 }
119 if (wsum > 1) {
120 d = 1.0/wsum;
121 scalecolor(csum, d);
122 }
123 }
124
125
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, offs;
139 register COLOR *scan;
140
141 wsum = FTINY;
142 setcolor(csum, 0.0, 0.0, 0.0);
143 for (y = ycent-yrad; y <= ycent+yrad; y++) {
144 if (y < 0) continue;
145 if (y >= yres) break;
146 dy = (y_r*(y+.5) - (r+.5))/rad;
147 scan = scanin[y%barsize];
148 for (x = xcent-xrad; x <= xcent+xrad; x++) {
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+offs]);
156 scalecolor(ctmp, weight);
157 addcolor(csum, ctmp);
158 }
159 }
160 weight = 1.0/wsum;
161 scalecolor(csum, weight);
162 }
163
164
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, offs;
175 register int c, x;
176 register float *gscan;
177 /* compute ring sums */
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 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+offs];
192 ringwt[x]++;
193 }
194 }
195 /* filter each subpixel */
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+.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 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 sumans(x, y, rcent, ccent,
210 pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
211 }
212 }
213 }
214
215
216 static double
217 pickfilt( /* find filter multiplier for p0 */
218 double p0
219 )
220 {
221 double m = 1.0;
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 = RSCA*CHECKRAD*rad*m + .5;
231 if (i > orad) i = orad;
232 avg = wsum = 0.0;
233 while (i--) {
234 t = (i+.5)/(m*rad);
235 t = lookgauss(t*t);
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 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(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 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 dy2, dx;
287 COLOR pval, ctmp;
288 int ksiz, r, offs;
289 double pc, pr, norm;
290 register int i, c;
291 register COLOR *scan;
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;
301 if (ksiz > orad) ksiz = orad;
302 /* compute normalization */
303 norm = 0.0;
304 i = 0;
305 for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
306 if (r < 0) continue;
307 if (r >= nrows) break;
308 dy2 = (pr - (r+.5))/(m*rad);
309 dy2 *= dy2;
310 for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
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 + dy2);
317 }
318 }
319 norm = 1.0/norm;
320 if (x_c < 1.0) norm *= x_c;
321 if (y_r < 1.0) norm *= y_r;
322 /* sum pixels */
323 i = 0;
324 for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
325 if (r < 0) continue;
326 if (r >= nrows) break;
327 scan = scoutbar[r%obarsize];
328 for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
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+offs], ctmp);
336 }
337 }
338 }