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 |
|
|
#define INFO "Fiber tractography, from one or more seeds" |
28 |
|
|
static const char *_tend_fiberInfoL = |
29 |
|
|
(INFO |
30 |
|
|
". A fairly complete command-line interface to the tenFiber API."); |
31 |
|
|
|
32 |
|
|
int |
33 |
|
|
tend_fiberMain(int argc, const char **argv, const char *me, |
34 |
|
|
hestParm *hparm) { |
35 |
|
|
int pret; |
36 |
|
2 |
hestOpt *hopt = NULL; |
37 |
|
1 |
char *perr, *err; |
38 |
|
|
airArray *mop; |
39 |
|
1 |
char *outS; |
40 |
|
|
|
41 |
|
|
tenFiberContext *tfx; |
42 |
|
|
tenFiberSingle *tfbs; |
43 |
|
1 |
NrrdKernelSpec *ksp; |
44 |
|
1 |
double start[3], step, *_stop, *stop; |
45 |
|
|
const airEnum *ftypeEnum; |
46 |
|
1 |
char *ftypeS; |
47 |
|
1 |
int E, intg, useDwi, allPaths, verbose, worldSpace, worldSpaceOut, |
48 |
|
|
ftype, ftypeDef; |
49 |
|
1 |
Nrrd *nin, *nseed, *nmat, *_nmat; |
50 |
|
1 |
unsigned int si, stopLen, whichPath; |
51 |
|
1 |
double matx[16]={1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1}; |
52 |
|
|
tenFiberMulti *tfml; |
53 |
|
|
limnPolyData *fiberPld; |
54 |
|
|
|
55 |
|
1 |
hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, "-", |
56 |
|
1 |
"input volume", NULL, NULL, nrrdHestNrrd); |
57 |
|
1 |
hestOptAdd(&hopt, "dwi", NULL, airTypeInt, 0, 0, &useDwi, NULL, |
58 |
|
|
"input volume is a DWI volume, not a single tensor volume"); |
59 |
|
1 |
hestOptAdd(&hopt, "s", "seed point", airTypeDouble, 3, 3, start, "0 0 0", |
60 |
|
|
"seed point for fiber; it will propogate in two opposite " |
61 |
|
|
"directions starting from here"); |
62 |
|
1 |
hestOptAdd(&hopt, "ns", "seed nrrd", airTypeOther, 1, 1, &nseed, "", |
63 |
|
1 |
"3-by-N nrrd of seedpoints", NULL, NULL, nrrdHestNrrd); |
64 |
|
1 |
hestOptAdd(&hopt, "wsp", NULL, airTypeInt, 0, 0, &worldSpace, NULL, |
65 |
|
|
"define seedpoint and output path in worldspace. Otherwise, " |
66 |
|
|
"(without using this option) everything is in index space"); |
67 |
|
1 |
hestOptAdd(&hopt, "t", "type", airTypeString, 1, 1, &ftypeS, "", |
68 |
|
|
"fiber type; defaults to something"); |
69 |
|
1 |
hestOptAdd(&hopt, "n", "intg", airTypeEnum, 1, 1, &intg, "rk4", |
70 |
|
|
"integration method for fiber tracking", |
71 |
|
1 |
NULL, tenFiberIntg); |
72 |
|
1 |
hestOptAdd(&hopt, "k", "kernel", airTypeOther, 1, 1, &ksp, "tent", |
73 |
|
|
"kernel for reconstructing tensor field", |
74 |
|
1 |
NULL, NULL, nrrdHestKernelSpec); |
75 |
|
1 |
hestOptAdd(&hopt, "wp", "which", airTypeUInt, 1, 1, &whichPath, "0", |
76 |
|
|
"when doing multi-tensor tracking, index of path to follow " |
77 |
|
|
"(made moot by \"-ap\")"); |
78 |
|
1 |
hestOptAdd(&hopt, "ap", "allpaths", airTypeInt, 0, 0, &allPaths, NULL, |
79 |
|
|
"follow all paths from (all) seedpoint(s), output will be " |
80 |
|
|
"polydata, rather than a single 3-by-N nrrd, even if only " |
81 |
|
|
"a single path is generated"); |
82 |
|
1 |
hestOptAdd(&hopt, "wspo", NULL, airTypeInt, 0, 0, &worldSpaceOut, NULL, |
83 |
|
|
"output should be in worldspace, even if input is not " |
84 |
|
|
"(this feature is unstable and/or confusing)"); |
85 |
|
1 |
hestOptAdd(&hopt, "step", "step size", airTypeDouble, 1, 1, &step, "0.01", |
86 |
|
|
"stepsize along fiber, in world space"); |
87 |
|
1 |
hestOptAdd(&hopt, "stop", "stop1", airTypeOther, 1, -1, &_stop, NULL, |
88 |
|
|
"the conditions that should signify the end of a fiber, or " |
89 |
|
|
"when to discard a fiber that is done propagating. " |
90 |
|
|
"Multiple stopping criteria are logically OR-ed and tested at " |
91 |
|
|
"every point along the fiber. Possibilities include:\n " |
92 |
|
|
"\b\bo \"aniso:<type>,<thresh>\": require anisotropy to be " |
93 |
|
|
"above the given threshold. Which anisotropy type is given " |
94 |
|
|
"as with \"tend anvol\" (see its usage info)\n " |
95 |
|
|
"\b\bo \"len:<length>\": limits the length, in world space, " |
96 |
|
|
"of each fiber half\n " |
97 |
|
|
"\b\bo \"steps:<N>\": the number of steps in each fiber half " |
98 |
|
|
"is capped at N\n " |
99 |
|
|
"\b\bo \"conf:<thresh>\": requires the tensor confidence value " |
100 |
|
|
"to be above the given thresh\n " |
101 |
|
|
"\b\bo \"radius:<thresh>\": requires that the radius of " |
102 |
|
|
"curvature of the fiber stay above given thr\n " |
103 |
|
|
"\b\bo \"frac:<F>\": in multi-tensor tracking, the fraction " |
104 |
|
|
"of the tracked component must stay above F\n " |
105 |
|
|
"\b\bo \"minlen:<len>\": discard fiber if its final whole " |
106 |
|
|
"length is below len (not really a termination criterion)\n " |
107 |
|
|
"\b\bo \"minsteps:<N>\": discard fiber if its final number of " |
108 |
|
|
"steps is below N (not really a termination criterion)", |
109 |
|
1 |
&stopLen, NULL, tendFiberStopCB); |
110 |
|
1 |
hestOptAdd(&hopt, "v", "verbose", airTypeInt, 1, 1, &verbose, "0", |
111 |
|
|
"verbosity level"); |
112 |
|
1 |
hestOptAdd(&hopt, "nmat", "transform", airTypeOther, 1, 1, &_nmat, "", |
113 |
|
|
"a 4x4 homogenous transform matrix (as a nrrd, or just a text " |
114 |
|
|
"file) given with this option will be applied to the output " |
115 |
|
|
"tractography vertices just prior to output", |
116 |
|
1 |
NULL, NULL, nrrdHestNrrd); |
117 |
|
1 |
hestOptAdd(&hopt, "o", "out", airTypeString, 1, 1, &outS, "-", |
118 |
|
|
"output fiber(s)"); |
119 |
|
|
|
120 |
|
1 |
mop = airMopNew(); |
121 |
|
1 |
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); |
122 |
✓✗ |
2 |
USAGE(_tend_fiberInfoL); |
123 |
|
|
PARSE(); |
124 |
|
|
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); |
125 |
|
|
|
126 |
|
|
tfbs = tenFiberSingleNew(); |
127 |
|
|
airMopAdd(mop, tfbs, (airMopper)tenFiberSingleNix, airMopAlways); |
128 |
|
|
|
129 |
|
|
if (_nmat) { |
130 |
|
|
if (!(2 == _nmat->dim |
131 |
|
|
&& 4 == _nmat->axis[0].size |
132 |
|
|
&& 4 == _nmat->axis[1].size)) { |
133 |
|
|
fprintf(stderr, "%s: transform matrix must be 2-D 4-by-4 array " |
134 |
|
|
"not %u-D %u-by-?\n", me, |
135 |
|
|
AIR_CAST(unsigned int, _nmat->dim), |
136 |
|
|
AIR_CAST(unsigned int, _nmat->axis[0].size)); |
137 |
|
|
airMopError(mop); return 1; |
138 |
|
|
} |
139 |
|
|
nmat = nrrdNew(); |
140 |
|
|
airMopAdd(mop, nmat, (airMopper)nrrdNuke, airMopAlways); |
141 |
|
|
if (nrrdConvert(nmat, _nmat, nrrdTypeDouble)) { |
142 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
143 |
|
|
fprintf(stderr, "%s: problem with transform matrix\n%s\n", me, err); |
144 |
|
|
airMopError(mop); return 1; |
145 |
|
|
} |
146 |
|
|
ELL_4M_COPY(matx, AIR_CAST(double *, nmat->data)); |
147 |
|
|
fprintf(stderr, "%s: transforming output by:\n", me); |
148 |
|
|
ell_4m_print_d(stderr, matx); |
149 |
|
|
} |
150 |
|
|
if (useDwi) { |
151 |
|
|
tfx = tenFiberContextDwiNew(nin, 50, 1, 1, |
152 |
|
|
tenEstimate1MethodLLS, |
153 |
|
|
tenEstimate2MethodPeled); |
154 |
|
|
} else { |
155 |
|
|
tfx = tenFiberContextNew(nin); |
156 |
|
|
} |
157 |
|
|
if (!tfx) { |
158 |
|
|
airMopAdd(mop, err = biffGetDone(TEN), airFree, airMopAlways); |
159 |
|
|
fprintf(stderr, "%s: failed to create the fiber context:\n%s\n", me, err); |
160 |
|
|
airMopError(mop); return 1; |
161 |
|
|
} |
162 |
|
|
airMopAdd(mop, tfx, (airMopper)tenFiberContextNix, airMopAlways); |
163 |
|
|
E = 0; |
164 |
|
|
for (si=0, stop=_stop; si<stopLen; si++, stop+=3) { |
165 |
|
|
int istop; /* not from Apple */ |
166 |
|
|
istop = AIR_CAST(int, stop[0]); |
167 |
|
|
switch(istop) { |
168 |
|
|
case tenFiberStopAniso: |
169 |
|
|
if (!E) E |= tenFiberStopSet(tfx, istop, |
170 |
|
|
AIR_CAST(int, stop[1]), stop[2]); |
171 |
|
|
break; |
172 |
|
|
case tenFiberStopNumSteps: |
173 |
|
|
case tenFiberStopMinNumSteps: |
174 |
|
|
if (!E) E |= tenFiberStopSet(tfx, istop, |
175 |
|
|
AIR_CAST(unsigned int, stop[1])); |
176 |
|
|
break; |
177 |
|
|
case tenFiberStopLength: |
178 |
|
|
case tenFiberStopConfidence: |
179 |
|
|
case tenFiberStopFraction: |
180 |
|
|
case tenFiberStopRadius: |
181 |
|
|
case tenFiberStopMinLength: |
182 |
|
|
if (!E) E |= tenFiberStopSet(tfx, istop, stop[1]); |
183 |
|
|
break; |
184 |
|
|
case tenFiberStopBounds: |
185 |
|
|
/* nothing to actually do */ |
186 |
|
|
break; |
187 |
|
|
default: |
188 |
|
|
fprintf(stderr, "%s: stop method %d not supported\n", me, |
189 |
|
|
istop); |
190 |
|
|
airMopError(mop); return 1; |
191 |
|
|
break; |
192 |
|
|
} |
193 |
|
|
} |
194 |
|
|
if (!E) { |
195 |
|
|
if (useDwi) { |
196 |
|
|
ftypeEnum = tenDwiFiberType; |
197 |
|
|
ftypeDef = tenDwiFiberType2Evec0; |
198 |
|
|
} else { |
199 |
|
|
ftypeEnum = tenFiberType; |
200 |
|
|
ftypeDef = tenFiberTypeEvec0; |
201 |
|
|
} |
202 |
|
|
if (airStrlen(ftypeS)) { |
203 |
|
|
ftype = airEnumVal(ftypeEnum, ftypeS); |
204 |
|
|
if (airEnumUnknown(ftypeEnum) == ftype) { |
205 |
|
|
fprintf(stderr, "%s: couldn't parse \"%s\" as a %s\n", me, |
206 |
|
|
ftypeS, ftypeEnum->name); |
207 |
|
|
airMopError(mop); return 1; |
208 |
|
|
} |
209 |
|
|
} else { |
210 |
|
|
ftype = ftypeDef; |
211 |
|
|
fprintf(stderr, "%s: (defaulting %s to %s)\n", me, |
212 |
|
|
ftypeEnum->name, airEnumStr(ftypeEnum, ftype)); |
213 |
|
|
} |
214 |
|
|
E |= tenFiberTypeSet(tfx, ftype); |
215 |
|
|
} |
216 |
|
|
if (!E) E |= tenFiberKernelSet(tfx, ksp->kernel, ksp->parm); |
217 |
|
|
if (!E) E |= tenFiberIntgSet(tfx, intg); |
218 |
|
|
if (!E) E |= tenFiberParmSet(tfx, tenFiberParmStepSize, step); |
219 |
|
|
if (!E) E |= tenFiberParmSet(tfx, tenFiberParmUseIndexSpace, |
220 |
|
|
worldSpace ? AIR_FALSE: AIR_TRUE); |
221 |
|
|
if (!E) E |= tenFiberUpdate(tfx); |
222 |
|
|
if (E) { |
223 |
|
|
airMopAdd(mop, err = biffGetDone(TEN), airFree, airMopAlways); |
224 |
|
|
fprintf(stderr, "%s: trouble:\n%s\n", me, err); |
225 |
|
|
airMopError(mop); return 1; |
226 |
|
|
} |
227 |
|
|
|
228 |
|
|
tenFiberVerboseSet(tfx, verbose); |
229 |
|
|
if (!allPaths) { |
230 |
|
|
if (tenFiberSingleTrace(tfx, tfbs, start, whichPath)) { |
231 |
|
|
airMopAdd(mop, err = biffGetDone(TEN), airFree, airMopAlways); |
232 |
|
|
fprintf(stderr, "%s: trouble:\n%s\n", me, err); |
233 |
|
|
airMopError(mop); return 1; |
234 |
|
|
} |
235 |
|
|
if (tenFiberStopUnknown == tfx->whyNowhere) { |
236 |
|
|
fprintf(stderr, "%s: steps[back,forw] = %u,%u; seedIdx = %u\n", me, |
237 |
|
|
tfbs->stepNum[0], tfbs->stepNum[1], tfbs->seedIdx); |
238 |
|
|
fprintf(stderr, "%s: whyStop[back,forw] = %s,%s\n", me, |
239 |
|
|
airEnumStr(tenFiberStop, tfbs->whyStop[0]), |
240 |
|
|
airEnumStr(tenFiberStop, tfbs->whyStop[1])); |
241 |
|
|
if (worldSpaceOut && !worldSpace) { |
242 |
|
|
/* have to convert output to worldspace */ |
243 |
|
|
fprintf(stderr, "%s: WARNING!!! output conversion " |
244 |
|
|
"to worldspace not done!!!\n", me); |
245 |
|
|
} |
246 |
|
|
if (_nmat) { |
247 |
|
|
fprintf(stderr, "%s: WARNING!!! output transform " |
248 |
|
|
"not done!!!\n", me); |
249 |
|
|
} |
250 |
|
|
if (nrrdSave(outS, tfbs->nvert, NULL)) { |
251 |
|
|
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); |
252 |
|
|
fprintf(stderr, "%s: trouble writing:\n%s\n", me, err); |
253 |
|
|
airMopError(mop); return 1; |
254 |
|
|
} |
255 |
|
|
} else { |
256 |
|
|
fprintf(stderr, "%s: fiber failed to start: %s.\n", |
257 |
|
|
me, airEnumDesc(tenFiberStop, tfx->whyNowhere)); |
258 |
|
|
} |
259 |
|
|
} else { |
260 |
|
|
if (!nseed) { |
261 |
|
|
fprintf(stderr, "%s: didn't get seed nrrd via \"-ns\"\n", me); |
262 |
|
|
airMopError(mop); return 1; |
263 |
|
|
} |
264 |
|
|
tfml = tenFiberMultiNew(); |
265 |
|
|
airMopAdd(mop, tfml, (airMopper)tenFiberMultiNix, airMopAlways); |
266 |
|
|
/* |
267 |
|
|
fiberArr = airArrayNew(AIR_CAST(void **, &fiber), NULL, |
268 |
|
|
sizeof(tenFiberSingle), 1024); |
269 |
|
|
airArrayStructCB(fiberArr, |
270 |
|
|
AIR_CAST(void (*)(void *), tenFiberSingleInit), |
271 |
|
|
AIR_CAST(void (*)(void *), tenFiberSingleDone)); |
272 |
|
|
airMopAdd(mop, fiberArr, (airMopper)airArrayNuke, airMopAlways); |
273 |
|
|
*/ |
274 |
|
|
fiberPld = limnPolyDataNew(); |
275 |
|
|
airMopAdd(mop, fiberPld, (airMopper)limnPolyDataNix, airMopAlways); |
276 |
|
|
if (tenFiberMultiTrace(tfx, tfml, nseed) |
277 |
|
|
|| tenFiberMultiPolyData(tfx, fiberPld, tfml)) { |
278 |
|
|
airMopAdd(mop, err = biffGetDone(TEN), airFree, airMopAlways); |
279 |
|
|
fprintf(stderr, "%s: trouble:\n%s\n", me, err); |
280 |
|
|
airMopError(mop); return 1; |
281 |
|
|
} |
282 |
|
|
if (worldSpaceOut && !worldSpace) { |
283 |
|
|
/* have to convert output to worldspace */ |
284 |
|
|
unsigned int ii; |
285 |
|
|
double indx[4], world[3]; |
286 |
|
|
for (ii=0; ii<fiberPld->xyzwNum; ii++) { |
287 |
|
|
ELL_4V_COPY(indx, fiberPld->xyzw + 4*ii); |
288 |
|
|
ELL_4V_HOMOG(indx, indx); |
289 |
|
|
gageShapeItoW(tfx->gtx->shape, world, indx); |
290 |
|
|
ELL_3V_COPY_TT(fiberPld->xyzw + 4*ii, float, world); |
291 |
|
|
(fiberPld->xyzw + 4*ii)[3] = 1.0; |
292 |
|
|
} |
293 |
|
|
} |
294 |
|
|
if (_nmat) { |
295 |
|
|
/* have to further transform output by given matrix */ |
296 |
|
|
unsigned int ii; |
297 |
|
|
double xxx[4], yyy[4]; |
298 |
|
|
for (ii=0; ii<fiberPld->xyzwNum; ii++) { |
299 |
|
|
ELL_4V_COPY(xxx, fiberPld->xyzw + 4*ii); |
300 |
|
|
ELL_4MV_MUL(yyy, matx, xxx); |
301 |
|
|
ELL_4V_HOMOG(yyy, yyy); |
302 |
|
|
ELL_4V_COPY_TT(fiberPld->xyzw + 4*ii, float, yyy); |
303 |
|
|
} |
304 |
|
|
} |
305 |
|
|
if (limnPolyDataSave(outS, fiberPld)) { |
306 |
|
|
airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways); |
307 |
|
|
fprintf(stderr, "%s: trouble saving:\n%s\n", me, err); |
308 |
|
|
airMopError(mop); return 1; |
309 |
|
|
} |
310 |
|
|
} |
311 |
|
|
|
312 |
|
|
airMopOkay(mop); |
313 |
|
|
return 0; |
314 |
|
1 |
} |
315 |
|
|
TEND_CMD(fiber, INFO); |