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 |
|
|
41 |
|
extern int obarsize; /* size of output scan bar */ |
42 |
|
extern int orad; /* output window radius */ |
43 |
|
|
44 |
+ |
extern int wrapfilt; /* wrap filter horizontally? */ |
45 |
+ |
|
46 |
|
extern char *progname; |
47 |
|
|
48 |
|
float *gausstable; /* gauss lookup table */ |
52 |
|
short *ringndx; /* ring index table */ |
53 |
|
float *warr; /* array of pixel weights */ |
54 |
|
|
55 |
+ |
extern double (*ourbright)(); /* brightness computation function */ |
56 |
+ |
|
57 |
|
double pickfilt(); |
58 |
|
|
59 |
|
#define lookgauss(x) gausstable[(int)(10.*(x)+.5)] |
66 |
|
double d; |
67 |
|
register int x; |
68 |
|
|
69 |
< |
gtabsiz = 150*CHECKRAD; |
69 |
> |
gtabsiz = 111*CHECKRAD*CHECKRAD; |
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++) |
116 |
|
int wsum; |
117 |
|
double d; |
118 |
|
int y; |
119 |
< |
register int x; |
119 |
> |
register int x, offs; |
120 |
|
register COLOR *scan; |
121 |
|
|
122 |
|
wsum = 0; |
124 |
|
for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) { |
125 |
|
if (y < 0) continue; |
126 |
|
if (y >= yres) break; |
127 |
< |
d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent); |
127 |
> |
d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent); |
128 |
|
if (d < -0.5) continue; |
129 |
|
if (d >= 0.5) break; |
130 |
|
scan = scanin[y%barsize]; |
131 |
|
for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) { |
132 |
< |
if (x < 0) continue; |
133 |
< |
if (x >= xres) break; |
134 |
< |
d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent); |
132 |
> |
offs = x < 0 ? xres : x >= xres ? -xres : 0; |
133 |
> |
if (offs && !wrapfilt) |
134 |
> |
continue; |
135 |
> |
d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent); |
136 |
|
if (d < -0.5) continue; |
137 |
|
if (d >= 0.5) break; |
138 |
|
wsum++; |
139 |
< |
addcolor(csum, scan[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 |
|
|
154 |
|
double dy, dx, weight, wsum; |
155 |
|
COLOR ctmp; |
156 |
|
int y; |
157 |
< |
register int x; |
157 |
> |
register int x, offs; |
158 |
|
register COLOR *scan; |
159 |
|
|
160 |
|
wsum = FTINY; |
165 |
|
dy = (y_r*(y+.5) - (r+.5))/rad; |
166 |
|
scan = scanin[y%barsize]; |
167 |
|
for (x = xcent-xrad; x <= xcent+xrad; x++) { |
168 |
< |
if (x < 0) continue; |
169 |
< |
if (x >= xres) break; |
168 |
> |
offs = x < 0 ? xres : x >= xres ? -xres : 0; |
169 |
> |
if (offs && !wrapfilt) |
170 |
> |
continue; |
171 |
|
dx = (x_c*(x+.5) - (c+.5))/rad; |
172 |
|
weight = lookgauss(dx*dx + dy*dy); |
173 |
|
wsum += weight; |
174 |
< |
copycolor(ctmp, scan[x]); |
174 |
> |
copycolor(ctmp, scan[x+offs]); |
175 |
|
scalecolor(ctmp, weight); |
176 |
|
addcolor(csum, ctmp); |
177 |
|
} |
178 |
|
} |
179 |
< |
scalecolor(csum, 1.0/wsum); |
179 |
> |
weight = 1.0/wsum; |
180 |
> |
scalecolor(csum, weight); |
181 |
|
} |
182 |
|
|
183 |
|
|
186 |
|
int ccent, rcent; |
187 |
|
{ |
188 |
|
double d; |
189 |
< |
int r, y; |
189 |
> |
int r, y, offs; |
190 |
|
register int c, x; |
191 |
|
register float *gscan; |
183 |
– |
#define pval gscan |
192 |
|
/* compute ring sums */ |
193 |
|
bzero((char *)ringsum, (orad+1)*sizeof(float)); |
194 |
|
bzero((char *)ringwt, (orad+1)*sizeof(short)); |
197 |
|
if (rcent+r >= nrows) break; |
198 |
|
gscan = greybar[(rcent+r)%obarsize]; |
199 |
|
for (c = -orad; c <= orad; c++) { |
200 |
< |
if (ccent+c < 0) continue; |
201 |
< |
if (ccent+c >= ncols) break; |
200 |
> |
offs = ccent+c < 0 ? ncols : |
201 |
> |
ccent+c >= ncols ? -ncols : 0; |
202 |
> |
if (offs && !wrapfilt) |
203 |
> |
continue; |
204 |
|
x = ringndx[c*c + r*r]; |
205 |
|
if (x < 0) continue; |
206 |
< |
ringsum[x] += gscan[ccent+c]; |
206 |
> |
ringsum[x] += gscan[ccent+c+offs]; |
207 |
|
ringwt[x]++; |
208 |
|
} |
209 |
|
} |
211 |
|
for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) { |
212 |
|
if (y < 0) continue; |
213 |
|
if (y >= yres) break; |
214 |
< |
d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent); |
214 |
> |
d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent); |
215 |
|
if (d < -0.5) continue; |
216 |
|
if (d >= 0.5) break; |
217 |
|
for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) { |
218 |
< |
if (x < 0) continue; |
219 |
< |
if (x >= xres) break; |
220 |
< |
d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent); |
218 |
> |
offs = x < 0 ? xres : x >= xres ? -xres : 0; |
219 |
> |
if (offs && !wrapfilt) |
220 |
> |
continue; |
221 |
> |
d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent); |
222 |
|
if (d < -0.5) continue; |
223 |
|
if (d >= 0.5) break; |
224 |
< |
pval = scanin[y%barsize][x]; |
225 |
< |
sumans(x, y, rcent, ccent, pickfilt(bright(pval))); |
224 |
> |
sumans(x, y, rcent, ccent, |
225 |
> |
pickfilt((*ourbright)(scanin[y%barsize][x+offs]))); |
226 |
|
} |
227 |
|
} |
217 |
– |
#undef pval |
228 |
|
} |
229 |
|
|
230 |
|
|
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; |
298 |
> |
int ksiz, r, offs; |
299 |
|
double pc, pr, norm; |
300 |
|
register int i, c; |
301 |
|
register COLOR *scan; |
302 |
< |
|
303 |
< |
copycolor(pval, scanin[py%barsize][px]); |
302 |
> |
/* |
303 |
> |
* This normalization method fails at the picture borders because |
304 |
> |
* a different number of input pixels contribute there. |
305 |
> |
*/ |
306 |
> |
scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0); |
307 |
> |
copycolor(pval, scan[px]); |
308 |
|
pc = x_c*(px+.5); |
309 |
|
pr = y_r*(py+.5); |
310 |
|
ksiz = CHECKRAD*m*rad + 1; |
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 (c < 0) continue; |
322 |
< |
if (c >= ncols) break; |
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; |
336 |
|
if (r >= nrows) break; |
337 |
|
scan = scoutbar[r%obarsize]; |
338 |
|
for (c = ccent-ksiz; c <= ccent+ksiz; c++) { |
339 |
< |
if (c < 0) continue; |
340 |
< |
if (c >= ncols) break; |
339 |
> |
offs = c < 0 ? ncols : c >= ncols ? -ncols : 0; |
340 |
> |
if (offs && !wrapfilt) |
341 |
> |
continue; |
342 |
|
copycolor(ctmp, pval); |
343 |
|
dx = norm*warr[i++]; |
344 |
|
scalecolor(ctmp, dx); |
345 |
< |
addcolor(scan[c], ctmp); |
345 |
> |
addcolor(scan[c+offs], ctmp); |
346 |
|
} |
347 |
|
} |
348 |
|
} |