| 84 |
|
return(ip); |
| 85 |
|
} |
| 86 |
|
|
| 87 |
+ |
/* Set minimum distance under which samples will start to merge */ |
| 88 |
+ |
void |
| 89 |
+ |
interp2_spacing(INTERP2 *ip, double mind) |
| 90 |
+ |
{ |
| 91 |
+ |
if (mind <= 0) |
| 92 |
+ |
return; |
| 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); |
| 97 |
+ |
ip->da = NULL; |
| 98 |
+ |
} |
| 99 |
+ |
ip->dmin = mind; |
| 100 |
+ |
} |
| 101 |
+ |
|
| 102 |
+ |
/* Modify smoothing parameter by the given factor */ |
| 103 |
+ |
void |
| 104 |
+ |
interp2_smooth(INTERP2 *ip, double sf) |
| 105 |
+ |
{ |
| 106 |
+ |
if ((ip->smf *= sf) < NI2DSMF) |
| 107 |
+ |
ip->smf = NI2DSMF; |
| 108 |
+ |
} |
| 109 |
+ |
|
| 110 |
|
/* private call-back to sort position index */ |
| 111 |
|
static int |
| 112 |
|
cmp_spos(const void *p1, const void *p2) |
| 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 */ |
