1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 */
23
24 #include "ten.h"
25 #include "privateTen.h"
26
27 #define INFO "Calculate one or more eigenvalues in a DT volume"
28 static const char *_tend_evalInfoL =
29 (INFO
30 ". ");
31
32 int
tend_evalMain(int argc,const char ** argv,const char * me,hestParm * hparm)33 tend_evalMain(int argc, const char **argv, const char *me,
34 hestParm *hparm) {
35 int pret, map[4];
36 hestOpt *hopt = NULL;
37 char *perr, *err;
38 airArray *mop;
39
40 int ret, *comp, compLen, cc;
41 Nrrd *nin, *nout;
42 char *outS;
43 float thresh, *edata, *tdata, eval[3], evec[9];
44 size_t N, I, sx, sy, sz;
45
46 hestOptAdd(&hopt, "c", "c0 ", airTypeInt, 1, 3, &comp, NULL,
47 "which eigenvalues should be saved out. \"0\" for the "
48 "largest, \"1\" for the middle, \"2\" for the smallest, "
49 "\"0 1\", \"1 2\", \"0 1 2\" or similar for more than one",
50 &compLen);
51 hestOptAdd(&hopt, "t", "thresh", airTypeFloat, 1, 1, &thresh, "0.5",
52 "confidence threshold");
53 hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, "-",
54 "input diffusion tensor volume", NULL, NULL, nrrdHestNrrd);
55 hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
56 "output image (floating point)");
57
58 mop = airMopNew();
59 airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
60 USAGE(_tend_evalInfoL);
61 PARSE();
62 airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
63
64 for (cc=0; cc<compLen; cc++) {
65 if (!AIR_IN_CL(0, comp[cc], 2)) {
66 fprintf(stderr, "%s: requested component %d (%d of 3) not in [0..2]\n",
67 me, comp[cc], cc+1);
68 airMopError(mop); return 1;
69 }
70 }
71 if (tenTensorCheck(nin, nrrdTypeFloat, AIR_TRUE, AIR_TRUE)) {
72 airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
73 fprintf(stderr, "%s: didn't get a valid DT volume:\n%s\n", me, err);
74 airMopError(mop); return 1;
75 }
76
77 sx = nin->axis[1].size;
78 sy = nin->axis[2].size;
79 sz = nin->axis[3].size;
80
81 nout = nrrdNew();
82 airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
83 if (1 == compLen) {
84 ret = nrrdMaybeAlloc_va(nout, nrrdTypeFloat, 3,
85 sx, sy, sz);
86 } else {
87 ret = nrrdMaybeAlloc_va(nout, nrrdTypeFloat, 4,
88 AIR_CAST(size_t, compLen), sx, sy, sz);
89 }
90 if (ret) {
91 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
92 fprintf(stderr, "%s: trouble allocating output:\n%s\n", me, err);
93 airMopError(mop); return 1;
94 }
95
96 N = sx*sy*sz;
97 edata = (float *)nout->data;
98 tdata = (float *)nin->data;
99 if (1 == compLen) {
100 ELL_3V_SET(map, 1, 2, 3);
101 for (I=0; I<N; I++) {
102 tenEigensolve_f(eval, evec, tdata);
103 edata[I] = (tdata[0] >= thresh)*eval[comp[0]];
104 tdata += 7;
105 }
106 } else {
107 ELL_4V_SET(map, 0, 1, 2, 3);
108 for (I=0; I<N; I++) {
109 tenEigensolve_f(eval, evec, tdata);
110 for (cc=0; cc<compLen; cc++)
111 edata[cc] = (tdata[0] >= thresh)*eval[comp[cc]];
112 edata += compLen;
113 tdata += 7;
114 }
115 }
116 if (nrrdAxisInfoCopy(nout, nin, map, NRRD_AXIS_INFO_SIZE_BIT)) {
117 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
118 fprintf(stderr, "%s: trouble:\n%s\n", me, err);
119 airMopError(mop); return 1;
120 }
121 if (nrrdBasicInfoCopy(nout, nin,
122 NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) {
123 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
124 fprintf(stderr, "%s: trouble:\n%s\n", me, err);
125 airMopError(mop); return 1;
126 }
127 if (1 != compLen) {
128 nout->axis[0].label = (char *)airFree(nout->axis[0].label);
129 nout->axis[0].kind = nrrdKindUnknown;
130 }
131
132 if (nrrdSave(outS, nout, NULL)) {
133 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
134 fprintf(stderr, "%s: trouble writing:\n%s\n", me, err);
135 airMopError(mop); return 1;
136 }
137
138 airMopOkay(mop);
139 return 0;
140 }
141 TEND_CMD(eval, INFO);
142