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