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 |
|
|
#include "ten.h" |
25 |
|
|
#include "privateTen.h" |
26 |
|
|
|
27 |
|
|
/* |
28 |
|
|
** learned: when it looks like good-old LLS estimation is producing |
29 |
|
|
** nothing but zero tensors, see if your tec->valueMin is larger |
30 |
|
|
** than (what are problably) floating-point DWIs |
31 |
|
|
*/ |
32 |
|
|
|
33 |
|
|
/* |
34 |
|
|
|
35 |
|
|
http://www.mathworks.com/access/helpdesk/help/toolbox/curvefit/ch_fitt5.html#40515 |
36 |
|
|
|
37 |
|
|
*/ |
38 |
|
|
|
39 |
|
|
/* ---------------------------------------------- */ |
40 |
|
|
|
41 |
|
|
int |
42 |
|
|
_tenGaussian(double *retP, double m, double t, double s) { |
43 |
|
|
static const char me[]="_tenGaussian"; |
44 |
|
|
double diff, earg, den; |
45 |
|
|
|
46 |
|
|
if (!retP) { |
47 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
48 |
|
|
return 1; |
49 |
|
|
} |
50 |
|
|
diff = (m-t)/2; |
51 |
|
|
earg = -diff*diff/2; |
52 |
|
|
den = s*sqrt(2*AIR_PI); |
53 |
|
|
*retP = exp(earg)/den; |
54 |
|
|
if (!AIR_EXISTS(*retP)) { |
55 |
|
|
biffAddf(TEN, "%s: m=%g, t=%g, s=%g", me, m, t, s); |
56 |
|
|
biffAddf(TEN, "%s: diff=%g, earg=%g, den=%g", me, diff, earg, den); |
57 |
|
|
biffAddf(TEN, "%s: failed with ret = exp(%g)/%g = %g/%g = %g", |
58 |
|
|
me, earg, den, exp(earg), den, *retP); |
59 |
|
|
*retP = AIR_NAN; return 1; |
60 |
|
|
} |
61 |
|
|
return 0; |
62 |
|
|
} |
63 |
|
|
|
64 |
|
|
int |
65 |
|
|
_tenRicianTrue(double *retP, |
66 |
|
|
double m /* measured */, |
67 |
|
|
double t /* truth */, |
68 |
|
|
double s /* sigma */) { |
69 |
|
|
static const char me[]="_tenRicianTrue"; |
70 |
|
|
double mos, moss, mos2, tos2, tos, ss, earg, barg; |
71 |
|
|
|
72 |
|
|
if (!retP) { |
73 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
74 |
|
|
return 1; |
75 |
|
|
} |
76 |
|
|
|
77 |
|
|
mos = m/s; |
78 |
|
|
moss = mos/s; |
79 |
|
|
tos = t/s; |
80 |
|
|
ss = s*s; |
81 |
|
|
mos2 = mos*mos; |
82 |
|
|
tos2 = tos*tos; |
83 |
|
|
earg = -(mos2 + tos2)/2; |
84 |
|
|
barg = mos*tos; |
85 |
|
|
*retP = exp(earg)*airBesselI0(barg)*moss; |
86 |
|
|
|
87 |
|
|
if (!AIR_EXISTS(*retP)) { |
88 |
|
|
biffAddf(TEN, "%s: m=%g, t=%g, s=%g", me, m, t, s); |
89 |
|
|
biffAddf(TEN, "%s: mos=%g, moss=%g, tos=%g, ss=%g", |
90 |
|
|
me, mos, moss, tos, ss); |
91 |
|
|
biffAddf(TEN, "%s: mos2=%g, tos2=%g, earg=%g, barg=%g", me, |
92 |
|
|
mos2, tos2, earg, barg); |
93 |
|
|
biffAddf(TEN, "%s: failed: ret=exp(%g)*bessi0(%g)*%g = %g * %g * %g = %g", |
94 |
|
|
me, earg, barg, moss, exp(earg), airBesselI0(barg), moss, *retP); |
95 |
|
|
*retP = AIR_NAN; return 1; |
96 |
|
|
} |
97 |
|
|
return 0; |
98 |
|
|
} |
99 |
|
|
|
100 |
|
|
int |
101 |
|
|
_tenRicianSafe(double *retP, double m, double t, double s) { |
102 |
|
|
static const char me[]="_tenRicianSafe"; |
103 |
|
|
double diff, ric, gau, neer=10, faar=20; |
104 |
|
|
int E; |
105 |
|
|
|
106 |
|
|
if (!retP) { |
107 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
108 |
|
|
return 1; |
109 |
|
|
} |
110 |
|
|
|
111 |
|
|
diff = AIR_ABS(m-t)/s; |
112 |
|
|
E = 0; |
113 |
|
|
if (diff < neer) { |
114 |
|
|
if (!E) E |= _tenRicianTrue(retP, m, t, s); |
115 |
|
|
} else if (diff < faar) { |
116 |
|
|
if (!E) E |= _tenRicianTrue(&ric, m, t, s); |
117 |
|
|
if (!E) E |= _tenGaussian(&gau, m, t, s); |
118 |
|
|
if (!E) *retP = AIR_AFFINE(neer, diff, faar, ric, gau); |
119 |
|
|
} else { |
120 |
|
|
if (!E) E |= _tenGaussian(retP, m, t, s); |
121 |
|
|
} |
122 |
|
|
if (E) { |
123 |
|
|
biffAddf(TEN, "%s: failed with m=%g, t=%g, s=%g -> diff=%g", |
124 |
|
|
me, m, t, s, diff); |
125 |
|
|
*retP = AIR_NAN; return 1; |
126 |
|
|
} |
127 |
|
|
return 0; |
128 |
|
|
} |
129 |
|
|
|
130 |
|
|
int |
131 |
|
|
_tenRician(double *retP, |
132 |
|
|
double m /* measured */, |
133 |
|
|
double t /* truth */, |
134 |
|
|
double s /* sigma */) { |
135 |
|
|
static const char me[]="_tenRician"; |
136 |
|
|
double tos, ric, gau, loSignal=4.0, hiSignal=8.0; |
137 |
|
|
int E; |
138 |
|
|
|
139 |
|
|
if (!retP) { |
140 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
141 |
|
|
return 1; |
142 |
|
|
} |
143 |
|
|
if (!( m >= 0 && t >= 0 && s > 0 )) { |
144 |
|
|
biffAddf(TEN, "%s: got bad args: m=%g t=%g s=%g", me, m, t, s); |
145 |
|
|
*retP = AIR_NAN; return 1; |
146 |
|
|
} |
147 |
|
|
|
148 |
|
|
tos = t/s; |
149 |
|
|
E = 0; |
150 |
|
|
if (tos < loSignal) { |
151 |
|
|
if (!E) E |= _tenRicianSafe(retP, m, t, s); |
152 |
|
|
} else if (tos < hiSignal) { |
153 |
|
|
if (!E) E |= _tenRicianSafe(&ric, m, t, s); |
154 |
|
|
if (!E) E |= _tenGaussian(&gau, m, t, s); |
155 |
|
|
if (!E) *retP = AIR_AFFINE(loSignal, tos, hiSignal, ric, gau); |
156 |
|
|
} else { |
157 |
|
|
if (!E) E |= _tenGaussian(retP, m, t, s); |
158 |
|
|
} |
159 |
|
|
if (E) { |
160 |
|
|
biffAddf(TEN, "%s: failed with m=%g, t=%g, s=%g -> tos=%g", |
161 |
|
|
me, m, t, s, tos); |
162 |
|
|
*retP = AIR_NAN; return 1; |
163 |
|
|
} |
164 |
|
|
return 0; |
165 |
|
|
} |
166 |
|
|
|
167 |
|
|
enum { |
168 |
|
|
flagUnknown, |
169 |
|
|
flagEstimateMethod, |
170 |
|
|
flagBInfo, |
171 |
|
|
flagAllNum, |
172 |
|
|
flagDwiNum, |
173 |
|
|
flagAllAlloc, |
174 |
|
|
flagDwiAlloc, |
175 |
|
|
flagAllSet, |
176 |
|
|
flagDwiSet, |
177 |
|
|
flagSkipSet, |
178 |
|
|
flagWght, |
179 |
|
|
flagEmat, |
180 |
|
|
flagLast |
181 |
|
|
}; |
182 |
|
|
|
183 |
|
|
void |
184 |
|
|
_tenEstimateOutputInit(tenEstimateContext *tec) { |
185 |
|
|
|
186 |
|
59464 |
tec->estimatedB0 = AIR_NAN; |
187 |
|
29732 |
TEN_T_SET(tec->ten, AIR_NAN, |
188 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN, |
189 |
|
|
AIR_NAN, AIR_NAN, |
190 |
|
|
AIR_NAN); |
191 |
|
29732 |
tec->conf = AIR_NAN; |
192 |
|
29732 |
tec->mdwi = AIR_NAN; |
193 |
|
29732 |
tec->time = AIR_NAN; |
194 |
|
29732 |
tec->errorDwi = AIR_NAN; |
195 |
|
29732 |
tec->errorLogDwi = AIR_NAN; |
196 |
|
29732 |
tec->likelihoodDwi = AIR_NAN; |
197 |
|
29732 |
} |
198 |
|
|
|
199 |
|
|
tenEstimateContext * |
200 |
|
|
tenEstimateContextNew() { |
201 |
|
|
tenEstimateContext *tec; |
202 |
|
|
unsigned int fi; |
203 |
|
|
airPtrPtrUnion appu; |
204 |
|
|
|
205 |
|
|
|
206 |
|
64 |
tec = AIR_CAST(tenEstimateContext *, malloc(sizeof(tenEstimateContext))); |
207 |
✓✗ |
32 |
if (tec) { |
208 |
|
32 |
tec->bValue = AIR_NAN; |
209 |
|
32 |
tec->valueMin = AIR_NAN; |
210 |
|
32 |
tec->sigma = AIR_NAN; |
211 |
|
32 |
tec->dwiConfThresh = AIR_NAN; |
212 |
|
32 |
tec->dwiConfSoft = AIR_NAN; |
213 |
|
32 |
tec->_ngrad = NULL; |
214 |
|
32 |
tec->_nbmat = NULL; |
215 |
|
32 |
tec->skipList = NULL; |
216 |
|
32 |
appu.ui = &(tec->skipList); |
217 |
|
32 |
tec->skipListArr = airArrayNew(appu.v, NULL, |
218 |
|
|
2*sizeof(unsigned int), 128); |
219 |
|
32 |
tec->skipListArr->noReallocWhenSmaller = AIR_TRUE; |
220 |
|
32 |
tec->all_f = NULL; |
221 |
|
32 |
tec->all_d = NULL; |
222 |
|
32 |
tec->simulate = AIR_FALSE; |
223 |
|
32 |
tec->estimate1Method = tenEstimate1MethodUnknown; |
224 |
|
32 |
tec->estimateB0 = AIR_TRUE; |
225 |
|
32 |
tec->recordTime = AIR_FALSE; |
226 |
|
32 |
tec->recordErrorDwi = AIR_FALSE; |
227 |
|
32 |
tec->recordErrorLogDwi = AIR_FALSE; |
228 |
|
32 |
tec->recordLikelihoodDwi = AIR_FALSE; |
229 |
|
32 |
tec->verbose = 0; |
230 |
|
32 |
tec->progress = AIR_FALSE; |
231 |
|
32 |
tec->WLSIterNum = 3; |
232 |
✓✓ |
768 |
for (fi=flagUnknown+1; fi<flagLast; fi++) { |
233 |
|
352 |
tec->flag[fi] = AIR_FALSE; |
234 |
|
|
} |
235 |
|
32 |
tec->allNum = 0; |
236 |
|
32 |
tec->dwiNum = 0; |
237 |
|
32 |
tec->nbmat = nrrdNew(); |
238 |
|
32 |
tec->nwght = nrrdNew(); |
239 |
|
32 |
tec->nemat = nrrdNew(); |
240 |
|
32 |
tec->knownB0 = AIR_NAN; |
241 |
|
32 |
tec->all = NULL; |
242 |
|
32 |
tec->bnorm = NULL; |
243 |
|
32 |
tec->allTmp = NULL; |
244 |
|
32 |
tec->dwiTmp = NULL; |
245 |
|
32 |
tec->dwi = NULL; |
246 |
|
32 |
tec->skipLut = NULL; |
247 |
|
32 |
_tenEstimateOutputInit(tec); |
248 |
|
32 |
} |
249 |
|
32 |
return tec; |
250 |
|
|
} |
251 |
|
|
|
252 |
|
|
tenEstimateContext * |
253 |
|
|
tenEstimateContextNix(tenEstimateContext *tec) { |
254 |
|
|
|
255 |
✓✗ |
64 |
if (tec) { |
256 |
|
32 |
nrrdNuke(tec->nbmat); |
257 |
|
32 |
nrrdNuke(tec->nwght); |
258 |
|
32 |
nrrdNuke(tec->nemat); |
259 |
|
32 |
airArrayNuke(tec->skipListArr); |
260 |
|
32 |
airFree(tec->all); |
261 |
|
32 |
airFree(tec->bnorm); |
262 |
|
32 |
airFree(tec->allTmp); |
263 |
|
32 |
airFree(tec->dwiTmp); |
264 |
|
32 |
airFree(tec->dwi); |
265 |
|
32 |
airFree(tec->skipLut); |
266 |
|
32 |
airFree(tec); |
267 |
|
32 |
} |
268 |
|
32 |
return NULL; |
269 |
|
|
} |
270 |
|
|
|
271 |
|
|
void |
272 |
|
|
tenEstimateVerboseSet(tenEstimateContext *tec, |
273 |
|
|
int verbose) { |
274 |
✓✗ |
64 |
if (tec) { |
275 |
|
32 |
tec->verbose = verbose; |
276 |
|
32 |
} |
277 |
|
32 |
return; |
278 |
|
|
} |
279 |
|
|
|
280 |
|
|
void |
281 |
|
|
tenEstimateNegEvalShiftSet(tenEstimateContext *tec, int doit) { |
282 |
|
|
|
283 |
✓✗ |
64 |
if (tec) { |
284 |
|
32 |
tec->negEvalShift = !!doit; |
285 |
|
32 |
} |
286 |
|
32 |
return; |
287 |
|
|
} |
288 |
|
|
|
289 |
|
|
int |
290 |
|
|
tenEstimateMethodSet(tenEstimateContext *tec, int estimateMethod) { |
291 |
|
|
static const char me[]="tenEstimateMethodSet"; |
292 |
|
|
|
293 |
✗✓ |
64 |
if (!tec) { |
294 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
295 |
|
|
return 1; |
296 |
|
|
} |
297 |
✗✓ |
32 |
if (airEnumValCheck(tenEstimate1Method, estimateMethod)) { |
298 |
|
|
biffAddf(TEN, "%s: estimateMethod %d not a valid %s", me, |
299 |
|
|
estimateMethod, tenEstimate1Method->name); |
300 |
|
|
return 1; |
301 |
|
|
} |
302 |
|
|
|
303 |
✓✗ |
32 |
if (tec->estimate1Method != estimateMethod) { |
304 |
|
32 |
tec->estimate1Method = estimateMethod; |
305 |
|
32 |
tec->flag[flagEstimateMethod] = AIR_TRUE; |
306 |
|
32 |
} |
307 |
|
|
|
308 |
|
32 |
return 0; |
309 |
|
32 |
} |
310 |
|
|
|
311 |
|
|
int |
312 |
|
|
tenEstimateSigmaSet(tenEstimateContext *tec, double sigma) { |
313 |
|
|
static const char me[]="tenEstimateSigmaSet"; |
314 |
|
|
|
315 |
|
|
if (!tec) { |
316 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
317 |
|
|
return 1; |
318 |
|
|
} |
319 |
|
|
if (!(AIR_EXISTS(sigma) && sigma >= 0.0)) { |
320 |
|
|
biffAddf(TEN, "%s: given sigma (%g) not existent and >= 0.0", me, sigma); |
321 |
|
|
return 1; |
322 |
|
|
} |
323 |
|
|
|
324 |
|
|
tec->sigma = sigma; |
325 |
|
|
|
326 |
|
|
return 0; |
327 |
|
|
} |
328 |
|
|
|
329 |
|
|
int |
330 |
|
|
tenEstimateValueMinSet(tenEstimateContext *tec, double valueMin) { |
331 |
|
|
static const char me[]="tenEstimateValueMinSet"; |
332 |
|
|
|
333 |
✗✓ |
64 |
if (!tec) { |
334 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
335 |
|
|
return 1; |
336 |
|
|
} |
337 |
✗✓ |
32 |
if (!(AIR_EXISTS(valueMin) && valueMin > 0.0)) { |
338 |
|
|
biffAddf(TEN, "%s: given valueMin (%g) not existent and > 0.0", |
339 |
|
|
me, valueMin); |
340 |
|
|
return 1; |
341 |
|
|
} |
342 |
|
|
|
343 |
|
32 |
tec->valueMin = valueMin; |
344 |
|
|
|
345 |
|
32 |
return 0; |
346 |
|
32 |
} |
347 |
|
|
|
348 |
|
|
int |
349 |
|
|
tenEstimateGradientsSet(tenEstimateContext *tec, |
350 |
|
|
const Nrrd *ngrad, double bValue, int estimateB0) { |
351 |
|
|
static const char me[]="tenEstimateGradientsSet"; |
352 |
|
|
|
353 |
✗✓ |
64 |
if (!(tec && ngrad)) { |
354 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
355 |
|
|
return 1; |
356 |
|
|
} |
357 |
✗✓ |
32 |
if (!AIR_EXISTS(bValue)) { |
358 |
|
|
biffAddf(TEN, "%s: given b value doesn't exist", me); |
359 |
|
|
return 1; |
360 |
|
|
} |
361 |
✗✓ |
32 |
if (tenGradientCheck(ngrad, nrrdTypeDefault, 7)) { |
362 |
|
|
biffAddf(TEN, "%s: problem with gradient list", me); |
363 |
|
|
return 1; |
364 |
|
|
} |
365 |
|
|
|
366 |
|
32 |
tec->bValue = bValue; |
367 |
|
32 |
tec->_ngrad = ngrad; |
368 |
|
32 |
tec->_nbmat = NULL; |
369 |
|
32 |
tec->estimateB0 = estimateB0; |
370 |
|
|
|
371 |
|
32 |
tec->flag[flagBInfo] = AIR_TRUE; |
372 |
|
32 |
return 0; |
373 |
|
32 |
} |
374 |
|
|
|
375 |
|
|
int |
376 |
|
|
tenEstimateBMatricesSet(tenEstimateContext *tec, |
377 |
|
|
const Nrrd *nbmat, double bValue, int estimateB0) { |
378 |
|
|
static const char me[]="tenEstimateBMatricesSet"; |
379 |
|
|
|
380 |
|
|
if (!(tec && nbmat)) { |
381 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
382 |
|
|
return 1; |
383 |
|
|
} |
384 |
|
|
if (!AIR_EXISTS(bValue)) { |
385 |
|
|
biffAddf(TEN, "%s: given b value doesn't exist", me); |
386 |
|
|
return 1; |
387 |
|
|
} |
388 |
|
|
if (tenBMatrixCheck(nbmat, nrrdTypeDefault, 7)) { |
389 |
|
|
biffAddf(TEN, "%s: problem with b-matrix list", me); |
390 |
|
|
return 1; |
391 |
|
|
} |
392 |
|
|
|
393 |
|
|
tec->bValue = bValue; |
394 |
|
|
tec->_ngrad = NULL; |
395 |
|
|
tec->_nbmat = nbmat; |
396 |
|
|
tec->estimateB0 = estimateB0; |
397 |
|
|
|
398 |
|
|
tec->flag[flagBInfo] = AIR_TRUE; |
399 |
|
|
return 0; |
400 |
|
|
} |
401 |
|
|
|
402 |
|
|
int |
403 |
|
|
tenEstimateSkipSet(tenEstimateContext *tec, |
404 |
|
|
unsigned int valIdx, int doSkip) { |
405 |
|
|
static const char me[]="tenEstimateSkipSet"; |
406 |
|
|
unsigned int skipIdx; |
407 |
|
|
|
408 |
|
|
if (!tec) { |
409 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
410 |
|
|
return 1; |
411 |
|
|
} |
412 |
|
|
|
413 |
|
|
skipIdx = airArrayLenIncr(tec->skipListArr, 1); |
414 |
|
|
tec->skipList[0 + 2*skipIdx] = valIdx; |
415 |
|
|
tec->skipList[1 + 2*skipIdx] = !!doSkip; |
416 |
|
|
|
417 |
|
|
tec->flag[flagSkipSet] = AIR_TRUE; |
418 |
|
|
return 0; |
419 |
|
|
} |
420 |
|
|
|
421 |
|
|
int |
422 |
|
|
tenEstimateSkipReset(tenEstimateContext *tec) { |
423 |
|
|
static const char me[]="tenEstimateSkipReset"; |
424 |
|
|
|
425 |
|
|
if (!tec) { |
426 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
427 |
|
|
return 1; |
428 |
|
|
} |
429 |
|
|
|
430 |
|
|
airArrayLenSet(tec->skipListArr, 0); |
431 |
|
|
|
432 |
|
|
tec->flag[flagSkipSet] = AIR_TRUE; |
433 |
|
|
return 0; |
434 |
|
|
} |
435 |
|
|
|
436 |
|
|
int |
437 |
|
|
tenEstimateThresholdFind(double *threshP, unsigned char *isB0, Nrrd *nin4d) { |
438 |
|
|
static const char me[]="tenEstimateThresholdFind"; |
439 |
|
|
Nrrd **ndwi; |
440 |
|
|
airArray *mop; |
441 |
|
|
unsigned int slIdx, slNum, dwiAx, dwiNum, |
442 |
|
|
rangeAxisNum, rangeAxisIdx[NRRD_DIM_MAX]; |
443 |
|
|
int dwiIdx; |
444 |
|
|
|
445 |
|
|
mop = airMopNew(); |
446 |
|
|
|
447 |
|
|
if (!(threshP && isB0 && nin4d)) { |
448 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
449 |
|
|
airMopError(mop); return 1; |
450 |
|
|
} |
451 |
|
|
|
452 |
|
|
/* HEY: copied from tenEpiRegister4D() */ |
453 |
|
|
rangeAxisNum = nrrdRangeAxesGet(nin4d, rangeAxisIdx); |
454 |
|
|
if (0 == rangeAxisNum) { |
455 |
|
|
/* we fall back on old behavior */ |
456 |
|
|
dwiAx = 0; |
457 |
|
|
} else if (1 == rangeAxisNum) { |
458 |
|
|
/* thankfully there's exactly one range axis */ |
459 |
|
|
dwiAx = rangeAxisIdx[0]; |
460 |
|
|
} else { |
461 |
|
|
biffAddf(TEN, "%s: have %u range axes instead of 1, don't know which " |
462 |
|
|
"is DWI axis", me, rangeAxisNum); |
463 |
|
|
airMopError(mop); return 1; |
464 |
|
|
} |
465 |
|
|
|
466 |
|
|
slNum = nin4d->axis[dwiAx].size; |
467 |
|
|
dwiNum = 0; |
468 |
|
|
for (slIdx=0; slIdx<slNum; slIdx++) { |
469 |
|
|
dwiNum += !isB0[slIdx]; |
470 |
|
|
} |
471 |
|
|
if (0 == dwiNum) { |
472 |
|
|
biffAddf(TEN, "%s: somehow got zero DWIs", me); |
473 |
|
|
airMopError(mop); return 1; |
474 |
|
|
} |
475 |
|
|
ndwi = AIR_CALLOC(dwiNum, Nrrd *); |
476 |
|
|
airMopAdd(mop, ndwi, (airMopper)airFree, airMopAlways); |
477 |
|
|
dwiIdx = -1; |
478 |
|
|
for (slIdx=0; slIdx<slNum; slIdx++) { |
479 |
|
|
if (!isB0[slIdx]) { |
480 |
|
|
dwiIdx++; |
481 |
|
|
ndwi[dwiIdx] = nrrdNew(); |
482 |
|
|
airMopAdd(mop, ndwi[dwiIdx], (airMopper)nrrdNuke, airMopAlways); |
483 |
|
|
if (nrrdSlice(ndwi[dwiIdx], nin4d, dwiAx, slIdx)) { |
484 |
|
|
biffMovef(TEN, NRRD, |
485 |
|
|
"%s: trouble slicing DWI at index %u", me, slIdx); |
486 |
|
|
airMopError(mop); return 1; |
487 |
|
|
} |
488 |
|
|
} |
489 |
|
|
} |
490 |
|
|
if (_tenEpiRegThresholdFind(threshP, ndwi, dwiNum, AIR_FALSE, 1.5)) { |
491 |
|
|
biffAddf(TEN, "%s: trouble finding thresh", me); |
492 |
|
|
airMopError(mop); return 1; |
493 |
|
|
} |
494 |
|
|
|
495 |
|
|
airMopOkay(mop); |
496 |
|
|
return 0; |
497 |
|
|
} |
498 |
|
|
|
499 |
|
|
int |
500 |
|
|
tenEstimateThresholdSet(tenEstimateContext *tec, |
501 |
|
|
double thresh, double soft) { |
502 |
|
|
static const char me[]="tenEstimateThresholdSet"; |
503 |
|
|
|
504 |
✗✓ |
64 |
if (!tec) { |
505 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
506 |
|
|
return 1; |
507 |
|
|
} |
508 |
✓✗✗✓
|
64 |
if (!(AIR_EXISTS(thresh) && AIR_EXISTS(soft))) { |
509 |
|
|
biffAddf(TEN, "%s: not both threshold (%g) and softness (%g) exist", me, |
510 |
|
|
thresh, soft); |
511 |
|
|
return 1; |
512 |
|
|
} |
513 |
|
|
|
514 |
|
32 |
tec->dwiConfThresh = thresh; |
515 |
|
32 |
tec->dwiConfSoft = soft; |
516 |
|
|
|
517 |
|
32 |
return 0; |
518 |
|
32 |
} |
519 |
|
|
|
520 |
|
|
int |
521 |
|
|
_tenEstimateCheck(tenEstimateContext *tec) { |
522 |
|
|
static const char me[]="_tenEstimateCheck"; |
523 |
|
|
|
524 |
✗✓ |
64 |
if (!tec) { |
525 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
526 |
|
|
return 1; |
527 |
|
|
} |
528 |
✓✗✗✓
|
64 |
if (!( AIR_EXISTS(tec->valueMin) && tec->valueMin > 0.0)) { |
529 |
|
|
biffAddf(TEN, "%s: need a positive valueMin set (not %g)", |
530 |
|
|
me, tec->valueMin); |
531 |
|
|
return 1; |
532 |
|
|
} |
533 |
✓✗ |
32 |
if (!tec->simulate) { |
534 |
✗✓ |
32 |
if (!AIR_EXISTS(tec->bValue)) { |
535 |
|
|
biffAddf(TEN, "%s: b-value not set", me); |
536 |
|
|
return 1; |
537 |
|
|
} |
538 |
✗✓ |
32 |
if (airEnumValCheck(tenEstimate1Method, tec->estimate1Method)) { |
539 |
|
|
biffAddf(TEN, "%s: estimation method not set", me); |
540 |
|
|
return 1; |
541 |
|
|
} |
542 |
✗✗ |
32 |
if (tenEstimate1MethodMLE == tec->estimate1Method |
543 |
✗✓✗✗
|
32 |
&& !(AIR_EXISTS(tec->sigma) && (tec->sigma >= 0.0)) |
544 |
|
|
) { |
545 |
|
|
biffAddf(TEN, "%s: can't do %s estim w/out non-negative sigma set", me, |
546 |
|
|
airEnumStr(tenEstimate1Method, tenEstimate1MethodMLE)); |
547 |
|
|
return 1; |
548 |
|
|
} |
549 |
✓✗✗✓
|
64 |
if (!(AIR_EXISTS(tec->dwiConfThresh) && AIR_EXISTS(tec->dwiConfSoft))) { |
550 |
|
|
biffAddf(TEN, "%s: not both threshold (%g) and softness (%g) exist", me, |
551 |
|
|
tec->dwiConfThresh, tec->dwiConfSoft); |
552 |
|
|
return 1; |
553 |
|
|
} |
554 |
|
|
} |
555 |
✗✓✗✗
|
32 |
if (!( tec->_ngrad || tec->_nbmat )) { |
556 |
|
|
biffAddf(TEN, "%s: need to set either gradients or B-matrices", me); |
557 |
|
|
return 1; |
558 |
|
|
} |
559 |
|
|
|
560 |
|
32 |
return 0; |
561 |
|
32 |
} |
562 |
|
|
|
563 |
|
|
/* |
564 |
|
|
** allNum includes the skipped images |
565 |
|
|
** dwiNum does not include the skipped images |
566 |
|
|
*/ |
567 |
|
|
int |
568 |
|
|
_tenEstimateNumUpdate(tenEstimateContext *tec) { |
569 |
|
|
static const char me[]="_tenEstimateNumUpdate"; |
570 |
|
|
unsigned int newAllNum, newDwiNum, allIdx, |
571 |
|
|
skipListIdx, skipIdx, skipDo, skipNotNum; |
572 |
|
|
double (*lup)(const void *, size_t), gg[3], bb[6]; |
573 |
|
|
|
574 |
✗✗ |
64 |
if (tec->flag[flagBInfo] |
575 |
✗✓ |
32 |
|| tec->flag[flagSkipSet]) { |
576 |
✓✗ |
32 |
if (tec->_ngrad) { |
577 |
|
32 |
newAllNum = AIR_CAST(unsigned int, tec->_ngrad->axis[1].size); |
578 |
|
32 |
lup = nrrdDLookup[tec->_ngrad->type]; |
579 |
|
32 |
} else { |
580 |
|
|
newAllNum = AIR_CAST(unsigned int, tec->_nbmat->axis[1].size); |
581 |
|
|
lup = nrrdDLookup[tec->_nbmat->type]; |
582 |
|
|
} |
583 |
✓✗ |
32 |
if (tec->allNum != newAllNum) { |
584 |
|
32 |
tec->allNum = newAllNum; |
585 |
|
32 |
tec->flag[flagAllNum] = AIR_TRUE; |
586 |
|
32 |
} |
587 |
|
|
|
588 |
|
|
/* HEY: this should probably be its own update function, but its very |
589 |
|
|
convenient to allocate these allNum-length arrays here, immediately */ |
590 |
|
32 |
airFree(tec->skipLut); |
591 |
|
32 |
tec->skipLut = AIR_CAST(unsigned char *, calloc(tec->allNum, |
592 |
|
|
sizeof(unsigned char))); |
593 |
|
32 |
airFree(tec->bnorm); |
594 |
|
32 |
tec->bnorm = AIR_CAST(double *, calloc(tec->allNum, sizeof(double))); |
595 |
✓✗✗✓
|
64 |
if (!(tec->skipLut && tec->bnorm)) { |
596 |
|
|
biffAddf(TEN, "%s: couldn't allocate skipLut, bnorm vectors length %u\n", |
597 |
|
|
me, tec->allNum); |
598 |
|
|
return 1; |
599 |
|
|
} |
600 |
|
|
|
601 |
✗✓ |
64 |
for (skipListIdx=0; skipListIdx<tec->skipListArr->len; skipListIdx++) { |
602 |
|
|
skipIdx = tec->skipList[0 + 2*skipListIdx]; |
603 |
|
|
skipDo = tec->skipList[1 + 2*skipListIdx]; |
604 |
|
|
if (!(skipIdx < tec->allNum)) { |
605 |
|
|
biffAddf(TEN, "%s: skipList entry %u value index %u not < # vals %u", |
606 |
|
|
me, skipListIdx, skipIdx, tec->allNum); |
607 |
|
|
return 1; |
608 |
|
|
} |
609 |
|
|
tec->skipLut[skipIdx] = skipDo; |
610 |
|
|
} |
611 |
|
|
skipNotNum = 0; |
612 |
✓✓ |
768 |
for (skipIdx=0; skipIdx<tec->allNum; skipIdx++) { |
613 |
|
352 |
skipNotNum += !tec->skipLut[skipIdx]; |
614 |
|
|
} |
615 |
✗✓ |
32 |
if (!(skipNotNum >= 7 )) { |
616 |
|
|
biffAddf(TEN, "%s: number of not-skipped (used) values %u < minimum 7", |
617 |
|
|
me, skipNotNum); |
618 |
|
|
return 1; |
619 |
|
|
} |
620 |
|
|
|
621 |
|
|
newDwiNum = 0; |
622 |
✓✓ |
768 |
for (allIdx=0; allIdx<tec->allNum; allIdx++) { |
623 |
✗✓ |
352 |
if (tec->skipLut[allIdx]) { |
624 |
|
|
tec->bnorm[allIdx] = AIR_NAN; |
625 |
|
|
} else { |
626 |
✓✗ |
352 |
if (tec->_ngrad) { |
627 |
|
352 |
gg[0] = lup(tec->_ngrad->data, 0 + 3*allIdx); |
628 |
|
352 |
gg[1] = lup(tec->_ngrad->data, 1 + 3*allIdx); |
629 |
|
352 |
gg[2] = lup(tec->_ngrad->data, 2 + 3*allIdx); |
630 |
|
352 |
bb[0] = gg[0]*gg[0]; |
631 |
|
352 |
bb[1] = gg[1]*gg[0]; |
632 |
|
352 |
bb[2] = gg[2]*gg[0]; |
633 |
|
352 |
bb[3] = gg[1]*gg[1]; |
634 |
|
352 |
bb[4] = gg[2]*gg[1]; |
635 |
|
352 |
bb[5] = gg[2]*gg[2]; |
636 |
|
352 |
} else { |
637 |
|
|
bb[0] = lup(tec->_nbmat->data, 0 + 6*allIdx); |
638 |
|
|
bb[1] = lup(tec->_nbmat->data, 1 + 6*allIdx); |
639 |
|
|
bb[2] = lup(tec->_nbmat->data, 2 + 6*allIdx); |
640 |
|
|
bb[3] = lup(tec->_nbmat->data, 3 + 6*allIdx); |
641 |
|
|
bb[4] = lup(tec->_nbmat->data, 4 + 6*allIdx); |
642 |
|
|
bb[5] = lup(tec->_nbmat->data, 5 + 6*allIdx); |
643 |
|
|
} |
644 |
|
704 |
tec->bnorm[allIdx] = sqrt(bb[0]*bb[0] + 2*bb[1]*bb[1] + 2*bb[2]*bb[2] |
645 |
|
352 |
+ bb[3]*bb[3] + 2*bb[4]*bb[4] |
646 |
|
352 |
+ bb[5]*bb[5]); |
647 |
✗✓ |
352 |
if (tec->estimateB0) { |
648 |
|
|
++newDwiNum; |
649 |
|
|
} else { |
650 |
|
352 |
newDwiNum += (0.0 != tec->bnorm[allIdx]); |
651 |
|
|
} |
652 |
|
|
} |
653 |
|
|
} |
654 |
✓✗ |
32 |
if (tec->dwiNum != newDwiNum) { |
655 |
|
32 |
tec->dwiNum = newDwiNum; |
656 |
|
32 |
tec->flag[flagDwiNum] = AIR_TRUE; |
657 |
|
32 |
} |
658 |
✓✗✗✓
|
64 |
if (!tec->estimateB0 && (tec->allNum == tec->dwiNum)) { |
659 |
|
|
biffAddf(TEN, "%s: don't want to estimate B0, but all values are DW", me); |
660 |
|
|
return 1; |
661 |
|
|
} |
662 |
|
|
} |
663 |
|
32 |
return 0; |
664 |
|
32 |
} |
665 |
|
|
|
666 |
|
|
int |
667 |
|
|
_tenEstimateAllAllocUpdate(tenEstimateContext *tec) { |
668 |
|
|
static const char me[]="_tenEstimateAllAllocUpdate"; |
669 |
|
|
|
670 |
✓✗ |
64 |
if (tec->flag[flagAllNum]) { |
671 |
|
32 |
airFree(tec->all); |
672 |
|
32 |
airFree(tec->allTmp); |
673 |
|
32 |
tec->all = AIR_CAST(double *, calloc(tec->allNum, sizeof(double))); |
674 |
|
32 |
tec->allTmp = AIR_CAST(double *, calloc(tec->allNum, sizeof(double))); |
675 |
✓✗✗✓
|
64 |
if (!( tec->all && tec->allTmp )) { |
676 |
|
|
biffAddf(TEN, "%s: couldn't allocate \"all\" arrays (length %u)", me, |
677 |
|
|
tec->allNum); |
678 |
|
|
return 1; |
679 |
|
|
} |
680 |
|
32 |
tec->flag[flagAllAlloc] = AIR_TRUE; |
681 |
|
32 |
} |
682 |
|
32 |
return 0; |
683 |
|
32 |
} |
684 |
|
|
|
685 |
|
|
int |
686 |
|
|
_tenEstimateDwiAllocUpdate(tenEstimateContext *tec) { |
687 |
|
|
static const char me[]="_tenEstimateDwiAllocUpdate"; |
688 |
|
64 |
size_t size[2]; |
689 |
|
|
int E; |
690 |
|
|
|
691 |
✓✗ |
32 |
if (tec->flag[flagDwiNum]) { |
692 |
|
32 |
airFree(tec->dwi); |
693 |
|
32 |
airFree(tec->dwiTmp); |
694 |
|
32 |
tec->dwi = AIR_CAST(double *, calloc(tec->dwiNum, sizeof(double))); |
695 |
|
32 |
tec->dwiTmp = AIR_CAST(double *, calloc(tec->dwiNum, sizeof(double))); |
696 |
✓✗✗✓
|
64 |
if (!(tec->dwi && tec->dwiTmp)) { |
697 |
|
|
biffAddf(TEN, "%s: couldn't allocate DWI arrays (length %u)", me, |
698 |
|
|
tec->dwiNum); |
699 |
|
|
return 1; |
700 |
|
|
} |
701 |
|
|
E = 0; |
702 |
✓✗ |
64 |
if (!E) size[0] = (tec->estimateB0 ? 7 : 6); |
703 |
✓✗ |
64 |
if (!E) size[1] = tec->dwiNum; |
704 |
✓✗ |
64 |
if (!E) E |= nrrdMaybeAlloc_nva(tec->nbmat, nrrdTypeDouble, 2, size); |
705 |
✓✗ |
64 |
if (!E) size[0] = tec->dwiNum; |
706 |
✓✗ |
64 |
if (!E) size[1] = tec->dwiNum; |
707 |
✓✗ |
64 |
if (!E) E |= nrrdMaybeAlloc_nva(tec->nwght, nrrdTypeDouble, 2, size); |
708 |
✗✓ |
32 |
if (E) { |
709 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate dwi nrrds", me); |
710 |
|
|
return 1; |
711 |
|
|
} |
712 |
|
|
/* nrrdSave("0-nbmat.txt", tec->nbmat, NULL); */ |
713 |
|
32 |
tec->flag[flagDwiAlloc] = AIR_TRUE; |
714 |
|
32 |
} |
715 |
|
32 |
return 0; |
716 |
|
32 |
} |
717 |
|
|
|
718 |
|
|
int |
719 |
|
|
_tenEstimateAllSetUpdate(tenEstimateContext *tec) { |
720 |
|
|
/* static const char me[]="_tenEstimateAllSetUpdate"; */ |
721 |
|
|
/* unsigned int skipListIdx, skipIdx, skip, dwiIdx */; |
722 |
|
|
|
723 |
|
|
if (tec->flag[flagAllAlloc] |
724 |
|
|
|| tec->flag[flagDwiNum]) { |
725 |
|
|
|
726 |
|
|
} |
727 |
|
64 |
return 0; |
728 |
|
|
} |
729 |
|
|
|
730 |
|
|
int |
731 |
|
|
_tenEstimateDwiSetUpdate(tenEstimateContext *tec) { |
732 |
|
|
/* static const char me[]="_tenEstimateDwiSetUpdate"; */ |
733 |
|
|
double (*lup)(const void *, size_t I), gg[3], *bmat; |
734 |
|
|
unsigned int allIdx, dwiIdx, bmIdx; |
735 |
|
|
|
736 |
✗✗ |
64 |
if (tec->flag[flagAllNum] |
737 |
✗✓ |
32 |
|| tec->flag[flagDwiAlloc]) { |
738 |
✓✗ |
32 |
if (tec->_ngrad) { |
739 |
|
32 |
lup = nrrdDLookup[tec->_ngrad->type]; |
740 |
|
32 |
} else { |
741 |
|
|
lup = nrrdDLookup[tec->_nbmat->type]; |
742 |
|
|
} |
743 |
|
|
dwiIdx = 0; |
744 |
|
32 |
bmat = AIR_CAST(double*, tec->nbmat->data); |
745 |
✓✓ |
768 |
for (allIdx=0; allIdx<tec->allNum; allIdx++) { |
746 |
✓✓ |
704 |
if (!tec->skipLut[allIdx] |
747 |
✓✗✓✗
|
1056 |
&& (tec->estimateB0 || tec->bnorm[allIdx])) { |
748 |
✓✗ |
320 |
if (tec->_ngrad) { |
749 |
|
320 |
gg[0] = lup(tec->_ngrad->data, 0 + 3*allIdx); |
750 |
|
320 |
gg[1] = lup(tec->_ngrad->data, 1 + 3*allIdx); |
751 |
|
320 |
gg[2] = lup(tec->_ngrad->data, 2 + 3*allIdx); |
752 |
|
320 |
bmat[0] = gg[0]*gg[0]; |
753 |
|
320 |
bmat[1] = gg[1]*gg[0]; |
754 |
|
320 |
bmat[2] = gg[2]*gg[0]; |
755 |
|
320 |
bmat[3] = gg[1]*gg[1]; |
756 |
|
320 |
bmat[4] = gg[2]*gg[1]; |
757 |
|
320 |
bmat[5] = gg[2]*gg[2]; |
758 |
|
320 |
} else { |
759 |
|
|
for (bmIdx=0; bmIdx<6; bmIdx++) { |
760 |
|
|
bmat[bmIdx] = lup(tec->_nbmat->data, bmIdx + 6*allIdx); |
761 |
|
|
} |
762 |
|
|
} |
763 |
|
320 |
bmat[1] *= 2.0; |
764 |
|
320 |
bmat[2] *= 2.0; |
765 |
|
320 |
bmat[4] *= 2.0; |
766 |
✗✓ |
320 |
if (tec->estimateB0) { |
767 |
|
|
bmat[6] = -1; |
768 |
|
|
} |
769 |
|
320 |
bmat += tec->nbmat->axis[0].size; |
770 |
|
320 |
dwiIdx++; |
771 |
|
320 |
} |
772 |
|
|
} |
773 |
|
|
} |
774 |
|
32 |
return 0; |
775 |
|
|
} |
776 |
|
|
|
777 |
|
|
int |
778 |
|
|
_tenEstimateWghtUpdate(tenEstimateContext *tec) { |
779 |
|
|
/* static const char me[]="_tenEstimateWghtUpdate"; */ |
780 |
|
|
unsigned int dwiIdx; |
781 |
|
|
double *wght; |
782 |
|
|
|
783 |
|
64 |
wght = AIR_CAST(double *, tec->nwght->data); |
784 |
✗✗ |
32 |
if (tec->flag[flagDwiAlloc] |
785 |
✗✓ |
32 |
|| tec->flag[flagEstimateMethod]) { |
786 |
|
|
|
787 |
|
|
/* HEY: this is only useful for linear LS, no? */ |
788 |
✓✓ |
704 |
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
789 |
|
320 |
wght[dwiIdx + tec->dwiNum*dwiIdx] = 1.0; |
790 |
|
|
} |
791 |
|
|
|
792 |
|
32 |
tec->flag[flagEstimateMethod] = AIR_FALSE; |
793 |
|
32 |
tec->flag[flagWght] = AIR_TRUE; |
794 |
|
32 |
} |
795 |
|
32 |
return 0; |
796 |
|
|
} |
797 |
|
|
|
798 |
|
|
int |
799 |
|
|
_tenEstimateEmatUpdate(tenEstimateContext *tec) { |
800 |
|
|
static const char me[]="tenEstimateEmatUpdate"; |
801 |
|
|
|
802 |
✓✗ |
96 |
if (tec->flag[flagDwiSet] |
803 |
✓✗ |
64 |
|| tec->flag[flagWght]) { |
804 |
|
|
|
805 |
✓✗ |
32 |
if (!tec->simulate) { |
806 |
|
|
/* HEY: ignores weights! */ |
807 |
✗✓ |
32 |
if (ell_Nm_pseudo_inv(tec->nemat, tec->nbmat)) { |
808 |
|
|
biffMovef(TEN, ELL, "%s: trouble pseudo-inverting %ux%u B-matrix", me, |
809 |
|
|
AIR_CAST(unsigned int, tec->nbmat->axis[1].size), |
810 |
|
|
AIR_CAST(unsigned int, tec->nbmat->axis[0].size)); |
811 |
|
|
return 1; |
812 |
|
|
} |
813 |
|
|
} |
814 |
|
|
|
815 |
|
32 |
tec->flag[flagDwiSet] = AIR_FALSE; |
816 |
|
32 |
tec->flag[flagWght] = AIR_FALSE; |
817 |
|
32 |
} |
818 |
|
32 |
return 0; |
819 |
|
32 |
} |
820 |
|
|
|
821 |
|
|
int |
822 |
|
|
tenEstimateUpdate(tenEstimateContext *tec) { |
823 |
|
|
static const char me[]="tenEstimateUpdate"; |
824 |
|
|
int EE; |
825 |
|
|
|
826 |
|
|
EE = 0; |
827 |
✓✗ |
96 |
if (!EE) EE |= _tenEstimateCheck(tec); |
828 |
✓✗ |
64 |
if (!EE) EE |= _tenEstimateNumUpdate(tec); |
829 |
✓✗ |
64 |
if (!EE) EE |= _tenEstimateAllAllocUpdate(tec); |
830 |
✓✗ |
64 |
if (!EE) EE |= _tenEstimateDwiAllocUpdate(tec); |
831 |
✓✗ |
64 |
if (!EE) EE |= _tenEstimateAllSetUpdate(tec); |
832 |
✓✗ |
64 |
if (!EE) EE |= _tenEstimateDwiSetUpdate(tec); |
833 |
✓✗ |
64 |
if (!EE) EE |= _tenEstimateWghtUpdate(tec); |
834 |
✓✗ |
64 |
if (!EE) EE |= _tenEstimateEmatUpdate(tec); |
835 |
✗✓ |
32 |
if (EE) { |
836 |
|
|
biffAddf(TEN, "%s: problem updating", me); |
837 |
|
|
return 1; |
838 |
|
|
} |
839 |
|
32 |
return 0; |
840 |
|
32 |
} |
841 |
|
|
|
842 |
|
|
/* |
843 |
|
|
** from given tec->all_f or tec->all_d (whichever is non-NULL), sets: |
844 |
|
|
** tec->all[], |
845 |
|
|
** tec->dwi[] |
846 |
|
|
** tec->knownB0, if !tec->estimateB0, |
847 |
|
|
** tec->mdwi, |
848 |
|
|
** tec->conf (from tec->mdwi) |
849 |
|
|
*/ |
850 |
|
|
void |
851 |
|
|
_tenEstimateValuesSet(tenEstimateContext *tec) { |
852 |
|
|
unsigned int allIdx, dwiIdx, B0Num; |
853 |
|
|
double normSum; |
854 |
|
|
|
855 |
✓✗ |
59400 |
if (!tec->estimateB0) { |
856 |
|
29700 |
tec->knownB0 = 0; |
857 |
|
29700 |
} else { |
858 |
|
|
tec->knownB0 = AIR_NAN; |
859 |
|
|
} |
860 |
|
|
normSum = 0; |
861 |
|
29700 |
tec->mdwi = 0; |
862 |
|
|
B0Num = 0; |
863 |
|
|
dwiIdx = 0; |
864 |
✓✓ |
712800 |
for (allIdx=0; allIdx<tec->allNum; allIdx++) { |
865 |
✓✗ |
326700 |
if (!tec->skipLut[allIdx]) { |
866 |
✗✓ |
980100 |
tec->all[allIdx] = (tec->all_f |
867 |
|
|
? tec->all_f[allIdx] |
868 |
|
326700 |
: tec->all_d[allIdx]); |
869 |
|
326700 |
tec->mdwi += tec->bnorm[allIdx]*tec->all[allIdx]; |
870 |
|
326700 |
normSum += tec->bnorm[allIdx]; |
871 |
✓✗✓✓
|
653400 |
if (tec->estimateB0 || tec->bnorm[allIdx]) { |
872 |
|
297000 |
tec->dwi[dwiIdx++] = tec->all[allIdx]; |
873 |
|
297000 |
} else { |
874 |
|
29700 |
tec->knownB0 += tec->all[allIdx]; |
875 |
|
29700 |
B0Num += 1; |
876 |
|
|
} |
877 |
|
|
} |
878 |
|
|
} |
879 |
✓✗ |
29700 |
if (!tec->estimateB0) { |
880 |
|
29700 |
tec->knownB0 /= B0Num; |
881 |
|
29700 |
} |
882 |
|
29700 |
tec->mdwi /= normSum; |
883 |
✗✓✗✓
|
59400 |
if (tec->dwiConfSoft > 0) { |
884 |
|
29700 |
tec->conf = AIR_AFFINE(-1, airErf((tec->mdwi - tec->dwiConfThresh) |
885 |
|
|
/tec->dwiConfSoft), 1, |
886 |
|
|
0, 1); |
887 |
|
|
} else { |
888 |
|
29700 |
tec->conf = tec->mdwi > tec->dwiConfThresh; |
889 |
|
|
} |
890 |
|
|
return; |
891 |
|
29700 |
} |
892 |
|
|
|
893 |
|
|
/* |
894 |
|
|
** ASSUMES THAT dwiTmp[] has been stuff with all values simulated from model |
895 |
|
|
*/ |
896 |
|
|
double |
897 |
|
|
_tenEstimateErrorDwi(tenEstimateContext *tec) { |
898 |
|
|
unsigned int dwiIdx; |
899 |
|
|
double err, diff; |
900 |
|
|
|
901 |
|
|
err = 0; |
902 |
|
|
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
903 |
|
|
diff = tec->dwi[dwiIdx] - tec->dwiTmp[dwiIdx]; |
904 |
|
|
/* |
905 |
|
|
avg = (tec->dwi[dwiIdx] + tec->dwiTmp[dwiIdx])/2; |
906 |
|
|
avg = AIR_ABS(avg); |
907 |
|
|
if (avg) { |
908 |
|
|
err += diff*diff/(avg*avg); |
909 |
|
|
} |
910 |
|
|
*/ |
911 |
|
|
err += diff*diff; |
912 |
|
|
} |
913 |
|
|
err /= tec->dwiNum; |
914 |
|
|
return sqrt(err); |
915 |
|
|
} |
916 |
|
|
double |
917 |
|
|
_tenEstimateErrorLogDwi(tenEstimateContext *tec) { |
918 |
|
|
unsigned int dwiIdx; |
919 |
|
|
double err, diff; |
920 |
|
|
|
921 |
|
|
err = 0; |
922 |
|
|
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
923 |
|
|
diff = (log(AIR_MAX(tec->valueMin, tec->dwi[dwiIdx])) |
924 |
|
|
- log(AIR_MAX(tec->valueMin, tec->dwiTmp[dwiIdx]))); |
925 |
|
|
err += diff*diff; |
926 |
|
|
} |
927 |
|
|
err /= tec->dwiNum; |
928 |
|
|
return sqrt(err); |
929 |
|
|
} |
930 |
|
|
|
931 |
|
|
/* |
932 |
|
|
** sets: |
933 |
|
|
** tec->dwiTmp[] |
934 |
|
|
** and sets of all of them, regardless of estimateB0 |
935 |
|
|
*/ |
936 |
|
|
int |
937 |
|
|
_tenEstimate1TensorSimulateSingle(tenEstimateContext *tec, |
938 |
|
|
double sigma, double bValue, double B0, |
939 |
|
|
const double ten[7]) { |
940 |
|
|
static const char me[]="_tenEstimate1TensorSimulateSingle"; |
941 |
|
|
unsigned int dwiIdx, jj; |
942 |
|
|
double nr, ni, vv; |
943 |
|
|
const double *bmat; |
944 |
|
|
|
945 |
|
|
if (!( ten && ten )) { |
946 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
947 |
|
|
return 1; |
948 |
|
|
} |
949 |
|
|
if (!( AIR_EXISTS(sigma) && sigma >= 0 |
950 |
|
|
&& AIR_EXISTS(bValue) && AIR_EXISTS(B0) )) { |
951 |
|
|
biffAddf(TEN, "%s: got bad args: sigma %g, bValue %g, B0 %g\n", me, |
952 |
|
|
sigma, bValue, B0); |
953 |
|
|
return 1; |
954 |
|
|
} |
955 |
|
|
|
956 |
|
|
bmat = AIR_CAST(const double *, tec->nbmat->data); |
957 |
|
|
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
958 |
|
|
vv = 0; |
959 |
|
|
for (jj=0; jj<6; jj++) { |
960 |
|
|
vv += bmat[jj]*ten[1+jj]; |
961 |
|
|
} |
962 |
|
|
/* |
963 |
|
|
fprintf(stderr, "!%s: sigma = %g, bValue = %g, B0 = %g\n", me, |
964 |
|
|
sigma, bValue, B0); |
965 |
|
|
fprintf(stderr, "!%s[%u]: bmat=(%g %g %g %g %g %g)." |
966 |
|
|
"ten=(%g %g %g %g %g %g)\n", |
967 |
|
|
me, dwiIdx, |
968 |
|
|
bmat[0], bmat[1], bmat[2], bmat[3], bmat[4], bmat[5], |
969 |
|
|
ten[1], ten[2], ten[3], ten[4], ten[5], ten[6]); |
970 |
|
|
fprintf(stderr, "!%s: %g * exp(- %g * %g) = %g * exp(%g) = " |
971 |
|
|
"%g * %g = ... \n", me, |
972 |
|
|
B0, bValue, vv, B0, -bValue*vv, B0, exp(-bValue*vv)); |
973 |
|
|
*/ |
974 |
|
|
/* need AIR_MAX(0, vv) because B:D might be negative */ |
975 |
|
|
vv = B0*exp(-bValue*AIR_MAX(0, vv)); |
976 |
|
|
/* |
977 |
|
|
fprintf(stderr, "!%s: vv = %g\n", me, vv); |
978 |
|
|
*/ |
979 |
|
|
if (sigma > 0) { |
980 |
|
|
airNormalRand(&nr, &ni); |
981 |
|
|
nr *= sigma; |
982 |
|
|
ni *= sigma; |
983 |
|
|
vv = sqrt((vv+nr)*(vv+nr) + ni*ni); |
984 |
|
|
} |
985 |
|
|
tec->dwiTmp[dwiIdx] = vv; |
986 |
|
|
if (!AIR_EXISTS(tec->dwiTmp[dwiIdx])) { |
987 |
|
|
fprintf(stderr, "**********************************\n"); |
988 |
|
|
|
989 |
|
|
} |
990 |
|
|
/* |
991 |
|
|
if (tec->verbose) { |
992 |
|
|
fprintf(stderr, "%s: dwi[%u] = %g\n", me, dwiIdx, tec->dwiTmp[dwiIdx]); |
993 |
|
|
} |
994 |
|
|
*/ |
995 |
|
|
bmat += tec->nbmat->axis[0].size; |
996 |
|
|
} |
997 |
|
|
return 0; |
998 |
|
|
} |
999 |
|
|
|
1000 |
|
|
int |
1001 |
|
|
tenEstimate1TensorSimulateSingle_f(tenEstimateContext *tec, |
1002 |
|
|
float *simval, |
1003 |
|
|
float sigma, float bValue, float B0, |
1004 |
|
|
const float _ten[7]) { |
1005 |
|
|
static const char me[]="tenEstimate1TensorSimulateSingle_f"; |
1006 |
|
|
unsigned int allIdx, dwiIdx; |
1007 |
|
|
double ten[7]; |
1008 |
|
|
|
1009 |
|
|
if (!(tec && simval && _ten)) { |
1010 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1011 |
|
|
return 1; |
1012 |
|
|
} |
1013 |
|
|
|
1014 |
|
|
TEN_T_COPY(ten, _ten); |
1015 |
|
|
if (_tenEstimate1TensorSimulateSingle(tec, sigma, bValue, B0, ten)) { |
1016 |
|
|
biffAddf(TEN, "%s: ", me); |
1017 |
|
|
return 1; |
1018 |
|
|
} |
1019 |
|
|
dwiIdx = 0; |
1020 |
|
|
for (allIdx=0; allIdx<tec->allNum; allIdx++) { |
1021 |
|
|
if (tec->estimateB0 || tec->bnorm[allIdx]) { |
1022 |
|
|
simval[allIdx] = AIR_CAST(float, tec->dwiTmp[dwiIdx++]); |
1023 |
|
|
} else { |
1024 |
|
|
simval[allIdx] = B0; |
1025 |
|
|
} |
1026 |
|
|
} |
1027 |
|
|
return 0; |
1028 |
|
|
} |
1029 |
|
|
|
1030 |
|
|
int |
1031 |
|
|
tenEstimate1TensorSimulateSingle_d(tenEstimateContext *tec, |
1032 |
|
|
double *simval, |
1033 |
|
|
double sigma, double bValue, double B0, |
1034 |
|
|
const double ten[7]) { |
1035 |
|
|
static const char me[]="tenEstimate1TensorSimulateSingle_d"; |
1036 |
|
|
unsigned int allIdx, dwiIdx; |
1037 |
|
|
|
1038 |
|
|
if (!(tec && simval && ten)) { |
1039 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1040 |
|
|
return 1; |
1041 |
|
|
} |
1042 |
|
|
if (!( AIR_EXISTS(sigma) && sigma >= 0 |
1043 |
|
|
&& AIR_EXISTS(bValue) && AIR_EXISTS(B0) )) { |
1044 |
|
|
biffAddf(TEN, "%s: got bad bargs sigma %g, bValue %g, B0 %g\n", me, |
1045 |
|
|
sigma, bValue, B0); |
1046 |
|
|
return 1; |
1047 |
|
|
} |
1048 |
|
|
|
1049 |
|
|
if (_tenEstimate1TensorSimulateSingle(tec, sigma, bValue, B0, ten)) { |
1050 |
|
|
biffAddf(TEN, "%s: ", me); |
1051 |
|
|
return 1; |
1052 |
|
|
} |
1053 |
|
|
dwiIdx = 0; |
1054 |
|
|
for (allIdx=0; allIdx<tec->allNum; allIdx++) { |
1055 |
|
|
if (tec->estimateB0 || tec->bnorm[allIdx]) { |
1056 |
|
|
simval[allIdx] = tec->dwiTmp[dwiIdx++]; |
1057 |
|
|
} else { |
1058 |
|
|
simval[allIdx] = B0; |
1059 |
|
|
} |
1060 |
|
|
} |
1061 |
|
|
return 0; |
1062 |
|
|
} |
1063 |
|
|
|
1064 |
|
|
int |
1065 |
|
|
tenEstimate1TensorSimulateVolume(tenEstimateContext *tec, |
1066 |
|
|
Nrrd *ndwi, |
1067 |
|
|
double sigma, double bValue, |
1068 |
|
|
const Nrrd *nB0, const Nrrd *nten, |
1069 |
|
|
int outType, int keyValueSet) { |
1070 |
|
|
static const char me[]="tenEstimate1TensorSimulateVolume"; |
1071 |
|
|
size_t sizeTen, sizeX, sizeY, sizeZ, NN, II; |
1072 |
|
|
double (*tlup)(const void *, size_t), (*blup)(const void *, size_t), |
1073 |
|
|
(*lup)(const void *, size_t), ten_d[7], *dwi_d, B0; |
1074 |
|
|
float *dwi_f, ten_f[7]; |
1075 |
|
|
unsigned int tt, allIdx; |
1076 |
|
|
int axmap[4], E; |
1077 |
|
|
airArray *mop; |
1078 |
|
|
char stmp[3][AIR_STRLEN_SMALL]; |
1079 |
|
|
|
1080 |
|
|
if (!(tec && ndwi && nB0 && nten)) { |
1081 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1082 |
|
|
return 1; |
1083 |
|
|
} |
1084 |
|
|
/* this should have been done by update(), but why not */ |
1085 |
|
|
if (_tenEstimateCheck(tec)) { |
1086 |
|
|
biffAddf(TEN, "%s: problem in given context", me); |
1087 |
|
|
return 1; |
1088 |
|
|
} |
1089 |
|
|
if (!(AIR_EXISTS(sigma) && sigma >= 0.0 |
1090 |
|
|
&& AIR_EXISTS(bValue) && bValue >= 0.0)) { |
1091 |
|
|
biffAddf(TEN, "%s: got invalid sigma (%g) or bValue (%g)\n", me, |
1092 |
|
|
sigma, bValue); |
1093 |
|
|
return 1; |
1094 |
|
|
} |
1095 |
|
|
if (airEnumValCheck(nrrdType, outType)) { |
1096 |
|
|
biffAddf(TEN, "%s: requested output type %d not valid", me, outType); |
1097 |
|
|
return 1; |
1098 |
|
|
} |
1099 |
|
|
if (!( nrrdTypeFloat == outType || nrrdTypeDouble == outType )) { |
1100 |
|
|
biffAddf(TEN, "%s: requested output type (%s) not %s or %s", me, |
1101 |
|
|
airEnumStr(nrrdType, outType), |
1102 |
|
|
airEnumStr(nrrdType, nrrdTypeFloat), |
1103 |
|
|
airEnumStr(nrrdType, nrrdTypeDouble)); |
1104 |
|
|
return 1; |
1105 |
|
|
} |
1106 |
|
|
|
1107 |
|
|
mop = airMopNew(); |
1108 |
|
|
|
1109 |
|
|
sizeTen = nrrdKindSize(nrrdKind3DMaskedSymMatrix); |
1110 |
|
|
sizeX = nten->axis[1].size; |
1111 |
|
|
sizeY = nten->axis[2].size; |
1112 |
|
|
sizeZ = nten->axis[3].size; |
1113 |
|
|
if (!(3 == nB0->dim && |
1114 |
|
|
sizeX == nB0->axis[0].size && |
1115 |
|
|
sizeY == nB0->axis[1].size && |
1116 |
|
|
sizeZ == nB0->axis[2].size)) { |
1117 |
|
|
biffAddf(TEN, "%s: given B0 (%u-D) volume not 3-D %sx%sx%s", me, nB0->dim, |
1118 |
|
|
airSprintSize_t(stmp[0], sizeX), |
1119 |
|
|
airSprintSize_t(stmp[1], sizeY), |
1120 |
|
|
airSprintSize_t(stmp[2], sizeZ)); |
1121 |
|
|
return 1; |
1122 |
|
|
} |
1123 |
|
|
if (nrrdMaybeAlloc_va(ndwi, outType, 4, |
1124 |
|
|
AIR_CAST(size_t, tec->allNum), sizeX, sizeY, sizeZ)) { |
1125 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate DWI output", me); |
1126 |
|
|
airMopError(mop); return 1; |
1127 |
|
|
} |
1128 |
|
|
NN = sizeX * sizeY * sizeZ; |
1129 |
|
|
tlup = nrrdDLookup[nten->type]; |
1130 |
|
|
blup = nrrdDLookup[nB0->type]; |
1131 |
|
|
dwi_d = AIR_CAST(double *, ndwi->data); |
1132 |
|
|
dwi_f = AIR_CAST(float *, ndwi->data); |
1133 |
|
|
E = 0; |
1134 |
|
|
for (II=0; !E && II<NN; II++) { |
1135 |
|
|
B0 = blup(nB0->data, II); |
1136 |
|
|
if (nrrdTypeDouble == outType) { |
1137 |
|
|
for (tt=0; tt<7; tt++) { |
1138 |
|
|
ten_d[tt] = tlup(nten->data, tt + sizeTen*II); |
1139 |
|
|
} |
1140 |
|
|
E = tenEstimate1TensorSimulateSingle_d(tec, dwi_d, sigma, |
1141 |
|
|
bValue, B0, ten_d); |
1142 |
|
|
dwi_d += tec->allNum; |
1143 |
|
|
} else { |
1144 |
|
|
for (tt=0; tt<7; tt++) { |
1145 |
|
|
ten_f[tt] = AIR_CAST(float, tlup(nten->data, tt + sizeTen*II)); |
1146 |
|
|
} |
1147 |
|
|
E = tenEstimate1TensorSimulateSingle_f(tec, dwi_f, |
1148 |
|
|
AIR_CAST(float, sigma), |
1149 |
|
|
AIR_CAST(float, bValue), |
1150 |
|
|
AIR_CAST(float, B0), |
1151 |
|
|
ten_f); |
1152 |
|
|
dwi_f += tec->allNum; |
1153 |
|
|
} |
1154 |
|
|
if (E) { |
1155 |
|
|
biffAddf(TEN, "%s: failed at sample %s", me, |
1156 |
|
|
airSprintSize_t(stmp[0], II)); |
1157 |
|
|
airMopError(mop); return 1; |
1158 |
|
|
} |
1159 |
|
|
} |
1160 |
|
|
|
1161 |
|
|
ELL_4V_SET(axmap, -1, 1, 2, 3); |
1162 |
|
|
nrrdAxisInfoCopy(ndwi, nten, axmap, NRRD_AXIS_INFO_NONE); |
1163 |
|
|
ndwi->axis[0].kind = nrrdKindList; |
1164 |
|
|
if (nrrdBasicInfoCopy(ndwi, nten, |
1165 |
|
|
NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) { |
1166 |
|
|
biffMovef(TEN, NRRD, "%s:", me); |
1167 |
|
|
airMopError(mop); return 1; |
1168 |
|
|
} |
1169 |
|
|
if (keyValueSet) { |
1170 |
|
|
char keystr[AIR_STRLEN_MED], valstr[AIR_STRLEN_MED]; |
1171 |
|
|
|
1172 |
|
|
nrrdKeyValueAdd(ndwi, tenDWMRIModalityKey, tenDWMRIModalityVal); |
1173 |
|
|
sprintf(valstr, "%g", bValue); |
1174 |
|
|
nrrdKeyValueAdd(ndwi, tenDWMRIBValueKey, valstr); |
1175 |
|
|
if (tec->_ngrad) { |
1176 |
|
|
lup = nrrdDLookup[tec->_ngrad->type]; |
1177 |
|
|
for (allIdx=0; allIdx<tec->allNum; allIdx++) { |
1178 |
|
|
sprintf(keystr, tenDWMRIGradKeyFmt, allIdx); |
1179 |
|
|
sprintf(valstr, "%g %g %g", |
1180 |
|
|
lup(tec->_ngrad->data, 0 + 3*allIdx), |
1181 |
|
|
lup(tec->_ngrad->data, 1 + 3*allIdx), |
1182 |
|
|
lup(tec->_ngrad->data, 2 + 3*allIdx)); |
1183 |
|
|
nrrdKeyValueAdd(ndwi, keystr, valstr); |
1184 |
|
|
} |
1185 |
|
|
} else { |
1186 |
|
|
lup = nrrdDLookup[tec->_nbmat->type]; |
1187 |
|
|
for (allIdx=0; allIdx<tec->allNum; allIdx++) { |
1188 |
|
|
sprintf(keystr, tenDWMRIBmatKeyFmt, allIdx); |
1189 |
|
|
sprintf(valstr, "%g %g %g %g %g %g", |
1190 |
|
|
lup(tec->_nbmat->data, 0 + 6*allIdx), |
1191 |
|
|
lup(tec->_nbmat->data, 1 + 6*allIdx), |
1192 |
|
|
lup(tec->_nbmat->data, 2 + 6*allIdx), |
1193 |
|
|
lup(tec->_nbmat->data, 3 + 6*allIdx), |
1194 |
|
|
lup(tec->_nbmat->data, 4 + 6*allIdx), |
1195 |
|
|
lup(tec->_nbmat->data, 5 + 6*allIdx)); |
1196 |
|
|
nrrdKeyValueAdd(ndwi, keystr, valstr); |
1197 |
|
|
} |
1198 |
|
|
} |
1199 |
|
|
} |
1200 |
|
|
airMopOkay(mop); |
1201 |
|
|
return 0; |
1202 |
|
|
} |
1203 |
|
|
|
1204 |
|
|
/* |
1205 |
|
|
** sets: |
1206 |
|
|
** tec->ten[1..6] |
1207 |
|
|
** tec->B0, if tec->estimateB0 |
1208 |
|
|
*/ |
1209 |
|
|
int |
1210 |
|
|
_tenEstimate1Tensor_LLS(tenEstimateContext *tec) { |
1211 |
|
|
static const char me[]="_tenEstimate1Tensor_LLS"; |
1212 |
|
|
double *emat, tmp, logB0; |
1213 |
|
|
unsigned int ii, jj; |
1214 |
|
|
|
1215 |
|
59400 |
emat = AIR_CAST(double *, tec->nemat->data); |
1216 |
✗✓ |
29700 |
if (tec->verbose) { |
1217 |
|
|
fprintf(stderr, "!%s: estimateB0 = %d\n", me, tec->estimateB0); |
1218 |
|
|
} |
1219 |
✗✓ |
29700 |
if (tec->estimateB0) { |
1220 |
|
|
for (ii=0; ii<tec->allNum; ii++) { |
1221 |
|
|
tmp = AIR_MAX(tec->valueMin, tec->all[ii]); |
1222 |
|
|
tec->allTmp[ii] = -log(tmp)/(tec->bValue); |
1223 |
|
|
} |
1224 |
|
|
for (jj=0; jj<7; jj++) { |
1225 |
|
|
tmp = 0; |
1226 |
|
|
for (ii=0; ii<tec->allNum; ii++) { |
1227 |
|
|
tmp += emat[ii + tec->allNum*jj]*tec->allTmp[ii]; |
1228 |
|
|
} |
1229 |
|
|
if (jj < 6) { |
1230 |
|
|
tec->ten[1+jj] = tmp; |
1231 |
|
|
if (!AIR_EXISTS(tmp)) { |
1232 |
|
|
biffAddf(TEN, "%s: estimated non-existent tensor coef (%u) %g", |
1233 |
|
|
me, jj, tmp); |
1234 |
|
|
return 1; |
1235 |
|
|
} |
1236 |
|
|
} else { |
1237 |
|
|
/* we're on seventh row, for finding B0 */ |
1238 |
|
|
tec->estimatedB0 = exp(tec->bValue*tmp); |
1239 |
|
|
tec->estimatedB0 = AIR_MIN(FLT_MAX, tec->estimatedB0); |
1240 |
|
|
if (!AIR_EXISTS(tec->estimatedB0)) { |
1241 |
|
|
biffAddf(TEN, "%s: estimated non-existent B0 %g (b=%g, tmp=%g)", |
1242 |
|
|
me, tec->estimatedB0, tec->bValue, tmp); |
1243 |
|
|
return 1; |
1244 |
|
|
} |
1245 |
|
|
} |
1246 |
|
|
} |
1247 |
|
|
} else { |
1248 |
✓✓ |
89100 |
logB0 = log(AIR_MAX(tec->valueMin, tec->knownB0)); |
1249 |
✓✓ |
653400 |
for (ii=0; ii<tec->dwiNum; ii++) { |
1250 |
✓✓ |
891000 |
tmp = AIR_MAX(tec->valueMin, tec->dwi[ii]); |
1251 |
|
297000 |
tec->dwiTmp[ii] = (logB0 - log(tmp))/(tec->bValue); |
1252 |
|
|
} |
1253 |
✓✓ |
415800 |
for (jj=0; jj<6; jj++) { |
1254 |
|
|
tmp = 0; |
1255 |
✓✓ |
3920400 |
for (ii=0; ii<tec->dwiNum; ii++) { |
1256 |
|
1782000 |
tmp += emat[ii + tec->dwiNum*jj]*tec->dwiTmp[ii]; |
1257 |
✗✓ |
1782000 |
if (tec->verbose > 5) { |
1258 |
|
|
fprintf(stderr, "%s: emat[(%u,%u)=%u]*dwi[%u] = %g*%g --> %g\n", me, |
1259 |
|
|
ii, jj, ii + tec->dwiNum*jj, ii, |
1260 |
|
|
emat[ii + tec->dwiNum*jj], tec->dwiTmp[ii], |
1261 |
|
|
tmp); |
1262 |
|
|
} |
1263 |
|
|
} |
1264 |
|
178200 |
tec->ten[1+jj] = tmp; |
1265 |
|
|
} |
1266 |
|
|
} |
1267 |
|
29700 |
return 0; |
1268 |
|
29700 |
} |
1269 |
|
|
|
1270 |
|
|
int |
1271 |
|
|
_tenEstimate1Tensor_WLS(tenEstimateContext *tec) { |
1272 |
|
|
static const char me[]="_tenEstimate1Tensor_WLS"; |
1273 |
|
|
unsigned int dwiIdx, iter; |
1274 |
|
|
double *wght, dwi, sum; |
1275 |
|
|
|
1276 |
|
|
if (!tec) { |
1277 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1278 |
|
|
return 1; |
1279 |
|
|
} |
1280 |
|
|
|
1281 |
|
|
wght = AIR_CAST(double *, tec->nwght->data); |
1282 |
|
|
sum = 0; |
1283 |
|
|
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
1284 |
|
|
dwi = tec->dwi[dwiIdx]; |
1285 |
|
|
dwi = AIR_MAX(tec->valueMin, dwi); |
1286 |
|
|
sum += dwi*dwi; |
1287 |
|
|
} |
1288 |
|
|
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
1289 |
|
|
dwi = tec->dwi[dwiIdx]; |
1290 |
|
|
dwi = AIR_MAX(tec->valueMin, dwi); |
1291 |
|
|
wght[dwiIdx + tec->dwiNum*dwiIdx] = dwi*dwi/sum; |
1292 |
|
|
} |
1293 |
|
|
if (ell_Nm_wght_pseudo_inv(tec->nemat, tec->nbmat, tec->nwght)) { |
1294 |
|
|
biffMovef(TEN, ELL, |
1295 |
|
|
"%s(1): trouble wght-pseudo-inverting %ux%u B-matrix", me, |
1296 |
|
|
AIR_CAST(unsigned int, tec->nbmat->axis[1].size), |
1297 |
|
|
AIR_CAST(unsigned int, tec->nbmat->axis[0].size)); |
1298 |
|
|
return 1; |
1299 |
|
|
} |
1300 |
|
|
/* |
1301 |
|
|
nrrdSave("wght.txt", tec->nwght, NULL); |
1302 |
|
|
nrrdSave("bmat.txt", tec->nbmat, NULL); |
1303 |
|
|
nrrdSave("emat.txt", tec->nemat, NULL); |
1304 |
|
|
*/ |
1305 |
|
|
if (_tenEstimate1Tensor_LLS(tec)) { |
1306 |
|
|
biffAddf(TEN, "%s: initial weighted LLS failed", me); |
1307 |
|
|
return 1; |
1308 |
|
|
} |
1309 |
|
|
|
1310 |
|
|
for (iter=0; iter<tec->WLSIterNum; iter++) { |
1311 |
|
|
/* |
1312 |
|
|
fprintf(stderr, "!%s: bValue = %g, B0 = %g, ten = %g %g %g %g %g %g\n", |
1313 |
|
|
me, |
1314 |
|
|
tec->bValue, (tec->estimateB0 ? tec->estimatedB0 : tec->knownB0), |
1315 |
|
|
tec->ten[1], tec->ten[2], tec->ten[3], |
1316 |
|
|
tec->ten[4], tec->ten[5], tec->ten[6]); |
1317 |
|
|
*/ |
1318 |
|
|
if (_tenEstimate1TensorSimulateSingle(tec, 0.0, tec->bValue, |
1319 |
|
|
(tec->estimateB0 ? |
1320 |
|
|
tec->estimatedB0 |
1321 |
|
|
: tec->knownB0), tec->ten)) { |
1322 |
|
|
biffAddf(TEN, "%s: iter %u", me, iter); |
1323 |
|
|
return 1; |
1324 |
|
|
} |
1325 |
|
|
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
1326 |
|
|
dwi = tec->dwiTmp[dwiIdx]; |
1327 |
|
|
if (!AIR_EXISTS(dwi)) { |
1328 |
|
|
biffAddf(TEN, "%s: bad simulated dwi[%u] == %g (iter %u)", |
1329 |
|
|
me, dwiIdx, dwi, iter); |
1330 |
|
|
return 1; |
1331 |
|
|
} |
1332 |
|
|
wght[dwiIdx + tec->dwiNum*dwiIdx] = AIR_MAX(FLT_MIN, dwi*dwi); |
1333 |
|
|
} |
1334 |
|
|
if (ell_Nm_wght_pseudo_inv(tec->nemat, tec->nbmat, tec->nwght)) { |
1335 |
|
|
biffMovef(TEN, ELL, "%s(2): trouble w/ %ux%u B-matrix (iter %u)", me, |
1336 |
|
|
AIR_CAST(unsigned int, tec->nbmat->axis[1].size), |
1337 |
|
|
AIR_CAST(unsigned int, tec->nbmat->axis[0].size), iter); |
1338 |
|
|
return 1; |
1339 |
|
|
} |
1340 |
|
|
_tenEstimate1Tensor_LLS(tec); |
1341 |
|
|
} |
1342 |
|
|
|
1343 |
|
|
return 0; |
1344 |
|
|
} |
1345 |
|
|
|
1346 |
|
|
int |
1347 |
|
|
_tenEstimate1TensorGradient(tenEstimateContext *tec, |
1348 |
|
|
double *gradB0P, double gradTen[7], |
1349 |
|
|
double B0, double ten[7], |
1350 |
|
|
double epsilon, |
1351 |
|
|
int (*gradientCB)(tenEstimateContext *tec, |
1352 |
|
|
double *gradB0P, double gTen[7], |
1353 |
|
|
double B0, double ten[7]), |
1354 |
|
|
int (*badnessCB)(tenEstimateContext *tec, |
1355 |
|
|
double *badP, |
1356 |
|
|
double B0, double ten[7])) { |
1357 |
|
|
static const char me[]="_tenEstimate1TensorGradper"; |
1358 |
|
|
double forwTen[7], backTen[7], forwBad, backBad; |
1359 |
|
|
unsigned int ti; |
1360 |
|
|
|
1361 |
|
|
if (!( tec && gradB0P && gradTen && badnessCB && ten)) { |
1362 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1363 |
|
|
return 1; |
1364 |
|
|
} |
1365 |
|
|
|
1366 |
|
|
if (gradientCB) { |
1367 |
|
|
if (gradientCB(tec, gradB0P, gradTen, B0, ten)) { |
1368 |
|
|
biffAddf(TEN, "%s: problem with grad callback", me); |
1369 |
|
|
return 1; |
1370 |
|
|
} |
1371 |
|
|
} else { |
1372 |
|
|
/* we find gradient manually */ |
1373 |
|
|
gradTen[0] = 0; |
1374 |
|
|
for (ti=0; ti<6; ti++) { |
1375 |
|
|
TEN_T_COPY(forwTen, ten); |
1376 |
|
|
TEN_T_COPY(backTen, ten); |
1377 |
|
|
forwTen[ti+1] += epsilon; |
1378 |
|
|
backTen[ti+1] -= epsilon; |
1379 |
|
|
if (badnessCB(tec, &forwBad, B0, forwTen) |
1380 |
|
|
|| badnessCB(tec, &backBad, B0, backTen)) { |
1381 |
|
|
biffAddf(TEN, "%s: trouble at ti=%u", me, ti); |
1382 |
|
|
return 1; |
1383 |
|
|
} |
1384 |
|
|
gradTen[ti+1] = (forwBad - backBad)/(2*epsilon); |
1385 |
|
|
} |
1386 |
|
|
} |
1387 |
|
|
|
1388 |
|
|
return 0; |
1389 |
|
|
} |
1390 |
|
|
|
1391 |
|
|
int |
1392 |
|
|
_tenEstimate1TensorDescent(tenEstimateContext *tec, |
1393 |
|
|
int (*gradientCB)(tenEstimateContext *tec, |
1394 |
|
|
double *gradB0, |
1395 |
|
|
double gradTen[7], |
1396 |
|
|
double B0, |
1397 |
|
|
double ten[7]), |
1398 |
|
|
int (*badnessCB)(tenEstimateContext *tec, |
1399 |
|
|
double *badP, |
1400 |
|
|
double B0, |
1401 |
|
|
double ten[7])) { |
1402 |
|
|
static const char me[]="_tenEstimate1TensorDescent"; |
1403 |
|
|
double currB0, lastB0, currTen[7], lastTen[7], gradB0=AIR_NAN, gradTen[7], |
1404 |
|
|
epsilon, |
1405 |
|
|
stepSize, badInit, bad, badDelta, stepSizeMin = 0.00000000001, badLast; |
1406 |
|
|
unsigned int iter, iterMax = 100000; |
1407 |
|
|
|
1408 |
|
|
/* start with WLS fit since its probably close */ |
1409 |
|
|
_tenEstimate1Tensor_WLS(tec); |
1410 |
|
|
if (tec->verbose) { |
1411 |
|
|
fprintf(stderr, "%s: WLS gave %g %g %g %g %g %g\n", me, |
1412 |
|
|
tec->ten[1], tec->ten[2], tec->ten[3], |
1413 |
|
|
tec->ten[4], tec->ten[5], tec->ten[6]); |
1414 |
|
|
} |
1415 |
|
|
|
1416 |
|
|
if (badnessCB(tec, &badInit, |
1417 |
|
|
(tec->estimateB0 ? tec->estimatedB0 : tec->knownB0), tec->ten) |
1418 |
|
|
|| !AIR_EXISTS(badInit)) { |
1419 |
|
|
biffAddf(TEN, "%s: problem getting initial bad", me); |
1420 |
|
|
return 1; |
1421 |
|
|
} |
1422 |
|
|
if (tec->verbose) { |
1423 |
|
|
fprintf(stderr, "\n%s: ________________________________________\n", me); |
1424 |
|
|
fprintf(stderr, "%s: start: badInit = %g ---------------\n", me, badInit); |
1425 |
|
|
} |
1426 |
|
|
|
1427 |
|
|
epsilon = 0.0000001; |
1428 |
|
|
newepsilon: |
1429 |
|
|
if (_tenEstimate1TensorGradient(tec, &gradB0, gradTen, |
1430 |
|
|
(tec->estimateB0 |
1431 |
|
|
? tec->estimatedB0 |
1432 |
|
|
: tec->knownB0), |
1433 |
|
|
tec->ten, epsilon, |
1434 |
|
|
gradientCB, badnessCB)) { |
1435 |
|
|
biffAddf(TEN, "%s: problem getting initial gradient", me); |
1436 |
|
|
return 1; |
1437 |
|
|
} |
1438 |
|
|
if (!( AIR_EXISTS(gradB0) || 0 <= TEN_T_NORM(gradTen) )) { |
1439 |
|
|
biffAddf(TEN, "%s: got bad gradB0 %g or zero-norm tensor grad", |
1440 |
|
|
me, gradB0); |
1441 |
|
|
return 1; |
1442 |
|
|
} |
1443 |
|
|
if (tec->verbose) { |
1444 |
|
|
fprintf(stderr, "%s: gradTen (%s) = %g %g %g %g %g %g\n", me, |
1445 |
|
|
gradientCB ? "analytic" : "cent-diff", |
1446 |
|
|
gradTen[1], gradTen[2], gradTen[3], |
1447 |
|
|
gradTen[4], gradTen[5], gradTen[6]); |
1448 |
|
|
} |
1449 |
|
|
|
1450 |
|
|
stepSize = 0.1; |
1451 |
|
|
do { |
1452 |
|
|
stepSize /= 10; |
1453 |
|
|
TEN_T_SCALE_ADD2(currTen, 1.0, tec->ten, -stepSize, gradTen); |
1454 |
|
|
if (tec->estimateB0) { |
1455 |
|
|
currB0 = tec->estimatedB0 + -stepSize*gradB0; |
1456 |
|
|
} else { |
1457 |
|
|
currB0 = tec->knownB0; |
1458 |
|
|
} |
1459 |
|
|
if (badnessCB(tec, &bad, currB0, currTen) |
1460 |
|
|
|| !AIR_EXISTS(bad)) { |
1461 |
|
|
biffAddf(TEN, "%s: problem getting badness for stepSize", me); |
1462 |
|
|
return 1; |
1463 |
|
|
} |
1464 |
|
|
if (tec->verbose) { |
1465 |
|
|
fprintf(stderr, "%s: ************ stepSize = %g --> bad = %g\n", |
1466 |
|
|
me, stepSize, bad); |
1467 |
|
|
} |
1468 |
|
|
} while (bad > badInit && stepSize > stepSizeMin); |
1469 |
|
|
|
1470 |
|
|
if (stepSize <= stepSizeMin) { |
1471 |
|
|
if (epsilon > FLT_MIN) { |
1472 |
|
|
epsilon /= 10; |
1473 |
|
|
fprintf(stderr, "%s: re-trying initial step w/ eps %g\n", me, epsilon); |
1474 |
|
|
goto newepsilon; |
1475 |
|
|
} else { |
1476 |
|
|
biffAddf(TEN, "%s: never found a usable step size", me); |
1477 |
|
|
return 1; |
1478 |
|
|
} |
1479 |
|
|
} else if (tec->verbose) { |
1480 |
|
|
biffAddf(TEN, "%s: using step size %g\n", me, stepSize); |
1481 |
|
|
} |
1482 |
|
|
|
1483 |
|
|
iter = 0; |
1484 |
|
|
badLast = bad; |
1485 |
|
|
do { |
1486 |
|
|
iter++; |
1487 |
|
|
TEN_T_COPY(lastTen, currTen); |
1488 |
|
|
lastB0 = currB0; |
1489 |
|
|
if (0 == (iter % 3)) { |
1490 |
|
|
if (_tenEstimate1TensorGradient(tec, &gradB0, gradTen, |
1491 |
|
|
currB0, currTen, stepSize/5, |
1492 |
|
|
gradientCB, badnessCB) |
1493 |
|
|
|| !AIR_EXISTS(gradB0)) { |
1494 |
|
|
biffAddf(TEN, "%s[%u]: problem getting iter grad", me, iter); |
1495 |
|
|
return 1; |
1496 |
|
|
} |
1497 |
|
|
} |
1498 |
|
|
TEN_T_SCALE_INCR(currTen, -stepSize, gradTen); |
1499 |
|
|
if (tec->estimateB0) { |
1500 |
|
|
currB0 -= stepSize*gradB0; |
1501 |
|
|
} |
1502 |
|
|
if (badnessCB(tec, &bad, currB0, currTen) |
1503 |
|
|
|| !AIR_EXISTS(bad)) { |
1504 |
|
|
biffAddf(TEN, "%s[%u]: problem getting badness during grad", me, iter); |
1505 |
|
|
return 1; |
1506 |
|
|
} |
1507 |
|
|
if (tec->verbose) { |
1508 |
|
|
fprintf(stderr, "%s: %u bad = %g\n", me, iter, bad); |
1509 |
|
|
} |
1510 |
|
|
badDelta = bad - badLast; |
1511 |
|
|
badLast = bad; |
1512 |
|
|
if (badDelta > 0) { |
1513 |
|
|
stepSize /= 10; |
1514 |
|
|
if (tec->verbose) { |
1515 |
|
|
fprintf(stderr, "%s: badDelta %g > 0 ---> stepSize = %g\n", |
1516 |
|
|
me, badDelta, stepSize); |
1517 |
|
|
} |
1518 |
|
|
badDelta = -1; /* bogus improvement for loop continuation */ |
1519 |
|
|
TEN_T_COPY(currTen, lastTen); |
1520 |
|
|
currB0 = lastB0; |
1521 |
|
|
} |
1522 |
|
|
} while (iter < iterMax && (iter < 2 || badDelta < -0.00005)); |
1523 |
|
|
if (iter >= iterMax) { |
1524 |
|
|
biffAddf(TEN, "%s: didn't converge after %u iterations", me, iter); |
1525 |
|
|
return 1; |
1526 |
|
|
} |
1527 |
|
|
if (tec->verbose) { |
1528 |
|
|
fprintf(stderr, "%s: finished\n", me); |
1529 |
|
|
} |
1530 |
|
|
|
1531 |
|
|
ELL_6V_COPY(tec->ten+1, currTen+1); |
1532 |
|
|
tec->estimatedB0 = currB0; |
1533 |
|
|
|
1534 |
|
|
return 0; |
1535 |
|
|
} |
1536 |
|
|
|
1537 |
|
|
int |
1538 |
|
|
_tenEstimate1Tensor_GradientNLS(tenEstimateContext *tec, |
1539 |
|
|
double *gradB0P, double gradTen[7], |
1540 |
|
|
double currB0, double currTen[7]) { |
1541 |
|
|
static const char me[]="_tenEstimate1Tensor_GradientNLS"; |
1542 |
|
|
double *bmat, dot, tmp, diff, scl; |
1543 |
|
|
unsigned int dwiIdx; |
1544 |
|
|
|
1545 |
|
|
if (!(tec && gradB0P && gradTen && currTen)) { |
1546 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1547 |
|
|
return 1; |
1548 |
|
|
} |
1549 |
|
|
*gradB0P = 0; |
1550 |
|
|
TEN_T_SET(gradTen, 0, 0, 0, 0, 0, 0, 0); |
1551 |
|
|
bmat = AIR_CAST(double *, tec->nbmat->data); |
1552 |
|
|
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
1553 |
|
|
dot = ELL_6V_DOT(bmat, currTen+1); |
1554 |
|
|
tmp = currB0*exp(-(tec->bValue)*dot); |
1555 |
|
|
diff = tec->dwi[dwiIdx] - tmp; |
1556 |
|
|
scl = 2*diff*tmp*(tec->bValue); |
1557 |
|
|
ELL_6V_SCALE_INCR(gradTen+1, scl, bmat); |
1558 |
|
|
bmat += tec->nbmat->axis[0].size; |
1559 |
|
|
/* HEY: increment *gradB0P */ |
1560 |
|
|
} |
1561 |
|
|
ELL_6V_SCALE_INCR(gradTen+1, 1.0/tec->dwiNum, gradTen+1); |
1562 |
|
|
return 0; |
1563 |
|
|
} |
1564 |
|
|
|
1565 |
|
|
int |
1566 |
|
|
_tenEstimate1Tensor_BadnessNLS(tenEstimateContext *tec, |
1567 |
|
|
double *retP, |
1568 |
|
|
double currB0, double currTen[7]) { |
1569 |
|
|
static const char me[]="_tenEstimate1Tensor_BadnessNLS"; |
1570 |
|
|
|
1571 |
|
|
if (!(retP && tec)) { |
1572 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1573 |
|
|
return 1; |
1574 |
|
|
} |
1575 |
|
|
if (_tenEstimate1TensorSimulateSingle(tec, 0.0, tec->bValue, |
1576 |
|
|
currB0, currTen)) { |
1577 |
|
|
biffAddf(TEN, "%s: ", me); |
1578 |
|
|
return 1; |
1579 |
|
|
} |
1580 |
|
|
if (tec->verbose > 2) { |
1581 |
|
|
unsigned int di; |
1582 |
|
|
fprintf(stderr, "%s: simdwi =", me); |
1583 |
|
|
for (di=0; di<tec->dwiNum; di++) { |
1584 |
|
|
fprintf(stderr, " %g", tec->dwiTmp[di]); |
1585 |
|
|
} |
1586 |
|
|
fprintf(stderr, "\n"); |
1587 |
|
|
} |
1588 |
|
|
*retP = _tenEstimateErrorDwi(tec); |
1589 |
|
|
if (tec->verbose > 2) { |
1590 |
|
|
fprintf(stderr, "!%s: badness(%g, (%g) %g %g %g %g %g %g) = %g\n", |
1591 |
|
|
me, currB0, currTen[0], |
1592 |
|
|
currTen[1], currTen[2], currTen[3], |
1593 |
|
|
currTen[4], currTen[5], |
1594 |
|
|
currTen[6], *retP); |
1595 |
|
|
} |
1596 |
|
|
return 0; |
1597 |
|
|
} |
1598 |
|
|
|
1599 |
|
|
int |
1600 |
|
|
_tenEstimate1Tensor_NLS(tenEstimateContext *tec) { |
1601 |
|
|
static const char me[]="_tenEstimate1Tensor_NLS"; |
1602 |
|
|
|
1603 |
|
|
if (_tenEstimate1TensorDescent(tec, |
1604 |
|
|
NULL |
1605 |
|
|
/* _tenEstimate1Tensor_GradientNLS */ |
1606 |
|
|
, |
1607 |
|
|
_tenEstimate1Tensor_BadnessNLS)) { |
1608 |
|
|
biffAddf(TEN, "%s: ", me); |
1609 |
|
|
return 1; |
1610 |
|
|
} |
1611 |
|
|
return 0; |
1612 |
|
|
} |
1613 |
|
|
|
1614 |
|
|
int |
1615 |
|
|
_tenEstimate1Tensor_GradientMLE(tenEstimateContext *tec, |
1616 |
|
|
double *gradB0P, double gradTen[7], |
1617 |
|
|
double currB0, double currTen[7]) { |
1618 |
|
|
static const char me[]="_tenEstimate1Tensor_GradientMLE"; |
1619 |
|
|
double *bmat, dot, barg, tmp, scl, dwi, sigma, bval; |
1620 |
|
|
unsigned int dwiIdx; |
1621 |
|
|
|
1622 |
|
|
if (!(tec && gradB0P && gradTen && currTen)) { |
1623 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1624 |
|
|
return 1; |
1625 |
|
|
} |
1626 |
|
|
if (tec->verbose) { |
1627 |
|
|
fprintf(stderr, "%s grad (currTen = %g %g %g %g %g %g)\n", me, |
1628 |
|
|
currTen[1], currTen[2], currTen[3], |
1629 |
|
|
currTen[4], currTen[5], |
1630 |
|
|
currTen[6]); |
1631 |
|
|
} |
1632 |
|
|
TEN_T_SET(gradTen, 0, 0, 0, 0, 0, 0, 0); |
1633 |
|
|
*gradB0P = 0; |
1634 |
|
|
sigma = tec->sigma; |
1635 |
|
|
bval = tec->bValue; |
1636 |
|
|
bmat = AIR_CAST(double *, tec->nbmat->data); |
1637 |
|
|
for (dwiIdx=0; dwiIdx<tec->dwiNum; dwiIdx++) { |
1638 |
|
|
dwi = tec->dwi[dwiIdx]; |
1639 |
|
|
dot = ELL_6V_DOT(bmat, currTen+1); |
1640 |
|
|
barg = exp(-bval*dot)*(dwi/sigma)*(currB0/sigma); |
1641 |
|
|
tmp = (exp(bval*dot)/sigma)*dwi/airBesselI0(barg); |
1642 |
|
|
if (tec->verbose) { |
1643 |
|
|
fprintf(stderr, "%s[%u]: dot = %g, barg = %g, tmp = %g\n", me, dwiIdx, |
1644 |
|
|
dot, barg, tmp); |
1645 |
|
|
} |
1646 |
|
|
if (tmp > DBL_MIN) { |
1647 |
|
|
tmp = currB0/sigma - tmp*airBesselI1(barg); |
1648 |
|
|
} else { |
1649 |
|
|
tmp = currB0/sigma; |
1650 |
|
|
} |
1651 |
|
|
if (tec->verbose) { |
1652 |
|
|
fprintf(stderr, " ---- tmp = %g\n", tmp); |
1653 |
|
|
} |
1654 |
|
|
scl = tmp*exp(-2*bval*dot)*bval*currB0/sigma; |
1655 |
|
|
ELL_6V_SCALE_INCR(gradTen+1, scl, bmat); |
1656 |
|
|
if (tec->verbose) { |
1657 |
|
|
fprintf(stderr, "%s[%u]: bmat = %g %g %g %g %g %g\n", |
1658 |
|
|
me, dwiIdx, |
1659 |
|
|
bmat[0], bmat[1], bmat[2], |
1660 |
|
|
bmat[3], bmat[4], |
1661 |
|
|
bmat[5]); |
1662 |
|
|
fprintf(stderr, "%s[%u]: scl = %g -> gradTen = %g %g %g %g %g %g\n", |
1663 |
|
|
me, dwiIdx, scl, |
1664 |
|
|
gradTen[1], gradTen[2], gradTen[3], |
1665 |
|
|
gradTen[4], gradTen[5], |
1666 |
|
|
gradTen[6]); |
1667 |
|
|
} |
1668 |
|
|
if (!AIR_EXISTS(scl)) { |
1669 |
|
|
TEN_T_SET(gradTen, AIR_NAN, |
1670 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN, |
1671 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
1672 |
|
|
*gradB0P = AIR_NAN; |
1673 |
|
|
biffAddf(TEN, "%s: scl = %g, very sorry", me, scl); |
1674 |
|
|
return 1; |
1675 |
|
|
} |
1676 |
|
|
bmat += tec->nbmat->axis[0].size; |
1677 |
|
|
/* HEY: increment gradB0 */ |
1678 |
|
|
} |
1679 |
|
|
ELL_6V_SCALE_INCR(gradTen+1, 1.0/tec->dwiNum, gradTen+1); |
1680 |
|
|
if (tec->verbose) { |
1681 |
|
|
fprintf(stderr, "%s: final gradTen = %g %g %g %g %g %g\n", me, |
1682 |
|
|
gradTen[1], gradTen[2], gradTen[3], |
1683 |
|
|
gradTen[4], gradTen[5], |
1684 |
|
|
gradTen[6]); |
1685 |
|
|
} |
1686 |
|
|
return 0; |
1687 |
|
|
} |
1688 |
|
|
|
1689 |
|
|
int |
1690 |
|
|
_tenEstimate1Tensor_BadnessMLE(tenEstimateContext *tec, |
1691 |
|
|
double *retP, |
1692 |
|
|
double currB0, double curt[7]) { |
1693 |
|
|
static const char me[]="_tenEstimate1Tensor_BadnessMLE"; |
1694 |
|
|
unsigned int dwiIdx; |
1695 |
|
|
double *bmat, sum, rice, logrice=0, mesdwi=0, simdwi=0, dot=0; |
1696 |
|
|
int E; |
1697 |
|
|
|
1698 |
|
|
E = 0; |
1699 |
|
|
sum = 0; |
1700 |
|
|
bmat = AIR_CAST(double *, tec->nbmat->data); |
1701 |
|
|
for (dwiIdx=0; !E && dwiIdx<tec->dwiNum; dwiIdx++) { |
1702 |
|
|
dot = ELL_6V_DOT(bmat, curt+1); |
1703 |
|
|
simdwi = currB0*exp(-(tec->bValue)*dot); |
1704 |
|
|
mesdwi = tec->dwi[dwiIdx]; |
1705 |
|
|
if (!E) E |= _tenRician(&rice, mesdwi, simdwi, tec->sigma); |
1706 |
|
|
if (!E) E |= !AIR_EXISTS(rice); |
1707 |
|
|
if (!E) logrice = log(rice + DBL_MIN); |
1708 |
|
|
if (!E) sum += logrice; |
1709 |
|
|
if (!E) E |= !AIR_EXISTS(sum); |
1710 |
|
|
if (!E) bmat += tec->nbmat->axis[0].size; |
1711 |
|
|
} |
1712 |
|
|
if (E) { |
1713 |
|
|
biffAddf(TEN, "%s[%u]: dot = (%g %g %g %g %g %g).(%g %g %g %g %g %g) = %g", |
1714 |
|
|
me, dwiIdx, |
1715 |
|
|
bmat[0], bmat[1], bmat[2], bmat[3], bmat[4], bmat[5], |
1716 |
|
|
curt[1], curt[2], curt[3], curt[4], curt[5], curt[6], dot); |
1717 |
|
|
biffAddf(TEN, "%s[%u]: simdwi = %g * exp(-%g * %g) = %g * exp(%g) " |
1718 |
|
|
"= %g * %g = %g", me, dwiIdx, |
1719 |
|
|
currB0, tec->bValue, dot, |
1720 |
|
|
currB0, -(tec->bValue)*dot, |
1721 |
|
|
currB0, exp(-(tec->bValue)*dot), |
1722 |
|
|
currB0*exp(-(tec->bValue)*dot)); |
1723 |
|
|
biffAddf(TEN, "%s[%u]: mesdwi = %g, simdwi = %g, sigma = %g", me, dwiIdx, |
1724 |
|
|
mesdwi, simdwi, tec->sigma); |
1725 |
|
|
biffAddf(TEN, "%s[%u]: rice = %g, logrice = %g, sum = %g", me, dwiIdx, |
1726 |
|
|
rice, logrice, sum); |
1727 |
|
|
*retP = AIR_NAN; |
1728 |
|
|
return 1; |
1729 |
|
|
} |
1730 |
|
|
*retP = -sum/tec->dwiNum; |
1731 |
|
|
return 0; |
1732 |
|
|
} |
1733 |
|
|
|
1734 |
|
|
int |
1735 |
|
|
_tenEstimate1Tensor_MLE(tenEstimateContext *tec) { |
1736 |
|
|
static const char me[]="_tenEstimate1Tensor_MLE"; |
1737 |
|
|
|
1738 |
|
|
if (_tenEstimate1TensorDescent(tec, NULL, |
1739 |
|
|
_tenEstimate1Tensor_BadnessMLE)) { |
1740 |
|
|
biffAddf(TEN, "%s: ", me); |
1741 |
|
|
return 1; |
1742 |
|
|
} |
1743 |
|
|
|
1744 |
|
|
return 0; |
1745 |
|
|
} |
1746 |
|
|
|
1747 |
|
|
/* |
1748 |
|
|
** sets: |
1749 |
|
|
** tec->ten[0] (from tec->conf) |
1750 |
|
|
** tec->time, if tec->recordTime |
1751 |
|
|
** tec->errorDwi, if tec->recordErrorDwi |
1752 |
|
|
** tec->errorLogDwi, if tec->recordErrorLogDwi |
1753 |
|
|
** tec->likelihoodDwi, if tec->recordLikelihoodDwi |
1754 |
|
|
*/ |
1755 |
|
|
int |
1756 |
|
|
_tenEstimate1TensorSingle(tenEstimateContext *tec) { |
1757 |
|
|
static const char me[]="_tenEstimate1TensorSingle"; |
1758 |
|
|
double time0, B0; |
1759 |
|
|
int E; |
1760 |
|
|
|
1761 |
|
59400 |
_tenEstimateOutputInit(tec); |
1762 |
✗✓ |
59400 |
time0 = tec->recordTime ? airTime() : 0; |
1763 |
|
29700 |
_tenEstimateValuesSet(tec); |
1764 |
|
29700 |
tec->ten[0] = tec->conf; |
1765 |
✓✗✗✗
✗ |
29700 |
switch(tec->estimate1Method) { |
1766 |
|
|
case tenEstimate1MethodLLS: |
1767 |
|
29700 |
E = _tenEstimate1Tensor_LLS(tec); |
1768 |
|
29700 |
break; |
1769 |
|
|
case tenEstimate1MethodWLS: |
1770 |
|
|
E = _tenEstimate1Tensor_WLS(tec); |
1771 |
|
|
break; |
1772 |
|
|
case tenEstimate1MethodNLS: |
1773 |
|
|
E = _tenEstimate1Tensor_NLS(tec); |
1774 |
|
|
break; |
1775 |
|
|
case tenEstimate1MethodMLE: |
1776 |
|
|
E = _tenEstimate1Tensor_MLE(tec); |
1777 |
|
|
break; |
1778 |
|
|
default: |
1779 |
|
|
biffAddf(TEN, "%s: estimation method %d unimplemented", |
1780 |
|
|
me, tec->estimate1Method); |
1781 |
|
|
return 1; |
1782 |
|
|
} |
1783 |
✗✓ |
59400 |
tec->time = tec->recordTime ? airTime() - time0 : 0; |
1784 |
✗✓ |
29700 |
if (tec->negEvalShift) { |
1785 |
|
|
double eval[3]; |
1786 |
|
|
tenEigensolve_d(eval, NULL, tec->ten); |
1787 |
|
|
if (eval[2] < 0) { |
1788 |
|
|
tec->ten[1] += -eval[2]; |
1789 |
|
|
tec->ten[4] += -eval[2]; |
1790 |
|
|
tec->ten[6] += -eval[2]; |
1791 |
|
|
} |
1792 |
|
|
} |
1793 |
✗✓ |
29700 |
if (E) { |
1794 |
|
|
TEN_T_SET(tec->ten, AIR_NAN, |
1795 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN, |
1796 |
|
|
AIR_NAN, AIR_NAN, AIR_NAN); |
1797 |
|
|
if (tec->estimateB0) { |
1798 |
|
|
tec->estimatedB0 = AIR_NAN; |
1799 |
|
|
} |
1800 |
|
|
biffAddf(TEN, "%s: estimation failed", me); |
1801 |
|
|
return 1; |
1802 |
|
|
} |
1803 |
|
|
|
1804 |
✗✓ |
59400 |
if (tec->recordErrorDwi |
1805 |
✓✗ |
59400 |
|| tec->recordErrorLogDwi) { |
1806 |
|
|
B0 = tec->estimateB0 ? tec->estimatedB0 : tec->knownB0; |
1807 |
|
|
if (_tenEstimate1TensorSimulateSingle(tec, 0.0, tec->bValue, |
1808 |
|
|
B0, tec->ten)) { |
1809 |
|
|
biffAddf(TEN, "%s: simulation failed", me); |
1810 |
|
|
return 1; |
1811 |
|
|
} |
1812 |
|
|
if (tec->recordErrorDwi) { |
1813 |
|
|
tec->errorDwi = _tenEstimateErrorDwi(tec); |
1814 |
|
|
} |
1815 |
|
|
if (tec->recordErrorLogDwi) { |
1816 |
|
|
tec->errorLogDwi = _tenEstimateErrorLogDwi(tec); |
1817 |
|
|
} |
1818 |
|
|
} |
1819 |
|
|
|
1820 |
|
|
/* HEY: record likelihood! */ |
1821 |
|
|
|
1822 |
|
29700 |
return 0; |
1823 |
|
29700 |
} |
1824 |
|
|
|
1825 |
|
|
int |
1826 |
|
|
tenEstimate1TensorSingle_f(tenEstimateContext *tec, |
1827 |
|
|
float ten[7], const float *all) { |
1828 |
|
|
static const char me[]="tenEstimate1TensorSingle_f"; |
1829 |
|
|
|
1830 |
|
|
if (!(tec && ten && all)) { |
1831 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1832 |
|
|
return 1; |
1833 |
|
|
} |
1834 |
|
|
|
1835 |
|
|
tec->all_f = all; |
1836 |
|
|
tec->all_d = NULL; |
1837 |
|
|
/* |
1838 |
|
|
fprintf(stderr, "!%s(%u): B0 = %g,%g\n", me, __LINE__, |
1839 |
|
|
tec->knownB0, tec->estimatedB0); |
1840 |
|
|
*/ |
1841 |
|
|
if (_tenEstimate1TensorSingle(tec)) { |
1842 |
|
|
biffAddf(TEN, "%s: ", me); |
1843 |
|
|
return 1; |
1844 |
|
|
} |
1845 |
|
|
/* |
1846 |
|
|
fprintf(stderr, "!%s(%u): B0 = %g,%g\n", me, __LINE__, |
1847 |
|
|
tec->knownB0, tec->estimatedB0); |
1848 |
|
|
*/ |
1849 |
|
|
TEN_T_COPY_TT(ten, float, tec->ten); |
1850 |
|
|
|
1851 |
|
|
return 0; |
1852 |
|
|
} |
1853 |
|
|
|
1854 |
|
|
int |
1855 |
|
|
tenEstimate1TensorSingle_d(tenEstimateContext *tec, |
1856 |
|
|
double ten[7], const double *all) { |
1857 |
|
|
static const char me[]="tenEstimate1TensorSingle_d"; |
1858 |
|
|
unsigned int ii; |
1859 |
|
|
|
1860 |
✗✓ |
59400 |
if (!(tec && ten && all)) { |
1861 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1862 |
|
|
return 1; |
1863 |
|
|
} |
1864 |
|
|
|
1865 |
|
29700 |
tec->all_f = NULL; |
1866 |
|
29700 |
tec->all_d = all; |
1867 |
✗✓ |
29700 |
if (tec->verbose) { |
1868 |
|
|
for (ii=0; ii<tec->allNum; ii++) { |
1869 |
|
|
fprintf(stderr, "%s: dwi[%u] = %g\n", me, ii, |
1870 |
|
|
tec->all_d ? tec->all_d[ii] : tec->all_f[ii]); |
1871 |
|
|
} |
1872 |
|
|
fprintf(stderr, "%s: will estimate by %d (%s) \n" |
1873 |
|
|
" estimateB0 %d; valueMin %g\n", me, |
1874 |
|
|
tec->estimate1Method, |
1875 |
|
|
airEnumStr(tenEstimate1Method, tec->estimate1Method), |
1876 |
|
|
tec->estimateB0, tec->valueMin); |
1877 |
|
|
} |
1878 |
✗✓ |
29700 |
if (_tenEstimate1TensorSingle(tec)) { |
1879 |
|
|
biffAddf(TEN, "%s: ", me); |
1880 |
|
|
return 1; |
1881 |
|
|
} |
1882 |
✗✓ |
29700 |
if (tec->verbose) { |
1883 |
|
|
fprintf(stderr, "%s: ten = %g %g %g %g %g %g %g\n", me, |
1884 |
|
|
tec->ten[0], |
1885 |
|
|
tec->ten[1], tec->ten[2], tec->ten[3], |
1886 |
|
|
tec->ten[4], tec->ten[5], |
1887 |
|
|
tec->ten[6]); |
1888 |
|
|
} |
1889 |
|
29700 |
TEN_T_COPY(ten, tec->ten); |
1890 |
|
29700 |
return 0; |
1891 |
|
29700 |
} |
1892 |
|
|
|
1893 |
|
|
int |
1894 |
|
|
tenEstimate1TensorVolume4D(tenEstimateContext *tec, |
1895 |
|
|
Nrrd *nten, Nrrd **nB0P, Nrrd **nterrP, |
1896 |
|
|
const Nrrd *ndwi, int outType) { |
1897 |
|
|
static const char me[]="tenEstimate1TensorVolume4D"; |
1898 |
|
|
char doneStr[20]; |
1899 |
|
|
size_t sizeTen, sizeX, sizeY, sizeZ, NN, II, tick; |
1900 |
|
|
double *all, ten[7], (*lup)(const void *, size_t), |
1901 |
|
|
(*ins)(void *v, size_t I, double d); |
1902 |
|
|
unsigned int dd; |
1903 |
|
|
airArray *mop; |
1904 |
|
|
int axmap[4]; |
1905 |
|
|
char stmp[AIR_STRLEN_SMALL]; |
1906 |
|
|
|
1907 |
|
|
#if 0 |
1908 |
|
|
#define NUM 800 |
1909 |
|
|
double val[NUM], minVal=0, maxVal=10, arg; |
1910 |
|
|
unsigned int valIdx; |
1911 |
|
|
Nrrd *nval; |
1912 |
|
|
for (valIdx=0; valIdx<NUM; valIdx++) { |
1913 |
|
|
arg = AIR_AFFINE(0, valIdx, NUM-1, minVal, maxVal); |
1914 |
|
|
if (_tenRician(val + valIdx, arg, 1, 1)) { |
1915 |
|
|
biffAddf(TEN, "%s: you are out of luck", me); |
1916 |
|
|
return 1; |
1917 |
|
|
} |
1918 |
|
|
} |
1919 |
|
|
nval = nrrdNew(); |
1920 |
|
|
nrrdWrap(nval, val, nrrdTypeDouble, 1, AIR_CAST(size_t, NUM)); |
1921 |
|
|
nrrdSave("nval.nrrd", nval, NULL); |
1922 |
|
|
nrrdNix(nval); |
1923 |
|
|
#endif |
1924 |
|
|
|
1925 |
|
|
if (!(tec && nten && ndwi)) { |
1926 |
|
|
/* nerrP and _NB0P can be NULL */ |
1927 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1928 |
|
|
return 1; |
1929 |
|
|
} |
1930 |
|
|
if (nrrdCheck(ndwi)) { |
1931 |
|
|
biffMovef(TEN, NRRD, "%s: DWI volume not valid", me); |
1932 |
|
|
return 1; |
1933 |
|
|
} |
1934 |
|
|
if (!( 4 == ndwi->dim && 7 <= ndwi->axis[0].size )) { |
1935 |
|
|
biffAddf(TEN, "%s: DWI volume should be 4-D with axis 0 size >= 7", me); |
1936 |
|
|
return 1; |
1937 |
|
|
} |
1938 |
|
|
if (tec->allNum != ndwi->axis[0].size) { |
1939 |
|
|
biffAddf(TEN, "%s: from %s info, expected %u values per sample, " |
1940 |
|
|
"but have %s in volume", me, |
1941 |
|
|
tec->_ngrad ? "gradient" : "B-matrix", tec->allNum, |
1942 |
|
|
airSprintSize_t(stmp, ndwi->axis[0].size)); |
1943 |
|
|
return 1; |
1944 |
|
|
} |
1945 |
|
|
if (nrrdTypeBlock == ndwi->type) { |
1946 |
|
|
biffAddf(TEN, "%s: DWI volume has non-scalar type %s", me, |
1947 |
|
|
airEnumStr(nrrdType, ndwi->type)); |
1948 |
|
|
return 1; |
1949 |
|
|
} |
1950 |
|
|
if (airEnumValCheck(nrrdType, outType)) { |
1951 |
|
|
biffAddf(TEN, "%s: requested output type %d not valid", me, outType); |
1952 |
|
|
return 1; |
1953 |
|
|
} |
1954 |
|
|
if (!( nrrdTypeFloat == outType || nrrdTypeDouble == outType )) { |
1955 |
|
|
biffAddf(TEN, "%s: requested output type (%s) not %s or %s", me, |
1956 |
|
|
airEnumStr(nrrdType, outType), |
1957 |
|
|
airEnumStr(nrrdType, nrrdTypeFloat), |
1958 |
|
|
airEnumStr(nrrdType, nrrdTypeDouble)); |
1959 |
|
|
return 1; |
1960 |
|
|
} |
1961 |
|
|
if (nterrP) { |
1962 |
|
|
int recE, recEL, recLK; |
1963 |
|
|
recE = !!(tec->recordErrorDwi); |
1964 |
|
|
recEL = !!(tec->recordErrorLogDwi); |
1965 |
|
|
recLK = !!(tec->recordLikelihoodDwi); |
1966 |
|
|
if (1 != recE + recEL + recLK) { |
1967 |
|
|
biffAddf(TEN, "%s: requested error volume but need exactly one of " |
1968 |
|
|
"recordErrorDwi, recordErrorLogDwi, recordLikelihoodDwi " |
1969 |
|
|
"to be set", me); |
1970 |
|
|
return 1; |
1971 |
|
|
} |
1972 |
|
|
} |
1973 |
|
|
|
1974 |
|
|
mop = airMopNew(); |
1975 |
|
|
|
1976 |
|
|
sizeTen = nrrdKindSize(nrrdKind3DMaskedSymMatrix); |
1977 |
|
|
sizeX = ndwi->axis[1].size; |
1978 |
|
|
sizeY = ndwi->axis[2].size; |
1979 |
|
|
sizeZ = ndwi->axis[3].size; |
1980 |
|
|
all = AIR_CAST(double *, calloc(tec->allNum, sizeof(double))); |
1981 |
|
|
if (!all) { |
1982 |
|
|
biffAddf(TEN, "%s: couldn't allocate length %u array", me, tec->allNum); |
1983 |
|
|
airMopError(mop); return 1; |
1984 |
|
|
} |
1985 |
|
|
airMopAdd(mop, all, airFree, airMopAlways); |
1986 |
|
|
|
1987 |
|
|
if (nrrdMaybeAlloc_va(nten, outType, 4, |
1988 |
|
|
sizeTen, sizeX, sizeY, sizeZ)) { |
1989 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate tensor output", me); |
1990 |
|
|
airMopError(mop); return 1; |
1991 |
|
|
} |
1992 |
|
|
if (nB0P) { |
1993 |
|
|
*nB0P = nrrdNew(); |
1994 |
|
|
if (nrrdMaybeAlloc_va(*nB0P, outType, 3, |
1995 |
|
|
sizeX, sizeY, sizeZ)) { |
1996 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate B0 output", me); |
1997 |
|
|
airMopError(mop); return 1; |
1998 |
|
|
} |
1999 |
|
|
airMopAdd(mop, *nB0P, (airMopper)nrrdNuke, airMopOnError); |
2000 |
|
|
airMopAdd(mop, nB0P, (airMopper)airSetNull, airMopOnError); |
2001 |
|
|
} |
2002 |
|
|
if (nterrP) { |
2003 |
|
|
*nterrP = nrrdNew(); |
2004 |
|
|
if (nrrdMaybeAlloc_va(*nterrP, outType, 3, |
2005 |
|
|
sizeX, sizeY, sizeZ) |
2006 |
|
|
|| nrrdBasicInfoCopy(*nterrP, ndwi, |
2007 |
|
|
NRRD_BASIC_INFO_DATA_BIT |
2008 |
|
|
| NRRD_BASIC_INFO_TYPE_BIT |
2009 |
|
|
| NRRD_BASIC_INFO_BLOCKSIZE_BIT |
2010 |
|
|
| NRRD_BASIC_INFO_DIMENSION_BIT |
2011 |
|
|
| NRRD_BASIC_INFO_CONTENT_BIT |
2012 |
|
|
| NRRD_BASIC_INFO_MEASUREMENTFRAME_BIT |
2013 |
|
|
| NRRD_BASIC_INFO_COMMENTS_BIT |
2014 |
|
|
| (nrrdStateKeyValuePairsPropagate |
2015 |
|
|
? 0 |
2016 |
|
|
: NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) { |
2017 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't creatting fitting error output", me); |
2018 |
|
|
airMopError(mop); return 1; |
2019 |
|
|
} |
2020 |
|
|
ELL_3V_SET(axmap, 1, 2, 3); |
2021 |
|
|
nrrdAxisInfoCopy(*nterrP, ndwi, axmap, NRRD_AXIS_INFO_NONE); |
2022 |
|
|
airMopAdd(mop, *nterrP, (airMopper)nrrdNuke, airMopOnError); |
2023 |
|
|
airMopAdd(mop, nterrP, (airMopper)airSetNull, airMopOnError); |
2024 |
|
|
} |
2025 |
|
|
NN = sizeX * sizeY * sizeZ; |
2026 |
|
|
lup = nrrdDLookup[ndwi->type]; |
2027 |
|
|
ins = nrrdDInsert[outType]; |
2028 |
|
|
if (tec->progress) { |
2029 |
|
|
fprintf(stderr, "%s: ", me); |
2030 |
|
|
} |
2031 |
|
|
fflush(stderr); |
2032 |
|
|
tick = NN / 200; |
2033 |
|
|
tick = AIR_MAX(1, tick); |
2034 |
|
|
for (II=0; II<NN; II++) { |
2035 |
|
|
if (tec->progress && 0 == II%tick) { |
2036 |
|
|
fprintf(stderr, "%s", airDoneStr(0, II, NN-1, doneStr)); |
2037 |
|
|
} |
2038 |
|
|
for (dd=0; dd<tec->allNum; dd++) { |
2039 |
|
|
all[dd] = lup(ndwi->data, dd + tec->allNum*II); |
2040 |
|
|
} |
2041 |
|
|
/* |
2042 |
|
|
tec->verbose = 10*(II == 42509); |
2043 |
|
|
*/ |
2044 |
|
|
if (tec->verbose) { |
2045 |
|
|
fprintf(stderr, "!%s: hello; II=%u\n", me, AIR_CAST(unsigned int, II)); |
2046 |
|
|
} |
2047 |
|
|
if (tenEstimate1TensorSingle_d(tec, ten, all)) { |
2048 |
|
|
biffAddf(TEN, "%s: failed at sample %s", me, |
2049 |
|
|
airSprintSize_t(stmp, II)); |
2050 |
|
|
airMopError(mop); return 1; |
2051 |
|
|
} |
2052 |
|
|
ins(nten->data, 0 + sizeTen*II, ten[0]); |
2053 |
|
|
ins(nten->data, 1 + sizeTen*II, ten[1]); |
2054 |
|
|
ins(nten->data, 2 + sizeTen*II, ten[2]); |
2055 |
|
|
ins(nten->data, 3 + sizeTen*II, ten[3]); |
2056 |
|
|
ins(nten->data, 4 + sizeTen*II, ten[4]); |
2057 |
|
|
ins(nten->data, 5 + sizeTen*II, ten[5]); |
2058 |
|
|
ins(nten->data, 6 + sizeTen*II, ten[6]); |
2059 |
|
|
if (nB0P) { |
2060 |
|
|
ins((*nB0P)->data, II, (tec->estimateB0 |
2061 |
|
|
? tec->estimatedB0 |
2062 |
|
|
: tec->knownB0)); |
2063 |
|
|
} |
2064 |
|
|
if (nterrP) { |
2065 |
|
|
/* this works because above we checked that only one of the |
2066 |
|
|
tec->record* flags is set */ |
2067 |
|
|
if (tec->recordErrorDwi) { |
2068 |
|
|
ins((*nterrP)->data, II, tec->errorDwi); |
2069 |
|
|
} else if (tec->recordErrorLogDwi) { |
2070 |
|
|
ins((*nterrP)->data, II, tec->errorLogDwi); |
2071 |
|
|
} else if (tec->recordLikelihoodDwi) { |
2072 |
|
|
ins((*nterrP)->data, II, tec->likelihoodDwi); |
2073 |
|
|
} |
2074 |
|
|
} |
2075 |
|
|
} |
2076 |
|
|
if (tec->progress) { |
2077 |
|
|
fprintf(stderr, "%s\n", airDoneStr(0, II, NN-1, doneStr)); |
2078 |
|
|
} |
2079 |
|
|
|
2080 |
|
|
ELL_4V_SET(axmap, -1, 1, 2, 3); |
2081 |
|
|
nrrdAxisInfoCopy(nten, ndwi, axmap, NRRD_AXIS_INFO_NONE); |
2082 |
|
|
nten->axis[0].kind = nrrdKind3DMaskedSymMatrix; |
2083 |
|
|
if (nrrdBasicInfoCopy(nten, ndwi, |
2084 |
|
|
NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) { |
2085 |
|
|
biffAddf(NRRD, "%s:", me); |
2086 |
|
|
return 1; |
2087 |
|
|
} |
2088 |
|
|
|
2089 |
|
|
airMopOkay(mop); |
2090 |
|
|
return 0; |
2091 |
|
|
} |