1 |
< |
/* Copyright (c) 1992 Regents of the University of California */ |
1 |
> |
/* Copyright (c) 1993 Regents of the University of California */ |
2 |
|
|
3 |
|
#ifndef lint |
4 |
|
static char SCCSid[] = "$SunId$ LBL"; |
14 |
|
|
15 |
|
#include "color.h" |
16 |
|
|
17 |
< |
#define TEPS 0.25 /* threshold proximity goal */ |
17 |
> |
#define TEPS 0.2 /* threshold proximity goal */ |
18 |
> |
#define REPS 0.1 /* radius proximity goal */ |
19 |
|
|
20 |
|
extern double CHECKRAD; /* radius over which gaussian is summed */ |
21 |
|
|
223 |
|
double p0; |
224 |
|
{ |
225 |
|
double m = 1.0; |
226 |
< |
double t, avg, wsum; |
227 |
< |
int ilimit = 5/TEPS; |
226 |
> |
double t, num, denom, avg, wsum; |
227 |
> |
double mlimit[2]; |
228 |
> |
int ilimit = 4/TEPS; |
229 |
|
register int i; |
230 |
|
/* iterative search for m */ |
231 |
+ |
mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD; |
232 |
|
do { |
233 |
|
/* compute grey weighted average */ |
234 |
|
i = 1.25*CHECKRAD*rad*m + .5; |
240 |
|
avg += t*ringsum[i]; |
241 |
|
wsum += t*ringwt[i]; |
242 |
|
} |
243 |
+ |
if (avg < 1e-20) /* zero inclusive average */ |
244 |
+ |
return(1.0); |
245 |
|
avg /= wsum; |
246 |
|
/* check threshold */ |
247 |
< |
t = p0 - avg; |
248 |
< |
if (t < 0.0) t = -t; |
249 |
< |
t *= gausstable[0]/(m*m*avg); |
250 |
< |
if (t <= thresh && (m <= 1.0+FTINY || |
251 |
< |
(thresh-t)/thresh <= TEPS)) |
252 |
< |
break; |
253 |
< |
if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY || |
254 |
< |
(t-thresh)/thresh <= TEPS)) |
255 |
< |
break; |
247 |
> |
denom = m*m/gausstable[0] - p0/avg; |
248 |
> |
if (denom <= FTINY) { /* zero exclusive average */ |
249 |
> |
if (m >= mlimit[1]-REPS) |
250 |
> |
break; |
251 |
> |
m = mlimit[1]; |
252 |
> |
continue; |
253 |
> |
} |
254 |
> |
num = p0/avg - 1.0; |
255 |
> |
if (num < 0.0) num = -num; |
256 |
> |
t = num/denom; |
257 |
> |
if (t <= thresh) { |
258 |
> |
if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS) |
259 |
> |
break; |
260 |
> |
} else { |
261 |
> |
if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS) |
262 |
> |
break; |
263 |
> |
} |
264 |
> |
t = m; /* remember current m */ |
265 |
|
/* next guesstimate */ |
266 |
< |
m *= sqrt(t/thresh); |
267 |
< |
if (m < 1.0) m = 1.0; |
268 |
< |
else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD; |
266 |
> |
m = sqrt(gausstable[0]*(num/thresh + p0/avg)); |
267 |
> |
if (m < t) { /* bound it */ |
268 |
> |
if (m <= mlimit[0]+FTINY) |
269 |
> |
m = 0.5*(mlimit[0] + t); |
270 |
> |
mlimit[1] = t; |
271 |
> |
} else { |
272 |
> |
if (m >= mlimit[1]-FTINY) |
273 |
> |
m = 0.5*(mlimit[1] + t); |
274 |
> |
mlimit[0] = t; |
275 |
> |
} |
276 |
|
} while (--ilimit > 0); |
277 |
|
return(m); |
278 |
|
} |