39 |
|
float dm; /* distance measure in this direction */ |
40 |
|
} SAMPORD; |
41 |
|
|
42 |
< |
/* Allocate a new set of interpolation samples */ |
42 |
> |
/* Allocate a new set of interpolation samples (caller assigns spt[] array) */ |
43 |
|
INTERP2 * |
44 |
|
interp2_alloc(int nsamps) |
45 |
|
{ |
60 |
|
return(nip); |
61 |
|
} |
62 |
|
|
63 |
+ |
/* Resize interpolation array (caller must assign any new values) */ |
64 |
+ |
INTERP2 * |
65 |
+ |
interp2_realloc(INTERP2 *ip, int nsamps) |
66 |
+ |
{ |
67 |
+ |
if (ip == NULL) |
68 |
+ |
return(interp2_alloc(nsamps)); |
69 |
+ |
if (nsamps <= 1) { |
70 |
+ |
interp2_free(ip); |
71 |
+ |
return(NULL); |
72 |
+ |
} |
73 |
+ |
if (nsamps == ip->ns); |
74 |
+ |
return(ip); |
75 |
+ |
if (ip->ra != NULL) { /* will need to recompute distribution */ |
76 |
+ |
free(ip->ra); |
77 |
+ |
ip->ra = NULL; |
78 |
+ |
} |
79 |
+ |
ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1)); |
80 |
+ |
if (ip == NULL) |
81 |
+ |
return(NULL); |
82 |
+ |
ip->ns = nsamps; |
83 |
+ |
return(ip); |
84 |
+ |
} |
85 |
+ |
|
86 |
|
/* private call-back to sort position index */ |
87 |
|
static int |
88 |
|
cmp_spos(const void *p1, const void *p2) |
97 |
|
return 0; |
98 |
|
} |
99 |
|
|
100 |
+ |
/* private routine to order samples in a particular direction */ |
101 |
+ |
static void |
102 |
+ |
sort_samples(SAMPORD *sord, const INTERP2 *ip, double ang) |
103 |
+ |
{ |
104 |
+ |
const double cosd = cos(ang); |
105 |
+ |
const double sind = sin(ang); |
106 |
+ |
int i; |
107 |
+ |
|
108 |
+ |
for (i = ip->ns; i--; ) { |
109 |
+ |
sord[i].si = i; |
110 |
+ |
sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1]; |
111 |
+ |
} |
112 |
+ |
qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos); |
113 |
+ |
} |
114 |
+ |
|
115 |
|
/* private routine to encode radius with range checks */ |
116 |
|
static int |
117 |
|
encode_radius(const INTERP2 *ip, double r) |
125 |
|
return(er); |
126 |
|
} |
127 |
|
|
128 |
< |
/* Compute anisotropic Gaussian basis function interpolant */ |
129 |
< |
static int |
130 |
< |
interp2_compute(INTERP2 *ip) |
128 |
> |
/* (Re)compute anisotropic basis function interpolant (normally automatic) */ |
129 |
> |
int |
130 |
> |
interp2_analyze(INTERP2 *ip) |
131 |
|
{ |
132 |
|
SAMPORD *sortord; |
133 |
< |
int *rightrndx, *leftrndx; |
133 |
> |
int *rightrndx, *leftrndx, *endrndx; |
134 |
|
int bd; |
135 |
|
/* sanity checks */ |
136 |
|
if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0)) |
146 |
|
sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns); |
147 |
|
rightrndx = (int *)malloc(sizeof(int)*ip->ns); |
148 |
|
leftrndx = (int *)malloc(sizeof(int)*ip->ns); |
149 |
< |
if ((sortord == NULL) | (rightrndx == NULL) | (leftrndx == NULL)) |
149 |
> |
endrndx = (int *)malloc(sizeof(int)*ip->ns); |
150 |
> |
if ((sortord == NULL) | (rightrndx == NULL) | |
151 |
> |
(leftrndx == NULL) | (endrndx == NULL)) |
152 |
|
return(0); |
153 |
|
/* run through bidirections */ |
154 |
|
for (bd = 0; bd < NI2DIR/2; bd++) { |
155 |
|
const double ang = 2.*PI/NI2DIR*bd; |
156 |
< |
double cosd, sind; |
156 |
> |
int *sptr; |
157 |
|
int i; |
158 |
|
/* create right reverse index */ |
159 |
< |
if (bd) { /* re-use from prev. iteration? */ |
160 |
< |
int *sptr = rightrndx; |
159 |
> |
if (bd) { /* re-use from previous iteration? */ |
160 |
> |
sptr = rightrndx; |
161 |
|
rightrndx = leftrndx; |
162 |
|
leftrndx = sptr; |
163 |
< |
} else { /* else compute it */ |
164 |
< |
cosd = cos(ang + (PI/2. - PI/NI2DIR)); |
165 |
< |
sind = sin(ang + (PI/2. - PI/NI2DIR)); |
126 |
< |
for (i = 0; i < ip->ns; i++) { |
127 |
< |
sortord[i].si = i; |
128 |
< |
sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1]; |
129 |
< |
} |
130 |
< |
qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos); |
131 |
< |
for (i = 0; i < ip->ns; i++) |
163 |
> |
} else { /* else sort first half-plane */ |
164 |
> |
sort_samples(sortord, ip, PI/2. - PI/NI2DIR); |
165 |
> |
for (i = ip->ns; i--; ) |
166 |
|
rightrndx[sortord[i].si] = i; |
167 |
+ |
/* & store reverse order for later */ |
168 |
+ |
for (i = ip->ns; i--; ) |
169 |
+ |
endrndx[sortord[i].si] = ip->ns-1 - i; |
170 |
|
} |
171 |
|
/* create new left reverse index */ |
172 |
< |
cosd = cos(ang + (PI/2. + PI/NI2DIR)); |
173 |
< |
sind = sin(ang + (PI/2. + PI/NI2DIR)); |
174 |
< |
for (i = 0; i < ip->ns; i++) { |
175 |
< |
sortord[i].si = i; |
176 |
< |
sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1]; |
177 |
< |
} |
178 |
< |
qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos); |
142 |
< |
for (i = 0; i < ip->ns; i++) |
172 |
> |
if (bd == NI2DIR/2 - 1) { /* use order from first iteration? */ |
173 |
> |
sptr = leftrndx; |
174 |
> |
leftrndx = endrndx; |
175 |
> |
endrndx = sptr; |
176 |
> |
} else { /* else compute new half-plane */ |
177 |
> |
sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR)); |
178 |
> |
for (i = ip->ns; i--; ) |
179 |
|
leftrndx[sortord[i].si] = i; |
144 |
– |
/* sort grid values in this direction */ |
145 |
– |
cosd = cos(ang); |
146 |
– |
sind = sin(ang); |
147 |
– |
for (i = 0; i < ip->ns; i++) { |
148 |
– |
sortord[i].si = i; |
149 |
– |
sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1]; |
180 |
|
} |
181 |
< |
qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos); |
181 |
> |
/* sort grid values in this direction */ |
182 |
> |
sort_samples(sortord, ip, ang); |
183 |
|
/* find nearest neighbors each side */ |
184 |
< |
for (i = 0; i < ip->ns; i++) { |
184 |
> |
for (i = ip->ns; i--; ) { |
185 |
|
const int rpos = rightrndx[sortord[i].si]; |
186 |
|
const int lpos = leftrndx[sortord[i].si]; |
187 |
|
int j; |
207 |
|
free(sortord); /* clean up */ |
208 |
|
free(rightrndx); |
209 |
|
free(leftrndx); |
210 |
+ |
free(endrndx); |
211 |
|
return(1); |
212 |
|
} |
213 |
|
|
242 |
|
if ((wtv == NULL) | (ip == NULL)) |
243 |
|
return(0); |
244 |
|
/* need to compute interpolant? */ |
245 |
< |
if (ip->ra == NULL && !interp2_compute(ip)) |
245 |
> |
if (ip->ra == NULL && !interp2_analyze(ip)) |
246 |
|
return(0); |
247 |
|
|
248 |
|
wnorm = 0; /* compute raw weights */ |
276 |
|
if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL)) |
277 |
|
return(0); |
278 |
|
/* need to compute interpolant? */ |
279 |
< |
if (ip->ra == NULL && !interp2_compute(ip)) |
279 |
> |
if (ip->ra == NULL && !interp2_analyze(ip)) |
280 |
|
return(0); |
281 |
|
/* identify top n weights */ |
282 |
|
for (i = ip->ns; i--; ) { |