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 |
|
* |
11 |
|
|
12 |
|
#include "color.h" |
13 |
|
|
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 |
|
|
66 |
|
double d; |
67 |
|
register int x; |
68 |
|
|
69 |
< |
gtabsiz = 150*CHECKRAD; |
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)/sqrt(d); |
74 |
> |
gausstable[0] = exp(-d); |
75 |
|
for (x = 1; x < gtabsiz; x++) |
76 |
< |
if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0]) |
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*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); |
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++) |
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 |
|
|
176 |
|
addcolor(csum, ctmp); |
177 |
|
} |
178 |
|
} |
179 |
< |
scalecolor(csum, 1.0/wsum); |
179 |
> |
weight = 1.0/wsum; |
180 |
> |
scalecolor(csum, weight); |
181 |
|
} |
182 |
|
|
183 |
|
|
235 |
|
double m = 1.0; |
236 |
|
double t, num, denom, avg, wsum; |
237 |
|
double mlimit[2]; |
238 |
< |
int ilimit = 4/TEPS; |
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 = 1.25*CHECKRAD*rad*m + .5; |
244 |
> |
i = RSCA*CHECKRAD*rad*m + .5; |
245 |
|
if (i > orad) i = orad; |
246 |
|
avg = wsum = 0.0; |
247 |
|
while (i--) { |
293 |
|
int rcent, ccent; |
294 |
|
double m; |
295 |
|
{ |
296 |
< |
double dy, dx; |
296 |
> |
double dy2, dx; |
297 |
|
COLOR pval, ctmp; |
298 |
|
int ksiz, r, offs; |
299 |
|
double pc, pr, norm; |
315 |
|
for (r = rcent-ksiz; r <= rcent+ksiz; r++) { |
316 |
|
if (r < 0) continue; |
317 |
|
if (r >= nrows) break; |
318 |
< |
dy = (pr - (r+.5))/(m*rad); |
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 + dy*dy); |
326 |
> |
norm += warr[i++] = lookgauss(dx*dx + dy2); |
327 |
|
} |
328 |
|
} |
329 |
|
norm = 1.0/norm; |