| 17 |
|
* each of NI2DIR surrounding directions. To speed this |
| 18 |
|
* calculation, we sort the data into half-planes and apply |
| 19 |
|
* simple tests to see which neighbor is closest in each |
| 20 |
< |
* direction. Once we have our approximate neighborhood |
| 20 |
> |
* angular slice. Once we have our approximate neighborhood |
| 21 |
|
* for a sample, we can use it in a modified Gaussian weighting |
| 22 |
|
* with allowing local anisotropy. Harmonic weighting is added |
| 23 |
|
* to reduce the influence of distant neighbors. This yields a |
| 90 |
|
{ |
| 91 |
|
if (mind <= 0) |
| 92 |
|
return; |
| 93 |
< |
if ((.998*ip->dmin <= mind) && (mind <= 1.002*ip->dmin)) |
| 93 |
> |
if ((.998*ip->dmin <= mind) & (mind <= 1.002*ip->dmin)) |
| 94 |
|
return; |
| 95 |
|
if (ip->da != NULL) { /* will need to recompute distribution */ |
| 96 |
|
free(ip->da); |
| 238 |
|
static double |
| 239 |
|
get_wt(const INTERP2 *ip, const int i, double x, double y) |
| 240 |
|
{ |
| 241 |
< |
double dir, rd, d2; |
| 241 |
> |
double dir, rd, r2, d2; |
| 242 |
|
int ri; |
| 243 |
|
/* get relative direction */ |
| 244 |
|
x -= ip->spt[i][0]; |
| 251 |
|
rd -= (double)ri; |
| 252 |
|
rd = (1.-rd)*ip->da[i][ri] + rd*ip->da[i][(ri+1)%NI2DIR]; |
| 253 |
|
rd = ip->smf * DECODE_DIA(ip, rd); |
| 254 |
+ |
r2 = 2.*rd*rd; |
| 255 |
|
d2 = x*x + y*y; |
| 256 |
+ |
if (d2 > 21.*r2) /* result would be < 1e-9 */ |
| 257 |
+ |
return(.0); |
| 258 |
|
/* Gaussian times harmonic weighting */ |
| 259 |
< |
return( exp(d2/(-2.*rd*rd)) * ip->dmin/(ip->dmin + sqrt(d2)) ); |
| 259 |
> |
return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) ); |
| 260 |
|
} |
| 261 |
|
|
| 262 |
|
/* Assign full set of normalized weights to interpolate the given position */ |
