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 |
|
|
int |
28 |
|
|
_tenEpiRegSave(const char *fname, Nrrd *nsingle, Nrrd **nmulti, |
29 |
|
|
int len, const char *desc) { |
30 |
|
|
static const char me[]="_tenEpiRegSave"; |
31 |
|
|
Nrrd *nout; |
32 |
|
|
airArray *mop; |
33 |
|
|
|
34 |
|
|
mop = airMopNew(); |
35 |
|
|
if (nsingle) { |
36 |
|
|
nout = nsingle; |
37 |
|
|
} else { |
38 |
|
|
airMopAdd(mop, nout=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
39 |
|
|
if (nrrdJoin(nout, (const Nrrd*const*)nmulti, len, 0, AIR_TRUE)) { |
40 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't join %s for output", me, desc); |
41 |
|
|
airMopError(mop); return 1; |
42 |
|
|
} |
43 |
|
|
} |
44 |
|
|
if (nrrdSave(fname, nout, NULL)) { |
45 |
|
|
biffMovef(TEN, NRRD, "%s: trouble saving %s to \"%s\"", me, desc, fname); |
46 |
|
|
airMopError(mop); return 1; |
47 |
|
|
} |
48 |
|
|
fprintf(stderr, "%s: saved %s to \"%s\"\n", me, desc, fname); |
49 |
|
|
airMopOkay(mop); |
50 |
|
|
return 0; |
51 |
|
|
} |
52 |
|
|
|
53 |
|
|
|
54 |
|
|
int |
55 |
|
|
_tenEpiRegCheck(Nrrd **nout, Nrrd **ndwi, unsigned int dwiLen, Nrrd *ngrad, |
56 |
|
|
int reference, |
57 |
|
|
double bwX, double bwY, double DWthr, |
58 |
|
|
const NrrdKernel *kern, double *kparm) { |
59 |
|
|
static const char me[]="_tenEpiRegCheck"; |
60 |
|
|
unsigned int ni; |
61 |
|
|
|
62 |
|
|
AIR_UNUSED(DWthr); |
63 |
|
|
if (!( nout && ndwi && ngrad && kern && kparm )) { |
64 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
65 |
|
|
return 1; |
66 |
|
|
} |
67 |
|
|
if (tenGradientCheck(ngrad, nrrdTypeDefault, 6 /* <-- HEY: not sure */)) { |
68 |
|
|
biffAddf(TEN, "%s: problem with given gradient list", me); |
69 |
|
|
return 1; |
70 |
|
|
} |
71 |
|
|
if (dwiLen != ngrad->axis[1].size) { |
72 |
|
|
char stmp[AIR_STRLEN_SMALL]; |
73 |
|
|
biffAddf(TEN, "%s: got %u DWIs, but %s gradient directions", me, |
74 |
|
|
dwiLen, airSprintSize_t(stmp, ngrad->axis[1].size)); |
75 |
|
|
return 1; |
76 |
|
|
} |
77 |
|
|
for (ni=0; ni<dwiLen; ni++) { |
78 |
|
|
if (!nout[ni]) { |
79 |
|
|
biffAddf(TEN, "%s: nout[%d] is NULL", me, ni); |
80 |
|
|
return 1; |
81 |
|
|
} |
82 |
|
|
if (nrrdCheck(ndwi[ni])) { |
83 |
|
|
biffMovef(TEN, NRRD, |
84 |
|
|
"%s: basic nrrd validity failed on ndwi[%d]", me, ni); |
85 |
|
|
return 1; |
86 |
|
|
} |
87 |
|
|
if (!nrrdSameSize(ndwi[0], ndwi[ni], AIR_TRUE)) { |
88 |
|
|
biffMovef(TEN, NRRD, "%s: ndwi[%d] is different from ndwi[0]", me, ni); |
89 |
|
|
return 1; |
90 |
|
|
} |
91 |
|
|
} |
92 |
|
|
if (!( 3 == ndwi[0]->dim )) { |
93 |
|
|
biffAddf(TEN, "%s: didn't get a set of 3-D arrays (got %d-D)", me, |
94 |
|
|
ndwi[0]->dim); |
95 |
|
|
return 1; |
96 |
|
|
} |
97 |
|
|
if (!( AIR_IN_CL(-1, reference, (int)dwiLen-1) )) { |
98 |
|
|
biffAddf(TEN, "%s: reference index %d not in valid range [-1,%d]", |
99 |
|
|
me, reference, dwiLen-1); |
100 |
|
|
return 1; |
101 |
|
|
} |
102 |
|
|
if (!( AIR_EXISTS(bwX) && AIR_EXISTS(bwY) )) { |
103 |
|
|
biffAddf(TEN, "%s: bwX, bwY don't both exist", me); |
104 |
|
|
return 1; |
105 |
|
|
} |
106 |
|
|
if (!( bwX >= 0 && bwY >= 0 )) { |
107 |
|
|
biffAddf(TEN, "%s: bwX (%g) and bwY (%g) are not both non-negative", |
108 |
|
|
me, bwX, bwY); |
109 |
|
|
return 1; |
110 |
|
|
} |
111 |
|
|
return 0; |
112 |
|
|
} |
113 |
|
|
|
114 |
|
|
/* |
115 |
|
|
** this assumes that all nblur[i] are valid nrrds, and does nothing |
116 |
|
|
** to manage them |
117 |
|
|
*/ |
118 |
|
|
int |
119 |
|
|
_tenEpiRegBlur(Nrrd **nblur, Nrrd **ndwi, unsigned int dwiLen, |
120 |
|
|
double bwX, double bwY, int verb) { |
121 |
|
|
static const char me[]="_tenEpiRegBlur"; |
122 |
|
|
NrrdResampleInfo *rinfo; |
123 |
|
|
airArray *mop; |
124 |
|
|
size_t ni, sx, sy, sz; |
125 |
|
|
double savemin[2], savemax[2]; |
126 |
|
|
|
127 |
|
|
if (!( bwX || bwY )) { |
128 |
|
|
if (verb) { |
129 |
|
|
fprintf(stderr, "%s:\n ", me); fflush(stderr); |
130 |
|
|
} |
131 |
|
|
for (ni=0; ni<dwiLen; ni++) { |
132 |
|
|
if (verb) { |
133 |
|
|
fprintf(stderr, "%2u ", (unsigned int)ni); fflush(stderr); |
134 |
|
|
} |
135 |
|
|
if (nrrdCopy(nblur[ni], ndwi[ni])) { |
136 |
|
|
biffMovef(TEN, NRRD, "%s: trouble copying ndwi[%u]", |
137 |
|
|
me, (unsigned int)ni); |
138 |
|
|
return 1; |
139 |
|
|
} |
140 |
|
|
} |
141 |
|
|
if (verb) { |
142 |
|
|
fprintf(stderr, "done\n"); |
143 |
|
|
} |
144 |
|
|
return 0; |
145 |
|
|
} |
146 |
|
|
/* else we need to blur */ |
147 |
|
|
sx = ndwi[0]->axis[0].size; |
148 |
|
|
sy = ndwi[0]->axis[1].size; |
149 |
|
|
sz = ndwi[0]->axis[2].size; |
150 |
|
|
mop = airMopNew(); |
151 |
|
|
rinfo = nrrdResampleInfoNew(); |
152 |
|
|
airMopAdd(mop, rinfo, (airMopper)nrrdResampleInfoNix, airMopAlways); |
153 |
|
|
if (bwX) { |
154 |
|
|
rinfo->kernel[0] = nrrdKernelGaussian; |
155 |
|
|
rinfo->parm[0][0] = bwX; |
156 |
|
|
rinfo->parm[0][1] = 3.0; /* how many stnd devs do we cut-off at */ |
157 |
|
|
} else { |
158 |
|
|
rinfo->kernel[0] = NULL; |
159 |
|
|
} |
160 |
|
|
if (bwY) { |
161 |
|
|
rinfo->kernel[1] = nrrdKernelGaussian; |
162 |
|
|
rinfo->parm[1][0] = bwY; |
163 |
|
|
rinfo->parm[1][1] = 3.0; /* how many stnd devs do we cut-off at */ |
164 |
|
|
} else { |
165 |
|
|
rinfo->kernel[1] = NULL; |
166 |
|
|
} |
167 |
|
|
rinfo->kernel[2] = NULL; |
168 |
|
|
ELL_3V_SET(rinfo->samples, sx, sy, sz); |
169 |
|
|
ELL_3V_SET(rinfo->min, 0, 0, 0); |
170 |
|
|
ELL_3V_SET(rinfo->max, sx-1, sy-1, sz-1); |
171 |
|
|
rinfo->boundary = nrrdBoundaryBleed; |
172 |
|
|
rinfo->type = nrrdTypeDefault; |
173 |
|
|
rinfo->renormalize = AIR_TRUE; |
174 |
|
|
rinfo->clamp = AIR_TRUE; |
175 |
|
|
if (verb) { |
176 |
|
|
fprintf(stderr, "%s:\n ", me); fflush(stderr); |
177 |
|
|
} |
178 |
|
|
for (ni=0; ni<dwiLen; ni++) { |
179 |
|
|
if (verb) { |
180 |
|
|
fprintf(stderr, "%2u ", (unsigned int)ni); fflush(stderr); |
181 |
|
|
} |
182 |
|
|
savemin[0] = ndwi[ni]->axis[0].min; savemax[0] = ndwi[ni]->axis[0].max; |
183 |
|
|
savemin[1] = ndwi[ni]->axis[1].min; savemax[1] = ndwi[ni]->axis[1].max; |
184 |
|
|
ndwi[ni]->axis[0].min = 0; ndwi[ni]->axis[0].max = sx-1; |
185 |
|
|
ndwi[ni]->axis[1].min = 0; ndwi[ni]->axis[1].max = sy-1; |
186 |
|
|
if (nrrdSpatialResample(nblur[ni], ndwi[ni], rinfo)) { |
187 |
|
|
biffMovef(TEN, NRRD, "%s: trouble blurring ndwi[%u]", |
188 |
|
|
me, (unsigned int)ni); |
189 |
|
|
airMopError(mop); return 1; |
190 |
|
|
} |
191 |
|
|
ndwi[ni]->axis[0].min = savemin[0]; ndwi[ni]->axis[0].max = savemax[0]; |
192 |
|
|
ndwi[ni]->axis[1].min = savemin[1]; ndwi[ni]->axis[1].max = savemax[1]; |
193 |
|
|
} |
194 |
|
|
if (verb) { |
195 |
|
|
fprintf(stderr, "done\n"); |
196 |
|
|
} |
197 |
|
|
airMopOkay(mop); |
198 |
|
|
return 0; |
199 |
|
|
} |
200 |
|
|
|
201 |
|
|
int |
202 |
|
|
_tenEpiRegThresholdFind(double *DWthrP, Nrrd **nin, int ninLen, |
203 |
|
|
int save, double expo) { |
204 |
|
|
static const char me[]="_tenEpiRegThresholdFind"; |
205 |
|
|
Nrrd *nhist, *ntmp; |
206 |
|
|
airArray *mop; |
207 |
|
|
int ni, bins, E; |
208 |
|
|
double min=0, max=0; |
209 |
|
|
NrrdRange *range; |
210 |
|
|
|
211 |
|
|
mop = airMopNew(); |
212 |
|
|
airMopAdd(mop, nhist=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
213 |
|
|
airMopAdd(mop, ntmp=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
214 |
|
|
|
215 |
|
|
for (ni=0; ni<ninLen; ni++) { |
216 |
|
|
range = nrrdRangeNewSet(nin[ni], nrrdBlind8BitRangeFalse); |
217 |
|
|
if (!ni) { |
218 |
|
|
min = range->min; |
219 |
|
|
max = range->max; |
220 |
|
|
} else { |
221 |
|
|
min = AIR_MIN(min, range->min); |
222 |
|
|
max = AIR_MAX(max, range->max); |
223 |
|
|
} |
224 |
|
|
range = nrrdRangeNix(range); |
225 |
|
|
} |
226 |
|
|
bins = AIR_MIN(1024, (int)(max - min + 1)); |
227 |
|
|
ntmp->axis[0].min = min; |
228 |
|
|
ntmp->axis[0].max = max; |
229 |
|
|
for (ni=0; ni<ninLen; ni++) { |
230 |
|
|
if (nrrdHisto(ntmp, nin[ni], NULL, NULL, bins, nrrdTypeFloat)) { |
231 |
|
|
biffMovef(TEN, NRRD, |
232 |
|
|
"%s: problem forming histogram of DWI %d", me, ni); |
233 |
|
|
airMopError(mop); return 1; |
234 |
|
|
} |
235 |
|
|
if (!ni) { |
236 |
|
|
E = nrrdCopy(nhist, ntmp); |
237 |
|
|
} else { |
238 |
|
|
E = nrrdArithBinaryOp(nhist, nrrdBinaryOpAdd, nhist, ntmp); |
239 |
|
|
} |
240 |
|
|
if (E) { |
241 |
|
|
biffMovef(TEN, NRRD, |
242 |
|
|
"%s: problem updating histogram sum on DWI %d", |
243 |
|
|
me, ni); |
244 |
|
|
airMopError(mop); return 1; |
245 |
|
|
} |
246 |
|
|
} |
247 |
|
|
if (save) { |
248 |
|
|
nrrdSave("regtmp-dwihist.nrrd", nhist, NULL); |
249 |
|
|
} |
250 |
|
|
/* |
251 |
|
|
if (_tenFindValley(DWthrP, nhist, 0.65, save)) { |
252 |
|
|
biffAddf(TEN, "%s: problem finding DWI histogram valley", me); |
253 |
|
|
airMopError(mop); return 1; |
254 |
|
|
} |
255 |
|
|
*/ |
256 |
|
|
if (nrrdHistoThresholdOtsu(DWthrP, nhist, expo)) { |
257 |
|
|
biffMovef(TEN, NRRD, "%s: problem finding DWI threshold", me); |
258 |
|
|
airMopError(mop); return 1; |
259 |
|
|
} |
260 |
|
|
|
261 |
|
|
airMopOkay(mop); |
262 |
|
|
return 0; |
263 |
|
|
} |
264 |
|
|
|
265 |
|
|
int |
266 |
|
|
_tenEpiRegThreshold(Nrrd **nthresh, Nrrd **nblur, unsigned int ninLen, |
267 |
|
|
double DWthr, int verb, int progress, double expo) { |
268 |
|
|
static const char me[]="_tenEpiRegThreshold"; |
269 |
|
|
airArray *mop; |
270 |
|
|
size_t I, sx, sy, sz, ni; |
271 |
|
|
float val; |
272 |
|
|
unsigned char *thr; |
273 |
|
|
|
274 |
|
|
if (!( AIR_EXISTS(DWthr) )) { |
275 |
|
|
if (_tenEpiRegThresholdFind(&DWthr, nblur, ninLen, progress, expo)) { |
276 |
|
|
biffAddf(TEN, "%s: trouble with automatic threshold determination", me); |
277 |
|
|
return 1; |
278 |
|
|
} |
279 |
|
|
fprintf(stderr, "%s: using %g for DWI threshold\n", me, DWthr); |
280 |
|
|
} |
281 |
|
|
|
282 |
|
|
mop = airMopNew(); |
283 |
|
|
if (verb) { |
284 |
|
|
fprintf(stderr, "%s:\n ", me); fflush(stderr); |
285 |
|
|
} |
286 |
|
|
sx = nblur[0]->axis[0].size; |
287 |
|
|
sy = nblur[0]->axis[1].size; |
288 |
|
|
sz = nblur[0]->axis[2].size; |
289 |
|
|
for (ni=0; ni<ninLen; ni++) { |
290 |
|
|
if (verb) { |
291 |
|
|
fprintf(stderr, "%2u ", (unsigned int)ni); fflush(stderr); |
292 |
|
|
} |
293 |
|
|
if (nrrdMaybeAlloc_va(nthresh[ni], nrrdTypeUChar, 3, |
294 |
|
|
sx, sy, sz)) { |
295 |
|
|
biffMovef(TEN, NRRD, "%s: trouble allocating threshold %u", |
296 |
|
|
me, (unsigned int)ni); |
297 |
|
|
airMopError(mop); return 1; |
298 |
|
|
} |
299 |
|
|
thr = (unsigned char *)(nthresh[ni]->data); |
300 |
|
|
for (I=0; I<sx*sy*sz; I++) { |
301 |
|
|
val = nrrdFLookup[nblur[ni]->type](nblur[ni]->data, I); |
302 |
|
|
val -= AIR_CAST(float, DWthr); |
303 |
|
|
thr[I] = (val >= 0 ? 1 : 0); |
304 |
|
|
} |
305 |
|
|
} |
306 |
|
|
if (verb) { |
307 |
|
|
fprintf(stderr, "done\n"); |
308 |
|
|
} |
309 |
|
|
|
310 |
|
|
airMopOkay(mop); |
311 |
|
|
return 0; |
312 |
|
|
} |
313 |
|
|
|
314 |
|
|
/* |
315 |
|
|
** _tenEpiRegBB: find the biggest bright CC |
316 |
|
|
*/ |
317 |
|
|
int |
318 |
|
|
_tenEpiRegBB(Nrrd *nval, Nrrd *nsize) { |
319 |
|
|
unsigned char *val; |
320 |
|
|
int *size, big; |
321 |
|
|
unsigned int ci; |
322 |
|
|
|
323 |
|
|
val = (unsigned char *)(nval->data); |
324 |
|
|
size = (int *)(nsize->data); |
325 |
|
|
big = 0; |
326 |
|
|
for (ci=0; ci<nsize->axis[0].size; ci++) { |
327 |
|
|
big = val[ci] ? AIR_MAX(big, size[ci]) : big; |
328 |
|
|
} |
329 |
|
|
return big; |
330 |
|
|
} |
331 |
|
|
|
332 |
|
|
int |
333 |
|
|
_tenEpiRegCC(Nrrd **nthr, int ninLen, int conny, int verb) { |
334 |
|
|
static const char me[]="_tenEpiRegCC"; |
335 |
|
|
Nrrd *nslc, *ncc, *nval, *nsize; |
336 |
|
|
airArray *mop; |
337 |
|
|
int ni, z, sz, big; |
338 |
|
|
|
339 |
|
|
mop = airMopNew(); |
340 |
|
|
airMopAdd(mop, nslc=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
341 |
|
|
airMopAdd(mop, nval=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
342 |
|
|
airMopAdd(mop, ncc=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
343 |
|
|
airMopAdd(mop, nsize=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
344 |
|
|
sz = nthr[0]->axis[2].size; |
345 |
|
|
if (verb) { |
346 |
|
|
fprintf(stderr, "%s:\n ", me); fflush(stderr); |
347 |
|
|
} |
348 |
|
|
for (ni=0; ni<ninLen; ni++) { |
349 |
|
|
if (verb) { |
350 |
|
|
fprintf(stderr, "%2d ", ni); fflush(stderr); |
351 |
|
|
} |
352 |
|
|
/* for each volume, we find the biggest bright 3-D CC, and merge |
353 |
|
|
down (to dark) all smaller bright pieces. Then, within each |
354 |
|
|
slice, we do 2-D CCs, find the biggest bright CC (size == big), |
355 |
|
|
and merge up (to bright) all small dark pieces, where |
356 |
|
|
(currently) small is big/2 */ |
357 |
|
|
big = -1; |
358 |
|
|
if (nrrdCCFind(ncc, &nval, nthr[ni], nrrdTypeDefault, conny) |
359 |
|
|
|| nrrdCCSize(nsize, ncc) |
360 |
|
|
|| !(big = _tenEpiRegBB(nval, nsize)) |
361 |
|
|
|| nrrdCCMerge(ncc, ncc, nval, -1, big-1, 0, conny) |
362 |
|
|
|| nrrdCCRevalue(nthr[ni], ncc, nval)) { |
363 |
|
|
if (big) { |
364 |
|
|
biffMovef(TEN, NRRD, |
365 |
|
|
"%s: trouble with 3-D processing nthr[%d]", me, ni); |
366 |
|
|
return 1; |
367 |
|
|
} else { |
368 |
|
|
biffAddf(TEN, "%s: got size 0 for biggest bright CC of nthr[%d]", |
369 |
|
|
me, ni); |
370 |
|
|
return 1; |
371 |
|
|
} |
372 |
|
|
} |
373 |
|
|
for (z=0; z<sz; z++) { |
374 |
|
|
big = -1; |
375 |
|
|
if ( nrrdSlice(nslc, nthr[ni], 2, z) |
376 |
|
|
|| nrrdCCFind(ncc, &nval, nslc, nrrdTypeDefault, conny) |
377 |
|
|
|| nrrdCCSize(nsize, ncc) |
378 |
|
|
|| !(big = _tenEpiRegBB(nval, nsize)) |
379 |
|
|
|| nrrdCCMerge(ncc, ncc, nval, 1, big/2, 0, conny) |
380 |
|
|
|| nrrdCCRevalue(nslc, ncc, nval) |
381 |
|
|
|| nrrdSplice(nthr[ni], nthr[ni], nslc, 2, z) ) { |
382 |
|
|
if (big) { |
383 |
|
|
biffMovef(TEN, NRRD, "%s: trouble processing slice %d of nthr[%d]", |
384 |
|
|
me, z, ni); |
385 |
|
|
return 1; |
386 |
|
|
} else { |
387 |
|
|
/* biggest bright CC on this slice had size 0 |
388 |
|
|
<==> there was no bright CC on this slice, move on */ |
389 |
|
|
} |
390 |
|
|
} |
391 |
|
|
} |
392 |
|
|
} |
393 |
|
|
if (verb) { |
394 |
|
|
fprintf(stderr, "done\n"); |
395 |
|
|
} |
396 |
|
|
|
397 |
|
|
airMopOkay(mop); |
398 |
|
|
return 0; |
399 |
|
|
} |
400 |
|
|
|
401 |
|
|
#define MEAN_X 0 |
402 |
|
|
#define MEAN_Y 1 |
403 |
|
|
#define M_02 2 |
404 |
|
|
#define M_11 3 |
405 |
|
|
#define M_20 4 |
406 |
|
|
|
407 |
|
|
/* |
408 |
|
|
** _tenEpiRegMoments() |
409 |
|
|
** |
410 |
|
|
** the moments are stored in (of course) a nrrd, one scanline per slice, |
411 |
|
|
** with each scanline containing: |
412 |
|
|
** |
413 |
|
|
** 0 1 2 3 4 |
414 |
|
|
** mean(x) mean(y) M_02 M_11 M_20 |
415 |
|
|
*/ |
416 |
|
|
int |
417 |
|
|
_tenEpiRegMoments(Nrrd **nmom, Nrrd **nthresh, unsigned int ninLen, |
418 |
|
|
int verb) { |
419 |
|
|
static const char me[]="_tenEpiRegMoments"; |
420 |
|
|
size_t sx, sy, sz, xi, yi, zi, ni; |
421 |
|
|
double N, mx, my, cx, cy, x, y, M02, M11, M20, *mom; |
422 |
|
|
float val; |
423 |
|
|
unsigned char *thr; |
424 |
|
|
|
425 |
|
|
sx = nthresh[0]->axis[0].size; |
426 |
|
|
sy = nthresh[0]->axis[1].size; |
427 |
|
|
sz = nthresh[0]->axis[2].size; |
428 |
|
|
if (verb) { |
429 |
|
|
fprintf(stderr, "%s:\n ", me); fflush(stderr); |
430 |
|
|
} |
431 |
|
|
for (ni=0; ni<ninLen; ni++) { |
432 |
|
|
if (verb) { |
433 |
|
|
fprintf(stderr, "%2u ", (unsigned int)ni); fflush(stderr); |
434 |
|
|
} |
435 |
|
|
if (nrrdMaybeAlloc_va(nmom[ni], nrrdTypeDouble, 2, |
436 |
|
|
AIR_CAST(size_t, 5), sz)) { |
437 |
|
|
biffMovef(TEN, NRRD, |
438 |
|
|
"%s: couldn't allocate nmom[%u]", me, (unsigned int)ni); |
439 |
|
|
return 1; |
440 |
|
|
} |
441 |
|
|
nrrdAxisInfoSet_va(nmom[ni], nrrdAxisInfoLabel, "mx,my,h,s,t", "z"); |
442 |
|
|
thr = (unsigned char *)(nthresh[ni]->data); |
443 |
|
|
mom = (double *)(nmom[ni]->data); |
444 |
|
|
for (zi=0; zi<sz; zi++) { |
445 |
|
|
/* ------ find mx, my */ |
446 |
|
|
N = 0; |
447 |
|
|
mx = my = 0.0; |
448 |
|
|
for (yi=0; yi<sy; yi++) { |
449 |
|
|
for (xi=0; xi<sx; xi++) { |
450 |
|
|
val = thr[xi + sx*yi]; |
451 |
|
|
N += val; |
452 |
|
|
mx += xi*val; |
453 |
|
|
my += yi*val; |
454 |
|
|
} |
455 |
|
|
} |
456 |
|
|
if (N == sx*sy) { |
457 |
|
|
biffAddf(TEN, "%s: saw only non-zero pixels in nthresh[%u]; " |
458 |
|
|
"DWI hreshold too low?", me, (unsigned int)ni); |
459 |
|
|
return 1; |
460 |
|
|
} |
461 |
|
|
if (N) { |
462 |
|
|
/* there were non-zero pixels */ |
463 |
|
|
mx /= N; |
464 |
|
|
my /= N; |
465 |
|
|
cx = sx/2.0; |
466 |
|
|
cy = sy/2.0; |
467 |
|
|
/* ------ find M02, M11, M20 */ |
468 |
|
|
M02 = M11 = M20 = 0.0; |
469 |
|
|
for (yi=0; yi<sy; yi++) { |
470 |
|
|
for (xi=0; xi<sx; xi++) { |
471 |
|
|
val = thr[xi + sx*yi]; |
472 |
|
|
x = xi - cx; |
473 |
|
|
y = yi - cy; |
474 |
|
|
M02 += y*y*val; |
475 |
|
|
M11 += x*y*val; |
476 |
|
|
M20 += x*x*val; |
477 |
|
|
} |
478 |
|
|
} |
479 |
|
|
M02 /= N; |
480 |
|
|
M11 /= N; |
481 |
|
|
M20 /= N; |
482 |
|
|
/* ------ set output */ |
483 |
|
|
mom[MEAN_X] = mx; |
484 |
|
|
mom[MEAN_Y] = my; |
485 |
|
|
mom[M_02] = M02; |
486 |
|
|
mom[M_11] = M11; |
487 |
|
|
mom[M_20] = M20; |
488 |
|
|
} else { |
489 |
|
|
/* there were no non-zero pixels */ |
490 |
|
|
mom[MEAN_X] = 0; |
491 |
|
|
mom[MEAN_Y] = 0; |
492 |
|
|
mom[M_02] = 0; |
493 |
|
|
mom[M_11] = 0; |
494 |
|
|
mom[M_20] = 0; |
495 |
|
|
} |
496 |
|
|
thr += sx*sy; |
497 |
|
|
mom += 5; |
498 |
|
|
} |
499 |
|
|
} |
500 |
|
|
if (verb) { |
501 |
|
|
fprintf(stderr, "done\n"); |
502 |
|
|
} |
503 |
|
|
|
504 |
|
|
return 0; |
505 |
|
|
} |
506 |
|
|
|
507 |
|
|
/* |
508 |
|
|
** _tenEpiRegPairXforms |
509 |
|
|
** |
510 |
|
|
** uses moment information to compute all pair-wise transforms, which are |
511 |
|
|
** stored in the 3 x ninLen x ninLen x sizeZ output. If xfr = npxfr->data, |
512 |
|
|
** xfr[0 + 3*(zi + sz*(A + ninLen*B))] is shear, |
513 |
|
|
** xfr[1 + " ] is scale, and |
514 |
|
|
** xfr[2 + " ] is translate in the transform |
515 |
|
|
** that maps slice zi from volume A to volume B. |
516 |
|
|
*/ |
517 |
|
|
int |
518 |
|
|
_tenEpiRegPairXforms(Nrrd *npxfr, Nrrd **nmom, int ninLen) { |
519 |
|
|
static const char me[]="_tenEpiRegPairXforms"; |
520 |
|
|
double *xfr, *A, *B, hh, ss, tt; |
521 |
|
|
int ai, bi, zi, sz; |
522 |
|
|
|
523 |
|
|
sz = nmom[0]->axis[1].size; |
524 |
|
|
if (nrrdMaybeAlloc_va(npxfr, nrrdTypeDouble, 4, |
525 |
|
|
AIR_CAST(size_t, 5), |
526 |
|
|
AIR_CAST(size_t, sz), |
527 |
|
|
AIR_CAST(size_t, ninLen), |
528 |
|
|
AIR_CAST(size_t, ninLen))) { |
529 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate transform nrrd", me); |
530 |
|
|
return 1; |
531 |
|
|
} |
532 |
|
|
nrrdAxisInfoSet_va(npxfr, nrrdAxisInfoLabel, |
533 |
|
|
"mx,my,h,s,t", "zi", "orig", "target"); |
534 |
|
|
xfr = (double *)(npxfr->data); |
535 |
|
|
for (bi=0; bi<ninLen; bi++) { |
536 |
|
|
for (ai=0; ai<ninLen; ai++) { |
537 |
|
|
for (zi=0; zi<sz; zi++) { |
538 |
|
|
A = (double*)(nmom[ai]->data) + 5*zi; |
539 |
|
|
B = (double*)(nmom[bi]->data) + 5*zi; |
540 |
|
|
ss = sqrt((A[M_20]*B[M_02] - B[M_11]*B[M_11]) / |
541 |
|
|
(A[M_20]*A[M_02] - A[M_11]*A[M_11])); |
542 |
|
|
hh = (B[M_11] - ss*A[M_11])/A[M_20]; |
543 |
|
|
tt = B[MEAN_Y] - A[MEAN_Y]; |
544 |
|
|
ELL_5V_SET(xfr + 5*(zi + sz*(ai + ninLen*bi)), |
545 |
|
|
A[MEAN_X], A[MEAN_Y], hh, ss, tt); |
546 |
|
|
} |
547 |
|
|
} |
548 |
|
|
} |
549 |
|
|
return 0; |
550 |
|
|
} |
551 |
|
|
|
552 |
|
|
#define SHEAR 2 |
553 |
|
|
#define SCALE 3 |
554 |
|
|
#define TRAN 4 |
555 |
|
|
|
556 |
|
|
int |
557 |
|
|
_tenEpiRegEstimHST(Nrrd *nhst, Nrrd *npxfr, int ninLen, Nrrd *ngrad) { |
558 |
|
|
static const char me[]="_tenEpiRegEstimHST"; |
559 |
|
|
double *hst, *grad, *mat, *vec, *ans, *pxfr, *gA, *gB; |
560 |
|
|
int z, sz, A, B, npairs, ri; |
561 |
|
|
Nrrd **nmat, *nvec, **ninv, *nans; |
562 |
|
|
airArray *mop; |
563 |
|
|
int order; |
564 |
|
|
|
565 |
|
|
order = 1; |
566 |
|
|
|
567 |
|
|
sz = npxfr->axis[1].size; |
568 |
|
|
npairs = ninLen*(ninLen-1); |
569 |
|
|
|
570 |
|
|
mop = airMopNew(); |
571 |
|
|
nmat = (Nrrd**)calloc(sz, sizeof(Nrrd*)); |
572 |
|
|
ninv = (Nrrd**)calloc(sz, sizeof(Nrrd*)); |
573 |
|
|
airMopAdd(mop, nmat, airFree, airMopAlways); |
574 |
|
|
airMopAdd(mop, ninv, airFree, airMopAlways); |
575 |
|
|
for (z=0; z<sz; z++) { |
576 |
|
|
airMopAdd(mop, nmat[z]=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
577 |
|
|
if (nrrdMaybeAlloc_va(nmat[z], nrrdTypeDouble, 2, |
578 |
|
|
AIR_CAST(size_t, (1 == order ? 3 : 9)), |
579 |
|
|
AIR_CAST(size_t, npairs))) { |
580 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate fitting matrices", me); |
581 |
|
|
airMopError(mop); return 1; |
582 |
|
|
} |
583 |
|
|
airMopAdd(mop, ninv[z]=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
584 |
|
|
} |
585 |
|
|
airMopAdd(mop, nvec=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
586 |
|
|
airMopAdd(mop, nans=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
587 |
|
|
if (nrrdMaybeAlloc_va(nhst, nrrdTypeDouble, 2, |
588 |
|
|
AIR_CAST(size_t, (1 == order ? 9 : 27)), |
589 |
|
|
AIR_CAST(size_t, sz)) |
590 |
|
|
|| nrrdMaybeAlloc_va(nvec, nrrdTypeDouble, 2, |
591 |
|
|
AIR_CAST(size_t, 1), |
592 |
|
|
AIR_CAST(size_t, npairs))) { |
593 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate HST nrrd", me); |
594 |
|
|
airMopError(mop); return 1; |
595 |
|
|
} |
596 |
|
|
nrrdAxisInfoSet_va(nhst, nrrdAxisInfoLabel, |
597 |
|
|
(1 == order |
598 |
|
|
? "Hx,Hy,Hz,Sx,Sy,Sz,Tx,Ty,Tz" |
599 |
|
|
: "HST parms"), "z"); |
600 |
|
|
|
601 |
|
|
/* ------ per-slice: compute model fitting matrix and its pseudo-inverse */ |
602 |
|
|
grad = (double *)(ngrad->data); |
603 |
|
|
for (z=0; z<sz; z++) { |
604 |
|
|
hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z; |
605 |
|
|
mat = (double *)(nmat[z]->data); |
606 |
|
|
ri = 0; |
607 |
|
|
for (A=0; A<ninLen; A++) { |
608 |
|
|
for (B=0; B<ninLen; B++) { |
609 |
|
|
if (A == B) continue; |
610 |
|
|
pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B)); |
611 |
|
|
gA = grad + 0 + 3*A; |
612 |
|
|
gB = grad + 0 + 3*B; |
613 |
|
|
if (1 == order) { |
614 |
|
|
ELL_3V_SET(mat + 3*ri, |
615 |
|
|
gB[0] - pxfr[SCALE]*gA[0], |
616 |
|
|
gB[1] - pxfr[SCALE]*gA[1], |
617 |
|
|
gB[2] - pxfr[SCALE]*gA[2]); |
618 |
|
|
} else { |
619 |
|
|
ELL_9V_SET(mat + 9*ri, |
620 |
|
|
gB[0] - pxfr[SCALE]*gA[0], |
621 |
|
|
gB[1] - pxfr[SCALE]*gA[1], |
622 |
|
|
gB[2] - pxfr[SCALE]*gA[2], |
623 |
|
|
gB[0]*gB[0]*gB[0] - pxfr[SCALE]*gA[0]*gA[0]*gA[0], |
624 |
|
|
gB[1]*gB[1]*gB[1] - pxfr[SCALE]*gA[1]*gA[1]*gA[1], |
625 |
|
|
gB[2]*gB[2]*gB[2] - pxfr[SCALE]*gA[2]*gA[2]*gA[2], |
626 |
|
|
gB[0]*gB[1]*gB[2] - pxfr[SCALE]*gA[0]*gA[1]*gA[2], |
627 |
|
|
1, 1); |
628 |
|
|
/* |
629 |
|
|
gB[0]*gB[1] - pxfr[SCALE]*gA[0]*gA[1], |
630 |
|
|
gB[0]*gB[2] - pxfr[SCALE]*gA[0]*gA[2], |
631 |
|
|
gB[1]*gB[2] - pxfr[SCALE]*gA[1]*gA[2]); |
632 |
|
|
*/ |
633 |
|
|
} |
634 |
|
|
ri += 1; |
635 |
|
|
} |
636 |
|
|
} |
637 |
|
|
if (nrrdHasNonExist(nmat[z])) { |
638 |
|
|
/* as happens if there were zero slices in the segmentation output */ |
639 |
|
|
if (1 == order) { |
640 |
|
|
ELL_3V_SET(hst + 0*3, 0, 0, 0); |
641 |
|
|
ELL_3V_SET(hst + 1*3, 0, 0, 0); |
642 |
|
|
ELL_3V_SET(hst + 2*3, 0, 0, 0); |
643 |
|
|
} else { |
644 |
|
|
ELL_9V_SET(hst + 0*9, 0, 0, 0, 0, 0, 0, 0, 0, 0); |
645 |
|
|
ELL_9V_SET(hst + 1*9, 0, 0, 0, 0, 0, 0, 0, 0, 0); |
646 |
|
|
ELL_9V_SET(hst + 2*9, 0, 0, 0, 0, 0, 0, 0, 0, 0); |
647 |
|
|
} |
648 |
|
|
} else { |
649 |
|
|
if (ell_Nm_pseudo_inv(ninv[z], nmat[z])) { |
650 |
|
|
biffMovef(TEN, ELL, "%s: trouble estimating model (slice %d)", me, z); |
651 |
|
|
airMopError(mop); return 1; |
652 |
|
|
} |
653 |
|
|
} |
654 |
|
|
} |
655 |
|
|
vec = (double *)(nvec->data); |
656 |
|
|
|
657 |
|
|
/* ------ find Sx, Sy, Sz per slice */ |
658 |
|
|
for (z=0; z<sz; z++) { |
659 |
|
|
if (nrrdHasNonExist(nmat[z])) { |
660 |
|
|
/* we've already zero-ed out this row in the HST nrrd */ |
661 |
|
|
continue; |
662 |
|
|
} |
663 |
|
|
hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z; |
664 |
|
|
ri = 0; |
665 |
|
|
for (A=0; A<ninLen; A++) { |
666 |
|
|
for (B=0; B<ninLen; B++) { |
667 |
|
|
if (A == B) continue; |
668 |
|
|
pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B)); |
669 |
|
|
vec[ri] = pxfr[SCALE] - 1; |
670 |
|
|
ri += 1; |
671 |
|
|
} |
672 |
|
|
} |
673 |
|
|
if (ell_Nm_mul(nans, ninv[z], nvec)) { |
674 |
|
|
biffMovef(TEN, ELL, |
675 |
|
|
"%s: trouble estimating model (slice %d): Sx, Sy, Sz", |
676 |
|
|
me, z); |
677 |
|
|
airMopError(mop); return 1; |
678 |
|
|
} |
679 |
|
|
ans = (double *)(nans->data); |
680 |
|
|
if (1 == order) { |
681 |
|
|
ELL_3V_COPY(hst + 1*3, ans); |
682 |
|
|
} else { |
683 |
|
|
ELL_9V_COPY(hst + 1*9, ans); |
684 |
|
|
} |
685 |
|
|
} |
686 |
|
|
|
687 |
|
|
/* ------ find Hx, Hy, Hz per slice */ |
688 |
|
|
for (z=0; z<sz; z++) { |
689 |
|
|
if (nrrdHasNonExist(nmat[z])) { |
690 |
|
|
/* we've already zero-ed out this row in the HST nrrd */ |
691 |
|
|
continue; |
692 |
|
|
} |
693 |
|
|
hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z; |
694 |
|
|
ri = 0; |
695 |
|
|
for (A=0; A<ninLen; A++) { |
696 |
|
|
for (B=0; B<ninLen; B++) { |
697 |
|
|
if (A == B) continue; |
698 |
|
|
pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B)); |
699 |
|
|
vec[ri] = pxfr[SHEAR]; |
700 |
|
|
ri += 1; |
701 |
|
|
} |
702 |
|
|
} |
703 |
|
|
if (ell_Nm_mul(nans, ninv[z], nvec)) { |
704 |
|
|
biffAddf(TEN, ELL, "%s: trouble estimating model (slice %d): Hx, Hy, Hz", |
705 |
|
|
me, z); |
706 |
|
|
airMopError(mop); return 1; |
707 |
|
|
} |
708 |
|
|
ans = (double *)(nans->data); |
709 |
|
|
if (1 == order) { |
710 |
|
|
ELL_3V_COPY(hst + 0*3, ans); |
711 |
|
|
} else { |
712 |
|
|
ELL_9V_COPY(hst + 0*9, ans); |
713 |
|
|
} |
714 |
|
|
} |
715 |
|
|
|
716 |
|
|
/* ------ find Tx, Ty, Tz per slice */ |
717 |
|
|
for (z=0; z<sz; z++) { |
718 |
|
|
if (nrrdHasNonExist(nmat[z])) { |
719 |
|
|
/* we've already zero-ed out this row in the HST nrrd */ |
720 |
|
|
continue; |
721 |
|
|
} |
722 |
|
|
hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z; |
723 |
|
|
ri = 0; |
724 |
|
|
for (A=0; A<ninLen; A++) { |
725 |
|
|
for (B=0; B<ninLen; B++) { |
726 |
|
|
if (A == B) continue; |
727 |
|
|
pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B)); |
728 |
|
|
vec[ri] = pxfr[TRAN]; |
729 |
|
|
ri += 1; |
730 |
|
|
} |
731 |
|
|
} |
732 |
|
|
if (ell_Nm_mul(nans, ninv[z], nvec)) { |
733 |
|
|
biffMovef(TEN, ELL, |
734 |
|
|
"%s: trouble estimating model (slice %d): Tx, Ty, Tz", |
735 |
|
|
me, z); |
736 |
|
|
airMopError(mop); return 1; |
737 |
|
|
} |
738 |
|
|
ans = (double *)(nans->data); |
739 |
|
|
if (1 == order) { |
740 |
|
|
ELL_3V_COPY(hst + 2*3, ans); |
741 |
|
|
} else { |
742 |
|
|
ELL_9V_COPY(hst + 2*9, ans); |
743 |
|
|
} |
744 |
|
|
} |
745 |
|
|
|
746 |
|
|
airMopOkay(mop); |
747 |
|
|
return 0; |
748 |
|
|
} |
749 |
|
|
|
750 |
|
|
int |
751 |
|
|
_tenEpiRegFitHST(Nrrd *nhst, Nrrd **_ncc, int ninLen, |
752 |
|
|
double goodFrac, int prog, int verb) { |
753 |
|
|
static const char me[]="_tenEpiRegFitHST"; |
754 |
|
|
airArray *mop; |
755 |
|
|
Nrrd *ncc, *ntA, *ntB, *nsd, *nl2; |
756 |
|
|
unsigned int cc, sz, zi, sh, hi; |
757 |
|
|
float *mess, *two, tmp; |
758 |
|
|
double *hst, x, y, xx, xy, mm, bb; |
759 |
|
|
|
760 |
|
|
mop = airMopNew(); |
761 |
|
|
airMopAdd(mop, ncc=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
762 |
|
|
airMopAdd(mop, ntA=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
763 |
|
|
airMopAdd(mop, ntB=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
764 |
|
|
airMopAdd(mop, nsd=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
765 |
|
|
airMopAdd(mop, nl2=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
766 |
|
|
/* do SD and L2 projections of the CCs along the DWI axis, |
767 |
|
|
integrate these over the X and Y axes of the slices, |
768 |
|
|
and define per-slice "messiness" as the quotient of the |
769 |
|
|
SD integral with the L2 integral */ |
770 |
|
|
if (verb) { |
771 |
|
|
fprintf(stderr, "%s: measuring segmentation uncertainty ... ", me); |
772 |
|
|
fflush(stderr); |
773 |
|
|
} |
774 |
|
|
if (nrrdJoin(ncc, (const Nrrd*const*)_ncc, ninLen, 0, AIR_TRUE) |
775 |
|
|
|| nrrdProject(ntA, ncc, 0, nrrdMeasureSD, nrrdTypeFloat) |
776 |
|
|
|| nrrdProject(ntB, ntA, 0, nrrdMeasureSum, nrrdTypeFloat) |
777 |
|
|
|| nrrdProject(nsd, ntB, 0, nrrdMeasureSum, nrrdTypeFloat) |
778 |
|
|
|| nrrdProject(ntA, ncc, 0, nrrdMeasureL2, nrrdTypeFloat) |
779 |
|
|
|| nrrdProject(ntB, ntA, 0, nrrdMeasureSum, nrrdTypeFloat) |
780 |
|
|
|| nrrdProject(nl2, ntB, 0, nrrdMeasureSum, nrrdTypeFloat) |
781 |
|
|
|| nrrdArithBinaryOp(ntA, nrrdBinaryOpDivide, nsd, nl2)) { |
782 |
|
|
biffMovef(TEN, NRRD, "%s: trouble doing CC projections", me); |
783 |
|
|
airMopError(mop); return 1; |
784 |
|
|
} |
785 |
|
|
if (verb) { |
786 |
|
|
fprintf(stderr, "done\n"); |
787 |
|
|
} |
788 |
|
|
if (prog && _tenEpiRegSave("regtmp-messy.txt", ntA, |
789 |
|
|
NULL, 0, "segmentation uncertainty")) { |
790 |
|
|
biffMovef(TEN, NRRD, "%s: EpiRegSave failed", me); |
791 |
|
|
airMopError(mop); return 1; |
792 |
|
|
} |
793 |
|
|
/* now ntA stores the per-slice messiness */ |
794 |
|
|
mess = AIR_CAST(float*, ntA->data); |
795 |
|
|
|
796 |
|
|
/* allocate an array of 2 floats per slice */ |
797 |
|
|
sz = ntA->axis[0].size; |
798 |
|
|
two = AIR_CAST(float*, calloc(2*sz, sizeof(float))); |
799 |
|
|
if (!two) { |
800 |
|
|
biffAddf(TEN, "%s: couldn't allocate tmp buffer", me); |
801 |
|
|
airMopError(mop); return 1; |
802 |
|
|
} |
803 |
|
|
airMopAdd(mop, two, airFree, airMopAlways); |
804 |
|
|
/* initial ordering: messiness, slice index */ |
805 |
|
|
for (zi=0; zi<sz; zi++) { |
806 |
|
|
two[0 + 2*zi] = (AIR_EXISTS(mess[zi]) |
807 |
|
|
? mess[zi] |
808 |
|
|
: 666); /* don't use empty slices */ |
809 |
|
|
two[1 + 2*zi] = AIR_CAST(float, zi); |
810 |
|
|
} |
811 |
|
|
/* sort into ascending messiness */ |
812 |
|
|
qsort(two, zi, 2*sizeof(float), nrrdValCompare[nrrdTypeFloat]); |
813 |
|
|
/* flip ordering while thresholding messiness into usability */ |
814 |
|
|
for (zi=0; zi<sz; zi++) { |
815 |
|
|
tmp = two[1 + 2*zi]; |
816 |
|
|
two[1 + 2*zi] = AIR_AFFINE(0, zi, sz-1, 0, 1) <= goodFrac ? 1.0f : 0.0f; |
817 |
|
|
two[0 + 2*zi] = tmp; |
818 |
|
|
} |
819 |
|
|
/* sort again, now into ascending slice order */ |
820 |
|
|
qsort(two, zi, 2*sizeof(float), nrrdValCompare[nrrdTypeFloat]); |
821 |
|
|
if (verb) { |
822 |
|
|
fprintf(stderr, "%s: using slices", me); |
823 |
|
|
for (zi=0; zi<sz; zi++) { |
824 |
|
|
if (two[1 + 2*zi]) { |
825 |
|
|
fprintf(stderr, " %u", zi); |
826 |
|
|
} |
827 |
|
|
} |
828 |
|
|
fprintf(stderr, " for fitting\n"); |
829 |
|
|
} |
830 |
|
|
|
831 |
|
|
/* perform fitting for each column in hst (regardless of |
832 |
|
|
whether we're using a 1st or 2nd order model */ |
833 |
|
|
hst = (double*)(nhst->data); |
834 |
|
|
sh = nhst->axis[0].size; |
835 |
|
|
for (hi=0; hi<sh; hi++) { |
836 |
|
|
x = y = xy = xx = 0; |
837 |
|
|
cc = 0; |
838 |
|
|
for (zi=0; zi<sz; zi++) { |
839 |
|
|
if (!two[1 + 2*zi]) |
840 |
|
|
continue; |
841 |
|
|
cc += 1; |
842 |
|
|
x += zi; |
843 |
|
|
xx += zi*zi; |
844 |
|
|
y += hst[hi + sh*zi]; |
845 |
|
|
xy += zi*hst[hi + sh*zi]; |
846 |
|
|
} |
847 |
|
|
x /= cc; xx /= cc; y /= cc; xy /= cc; |
848 |
|
|
mm = (xy - x*y)/(xx - x*x); |
849 |
|
|
bb = y - mm*x; |
850 |
|
|
for (zi=0; zi<sz; zi++) { |
851 |
|
|
hst[hi + sh*zi] = mm*zi + bb; |
852 |
|
|
} |
853 |
|
|
} |
854 |
|
|
|
855 |
|
|
airMopOkay(mop); |
856 |
|
|
return 0; |
857 |
|
|
} |
858 |
|
|
|
859 |
|
|
int |
860 |
|
|
_tenEpiRegGetHST(double *hhP, double *ssP, double *ttP, |
861 |
|
|
int reference, int ni, int zi, |
862 |
|
|
Nrrd *npxfr, Nrrd *nhst, Nrrd *ngrad) { |
863 |
|
|
double *xfr, *hst, *_g, grad[9]; /* big enough for 2nd order */ |
864 |
|
|
int sz, ninLen; |
865 |
|
|
|
866 |
|
|
int order; |
867 |
|
|
|
868 |
|
|
order = 1; |
869 |
|
|
|
870 |
|
|
/* these could also have been passed to us, but we can also discover them */ |
871 |
|
|
sz = npxfr->axis[1].size; |
872 |
|
|
ninLen = npxfr->axis[2].size; |
873 |
|
|
|
874 |
|
|
if (-1 == reference) { |
875 |
|
|
/* we use the estimated H,S,T vectors to determine distortion |
876 |
|
|
as a function of gradient direction, and then invert this */ |
877 |
|
|
_g = (double*)(ngrad->data) + 0 + 3*ni; |
878 |
|
|
if (1 == order) { |
879 |
|
|
hst = (double*)(nhst->data) + 0 + 9*zi; |
880 |
|
|
*hhP = ELL_3V_DOT(_g, hst + 0*3); |
881 |
|
|
*ssP = 1 + ELL_3V_DOT(_g, hst + 1*3); |
882 |
|
|
*ttP = ELL_3V_DOT(_g, hst + 2*3); |
883 |
|
|
} else { |
884 |
|
|
hst = (double*)(nhst->data) + 0 + 27*zi; |
885 |
|
|
ELL_9V_SET(grad, _g[0], _g[1], _g[2], |
886 |
|
|
_g[0]*_g[0]*_g[0], _g[1]*_g[1]*_g[1], _g[2]*_g[2]*_g[2], |
887 |
|
|
_g[0]*_g[1]*_g[2], |
888 |
|
|
0, 0); |
889 |
|
|
/* |
890 |
|
|
_g[0]*_g[1], _g[0]*_g[2], _g[1]*_g[2]); |
891 |
|
|
*/ |
892 |
|
|
*hhP = ELL_9V_DOT(grad, hst + 0*9); |
893 |
|
|
*ssP = 1 + ELL_9V_DOT(grad, hst + 1*9); |
894 |
|
|
*ttP = ELL_9V_DOT(grad, hst + 2*9); |
895 |
|
|
} |
896 |
|
|
} else { |
897 |
|
|
/* we register against a specific DWI */ |
898 |
|
|
xfr = (double*)(npxfr->data) + 0 + 5*(zi + sz*(reference + ninLen*ni)); |
899 |
|
|
*hhP = xfr[2]; |
900 |
|
|
*ssP = xfr[3]; |
901 |
|
|
*ttP = xfr[4]; |
902 |
|
|
} |
903 |
|
|
return 0; |
904 |
|
|
} |
905 |
|
|
|
906 |
|
|
/* |
907 |
|
|
** _tenEpiRegSliceWarp |
908 |
|
|
** |
909 |
|
|
** Apply [hh,ss,tt] transform to nin, putting results in nout, but with |
910 |
|
|
** some trickiness: |
911 |
|
|
** - nwght and nidx are already allocated to the the weights (type float) |
912 |
|
|
** and indices for resampling nin with "kern" and "kparm" |
913 |
|
|
** - nout is already allocated to the correct size and type |
914 |
|
|
** - nin is type float, but output must be type nout->type |
915 |
|
|
** - nin is been transposed to have the resampled axis fastest in memory, |
916 |
|
|
** but nout output will not be transposed |
917 |
|
|
*/ |
918 |
|
|
int |
919 |
|
|
_tenEpiRegSliceWarp(Nrrd *nout, Nrrd *nin, Nrrd *nwght, Nrrd *nidx, |
920 |
|
|
const NrrdKernel *kern, double *kparm, |
921 |
|
|
double hh, double ss, double tt, double cx, double cy) { |
922 |
|
|
float *wght, *in, pp, pf, tmp; |
923 |
|
|
int *idx; |
924 |
|
|
unsigned int supp; |
925 |
|
|
size_t sx, sy, xi, yi, pb, pi; |
926 |
|
|
double (*ins)(void *, size_t, double), (*clamp)(double); |
927 |
|
|
|
928 |
|
|
sy = nin->axis[0].size; |
929 |
|
|
sx = nin->axis[1].size; |
930 |
|
|
supp = AIR_CAST(unsigned int, kern->support(kparm)); |
931 |
|
|
ins = nrrdDInsert[nout->type]; |
932 |
|
|
clamp = nrrdDClamp[nout->type]; |
933 |
|
|
|
934 |
|
|
in = AIR_CAST(float*, nin->data); |
935 |
|
|
for (xi=0; xi<sx; xi++) { |
936 |
|
|
idx = AIR_CAST(int*, nidx->data); |
937 |
|
|
wght = AIR_CAST(float*, nwght->data); |
938 |
|
|
for (yi=0; yi<sy; yi++) { |
939 |
|
|
pp = AIR_CAST(float, hh*(xi - cx) + ss*(yi - cy) + tt + cy); |
940 |
|
|
pb = AIR_CAST(size_t, floor(pp)); |
941 |
|
|
pf = pp - pb; |
942 |
|
|
for (pi=0; pi<2*supp; pi++) { |
943 |
|
|
idx[pi] = AIR_MIN(pb + pi - (supp-1), sy-1); |
944 |
|
|
wght[pi] = pi - (supp-1) - pf; |
945 |
|
|
} |
946 |
|
|
idx += 2*supp; |
947 |
|
|
wght += 2*supp; |
948 |
|
|
} |
949 |
|
|
idx = AIR_CAST(int*, nidx->data); |
950 |
|
|
wght = AIR_CAST(float*, nwght->data); |
951 |
|
|
kern->evalN_f(wght, wght, 2*supp*sy, kparm); |
952 |
|
|
for (yi=0; yi<sy; yi++) { |
953 |
|
|
tmp = 0; |
954 |
|
|
for (pi=0; pi<2*supp; pi++) { |
955 |
|
|
tmp += in[idx[pi]]*wght[pi]; |
956 |
|
|
} |
957 |
|
|
ins(nout->data, xi + sx*yi, clamp(ss*tmp)); |
958 |
|
|
idx += 2*supp; |
959 |
|
|
wght += 2*supp; |
960 |
|
|
} |
961 |
|
|
in += sy; |
962 |
|
|
} |
963 |
|
|
|
964 |
|
|
return 0; |
965 |
|
|
} |
966 |
|
|
|
967 |
|
|
/* |
968 |
|
|
** _tenEpiRegWarp() |
969 |
|
|
** |
970 |
|
|
*/ |
971 |
|
|
int |
972 |
|
|
_tenEpiRegWarp(Nrrd **ndone, Nrrd *npxfr, Nrrd *nhst, Nrrd *ngrad, |
973 |
|
|
Nrrd **nin, int ninLen, |
974 |
|
|
int reference, const NrrdKernel *kern, double *kparm, |
975 |
|
|
int verb) { |
976 |
|
|
static const char me[]="_tenEpiRegWarp"; |
977 |
|
|
Nrrd *ntmp, *nfin, *nslcA, *nslcB, *nwght, *nidx; |
978 |
|
|
airArray *mop; |
979 |
|
|
int sx, sy, sz, ni, zi, supp; |
980 |
|
|
double hh, ss, tt, cx, cy; |
981 |
|
|
|
982 |
|
|
mop = airMopNew(); |
983 |
|
|
airMopAdd(mop, ntmp=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
984 |
|
|
airMopAdd(mop, nfin=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
985 |
|
|
airMopAdd(mop, nslcA=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
986 |
|
|
airMopAdd(mop, nslcB=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
987 |
|
|
airMopAdd(mop, nwght=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
988 |
|
|
airMopAdd(mop, nidx=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
989 |
|
|
|
990 |
|
|
if (verb) { |
991 |
|
|
fprintf(stderr, "%s:\n ", me); fflush(stderr); |
992 |
|
|
} |
993 |
|
|
sx = nin[0]->axis[0].size; |
994 |
|
|
sy = nin[0]->axis[1].size; |
995 |
|
|
sz = nin[0]->axis[2].size; |
996 |
|
|
cx = sx/2.0; |
997 |
|
|
cy = sy/2.0; |
998 |
|
|
supp = (int)kern->support(kparm); |
999 |
|
|
if (nrrdMaybeAlloc_va(nwght, nrrdTypeFloat, 2, |
1000 |
|
|
AIR_CAST(size_t, 2*supp), |
1001 |
|
|
AIR_CAST(size_t, sy)) |
1002 |
|
|
|| nrrdMaybeAlloc_va(nidx, nrrdTypeInt, 2, |
1003 |
|
|
AIR_CAST(size_t, 2*supp), |
1004 |
|
|
AIR_CAST(size_t, sy))) { |
1005 |
|
|
biffMovef(TEN, NRRD, "%s: trouble allocating buffers", me); |
1006 |
|
|
airMopError(mop); return 1; |
1007 |
|
|
} |
1008 |
|
|
for (ni=0; ni<ninLen; ni++) { |
1009 |
|
|
if (verb) { |
1010 |
|
|
fprintf(stderr, "%2d ", ni); fflush(stderr); |
1011 |
|
|
} |
1012 |
|
|
if (nrrdCopy(ndone[ni], nin[ni]) |
1013 |
|
|
|| ((!ni) && nrrdSlice(nslcB, ndone[ni], 2, 0)) /* only when 0==ni */ |
1014 |
|
|
|| nrrdAxesSwap(ntmp, nin[ni], 0, 1) |
1015 |
|
|
|| nrrdConvert(nfin, ntmp, nrrdTypeFloat)) { |
1016 |
|
|
biffMovef(TEN, NRRD, "%s: trouble prepping at ni=%d", me, ni); |
1017 |
|
|
airMopError(mop); return 1; |
1018 |
|
|
} |
1019 |
|
|
for (zi=0; zi<sz; zi++) { |
1020 |
|
|
if (_tenEpiRegGetHST(&hh, &ss, &tt, reference, |
1021 |
|
|
ni, zi, npxfr, nhst, ngrad) |
1022 |
|
|
|| nrrdSlice(nslcA, nfin, 2, zi) |
1023 |
|
|
|| _tenEpiRegSliceWarp(nslcB, nslcA, nwght, nidx, kern, kparm, |
1024 |
|
|
hh, ss, tt, cx, cy) |
1025 |
|
|
|| nrrdSplice(ndone[ni], ndone[ni], nslcB, 2, zi)) { |
1026 |
|
|
biffMovef(TEN, NRRD, "%s: trouble on slice %d if ni=%d", me, zi, ni); |
1027 |
|
|
/* because the _tenEpiReg calls above don't use biff */ |
1028 |
|
|
airMopError(mop); return 1; |
1029 |
|
|
} |
1030 |
|
|
} |
1031 |
|
|
} |
1032 |
|
|
if (verb) { |
1033 |
|
|
fprintf(stderr, "done\n"); |
1034 |
|
|
} |
1035 |
|
|
|
1036 |
|
|
airMopOkay(mop); |
1037 |
|
|
return 0; |
1038 |
|
|
} |
1039 |
|
|
|
1040 |
|
|
int |
1041 |
|
|
tenEpiRegister3D(Nrrd **nout, Nrrd **nin, unsigned int ninLen, Nrrd *_ngrad, |
1042 |
|
|
int reference, |
1043 |
|
|
double bwX, double bwY, double fitFrac, |
1044 |
|
|
double DWthr, int doCC, |
1045 |
|
|
const NrrdKernel *kern, double *kparm, |
1046 |
|
|
int progress, int verbose) { |
1047 |
|
|
static const char me[]="tenEpiRegister3D"; |
1048 |
|
|
airArray *mop; |
1049 |
|
|
Nrrd **nbuffA, **nbuffB, *npxfr, *nhst, *ngrad; |
1050 |
|
|
int hack1, hack2; |
1051 |
|
|
unsigned int i; |
1052 |
|
|
|
1053 |
|
|
hack1 = nrrdStateAlwaysSetContent; |
1054 |
|
|
hack2 = nrrdStateDisableContent; |
1055 |
|
|
nrrdStateAlwaysSetContent = AIR_FALSE; |
1056 |
|
|
nrrdStateDisableContent = AIR_TRUE; |
1057 |
|
|
|
1058 |
|
|
mop = airMopNew(); |
1059 |
|
|
if (_tenEpiRegCheck(nout, nin, ninLen, _ngrad, reference, |
1060 |
|
|
bwX, bwY, DWthr, |
1061 |
|
|
kern, kparm)) { |
1062 |
|
|
biffAddf(TEN, "%s: trouble with input", me); |
1063 |
|
|
airMopError(mop); return 1; |
1064 |
|
|
} |
1065 |
|
|
|
1066 |
|
|
nbuffA = (Nrrd **)calloc(ninLen, sizeof(Nrrd*)); |
1067 |
|
|
nbuffB = (Nrrd **)calloc(ninLen, sizeof(Nrrd*)); |
1068 |
|
|
if (!( nbuffA && nbuffB )) { |
1069 |
|
|
biffAddf(TEN, "%s: couldn't allocate tmp nrrd pointer arrays", me); |
1070 |
|
|
airMopError(mop); return 1; |
1071 |
|
|
} |
1072 |
|
|
airMopAdd(mop, nbuffA, airFree, airMopAlways); |
1073 |
|
|
airMopAdd(mop, nbuffB, airFree, airMopAlways); |
1074 |
|
|
for (i=0; i<ninLen; i++) { |
1075 |
|
|
airMopAdd(mop, nbuffA[i] = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
1076 |
|
|
airMopAdd(mop, nbuffB[i] = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
1077 |
|
|
nrrdAxisInfoCopy(nout[i], nin[i], NULL, NRRD_AXIS_INFO_NONE); |
1078 |
|
|
} |
1079 |
|
|
airMopAdd(mop, npxfr = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
1080 |
|
|
airMopAdd(mop, nhst = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
1081 |
|
|
airMopAdd(mop, ngrad = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
1082 |
|
|
if (nrrdConvert(ngrad, _ngrad, nrrdTypeDouble)) { |
1083 |
|
|
biffMovef(TEN, NRRD, "%s: trouble converting gradients to doubles", me); |
1084 |
|
|
airMopError(mop); return 1; |
1085 |
|
|
} |
1086 |
|
|
|
1087 |
|
|
/* ------ blur */ |
1088 |
|
|
if (_tenEpiRegBlur(nbuffA, nin, ninLen, bwX, bwY, verbose)) { |
1089 |
|
|
biffAddf(TEN, "%s: trouble %s", me, (bwX || bwY) ? "blurring" : "copying"); |
1090 |
|
|
airMopError(mop); return 1; |
1091 |
|
|
} |
1092 |
|
|
if (progress && _tenEpiRegSave("regtmp-blur.nrrd", NULL, |
1093 |
|
|
nbuffA, ninLen, "blurred DWIs")) { |
1094 |
|
|
airMopError(mop); return 1; |
1095 |
|
|
} |
1096 |
|
|
|
1097 |
|
|
/* ------ threshold */ |
1098 |
|
|
if (_tenEpiRegThreshold(nbuffB, nbuffA, ninLen, |
1099 |
|
|
DWthr, verbose, progress, 1.5)) { |
1100 |
|
|
biffAddf(TEN, "%s: trouble thresholding", me); |
1101 |
|
|
airMopError(mop); return 1; |
1102 |
|
|
} |
1103 |
|
|
if (progress && _tenEpiRegSave("regtmp-thresh.nrrd", NULL, |
1104 |
|
|
nbuffB, ninLen, "thresholded DWIs")) { |
1105 |
|
|
airMopError(mop); return 1; |
1106 |
|
|
} |
1107 |
|
|
|
1108 |
|
|
/* ------ connected components */ |
1109 |
|
|
if (doCC) { |
1110 |
|
|
if (_tenEpiRegCC(nbuffB, ninLen, 1, verbose)) { |
1111 |
|
|
biffAddf(TEN, "%s: trouble doing connected components", me); |
1112 |
|
|
airMopError(mop); return 1; |
1113 |
|
|
} |
1114 |
|
|
if (progress && _tenEpiRegSave("regtmp-ccs.nrrd", NULL, |
1115 |
|
|
nbuffB, ninLen, "connected components")) { |
1116 |
|
|
airMopError(mop); return 1; |
1117 |
|
|
} |
1118 |
|
|
} |
1119 |
|
|
|
1120 |
|
|
/* ------ moments */ |
1121 |
|
|
if (_tenEpiRegMoments(nbuffA, nbuffB, ninLen, verbose)) { |
1122 |
|
|
biffAddf(TEN, "%s: trouble finding moments", me); |
1123 |
|
|
airMopError(mop); return 1; |
1124 |
|
|
} |
1125 |
|
|
if (progress && _tenEpiRegSave("regtmp-mom.nrrd", NULL, |
1126 |
|
|
nbuffA, ninLen, "moments")) { |
1127 |
|
|
airMopError(mop); return 1; |
1128 |
|
|
} |
1129 |
|
|
|
1130 |
|
|
/* ------ transforms */ |
1131 |
|
|
if (_tenEpiRegPairXforms(npxfr, nbuffA, ninLen)) { |
1132 |
|
|
biffAddf(TEN, "%s: trouble calculating transforms", me); |
1133 |
|
|
airMopError(mop); return 1; |
1134 |
|
|
} |
1135 |
|
|
if (progress && _tenEpiRegSave("regtmp-pxfr.nrrd", npxfr, |
1136 |
|
|
NULL, 0, "pair-wise xforms")) { |
1137 |
|
|
airMopError(mop); return 1; |
1138 |
|
|
} |
1139 |
|
|
|
1140 |
|
|
if (-1 == reference) { |
1141 |
|
|
/* ------ HST estimation */ |
1142 |
|
|
if (_tenEpiRegEstimHST(nhst, npxfr, ninLen, ngrad)) { |
1143 |
|
|
biffAddf(TEN, "%s: trouble estimating HST", me); |
1144 |
|
|
airMopError(mop); return 1; |
1145 |
|
|
} |
1146 |
|
|
if (progress && _tenEpiRegSave("regtmp-hst.txt", nhst, |
1147 |
|
|
NULL, 0, "HST estimates")) { |
1148 |
|
|
airMopError(mop); return 1; |
1149 |
|
|
} |
1150 |
|
|
|
1151 |
|
|
if (fitFrac) { |
1152 |
|
|
/* ------ HST parameter fitting */ |
1153 |
|
|
if (_tenEpiRegFitHST(nhst, nbuffB, ninLen, fitFrac, progress, verbose)) { |
1154 |
|
|
biffAddf(TEN, "%s: trouble fitting HST", me); |
1155 |
|
|
airMopError(mop); return 1; |
1156 |
|
|
} |
1157 |
|
|
if (progress && _tenEpiRegSave("regtmp-fit-hst.txt", nhst, |
1158 |
|
|
NULL, 0, "fitted HST")) { |
1159 |
|
|
airMopError(mop); return 1; |
1160 |
|
|
} |
1161 |
|
|
} |
1162 |
|
|
} |
1163 |
|
|
|
1164 |
|
|
/* ------ doit */ |
1165 |
|
|
if (_tenEpiRegWarp(nout, npxfr, nhst, ngrad, nin, ninLen, |
1166 |
|
|
reference, kern, kparm, verbose)) { |
1167 |
|
|
biffAddf(TEN, "%s: trouble performing final registration", me); |
1168 |
|
|
airMopError(mop); return 1; |
1169 |
|
|
} |
1170 |
|
|
|
1171 |
|
|
airMopOkay(mop); |
1172 |
|
|
nrrdStateAlwaysSetContent = hack1; |
1173 |
|
|
nrrdStateDisableContent = hack2; |
1174 |
|
|
return 0; |
1175 |
|
|
} |
1176 |
|
|
|
1177 |
|
|
int |
1178 |
|
|
tenEpiRegister4D(Nrrd *_nout, Nrrd *_nin, Nrrd *_ngrad, |
1179 |
|
|
int reference, |
1180 |
|
|
double bwX, double bwY, double fitFrac, |
1181 |
|
|
double DWthr, int doCC, |
1182 |
|
|
const NrrdKernel *kern, double *kparm, |
1183 |
|
|
int progress, int verbose) { |
1184 |
|
|
static const char me[]="tenEpiRegister4D"; |
1185 |
|
|
unsigned int ninIdx, ninLen, |
1186 |
|
|
dwiAx, rangeAxisNum, rangeAxisIdx[NRRD_DIM_MAX]; |
1187 |
|
|
int dwiIdx; |
1188 |
|
|
Nrrd **nout, **nin, **ndwi, **ndwiOut, *ngrad, *ndwigrad; |
1189 |
|
|
airArray *mop; |
1190 |
|
|
double *grad, *dwigrad, glen; |
1191 |
|
|
|
1192 |
|
|
if (!(_nout && _nin)) { |
1193 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
1194 |
|
|
return 1; |
1195 |
|
|
} |
1196 |
|
|
if (4 != _nin->dim) { |
1197 |
|
|
biffAddf(TEN, "%s: need a 4-D input array, not %d-D", me, _nin->dim); |
1198 |
|
|
return 1; |
1199 |
|
|
} |
1200 |
|
|
if (tenGradientCheck(_ngrad, nrrdTypeDefault, 6 /* <-- HEY: not sure */)) { |
1201 |
|
|
biffAddf(TEN, "%s: problem with given gradient list", me); |
1202 |
|
|
return 1; |
1203 |
|
|
} |
1204 |
|
|
|
1205 |
|
|
rangeAxisNum = nrrdRangeAxesGet(_nin, rangeAxisIdx); |
1206 |
|
|
if (0 == rangeAxisNum) { |
1207 |
|
|
/* we fall back on old behavior */ |
1208 |
|
|
dwiAx = 0; |
1209 |
|
|
} else if (1 == rangeAxisNum) { |
1210 |
|
|
/* thankfully there's exactly one range axis */ |
1211 |
|
|
dwiAx = rangeAxisIdx[0]; |
1212 |
|
|
} else { |
1213 |
|
|
biffAddf(TEN, "%s: have %u range axes instead of 1, don't know which " |
1214 |
|
|
"is DWI axis", me, rangeAxisNum); |
1215 |
|
|
return 1; |
1216 |
|
|
} |
1217 |
|
|
|
1218 |
|
|
ninLen = _nin->axis[dwiAx].size; |
1219 |
|
|
/* outdated |
1220 |
|
|
if (!( AIR_IN_CL(6, ninLen, 120) )) { |
1221 |
|
|
biffAddf(TEN, "%s: %u (size of axis %u, and # DWIs) is unreasonable", |
1222 |
|
|
me, ninLen, dwiAx); |
1223 |
|
|
return 1; |
1224 |
|
|
} |
1225 |
|
|
*/ |
1226 |
|
|
if (ninLen != _ngrad->axis[1].size) { |
1227 |
|
|
char stmp[AIR_STRLEN_SMALL]; |
1228 |
|
|
biffAddf(TEN, "%s: ninLen %u != # grads %s", me, ninLen, |
1229 |
|
|
airSprintSize_t(stmp, _ngrad->axis[1].size)); |
1230 |
|
|
return 1; |
1231 |
|
|
} |
1232 |
|
|
|
1233 |
|
|
mop = airMopNew(); |
1234 |
|
|
ngrad = nrrdNew(); |
1235 |
|
|
ndwigrad = nrrdNew(); |
1236 |
|
|
airMopAdd(mop, ngrad, (airMopper)nrrdNuke, airMopAlways); |
1237 |
|
|
airMopAdd(mop, ndwigrad, (airMopper)nrrdNuke, airMopAlways); |
1238 |
|
|
if (nrrdConvert(ngrad, _ngrad, nrrdTypeDouble) |
1239 |
|
|
|| nrrdConvert(ndwigrad, _ngrad, nrrdTypeDouble)) { /* HACK applies */ |
1240 |
|
|
biffMovef(TEN, NRRD, "%s: trouble converting gradients to doubles", me); |
1241 |
|
|
airMopError(mop); return 1; |
1242 |
|
|
} |
1243 |
|
|
|
1244 |
|
|
nin = (Nrrd **)calloc(ninLen, sizeof(Nrrd*)); |
1245 |
|
|
ndwi = (Nrrd **)calloc(ninLen, sizeof(Nrrd*)); |
1246 |
|
|
nout = (Nrrd **)calloc(ninLen, sizeof(Nrrd*)); |
1247 |
|
|
ndwiOut = (Nrrd **)calloc(ninLen, sizeof(Nrrd*)); |
1248 |
|
|
if (!(nin && ndwi && nout && ndwiOut)) { |
1249 |
|
|
biffAddf(TEN, "%s: trouble allocating local arrays", me); |
1250 |
|
|
airMopError(mop); return 1; |
1251 |
|
|
} |
1252 |
|
|
airMopAdd(mop, nin, airFree, airMopAlways); |
1253 |
|
|
airMopAdd(mop, ndwi, airFree, airMopAlways); |
1254 |
|
|
airMopAdd(mop, nout, airFree, airMopAlways); |
1255 |
|
|
airMopAdd(mop, ndwiOut, airFree, airMopAlways); |
1256 |
|
|
dwiIdx = -1; |
1257 |
|
|
grad = AIR_CAST(double *, ngrad->data); |
1258 |
|
|
dwigrad = AIR_CAST(double *, ndwigrad->data); |
1259 |
|
|
for (ninIdx=0; ninIdx<ninLen; ninIdx++) { |
1260 |
|
|
glen = ELL_3V_LEN(grad + 3*ninIdx); |
1261 |
|
|
if (-1 != reference || glen > 0.0) { |
1262 |
|
|
/* we're not allowed to use only those images that are truly DWIs, |
1263 |
|
|
or, |
1264 |
|
|
(we are so allowed and) this is a true DWI */ |
1265 |
|
|
dwiIdx++; |
1266 |
|
|
ndwi[dwiIdx] = nrrdNew(); |
1267 |
|
|
ndwiOut[dwiIdx] = nrrdNew(); |
1268 |
|
|
airMopAdd(mop, ndwi[dwiIdx], (airMopper)nrrdNuke, airMopAlways); |
1269 |
|
|
airMopAdd(mop, ndwiOut[dwiIdx], (airMopper)nrrdNuke, airMopAlways); |
1270 |
|
|
if (nrrdSlice(ndwi[dwiIdx], _nin, dwiAx, |
1271 |
|
|
AIR_CAST(unsigned int, ninIdx))) { |
1272 |
|
|
biffMovef(TEN, NRRD, "%s: trouble slicing at %d on axis %u", |
1273 |
|
|
me, ninIdx, dwiAx); |
1274 |
|
|
airMopError(mop); return 1; |
1275 |
|
|
} |
1276 |
|
|
/* NOTE: this works because dwiIdx <= ninIdx */ |
1277 |
|
|
ELL_3V_COPY(dwigrad + 3*dwiIdx, grad + 3*ninIdx); |
1278 |
|
|
} else { |
1279 |
|
|
/* we are only looking at true DWIs, and this isn't one of them */ |
1280 |
|
|
continue; |
1281 |
|
|
} |
1282 |
|
|
} |
1283 |
|
|
if (-1 == dwiIdx) { |
1284 |
|
|
biffAddf(TEN, "%s: somehow got no DWIs", me); |
1285 |
|
|
airMopError(mop); return 1; |
1286 |
|
|
} |
1287 |
|
|
/* HEY: HACK! */ |
1288 |
|
|
ndwigrad->axis[1].size = 1 + AIR_CAST(unsigned int, dwiIdx); |
1289 |
|
|
if (tenEpiRegister3D(ndwiOut, ndwi, ndwigrad->axis[1].size, ndwigrad, |
1290 |
|
|
reference, |
1291 |
|
|
bwX, bwY, fitFrac, DWthr, |
1292 |
|
|
doCC, |
1293 |
|
|
kern, kparm, |
1294 |
|
|
progress, verbose)) { |
1295 |
|
|
biffAddf(TEN, "%s: trouble", me); |
1296 |
|
|
airMopError(mop); return 1; |
1297 |
|
|
} |
1298 |
|
|
dwiIdx = -1; |
1299 |
|
|
for (ninIdx=0; ninIdx<ninLen; ninIdx++) { |
1300 |
|
|
glen = ELL_3V_LEN(grad + 3*ninIdx); |
1301 |
|
|
if (-1 != reference || glen > 0.0) { |
1302 |
|
|
/* we're not allowed to use only those images that are truly DWIs, |
1303 |
|
|
or, |
1304 |
|
|
(we are so allowed and) this is a true DWI */ |
1305 |
|
|
dwiIdx++; |
1306 |
|
|
nout[ninIdx] = ndwiOut[dwiIdx]; |
1307 |
|
|
/* |
1308 |
|
|
fprintf(stderr, "!%s: (A) nout[%u] = %p\n", me, ninIdx, nout[ninIdx]); |
1309 |
|
|
*/ |
1310 |
|
|
} else { |
1311 |
|
|
/* we are only looking at true DWIs, and this isn't one of them */ |
1312 |
|
|
nout[ninIdx] = nrrdNew(); |
1313 |
|
|
airMopAdd(mop, nout[ninIdx], (airMopper)nrrdNuke, airMopAlways); |
1314 |
|
|
if (nrrdSlice(nout[ninIdx], _nin, dwiAx, ninIdx)) { |
1315 |
|
|
biffMovef(TEN, NRRD, "%s: trouble slicing at %d on axis %u", |
1316 |
|
|
me, dwiIdx, dwiAx); |
1317 |
|
|
airMopError(mop); return 1; |
1318 |
|
|
} |
1319 |
|
|
/* |
1320 |
|
|
fprintf(stderr, "!%s: (B) nout[%u] = %p\n", me, ninIdx, nout[ninIdx]); |
1321 |
|
|
*/ |
1322 |
|
|
} |
1323 |
|
|
} |
1324 |
|
|
if (nrrdJoin(_nout, (const Nrrd*const*)nout, ninLen, dwiAx, AIR_TRUE)) { |
1325 |
|
|
biffMovef(TEN, NRRD, "%s: trouble joining output", me); |
1326 |
|
|
airMopError(mop); return 1; |
1327 |
|
|
} |
1328 |
|
|
nrrdAxisInfoCopy(_nout, _nin, NULL, NRRD_AXIS_INFO_NONE); |
1329 |
|
|
if (nrrdBasicInfoCopy(_nout, _nin, |
1330 |
|
|
NRRD_BASIC_INFO_DATA_BIT |
1331 |
|
|
| NRRD_BASIC_INFO_TYPE_BIT |
1332 |
|
|
| NRRD_BASIC_INFO_BLOCKSIZE_BIT |
1333 |
|
|
| NRRD_BASIC_INFO_DIMENSION_BIT |
1334 |
|
|
| NRRD_BASIC_INFO_CONTENT_BIT |
1335 |
|
|
| NRRD_BASIC_INFO_COMMENTS_BIT)) { |
1336 |
|
|
/* note that we're ALWAYS copying the key/value pairs- its just |
1337 |
|
|
too annoying to have to set nrrdStateKeyValuePairsPropagate |
1338 |
|
|
in order for the DWI-specific key/value pairs to be set */ |
1339 |
|
|
biffMovef(TEN, NRRD, "%s:", me); |
1340 |
|
|
airMopError(mop); return 1; |
1341 |
|
|
} |
1342 |
|
|
|
1343 |
|
|
airMopOkay(mop); |
1344 |
|
|
return 0; |
1345 |
|
|
} |