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"; |
15 |
|
#include "color.h" |
16 |
|
|
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; |
226 |
> |
double t, num, denom, avg, wsum; |
227 |
> |
double maxm = orad/rad/CHECKRAD; |
228 |
> |
double mlimit[2]; |
229 |
|
int ilimit = 4/TEPS; |
230 |
|
register int i; |
231 |
|
/* iterative search for m */ |
232 |
+ |
mlimit[0] = 1.0; mlimit[1] = maxm; |
233 |
|
do { |
234 |
|
/* compute grey weighted average */ |
235 |
|
i = 1.25*CHECKRAD*rad*m + .5; |
241 |
|
avg += t*ringsum[i]; |
242 |
|
wsum += t*ringwt[i]; |
243 |
|
} |
244 |
+ |
if (avg < 1e-20) /* zero inclusive average */ |
245 |
+ |
return(1.0); |
246 |
|
avg /= wsum; |
247 |
|
/* check threshold */ |
248 |
< |
t = p0 - avg; |
249 |
< |
if (t < 0.0) t = -t; |
250 |
< |
t *= gausstable[0]/(m*m*avg); |
251 |
< |
if (t <= thresh && (m <= 1.0+FTINY || |
252 |
< |
(thresh-t)/thresh <= TEPS)) |
253 |
< |
break; |
254 |
< |
if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY || |
255 |
< |
(t-thresh)/thresh <= TEPS)) |
256 |
< |
break; |
248 |
> |
denom = m*m/gausstable[0] - p0/avg; |
249 |
> |
if (denom <= FTINY) /* zero exclusive average */ |
250 |
> |
return(maxm); |
251 |
> |
num = p0/avg - 1.0; |
252 |
> |
if (num < 0.0) num = -num; |
253 |
> |
t = num/denom; |
254 |
> |
if (t <= thresh) { |
255 |
> |
if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS) |
256 |
> |
break; |
257 |
> |
} else { |
258 |
> |
if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS) |
259 |
> |
break; |
260 |
> |
} |
261 |
> |
t = m; /* remember current m */ |
262 |
|
/* next guesstimate */ |
263 |
< |
m *= sqrt(t/thresh); |
264 |
< |
if (m < 1.0) m = 1.0; |
265 |
< |
else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD; |
263 |
> |
m = sqrt(gausstable[0]*(num/thresh + p0/avg)); |
264 |
> |
if (m < t) { /* bound it */ |
265 |
> |
if (m <= mlimit[0]+FTINY) |
266 |
> |
m = 0.5*(mlimit[0] + t); |
267 |
> |
mlimit[1] = t; |
268 |
> |
} else { |
269 |
> |
if (m >= mlimit[1]-FTINY) |
270 |
> |
m = 0.5*(mlimit[1] + t); |
271 |
> |
mlimit[0] = t; |
272 |
> |
} |
273 |
|
} while (--ilimit > 0); |
274 |
|
return(m); |
275 |
|
} |