Bug Summary

File:src/ten/path.c
Location:line 135, column 7
Description:Potential leak of memory pointed to by 'newtip'

Annotated Source Code

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
28tenInterpParm *
29tenInterpParmNew(void) {
30 tenInterpParm *tip;
31
32 tip = AIR_CAST(tenInterpParm *, malloc(sizeof(tenInterpParm)))((tenInterpParm *)(malloc(sizeof(tenInterpParm))));
2
←
Within the expansion of the macro 'AIR_CAST':
→
a
Memory is allocated
33 if (tip) {
3
←
Assuming 'tip' is non-null
→
4
←
Taking true branch
→
34 tip->verbose = AIR_FALSE0;
35 tip->convStep = 0.2;
36 tip->minNorm = 0.0;
37 tip->convEps = 0.0000000001;
38 tip->wghtSumEps = 0.0000001;
39 tip->enableRecurse = AIR_TRUE1;
40 tip->maxIter = 20;
41 tip->numSteps = 100;
42 tip->lengthFancy = AIR_FALSE0;
43
44 tip->allocLen = 0;
45 tip->eval = NULL((void*)0);
46 tip->evec = NULL((void*)0);
47 tip->rtIn = NULL((void*)0);
48 tip->rtLog = NULL((void*)0);
49 tip->qIn = NULL((void*)0);
50 tip->qBuff = NULL((void*)0);
51 tip->qInter = NULL((void*)0);
52
53 tip->numIter = 0;
54 tip->convFinal = AIR_NAN(airFloatQNaN.f);
55 tip->lengthShape = AIR_NAN(airFloatQNaN.f);
56 tip->lengthOrient = AIR_NAN(airFloatQNaN.f);
57 }
58 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*/
68int
69tenInterpParmBufferAlloc(tenInterpParm *tip, unsigned int num) {
70 static const char me[]="tenInterpParmBufferAlloc";
71
72 if (0 == num) {
73 /* user wants to free buffers for some reason */
74 airFree(tip->eval); tip->eval = NULL((void*)0);
75 airFree(tip->evec); tip->evec = NULL((void*)0);
76 airFree(tip->rtIn); tip->rtIn = NULL((void*)0);
77 airFree(tip->rtLog); tip->rtLog = NULL((void*)0);
78 airFree(tip->qIn); tip->qIn = NULL((void*)0);
79 airFree(tip->qBuff); tip->qBuff = NULL((void*)0);
80 airFree(tip->qInter); tip->qInter = NULL((void*)0);
81 tip->allocLen = 0;
82 } else if (1 == num) {
83 biffAddf(TENtenBiffKey, "%s: need num >= 2 (not %u)", me, num);
84 return 1;
85 } else if (num != tip->allocLen) {
86 airFree(tip->eval); tip->eval = NULL((void*)0);
87 airFree(tip->evec); tip->evec = NULL((void*)0);
88 airFree(tip->rtIn); tip->rtIn = NULL((void*)0);
89 airFree(tip->rtLog); tip->rtLog = NULL((void*)0);
90 airFree(tip->qIn); tip->qIn = NULL((void*)0);
91 airFree(tip->qBuff); tip->qBuff = NULL((void*)0);
92 airFree(tip->qInter); tip->qInter = NULL((void*)0);
93 tip->eval = AIR_CALLOC(3*num, double)(double*)(calloc((3*num), sizeof(double)));
94 tip->evec = AIR_CALLOC(9*num, double)(double*)(calloc((9*num), sizeof(double)));
95 tip->rtIn = AIR_CALLOC(3*num, double)(double*)(calloc((3*num), sizeof(double)));
96 tip->rtLog = AIR_CALLOC(3*num, double)(double*)(calloc((3*num), sizeof(double)));
97 tip->qIn = AIR_CALLOC(4*num, double)(double*)(calloc((4*num), sizeof(double)));
98 tip->qBuff = AIR_CALLOC(4*num, double)(double*)(calloc((4*num), sizeof(double)));
99 tip->qInter = AIR_CALLOC(num*num, double)(double*)(calloc((num*num), sizeof(double)));
100 if (!(tip->eval && tip->evec &&
101 tip->rtIn && tip->rtLog &&
102 tip->qIn && tip->qBuff && tip->qInter)) {
103 biffAddf(TENtenBiffKey, "%s: didn't alloc buffers (%p,%p,%p %p %p %p %p)", me,
104 AIR_VOIDP(tip->eval)((void *)(tip->eval)), AIR_VOIDP(tip->evec)((void *)(tip->evec)),
105 AIR_VOIDP(tip->rtIn)((void *)(tip->rtIn)), AIR_VOIDP(tip->rtLog)((void *)(tip->rtLog)),
106 AIR_VOIDP(tip->qIn)((void *)(tip->qIn)), AIR_VOIDP(tip->qBuff)((void *)(tip->qBuff)),
107 AIR_VOIDP(tip->qInter)((void *)(tip->qInter)));
108 return 1;
109 }
110 tip->allocLen = num;
111 }
112 return 0;
113}
114
115tenInterpParm *
116tenInterpParmCopy(tenInterpParm *tip) {
117 static const char me[]="tenInterpParmCopy";
118 tenInterpParm *newtip;
119 unsigned int num;
120
121 num = tip->allocLen;
122 newtip = tenInterpParmNew();
1
Calling 'tenInterpParmNew'
→
5
←
Returned allocated memory
→
123 if (newtip) {
6
←
Taking true branch
→
124 memcpy(newtip, tip, sizeof(tenInterpParm))__builtin___memcpy_chk (newtip, tip, sizeof(tenInterpParm), __builtin_object_size
(newtip, 0));
125 /* manually set all pointers */
126 newtip->allocLen = 0;
127 newtip->eval = NULL((void*)0);
128 newtip->evec = NULL((void*)0);
129 newtip->rtIn = NULL((void*)0);
130 newtip->rtLog = NULL((void*)0);
131 newtip->qIn = NULL((void*)0);
132 newtip->qBuff = NULL((void*)0);
133 newtip->qInter = NULL((void*)0);
134 if (tenInterpParmBufferAlloc(newtip, num)) {
7
←
Taking true branch
→
135 biffAddf(TENtenBiffKey, "%s: trouble allocating output", me);
8
←
Potential leak of memory pointed to by 'newtip'
136 return NULL((void*)0);
137 }
138 memcpy(newtip->eval, tip->eval, 3*num*sizeof(double))__builtin___memcpy_chk (newtip->eval, tip->eval, 3*num*
sizeof(double), __builtin_object_size (newtip->eval, 0));
139 memcpy(newtip->evec, tip->evec, 9*num*sizeof(double))__builtin___memcpy_chk (newtip->evec, tip->evec, 9*num*
sizeof(double), __builtin_object_size (newtip->evec, 0));
140 memcpy(newtip->rtIn, tip->rtIn, 3*num*sizeof(double))__builtin___memcpy_chk (newtip->rtIn, tip->rtIn, 3*num*
sizeof(double), __builtin_object_size (newtip->rtIn, 0));
141 memcpy(newtip->rtLog, tip->rtLog, 3*num*sizeof(double))__builtin___memcpy_chk (newtip->rtLog, tip->rtLog, 3*num
*sizeof(double), __builtin_object_size (newtip->rtLog, 0));
142 memcpy(newtip->qIn, tip->qIn, 4*num*sizeof(double))__builtin___memcpy_chk (newtip->qIn, tip->qIn, 4*num*sizeof
(double), __builtin_object_size (newtip->qIn, 0));
143 memcpy(newtip->qBuff, tip->qBuff, 4*num*sizeof(double))__builtin___memcpy_chk (newtip->qBuff, tip->qBuff, 4*num
*sizeof(double), __builtin_object_size (newtip->qBuff, 0));
144 memcpy(newtip->qInter, tip->qInter, num*num*sizeof(double))__builtin___memcpy_chk (newtip->qInter, tip->qInter, num
*num*sizeof(double), __builtin_object_size (newtip->qInter
, 0));
145 }
146 return newtip;
147}
148
149tenInterpParm *
150tenInterpParmNix(tenInterpParm *tip) {
151
152 if (tip) {
153 airFree(tip->eval);
154 airFree(tip->evec);
155 airFree(tip->rtIn);
156 airFree(tip->rtLog);
157 airFree(tip->qIn);
158 airFree(tip->qBuff);
159 airFree(tip->qInter);
160 free(tip);
161 }
162 return NULL((void*)0);
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*/
178void
179tenInterpTwo_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,( (oten)[0] = ((airFloatQNaN.f)), (oten)[1] = ((airFloatQNaN.
f)), (oten)[2] = ((airFloatQNaN.f)), (oten)[3] = ((airFloatQNaN
.f)), (oten)[4] = ((airFloatQNaN.f)), (oten)[5] = ((airFloatQNaN
.f)), (oten)[6] = ((airFloatQNaN.f)) )
194 AIR_NAN, AIR_NAN, AIR_NAN)( (oten)[0] = ((airFloatQNaN.f)), (oten)[1] = ((airFloatQNaN.
f)), (oten)[2] = ((airFloatQNaN.f)), (oten)[3] = ((airFloatQNaN
.f)), (oten)[4] = ((airFloatQNaN.f)), (oten)[5] = ((airFloatQNaN
.f)), (oten)[6] = ((airFloatQNaN.f)) );
195 }
196 return;
197 }
198
199 switch (ptype) {
200 case tenInterpTypeLinear:
201 TEN_T_LERP(oten, aa, tenA, tenB)( (oten)[0] = (((aa))*(((tenB)[0]) - ((tenA)[0])) + ((tenA)[0
])), (oten)[1] = (((aa))*(((tenB)[1]) - ((tenA)[1])) + ((tenA
)[1])), (oten)[2] = (((aa))*(((tenB)[2]) - ((tenA)[2])) + ((tenA
)[2])), (oten)[3] = (((aa))*(((tenB)[3]) - ((tenA)[3])) + ((tenA
)[3])), (oten)[4] = (((aa))*(((tenB)[4]) - ((tenA)[4])) + ((tenA
)[4])), (oten)[5] = (((aa))*(((tenB)[5]) - ((tenA)[5])) + ((tenA
)[5])), (oten)[6] = (((aa))*(((tenB)[6]) - ((tenA)[6])) + ((tenA
)[6])));
202 break;
203 case tenInterpTypeLogLinear:
204 tenLogSingle_d(logA, tenA);
205 tenLogSingle_d(logB, tenB);
206 TEN_T_LERP(logMean, aa, logA, logB)( (logMean)[0] = (((aa))*(((logB)[0]) - ((logA)[0])) + ((logA
)[0])), (logMean)[1] = (((aa))*(((logB)[1]) - ((logA)[1])) + (
(logA)[1])), (logMean)[2] = (((aa))*(((logB)[2]) - ((logA)[2]
)) + ((logA)[2])), (logMean)[3] = (((aa))*(((logB)[3]) - ((logA
)[3])) + ((logA)[3])), (logMean)[4] = (((aa))*(((logB)[4]) - (
(logA)[4])) + ((logA)[4])), (logMean)[5] = (((aa))*(((logB)[5
]) - ((logA)[5])) + ((logA)[5])), (logMean)[6] = (((aa))*(((logB
)[6]) - ((logA)[6])) + ((logA)[6])));
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)( (mat1)[0] = (tenB)[1], (mat1)[1] = (tenB)[2], (mat1)[2] = (
tenB)[3], (mat1)[3] = (tenB)[2], (mat1)[4] = (tenB)[4], (mat1
)[5] = (tenB)[5], (mat1)[6] = (tenB)[3], (mat1)[7] = (tenB)[5
], (mat1)[8] = (tenB)[6] );
213 TEN_T2M(mat2, isqrtA)( (mat2)[0] = (isqrtA)[1], (mat2)[1] = (isqrtA)[2], (mat2)[2]
= (isqrtA)[3], (mat2)[3] = (isqrtA)[2], (mat2)[4] = (isqrtA)
[4], (mat2)[5] = (isqrtA)[5], (mat2)[6] = (isqrtA)[3], (mat2)
[7] = (isqrtA)[5], (mat2)[8] = (isqrtA)[6] );
214 ELL_3M_MUL(mat3, mat1, mat2)((mat3)[0] = (mat1)[0]*(mat2)[0] + (mat1)[1]*(mat2)[3] + (mat1
)[2]*(mat2)[6], (mat3)[1] = (mat1)[0]*(mat2)[1] + (mat1)[1]*(
mat2)[4] + (mat1)[2]*(mat2)[7], (mat3)[2] = (mat1)[0]*(mat2)[
2] + (mat1)[1]*(mat2)[5] + (mat1)[2]*(mat2)[8], (mat3)[3] = (
mat1)[3]*(mat2)[0] + (mat1)[4]*(mat2)[3] + (mat1)[5]*(mat2)[6
], (mat3)[4] = (mat1)[3]*(mat2)[1] + (mat1)[4]*(mat2)[4] + (mat1
)[5]*(mat2)[7], (mat3)[5] = (mat1)[3]*(mat2)[2] + (mat1)[4]*(
mat2)[5] + (mat1)[5]*(mat2)[8], (mat3)[6] = (mat1)[6]*(mat2)[
0] + (mat1)[7]*(mat2)[3] + (mat1)[8]*(mat2)[6], (mat3)[7] = (
mat1)[6]*(mat2)[1] + (mat1)[7]*(mat2)[4] + (mat1)[8]*(mat2)[7
], (mat3)[8] = (mat1)[6]*(mat2)[2] + (mat1)[7]*(mat2)[5] + (mat1
)[8]*(mat2)[8]); /* B * is(A) */
215 ELL_3M_MUL(mat1, mat2, mat3)((mat1)[0] = (mat2)[0]*(mat3)[0] + (mat2)[1]*(mat3)[3] + (mat2
)[2]*(mat3)[6], (mat1)[1] = (mat2)[0]*(mat3)[1] + (mat2)[1]*(
mat3)[4] + (mat2)[2]*(mat3)[7], (mat1)[2] = (mat2)[0]*(mat3)[
2] + (mat2)[1]*(mat3)[5] + (mat2)[2]*(mat3)[8], (mat1)[3] = (
mat2)[3]*(mat3)[0] + (mat2)[4]*(mat3)[3] + (mat2)[5]*(mat3)[6
], (mat1)[4] = (mat2)[3]*(mat3)[1] + (mat2)[4]*(mat3)[4] + (mat2
)[5]*(mat3)[7], (mat1)[5] = (mat2)[3]*(mat3)[2] + (mat2)[4]*(
mat3)[5] + (mat2)[5]*(mat3)[8], (mat1)[6] = (mat2)[6]*(mat3)[
0] + (mat2)[7]*(mat3)[3] + (mat2)[8]*(mat3)[6], (mat1)[7] = (
mat2)[6]*(mat3)[1] + (mat2)[7]*(mat3)[4] + (mat2)[8]*(mat3)[7
], (mat1)[8] = (mat2)[6]*(mat3)[2] + (mat2)[7]*(mat3)[5] + (mat2
)[8]*(mat3)[8]); /* is(A) * B * is(A) */
216 TEN_M2T(tmp2, mat1)( (tmp2)[1] = (mat1)[0], (tmp2)[2] = ((mat1)[1]+(mat1)[3])/2.0
, (tmp2)[3] = ((mat1)[2]+(mat1)[6])/2.0, (tmp2)[4] = (mat1)[4
], (tmp2)[5] = ((mat1)[5]+(mat1)[7])/2.0, (tmp2)[6] = (mat1)[
8] );
217 tenPowSingle_d(tmp1, tmp2, aa); /* m = (is(A) * B * is(A))^aa */
218 TEN_T2M(mat1, tmp1)( (mat1)[0] = (tmp1)[1], (mat1)[1] = (tmp1)[2], (mat1)[2] = (
tmp1)[3], (mat1)[3] = (tmp1)[2], (mat1)[4] = (tmp1)[4], (mat1
)[5] = (tmp1)[5], (mat1)[6] = (tmp1)[3], (mat1)[7] = (tmp1)[5
], (mat1)[8] = (tmp1)[6] );
219 TEN_T2M(mat2, sqrtA)( (mat2)[0] = (sqrtA)[1], (mat2)[1] = (sqrtA)[2], (mat2)[2] =
(sqrtA)[3], (mat2)[3] = (sqrtA)[2], (mat2)[4] = (sqrtA)[4], (
mat2)[5] = (sqrtA)[5], (mat2)[6] = (sqrtA)[3], (mat2)[7] = (sqrtA
)[5], (mat2)[8] = (sqrtA)[6] );
220 ELL_3M_MUL(mat3, mat1, mat2)((mat3)[0] = (mat1)[0]*(mat2)[0] + (mat1)[1]*(mat2)[3] + (mat1
)[2]*(mat2)[6], (mat3)[1] = (mat1)[0]*(mat2)[1] + (mat1)[1]*(
mat2)[4] + (mat1)[2]*(mat2)[7], (mat3)[2] = (mat1)[0]*(mat2)[
2] + (mat1)[1]*(mat2)[5] + (mat1)[2]*(mat2)[8], (mat3)[3] = (
mat1)[3]*(mat2)[0] + (mat1)[4]*(mat2)[3] + (mat1)[5]*(mat2)[6
], (mat3)[4] = (mat1)[3]*(mat2)[1] + (mat1)[4]*(mat2)[4] + (mat1
)[5]*(mat2)[7], (mat3)[5] = (mat1)[3]*(mat2)[2] + (mat1)[4]*(
mat2)[5] + (mat1)[5]*(mat2)[8], (mat3)[6] = (mat1)[6]*(mat2)[
0] + (mat1)[7]*(mat2)[3] + (mat1)[8]*(mat2)[6], (mat3)[7] = (
mat1)[6]*(mat2)[1] + (mat1)[7]*(mat2)[4] + (mat1)[8]*(mat2)[7
], (mat3)[8] = (mat1)[6]*(mat2)[2] + (mat1)[7]*(mat2)[5] + (mat1
)[8]*(mat2)[8]); /* m * sqrt(A) */
221 ELL_3M_MUL(mat1, mat2, mat3)((mat1)[0] = (mat2)[0]*(mat3)[0] + (mat2)[1]*(mat3)[3] + (mat2
)[2]*(mat3)[6], (mat1)[1] = (mat2)[0]*(mat3)[1] + (mat2)[1]*(
mat3)[4] + (mat2)[2]*(mat3)[7], (mat1)[2] = (mat2)[0]*(mat3)[
2] + (mat2)[1]*(mat3)[5] + (mat2)[2]*(mat3)[8], (mat1)[3] = (
mat2)[3]*(mat3)[0] + (mat2)[4]*(mat3)[3] + (mat2)[5]*(mat3)[6
], (mat1)[4] = (mat2)[3]*(mat3)[1] + (mat2)[4]*(mat3)[4] + (mat2
)[5]*(mat3)[7], (mat1)[5] = (mat2)[3]*(mat3)[2] + (mat2)[4]*(
mat3)[5] + (mat2)[5]*(mat3)[8], (mat1)[6] = (mat2)[6]*(mat3)[
0] + (mat2)[7]*(mat3)[3] + (mat2)[8]*(mat3)[6], (mat1)[7] = (
mat2)[6]*(mat3)[1] + (mat2)[7]*(mat3)[4] + (mat2)[8]*(mat3)[7
], (mat1)[8] = (mat2)[6]*(mat3)[2] + (mat2)[7]*(mat3)[5] + (mat2
)[8]*(mat3)[8]); /* sqrt(A) * m * sqrt(A) */
222 TEN_M2T(oten, mat1)( (oten)[1] = (mat1)[0], (oten)[2] = ((mat1)[1]+(mat1)[3])/2.0
, (oten)[3] = ((mat1)[2]+(mat1)[6])/2.0, (oten)[4] = (mat1)[4
], (oten)[5] = ((mat1)[5]+(mat1)[7])/2.0, (oten)[6] = (mat1)[
8] );
223 oten[0] = AIR_LERP(aa, tenA[0], tenB[0])((aa)*((tenB[0]) - (tenA[0])) + (tenA[0]));
224 if (tip->verbose) {
225 fprintf(stderr__stderrp, "%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)( (mean)[0] = (((aa))*(((tenB)[0]) - ((tenA)[0])) + ((tenA)[0
])), (mean)[1] = (((aa))*(((tenB)[1]) - ((tenA)[1])) + ((tenA
)[1])), (mean)[2] = (((aa))*(((tenB)[2]) - ((tenA)[2])) + ((tenA
)[2])), (mean)[3] = (((aa))*(((tenB)[3]) - ((tenA)[3])) + ((tenA
)[3])), (mean)[4] = (((aa))*(((tenB)[4]) - ((tenA)[4])) + ((tenA
)[4])), (mean)[5] = (((aa))*(((tenB)[5]) - ((tenA)[5])) + ((tenA
)[5])), (mean)[6] = (((aa))*(((tenB)[6]) - ((tenA)[6])) + ((tenA
)[6]))); /* "A" = mean */
238 tenLogSingle_d(logA, tenA);
239 tenLogSingle_d(logB, tenB);
240 TEN_T_LERP(logMean, aa, logA, logB)( (logMean)[0] = (((aa))*(((logB)[0]) - ((logA)[0])) + ((logA
)[0])), (logMean)[1] = (((aa))*(((logB)[1]) - ((logA)[1])) + (
(logA)[1])), (logMean)[2] = (((aa))*(((logB)[2]) - ((logA)[2]
)) + ((logA)[2])), (logMean)[3] = (((aa))*(((logB)[3]) - ((logA
)[3])) + ((logA)[3])), (logMean)[4] = (((aa))*(((logB)[4]) - (
(logA)[4])) + ((logA)[4])), (logMean)[5] = (((aa))*(((logB)[5
]) - ((logA)[5])) + ((logA)[5])), (logMean)[6] = (((aa))*(((logB
)[6]) - ((logA)[6])) + ((logA)[6]))); /* "B" = logMean */
241 tenSqrtSingle_d(sqrtB, logMean);
242 tenInv_d(isqrtB, sqrtB);
243 TEN_T2M(mat1, mean)( (mat1)[0] = (mean)[1], (mat1)[1] = (mean)[2], (mat1)[2] = (
mean)[3], (mat1)[3] = (mean)[2], (mat1)[4] = (mean)[4], (mat1
)[5] = (mean)[5], (mat1)[6] = (mean)[3], (mat1)[7] = (mean)[5
], (mat1)[8] = (mean)[6] );
244 TEN_T2M(mat2, isqrtB)( (mat2)[0] = (isqrtB)[1], (mat2)[1] = (isqrtB)[2], (mat2)[2]
= (isqrtB)[3], (mat2)[3] = (isqrtB)[2], (mat2)[4] = (isqrtB)
[4], (mat2)[5] = (isqrtB)[5], (mat2)[6] = (isqrtB)[3], (mat2)
[7] = (isqrtB)[5], (mat2)[8] = (isqrtB)[6] );
245 ELL_3M_MUL(mat3, mat1, mat2)((mat3)[0] = (mat1)[0]*(mat2)[0] + (mat1)[1]*(mat2)[3] + (mat1
)[2]*(mat2)[6], (mat3)[1] = (mat1)[0]*(mat2)[1] + (mat1)[1]*(
mat2)[4] + (mat1)[2]*(mat2)[7], (mat3)[2] = (mat1)[0]*(mat2)[
2] + (mat1)[1]*(mat2)[5] + (mat1)[2]*(mat2)[8], (mat3)[3] = (
mat1)[3]*(mat2)[0] + (mat1)[4]*(mat2)[3] + (mat1)[5]*(mat2)[6
], (mat3)[4] = (mat1)[3]*(mat2)[1] + (mat1)[4]*(mat2)[4] + (mat1
)[5]*(mat2)[7], (mat3)[5] = (mat1)[3]*(mat2)[2] + (mat1)[4]*(
mat2)[5] + (mat1)[5]*(mat2)[8], (mat3)[6] = (mat1)[6]*(mat2)[
0] + (mat1)[7]*(mat2)[3] + (mat1)[8]*(mat2)[6], (mat3)[7] = (
mat1)[6]*(mat2)[1] + (mat1)[7]*(mat2)[4] + (mat1)[8]*(mat2)[7
], (mat3)[8] = (mat1)[6]*(mat2)[2] + (mat1)[7]*(mat2)[5] + (mat1
)[8]*(mat2)[8]);
246 ELL_3M_MUL(mat1, mat2, mat3)((mat1)[0] = (mat2)[0]*(mat3)[0] + (mat2)[1]*(mat3)[3] + (mat2
)[2]*(mat3)[6], (mat1)[1] = (mat2)[0]*(mat3)[1] + (mat2)[1]*(
mat3)[4] + (mat2)[2]*(mat3)[7], (mat1)[2] = (mat2)[0]*(mat3)[
2] + (mat2)[1]*(mat3)[5] + (mat2)[2]*(mat3)[8], (mat1)[3] = (
mat2)[3]*(mat3)[0] + (mat2)[4]*(mat3)[3] + (mat2)[5]*(mat3)[6
], (mat1)[4] = (mat2)[3]*(mat3)[1] + (mat2)[4]*(mat3)[4] + (mat2
)[5]*(mat3)[7], (mat1)[5] = (mat2)[3]*(mat3)[2] + (mat2)[4]*(
mat3)[5] + (mat2)[5]*(mat3)[8], (mat1)[6] = (mat2)[6]*(mat3)[
0] + (mat2)[7]*(mat3)[3] + (mat2)[8]*(mat3)[6], (mat1)[7] = (
mat2)[6]*(mat3)[1] + (mat2)[7]*(mat3)[4] + (mat2)[8]*(mat3)[7
], (mat1)[8] = (mat2)[6]*(mat3)[2] + (mat2)[7]*(mat3)[5] + (mat2
)[8]*(mat3)[8]);
247 TEN_M2T(tmp1, mat1)( (tmp1)[1] = (mat1)[0], (tmp1)[2] = ((mat1)[1]+(mat1)[3])/2.0
, (tmp1)[3] = ((mat1)[2]+(mat1)[6])/2.0, (tmp1)[4] = (mat1)[4
], (tmp1)[5] = ((mat1)[5]+(mat1)[7])/2.0, (tmp1)[6] = (mat1)[
8] );
248 tenSqrtSingle_d(oten, tmp1);
249 oten[0] = AIR_LERP(aa, tenA[0], tenB[0])((aa)*((tenB[0]) - (tenA[0])) + (tenA[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])((aa)*((tenB[0]) - (tenA[0])) + (tenA[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,( (oten)[0] = ((airFloatQNaN.f)), (oten)[1] = ((airFloatQNaN.
f)), (oten)[2] = ((airFloatQNaN.f)), (oten)[3] = ((airFloatQNaN
.f)), (oten)[4] = ((airFloatQNaN.f)), (oten)[5] = ((airFloatQNaN
.f)), (oten)[6] = ((airFloatQNaN.f)) )
269 AIR_NAN, AIR_NAN, AIR_NAN)( (oten)[0] = ((airFloatQNaN.f)), (oten)[1] = ((airFloatQNaN.
f)), (oten)[2] = ((airFloatQNaN.f)), (oten)[3] = ((airFloatQNaN
.f)), (oten)[4] = ((airFloatQNaN.f)), (oten)[5] = ((airFloatQNaN
.f)), (oten)[6] = ((airFloatQNaN.f)) );
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((void*)0), tenA);
275 tenTripleConvertSingle_d(rtpA, tenTripleTypeRThetaPhi,
276 eval, tenTripleTypeEigenvalue);
277 tenEigensolve_d(eval, NULL((void*)0), tenB);
278 tenTripleConvertSingle_d(rtpB, tenTripleTypeRThetaPhi,
279 eval, tenTripleTypeEigenvalue);
280 TEN_T_LERP(tenM, aa, tenA, tenB)( (tenM)[0] = (((aa))*(((tenB)[0]) - ((tenA)[0])) + ((tenA)[0
])), (tenM)[1] = (((aa))*(((tenB)[1]) - ((tenA)[1])) + ((tenA
)[1])), (tenM)[2] = (((aa))*(((tenB)[2]) - ((tenA)[2])) + ((tenA
)[2])), (tenM)[3] = (((aa))*(((tenB)[3]) - ((tenA)[3])) + ((tenA
)[3])), (tenM)[4] = (((aa))*(((tenB)[4]) - ((tenA)[4])) + ((tenA
)[4])), (tenM)[5] = (((aa))*(((tenB)[5]) - ((tenA)[5])) + ((tenA
)[5])), (tenM)[6] = (((aa))*(((tenB)[6]) - ((tenA)[6])) + ((tenA
)[6])));
281 tenEigensolve_d(eval, oevec, tenM);
282 ELL_3V_LERP(rtpM, aa, rtpA, rtpB)((rtpM)[0] = (((aa))*(((rtpB)[0]) - ((rtpA)[0])) + ((rtpA)[0]
)), (rtpM)[1] = (((aa))*(((rtpB)[1]) - ((rtpA)[1])) + ((rtpA)
[1])), (rtpM)[2] = (((aa))*(((rtpB)[2]) - ((rtpA)[2])) + ((rtpA
)[2])));
283 tenTripleConvertSingle_d(oeval, tenTripleTypeEigenvalue,
284 rtpM, tenTripleTypeRThetaPhi);
285 }
286 tenMakeSingle_d(oten, AIR_LERP(aa, tenA[0], tenB[0])((aa)*((tenB[0]) - (tenA[0])) + (tenA[0])), oeval, oevec);
287 break;
288 default:
289 /* error */
290 TEN_T_SET(oten, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN,( (oten)[0] = ((airFloatQNaN.f)), (oten)[1] = ((airFloatQNaN.
f)), (oten)[2] = ((airFloatQNaN.f)), (oten)[3] = ((airFloatQNaN
.f)), (oten)[4] = ((airFloatQNaN.f)), (oten)[5] = ((airFloatQNaN
.f)), (oten)[6] = ((airFloatQNaN.f)) )
291 AIR_NAN, AIR_NAN, AIR_NAN)( (oten)[0] = ((airFloatQNaN.f)), (oten)[1] = ((airFloatQNaN.
f)), (oten)[2] = ((airFloatQNaN.f)), (oten)[3] = ((airFloatQNaN
.f)), (oten)[4] = ((airFloatQNaN.f)), (oten)[5] = ((airFloatQNaN
.f)), (oten)[6] = ((airFloatQNaN.f)) );
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*/
303int
304tenInterpN_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,( (tenErr)[0] = ((airFloatQNaN.f)), (tenErr)[1] = ((airFloatQNaN
.f)), (tenErr)[2] = ((airFloatQNaN.f)), (tenErr)[3] = ((airFloatQNaN
.f)), (tenErr)[4] = ((airFloatQNaN.f)), (tenErr)[5] = ((airFloatQNaN
.f)), (tenErr)[6] = ((airFloatQNaN.f)) )
312 AIR_NAN, AIR_NAN, AIR_NAN)( (tenErr)[0] = ((airFloatQNaN.f)), (tenErr)[1] = ((airFloatQNaN
.f)), (tenErr)[2] = ((airFloatQNaN.f)), (tenErr)[3] = ((airFloatQNaN
.f)), (tenErr)[4] = ((airFloatQNaN.f)), (tenErr)[5] = ((airFloatQNaN
.f)), (tenErr)[6] = ((airFloatQNaN.f)) );
313 /* wght can be NULL ==> equal 1/num weight for all */
314 if (!(tenOut && tenIn && tip)) {
315 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
316 return 1;
317 }
318 if (!( num >= 2 )) {
319 biffAddf(TENtenBiffKey, "%s: need num >= 2 (not %u)", me, num);
320 TEN_T_COPY(tenOut, tenErr)( (tenOut)[0] = (tenErr)[0], (tenOut)[1] = (tenErr)[1], (tenOut
)[2] = (tenErr)[2], (tenOut)[3] = (tenErr)[3], (tenOut)[4] = (
tenErr)[4], (tenOut)[5] = (tenErr)[5], (tenOut)[6] = (tenErr)
[6] ); return 1;
321 }
322 if (airEnumValCheck(tenInterpType, ptype)) {
323 biffAddf(TENtenBiffKey, "%s: invalid %s %d", me, tenInterpType->name, ptype);
324 TEN_T_COPY(tenOut, tenErr)( (tenOut)[0] = (tenErr)[0], (tenOut)[1] = (tenErr)[1], (tenOut
)[2] = (tenErr)[2], (tenOut)[3] = (tenErr)[3], (tenOut)[4] = (
tenErr)[4], (tenOut)[5] = (tenErr)[5], (tenOut)[6] = (tenErr)
[6] ); 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)((1 - tip->wghtSumEps) <= (wghtSum) && (wghtSum
) <= (1 + tip->wghtSumEps)) )) {
332 biffAddf(TENtenBiffKey, "%s: wght sum %g not within %g of 1.0", me,
333 wghtSum, tip->wghtSumEps);
334 TEN_T_COPY(tenOut, tenErr)( (tenOut)[0] = (tenErr)[0], (tenOut)[1] = (tenErr)[1], (tenOut
)[2] = (tenErr)[2], (tenOut)[3] = (tenErr)[3], (tenOut)[4] = (
tenErr)[4], (tenOut)[5] = (tenErr)[5], (tenOut)[6] = (tenErr)
[6] ); return 1;
335 }
336
337 switch (ptype) {
338 case tenInterpTypeLinear:
339 TEN_T_SET(tenOut, 0, 0, 0, 0, 0, 0, 0)( (tenOut)[0] = (0), (tenOut)[1] = (0), (tenOut)[2] = (0), (tenOut
)[3] = (0), (tenOut)[4] = (0), (tenOut)[5] = (0), (tenOut)[6]
= (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)( (tenOut)[0] = (tenIn + 7*ii)[0], (tenOut)[1] += (ww)*(tenIn
+ 7*ii)[1], (tenOut)[2] += (ww)*(tenIn + 7*ii)[2], (tenOut)[
3] += (ww)*(tenIn + 7*ii)[3], (tenOut)[4] += (ww)*(tenIn + 7*
ii)[4], (tenOut)[5] += (ww)*(tenIn + 7*ii)[5], (tenOut)[6] +=
(ww)*(tenIn + 7*ii)[6]);
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)( (tenOut)[0] = (0), (tenOut)[1] = (0), (tenOut)[2] = (0), (tenOut
)[3] = (0), (tenOut)[4] = (0), (tenOut)[5] = (0), (tenOut)[6]
= (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)( (tenOut)[0] = (tmpTen)[0], (tenOut)[1] += (ww)*(tmpTen)[1],
(tenOut)[2] += (ww)*(tmpTen)[2], (tenOut)[3] += (ww)*(tmpTen
)[3], (tenOut)[4] += (ww)*(tmpTen)[4], (tenOut)[5] += (ww)*(tmpTen
)[5], (tenOut)[6] += (ww)*(tmpTen)[6]);
355 cc += ww*(tenIn + 7*ii)[0];
356 }
357 tenOut[0] = cc;
358 TEN_T_COPY(tmpTen, tenOut)( (tmpTen)[0] = (tenOut)[0], (tmpTen)[1] = (tenOut)[1], (tmpTen
)[2] = (tenOut)[2], (tmpTen)[3] = (tenOut)[3], (tmpTen)[4] = (
tenOut)[4], (tmpTen)[5] = (tenOut)[5], (tmpTen)[6] = (tenOut)
[6] );
359 tenExpSingle_d(tenOut, tmpTen);
360 break;
361 case tenInterpTypeAffineInvariant:
362 case tenInterpTypeWang:
363 biffAddf(TENtenBiffKey, "%s: sorry, not implemented", me);
364 TEN_T_COPY(tenOut, tenErr)( (tenOut)[0] = (tenErr)[0], (tenOut)[1] = (tenErr)[1], (tenOut
)[2] = (tenErr)[2], (tenOut)[3] = (tenErr)[3], (tenOut)[4] = (
tenErr)[4], (tenOut)[5] = (tenErr)[5], (tenOut)[6] = (tenErr)
[6] ); return 1;
365 break;
366 case tenInterpTypeGeoLoxK:
367 case tenInterpTypeGeoLoxR:
368 case tenInterpTypeLoxK:
369 case tenInterpTypeLoxR:
370 biffAddf(TENtenBiffKey, "%s: %s doesn't support averaging multiple tensors", me,
371 airEnumStr(tenInterpType, ptype));
372 TEN_T_COPY(tenOut, tenErr)( (tenOut)[0] = (tenErr)[0], (tenOut)[1] = (tenErr)[1], (tenOut
)[2] = (tenErr)[2], (tenOut)[3] = (tenErr)[3], (tenOut)[4] = (
tenErr)[4], (tenOut)[5] = (tenErr)[5], (tenOut)[6] = (tenErr)
[6] ); return 1;
373 break;
374 case tenInterpTypeQuatGeoLoxK:
375 case tenInterpTypeQuatGeoLoxR:
376 if (tenInterpParmBufferAlloc(tip, num)) {
377 biffAddf(TENtenBiffKey, "%s: trouble getting buffers", me);
378 TEN_T_COPY(tenOut, tenErr)( (tenOut)[0] = (tenErr)[0], (tenOut)[1] = (tenErr)[1], (tenOut
)[2] = (tenErr)[2], (tenOut)[3] = (tenErr)[3], (tenOut)[4] = (
tenErr)[4], (tenOut)[5] = (tenErr)[5], (tenOut)[6] = (tenErr)
[6] ); 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(TENtenBiffKey, "%s: trouble computing", me);
389 TEN_T_COPY(tenOut, tenErr)( (tenOut)[0] = (tenErr)[0], (tenOut)[1] = (tenErr)[1], (tenOut
)[2] = (tenErr)[2], (tenOut)[3] = (tenErr)[3], (tenOut)[4] = (
tenErr)[4], (tenOut)[5] = (tenErr)[5], (tenOut)[6] = (tenErr)
[6] ); 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)( (tmpTen)[0] = (0), (tmpTen)[1] = (0), (tmpTen)[2] = (0), (tmpTen
)[3] = (0), (tmpTen)[4] = (0), (tmpTen)[5] = (0), (tmpTen)[6]
= (0) );
396 ELL_3V_SET(rtp, 0, 0, 0)((rtp)[0] = (0), (rtp)[1] = (0), (rtp)[2] = (0));
397 for (ii=0; ii<num; ii++) {
398 double tmpeval[3], tmprtp[3];
399 tenEigensolve_d(tmpeval, NULL((void*)0), 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)( (tmpTen)[0] = (tenIn + 7*ii)[0], (tmpTen)[1] += (ww)*(tenIn
+ 7*ii)[1], (tmpTen)[2] += (ww)*(tenIn + 7*ii)[2], (tmpTen)[
3] += (ww)*(tenIn + 7*ii)[3], (tmpTen)[4] += (ww)*(tenIn + 7*
ii)[4], (tmpTen)[5] += (ww)*(tenIn + 7*ii)[5], (tmpTen)[6] +=
(ww)*(tenIn + 7*ii)[6]);
404 ELL_3V_SCALE_INCR(rtp, ww, tmprtp)((rtp)[0] += (ww)*(tmprtp)[0], (rtp)[1] += (ww)*(tmprtp)[1], (
rtp)[2] += (ww)*(tmprtp)[2]);
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(TENtenBiffKey, "%s: sorry, interp type %s (%d) not implemented",
413 me, airEnumStr(tenInterpType, ptype), ptype);
414 TEN_T_COPY(tenOut, tenErr)( (tenOut)[0] = (tenErr)[0], (tenOut)[1] = (tenErr)[1], (tenOut
)[2] = (tenErr)[2], (tenOut)[3] = (tenErr)[3], (tenOut)[4] = (
tenErr)[4], (tenOut)[5] = (tenErr)[5], (tenOut)[6] = (tenErr)
[6] );
415 return 1;
416 }
417
418 return 0;
419}
420
421int
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__stderrp, "---- %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)((double *)(ntdata->data));
435 odata = AIR_CAST(double *, nodata->data)((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)((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])( (igrt[1][jj])[1]*(igrt[2][jj])[1] + 2*(igrt[1][jj])[2]*(igrt
[2][jj])[2] + 2*(igrt[1][jj])[3]*(igrt[2][jj])[3] + (igrt[1][
jj])[4]*(igrt[2][jj])[4] + 2*(igrt[1][jj])[5]*(igrt[2][jj])[5
] + (igrt[1][jj])[6]*(igrt[2][jj])[6] ) < 0) {
455 TEN_T_SCALE(igrt[1][jj], -1, igrt[1][jj])( (igrt[1][jj])[0] = (igrt[1][jj])[0], (igrt[1][jj])[1] = (-1
)*(igrt[1][jj])[1], (igrt[1][jj])[2] = (-1)*(igrt[1][jj])[2],
(igrt[1][jj])[3] = (-1)*(igrt[1][jj])[3], (igrt[1][jj])[4] =
(-1)*(igrt[1][jj])[4], (igrt[1][jj])[5] = (-1)*(igrt[1][jj])
[5], (igrt[1][jj])[6] = (-1)*(igrt[1][jj])[6]);
456 }
457 if (TEN_T_DOT(igrt[3][jj], igrt[2][jj])( (igrt[3][jj])[1]*(igrt[2][jj])[1] + 2*(igrt[3][jj])[2]*(igrt
[2][jj])[2] + 2*(igrt[3][jj])[3]*(igrt[2][jj])[3] + (igrt[3][
jj])[4]*(igrt[2][jj])[4] + 2*(igrt[3][jj])[5]*(igrt[2][jj])[5
] + (igrt[3][jj])[6]*(igrt[2][jj])[6] ) < 0) {
458 TEN_T_SCALE(igrt[3][jj], -1, igrt[1][jj])( (igrt[3][jj])[0] = (igrt[1][jj])[0], (igrt[3][jj])[1] = (-1
)*(igrt[1][jj])[1], (igrt[3][jj])[2] = (-1)*(igrt[1][jj])[2],
(igrt[3][jj])[3] = (-1)*(igrt[1][jj])[3], (igrt[3][jj])[4] =
(-1)*(igrt[1][jj])[4], (igrt[3][jj])[5] = (-1)*(igrt[1][jj])
[5], (igrt[3][jj])[6] = (-1)*(igrt[1][jj])[6]);
459 }
460 }
461
462 TEN_T_SUB(tng, tt[4], tt[0])( (tng)[0] = ((tt[4])[0] + (tt[0])[0])/2.0, (tng)[1] = (tt[4]
)[1] - (tt[0])[1], (tng)[2] = (tt[4])[2] - (tt[0])[2], (tng)[
3] = (tt[4])[3] - (tt[0])[3], (tng)[4] = (tt[4])[4] - (tt[0])
[4], (tng)[5] = (tt[4])[5] - (tt[0])[5], (tng)[6] = (tt[4])[6
] - (tt[0])[6]);
463 tmp = 1.0/TEN_T_NORM(tng)(sqrt(( (tng)[1]*(tng)[1] + 2*(tng)[2]*(tng)[2] + 2*(tng)[3]*
(tng)[3] + (tng)[4]*(tng)[4] + 2*(tng)[5]*(tng)[5] + (tng)[6]
*(tng)[6] )));
464 TEN_T_SCALE(tng, tmp, tng)( (tng)[0] = (tng)[0], (tng)[1] = (tmp)*(tng)[1], (tng)[2] = (
tmp)*(tng)[2], (tng)[3] = (tmp)*(tng)[3], (tng)[4] = (tmp)*(tng
)[4], (tng)[5] = (tmp)*(tng)[5], (tng)[6] = (tmp)*(tng)[6]);
465
466 TEN_T_SUB(d02, tt[2], tt[0])( (d02)[0] = ((tt[2])[0] + (tt[0])[0])/2.0, (d02)[1] = (tt[2]
)[1] - (tt[0])[1], (d02)[2] = (tt[2])[2] - (tt[0])[2], (d02)[
3] = (tt[2])[3] - (tt[0])[3], (d02)[4] = (tt[2])[4] - (tt[0])
[4], (d02)[5] = (tt[2])[5] - (tt[0])[5], (d02)[6] = (tt[2])[6
] - (tt[0])[6]);
467 TEN_T_SUB(d24, tt[4], tt[2])( (d24)[0] = ((tt[4])[0] + (tt[2])[0])/2.0, (d24)[1] = (tt[4]
)[1] - (tt[2])[1], (d24)[2] = (tt[4])[2] - (tt[2])[2], (d24)[
3] = (tt[4])[3] - (tt[2])[3], (d24)[4] = (tt[4])[4] - (tt[2])
[4], (d24)[5] = (tt[4])[5] - (tt[2])[5], (d24)[6] = (tt[4])[6
] - (tt[2])[6]);
468 TEN_T_SET(update, 1, 0, 0, 0, 0, 0, 0)( (update)[0] = (1), (update)[1] = (0), (update)[2] = (0), (update
)[3] = (0), (update)[4] = (0), (update)[5] = (0), (update)[6]
= (0) );
469 for (jj=0; jj<(rotnoop ? 3u : 6u); jj++) {
470 len02 = TEN_T_DOT(igrt[1][jj], d02)( (igrt[1][jj])[1]*(d02)[1] + 2*(igrt[1][jj])[2]*(d02)[2] + 2
*(igrt[1][jj])[3]*(d02)[3] + (igrt[1][jj])[4]*(d02)[4] + 2*(igrt
[1][jj])[5]*(d02)[5] + (igrt[1][jj])[6]*(d02)[6] );
471 len24 = TEN_T_DOT(igrt[3][jj], d24)( (igrt[3][jj])[1]*(d24)[1] + 2*(igrt[3][jj])[2]*(d24)[2] + 2
*(igrt[3][jj])[3]*(d24)[3] + (igrt[3][jj])[4]*(d24)[4] + 2*(igrt
[3][jj])[5]*(d24)[5] + (igrt[3][jj])[6]*(d24)[6] );
472 correct = (len24 - len02)/2;
473 TEN_T_SCALE_INCR(update, correct*scl, igrt[2][jj])( (update)[0] = (igrt[2][jj])[0], (update)[1] += (correct*scl
)*(igrt[2][jj])[1], (update)[2] += (correct*scl)*(igrt[2][jj]
)[2], (update)[3] += (correct*scl)*(igrt[2][jj])[3], (update)
[4] += (correct*scl)*(igrt[2][jj])[4], (update)[5] += (correct
*scl)*(igrt[2][jj])[5], (update)[6] += (correct*scl)*(igrt[2]
[jj])[6]);
474 if (tip->verbose) {
475 fprintf(stderr__stderrp, "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__stderrp, "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__stderrp, "(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)( (igrt[2][0])[1]*(update)[1] + 2*(igrt[2][0])[2]*(update)[2]
+ 2*(igrt[2][0])[3]*(update)[3] + (igrt[2][0])[4]*(update)[4
] + 2*(igrt[2][0])[5]*(update)[5] + (igrt[2][0])[6]*(update)[
6] ),
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])( (avg)[0] = (((0.5))*(((tt[4])[0]) - ((tt[0])[0])) + ((tt[0]
)[0])), (avg)[1] = (((0.5))*(((tt[4])[1]) - ((tt[0])[1])) + (
(tt[0])[1])), (avg)[2] = (((0.5))*(((tt[4])[2]) - ((tt[0])[2]
)) + ((tt[0])[2])), (avg)[3] = (((0.5))*(((tt[4])[3]) - ((tt[
0])[3])) + ((tt[0])[3])), (avg)[4] = (((0.5))*(((tt[4])[4]) -
((tt[0])[4])) + ((tt[0])[4])), (avg)[5] = (((0.5))*(((tt[4])
[5]) - ((tt[0])[5])) + ((tt[0])[5])), (avg)[6] = (((0.5))*(((
tt[4])[6]) - ((tt[0])[6])) + ((tt[0])[6])));
492 TEN_T_SUB(diff, avg, tt[2])( (diff)[0] = ((avg)[0] + (tt[2])[0])/2.0, (diff)[1] = (avg)[
1] - (tt[2])[1], (diff)[2] = (avg)[2] - (tt[2])[2], (diff)[3]
= (avg)[3] - (tt[2])[3], (diff)[4] = (avg)[4] - (tt[2])[4], (
diff)[5] = (avg)[5] - (tt[2])[5], (diff)[6] = (avg)[6] - (tt[
2])[6]);
493 for (jj=0; jj<3; jj++) {
494 len = TEN_T_DOT(igrt[2][jj], diff)( (igrt[2][jj])[1]*(diff)[1] + 2*(igrt[2][jj])[2]*(diff)[2] +
2*(igrt[2][jj])[3]*(diff)[3] + (igrt[2][jj])[4]*(diff)[4] + 2
*(igrt[2][jj])[5]*(diff)[5] + (igrt[2][jj])[6]*(diff)[6] );
495 TEN_T_SCALE_INCR(diff, -len, igrt[2][jj])( (diff)[0] = (igrt[2][jj])[0], (diff)[1] += (-len)*(igrt[2][
jj])[1], (diff)[2] += (-len)*(igrt[2][jj])[2], (diff)[3] += (
-len)*(igrt[2][jj])[3], (diff)[4] += (-len)*(igrt[2][jj])[4],
(diff)[5] += (-len)*(igrt[2][jj])[5], (diff)[6] += (-len)*(igrt
[2][jj])[6]);
496 }
497 TEN_T_SCALE_INCR(update, scl*0.2, diff)( (update)[0] = (diff)[0], (update)[1] += (scl*0.2)*(diff)[1]
, (update)[2] += (scl*0.2)*(diff)[2], (update)[3] += (scl*0.2
)*(diff)[3], (update)[4] += (scl*0.2)*(diff)[4], (update)[5] +=
(scl*0.2)*(diff)[5], (update)[6] += (scl*0.2)*(diff)[6]); /* HEY: scaling is a hack */
498 if (tip->verbose) {
499 fprintf(stderr__stderrp, "(rotnoop) (d = %g) "
500 "update = %g %g %g %g %g %g\n",
501 TEN_T_DOT(igrt[2][0], update)( (igrt[2][0])[1]*(update)[1] + 2*(igrt[2][0])[2]*(update)[2]
+ 2*(igrt[2][0])[3]*(update)[3] + (igrt[2][0])[4]*(update)[4
] + 2*(igrt[2][0])[5]*(update)[5] + (igrt[2][0])[6]*(update)[
6] ),
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)( (((int)(!(((update)[0]) - ((update)[0]))))) && (((int
)(!(((update)[1]) - ((update)[1]))))) && (((int)(!(((
update)[2]) - ((update)[2]))))) && (((int)(!(((update
)[3]) - ((update)[3]))))) && (((int)(!(((update)[4]) -
((update)[4]))))) && (((int)(!(((update)[5]) - ((update
)[5]))))) && (((int)(!(((update)[6]) - ((update)[6]))
))) )) {
516 biffAddf(TENtenBiffKey, "%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)( (odata + 7*(2*ii + 0))[0] = ((tt[2])[0] + (update)[0])/2.0,
(odata + 7*(2*ii + 0))[1] = (tt[2])[1] + (update)[1], (odata
+ 7*(2*ii + 0))[2] = (tt[2])[2] + (update)[2], (odata + 7*(2
*ii + 0))[3] = (tt[2])[3] + (update)[3], (odata + 7*(2*ii + 0
))[4] = (tt[2])[4] + (update)[4], (odata + 7*(2*ii + 0))[5] =
(tt[2])[5] + (update)[5], (odata + 7*(2*ii + 0))[6] = (tt[2]
)[6] + (update)[6]);
521
522 return 0;
523}
524
525void
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,( (igrt + 7*3)[0] = (1), (igrt + 7*3)[1] = ((airFloatQNaN.f))
, (igrt + 7*3)[2] = ((airFloatQNaN.f)), (igrt + 7*3)[3] = ((airFloatQNaN
.f)), (igrt + 7*3)[4] = ((airFloatQNaN.f)), (igrt + 7*3)[5] =
((airFloatQNaN.f)), (igrt + 7*3)[6] = ((airFloatQNaN.f)) )
539 AIR_NAN, AIR_NAN, AIR_NAN)( (igrt + 7*3)[0] = (1), (igrt + 7*3)[1] = ((airFloatQNaN.f))
, (igrt + 7*3)[2] = ((airFloatQNaN.f)), (igrt + 7*3)[3] = ((airFloatQNaN
.f)), (igrt + 7*3)[4] = ((airFloatQNaN.f)), (igrt + 7*3)[5] =
((airFloatQNaN.f)), (igrt + 7*3)[6] = ((airFloatQNaN.f)) );
540 TEN_T_SET(igrt + 7*4, 1, AIR_NAN, AIR_NAN, AIR_NAN,( (igrt + 7*4)[0] = (1), (igrt + 7*4)[1] = ((airFloatQNaN.f))
, (igrt + 7*4)[2] = ((airFloatQNaN.f)), (igrt + 7*4)[3] = ((airFloatQNaN
.f)), (igrt + 7*4)[4] = ((airFloatQNaN.f)), (igrt + 7*4)[5] =
((airFloatQNaN.f)), (igrt + 7*4)[6] = ((airFloatQNaN.f)) )
541 AIR_NAN, AIR_NAN, AIR_NAN)( (igrt + 7*4)[0] = (1), (igrt + 7*4)[1] = ((airFloatQNaN.f))
, (igrt + 7*4)[2] = ((airFloatQNaN.f)), (igrt + 7*4)[3] = ((airFloatQNaN
.f)), (igrt + 7*4)[4] = ((airFloatQNaN.f)), (igrt + 7*4)[5] =
((airFloatQNaN.f)), (igrt + 7*4)[6] = ((airFloatQNaN.f)) );
542 TEN_T_SET(igrt + 7*5, 1, AIR_NAN, AIR_NAN, AIR_NAN,( (igrt + 7*5)[0] = (1), (igrt + 7*5)[1] = ((airFloatQNaN.f))
, (igrt + 7*5)[2] = ((airFloatQNaN.f)), (igrt + 7*5)[3] = ((airFloatQNaN
.f)), (igrt + 7*5)[4] = ((airFloatQNaN.f)), (igrt + 7*5)[5] =
((airFloatQNaN.f)), (igrt + 7*5)[6] = ((airFloatQNaN.f)) )
543 AIR_NAN, AIR_NAN, AIR_NAN)( (igrt + 7*5)[0] = (1), (igrt + 7*5)[1] = ((airFloatQNaN.f))
, (igrt + 7*5)[2] = ((airFloatQNaN.f)), (igrt + 7*5)[3] = ((airFloatQNaN
.f)), (igrt + 7*5)[4] = ((airFloatQNaN.f)), (igrt + 7*5)[5] =
((airFloatQNaN.f)), (igrt + 7*5)[6] = ((airFloatQNaN.f)) );
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*/
558double
559tenInterpPathLength(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)((double *)(ntt->data));
564 if (doubleVerts) {
565 NN = AIR_CAST(unsigned int, (ntt->axis[1].size-1)/2)((unsigned int)((ntt->axis[1].size-1)/2));
566 } else {
567 NN = AIR_CAST(unsigned int, ntt->axis[1].size-1)((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)( (diff)[0] = ((tenA)[0] + (tenB)[0])/2.0, (diff)[1] = (tenA)
[1] - (tenB)[1], (diff)[2] = (tenA)[2] - (tenB)[2], (diff)[3]
= (tenA)[3] - (tenB)[3], (diff)[4] = (tenA)[4] - (tenB)[4], (
diff)[5] = (tenA)[5] - (tenB)[5], (diff)[6] = (tenA)[6] - (tenB
)[6]);
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)( (mean)[0] = (((0.5))*(((tenB)[0]) - ((tenA)[0])) + ((tenA)[
0])), (mean)[1] = (((0.5))*(((tenB)[1]) - ((tenA)[1])) + ((tenA
)[1])), (mean)[2] = (((0.5))*(((tenB)[2]) - ((tenA)[2])) + ((
tenA)[2])), (mean)[3] = (((0.5))*(((tenB)[3]) - ((tenA)[3])) +
((tenA)[3])), (mean)[4] = (((0.5))*(((tenB)[4]) - ((tenA)[4]
)) + ((tenA)[4])), (mean)[5] = (((0.5))*(((tenB)[5]) - ((tenA
)[5])) + ((tenA)[5])), (mean)[6] = (((0.5))*(((tenB)[6]) - ((
tenA)[6])) + ((tenA)[6])));
584 _tenInterpGeoLoxIGRT(igrt, mean, AIR_FALSE0, AIR_FALSE0, 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)( (igrt + 7*ii)[1]*(diff)[1] + 2*(igrt + 7*ii)[2]*(diff)[2] +
2*(igrt + 7*ii)[3]*(diff)[3] + (igrt + 7*ii)[4]*(diff)[4] + 2
*(igrt + 7*ii)[5]*(diff)[5] + (igrt + 7*ii)[6]*(diff)[6] );
595 incr += dot*dot;
596 }
597 len += sqrt(incr);
598 } else {
599 len += TEN_T_NORM(diff)(sqrt(( (diff)[1]*(diff)[1] + 2*(diff)[2]*(diff)[2] + 2*(diff
)[3]*(diff)[3] + (diff)[4]*(diff)[4] + 2*(diff)[5]*(diff)[5] +
(diff)[6]*(diff)[6] )));
600 }
601 }
602 return len;
603}
604
605double
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)((double *)(nin->data));
618 out = AIR_CAST(double *, nout->data)((double *)(nout->data));
619 NN = (nin->axis[1].size-1)/2;
620 lenTotal = tenInterpPathLength(nin, AIR_TRUE1, AIR_FALSE0, AIR_FALSE0);
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))( (out + 7*2*(0 + 0))[0] = (in + 7*2*(0 + 0))[0], (out + 7*2*
(0 + 0))[1] = (in + 7*2*(0 + 0))[1], (out + 7*2*(0 + 0))[2] =
(in + 7*2*(0 + 0))[2], (out + 7*2*(0 + 0))[3] = (in + 7*2*(0
+ 0))[3], (out + 7*2*(0 + 0))[4] = (in + 7*2*(0 + 0))[4], (out
+ 7*2*(0 + 0))[5] = (in + 7*2*(0 + 0))[5], (out + 7*2*(0 + 0
))[6] = (in + 7*2*(0 + 0))[6] );
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)( (diff)[0] = ((tenNext)[0] + (tenHere)[0])/2.0, (diff)[1] = (
tenNext)[1] - (tenHere)[1], (diff)[2] = (tenNext)[2] - (tenHere
)[2], (diff)[3] = (tenNext)[3] - (tenHere)[3], (diff)[4] = (tenNext
)[4] - (tenHere)[4], (diff)[5] = (tenNext)[5] - (tenHere)[5],
(diff)[6] = (tenNext)[6] - (tenHere)[6]);
633 lenHere = TEN_T_NORM(diff)(sqrt(( (diff)[1]*(diff)[1] + 2*(diff)[2]*(diff)[2] + 2*(diff
)[3]*(diff)[3] + (diff)[4]*(diff)[4] + 2*(diff)[5]*(diff)[5] +
(diff)[6]*(diff)[6] )));
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),( (out + 7*(2*idxOut + 0))[0] = ( ((double)((tenNext)[0])-((tenHere
)[0]))*((double)((len))-((lenHere))) / ((double)((0))-((lenHere
))) + ((tenHere)[0])), (out + 7*(2*idxOut + 0))[1] = ( ((double
)((tenNext)[1])-((tenHere)[1]))*((double)((len))-((lenHere)))
/ ((double)((0))-((lenHere))) + ((tenHere)[1])), (out + 7*(2
*idxOut + 0))[2] = ( ((double)((tenNext)[2])-((tenHere)[2]))*
((double)((len))-((lenHere))) / ((double)((0))-((lenHere))) +
((tenHere)[2])), (out + 7*(2*idxOut + 0))[3] = ( ((double)((
tenNext)[3])-((tenHere)[3]))*((double)((len))-((lenHere))) / (
(double)((0))-((lenHere))) + ((tenHere)[3])), (out + 7*(2*idxOut
+ 0))[4] = ( ((double)((tenNext)[4])-((tenHere)[4]))*((double
)((len))-((lenHere))) / ((double)((0))-((lenHere))) + ((tenHere
)[4])), (out + 7*(2*idxOut + 0))[5] = ( ((double)((tenNext)[5
])-((tenHere)[5]))*((double)((len))-((lenHere))) / ((double)(
(0))-((lenHere))) + ((tenHere)[5])), (out + 7*(2*idxOut + 0))
[6] = ( ((double)((tenNext)[6])-((tenHere)[6]))*((double)((len
))-((lenHere))) / ((double)((0))-((lenHere))) + ((tenHere)[6]
)))
649 lenHere, len, 0, tenHere, tenNext)( (out + 7*(2*idxOut + 0))[0] = ( ((double)((tenNext)[0])-((tenHere
)[0]))*((double)((len))-((lenHere))) / ((double)((0))-((lenHere
))) + ((tenHere)[0])), (out + 7*(2*idxOut + 0))[1] = ( ((double
)((tenNext)[1])-((tenHere)[1]))*((double)((len))-((lenHere)))
/ ((double)((0))-((lenHere))) + ((tenHere)[1])), (out + 7*(2
*idxOut + 0))[2] = ( ((double)((tenNext)[2])-((tenHere)[2]))*
((double)((len))-((lenHere))) / ((double)((0))-((lenHere))) +
((tenHere)[2])), (out + 7*(2*idxOut + 0))[3] = ( ((double)((
tenNext)[3])-((tenHere)[3]))*((double)((len))-((lenHere))) / (
(double)((0))-((lenHere))) + ((tenHere)[3])), (out + 7*(2*idxOut
+ 0))[4] = ( ((double)((tenNext)[4])-((tenHere)[4]))*((double
)((len))-((lenHere))) / ((double)((0))-((lenHere))) + ((tenHere
)[4])), (out + 7*(2*idxOut + 0))[5] = ( ((double)((tenNext)[5
])-((tenHere)[5]))*((double)((len))-((lenHere))) / ((double)(
(0))-((lenHere))) + ((tenHere)[5])), (out + 7*(2*idxOut + 0))
[6] = ( ((double)((tenNext)[6])-((tenHere)[6]))*((double)((len
))-((lenHere))) / ((double)((0))-((lenHere))) + ((tenHere)[6]
)));
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))( (out + 7*2*(NN + 0))[0] = (in + 7*2*(NN + 0))[0], (out + 7*
2*(NN + 0))[1] = (in + 7*2*(NN + 0))[1], (out + 7*2*(NN + 0))
[2] = (in + 7*2*(NN + 0))[2], (out + 7*2*(NN + 0))[3] = (in +
7*2*(NN + 0))[3], (out + 7*2*(NN + 0))[4] = (in + 7*2*(NN + 0
))[4], (out + 7*2*(NN + 0))[5] = (in + 7*2*(NN + 0))[5], (out
+ 7*2*(NN + 0))[6] = (in + 7*2*(NN + 0))[6] );
669
670 /* fill in vertex mid-points */
671 for (idxOut=0; idxOut<NN; idxOut++) {
672 TEN_T_LERP(out + 7*(2*idxOut + 1),( (out + 7*(2*idxOut + 1))[0] = (((0.5))*(((out + 7*(2*idxOut
+ 2))[0]) - ((out + 7*(2*idxOut + 0))[0])) + ((out + 7*(2*idxOut
+ 0))[0])), (out + 7*(2*idxOut + 1))[1] = (((0.5))*(((out + 7
*(2*idxOut + 2))[1]) - ((out + 7*(2*idxOut + 0))[1])) + ((out
+ 7*(2*idxOut + 0))[1])), (out + 7*(2*idxOut + 1))[2] = (((0.5
))*(((out + 7*(2*idxOut + 2))[2]) - ((out + 7*(2*idxOut + 0))
[2])) + ((out + 7*(2*idxOut + 0))[2])), (out + 7*(2*idxOut + 1
))[3] = (((0.5))*(((out + 7*(2*idxOut + 2))[3]) - ((out + 7*(
2*idxOut + 0))[3])) + ((out + 7*(2*idxOut + 0))[3])), (out + 7
*(2*idxOut + 1))[4] = (((0.5))*(((out + 7*(2*idxOut + 2))[4])
- ((out + 7*(2*idxOut + 0))[4])) + ((out + 7*(2*idxOut + 0))
[4])), (out + 7*(2*idxOut + 1))[5] = (((0.5))*(((out + 7*(2*idxOut
+ 2))[5]) - ((out + 7*(2*idxOut + 0))[5])) + ((out + 7*(2*idxOut
+ 0))[5])), (out + 7*(2*idxOut + 1))[6] = (((0.5))*(((out + 7
*(2*idxOut + 2))[6]) - ((out + 7*(2*idxOut + 0))[6])) + ((out
+ 7*(2*idxOut + 0))[6])))
673 0.5, out + 7*(2*idxOut + 0), out + 7*(2*idxOut + 2))( (out + 7*(2*idxOut + 1))[0] = (((0.5))*(((out + 7*(2*idxOut
+ 2))[0]) - ((out + 7*(2*idxOut + 0))[0])) + ((out + 7*(2*idxOut
+ 0))[0])), (out + 7*(2*idxOut + 1))[1] = (((0.5))*(((out + 7
*(2*idxOut + 2))[1]) - ((out + 7*(2*idxOut + 0))[1])) + ((out
+ 7*(2*idxOut + 0))[1])), (out + 7*(2*idxOut + 1))[2] = (((0.5
))*(((out + 7*(2*idxOut + 2))[2]) - ((out + 7*(2*idxOut + 0))
[2])) + ((out + 7*(2*idxOut + 0))[2])), (out + 7*(2*idxOut + 1
))[3] = (((0.5))*(((out + 7*(2*idxOut + 2))[3]) - ((out + 7*(
2*idxOut + 0))[3])) + ((out + 7*(2*idxOut + 0))[3])), (out + 7
*(2*idxOut + 1))[4] = (((0.5))*(((out + 7*(2*idxOut + 2))[4])
- ((out + 7*(2*idxOut + 0))[4])) + ((out + 7*(2*idxOut + 0))
[4])), (out + 7*(2*idxOut + 1))[5] = (((0.5))*(((out + 7*(2*idxOut
+ 2))[5]) - ((out + 7*(2*idxOut + 0))[5])) + ((out + 7*(2*idxOut
+ 0))[5])), (out + 7*(2*idxOut + 1))[6] = (((0.5))*(((out + 7
*(2*idxOut + 2))[6]) - ((out + 7*(2*idxOut + 0))[6])) + ((out
+ 7*(2*idxOut + 0))[6])));
674 }
675 return lenTotal;
676}
677
678int
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(TENtenBiffKey, "%s: got NULL pointer", me);
691 return 1;
692 }
693 if (!(NN >= 2)) {
694 biffAddf(TENtenBiffKey, "%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)((size_t)(7)),
709 AIR_CAST(size_t, NN+1)((size_t)(NN+1)))
710 || nrrdMaybeAlloc_va(ntt, nrrdTypeDouble, 2,
711 AIR_CAST(size_t, 7)((size_t)(7)),
712 AIR_CAST(size_t, 2*NN + 1)((size_t)(2*NN + 1)))
713 || nrrdMaybeAlloc_va(nigrt, nrrdTypeDouble, 3,
714 AIR_CAST(size_t, 7)((size_t)(7)),
715 AIR_CAST(size_t, 6)((size_t)(6)),
716 AIR_CAST(size_t, 2*NN + 1)((size_t)(2*NN + 1)))) {
717 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output", me);
718 airMopError(mop); return 1;
719 }
720 geod = AIR_CAST(double *, ngeod->data)((double *)(ngeod->data));
721 tt = AIR_CAST(double *, ntt->data)((double *)(ntt->data));
722 igrt = AIR_CAST(double *, nigrt->data)((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(TENtenBiffKey, "%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_FALSE0;
740 if (!E) E |= nrrdResampleDefaultCenterSet(rsmc, nrrdCenterNode);
741 if (!E) E |= nrrdResampleInputSet(rsmc, nsub);
742 if (!E) E |= nrrdResampleKernelSet(rsmc, 0, NULL((void*)0), NULL((void*)0));
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_TRUE1);
749 if (!E) E |= nrrdResampleExecute(rsmc, ntt);
750 if (E) {
751 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%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)( (tt + 7*ii)[0] = ( ((double)((tenB)[0])-((tenA)[0]))*((double
)((ii))-((0))) / ((double)((2*NN))-((0))) + ((tenA)[0])), (tt
+ 7*ii)[1] = ( ((double)((tenB)[1])-((tenA)[1]))*((double)((
ii))-((0))) / ((double)((2*NN))-((0))) + ((tenA)[1])), (tt + 7
*ii)[2] = ( ((double)((tenB)[2])-((tenA)[2]))*((double)((ii))
-((0))) / ((double)((2*NN))-((0))) + ((tenA)[2])), (tt + 7*ii
)[3] = ( ((double)((tenB)[3])-((tenA)[3]))*((double)((ii))-((
0))) / ((double)((2*NN))-((0))) + ((tenA)[3])), (tt + 7*ii)[4
] = ( ((double)((tenB)[4])-((tenA)[4]))*((double)((ii))-((0))
) / ((double)((2*NN))-((0))) + ((tenA)[4])), (tt + 7*ii)[5] =
( ((double)((tenB)[5])-((tenA)[5]))*((double)((ii))-((0))) /
((double)((2*NN))-((0))) + ((tenA)[5])), (tt + 7*ii)[6] = ( (
(double)((tenB)[6])-((tenA)[6]))*((double)((ii))-((0))) / ((double
)((2*NN))-((0))) + ((tenA)[6])));
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_TRUE1, AIR_FALSE0, AIR_FALSE0);
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__stderrp, "%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(TENtenBiffKey, "%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) > 0.0f ? (newlen - len) : -(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)( (geod + 7*ii)[0] = (tt + 7*2*ii)[0], (geod + 7*ii)[1] = (tt
+ 7*2*ii)[1], (geod + 7*ii)[2] = (tt + 7*2*ii)[2], (geod + 7
*ii)[3] = (tt + 7*2*ii)[3], (geod + 7*ii)[4] = (tt + 7*2*ii)[
4], (geod + 7*ii)[5] = (tt + 7*2*ii)[5], (geod + 7*ii)[6] = (
tt + 7*2*ii)[6] );
806 }
807 /* values from outer-most recursion will stick */
808 tip->numIter = *numIter;
809 tip->convFinal = 2*AIR_ABS(newlen - len)((newlen - len) > 0.0f ? (newlen - len) : -(newlen - len))/(newlen + len);
810
811 airMopOkay(mop);
812 return 0;
813}
814
815int
816tenInterpTwoDiscrete_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(TENtenBiffKey, "%s: got NULL pointer", me);
828 return 1;
829 }
830 if (airEnumValCheck(tenInterpType, ptype)) {
831 biffAddf(TENtenBiffKey, "%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(TENtenBiffKey, "%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)((size_t)(7)),
849 AIR_CAST(size_t, num)((size_t)(num)))) {
850 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble allocating output", me);
851 airMopError(mop); return 1;
852 }
853 out = AIR_CAST(double *, nout->data)((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__stderrp, "!%s: useK = %d, rotnoop = %d\n", me, useK, rotnoop);
882 if (_tenInterpGeoLoxPolyLine(nout, &numIter,
883 tenA, tenB,
884 num, useK, rotnoop, tip)) {
885 biffAddf(TENtenBiffKey, "%s: trouble finding path", me);
886 airMopError(mop); return 1;
887 }
888 } else {
889 biffAddf(TENtenBiffKey, "%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
898double
899tenInterpDistanceTwo_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(airFloatQNaN.f);
911 }
912
913 mop = airMopNew();
914 switch (ptype) {
915 case tenInterpTypeLinear:
916 TEN_T_SUB(diff, tenA, tenB)( (diff)[0] = ((tenA)[0] + (tenB)[0])/2.0, (diff)[1] = (tenA)
[1] - (tenB)[1], (diff)[2] = (tenA)[2] - (tenB)[2], (diff)[3]
= (tenA)[3] - (tenB)[3], (diff)[4] = (tenA)[4] - (tenB)[4], (
diff)[5] = (tenA)[5] - (tenB)[5], (diff)[6] = (tenA)[6] - (tenB
)[6]);
917 ret = TEN_T_NORM(diff)(sqrt(( (diff)[1]*(diff)[1] + 2*(diff)[2]*(diff)[2] + 2*(diff
)[3]*(diff)[3] + (diff)[4]*(diff)[4] + 2*(diff)[5]*(diff)[5] +
(diff)[6]*(diff)[6] )));
918 break;
919 case tenInterpTypeLogLinear:
920 tenLogSingle_d(logA, tenA);
921 tenLogSingle_d(logB, tenB);
922 TEN_T_SUB(diff, logA, logB)( (diff)[0] = ((logA)[0] + (logB)[0])/2.0, (diff)[1] = (logA)
[1] - (logB)[1], (diff)[2] = (logA)[2] - (logB)[2], (diff)[3]
= (logA)[3] - (logB)[3], (diff)[4] = (logA)[4] - (logB)[4], (
diff)[5] = (logA)[5] - (logB)[5], (diff)[6] = (logA)[6] - (logB
)[6]);
923 ret = TEN_T_NORM(diff)(sqrt(( (diff)[1]*(diff)[1] + 2*(diff)[2]*(diff)[2] + 2*(diff
)[3]*(diff)[3] + (diff)[4]*(diff)[4] + 2*(diff)[5]*(diff)[5] +
(diff)[6]*(diff)[6] )));
924 break;
925 case tenInterpTypeAffineInvariant:
926 TEN_T_INV(invA, tenA, det)((det) = ( (tenA)[1]*((tenA)[4]*(tenA)[6] - (tenA)[5]*(tenA)[
5]) + (tenA)[2]*((tenA)[5]*(tenA)[3] - (tenA)[2]*(tenA)[6]) +
(tenA)[3]*((tenA)[2]*(tenA)[5] - (tenA)[3]*(tenA)[4])), (invA
)[0] = (tenA)[0], (invA)[1] = (((tenA)[4])*((tenA)[6]) - ((tenA
)[5])*((tenA)[5]))/(det), (invA)[2] = -(((tenA)[2])*((tenA)[6
]) - ((tenA)[5])*((tenA)[3]))/(det), (invA)[3] = (((tenA)[2])
*((tenA)[5]) - ((tenA)[4])*((tenA)[3]))/(det), (invA)[4] = ((
(tenA)[1])*((tenA)[6]) - ((tenA)[3])*((tenA)[3]))/(det), (invA
)[5] = -(((tenA)[1])*((tenA)[5]) - ((tenA)[2])*((tenA)[3]))/(
det), (invA)[6] = (((tenA)[1])*((tenA)[4]) - ((tenA)[2])*((tenA
)[2]))/(det));
927 tenSqrtSingle_d(siA, invA);
928 TEN_T2M(mat1, tenB)( (mat1)[0] = (tenB)[1], (mat1)[1] = (tenB)[2], (mat1)[2] = (
tenB)[3], (mat1)[3] = (tenB)[2], (mat1)[4] = (tenB)[4], (mat1
)[5] = (tenB)[5], (mat1)[6] = (tenB)[3], (mat1)[7] = (tenB)[5
], (mat1)[8] = (tenB)[6] );
929 TEN_T2M(mat2, siA)( (mat2)[0] = (siA)[1], (mat2)[1] = (siA)[2], (mat2)[2] = (siA
)[3], (mat2)[3] = (siA)[2], (mat2)[4] = (siA)[4], (mat2)[5] =
(siA)[5], (mat2)[6] = (siA)[3], (mat2)[7] = (siA)[5], (mat2)
[8] = (siA)[6] );
930 ell_3m_mul_d(mat3, mat1, mat2);
931 ell_3m_mul_d(mat1, mat2, mat3);
932 TEN_M2T(diff, mat1)( (diff)[1] = (mat1)[0], (diff)[2] = ((mat1)[1]+(mat1)[3])/2.0
, (diff)[3] = ((mat1)[2]+(mat1)[6])/2.0, (diff)[4] = (mat1)[4
], (diff)[5] = ((mat1)[5]+(mat1)[7])/2.0, (diff)[6] = (mat1)[
8] );
933 tenLogSingle_d(logDiff, diff);
934 ret = TEN_T_NORM(logDiff)(sqrt(( (logDiff)[1]*(logDiff)[1] + 2*(logDiff)[2]*(logDiff)[
2] + 2*(logDiff)[3]*(logDiff)[3] + (logDiff)[4]*(logDiff)[4] +
2*(logDiff)[5]*(logDiff)[5] + (logDiff)[6]*(logDiff)[6] )));
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(TENtenBiffKey), airFree, airMopAlways);
953 fprintf(stderr__stderrp, "%s: trouble computing path:\n%s\n", me, err);
954 airMopError(mop); return AIR_NAN(airFloatQNaN.f);
955 }
956 ret = tenInterpPathLength(npath, AIR_FALSE0, AIR_FALSE0, AIR_FALSE0);
957 if (tip->lengthFancy) {
958 tip->lengthShape = tenInterpPathLength(npath, AIR_FALSE0,
959 AIR_TRUE1, AIR_TRUE1);
960 tip->lengthOrient = tenInterpPathLength(npath, AIR_FALSE0,
961 AIR_TRUE1, AIR_FALSE0);
962 }
963 break;
964 case tenInterpTypeWang:
965 default:
966 fprintf(stderr__stderrp, "%s: unimplemented %s %d!!!!\n", me,
967 tenInterpType->name, ptype);
968 ret = AIR_NAN(airFloatQNaN.f);
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*/
979int
980tenInterpMulti3D(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(TENtenBiffKey, "%s: got NULL pointer", me);
993 return 1;
994 }
995 if (!( ninLen > 0 )) {
996 biffAddf(TENtenBiffKey, "%s: need at least 1 nin, not 0", me);
997 return 1;
998 }
999 if (airEnumValCheck(tenInterpType, ptype)) {
1000 biffAddf(TENtenBiffKey, "%s: invalid %s %d", me,
1001 tenInterpType->name, ptype);
1002 return 1;
1003 }
1004
1005 if (tenTensorCheck(nin[0], nrrdTypeDefault, AIR_FALSE0, AIR_TRUE1)) {
1006 biffAddf(TENtenBiffKey, "%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(TENtenBiffKey, "%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_FALSE0, AIR_TRUE1)) {
1019 biffAddf(TENtenBiffKey, "%s: nin[%u] not a tensor array", me, ninIdx);
1020 return 1;
1021 }
1022 if (!nrrdSameSize(nin[0], nin[ninIdx], AIR_TRUE1)) {
1023 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: nin[0] doesn't match nin[%u]", me, ninIdx);
1024 return 1;
1025 }
1026 if (nin[0]->type != nin[ninIdx]->type) {
1027 biffAddf(TENtenBiffKey, "%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(TENtenBiffKey, NRRDnrrdBiffKey, "%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)))((double *)(calloc(7*ninLen, sizeof(double))));
1046 if (!tbuff) {
1047 biffAddf(TENtenBiffKey, "%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(128+1)];
1063 biffAddf(TENtenBiffKey, "%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}