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 "Register diffusion-weighted echo-planar images" |
28 |
|
|
static const char *_tend_epiregInfoL = |
29 |
|
|
(INFO |
30 |
|
|
". This registration corrects the shear, scale, and translate along " |
31 |
|
|
"the phase encoding direction (assumed to be the Y (second) axis of " |
32 |
|
|
"the image) caused by eddy currents from the diffusion-encoding " |
33 |
|
|
"gradients with echo-planar imaging. The method is based on calculating " |
34 |
|
|
"moments of segmented images, where the segmentation is a simple " |
35 |
|
|
"procedure based on blurring (optional), thresholding and " |
36 |
|
|
"connected component analysis. " |
37 |
|
|
"The registered DWIs are resampled with the " |
38 |
|
|
"chosen kernel, with the separate DWIs stacked along axis 0."); |
39 |
|
|
|
40 |
|
|
int |
41 |
|
|
tend_epiregMain(int argc, const char **argv, const char *me, |
42 |
|
|
hestParm *hparm) { |
43 |
|
|
int pret, rret; |
44 |
|
2 |
hestOpt *hopt = NULL; |
45 |
|
1 |
char *perr, *err; |
46 |
|
|
airArray *mop; |
47 |
|
1 |
char *outS, *buff; |
48 |
|
|
|
49 |
|
1 |
char *gradS; |
50 |
|
1 |
NrrdKernelSpec *ksp; |
51 |
|
1 |
Nrrd **nin, **nout3D, *nout4D, *ngrad, *ngradKVP, *nbmatKVP; |
52 |
|
1 |
unsigned int ni, ninLen, *skip, skipNum; |
53 |
|
1 |
int ref, noverbose, progress, nocc, baseNum; |
54 |
|
1 |
float bw[2], thr, fitFrac; |
55 |
|
1 |
double bvalue; |
56 |
|
|
|
57 |
|
1 |
hestOptAdd(&hopt, "i", "dwi0 dwi1", airTypeOther, 1, -1, &nin, NULL, |
58 |
|
|
"all the diffusion-weighted images (DWIs), as separate 3D nrrds, " |
59 |
|
|
"**OR**: one 4D nrrd of all DWIs stacked along axis 0", |
60 |
|
1 |
&ninLen, NULL, nrrdHestNrrd); |
61 |
|
1 |
hestOptAdd(&hopt, "g", "grads", airTypeString, 1, 1, &gradS, NULL, |
62 |
|
|
"array of gradient directions, in the same order as the " |
63 |
|
|
"associated DWIs were given to \"-i\", " |
64 |
|
|
"**OR** \"-g kvp\" signifies that gradient directions should " |
65 |
|
|
"be read from the key/value pairs of the DWI", |
66 |
|
1 |
NULL, NULL, nrrdHestNrrd); |
67 |
|
1 |
hestOptAdd(&hopt, "r", "reference", airTypeInt, 1, 1, &ref, "-1", |
68 |
|
|
"which of the DW volumes (zero-based numbering) should be used " |
69 |
|
|
"as the standard, to which all other images are transformed. " |
70 |
|
|
"Using -1 (the default) means that 9 intrinsic parameters " |
71 |
|
|
"governing the relationship between the gradient direction " |
72 |
|
|
"and the resulting distortion are estimated and fitted, " |
73 |
|
|
"ensuring good registration with the non-diffusion-weighted " |
74 |
|
|
"T2 image (which is never explicitly used in registration). " |
75 |
|
|
"Otherwise, by picking a specific DWI, no distortion parameter " |
76 |
|
|
"estimation is done. "); |
77 |
|
1 |
hestOptAdd(&hopt, "nv", NULL, airTypeInt, 0, 0, &noverbose, NULL, |
78 |
|
|
"turn OFF verbose mode, and " |
79 |
|
|
"have no idea what stage processing is at."); |
80 |
|
1 |
hestOptAdd(&hopt, "p", NULL, airTypeInt, 0, 0, &progress, NULL, |
81 |
|
|
"save out intermediate steps of processing"); |
82 |
|
1 |
hestOptAdd(&hopt, "bw", "x,y blur", airTypeFloat, 2, 2, bw, "1.0 2.0", |
83 |
|
|
"standard devs in X and Y directions of gaussian filter used " |
84 |
|
|
"to blur the DWIs prior to doing segmentation. This blurring " |
85 |
|
|
"does not effect the final resampling of registered DWIs. " |
86 |
|
|
"Use \"0.0 0.0\" to say \"no blurring\""); |
87 |
|
1 |
hestOptAdd(&hopt, "t", "DWI thresh", airTypeFloat, 1, 1, &thr, "nan", |
88 |
|
|
"Threshold value to use on DWIs, " |
89 |
|
|
"to do initial separation of brain and non-brain. By default, " |
90 |
|
|
"the threshold is determined automatically by histogram " |
91 |
|
|
"analysis. "); |
92 |
|
1 |
hestOptAdd(&hopt, "ncc", NULL, airTypeInt, 0, 0, &nocc, NULL, |
93 |
|
|
"do *NOT* do connected component (CC) analysis, after " |
94 |
|
|
"thresholding and before moment calculation. Doing CC analysis " |
95 |
|
|
"usually gives better results because it converts the " |
96 |
|
|
"thresholding output into something much closer to a " |
97 |
|
|
"real segmentation"); |
98 |
|
1 |
hestOptAdd(&hopt, "f", "fit frac", airTypeFloat, 1, 1, &fitFrac, "0.70", |
99 |
|
|
"(only meaningful with \"-r -1\") When doing linear fitting " |
100 |
|
|
"of the intrinsic distortion parameters, it is good " |
101 |
|
|
"to ignore the slices for which the segmentation was poor. A " |
102 |
|
|
"heuristic is used to rank the slices according to segmentation " |
103 |
|
|
"quality. This option controls how many of the (best) slices " |
104 |
|
|
"contribute to the fitting. Use \"0\" to disable distortion " |
105 |
|
|
"parameter fitting. "); |
106 |
|
1 |
hestOptAdd(&hopt, "k", "kernel", airTypeOther, 1, 1, &ksp, "cubic:0,0.5", |
107 |
|
|
"kernel for resampling DWIs along the phase-encoding " |
108 |
|
|
"direction during final registration stage", |
109 |
|
1 |
NULL, NULL, nrrdHestKernelSpec); |
110 |
|
1 |
hestOptAdd(&hopt, "s", "start #", airTypeInt, 1, 1, &baseNum, "1", |
111 |
|
|
"first number to use in numbered sequence of output files."); |
112 |
|
1 |
hestOptAdd(&hopt, "o", "output/prefix", airTypeString, 1, 1, &outS, "-", |
113 |
|
|
"For separate 3D DWI volume inputs: prefix for output filenames; " |
114 |
|
|
"will save out one (registered) " |
115 |
|
|
"DWI for each input DWI, using the same type as the input. " |
116 |
|
|
"**OR**: For single 4D DWI input: output file name. "); |
117 |
|
|
|
118 |
|
1 |
mop = airMopNew(); |
119 |
|
1 |
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); |
120 |
✓✗ |
2 |
USAGE(_tend_epiregInfoL); |
121 |
|
|
JUSTPARSE(); |
122 |
|
|
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); |
123 |
|
|
|
124 |
|
|
if (strcmp("kvp", gradS)) { |
125 |
|
|
/* they're NOT coming from key/value pairs */ |
126 |
|
|
if (nrrdLoad(ngrad=nrrdNew(), gradS, NULL)) { |
127 |
|
|
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); |
128 |
|
|
fprintf(stderr, "%s: trouble loading gradient list:\n%s\n", me, err); |
129 |
|
|
airMopError(mop); return 1; |
130 |
|
|
} |
131 |
|
|
} else { |
132 |
|
|
if (1 != ninLen) { |
133 |
|
|
fprintf(stderr, "%s: can do key/value pairs only from single nrrd", me); |
134 |
|
|
airMopError(mop); return 1; |
135 |
|
|
} |
136 |
|
|
/* they are coming from key/value pairs */ |
137 |
|
|
if (tenDWMRIKeyValueParse(&ngradKVP, &nbmatKVP, &bvalue, |
138 |
|
|
&skip, &skipNum, nin[0])) { |
139 |
|
|
airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); |
140 |
|
|
fprintf(stderr, "%s: trouble parsing gradient list:\n%s\n", me, err); |
141 |
|
|
airMopError(mop); return 1; |
142 |
|
|
} |
143 |
|
|
if (nbmatKVP) { |
144 |
|
|
fprintf(stderr, "%s: sorry, can only use gradients, not b-matrices", me); |
145 |
|
|
airMopError(mop); return 1; |
146 |
|
|
} |
147 |
|
|
ngrad = ngradKVP; |
148 |
|
|
} |
149 |
|
|
airMopAdd(mop, ngrad, (airMopper)nrrdNuke, airMopAlways); |
150 |
|
|
|
151 |
|
|
nout3D = AIR_CALLOC(ninLen, Nrrd *); |
152 |
|
|
airMopAdd(mop, nout3D, airFree, airMopAlways); |
153 |
|
|
nout4D = nrrdNew(); |
154 |
|
|
airMopAdd(mop, nout4D, (airMopper)nrrdNuke, airMopAlways); |
155 |
|
|
buff = AIR_CALLOC(airStrlen(outS) + 10, char); |
156 |
|
|
airMopAdd(mop, buff, airFree, airMopAlways); |
157 |
|
|
if (!( nout3D && nout4D && buff )) { |
158 |
|
|
fprintf(stderr, "%s: couldn't allocate buffers", me); |
159 |
|
|
airMopError(mop); return 1; |
160 |
|
|
} |
161 |
|
|
for (ni=0; ni<ninLen; ni++) { |
162 |
|
|
nout3D[ni]=nrrdNew(); |
163 |
|
|
airMopAdd(mop, nout3D[ni], (airMopper)nrrdNuke, airMopAlways); |
164 |
|
|
} |
165 |
|
|
if (1 == ninLen) { |
166 |
|
|
rret = tenEpiRegister4D(nout4D, nin[0], ngrad, |
167 |
|
|
ref, |
168 |
|
|
bw[0], bw[1], fitFrac, thr, !nocc, |
169 |
|
|
ksp->kernel, ksp->parm, |
170 |
|
|
progress, !noverbose); |
171 |
|
|
} else { |
172 |
|
|
rret = tenEpiRegister3D(nout3D, nin, ninLen, ngrad, |
173 |
|
|
ref, |
174 |
|
|
bw[0], bw[1], fitFrac, thr, !nocc, |
175 |
|
|
ksp->kernel, ksp->parm, |
176 |
|
|
progress, !noverbose); |
177 |
|
|
} |
178 |
|
|
if (rret) { |
179 |
|
|
airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); |
180 |
|
|
fprintf(stderr, "%s: trouble doing epireg:\n%s\n", me, err); |
181 |
|
|
airMopError(mop); return 1; |
182 |
|
|
} |
183 |
|
|
|
184 |
|
|
if (1 == ninLen) { |
185 |
|
|
if (nrrdSave(outS, nout4D, NULL)) { |
186 |
|
|
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); |
187 |
|
|
fprintf(stderr, "%s: trouble writing \"%s\":\n%s\n", me, outS, err); |
188 |
|
|
airMopError(mop); return 1; |
189 |
|
|
} |
190 |
|
|
} else { |
191 |
|
|
for (ni=0; ni<ninLen; ni++) { |
192 |
|
|
if (ninLen+baseNum > 99) { |
193 |
|
|
sprintf(buff, "%s%05d.nrrd", outS, ni+baseNum); |
194 |
|
|
} else if (ninLen+baseNum > 9) { |
195 |
|
|
sprintf(buff, "%s%02d.nrrd", outS, ni+baseNum); |
196 |
|
|
} else { |
197 |
|
|
sprintf(buff, "%s%d.nrrd", outS, ni+baseNum); |
198 |
|
|
} |
199 |
|
|
if (nrrdSave(buff, nout3D[ni], NULL)) { |
200 |
|
|
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); |
201 |
|
|
fprintf(stderr, "%s: trouble writing \"%s\":\n%s\n", me, buff, err); |
202 |
|
|
airMopError(mop); return 1; |
203 |
|
|
} |
204 |
|
|
} |
205 |
|
|
} |
206 |
|
|
|
207 |
|
|
airMopOkay(mop); |
208 |
|
|
return 0; |
209 |
|
1 |
} |
210 |
|
|
TEND_CMD(epireg, INFO); |