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 |
|
* |
9 |
|
|
10 |
|
#include "standard.h" |
11 |
|
|
12 |
+ |
#include <string.h> |
13 |
+ |
|
14 |
|
#include "color.h" |
15 |
|
|
16 |
+ |
#define RSCA 1.13 /* square-radius multiplier: sqrt(4/PI) */ |
17 |
|
#define TEPS 0.2 /* threshold proximity goal */ |
18 |
|
#define REPS 0.1 /* radius proximity goal */ |
19 |
|
|
68 |
|
double d; |
69 |
|
register int x; |
70 |
|
|
71 |
< |
gtabsiz = 150*CHECKRAD; |
71 |
> |
gtabsiz = 111*CHECKRAD*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); |
76 |
> |
gausstable[0] = exp(-d); |
77 |
|
for (x = 1; x < gtabsiz; x++) |
78 |
< |
if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0]) |
78 |
> |
if (x*0.1 <= d) |
79 |
|
gausstable[x] = gausstable[0]; |
80 |
+ |
else |
81 |
+ |
gausstable[x] = exp(-x*0.1); |
82 |
|
if (obarsize == 0) |
83 |
|
return; |
84 |
|
/* compute integral of filter */ |
85 |
< |
gaussN = PI*sqrt(d)*exp(-d); /* plateau */ |
86 |
< |
for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1) |
87 |
< |
gaussN += 0.1*2.0*PI*exp(-d*d); |
85 |
> |
gaussN = PI*d*exp(-d); /* plateau */ |
86 |
> |
for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1) |
87 |
> |
gaussN += 0.1*2.0*PI*d*exp(-d*d); |
88 |
|
/* normalize filter */ |
89 |
|
gaussN = x_c*y_r/(rad*rad*gaussN); |
90 |
|
for (x = 0; x < gtabsiz; x++) |
141 |
|
addcolor(csum, scan[x+offs]); |
142 |
|
} |
143 |
|
} |
144 |
< |
if (wsum > 1) |
145 |
< |
scalecolor(csum, 1.0/wsum); |
144 |
> |
if (wsum > 1) { |
145 |
> |
d = 1.0/wsum; |
146 |
> |
scalecolor(csum, d); |
147 |
> |
} |
148 |
|
} |
149 |
|
|
150 |
|
|
178 |
|
addcolor(csum, ctmp); |
179 |
|
} |
180 |
|
} |
181 |
< |
scalecolor(csum, 1.0/wsum); |
181 |
> |
weight = 1.0/wsum; |
182 |
> |
scalecolor(csum, weight); |
183 |
|
} |
184 |
|
|
185 |
|
|
192 |
|
register int c, x; |
193 |
|
register float *gscan; |
194 |
|
/* compute ring sums */ |
195 |
< |
bzero((char *)ringsum, (orad+1)*sizeof(float)); |
196 |
< |
bzero((char *)ringwt, (orad+1)*sizeof(short)); |
195 |
> |
memset((char *)ringsum, '\0', (orad+1)*sizeof(float)); |
196 |
> |
memset((char *)ringwt, '\0', (orad+1)*sizeof(short)); |
197 |
|
for (r = -orad; r <= orad; r++) { |
198 |
|
if (rcent+r < 0) continue; |
199 |
|
if (rcent+r >= nrows) break; |
237 |
|
double m = 1.0; |
238 |
|
double t, num, denom, avg, wsum; |
239 |
|
double mlimit[2]; |
240 |
< |
int ilimit = 4/TEPS; |
240 |
> |
int ilimit = 4.0/TEPS; |
241 |
|
register int i; |
242 |
|
/* iterative search for m */ |
243 |
|
mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD; |
244 |
|
do { |
245 |
|
/* compute grey weighted average */ |
246 |
< |
i = 1.25*CHECKRAD*rad*m + .5; |
246 |
> |
i = RSCA*CHECKRAD*rad*m + .5; |
247 |
|
if (i > orad) i = orad; |
248 |
|
avg = wsum = 0.0; |
249 |
|
while (i--) { |
295 |
|
int rcent, ccent; |
296 |
|
double m; |
297 |
|
{ |
298 |
< |
double dy, dx; |
298 |
> |
double dy2, dx; |
299 |
|
COLOR pval, ctmp; |
300 |
|
int ksiz, r, offs; |
301 |
|
double pc, pr, norm; |
317 |
|
for (r = rcent-ksiz; r <= rcent+ksiz; r++) { |
318 |
|
if (r < 0) continue; |
319 |
|
if (r >= nrows) break; |
320 |
< |
dy = (pr - (r+.5))/(m*rad); |
320 |
> |
dy2 = (pr - (r+.5))/(m*rad); |
321 |
> |
dy2 *= dy2; |
322 |
|
for (c = ccent-ksiz; c <= ccent+ksiz; c++) { |
323 |
|
if (!wrapfilt) { |
324 |
|
if (c < 0) continue; |
325 |
|
if (c >= ncols) break; |
326 |
|
} |
327 |
|
dx = (pc - (c+.5))/(m*rad); |
328 |
< |
norm += warr[i++] = lookgauss(dx*dx + dy*dy); |
328 |
> |
norm += warr[i++] = lookgauss(dx*dx + dy2); |
329 |
|
} |
330 |
|
} |
331 |
|
norm = 1.0/norm; |