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 eigenvectors in a DT volume"
28 static const char *_tend_evecInfoL =
29 (INFO
30 ". ");
31
32 int
tend_evecMain(int argc,const char ** argv,const char * me,hestParm * hparm)33 tend_evecMain(int argc, const char **argv, const char *me,
34 hestParm *hparm) {
35 int pret;
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], scl;
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_evecInfoL);
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 ret = nrrdMaybeAlloc_va(nout, nrrdTypeFloat, 4,
84 AIR_CAST(size_t, 3*compLen), sx, sy, sz);
85 if (ret) {
86 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
87 fprintf(stderr, "%s: trouble allocating output:\n%s\n", me, err);
88 airMopError(mop); return 1;
89 }
90
91 N = sx*sy*sz;
92 edata = (float *)nout->data;
93 tdata = (float *)nin->data;
94 if (1 == compLen) {
95 for (I=0; I<N; I++) {
96 tenEigensolve_f(eval, evec, tdata);
97 scl = AIR_CAST(float, tdata[0] >= thresh);
98 ELL_3V_SCALE(edata, scl, evec+3*comp[0]);
99 edata += 3;
100 tdata += 7;
101 }
102 } else {
103 for (I=0; I<N; I++) {
104 tenEigensolve_f(eval, evec, tdata);
105 scl = AIR_CAST(float, tdata[0] >= thresh);
106 for (cc=0; cc<compLen; cc++) {
107 ELL_3V_SCALE(edata+3*cc, scl, evec+3*comp[cc]);
108 }
109 edata += 3*compLen;
110 tdata += 7;
111 }
112 }
113 if (nrrdAxisInfoCopy(nout, nin, NULL, NRRD_AXIS_INFO_SIZE_BIT)) {
114 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
115 fprintf(stderr, "%s: trouble:\n%s\n", me, err);
116 airMopError(mop); return 1;
117 }
118 if (nrrdBasicInfoCopy(nout, nin,
119 NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) {
120 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
121 fprintf(stderr, "%s: trouble:\n%s\n", me, err);
122 airMopError(mop); return 1;
123 }
124 nout->axis[0].label = (char *)airFree(nout->axis[0].label);
125 nout->axis[0].kind = nrrdKindUnknown;
126
127 if (nrrdSave(outS, nout, NULL)) {
128 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
129 fprintf(stderr, "%s: trouble writing:\n%s\n", me, err);
130 airMopError(mop); return 1;
131 }
132
133 airMopOkay(mop);
134 return 0;
135 }
136 TEND_CMD(evec, INFO);
137