| 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; |
