GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/ten/tendFiber.c Lines: 36 154 23.4 %
Date: 2017-05-26 Branches: 1 102 1.0 %

Line Branch Exec Source
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);