1 |
|
|
/* |
2 |
|
|
Teem: Tools to process and visualize scientific data and images . |
3 |
|
|
Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago |
4 |
|
|
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann |
5 |
|
|
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah |
6 |
|
|
|
7 |
|
|
This library is free software; you can redistribute it and/or |
8 |
|
|
modify it under the terms of the GNU Lesser General Public License |
9 |
|
|
(LGPL) as published by the Free Software Foundation; either |
10 |
|
|
version 2.1 of the License, or (at your option) any later version. |
11 |
|
|
The terms of redistributing and/or modifying this software also |
12 |
|
|
include exceptions to the LGPL that facilitate static linking. |
13 |
|
|
|
14 |
|
|
This library is distributed in the hope that it will be useful, |
15 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
17 |
|
|
Lesser General Public License for more details. |
18 |
|
|
|
19 |
|
|
You should have received a copy of the GNU Lesser General Public License |
20 |
|
|
along with this library; if not, write to Free Software Foundation, Inc., |
21 |
|
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
22 |
|
|
*/ |
23 |
|
|
|
24 |
|
|
|
25 |
|
|
#include "ten.h" |
26 |
|
|
#include "privateTen.h" |
27 |
|
|
|
28 |
|
|
tenInterpParm * |
29 |
|
|
tenInterpParmNew(void) { |
30 |
|
|
tenInterpParm *tip; |
31 |
|
|
|
32 |
|
32 |
tip = AIR_CAST(tenInterpParm *, malloc(sizeof(tenInterpParm))); |
33 |
✓✗ |
16 |
if (tip) { |
34 |
|
16 |
tip->verbose = AIR_FALSE; |
35 |
|
16 |
tip->convStep = 0.2; |
36 |
|
16 |
tip->minNorm = 0.0; |
37 |
|
16 |
tip->convEps = 0.0000000001; |
38 |
|
16 |
tip->wghtSumEps = 0.0000001; |
39 |
|
16 |
tip->enableRecurse = AIR_TRUE; |
40 |
|
16 |
tip->maxIter = 20; |
41 |
|
16 |
tip->numSteps = 100; |
42 |
|
16 |
tip->lengthFancy = AIR_FALSE; |
43 |
|
|
|
44 |
|
16 |
tip->allocLen = 0; |
45 |
|
16 |
tip->eval = NULL; |
46 |
|
16 |
tip->evec = NULL; |
47 |
|
16 |
tip->rtIn = NULL; |
48 |
|
16 |
tip->rtLog = NULL; |
49 |
|
16 |
tip->qIn = NULL; |
50 |
|
16 |
tip->qBuff = NULL; |
51 |
|
16 |
tip->qInter = NULL; |
52 |
|
|
|
53 |
|
16 |
tip->numIter = 0; |
54 |
|
16 |
tip->convFinal = AIR_NAN; |
55 |
|
16 |
tip->lengthShape = AIR_NAN; |
56 |
|
16 |
tip->lengthOrient = AIR_NAN; |
57 |
|
16 |
} |
58 |
|
16 |
return tip; |
59 |
|
|
} |
60 |
|
|
|
61 |
|
|
|
62 |
|
|
|
63 |
|
|
/* |
64 |
|
|
** handles allocating all the various buffers that are needed for QGL |
65 |
|
|
** interpolation, so that they are repeatedly allocated and freed |
66 |
|
|
** between calls |
67 |
|
|
*/ |
68 |
|
|
int |
69 |
|
|
tenInterpParmBufferAlloc(tenInterpParm *tip, unsigned int num) { |
70 |
|
|
static const char me[]="tenInterpParmBufferAlloc"; |
71 |
|
|
|
72 |
✗✓ |
32 |
if (0 == num) { |
73 |
|
|
/* user wants to free buffers for some reason */ |
74 |
|
|
airFree(tip->eval); tip->eval = NULL; |
75 |
|
|
airFree(tip->evec); tip->evec = NULL; |
76 |
|
|
airFree(tip->rtIn); tip->rtIn = NULL; |
77 |
|
|
airFree(tip->rtLog); tip->rtLog = NULL; |
78 |
|
|
airFree(tip->qIn); tip->qIn = NULL; |
79 |
|
|
airFree(tip->qBuff); tip->qBuff = NULL; |
80 |
|
|
airFree(tip->qInter); tip->qInter = NULL; |
81 |
|
|
tip->allocLen = 0; |
82 |
✗✓ |
16 |
} else if (1 == num) { |
83 |
|
|
biffAddf(TEN, "%s: need num >= 2 (not %u)", me, num); |
84 |
|
|
return 1; |
85 |
✓✗ |
16 |
} else if (num != tip->allocLen) { |
86 |
|
16 |
airFree(tip->eval); tip->eval = NULL; |
87 |
|
16 |
airFree(tip->evec); tip->evec = NULL; |
88 |
|
16 |
airFree(tip->rtIn); tip->rtIn = NULL; |
89 |
|
16 |
airFree(tip->rtLog); tip->rtLog = NULL; |
90 |
|
16 |
airFree(tip->qIn); tip->qIn = NULL; |
91 |
|
16 |
airFree(tip->qBuff); tip->qBuff = NULL; |
92 |
|
16 |
airFree(tip->qInter); tip->qInter = NULL; |
93 |
|
16 |
tip->eval = AIR_CALLOC(3*num, double); |
94 |
|
16 |
tip->evec = AIR_CALLOC(9*num, double); |
95 |
|
16 |
tip->rtIn = AIR_CALLOC(3*num, double); |
96 |
|
16 |
tip->rtLog = AIR_CALLOC(3*num, double); |
97 |
|
16 |
tip->qIn = AIR_CALLOC(4*num, double); |
98 |
|
16 |
tip->qBuff = AIR_CALLOC(4*num, double); |
99 |
|
16 |
tip->qInter = AIR_CALLOC(num*num, double); |
100 |
✓✗✓✗
✗✓ |
48 |
if (!(tip->eval && tip->evec && |
101 |
✓✗✓✗
|
32 |
tip->rtIn && tip->rtLog && |
102 |
✓✗✓✗
|
48 |
tip->qIn && tip->qBuff && tip->qInter)) { |
103 |
|
|
biffAddf(TEN, "%s: didn't alloc buffers (%p,%p,%p %p %p %p %p)", me, |
104 |
|
|
AIR_VOIDP(tip->eval), AIR_VOIDP(tip->evec), |
105 |
|
|
AIR_VOIDP(tip->rtIn), AIR_VOIDP(tip->rtLog), |
106 |
|
|
AIR_VOIDP(tip->qIn), AIR_VOIDP(tip->qBuff), |
107 |
|
|
AIR_VOIDP(tip->qInter)); |
108 |
|
|
return 1; |
109 |
|
|
} |
110 |
|
16 |
tip->allocLen = num; |
111 |
|
16 |
} |
112 |
|
16 |
return 0; |
113 |
|
16 |
} |
114 |
|
|
|
115 |
|
|
tenInterpParm * |
116 |
|
|
tenInterpParmCopy(tenInterpParm *tip) { |
117 |
|
|
static const char me[]="tenInterpParmCopy"; |
118 |
|
|
tenInterpParm *newtip; |
119 |
|
|
unsigned int num; |
120 |
|
|
|
121 |
|
16 |
num = tip->allocLen; |
122 |
|
8 |
newtip = tenInterpParmNew(); |
123 |
✓✗ |
8 |
if (newtip) { |
124 |
|
8 |
memcpy(newtip, tip, sizeof(tenInterpParm)); |
125 |
|
|
/* manually set all pointers */ |
126 |
|
8 |
newtip->allocLen = 0; |
127 |
|
8 |
newtip->eval = NULL; |
128 |
|
8 |
newtip->evec = NULL; |
129 |
|
8 |
newtip->rtIn = NULL; |
130 |
|
8 |
newtip->rtLog = NULL; |
131 |
|
8 |
newtip->qIn = NULL; |
132 |
|
8 |
newtip->qBuff = NULL; |
133 |
|
8 |
newtip->qInter = NULL; |
134 |
✗✓ |
8 |
if (tenInterpParmBufferAlloc(newtip, num)) { |
135 |
|
|
biffAddf(TEN, "%s: trouble allocating output", me); |
136 |
|
|
return NULL; |
137 |
|
|
} |
138 |
|
8 |
memcpy(newtip->eval, tip->eval, 3*num*sizeof(double)); |
139 |
|
8 |
memcpy(newtip->evec, tip->evec, 9*num*sizeof(double)); |
140 |
|
8 |
memcpy(newtip->rtIn, tip->rtIn, 3*num*sizeof(double)); |
141 |
|
8 |
memcpy(newtip->rtLog, tip->rtLog, 3*num*sizeof(double)); |
142 |
|
8 |
memcpy(newtip->qIn, tip->qIn, 4*num*sizeof(double)); |
143 |
|
8 |
memcpy(newtip->qBuff, tip->qBuff, 4*num*sizeof(double)); |
144 |
|
8 |
memcpy(newtip->qInter, tip->qInter, num*num*sizeof(double)); |
145 |
|
8 |
} |
146 |
|
8 |
return newtip; |
147 |
|
8 |
} |
148 |
|
|
|
149 |
|
|
tenInterpParm * |
150 |
|
|
tenInterpParmNix(tenInterpParm *tip) { |
151 |
|
|
|
152 |
✓✗ |
32 |
if (tip) { |
153 |
|
16 |
airFree(tip->eval); |
154 |
|
16 |
airFree(tip->evec); |
155 |
|
16 |
airFree(tip->rtIn); |
156 |
|
16 |
airFree(tip->rtLog); |
157 |
|
16 |
airFree(tip->qIn); |
158 |
|
16 |
airFree(tip->qBuff); |
159 |
|
16 |
airFree(tip->qInter); |
160 |
|
16 |
free(tip); |
161 |
|
16 |
} |
162 |
|
16 |
return NULL; |
163 |
|
|
} |
164 |
|
|
|
165 |
|
|
/* |
166 |
|
|
******** tenInterpTwo_d |
167 |
|
|
** |
168 |
|
|
** interpolates between two tensors, in various ways |
169 |
|
|
** |
170 |
|
|
** this is really only used for demo purposes; its not useful for |
171 |
|
|
** doing real work in DTI fields. So: its okay that its slow |
172 |
|
|
** (e.g. for tenInterpTypeQuatGeoLox{K,R}, it recomputes the |
173 |
|
|
** eigensystems at the endpoints for every call, even though they are |
174 |
|
|
** apt to be the same between calls. |
175 |
|
|
** |
176 |
|
|
** this |
177 |
|
|
*/ |
178 |
|
|
void |
179 |
|
|
tenInterpTwo_d(double oten[7], |
180 |
|
|
const double tenA[7], const double tenB[7], |
181 |
|
|
int ptype, double aa, |
182 |
|
|
tenInterpParm *tip) { |
183 |
|
|
static const char me[]="tenInterpTwo_d"; |
184 |
|
|
double logA[7], logB[7], tmp1[7], tmp2[7], logMean[7], |
185 |
|
|
mat1[9], mat2[9], mat3[9], sqrtA[7], isqrtA[7], |
186 |
|
|
mean[7], sqrtB[7], isqrtB[7], |
187 |
|
|
oeval[3], evalA[3], evalB[3], oevec[9], evecA[9], evecB[9]; |
188 |
|
|
|
189 |
|
|
|
190 |
|
|
if (!( oten && tenA && tenB )) { |
191 |
|
|
/* got NULL pointer, but not using biff */ |
192 |
|
|
if (oten) { |
193 |
|
|
TEN_T_SET(oten, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN, |
194 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
195 |
|
|
} |
196 |
|
|
return; |
197 |
|
|
} |
198 |
|
|
|
199 |
|
|
switch (ptype) { |
200 |
|
|
case tenInterpTypeLinear: |
201 |
|
|
TEN_T_LERP(oten, aa, tenA, tenB); |
202 |
|
|
break; |
203 |
|
|
case tenInterpTypeLogLinear: |
204 |
|
|
tenLogSingle_d(logA, tenA); |
205 |
|
|
tenLogSingle_d(logB, tenB); |
206 |
|
|
TEN_T_LERP(logMean, aa, logA, logB); |
207 |
|
|
tenExpSingle_d(oten, logMean); |
208 |
|
|
break; |
209 |
|
|
case tenInterpTypeAffineInvariant: |
210 |
|
|
tenSqrtSingle_d(sqrtA, tenA); |
211 |
|
|
tenInv_d(isqrtA, sqrtA); |
212 |
|
|
TEN_T2M(mat1, tenB); |
213 |
|
|
TEN_T2M(mat2, isqrtA); |
214 |
|
|
ELL_3M_MUL(mat3, mat1, mat2); /* B * is(A) */ |
215 |
|
|
ELL_3M_MUL(mat1, mat2, mat3); /* is(A) * B * is(A) */ |
216 |
|
|
TEN_M2T(tmp2, mat1); |
217 |
|
|
tenPowSingle_d(tmp1, tmp2, aa); /* m = (is(A) * B * is(A))^aa */ |
218 |
|
|
TEN_T2M(mat1, tmp1); |
219 |
|
|
TEN_T2M(mat2, sqrtA); |
220 |
|
|
ELL_3M_MUL(mat3, mat1, mat2); /* m * sqrt(A) */ |
221 |
|
|
ELL_3M_MUL(mat1, mat2, mat3); /* sqrt(A) * m * sqrt(A) */ |
222 |
|
|
TEN_M2T(oten, mat1); |
223 |
|
|
oten[0] = AIR_LERP(aa, tenA[0], tenB[0]); |
224 |
|
|
if (tip->verbose) { |
225 |
|
|
fprintf(stderr, "%s:\nA= %g %g %g %g %g %g\n" |
226 |
|
|
"B = %g %g %g %g %g %g\n" |
227 |
|
|
"foo = %g %g %g %g %g %g\n" |
228 |
|
|
"bar(%g) = %g %g %g %g %g %g\n", me, |
229 |
|
|
tenA[1], tenA[2], tenA[3], tenA[4], tenA[5], tenA[6], |
230 |
|
|
tenB[1], tenB[2], tenB[3], tenB[4], tenB[5], tenB[6], |
231 |
|
|
tmp1[1], tmp1[2], tmp1[3], tmp1[4], tmp1[5], tmp1[6], |
232 |
|
|
aa, oten[1], oten[2], oten[3], oten[4], oten[5], oten[6]); |
233 |
|
|
} |
234 |
|
|
break; |
235 |
|
|
case tenInterpTypeWang: |
236 |
|
|
/* HEY: this seems to be broken */ |
237 |
|
|
TEN_T_LERP(mean, aa, tenA, tenB); /* "A" = mean */ |
238 |
|
|
tenLogSingle_d(logA, tenA); |
239 |
|
|
tenLogSingle_d(logB, tenB); |
240 |
|
|
TEN_T_LERP(logMean, aa, logA, logB); /* "B" = logMean */ |
241 |
|
|
tenSqrtSingle_d(sqrtB, logMean); |
242 |
|
|
tenInv_d(isqrtB, sqrtB); |
243 |
|
|
TEN_T2M(mat1, mean); |
244 |
|
|
TEN_T2M(mat2, isqrtB); |
245 |
|
|
ELL_3M_MUL(mat3, mat1, mat2); |
246 |
|
|
ELL_3M_MUL(mat1, mat2, mat3); |
247 |
|
|
TEN_M2T(tmp1, mat1); |
248 |
|
|
tenSqrtSingle_d(oten, tmp1); |
249 |
|
|
oten[0] = AIR_LERP(aa, tenA[0], tenB[0]); |
250 |
|
|
break; |
251 |
|
|
case tenInterpTypeQuatGeoLoxK: |
252 |
|
|
case tenInterpTypeQuatGeoLoxR: |
253 |
|
|
tenEigensolve_d(evalA, evecA, tenA); |
254 |
|
|
tenEigensolve_d(evalB, evecB, tenB); |
255 |
|
|
if (tenInterpTypeQuatGeoLoxK == ptype) { |
256 |
|
|
tenQGLInterpTwoEvalK(oeval, evalA, evalB, aa); |
257 |
|
|
} else { |
258 |
|
|
tenQGLInterpTwoEvalR(oeval, evalA, evalB, aa); |
259 |
|
|
} |
260 |
|
|
tenQGLInterpTwoEvec(oevec, evecA, evecB, aa); |
261 |
|
|
tenMakeSingle_d(oten, AIR_LERP(aa, tenA[0], tenB[0]), oeval, oevec); |
262 |
|
|
break; |
263 |
|
|
case tenInterpTypeGeoLoxK: |
264 |
|
|
case tenInterpTypeGeoLoxR: |
265 |
|
|
case tenInterpTypeLoxK: |
266 |
|
|
case tenInterpTypeLoxR: |
267 |
|
|
/* (currently) no closed-form expression for these */ |
268 |
|
|
TEN_T_SET(oten, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN, |
269 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
270 |
|
|
break; |
271 |
|
|
case tenInterpTypeRThetaPhiLinear: |
272 |
|
|
if (1) { |
273 |
|
|
double rtpA[3], rtpB[3], rtpM[3], eval[3], tenM[7]; |
274 |
|
|
tenEigensolve_d(eval, NULL, tenA); |
275 |
|
|
tenTripleConvertSingle_d(rtpA, tenTripleTypeRThetaPhi, |
276 |
|
|
eval, tenTripleTypeEigenvalue); |
277 |
|
|
tenEigensolve_d(eval, NULL, tenB); |
278 |
|
|
tenTripleConvertSingle_d(rtpB, tenTripleTypeRThetaPhi, |
279 |
|
|
eval, tenTripleTypeEigenvalue); |
280 |
|
|
TEN_T_LERP(tenM, aa, tenA, tenB); |
281 |
|
|
tenEigensolve_d(eval, oevec, tenM); |
282 |
|
|
ELL_3V_LERP(rtpM, aa, rtpA, rtpB); |
283 |
|
|
tenTripleConvertSingle_d(oeval, tenTripleTypeEigenvalue, |
284 |
|
|
rtpM, tenTripleTypeRThetaPhi); |
285 |
|
|
} |
286 |
|
|
tenMakeSingle_d(oten, AIR_LERP(aa, tenA[0], tenB[0]), oeval, oevec); |
287 |
|
|
break; |
288 |
|
|
default: |
289 |
|
|
/* error */ |
290 |
|
|
TEN_T_SET(oten, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN, |
291 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
292 |
|
|
break; |
293 |
|
|
} |
294 |
|
|
return; |
295 |
|
|
} |
296 |
|
|
|
297 |
|
|
/* |
298 |
|
|
** this NON-optionally uses biff |
299 |
|
|
** |
300 |
|
|
** for simplicity, a pre-allocated tenInterpParm MUST be passed, |
301 |
|
|
** regardless of the interpolation requested |
302 |
|
|
*/ |
303 |
|
|
int |
304 |
|
|
tenInterpN_d(double tenOut[7], |
305 |
|
|
const double *tenIn, const double *wght, |
306 |
|
|
unsigned int num, int ptype, tenInterpParm *tip) { |
307 |
|
|
static const char me[]="tenInterpN_d"; |
308 |
|
|
unsigned int ii; |
309 |
|
|
double ww, cc, tenErr[7], tmpTen[7], wghtSum, eval[3], evec[9], rtp[3]; |
310 |
|
|
|
311 |
|
|
TEN_T_SET(tenErr, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN, |
312 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
313 |
|
|
/* wght can be NULL ==> equal 1/num weight for all */ |
314 |
|
|
if (!(tenOut && tenIn && tip)) { |
315 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
316 |
|
|
return 1; |
317 |
|
|
} |
318 |
|
|
if (!( num >= 2 )) { |
319 |
|
|
biffAddf(TEN, "%s: need num >= 2 (not %u)", me, num); |
320 |
|
|
TEN_T_COPY(tenOut, tenErr); return 1; |
321 |
|
|
} |
322 |
|
|
if (airEnumValCheck(tenInterpType, ptype)) { |
323 |
|
|
biffAddf(TEN, "%s: invalid %s %d", me, tenInterpType->name, ptype); |
324 |
|
|
TEN_T_COPY(tenOut, tenErr); return 1; |
325 |
|
|
} |
326 |
|
|
wghtSum = 0; |
327 |
|
|
for (ii=0; ii<num; ii++) { |
328 |
|
|
ww = wght ? wght[ii] : 1.0/num; |
329 |
|
|
wghtSum += ww; |
330 |
|
|
} |
331 |
|
|
if (!( AIR_IN_CL(1 - tip->wghtSumEps, wghtSum, 1 + tip->wghtSumEps) )) { |
332 |
|
|
biffAddf(TEN, "%s: wght sum %g not within %g of 1.0", me, |
333 |
|
|
wghtSum, tip->wghtSumEps); |
334 |
|
|
TEN_T_COPY(tenOut, tenErr); return 1; |
335 |
|
|
} |
336 |
|
|
|
337 |
|
|
switch (ptype) { |
338 |
|
|
case tenInterpTypeLinear: |
339 |
|
|
TEN_T_SET(tenOut, 0, 0, 0, 0, 0, 0, 0); |
340 |
|
|
cc = 0; |
341 |
|
|
for (ii=0; ii<num; ii++) { |
342 |
|
|
ww = wght ? wght[ii] : 1.0/num; |
343 |
|
|
TEN_T_SCALE_INCR(tenOut, ww, tenIn + 7*ii); |
344 |
|
|
cc += ww*(tenIn + 7*ii)[0]; |
345 |
|
|
} |
346 |
|
|
tenOut[0] = cc; |
347 |
|
|
break; |
348 |
|
|
case tenInterpTypeLogLinear: |
349 |
|
|
TEN_T_SET(tenOut, 0, 0, 0, 0, 0, 0, 0); |
350 |
|
|
cc = 0; |
351 |
|
|
for (ii=0; ii<num; ii++) { |
352 |
|
|
ww = wght ? wght[ii] : 1.0/num; |
353 |
|
|
tenLogSingle_d(tmpTen, tenIn + 7*ii); |
354 |
|
|
TEN_T_SCALE_INCR(tenOut, ww, tmpTen); |
355 |
|
|
cc += ww*(tenIn + 7*ii)[0]; |
356 |
|
|
} |
357 |
|
|
tenOut[0] = cc; |
358 |
|
|
TEN_T_COPY(tmpTen, tenOut); |
359 |
|
|
tenExpSingle_d(tenOut, tmpTen); |
360 |
|
|
break; |
361 |
|
|
case tenInterpTypeAffineInvariant: |
362 |
|
|
case tenInterpTypeWang: |
363 |
|
|
biffAddf(TEN, "%s: sorry, not implemented", me); |
364 |
|
|
TEN_T_COPY(tenOut, tenErr); return 1; |
365 |
|
|
break; |
366 |
|
|
case tenInterpTypeGeoLoxK: |
367 |
|
|
case tenInterpTypeGeoLoxR: |
368 |
|
|
case tenInterpTypeLoxK: |
369 |
|
|
case tenInterpTypeLoxR: |
370 |
|
|
biffAddf(TEN, "%s: %s doesn't support averaging multiple tensors", me, |
371 |
|
|
airEnumStr(tenInterpType, ptype)); |
372 |
|
|
TEN_T_COPY(tenOut, tenErr); return 1; |
373 |
|
|
break; |
374 |
|
|
case tenInterpTypeQuatGeoLoxK: |
375 |
|
|
case tenInterpTypeQuatGeoLoxR: |
376 |
|
|
if (tenInterpParmBufferAlloc(tip, num)) { |
377 |
|
|
biffAddf(TEN, "%s: trouble getting buffers", me); |
378 |
|
|
TEN_T_COPY(tenOut, tenErr); return 1; |
379 |
|
|
} else { |
380 |
|
|
cc = 0; |
381 |
|
|
for (ii=0; ii<num; ii++) { |
382 |
|
|
tenEigensolve_d(tip->eval + 3*ii, tip->evec + 9*ii, tenIn + 7*ii); |
383 |
|
|
ww = wght ? wght[ii] : 1.0/num; |
384 |
|
|
cc += ww*(tenIn + 7*ii)[0]; |
385 |
|
|
} |
386 |
|
|
if (_tenQGLInterpNEval(eval, tip->eval, wght, num, ptype, tip) |
387 |
|
|
|| _tenQGLInterpNEvec(evec, tip->evec, wght, num, tip)) { |
388 |
|
|
biffAddf(TEN, "%s: trouble computing", me); |
389 |
|
|
TEN_T_COPY(tenOut, tenErr); return 1; |
390 |
|
|
} |
391 |
|
|
tenMakeSingle_d(tenOut, cc, eval, evec); |
392 |
|
|
} |
393 |
|
|
break; |
394 |
|
|
case tenInterpTypeRThetaPhiLinear: |
395 |
|
|
TEN_T_SET(tmpTen, 0, 0, 0, 0, 0, 0, 0); |
396 |
|
|
ELL_3V_SET(rtp, 0, 0, 0); |
397 |
|
|
for (ii=0; ii<num; ii++) { |
398 |
|
|
double tmpeval[3], tmprtp[3]; |
399 |
|
|
tenEigensolve_d(tmpeval, NULL, tenIn + 7*ii); |
400 |
|
|
tenTripleConvertSingle_d(tmprtp, tenTripleTypeRThetaPhi, |
401 |
|
|
tmpeval, tenTripleTypeEigenvalue); |
402 |
|
|
ww = wght ? wght[ii] : 1.0/num; |
403 |
|
|
TEN_T_SCALE_INCR(tmpTen, ww, tenIn + 7*ii); |
404 |
|
|
ELL_3V_SCALE_INCR(rtp, ww, tmprtp); |
405 |
|
|
} |
406 |
|
|
tenEigensolve_d(eval, evec, tmpTen); /* only care about evec */ |
407 |
|
|
tenTripleConvertSingle_d(eval, tenTripleTypeEigenvalue, |
408 |
|
|
rtp, tenTripleTypeRThetaPhi); |
409 |
|
|
tenMakeSingle_d(tenOut, tmpTen[0], eval, evec); |
410 |
|
|
break; |
411 |
|
|
default: |
412 |
|
|
biffAddf(TEN, "%s: sorry, interp type %s (%d) not implemented", |
413 |
|
|
me, airEnumStr(tenInterpType, ptype), ptype); |
414 |
|
|
TEN_T_COPY(tenOut, tenErr); |
415 |
|
|
return 1; |
416 |
|
|
} |
417 |
|
|
|
418 |
|
|
return 0; |
419 |
|
|
} |
420 |
|
|
|
421 |
|
|
int |
422 |
|
|
_tenInterpGeoLoxRelaxOne(Nrrd *nodata, Nrrd *ntdata, Nrrd *nigrtdata, |
423 |
|
|
unsigned int ii, int rotnoop, double scl, |
424 |
|
|
tenInterpParm *tip) { |
425 |
|
|
static const char me[]="_tenInterpGeoLoxRelaxOne"; |
426 |
|
|
double *tdata, *odata, *igrtdata, *tt[5], *igrt[5][6], d02[7], d24[7], |
427 |
|
|
len02, len24, tmp, tng[7], correct, update[7]; |
428 |
|
|
unsigned int jj; |
429 |
|
|
|
430 |
|
|
if (tip->verbose) { |
431 |
|
|
fprintf(stderr, "---- %u --> %u %u %u %u %u\n", ii, |
432 |
|
|
2*ii - 2, 2*ii - 1, 2*ii, 2*ii + 1, 2*ii + 2); |
433 |
|
|
} |
434 |
|
|
tdata = AIR_CAST(double *, ntdata->data); |
435 |
|
|
odata = AIR_CAST(double *, nodata->data); |
436 |
|
|
tt[0] = tdata + 7*(2*ii - 2); |
437 |
|
|
tt[1] = tdata + 7*(2*ii - 1); /* unused */ |
438 |
|
|
tt[2] = tdata + 7*(2*ii + 0); |
439 |
|
|
tt[3] = tdata + 7*(2*ii + 1); /* unused */ |
440 |
|
|
tt[4] = tdata + 7*(2*ii + 2); |
441 |
|
|
igrtdata = AIR_CAST(double *, nigrtdata->data); |
442 |
|
|
for (jj=0; jj<6; jj++) { |
443 |
|
|
igrt[0][jj] = igrtdata + 7*(jj + 6*(2*ii - 2)); /* unused */ |
444 |
|
|
igrt[1][jj] = igrtdata + 7*(jj + 6*(2*ii - 1)); |
445 |
|
|
igrt[2][jj] = igrtdata + 7*(jj + 6*(2*ii + 0)); |
446 |
|
|
igrt[3][jj] = igrtdata + 7*(jj + 6*(2*ii + 1)); |
447 |
|
|
igrt[4][jj] = igrtdata + 7*(jj + 6*(2*ii + 2)); /* unused */ |
448 |
|
|
} |
449 |
|
|
|
450 |
|
|
/* re-align [1] and [3] bases relative to [2] */ |
451 |
|
|
/* HEY: should I be worrying about aligning the mode normal |
452 |
|
|
when it had to be computed from eigenvectors? */ |
453 |
|
|
for (jj=3; jj<6; jj++) { |
454 |
|
|
if (TEN_T_DOT(igrt[1][jj], igrt[2][jj]) < 0) { |
455 |
|
|
TEN_T_SCALE(igrt[1][jj], -1, igrt[1][jj]); |
456 |
|
|
} |
457 |
|
|
if (TEN_T_DOT(igrt[3][jj], igrt[2][jj]) < 0) { |
458 |
|
|
TEN_T_SCALE(igrt[3][jj], -1, igrt[1][jj]); |
459 |
|
|
} |
460 |
|
|
} |
461 |
|
|
|
462 |
|
|
TEN_T_SUB(tng, tt[4], tt[0]); |
463 |
|
|
tmp = 1.0/TEN_T_NORM(tng); |
464 |
|
|
TEN_T_SCALE(tng, tmp, tng); |
465 |
|
|
|
466 |
|
|
TEN_T_SUB(d02, tt[2], tt[0]); |
467 |
|
|
TEN_T_SUB(d24, tt[4], tt[2]); |
468 |
|
|
TEN_T_SET(update, 1, 0, 0, 0, 0, 0, 0); |
469 |
|
|
for (jj=0; jj<(rotnoop ? 3u : 6u); jj++) { |
470 |
|
|
len02 = TEN_T_DOT(igrt[1][jj], d02); |
471 |
|
|
len24 = TEN_T_DOT(igrt[3][jj], d24); |
472 |
|
|
correct = (len24 - len02)/2; |
473 |
|
|
TEN_T_SCALE_INCR(update, correct*scl, igrt[2][jj]); |
474 |
|
|
if (tip->verbose) { |
475 |
|
|
fprintf(stderr, "igrt[1][%u] = %g %g %g %g %g %g\n", jj, |
476 |
|
|
igrt[1][jj][1], igrt[1][jj][2], igrt[1][jj][3], |
477 |
|
|
igrt[1][jj][4], igrt[1][jj][5], igrt[1][jj][6]); |
478 |
|
|
fprintf(stderr, "igrt[3][%u] = %g %g %g %g %g %g\n", jj, |
479 |
|
|
igrt[3][jj][1], igrt[3][jj][2], igrt[3][jj][3], |
480 |
|
|
igrt[3][jj][4], igrt[3][jj][5], igrt[3][jj][6]); |
481 |
|
|
fprintf(stderr, "(jj=%u) len = %g %g --> (d = %g) " |
482 |
|
|
"update = %g %g %g %g %g %g\n", |
483 |
|
|
jj, len02, len24, |
484 |
|
|
TEN_T_DOT(igrt[2][0], update), |
485 |
|
|
update[1], update[2], update[3], |
486 |
|
|
update[4], update[5], update[6]); |
487 |
|
|
} |
488 |
|
|
} |
489 |
|
|
if (rotnoop) { |
490 |
|
|
double avg[7], diff[7], len; |
491 |
|
|
TEN_T_LERP(avg, 0.5, tt[0], tt[4]); |
492 |
|
|
TEN_T_SUB(diff, avg, tt[2]); |
493 |
|
|
for (jj=0; jj<3; jj++) { |
494 |
|
|
len = TEN_T_DOT(igrt[2][jj], diff); |
495 |
|
|
TEN_T_SCALE_INCR(diff, -len, igrt[2][jj]); |
496 |
|
|
} |
497 |
|
|
TEN_T_SCALE_INCR(update, scl*0.2, diff); /* HEY: scaling is a hack */ |
498 |
|
|
if (tip->verbose) { |
499 |
|
|
fprintf(stderr, "(rotnoop) (d = %g) " |
500 |
|
|
"update = %g %g %g %g %g %g\n", |
501 |
|
|
TEN_T_DOT(igrt[2][0], update), |
502 |
|
|
update[1], update[2], update[3], |
503 |
|
|
update[4], update[5], update[6]); |
504 |
|
|
} |
505 |
|
|
} |
506 |
|
|
/* |
507 |
|
|
TEN_T_SUB(d02, tt[2], tt[0]); |
508 |
|
|
TEN_T_SUB(d24, tt[4], tt[2]); |
509 |
|
|
len02 = TEN_T_DOT(tng, d02); |
510 |
|
|
len24 = TEN_T_DOT(tng, d24); |
511 |
|
|
correct = (len24 - len02); |
512 |
|
|
TEN_T_SCALE_INCR(update, scl*correct, tng); |
513 |
|
|
*/ |
514 |
|
|
|
515 |
|
|
if (!TEN_T_EXISTS(update)) { |
516 |
|
|
biffAddf(TEN, "%s: computed non-existent update (step-size too big?)", me); |
517 |
|
|
return 1; |
518 |
|
|
} |
519 |
|
|
|
520 |
|
|
TEN_T_ADD(odata + 7*(2*ii + 0), tt[2], update); |
521 |
|
|
|
522 |
|
|
return 0; |
523 |
|
|
} |
524 |
|
|
|
525 |
|
|
void |
526 |
|
|
_tenInterpGeoLoxIGRT(double *igrt, double *ten, int useK, int rotNoop, |
527 |
|
|
double minnorm) { |
528 |
|
|
/* static const char me[]="_tenInterpGeoLoxIGRT"; */ |
529 |
|
|
double eval[3], evec[9]; |
530 |
|
|
|
531 |
|
|
if (useK) { |
532 |
|
|
tenInvariantGradientsK_d(igrt + 7*0, igrt + 7*1, igrt + 7*2, ten, minnorm); |
533 |
|
|
} else { |
534 |
|
|
tenInvariantGradientsR_d(igrt + 7*0, igrt + 7*1, igrt + 7*2, ten, minnorm); |
535 |
|
|
} |
536 |
|
|
if (rotNoop) { |
537 |
|
|
/* these shouldn't be used */ |
538 |
|
|
TEN_T_SET(igrt + 7*3, 1, AIR_NAN, AIR_NAN, AIR_NAN, |
539 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
540 |
|
|
TEN_T_SET(igrt + 7*4, 1, AIR_NAN, AIR_NAN, AIR_NAN, |
541 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
542 |
|
|
TEN_T_SET(igrt + 7*5, 1, AIR_NAN, AIR_NAN, AIR_NAN, |
543 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
544 |
|
|
} else { |
545 |
|
|
tenEigensolve_d(eval, evec, ten); |
546 |
|
|
tenRotationTangents_d(igrt + 7*3, igrt + 7*4, igrt + 7*5, evec); |
547 |
|
|
} |
548 |
|
|
return; |
549 |
|
|
} |
550 |
|
|
|
551 |
|
|
/* |
552 |
|
|
** if "doubling" is non-zero, this assumes that the real |
553 |
|
|
** vertices are on the even-numbered indices: |
554 |
|
|
** (0 1 2 3 4) |
555 |
|
|
** 0 2 4 6 8 --> size=9 --> NN=4 |
556 |
|
|
** 1 3 5 7 |
557 |
|
|
*/ |
558 |
|
|
double |
559 |
|
|
tenInterpPathLength(Nrrd *ntt, int doubleVerts, int fancy, int shape) { |
560 |
|
|
double *tt, len, diff[7], *tenA, *tenB; |
561 |
|
|
unsigned int ii, NN; |
562 |
|
|
|
563 |
|
|
tt = AIR_CAST(double *, ntt->data); |
564 |
|
|
if (doubleVerts) { |
565 |
|
|
NN = AIR_CAST(unsigned int, (ntt->axis[1].size-1)/2); |
566 |
|
|
} else { |
567 |
|
|
NN = AIR_CAST(unsigned int, ntt->axis[1].size-1); |
568 |
|
|
} |
569 |
|
|
len = 0; |
570 |
|
|
for (ii=0; ii<NN; ii++) { |
571 |
|
|
if (doubleVerts) { |
572 |
|
|
tenA = tt + 7*2*(ii + 1); |
573 |
|
|
tenB = tt + 7*2*(ii + 0); |
574 |
|
|
} else { |
575 |
|
|
tenA = tt + 7*(ii + 1); |
576 |
|
|
tenB = tt + 7*(ii + 0); |
577 |
|
|
} |
578 |
|
|
TEN_T_SUB(diff, tenA, tenB); |
579 |
|
|
if (fancy) { |
580 |
|
|
double mean[7], igrt[7*6], dot, incr; |
581 |
|
|
unsigned int lo, hi; |
582 |
|
|
|
583 |
|
|
TEN_T_LERP(mean, 0.5, tenA, tenB); |
584 |
|
|
_tenInterpGeoLoxIGRT(igrt, mean, AIR_FALSE, AIR_FALSE, 0.0); |
585 |
|
|
if (shape) { |
586 |
|
|
lo = 0; |
587 |
|
|
hi = 2; |
588 |
|
|
} else { |
589 |
|
|
lo = 3; |
590 |
|
|
hi = 5; |
591 |
|
|
} |
592 |
|
|
incr = 0; |
593 |
|
|
for (ii=lo; ii<=hi; ii++) { |
594 |
|
|
dot = TEN_T_DOT(igrt + 7*ii, diff); |
595 |
|
|
incr += dot*dot; |
596 |
|
|
} |
597 |
|
|
len += sqrt(incr); |
598 |
|
|
} else { |
599 |
|
|
len += TEN_T_NORM(diff); |
600 |
|
|
} |
601 |
|
|
} |
602 |
|
|
return len; |
603 |
|
|
} |
604 |
|
|
|
605 |
|
|
double |
606 |
|
|
_tenPathSpacingEqualize(Nrrd *nout, Nrrd *nin) { |
607 |
|
|
/* static const char me[]="_tenPathSpacingEqualize"; */ |
608 |
|
|
double *in, *out, len, diff[7], |
609 |
|
|
lenTotal, /* total length of input */ |
610 |
|
|
lenStep, /* correct length on input polyline between output vertices */ |
611 |
|
|
lenIn, /* length along input processed so far */ |
612 |
|
|
lenHere, /* length of segment associated with current input index */ |
613 |
|
|
lenRmn, /* length along past input segments as yet unmapped to output */ |
614 |
|
|
*tenHere, *tenNext; |
615 |
|
|
unsigned int idxIn, idxOut, NN; |
616 |
|
|
|
617 |
|
|
in = AIR_CAST(double *, nin->data); |
618 |
|
|
out = AIR_CAST(double *, nout->data); |
619 |
|
|
NN = (nin->axis[1].size-1)/2; |
620 |
|
|
lenTotal = tenInterpPathLength(nin, AIR_TRUE, AIR_FALSE, AIR_FALSE); |
621 |
|
|
lenStep = lenTotal/NN; |
622 |
|
|
/* |
623 |
|
|
fprintf(stderr, "!%s: lenTotal/NN = %g/%u = %g = lenStep\n", me, |
624 |
|
|
lenTotal, NN, lenStep); |
625 |
|
|
*/ |
626 |
|
|
TEN_T_COPY(out + 7*2*(0 + 0), in + 7*2*(0 + 0)); |
627 |
|
|
lenIn = lenRmn = 0; |
628 |
|
|
idxOut = 1; |
629 |
|
|
for (idxIn=0; idxIn<NN; idxIn++) { |
630 |
|
|
tenNext = in + 7*2*(idxIn + 1); |
631 |
|
|
tenHere = in + 7*2*(idxIn + 0); |
632 |
|
|
TEN_T_SUB(diff, tenNext, tenHere); |
633 |
|
|
lenHere = TEN_T_NORM(diff); |
634 |
|
|
/* |
635 |
|
|
fprintf(stderr, "!%s(%u): %g + %g >(%s)= %g\n", me, idxIn, |
636 |
|
|
lenRmn, lenHere, |
637 |
|
|
(lenRmn + lenHere >= lenStep ? "yes" : "no"), |
638 |
|
|
lenStep); |
639 |
|
|
*/ |
640 |
|
|
if (lenRmn + lenHere >= lenStep) { |
641 |
|
|
len = lenRmn + lenHere; |
642 |
|
|
while (len > lenStep) { |
643 |
|
|
len -= lenStep; |
644 |
|
|
/* |
645 |
|
|
fprintf(stderr, "!%s(%u): len = %g -> %g\n", me, idxIn, |
646 |
|
|
len + lenStep, len); |
647 |
|
|
*/ |
648 |
|
|
TEN_T_AFFINE(out + 7*(2*idxOut + 0), |
649 |
|
|
lenHere, len, 0, tenHere, tenNext); |
650 |
|
|
/* |
651 |
|
|
fprintf(stderr, "!%s(%u): out[%u] ~ %g\n", me, idxIn, idxOut, |
652 |
|
|
AIR_AFFINE(lenHere, len, 0, 0, 1)); |
653 |
|
|
*/ |
654 |
|
|
idxOut++; |
655 |
|
|
} |
656 |
|
|
lenRmn = len; |
657 |
|
|
} else { |
658 |
|
|
lenRmn += lenHere; |
659 |
|
|
/* |
660 |
|
|
fprintf(stderr, "!%s(%u): (==> lenRmn = %g -> %g)\n", me, idxIn, |
661 |
|
|
lenRmn - lenHere, lenRmn); |
662 |
|
|
*/ |
663 |
|
|
} |
664 |
|
|
/* now lenRmn < lenStep */ |
665 |
|
|
lenIn += lenHere; |
666 |
|
|
} |
667 |
|
|
/* copy very last one in case we didn't get to it somehow */ |
668 |
|
|
TEN_T_COPY(out + 7*2*(NN + 0), in + 7*2*(NN + 0)); |
669 |
|
|
|
670 |
|
|
/* fill in vertex mid-points */ |
671 |
|
|
for (idxOut=0; idxOut<NN; idxOut++) { |
672 |
|
|
TEN_T_LERP(out + 7*(2*idxOut + 1), |
673 |
|
|
0.5, out + 7*(2*idxOut + 0), out + 7*(2*idxOut + 2)); |
674 |
|
|
} |
675 |
|
|
return lenTotal; |
676 |
|
|
} |
677 |
|
|
|
678 |
|
|
int |
679 |
|
|
_tenInterpGeoLoxPolyLine(Nrrd *ngeod, unsigned int *numIter, |
680 |
|
|
const double tenA[7], const double tenB[7], |
681 |
|
|
unsigned int NN, int useK, int rotnoop, |
682 |
|
|
tenInterpParm *tip) { |
683 |
|
|
static const char me[]="_tenInterpGeoLoxPolyLine"; |
684 |
|
|
Nrrd *nigrt, *ntt, *nss, *nsub; |
685 |
|
|
double *igrt, *geod, *tt, len, newlen; |
686 |
|
|
unsigned int ii; |
687 |
|
|
airArray *mop; |
688 |
|
|
|
689 |
|
|
if (!(ngeod && numIter && tenA && tenB)) { |
690 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
691 |
|
|
return 1; |
692 |
|
|
} |
693 |
|
|
if (!(NN >= 2)) { |
694 |
|
|
biffAddf(TEN, "%s: # steps %u too small", me, NN); |
695 |
|
|
return 1; |
696 |
|
|
} |
697 |
|
|
|
698 |
|
|
mop = airMopNew(); |
699 |
|
|
ntt = nrrdNew(); |
700 |
|
|
airMopAdd(mop, ntt, (airMopper)nrrdNuke, airMopAlways); |
701 |
|
|
nss = nrrdNew(); |
702 |
|
|
airMopAdd(mop, nss, (airMopper)nrrdNuke, airMopAlways); |
703 |
|
|
nigrt = nrrdNew(); |
704 |
|
|
airMopAdd(mop, nigrt, (airMopper)nrrdNuke, airMopAlways); |
705 |
|
|
nsub = nrrdNew(); |
706 |
|
|
airMopAdd(mop, nsub, (airMopper)nrrdNuke, airMopAlways); |
707 |
|
|
if (nrrdMaybeAlloc_va(ngeod, nrrdTypeDouble, 2, |
708 |
|
|
AIR_CAST(size_t, 7), |
709 |
|
|
AIR_CAST(size_t, NN+1)) |
710 |
|
|
|| nrrdMaybeAlloc_va(ntt, nrrdTypeDouble, 2, |
711 |
|
|
AIR_CAST(size_t, 7), |
712 |
|
|
AIR_CAST(size_t, 2*NN + 1)) |
713 |
|
|
|| nrrdMaybeAlloc_va(nigrt, nrrdTypeDouble, 3, |
714 |
|
|
AIR_CAST(size_t, 7), |
715 |
|
|
AIR_CAST(size_t, 6), |
716 |
|
|
AIR_CAST(size_t, 2*NN + 1))) { |
717 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate output", me); |
718 |
|
|
airMopError(mop); return 1; |
719 |
|
|
} |
720 |
|
|
geod = AIR_CAST(double *, ngeod->data); |
721 |
|
|
tt = AIR_CAST(double *, ntt->data); |
722 |
|
|
igrt = AIR_CAST(double *, nigrt->data); |
723 |
|
|
|
724 |
|
|
*numIter = 0; |
725 |
|
|
if (NN > 14 && tip->enableRecurse) { |
726 |
|
|
unsigned int subIter; |
727 |
|
|
int E; |
728 |
|
|
NrrdResampleContext *rsmc; |
729 |
|
|
double kparm[3] = {1.0, 0.0, 0.5}; |
730 |
|
|
/* recurse and find geodesic with smaller number of vertices */ |
731 |
|
|
if (_tenInterpGeoLoxPolyLine(nsub, &subIter, tenA, tenB, |
732 |
|
|
NN/2, useK, rotnoop, tip)) { |
733 |
|
|
biffAddf(TEN, "%s: problem with recursive call", me); |
734 |
|
|
airMopError(mop); return 1; |
735 |
|
|
} |
736 |
|
|
/* upsample coarse geodesic to higher resolution */ |
737 |
|
|
rsmc = nrrdResampleContextNew(); |
738 |
|
|
airMopAdd(mop, rsmc, (airMopper)nrrdResampleContextNix, airMopAlways); |
739 |
|
|
E = AIR_FALSE; |
740 |
|
|
if (!E) E |= nrrdResampleDefaultCenterSet(rsmc, nrrdCenterNode); |
741 |
|
|
if (!E) E |= nrrdResampleInputSet(rsmc, nsub); |
742 |
|
|
if (!E) E |= nrrdResampleKernelSet(rsmc, 0, NULL, NULL); |
743 |
|
|
if (!E) E |= nrrdResampleKernelSet(rsmc, 1, nrrdKernelTent, kparm); |
744 |
|
|
if (!E) E |= nrrdResampleSamplesSet(rsmc, 1, 2*NN + 1); |
745 |
|
|
if (!E) E |= nrrdResampleRangeFullSet(rsmc, 1); |
746 |
|
|
if (!E) E |= nrrdResampleBoundarySet(rsmc, nrrdBoundaryBleed); |
747 |
|
|
if (!E) E |= nrrdResampleTypeOutSet(rsmc, nrrdTypeDefault); |
748 |
|
|
if (!E) E |= nrrdResampleRenormalizeSet(rsmc, AIR_TRUE); |
749 |
|
|
if (!E) E |= nrrdResampleExecute(rsmc, ntt); |
750 |
|
|
if (E) { |
751 |
|
|
biffMovef(TEN, NRRD, "%s: problem upsampling course solution", me); |
752 |
|
|
airMopError(mop); return 1; |
753 |
|
|
} |
754 |
|
|
*numIter += subIter; |
755 |
|
|
} else { |
756 |
|
|
/* initialize the path, including all the segment midpoints */ |
757 |
|
|
for (ii=0; ii<=2*NN; ii++) { |
758 |
|
|
TEN_T_AFFINE(tt + 7*ii, 0, ii, 2*NN, tenA, tenB); |
759 |
|
|
} |
760 |
|
|
} |
761 |
|
|
for (ii=0; ii<=2*NN; ii++) { |
762 |
|
|
_tenInterpGeoLoxIGRT(igrt + 7*6*ii, tt + 7*ii, useK, rotnoop, |
763 |
|
|
tip->minNorm); |
764 |
|
|
} |
765 |
|
|
nrrdCopy(nss, ntt); |
766 |
|
|
|
767 |
|
|
newlen = tenInterpPathLength(ntt, AIR_TRUE, AIR_FALSE, AIR_FALSE); |
768 |
|
|
do { |
769 |
|
|
unsigned int lo, hi; |
770 |
|
|
int dd; |
771 |
|
|
len = newlen; |
772 |
|
|
if (0 == *numIter % 2) { |
773 |
|
|
lo = 1; |
774 |
|
|
hi = NN; |
775 |
|
|
dd = 1; |
776 |
|
|
} else { |
777 |
|
|
lo = NN-1; |
778 |
|
|
hi = 0; |
779 |
|
|
dd = -1; |
780 |
|
|
} |
781 |
|
|
if (tip->verbose) { |
782 |
|
|
fprintf(stderr, "%s: ======= iter = %u (NN=%u)\n", me, *numIter, NN); |
783 |
|
|
} |
784 |
|
|
for (ii=lo; ii!=hi; ii+=dd) { |
785 |
|
|
double sclHack; |
786 |
|
|
sclHack = ii*4.0/NN - ii*ii*4.0/NN/NN; |
787 |
|
|
if (_tenInterpGeoLoxRelaxOne(nss, ntt, nigrt, ii, rotnoop, |
788 |
|
|
sclHack*tip->convStep, tip)) { |
789 |
|
|
biffAddf(TEN, "%s: problem on vert %u, iter %u\n", me, ii, *numIter); |
790 |
|
|
return 1; |
791 |
|
|
} |
792 |
|
|
} |
793 |
|
|
newlen = _tenPathSpacingEqualize(ntt, nss); |
794 |
|
|
/* try doing this less often */ |
795 |
|
|
for (ii=0; ii<=2*NN; ii++) { |
796 |
|
|
_tenInterpGeoLoxIGRT(igrt + 7*6*ii, tt + 7*ii, useK, rotnoop, |
797 |
|
|
tip->minNorm); |
798 |
|
|
} |
799 |
|
|
*numIter += 1; |
800 |
|
|
} while ((0 == tip->maxIter || *numIter < tip->maxIter) |
801 |
|
|
&& 2*AIR_ABS(newlen - len)/(newlen + len) > tip->convEps); |
802 |
|
|
|
803 |
|
|
/* copy final result to output */ |
804 |
|
|
for (ii=0; ii<=NN; ii++) { |
805 |
|
|
TEN_T_COPY(geod + 7*ii, tt + 7*2*ii); |
806 |
|
|
} |
807 |
|
|
/* values from outer-most recursion will stick */ |
808 |
|
|
tip->numIter = *numIter; |
809 |
|
|
tip->convFinal = 2*AIR_ABS(newlen - len)/(newlen + len); |
810 |
|
|
|
811 |
|
|
airMopOkay(mop); |
812 |
|
|
return 0; |
813 |
|
|
} |
814 |
|
|
|
815 |
|
|
int |
816 |
|
|
tenInterpTwoDiscrete_d(Nrrd *nout, |
817 |
|
|
const double tenA[7], const double tenB[7], |
818 |
|
|
int ptype, unsigned int num, |
819 |
|
|
tenInterpParm *_tip) { |
820 |
|
|
static const char me[]="tenInterpTwoDiscrete_d"; |
821 |
|
|
double *out; |
822 |
|
|
unsigned int ii; |
823 |
|
|
airArray *mop; |
824 |
|
|
tenInterpParm *tip; |
825 |
|
|
|
826 |
|
|
if (!nout) { |
827 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
828 |
|
|
return 1; |
829 |
|
|
} |
830 |
|
|
if (airEnumValCheck(tenInterpType, ptype)) { |
831 |
|
|
biffAddf(TEN, "%s: path type %d not a valid %s", me, ptype, |
832 |
|
|
tenInterpType->name); |
833 |
|
|
return 1; |
834 |
|
|
} |
835 |
|
|
|
836 |
|
|
mop = airMopNew(); |
837 |
|
|
if (_tip) { |
838 |
|
|
tip = _tip; |
839 |
|
|
} else { |
840 |
|
|
tip = tenInterpParmNew(); |
841 |
|
|
airMopAdd(mop, tip, (airMopper)tenInterpParmNix, airMopAlways); |
842 |
|
|
} |
843 |
|
|
if (!( num >= 2 )) { |
844 |
|
|
biffAddf(TEN, "%s: need num >= 2 (not %u)", me, num); |
845 |
|
|
airMopError(mop); return 1; |
846 |
|
|
} |
847 |
|
|
if (nrrdMaybeAlloc_va(nout, nrrdTypeDouble, 2, |
848 |
|
|
AIR_CAST(size_t, 7), |
849 |
|
|
AIR_CAST(size_t, num))) { |
850 |
|
|
biffMovef(TEN, NRRD, "%s: trouble allocating output", me); |
851 |
|
|
airMopError(mop); return 1; |
852 |
|
|
} |
853 |
|
|
out = AIR_CAST(double *, nout->data); |
854 |
|
|
|
855 |
|
|
if (ptype == tenInterpTypeLinear |
856 |
|
|
|| ptype == tenInterpTypeLogLinear |
857 |
|
|
|| ptype == tenInterpTypeAffineInvariant |
858 |
|
|
|| ptype == tenInterpTypeWang |
859 |
|
|
|| ptype == tenInterpTypeQuatGeoLoxK |
860 |
|
|
|| ptype == tenInterpTypeQuatGeoLoxR |
861 |
|
|
|| ptype == tenInterpTypeRThetaPhiLinear) { |
862 |
|
|
/* we have fast ways of doing interpolation |
863 |
|
|
between two tensors for these path types */ |
864 |
|
|
for (ii=0; ii<num; ii++) { |
865 |
|
|
/* yes, this is often doing a lot of needless recomputations. */ |
866 |
|
|
tenInterpTwo_d(out + 7*ii, tenA, tenB, |
867 |
|
|
ptype, (double)ii/(num-1), tip); |
868 |
|
|
} |
869 |
|
|
} else if (ptype == tenInterpTypeGeoLoxK |
870 |
|
|
|| ptype == tenInterpTypeGeoLoxR |
871 |
|
|
|| ptype == tenInterpTypeLoxK |
872 |
|
|
|| ptype == tenInterpTypeLoxR) { |
873 |
|
|
/* we have slow iterative code for these */ |
874 |
|
|
unsigned int numIter; |
875 |
|
|
int useK, rotnoop; |
876 |
|
|
|
877 |
|
|
useK = (tenInterpTypeGeoLoxK == ptype |
878 |
|
|
|| tenInterpTypeLoxK == ptype); |
879 |
|
|
rotnoop = (tenInterpTypeGeoLoxK == ptype |
880 |
|
|
|| tenInterpTypeGeoLoxR == ptype); |
881 |
|
|
fprintf(stderr, "!%s: useK = %d, rotnoop = %d\n", me, useK, rotnoop); |
882 |
|
|
if (_tenInterpGeoLoxPolyLine(nout, &numIter, |
883 |
|
|
tenA, tenB, |
884 |
|
|
num, useK, rotnoop, tip)) { |
885 |
|
|
biffAddf(TEN, "%s: trouble finding path", me); |
886 |
|
|
airMopError(mop); return 1; |
887 |
|
|
} |
888 |
|
|
} else { |
889 |
|
|
biffAddf(TEN, "%s: sorry, interp for path %s not implemented", me, |
890 |
|
|
airEnumStr(tenInterpType, ptype)); |
891 |
|
|
airMopError(mop); return 1; |
892 |
|
|
} |
893 |
|
|
|
894 |
|
|
airMopOkay(mop); |
895 |
|
|
return 0; |
896 |
|
|
} |
897 |
|
|
|
898 |
|
|
double |
899 |
|
|
tenInterpDistanceTwo_d(const double tenA[7], const double tenB[7], |
900 |
|
|
int ptype, tenInterpParm *_tip) { |
901 |
|
|
static const char me[]="tenInterpDistanceTwo_d"; |
902 |
|
|
char *err; |
903 |
|
|
tenInterpParm *tip; |
904 |
|
|
airArray *mop; |
905 |
|
|
double ret, diff[7], logA[7], logB[7], invA[7], det, siA[7], |
906 |
|
|
mat1[9], mat2[9], mat3[9], logDiff[7]; |
907 |
|
|
Nrrd *npath; |
908 |
|
|
|
909 |
|
|
if (!( tenA && tenB && !airEnumValCheck(tenInterpType, ptype) )) { |
910 |
|
|
return AIR_NAN; |
911 |
|
|
} |
912 |
|
|
|
913 |
|
|
mop = airMopNew(); |
914 |
|
|
switch (ptype) { |
915 |
|
|
case tenInterpTypeLinear: |
916 |
|
|
TEN_T_SUB(diff, tenA, tenB); |
917 |
|
|
ret = TEN_T_NORM(diff); |
918 |
|
|
break; |
919 |
|
|
case tenInterpTypeLogLinear: |
920 |
|
|
tenLogSingle_d(logA, tenA); |
921 |
|
|
tenLogSingle_d(logB, tenB); |
922 |
|
|
TEN_T_SUB(diff, logA, logB); |
923 |
|
|
ret = TEN_T_NORM(diff); |
924 |
|
|
break; |
925 |
|
|
case tenInterpTypeAffineInvariant: |
926 |
|
|
TEN_T_INV(invA, tenA, det); |
927 |
|
|
tenSqrtSingle_d(siA, invA); |
928 |
|
|
TEN_T2M(mat1, tenB); |
929 |
|
|
TEN_T2M(mat2, siA); |
930 |
|
|
ell_3m_mul_d(mat3, mat1, mat2); |
931 |
|
|
ell_3m_mul_d(mat1, mat2, mat3); |
932 |
|
|
TEN_M2T(diff, mat1); |
933 |
|
|
tenLogSingle_d(logDiff, diff); |
934 |
|
|
ret = TEN_T_NORM(logDiff); |
935 |
|
|
break; |
936 |
|
|
case tenInterpTypeGeoLoxK: |
937 |
|
|
case tenInterpTypeGeoLoxR: |
938 |
|
|
case tenInterpTypeLoxK: |
939 |
|
|
case tenInterpTypeLoxR: |
940 |
|
|
case tenInterpTypeQuatGeoLoxK: |
941 |
|
|
case tenInterpTypeQuatGeoLoxR: |
942 |
|
|
npath = nrrdNew(); |
943 |
|
|
airMopAdd(mop, npath, (airMopper)nrrdNuke, airMopAlways); |
944 |
|
|
if (_tip) { |
945 |
|
|
tip = _tip; |
946 |
|
|
} else { |
947 |
|
|
tip = tenInterpParmNew(); |
948 |
|
|
airMopAdd(mop, tip, (airMopper)tenInterpParmNix, airMopAlways); |
949 |
|
|
} |
950 |
|
|
if (tenInterpTwoDiscrete_d(npath, tenA, tenB, ptype, |
951 |
|
|
tip->numSteps, tip)) { |
952 |
|
|
airMopAdd(mop, err = biffGetDone(TEN), airFree, airMopAlways); |
953 |
|
|
fprintf(stderr, "%s: trouble computing path:\n%s\n", me, err); |
954 |
|
|
airMopError(mop); return AIR_NAN; |
955 |
|
|
} |
956 |
|
|
ret = tenInterpPathLength(npath, AIR_FALSE, AIR_FALSE, AIR_FALSE); |
957 |
|
|
if (tip->lengthFancy) { |
958 |
|
|
tip->lengthShape = tenInterpPathLength(npath, AIR_FALSE, |
959 |
|
|
AIR_TRUE, AIR_TRUE); |
960 |
|
|
tip->lengthOrient = tenInterpPathLength(npath, AIR_FALSE, |
961 |
|
|
AIR_TRUE, AIR_FALSE); |
962 |
|
|
} |
963 |
|
|
break; |
964 |
|
|
case tenInterpTypeWang: |
965 |
|
|
default: |
966 |
|
|
fprintf(stderr, "%s: unimplemented %s %d!!!!\n", me, |
967 |
|
|
tenInterpType->name, ptype); |
968 |
|
|
ret = AIR_NAN; |
969 |
|
|
break; |
970 |
|
|
} |
971 |
|
|
|
972 |
|
|
airMopOkay(mop); |
973 |
|
|
return ret; |
974 |
|
|
} |
975 |
|
|
|
976 |
|
|
/* |
977 |
|
|
** actually, the input nrrds don't have to be 3D ... |
978 |
|
|
*/ |
979 |
|
|
int |
980 |
|
|
tenInterpMulti3D(Nrrd *nout, const Nrrd *const *nin, const double *wght, |
981 |
|
|
unsigned int ninLen, int ptype, tenInterpParm *_tip) { |
982 |
|
|
static const char me[]="tenInterpMulti3D"; |
983 |
|
|
unsigned int ninIdx; |
984 |
|
|
size_t II, NN; |
985 |
|
|
double (*lup)(const void *, size_t), (*ins)(void *, size_t, double), |
986 |
|
|
*tbuff; |
987 |
|
|
tenInterpParm *tip; |
988 |
|
|
airArray *mop; |
989 |
|
|
|
990 |
|
|
/* allow NULL wght, to signify equal weighting */ |
991 |
|
|
if (!(nout && nin)) { |
992 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
993 |
|
|
return 1; |
994 |
|
|
} |
995 |
|
|
if (!( ninLen > 0 )) { |
996 |
|
|
biffAddf(TEN, "%s: need at least 1 nin, not 0", me); |
997 |
|
|
return 1; |
998 |
|
|
} |
999 |
|
|
if (airEnumValCheck(tenInterpType, ptype)) { |
1000 |
|
|
biffAddf(TEN, "%s: invalid %s %d", me, |
1001 |
|
|
tenInterpType->name, ptype); |
1002 |
|
|
return 1; |
1003 |
|
|
} |
1004 |
|
|
|
1005 |
|
|
if (tenTensorCheck(nin[0], nrrdTypeDefault, AIR_FALSE, AIR_TRUE)) { |
1006 |
|
|
biffAddf(TEN, "%s: first nrrd not a tensor array", me); |
1007 |
|
|
return 1; |
1008 |
|
|
} |
1009 |
|
|
if (!( nrrdTypeFloat == nin[0]->type || |
1010 |
|
|
nrrdTypeDouble == nin[0]->type )) { |
1011 |
|
|
biffAddf(TEN, "%s: need type %s or %s (not %s) in first nrrd", me, |
1012 |
|
|
airEnumStr(nrrdType, nrrdTypeFloat), |
1013 |
|
|
airEnumStr(nrrdType, nrrdTypeDouble), |
1014 |
|
|
airEnumStr(nrrdType, nin[0]->type)); |
1015 |
|
|
return 1; |
1016 |
|
|
} |
1017 |
|
|
for (ninIdx=1; ninIdx<ninLen; ninIdx++) { |
1018 |
|
|
if (tenTensorCheck(nin[ninIdx], nrrdTypeDefault, AIR_FALSE, AIR_TRUE)) { |
1019 |
|
|
biffAddf(TEN, "%s: nin[%u] not a tensor array", me, ninIdx); |
1020 |
|
|
return 1; |
1021 |
|
|
} |
1022 |
|
|
if (!nrrdSameSize(nin[0], nin[ninIdx], AIR_TRUE)) { |
1023 |
|
|
biffMovef(TEN, NRRD, "%s: nin[0] doesn't match nin[%u]", me, ninIdx); |
1024 |
|
|
return 1; |
1025 |
|
|
} |
1026 |
|
|
if (nin[0]->type != nin[ninIdx]->type) { |
1027 |
|
|
biffAddf(TEN, "%s: nin[0] type (%s) != nin[%u] type (%s)", me, |
1028 |
|
|
airEnumStr(nrrdType, nin[0]->type), |
1029 |
|
|
ninIdx, airEnumStr(nrrdType, nin[ninIdx]->type)); |
1030 |
|
|
return 1; |
1031 |
|
|
} |
1032 |
|
|
} |
1033 |
|
|
|
1034 |
|
|
mop = airMopNew(); |
1035 |
|
|
if (nrrdCopy(nout, nin[0])) { |
1036 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't initialize output", me); |
1037 |
|
|
airMopError(mop); return 1; |
1038 |
|
|
} |
1039 |
|
|
if (_tip) { |
1040 |
|
|
tip = _tip; |
1041 |
|
|
} else { |
1042 |
|
|
tip = tenInterpParmNew(); |
1043 |
|
|
airMopAdd(mop, tip, (airMopper)tenInterpParmNix, airMopAlways); |
1044 |
|
|
} |
1045 |
|
|
tbuff = AIR_CAST(double *, calloc(7*ninLen, sizeof(double))); |
1046 |
|
|
if (!tbuff) { |
1047 |
|
|
biffAddf(TEN, "%s: couldn't allocate tensor buff", me); |
1048 |
|
|
airMopError(mop); return 1; |
1049 |
|
|
} |
1050 |
|
|
ins = nrrdDInsert[nin[0]->type]; |
1051 |
|
|
lup = nrrdDLookup[nin[0]->type]; |
1052 |
|
|
NN = nrrdElementNumber(nin[0])/7; |
1053 |
|
|
for (II=0; II<NN; II++) { |
1054 |
|
|
double tenOut[7]; |
1055 |
|
|
unsigned int tt; |
1056 |
|
|
for (ninIdx=0; ninIdx<ninLen; ninIdx++) { |
1057 |
|
|
for (tt=0; tt<7; tt++) { |
1058 |
|
|
tbuff[tt + 7*ninIdx] = lup(nin[ninIdx]->data, tt + 7*II); |
1059 |
|
|
} |
1060 |
|
|
} |
1061 |
|
|
if (tenInterpN_d(tenOut, tbuff, wght, ninLen, ptype, tip)) { |
1062 |
|
|
char stmp[AIR_STRLEN_SMALL]; |
1063 |
|
|
biffAddf(TEN, "%s: trouble on sample %s", me, |
1064 |
|
|
airSprintSize_t(stmp, II)); |
1065 |
|
|
airMopError(mop); return 1; |
1066 |
|
|
} |
1067 |
|
|
for (tt=0; tt<7; tt++) { |
1068 |
|
|
ins(nout->data, tt + 7*II, tenOut[tt]); |
1069 |
|
|
} |
1070 |
|
|
} |
1071 |
|
|
|
1072 |
|
|
airMopOkay(mop); |
1073 |
|
|
return 0; |
1074 |
|
|
} |