1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 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 #include "ten.h"
25 #include "privateTen.h"
26
27 #if TEEM_LEVMAR
28 #include <levmar.h>
29 #endif
30
31 /* --------------------------------------------------------------------- */
32
33 const char *
34 _tenDwiGageStr[] = {
35 "(unknown tenDwiGage)",
36 "all",
37 "b0",
38 "jdwi",
39 "adc",
40 "mdwi",
41 "tlls",
42 "tllserr",
43 "tllserrlog",
44 "tllslike",
45 "twls",
46 "twlserr",
47 "twlserrlog",
48 "twlslike",
49 "tnls",
50 "tnlserr",
51 "tnlserrlog",
52 "tnlslike",
53 "tmle",
54 "tmleerr",
55 "tmleerrlog",
56 "tmlelike",
57 "t",
58 "terr",
59 "terrlog",
60 "tlike",
61 "c",
62 "fa",
63 "adwie",
64 "2qs",
65 "2qserr",
66 "2qsnerr",
67 "2peled",
68 "2pelederr",
69 "2pelednerr",
70 "2peledlminfo",
71 };
72
73 const int
74 _tenDwiGageVal[] = {
75 tenDwiGageUnknown,
76 tenDwiGageAll,
77 tenDwiGageB0,
78 tenDwiGageJustDWI,
79 tenDwiGageADC,
80 tenDwiGageMeanDWIValue,
81 tenDwiGageTensorLLS,
82 tenDwiGageTensorLLSError,
83 tenDwiGageTensorLLSErrorLog,
84 tenDwiGageTensorLLSLikelihood,
85 tenDwiGageTensorWLS,
86 tenDwiGageTensorWLSError,
87 tenDwiGageTensorWLSErrorLog,
88 tenDwiGageTensorWLSLikelihood,
89 tenDwiGageTensorNLS,
90 tenDwiGageTensorNLSError,
91 tenDwiGageTensorNLSErrorLog,
92 tenDwiGageTensorNLSLikelihood,
93 tenDwiGageTensorMLE,
94 tenDwiGageTensorMLEError,
95 tenDwiGageTensorMLEErrorLog,
96 tenDwiGageTensorMLELikelihood,
97 tenDwiGageTensor,
98 tenDwiGageTensorError,
99 tenDwiGageTensorErrorLog,
100 tenDwiGageTensorLikelihood,
101 tenDwiGageConfidence,
102 tenDwiGageFA,
103 tenDwiGageTensorAllDWIError,
104 tenDwiGage2TensorQSeg,
105 tenDwiGage2TensorQSegError,
106 tenDwiGage2TensorQSegAndError,
107 tenDwiGage2TensorPeled,
108 tenDwiGage2TensorPeledError,
109 tenDwiGage2TensorPeledAndError,
110 tenDwiGage2TensorPeledLevmarInfo
111 };
112
113 const airEnum
114 _tenDwiGage = {
115 "tenDwiGage",
116 TEN_DWI_GAGE_ITEM_MAX,
117 _tenDwiGageStr, _tenDwiGageVal,
118 NULL,
119 NULL, NULL,
120 AIR_FALSE
121 };
122 const airEnum *const
123 tenDwiGage = &_tenDwiGage;
124
125 /* --------------------------------------------------------------------- */
126
127 gageItemEntry
128 _tenDwiGageTable[TEN_DWI_GAGE_ITEM_MAX+1] = {
129 /* enum value len,deriv, prereqs, parent item, parent index, needData */
130 {tenDwiGageUnknown, 0, 0, {0}, 0, 0, AIR_TRUE},
131 /* len == 0 for tenDwiGage{All,JustDWI,ADC} means "learn later at run-time" */
132 {tenDwiGageAll, 0, 0, {0}, 0, 0, AIR_TRUE},
133 {tenDwiGageB0, 1, 0, {tenDwiGageAll}, tenDwiGageAll, 0, AIR_TRUE},
134 {tenDwiGageJustDWI, 0, 0, {tenDwiGageAll}, tenDwiGageAll, 1, AIR_TRUE},
135 {tenDwiGageADC, 0, 0, {tenDwiGageB0, tenDwiGageJustDWI}, 0, 0, AIR_TRUE},
136
137 {tenDwiGageMeanDWIValue, 1, 0, {tenDwiGageAll}, 0, 0, AIR_TRUE},
138
139 {tenDwiGageTensorLLS, 7, 0, {tenDwiGageAll, tenDwiGageMeanDWIValue}, 0, 0, AIR_TRUE},
140 {tenDwiGageTensorLLSError, 1, 0, {tenDwiGageTensorLLS}, 0, 0, AIR_TRUE},
141 {tenDwiGageTensorLLSErrorLog, 1, 0, {tenDwiGageTensorLLS}, 0, 0, AIR_TRUE},
142 {tenDwiGageTensorLLSLikelihood, 1, 0, {tenDwiGageTensorLLS}, 0, 0, AIR_TRUE},
143
144 {tenDwiGageTensorWLS, 7, 0, {tenDwiGageAll, tenDwiGageMeanDWIValue}, 0, 0, AIR_TRUE},
145 {tenDwiGageTensorWLSError, 1, 0, {tenDwiGageTensorWLS}, 0, 0, AIR_TRUE},
146 {tenDwiGageTensorWLSErrorLog, 1, 0, {tenDwiGageTensorWLS}, 0, 0, AIR_TRUE},
147 {tenDwiGageTensorWLSLikelihood, 1, 0, {tenDwiGageTensorWLS}, 0, 0, AIR_TRUE},
148
149 {tenDwiGageTensorNLS, 7, 0, {tenDwiGageAll, tenDwiGageMeanDWIValue}, 0, 0, AIR_TRUE},
150 {tenDwiGageTensorNLSError, 1, 0, {tenDwiGageTensorNLS}, 0, 0, AIR_TRUE},
151 {tenDwiGageTensorNLSErrorLog, 1, 0, {tenDwiGageTensorNLS}, 0, 0, AIR_TRUE},
152 {tenDwiGageTensorNLSLikelihood, 1, 0, {tenDwiGageTensorNLS}, 0, 0, AIR_TRUE},
153
154 {tenDwiGageTensorMLE, 7, 0, {tenDwiGageAll, tenDwiGageMeanDWIValue}, 0, 0, AIR_TRUE},
155 {tenDwiGageTensorMLEError, 1, 0, {tenDwiGageTensorMLE}, 0, 0, AIR_TRUE},
156 {tenDwiGageTensorMLEErrorLog, 1, 0, {tenDwiGageTensorMLE}, 0, 0, AIR_TRUE},
157 {tenDwiGageTensorMLELikelihood, 1, 0, {tenDwiGageTensorMLE}, 0, 0, AIR_TRUE},
158
159 /* these are NOT sub-items: they are copies, as controlled by the
160 kind->data, but not the query: the query can't capture the kind
161 of dependency implemented by having a dynamic kind */
162 {tenDwiGageTensor, 7, 0, {0}, /* 0 == "learn later at run time" */ 0, 0, AIR_TRUE},
163 {tenDwiGageTensorError, 1, 0, {0}, 0, 0, AIR_TRUE},
164 {tenDwiGageTensorErrorLog, 1, 0, {0}, 0, 0, AIR_TRUE},
165 {tenDwiGageTensorLikelihood, 1, 0, {0}, 0, 0, AIR_TRUE},
166
167 /* back to normal non-run-time items */
168 {tenDwiGageConfidence, 1, 0, {tenDwiGageTensor}, tenDwiGageTensor, 0, AIR_TRUE},
169 {tenDwiGageFA, 1, 0, {tenDwiGageTensor}, 0, 0, AIR_TRUE},
170
171 {tenDwiGageTensorAllDWIError, 0, 0, {tenDwiGageTensor, tenDwiGageJustDWI}, 0, 0, AIR_TRUE},
172
173 /* it actually doesn't make sense for tenDwiGage2TensorQSegAndError to be the parent,
174 because of the situations where you want the q-seg result, but don't care about error */
175 {tenDwiGage2TensorQSeg, 14, 0, {tenDwiGageAll}, 0, 0, AIR_TRUE},
176 {tenDwiGage2TensorQSegError, 1, 0, {tenDwiGageAll, tenDwiGage2TensorQSeg}, 0, 0, AIR_TRUE},
177 {tenDwiGage2TensorQSegAndError, 15, 0, {tenDwiGage2TensorQSeg, tenDwiGage2TensorQSegError}, 0, 0, AIR_TRUE},
178
179 {tenDwiGage2TensorPeled, 14, 0, {tenDwiGageAll}, 0, 0, AIR_TRUE},
180 {tenDwiGage2TensorPeledError, 1, 0, {tenDwiGageAll, tenDwiGage2TensorPeled}, 0, 0, AIR_TRUE},
181 {tenDwiGage2TensorPeledAndError, 15, 0, {tenDwiGage2TensorPeled, tenDwiGage2TensorPeledError}, 0, 0, AIR_TRUE},
182
183 {tenDwiGage2TensorPeledLevmarInfo, 5, 0, {tenDwiGage2TensorPeled}, 0, 0, AIR_TRUE}
184 };
185
186 void
_tenDwiGageIv3Print(FILE * file,gageContext * ctx,gagePerVolume * pvl)187 _tenDwiGageIv3Print(FILE *file, gageContext *ctx, gagePerVolume *pvl) {
188 static const char me[]="_tenDwiGageIv3Print";
189
190 AIR_UNUSED(ctx);
191 AIR_UNUSED(pvl);
192 fprintf(file, "%s: sorry, unimplemented\n", me);
193 return;
194 }
195
196 void
_tenDwiGageFilter(gageContext * ctx,gagePerVolume * pvl)197 _tenDwiGageFilter(gageContext *ctx, gagePerVolume *pvl) {
198 static const char me[]="_tenDwiGageFilter";
199 double *fw00, *fw11, *fw22, *dwi;
200 int fd, needD[3]={AIR_TRUE, AIR_FALSE, AIR_FALSE};
201 /* tenDwiGageKindData *kindData; */
202 gageScl3PFilter_t *filter[5] = {NULL, gageScl3PFilter2, gageScl3PFilter4,
203 gageScl3PFilter6, gageScl3PFilter8};
204 unsigned int J, dwiNum;
205
206 fd = 2*ctx->radius;
207 dwi = pvl->directAnswer[tenDwiGageAll];
208 /* kindData = AIR_CAST(tenDwiGageKindData *, pvl->kind->data); */
209 dwiNum = pvl->kind->valLen;
210 if (!ctx->parm.k3pack) {
211 fprintf(stderr, "%s: sorry, 6pack filtering not implemented\n", me);
212 return;
213 }
214 fw00 = ctx->fw + fd*3*gageKernel00;
215 fw11 = ctx->fw + fd*3*gageKernel11;
216 fw22 = ctx->fw + fd*3*gageKernel22;
217 /* HEY: these will have to be updated if there is ever any use for
218 derivatives in DWIs: can't pass NULL pointers for gradient info.
219 The unusual use of a hard-coded local needD is because there
220 currently isn't allocated space in the tenDwiGage kind (which is
221 unusual for its dynamic allocation) for DWI derivatives */
222 if (fd <= 8) {
223 for (J=0; J<dwiNum; J++) {
224 filter[ctx->radius](ctx->shape, pvl->iv3 + J*fd*fd*fd,
225 pvl->iv2 + J*fd*fd,
226 pvl->iv1 + J*fd,
227 fw00, fw11, fw22,
228 dwi + J, NULL, NULL,
229 needD);
230 }
231 } else {
232 for (J=0; J<dwiNum; J++) {
233 gageScl3PFilterN(ctx->shape, fd, pvl->iv3 + J*fd*fd*fd,
234 pvl->iv2 + J*fd*fd, pvl->iv1 + J*fd,
235 fw00, fw11, fw22,
236 dwi + J, NULL, NULL,
237 needD);
238 }
239 }
240
241 return;
242 }
243
244 /* Returns the Akaike Information Criterion */
245
246 /*
247 ** residual: is the variance
248 ** n: number of observations: number of DWI's in our case
249 ** k: number of parameters: number of tensor components in our case
250 */
251 double
_tenComputeAIC(double residual,int n,int k)252 _tenComputeAIC(double residual, int n, int k) {
253 double AIC = 0;
254
255 if (residual == 0) {
256 return 0;
257 }
258
259 /* AIC, RSS used when doing regression */
260 AIC = 2*k + n*log(residual);
261 /* Always use bias adjustment */
262 /* if (n/k < 40) { */
263 AIC = AIC + ((2*k*(k + 1))/(n - k - 1));
264 /* } */
265
266 return AIC;
267 }
268
269 /* Form a 2D tensor from the parameters */
270 void
_tenPeledRotate2D(double ten[7],double lam1,double lam3,double phi)271 _tenPeledRotate2D(double ten[7], double lam1, double lam3, double phi) {
272 double cc, ss, d3, d1, d2;
273
274 cc = cos(phi);
275 ss = sin(phi);
276 d1 = cc*cc*lam1 + ss*ss*lam3;
277 d3 = cc*ss*(lam1 - lam3);
278 d2 = ss*ss*lam1 + cc*cc*lam3;
279
280 TEN_T_SET(ten, 1.0, d1, d3, 0, d2, 0, lam3);
281 return;
282 }
283
284 /* The main callback function that is iterated during levmar */
285
286 /* vector pp of parameters is as follows:
287 ** pp[0]: principal eigenvalue
288 ** pp[1]: fraction of 1st tensor
289 ** pp[2]: phi for 1st tensor
290 ** pp[3]: phi for 2nd tensor
291 */
292 void
_tenLevmarPeledCB(double * pp,double * xx,int mm,int nn,void * _pvlData)293 _tenLevmarPeledCB(double *pp, double *xx, int mm, int nn, void *_pvlData) {
294 /* static const char me[]="_tenLevmarPeledCB"; */
295 double tenA[7], tenB[7];
296 int ii;
297 tenDwiGagePvlData *pvlData;
298 double *egrad;
299
300 AIR_UNUSED(mm);
301 pvlData = AIR_CAST(tenDwiGagePvlData *, _pvlData);
302
303 /* Form the tensors using the estimated parms */
304 _tenPeledRotate2D(tenA, pp[0], pvlData->ten1Eval[2], pp[2]);
305 _tenPeledRotate2D(tenB, pp[0], pvlData->ten1Eval[2], pp[3]);
306
307 egrad = AIR_CAST(double *, pvlData->nten1EigenGrads->data);
308 /* skip past b0 gradient, HEY: not general purpose */
309 egrad += 3;
310 for (ii=0; ii<nn; ii++) {
311 double argA, argB, sigA, sigB;
312 argA = -pvlData->tec1->bValue*TEN_T3V_CONTR(tenA, egrad + 3*ii);
313 argB = -pvlData->tec1->bValue*TEN_T3V_CONTR(tenB, egrad + 3*ii);
314 if (pvlData->levmarUseFastExp) {
315 sigA = airFastExp(argA);
316 sigB = airFastExp(argB);
317 } else {
318 sigA = exp(argA);
319 sigB = exp(argB);
320 }
321 xx[ii] = pvlData->tec1->knownB0*(pp[1]*sigA + (1-pp[1])*sigB);
322 }
323 return;
324 }
325
326 void
_tenDwiGageAnswer(gageContext * ctx,gagePerVolume * pvl)327 _tenDwiGageAnswer(gageContext *ctx, gagePerVolume *pvl) {
328 static const char me[]="_tenDwiGageAnswer";
329 unsigned int dwiIdx;
330 tenDwiGageKindData *kindData;
331 tenDwiGagePvlData *pvlData;
332 double *dwiAll, dwiMean=0, tentmp[7];
333
334 kindData = AIR_CAST(tenDwiGageKindData *, pvl->kind->data);
335 pvlData = AIR_CAST(tenDwiGagePvlData *, pvl->data);
336
337 dwiAll = pvl->directAnswer[tenDwiGageAll];
338 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageAll)) {
339 /* done if doV */
340 if (ctx->verbose) {
341 for (dwiIdx=0; dwiIdx<pvl->kind->valLen; dwiIdx++) {
342 fprintf(stderr, "%s(%d+%g,%d+%g,%d+%g): dwi[%u] = %g\n", me,
343 ctx->point.idx[0], ctx->point.frac[0],
344 ctx->point.idx[1], ctx->point.frac[1],
345 ctx->point.idx[2], ctx->point.frac[2],
346 dwiIdx, dwiAll[dwiIdx]);
347 }
348 fprintf(stderr, "%s: type(ngrad) = %d = %s\n", me,
349 kindData->ngrad->type,
350 airEnumStr(nrrdType, kindData->ngrad->type));
351 }
352 }
353
354 /*
355 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageB0)) {
356 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageJustDWI)) {
357 done if doV
358 }
359 */
360 /* HEY this isn't valid for multiple b-values */
361 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageADC)) {
362 double logdwi, logb0;
363 logb0 = log(AIR_MAX(kindData->valueMin,
364 pvl->directAnswer[tenDwiGageB0][0]));
365 for (dwiIdx=1; dwiIdx<pvl->kind->valLen; dwiIdx++) {
366 logdwi = log(AIR_MAX(kindData->valueMin,
367 pvl->directAnswer[tenDwiGageJustDWI][dwiIdx-1]));
368 pvl->directAnswer[tenDwiGageADC][dwiIdx-1]
369 = (logb0 - logdwi)/kindData->bval;
370 }
371 }
372 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageMeanDWIValue)) {
373 dwiMean = 0;
374 for (dwiIdx=1; dwiIdx<pvl->kind->valLen; dwiIdx++) {
375 dwiMean += dwiAll[dwiIdx];
376 }
377 dwiMean /= pvl->kind->valLen;
378 pvl->directAnswer[tenDwiGageMeanDWIValue][0] = dwiMean;
379 }
380
381 /* note: the gage interface to tenEstimate functionality
382 allows you exactly one kind of tensor estimation (per kind),
383 so the function call to do the estimation is actually
384 repeated over and over again; the copy into the answer
385 buffer is what changes... */
386 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorLLS)) {
387 tenEstimate1TensorSingle_d(pvlData->tec1, tentmp, dwiAll);
388 TEN_T_COPY(pvl->directAnswer[tenDwiGageTensorLLS], tentmp);
389 }
390 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorLLSError)) {
391 pvl->directAnswer[tenDwiGageTensorLLSError][0] = pvlData->tec1->errorDwi;
392 }
393 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorLLSErrorLog)) {
394 pvl->directAnswer[tenDwiGageTensorLLSErrorLog][0]
395 = pvlData->tec1->errorLogDwi;
396 }
397 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorWLS)) {
398 tenEstimate1TensorSingle_d(pvlData->tec1, tentmp, dwiAll);
399 TEN_T_COPY(pvl->directAnswer[tenDwiGageTensorWLS], tentmp);
400 }
401 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorNLS)) {
402 tenEstimate1TensorSingle_d(pvlData->tec1, tentmp, dwiAll);
403 TEN_T_COPY(pvl->directAnswer[tenDwiGageTensorNLS], tentmp);
404 }
405 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorMLE)) {
406 tenEstimate1TensorSingle_d(pvlData->tec1, tentmp, dwiAll);
407 TEN_T_COPY(pvl->directAnswer[tenDwiGageTensorMLE], tentmp);
408 }
409 /* HEY: have to implement all the different kinds of errors */
410
411 /* BEGIN sneakiness ........ */
412 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensor)) {
413 gageItemEntry *item;
414 item = pvl->kind->table + tenDwiGageTensor;
415 TEN_T_COPY(pvl->directAnswer[tenDwiGageTensor],
416 pvl->directAnswer[item->prereq[0]]);
417 }
418 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorError)) {
419 gageItemEntry *item;
420 item = pvl->kind->table + tenDwiGageTensorError;
421 pvl->directAnswer[tenDwiGageTensorError][0]
422 = pvl->directAnswer[item->prereq[0]][0];
423 }
424 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorErrorLog)) {
425 gageItemEntry *item;
426 item = pvl->kind->table + tenDwiGageTensorErrorLog;
427 pvl->directAnswer[tenDwiGageTensorErrorLog][0]
428 = pvl->directAnswer[item->prereq[0]][0];
429 }
430 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorLikelihood)) {
431 gageItemEntry *item;
432 item = pvl->kind->table + tenDwiGageTensorLikelihood;
433 pvl->directAnswer[tenDwiGageTensorLikelihood][0]
434 = pvl->directAnswer[item->prereq[0]][0];
435 }
436 /* END sneakiness ........ */
437
438 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageFA)) {
439 pvl->directAnswer[tenDwiGageFA][0]
440 = pvl->directAnswer[tenDwiGageTensor][0]
441 * tenAnisoTen_d(pvl->directAnswer[tenDwiGageTensor],
442 tenAniso_FA);
443 }
444
445 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorAllDWIError)) {
446 const double *grads;
447 int gradcount;
448 double *ten, d;
449 int i;
450
451 /* HEY: should switch to tenEstimate-based DWI simulation */
452 ten = pvl->directAnswer[tenDwiGageTensor];
453 gradcount = pvl->kind->valLen -1; /* Dont count b0 */
454 grads = ((const double*) kindData->ngrad->data) +3; /* Ignore b0 grad */
455 for( i=0; i < gradcount; i++ ) {
456 d = dwiAll[0]*exp(- pvlData->tec1->bValue
457 * TEN_T3V_CONTR(ten, grads + 3*i));
458 pvl->directAnswer[tenDwiGageTensorAllDWIError][i] = dwiAll[i+1] - d;
459 }
460 }
461
462 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGage2TensorQSeg)) {
463 const double *grads;
464 int gradcount;
465 double *twoten;
466 unsigned int valIdx, E;
467
468 twoten = pvl->directAnswer[tenDwiGage2TensorQSeg];
469
470 gradcount = pvl->kind->valLen -1; /* Dont count b0 */
471 grads = ((const double*) kindData->ngrad->data) +3; /* Ignore b0 grad */
472 if (dwiAll[0] != 0) { /* S0 = 0 */
473 _tenQball(pvlData->tec2->bValue, gradcount, dwiAll, grads,
474 pvlData->qvals);
475 _tenQvals2points(gradcount, pvlData->qvals, grads, pvlData->qpoints);
476 _tenSegsamp2(gradcount, pvlData->qvals, grads, pvlData->qpoints,
477 pvlData->wght + 1, pvlData->dists );
478 } else {
479 /* stupid; should really return right here since data is garbage */
480 for (valIdx=1; valIdx < AIR_CAST(unsigned int, gradcount+1); valIdx++) {
481 pvlData->wght[valIdx] = valIdx % 2;
482 }
483 }
484
485 E = 0;
486 for (valIdx=1; valIdx<pvl->kind->valLen; valIdx++) {
487 if (!E) E |= tenEstimateSkipSet(pvlData->tec2, valIdx,
488 pvlData->wght[valIdx]);
489 }
490 if (!E) E |= tenEstimateUpdate(pvlData->tec2);
491 if (!E) E |= tenEstimate1TensorSingle_d(pvlData->tec2,
492 twoten + 0, dwiAll);
493 for (valIdx=1; valIdx<pvl->kind->valLen; valIdx++) {
494 if (!E) E |= tenEstimateSkipSet(pvlData->tec2, valIdx,
495 1 - pvlData->wght[valIdx]);
496 }
497 if (!E) E |= tenEstimateUpdate(pvlData->tec2);
498 if (!E) E |= tenEstimate1TensorSingle_d(pvlData->tec2,
499 twoten + 7, dwiAll);
500 if (E) {
501 char *terr;
502 terr = biffGetDone(TEN);
503 fprintf(stderr, "%s: (trouble) %s\n", me, terr);
504 free(terr);
505 }
506
507 /* hack: confidence for two-tensor fit */
508 twoten[0] = (twoten[0] + twoten[7])/2;
509 twoten[7] = 0.5; /* fraction that is the first tensor (initial value) */
510 /* twoten[1 .. 6] = first tensor */
511 /* twoten[8 .. 13] = second tensor */
512
513 /* Compute fraction between tensors if not garbage in this voxel */
514 if (twoten[0] > 0.5) {
515 double exp0,exp1,d,e=0,g=0, a=0,b=0;
516 int i;
517
518 for( i=0; i < gradcount; i++ ) {
519 exp0 = exp(-pvlData->tec2->bValue * TEN_T3V_CONTR(twoten + 0,
520 grads + 3*i));
521 exp1 = exp(-pvlData->tec2->bValue * TEN_T3V_CONTR(twoten + 7,
522 grads + 3*i));
523
524 d = dwiAll[i+1] / dwiAll[0];
525 e = exp0 - exp1;
526 g = d - exp1;
527
528 a += .5*e*e;
529 b += e*g;
530 }
531
532 twoten[7] = AIR_CLAMP(0, 0.5*(b/a), 1);
533 }
534 }
535 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGage2TensorQSegError)) {
536 const double *grads;
537 int gradcount;
538 double *twoten, d;
539 int i;
540
541 /* HEY: should switch to tenEstimate-based DWI simulation */
542 if (dwiAll[0] != 0) { /* S0 = 0 */
543 twoten = pvl->directAnswer[tenDwiGage2TensorQSeg];
544 gradcount = pvl->kind->valLen -1; /* Dont count b0 */
545 grads = ((const double*) kindData->ngrad->data) +3; /* Ignore b0 grad */
546
547 pvl->directAnswer[tenDwiGage2TensorQSegError][0] = 0;
548 for( i=0; i < gradcount; i++ ) {
549 d = twoten[7]*exp(-pvlData->tec2->bValue * TEN_T3V_CONTR(twoten + 0,
550 grads + 3*i));
551 d += (1 - twoten[7])*exp(-pvlData->tec2->bValue
552 *TEN_T3V_CONTR(twoten + 7, grads + 3*i));
553 d = dwiAll[i+1]/dwiAll[0] - d;
554 pvl->directAnswer[tenDwiGage2TensorQSegError][0] += d*d;
555 }
556 pvl->directAnswer[tenDwiGage2TensorQSegError][0] =
557 sqrt( pvl->directAnswer[tenDwiGage2TensorQSegError][0] );
558 } else {
559 /* HEY: COMPLETELY WRONG!! An error is not defined! */
560 pvl->directAnswer[tenDwiGage2TensorQSegError][0] = 0;
561 }
562 /* printf("%f\n",pvl->directAnswer[tenDwiGage2TensorQSegError][0]); */
563 }
564 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGage2TensorQSegAndError)) {
565 double *twoten, *err, *twotenerr;
566
567 twoten = pvl->directAnswer[tenDwiGage2TensorQSeg];
568 err = pvl->directAnswer[tenDwiGage2TensorQSegError];
569 twotenerr = pvl->directAnswer[tenDwiGage2TensorQSegAndError];
570 TEN_T_COPY(twotenerr + 0, twoten + 0);
571 TEN_T_COPY(twotenerr + 7, twoten + 7);
572 twotenerr[14] = err[0];
573 }
574
575 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGage2TensorPeled)) {
576 #if TEEM_LEVMAR
577 #define PARAMS 4
578 double *twoTen, Cp /* , residual, AICSingFit, AICTwoFit */;
579 /* Vars for the NLLS */
580 double guess[PARAMS], loBnd[PARAMS], upBnd[PARAMS],
581 opts[LM_OPTS_SZ], *grad, *egrad, tenA[7], tenB[7],
582 matA[9], matB[9], matTmp[9], rott[9];
583 unsigned int gi;
584 int lmret;
585
586 /* Pointer to the location where the two tensor will be written */
587 twoTen = pvl->directAnswer[tenDwiGage2TensorPeled];
588 /* Estimate the DWI error, error is given as standard deviation */
589 pvlData->tec1->recordErrorDwi = AIR_FALSE;
590 /* Estimate the single tensor */
591 tenEstimate1TensorSingle_d(pvlData->tec1, pvlData->ten1, dwiAll);
592 /* Get the eigenValues and eigen vectors for this tensor */
593 tenEigensolve_d(pvlData->ten1Eval, pvlData->ten1Evec, pvlData->ten1);
594 /* Get westins Cp */
595 Cp = tenAnisoEval_d(pvlData->ten1Eval, tenAniso_Cp1);
596
597 /* Only do two-tensor fitting if CP is greater or equal to than a
598 user-defined threshold */
599 if (Cp >= pvlData->levmarMinCp) {
600 /* Calculate the residual, need the variance to sqr it */
601 /* residual = pvlData->tec1->errorDwi*pvlData->tec1->errorDwi; */
602 /* Calculate the AIC for single tensor fit */
603 /* AICSingFit = _tenComputeAIC(residual, pvlData->tec1->dwiNum, 6); */
604
605 /* the CP-based test is gone; caller's responsibility */
606
607 /* rotate DW gradients by inverse of eigenvector column matrix
608 and place into pvlData->nten1EigenGrads (which has been
609 allocated by _tenDwiGagePvlDataNew()) */
610 grad = AIR_CAST(double *, kindData->ngrad->data);
611 egrad = AIR_CAST(double *, pvlData->nten1EigenGrads->data);
612 for (gi=0; gi<kindData->ngrad->axis[1].size; gi++) {
613 /* yes, this is also transforming some zero-length (B0) gradients;
614 that's harmless */
615 ELL_3MV_MUL(egrad, pvlData->ten1Evec, grad);
616 grad += 3;
617 egrad += 3;
618 }
619
620 /* Lower and upper bounds for the NLLS routine */
621 loBnd[0] = 0.0;
622 loBnd[1] = 0.0;
623 loBnd[2] = -AIR_PI/2;
624 loBnd[3] = -AIR_PI/2;
625 upBnd[0] = pvlData->ten1Eval[0]*5;
626 upBnd[1] = 1.0;
627 upBnd[2] = AIR_PI/2;
628 upBnd[3] = AIR_PI/2;
629 /* Starting point for the NLLS */
630 guess[0] = pvlData->ten1Eval[0];
631 guess[1] = 0.5;
632
633 guess[2] = AIR_PI/4;
634 guess[3] = -AIR_PI/4;
635 /*
636 guess[2] = AIR_AFFINE(0, airDrandMT_r(pvlData->randState), 1,
637 AIR_PI/6, AIR_PI/3);
638 guess[3] = AIR_AFFINE(0, airDrandMT_r(pvlData->randState), 1,
639 -AIR_PI/6, -AIR_PI/3);
640 */
641 /* Fill in the constraints for the LM optimization,
642 the threshold of error difference */
643 opts[0] = pvlData->levmarTau;
644 opts[1] = pvlData->levmarEps1;
645 opts[2] = pvlData->levmarEps2;
646 opts[3] = pvlData->levmarEps3;
647 /* Very imp to set this opt, note that only forward
648 differences are used to approx Jacobian */
649 opts[4] = pvlData->levmarDelta;
650
651 /* run NLLS, results are stored back into guess[] */
652 pvlData->levmarUseFastExp = AIR_FALSE;
653 lmret = dlevmar_bc_dif(_tenLevmarPeledCB, guess, pvlData->tec1->dwi,
654 PARAMS, pvlData->tec1->dwiNum, loBnd, upBnd,
655 NULL, pvlData->levmarMaxIter, opts,
656 pvlData->levmarInfo,
657 NULL, NULL, pvlData);
658 if (-1 == lmret) {
659 ctx->errNum = 1;
660 sprintf(ctx->errStr, "%s: dlevmar_bc_dif() failed!", me);
661 } else {
662 /* Get the AIC for the two tensor fit, use the levmarinfo
663 to get the residual */
664 /*
665 residual = pvlData->levmarInfo[1]/pvlData->tec1->dwiNum;
666 AICTwoFit = _tenComputeAIC(residual, pvlData->tec1->dwiNum, 12);
667 */
668 /* Form the tensors using the estimated pp, returned in guess */
669 _tenPeledRotate2D(tenA, guess[0], pvlData->ten1Eval[2], guess[2]);
670 _tenPeledRotate2D(tenB, guess[0], pvlData->ten1Eval[2], guess[3]);
671 TEN_T2M(matA, tenA);
672 TEN_T2M(matB, tenB);
673
674 ELL_3M_TRANSPOSE(rott, pvlData->ten1Evec);
675 ELL_3M_MUL(matTmp, matA, pvlData->ten1Evec);
676 ELL_3M_MUL(matA, rott, matTmp);
677 ELL_3M_MUL(matTmp, matB, pvlData->ten1Evec);
678 ELL_3M_MUL(matB, rott, matTmp);
679
680 /* Copy two two tensors */
681 /* guess[1] is population fraction of first tensor */
682 if (guess[1] > 0.5) {
683 twoTen[7] = guess[1];
684 TEN_M2T(twoTen + 0, matA);
685 TEN_M2T(twoTen + 7, matB);
686 } else {
687 twoTen[7] = 1 - guess[1];
688 TEN_M2T(twoTen + 0, matB);
689 TEN_M2T(twoTen + 7, matA);
690 }
691 twoTen[0] = 1;
692 }
693 } else {
694 /* its too planar- just do single tensor fit */
695 TEN_T_COPY(twoTen, pvlData->ten1);
696 TEN_T_SET(twoTen + 7, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
697 }
698 #undef PARAMS
699 #else
700 double *twoTen;
701 twoTen = pvl->directAnswer[tenDwiGage2TensorPeled];
702 TEN_T_SET(twoTen + 0, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN,
703 AIR_NAN, AIR_NAN, AIR_NAN);
704 TEN_T_SET(twoTen + 7, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN,
705 AIR_NAN, AIR_NAN, AIR_NAN);
706 fprintf(stderr, "%s: sorry, not compiled with TEEM_LEVMAR\n", me);
707 #endif
708 }
709 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGage2TensorPeledError)) {
710 double *info;
711 info = pvlData->levmarInfo;
712 pvl->directAnswer[tenDwiGage2TensorPeledError][0] = 0;
713
714 if (info[1] > 0) {
715 /* Returning the standard deviation */
716 pvl->directAnswer[tenDwiGage2TensorPeledError][0] =
717 sqrt(info[1]/pvlData->tec1->dwiNum);
718 }
719 }
720 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGage2TensorPeledAndError)) {
721 double *twoten, *err, *twotenerr;
722 /* HEY cut and paste */
723 twoten = pvl->directAnswer[tenDwiGage2TensorPeled];
724 err = pvl->directAnswer[tenDwiGage2TensorPeledError];
725 twotenerr = pvl->directAnswer[tenDwiGage2TensorPeledAndError];
726 TEN_T_COPY(twotenerr + 0, twoten + 0);
727 TEN_T_COPY(twotenerr + 7, twoten + 7);
728 twotenerr[14] = err[0];
729 }
730
731 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGage2TensorPeledLevmarInfo)) {
732 double *info;
733 unsigned int ii, alen;
734 alen = gageKindAnswerLength(pvl->kind, tenDwiGage2TensorPeledLevmarInfo);
735 info = pvl->directAnswer[tenDwiGage2TensorPeledLevmarInfo];
736 for (ii=0; ii<alen; ii++) {
737 info[ii] = pvlData->levmarInfo[ii];
738 }
739 }
740
741 return;
742 }
743
744 /* --------------------- pvlData */
745
746 /* note use of the GAGE biff key */
747 void *
_tenDwiGagePvlDataNew(const gageKind * kind)748 _tenDwiGagePvlDataNew(const gageKind *kind) {
749 static const char me[]="_tenDwiGagePvlDataNew";
750 tenDwiGagePvlData *pvlData;
751 tenDwiGageKindData *kindData;
752 const int segcount = 2;
753 unsigned int num;
754 int E;
755
756 if (tenDwiGageKindCheck(kind)) {
757 biffMovef(GAGE, TEN, "%s: kindData not ready for use", me);
758 return NULL;
759 }
760 kindData = AIR_CAST(tenDwiGageKindData *, kind->data);
761
762 pvlData = AIR_CALLOC(1, tenDwiGagePvlData);
763 if (!pvlData) {
764 biffAddf(GAGE, "%s: couldn't allocate pvl data!", me);
765 return NULL;
766 }
767 pvlData->tec1 = tenEstimateContextNew();
768 pvlData->tec2 = tenEstimateContextNew();
769 for (num=1; num<=2; num++) {
770 tenEstimateContext *tec;
771 tec = (1 == num ? pvlData->tec1 : pvlData->tec2);
772 E = 0;
773 if (!E) tenEstimateVerboseSet(tec, 0);
774 if (!E) tenEstimateNegEvalShiftSet(tec, AIR_FALSE);
775 if (!E) E |= tenEstimateMethodSet(tec, 1 == num
776 ? kindData->est1Method
777 : kindData->est2Method);
778 if (!E) E |= tenEstimateValueMinSet(tec, kindData->valueMin);
779 if (kindData->ngrad->data) {
780 if (!E) E |= tenEstimateGradientsSet(tec, kindData->ngrad,
781 kindData->bval, AIR_FALSE);
782 } else {
783 if (!E) E |= tenEstimateBMatricesSet(tec, kindData->nbmat,
784 kindData->bval, AIR_FALSE);
785 }
786 if (!E) E |= tenEstimateThresholdSet(tec,
787 kindData->thresh, kindData->soft);
788 if (!E) E |= tenEstimateUpdate(tec);
789 if (E) {
790 biffMovef(GAGE, TEN, "%s: trouble setting %u estimation", me, num);
791 return NULL;
792 }
793 }
794 pvlData->vbuf = AIR_CALLOC(kind->valLen, double);
795 pvlData->wght = AIR_CALLOC(kind->valLen, unsigned int);
796 /* HEY: this is where we act on the the assumption about
797 having val[0] be T2 baseline and all subsequent val[i] be DWIs */
798 pvlData->wght[0] = 1;
799 pvlData->qvals = AIR_CALLOC(kind->valLen-1, double);
800 pvlData->qpoints = AIR_CALLOC(3*(kind->valLen-1), double);
801 pvlData->dists = AIR_CALLOC(segcount*(kind->valLen-1), double);
802 pvlData->weights = AIR_CALLOC(segcount*(kind->valLen-1), double);
803
804 if (kindData->ngrad->data) {
805 pvlData->nten1EigenGrads = nrrdNew();
806 /* this is for allocation only; values will get over-written */
807 nrrdCopy(pvlData->nten1EigenGrads, kindData->ngrad);
808 } else {
809 /* HEY: currently don't handle general B-matrices here */
810 pvlData->nten1EigenGrads = NULL;
811 }
812
813 pvlData->randSeed = kindData->randSeed;
814 pvlData->randState = airRandMTStateNew(pvlData->randSeed);
815
816 /* initialize single-tensor info to all NaNs */
817 TEN_T_SET(pvlData->ten1, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN,
818 AIR_NAN, AIR_NAN, AIR_NAN);
819 ELL_3V_SET(pvlData->ten1Evec + 0, AIR_NAN, AIR_NAN, AIR_NAN);
820 ELL_3V_SET(pvlData->ten1Evec + 3, AIR_NAN, AIR_NAN, AIR_NAN);
821 ELL_3V_SET(pvlData->ten1Evec + 6, AIR_NAN, AIR_NAN, AIR_NAN);
822 ELL_3V_SET(pvlData->ten1Eval, AIR_NAN, AIR_NAN, AIR_NAN);
823
824 /* here's an okay spot to check our compile-time assumptions
825 about the levmar library */
826 #if TEEM_LEVMAR
827 /* this is needed to make sure that the tenDwiGage2TensorPeledLevmarInfo
828 item definition above is valid */
829 if (5 != LM_OPTS_SZ) {
830 biffAddf(GAGE, "%s: LM_OPTS_SZ (%d) != expected 5\n", me, LM_OPTS_SZ);
831 return NULL;
832 }
833 #endif
834 pvlData->levmarUseFastExp = AIR_FALSE;
835 pvlData->levmarMaxIter = 200;
836 pvlData->levmarTau = 1E-03; /* LM_INIT_MU; */
837 pvlData->levmarEps1 = 1E-8;
838 pvlData->levmarEps2 = 1E-8;
839 pvlData->levmarEps3 = 1E-8;
840 pvlData->levmarDelta = 1E-8;
841 pvlData->levmarMinCp = 0.1;
842
843 /* pvlData->levmarInfo[] is output; not initialized */
844
845 return AIR_CAST(void *, pvlData);
846 }
847
848 void *
_tenDwiGagePvlDataCopy(const gageKind * kind,const void * _pvlDataOld)849 _tenDwiGagePvlDataCopy(const gageKind *kind, const void *_pvlDataOld) {
850 tenDwiGagePvlData *pvlDataOld, *pvlDataNew;
851
852 pvlDataOld = AIR_CAST(tenDwiGagePvlData *, _pvlDataOld);
853 pvlDataNew = AIR_CAST(tenDwiGagePvlData *, _tenDwiGagePvlDataNew(kind));
854
855 /* HEY: no error checking? */
856 if (pvlDataOld->nten1EigenGrads) {
857 nrrdCopy(pvlDataNew->nten1EigenGrads, pvlDataOld->nten1EigenGrads);
858 }
859 /* need to copy randState or randSeed? */
860
861 TEN_T_COPY(pvlDataNew->ten1, pvlDataOld->ten1);
862 ELL_3M_COPY(pvlDataNew->ten1Evec, pvlDataOld->ten1Evec);
863 ELL_3V_COPY(pvlDataNew->ten1Eval, pvlDataOld->ten1Eval);
864
865 pvlDataNew->levmarUseFastExp = pvlDataOld->levmarUseFastExp;
866 pvlDataNew->levmarMaxIter = pvlDataOld->levmarMaxIter;
867 pvlDataNew->levmarTau = pvlDataOld->levmarTau;
868 pvlDataNew->levmarEps1 = pvlDataOld->levmarEps1;
869 pvlDataNew->levmarEps2 = pvlDataOld->levmarEps2;
870 pvlDataNew->levmarEps3 = pvlDataOld->levmarEps3;
871 pvlDataNew->levmarDelta = pvlDataOld->levmarDelta;
872 pvlDataNew->levmarMinCp = pvlDataOld->levmarMinCp;
873 /* pvlData->levmarInfo[] is output; not copied */
874
875 return pvlDataNew;
876 }
877
878 int
_tenDwiGagePvlDataUpdate(const gageKind * kind,const gageContext * ctx,const gagePerVolume * pvl,const void * _pvlData)879 _tenDwiGagePvlDataUpdate(const gageKind *kind,
880 const gageContext *ctx,
881 const gagePerVolume *pvl, const void *_pvlData) {
882 /* static const char me[]="_tenDwiGagePvlDataUpdate"; */
883 tenDwiGagePvlData *pvlData;
884
885 AIR_UNUSED(ctx);
886 pvlData = AIR_CAST(tenDwiGagePvlData *, _pvlData);
887 AIR_UNUSED(kind);
888 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorLLSError)
889 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorWLSError)
890 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorNLSError)
891 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorMLEError)) {
892 pvlData->tec1->recordErrorDwi = AIR_TRUE;
893 } else {
894 pvlData->tec1->recordErrorDwi = AIR_FALSE;
895 }
896 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorLLSErrorLog)
897 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorWLSErrorLog)
898 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorNLSErrorLog)
899 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorMLEErrorLog)) {
900 pvlData->tec1->recordErrorLogDwi = AIR_TRUE;
901 } else {
902 pvlData->tec1->recordErrorLogDwi = AIR_FALSE;
903 }
904 if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorLLSLikelihood)
905 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorWLSLikelihood)
906 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorNLSLikelihood)
907 || GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGageTensorMLELikelihood)) {
908 pvlData->tec1->recordLikelihoodDwi = AIR_TRUE;
909 } else {
910 pvlData->tec1->recordLikelihoodDwi = AIR_FALSE;
911 }
912 /*
913 fprintf(stderr, "%s: record %d %d %d\n", me,
914 pvlData->tec1->recordErrorDwi,
915 pvlData->tec1->recordErrorLogDwi,
916 pvlData->tec1->recordLikelihoodDwi);
917 */
918 return 0;
919 }
920
921 void *
_tenDwiGagePvlDataNix(const gageKind * kind,void * _pvlData)922 _tenDwiGagePvlDataNix(const gageKind *kind, void *_pvlData) {
923 tenDwiGagePvlData *pvlData;
924
925 AIR_UNUSED(kind);
926 pvlData = AIR_CAST(tenDwiGagePvlData *, _pvlData);
927 if (pvlData) {
928 tenEstimateContextNix(pvlData->tec1);
929 tenEstimateContextNix(pvlData->tec2);
930 airFree(pvlData->vbuf);
931 airFree(pvlData->wght);
932 airFree(pvlData->qvals);
933 airFree(pvlData->qpoints);
934 airFree(pvlData->dists);
935 airFree(pvlData->weights);
936 nrrdNuke(pvlData->nten1EigenGrads);
937 airRandMTStateNix(pvlData->randState);
938 airFree(pvlData);
939 }
940 return NULL;
941 }
942
943 /* --------------------- kindData */
944
945 tenDwiGageKindData*
tenDwiGageKindDataNew(void)946 tenDwiGageKindDataNew(void) {
947 tenDwiGageKindData *ret;
948
949 ret = AIR_CALLOC(1, tenDwiGageKindData);
950 if (ret) {
951 /* it may be that only one of these is actually filled */
952 ret->ngrad = nrrdNew();
953 ret->nbmat = nrrdNew();
954 ret->thresh = ret->soft = ret->bval = AIR_NAN;
955 ret->est1Method = tenEstimate1MethodUnknown;
956 ret->est2Method = tenEstimate2MethodUnknown;
957 ret->randSeed = 42;
958 }
959 return ret;
960 }
961
962 tenDwiGageKindData*
tenDwiGageKindDataNix(tenDwiGageKindData * kindData)963 tenDwiGageKindDataNix(tenDwiGageKindData *kindData) {
964
965 if (kindData) {
966 nrrdNuke(kindData->ngrad);
967 nrrdNuke(kindData->nbmat);
968 airFree(kindData);
969 }
970 return NULL;
971 }
972
973 /* --------------------- dwiKind, and dwiKind->data setting*/
974
975 /*
976 ** Because this kind has to be dynamically allocated,
977 ** this is not the kind, but just the template for it
978 ** HEY: having a const public version of this could be a
979 ** nice way of having a way of referring to the dwiKind
980 ** without having to allocate it each time
981 */
982 gageKind
983 _tenDwiGageKindTmpl = {
984 AIR_TRUE, /* dynamically allocated */
985 TEN_DWI_GAGE_KIND_NAME,
986 &_tenDwiGage,
987 1,
988 0, /* NOT: set later by tenDwiGageKindSet() */
989 TEN_DWI_GAGE_ITEM_MAX,
990 NULL, /* NOT: modified copy of _tenDwiGageTable,
991 allocated by tenDwiGageKindNew(), and
992 set by _tenDwiGageKindSet() */
993 _tenDwiGageIv3Print,
994 _tenDwiGageFilter,
995 _tenDwiGageAnswer,
996 _tenDwiGagePvlDataNew,
997 _tenDwiGagePvlDataCopy,
998 _tenDwiGagePvlDataNix,
999 _tenDwiGagePvlDataUpdate,
1000 NULL /* NOT: allocated by tenDwiGageKindNew(),
1001 insides set by tenDwiGageKindSet() */
1002 };
1003
1004 gageKind *
tenDwiGageKindNew()1005 tenDwiGageKindNew() {
1006 gageKind *kind;
1007
1008 kind = AIR_CALLOC(1, gageKind);
1009 if (kind) {
1010 memcpy(kind, &_tenDwiGageKindTmpl, sizeof(gageKind));
1011 kind->valLen = 0; /* still has to be set later */
1012 kind->table = AIR_CAST(gageItemEntry *,
1013 malloc(sizeof(_tenDwiGageTable)));
1014 memcpy(kind->table, _tenDwiGageTable, sizeof(_tenDwiGageTable));
1015 kind->data = AIR_CAST(void *, tenDwiGageKindDataNew());
1016 }
1017 return kind;
1018 }
1019
1020 gageKind *
tenDwiGageKindNix(gageKind * kind)1021 tenDwiGageKindNix(gageKind *kind) {
1022
1023 if (kind) {
1024 airFree(kind->table);
1025 tenDwiGageKindDataNix(AIR_CAST(tenDwiGageKindData *, kind->data));
1026 airFree(kind);
1027 }
1028 return NULL;
1029 }
1030
1031 /*
1032 ** NOTE: this sets information in both the kind and kindData
1033 */
1034 int
tenDwiGageKindSet(gageKind * dwiKind,double thresh,double soft,double bval,double valueMin,const Nrrd * ngrad,const Nrrd * nbmat,int e1method,int e2method,unsigned int randSeed)1035 tenDwiGageKindSet(gageKind *dwiKind,
1036 double thresh, double soft, double bval, double valueMin,
1037 const Nrrd *ngrad,
1038 const Nrrd *nbmat,
1039 int e1method, int e2method,
1040 unsigned int randSeed) {
1041 static const char me[]="tenDwiGageKindSet";
1042 tenDwiGageKindData *kindData;
1043 double grad[3], (*lup)(const void *, size_t);
1044 unsigned int gi;
1045
1046 if (!dwiKind) {
1047 biffAddf(TEN, "%s: got NULL pointer", me);
1048 return 0;
1049 }
1050 if (!( !!(ngrad) ^ !!(nbmat) )) {
1051 biffAddf(TEN, "%s: need exactly one non-NULL in {ngrad,nbmat}", me);
1052 return 1;
1053 }
1054 if (nbmat) {
1055 biffAddf(TEN, "%s: sorry, B-matrices temporarily disabled", me);
1056 return 1;
1057 }
1058 /* (used for detecting errors in losslessly writing/reading
1059 a gradient set)
1060 {
1061 fprintf(stderr, "!%s: saving ngrad.nrrd\n", me);
1062 if (ngrad) {
1063 nrrdSave("ngrad.nrrd", ngrad, NULL);
1064 }
1065 }
1066 */
1067
1068 if (tenGradientCheck(ngrad, nrrdTypeDefault, 7)) {
1069 biffAddf(TEN, "%s: problem with given gradients", me);
1070 return 1;
1071 }
1072 /* make sure that gradient lengths are as expected */
1073 lup = nrrdDLookup[ngrad->type];
1074 grad[0] = lup(ngrad->data, 0);
1075 grad[1] = lup(ngrad->data, 1);
1076 grad[2] = lup(ngrad->data, 2);
1077 if (0.0 != ELL_3V_LEN(grad)) {
1078 biffAddf(TEN, "%s: sorry, currently need grad[0] to be len 0 (not %g)",
1079 me, ELL_3V_LEN(grad));
1080 return 1;
1081 }
1082 for (gi=1; gi<ngrad->axis[1].size; gi++) {
1083 grad[0] = lup(ngrad->data, 0 + 3*gi);
1084 grad[1] = lup(ngrad->data, 1 + 3*gi);
1085 grad[2] = lup(ngrad->data, 2 + 3*gi);
1086 if (0.0 == ELL_3V_LEN(grad)) {
1087 biffAddf(TEN, "%s: sorry, all but first gradient must be non-zero "
1088 "(%u is zero)", me, gi);
1089 return 1;
1090 }
1091 }
1092 if (airEnumValCheck(tenEstimate1Method, e1method)) {
1093 biffAddf(TEN, "%s: e1method %d is not a valid %s", me,
1094 e1method, tenEstimate1Method->name);
1095 return 1;
1096 }
1097 if (airEnumValCheck(tenEstimate2Method, e2method)) {
1098 biffAddf(TEN, "%s: emethod %d is not a valid %s", me,
1099 e2method, tenEstimate2Method->name);
1100 return 1;
1101 }
1102
1103 kindData = AIR_CAST(tenDwiGageKindData *, dwiKind->data);
1104 if (nrrdConvert(kindData->ngrad, ngrad, nrrdTypeDouble)) {
1105 biffMovef(TEN, NRRD, "%s: trouble converting", me);
1106 return 1;
1107 }
1108 dwiKind->valLen = kindData->ngrad->axis[1].size;
1109
1110 /* fixing up the item table ... */
1111 dwiKind->table[tenDwiGageAll].answerLength = dwiKind->valLen;
1112 dwiKind->table[tenDwiGageJustDWI].answerLength = dwiKind->valLen - 1;
1113 dwiKind->table[tenDwiGageADC].answerLength = dwiKind->valLen - 1;
1114 dwiKind->table[tenDwiGageTensorAllDWIError].answerLength =
1115 dwiKind->valLen - 1;
1116 switch (e1method) {
1117 case tenEstimate1MethodLLS:
1118 dwiKind->table[tenDwiGageTensor].prereq[0]
1119 = tenDwiGageTensorLLS;
1120 dwiKind->table[tenDwiGageTensorError].prereq[0]
1121 = tenDwiGageTensorLLSError;
1122 dwiKind->table[tenDwiGageTensorErrorLog].prereq[0]
1123 = tenDwiGageTensorLLSErrorLog;
1124 dwiKind->table[tenDwiGageTensorLikelihood].prereq[0]
1125 = tenDwiGageTensorLLSLikelihood;
1126 break;
1127 case tenEstimate1MethodWLS:
1128 dwiKind->table[tenDwiGageTensor].prereq[0]
1129 = tenDwiGageTensorWLS;
1130 dwiKind->table[tenDwiGageTensorError].prereq[0]
1131 = tenDwiGageTensorWLSError;
1132 dwiKind->table[tenDwiGageTensorErrorLog].prereq[0]
1133 = tenDwiGageTensorWLSErrorLog;
1134 dwiKind->table[tenDwiGageTensorLikelihood].prereq[0]
1135 = tenDwiGageTensorWLSLikelihood;
1136 break;
1137 case tenEstimate1MethodNLS:
1138 dwiKind->table[tenDwiGageTensor].prereq[0]
1139 = tenDwiGageTensorNLS;
1140 dwiKind->table[tenDwiGageTensorError].prereq[0]
1141 = tenDwiGageTensorNLSError;
1142 dwiKind->table[tenDwiGageTensorErrorLog].prereq[0]
1143 = tenDwiGageTensorNLSErrorLog;
1144 dwiKind->table[tenDwiGageTensorLikelihood].prereq[0]
1145 = tenDwiGageTensorNLSLikelihood;
1146 break;
1147 case tenEstimate1MethodMLE:
1148 dwiKind->table[tenDwiGageTensor].prereq[0]
1149 = tenDwiGageTensorMLE;
1150 dwiKind->table[tenDwiGageTensorError].prereq[0]
1151 = tenDwiGageTensorMLEError;
1152 dwiKind->table[tenDwiGageTensorErrorLog].prereq[0]
1153 = tenDwiGageTensorMLEErrorLog;
1154 dwiKind->table[tenDwiGageTensorLikelihood].prereq[0]
1155 = tenDwiGageTensorMLELikelihood;
1156 break;
1157 default:
1158 biffAddf(TEN, "%s: unimplemented %s: %s (%d)", me,
1159 tenEstimate1Method->name,
1160 airEnumStr(tenEstimate1Method, e1method), e1method);
1161 return 1;
1162 break;
1163 }
1164 kindData->thresh = thresh;
1165 kindData->soft = soft;
1166 kindData->bval = bval;
1167 kindData->valueMin = valueMin;
1168 kindData->est1Method = e1method;
1169 kindData->est2Method = e2method;
1170 kindData->randSeed = randSeed;
1171 return 0;
1172 }
1173
1174 int
tenDwiGageKindCheck(const gageKind * kind)1175 tenDwiGageKindCheck(const gageKind *kind) {
1176 static const char me[]="tenDwiGageKindCheck";
1177
1178 if (!kind) {
1179 biffAddf(TEN, "%s: got NULL pointer", me);
1180 return 1;
1181 }
1182 if (strcmp(kind->name, TEN_DWI_GAGE_KIND_NAME)) {
1183 biffAddf(TEN, "%s: got \"%s\" kind, not \"%s\"", me,
1184 kind->name, TEN_DWI_GAGE_KIND_NAME);
1185 return 1;
1186 }
1187 if (0 == kind->valLen) {
1188 biffAddf(TEN, "%s: don't yet know valLen", me);
1189 return 1;
1190 }
1191 if (!kind->data) {
1192 biffAddf(TEN, "%s: kind->data is NULL", me);
1193 return 1;
1194 }
1195 return 0;
1196 }
1197