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