1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 /* bodily include the relevant interpolation and warp functions, for speed */
8 #include "mri_genalign_util.c"
9 #include "mri_genalign.c"
10 #include "mri_nwarp.c"
11
12 /*----------------------------------------------------------------------------*/
13
NWA_help(void)14 void NWA_help(void)
15 {
16 printf("\n"
17 "Usage: 3dNwarpApply [options]\n"
18 "\n"
19 "Program to apply a nonlinear 3D warp saved from 3dQwarp (or 3dNwarpCat, etc.)\n"
20 "to a 3D dataset, to produce a warped version of the source dataset.\n"
21 "\n"
22 "The '-nwarp' and '-source' options are MANDATORY. For both of these options,\n"
23 "as well as '-prefix', the input arguments after the option name are applied up\n"
24 "until an argument starts with the '-' character, or until the arguments run out.\n"
25 "\n"
26 "This program has been heavily modified [01 Dec 2014], including the following\n"
27 "major improvements:\n"
28 "(1) Allow catenation of warps with different grid spacings -- the functions\n"
29 " that deal with the '-nwarp' option will automatically deal with the grids.\n"
30 "(2) Allow input of affine warps with multiple time points, so that 3D+time\n"
31 " datasets can be warped with a time dependent '-nwarp' list.\n"
32 "(3) Allow input of multiple source datasets, so that several datasets can be\n"
33 " warped the same way at once. This operation is more efficient than running\n"
34 " 3dNwarpApply several times, since the auto-regridding and auto-catenation\n"
35 " in '-nwarp' will only have to be done once.\n"
36 " * Specification of the output dataset names can be done via multiple\n"
37 " arguments to the '-prefix' option, or via the new '-suffix' option.\n"
38 "\n"
39 "New Feature [28 Mar 2018]:\n"
40 "(4) If a source dataset contains complex numbers, then 3dNwarpApply will warp\n"
41 " the real and imaginary parts separately, combine them, and produce a\n"
42 " complex-valued dataset as output.\n"
43 " * Previously, the program would have warped the magnitude of the input\n"
44 " dataset and written out a float-valued dataset.\n"
45 " * No special option is needed to warp complex-valued datasets.\n"
46 " * If you WANT to warp the magnitude of a complex-valued dataset, you will\n"
47 " have to convert the dataset to a float dataset via 3dcalc, then use\n"
48 " 3dNwarpApply on THAT dataset instead.\n"
49 " * You cannot use option '-short' with complex-valued source datasets!\n"
50 " More precisely, you can try to use this option, but it will be ignored.\n"
51 " * This ability is added for those of you who deal with complex-valued\n"
52 " EPI datasets (I'm looking at YOU, O International Man of Mystery).\n"
53 "\n"
54 "OPTIONS:\n"
55 "--------\n"
56 " -nwarp www = 'www' is the name of the 3D warp dataset\n"
57 " (this is a mandatory option!)\n"
58 " ++ Multiple warps can be catenated here.\n"
59 " -->> Please see the lengthier discussion below on this feature!\n"
60 " -->> Also see the help for 3dNwarpCat for some more information\n"
61 " on the formats allowed for inputting warp fields; for\n"
62 " example, warping in one direction only (e.g., 'AP') is\n"
63 " possible.\n"
64 " ++ NOTE WELL: The interpretation of this option has changed somewhat,\n"
65 " as of 01 Dec 2014. In particular, this option is\n"
66 " generalized from the version in other programs, including\n"
67 " 3dNwarpCat, 3dNwarpFuncs, and 3dNwarpXYZ. The major\n"
68 " change is that multi-line matrix files are allowed to\n"
69 " be included in the 'www' mixture, so that the nonlinear\n"
70 " warp being calculated can be time-dependent.\n"
71 " In addition, the warps supplied need not all be on the\n"
72 " same 3D grid -- this ability lets you catenate a warp\n"
73 " defined on the EPI data grid with a warp defined on the\n"
74 " structural data grid (e.g.).\n"
75 "\n"
76 " -iwarp = After the warp specified in '-nwarp' is computed,\n"
77 " invert it. If the input warp would take a dataset\n"
78 " from space A to B, then the inverted warp will do\n"
79 " the reverse.\n"
80 " ++ The combination \"-iwarp -nwarp 'A B C'\" is equivalent\n"
81 " to \"-nwarp 'INV(C) INV(B) INV(A)'\" -- that is, inverting\n"
82 " each warp/matrix in the list *and* reversing their order.\n"
83 " ++ The '-iwarp' option is provided for convenience, and\n"
84 " may prove to be very slow for time-dependent '-nwarp' inputs.\n"
85 "\n"
86 " -affter aaa = *** THIS OPTION IS NO LONGER AVAILABLE ***\n"
87 " See the discussion of the new '-nwarp' option above to see\n"
88 " how to do include time-dependent matrix transformations\n"
89 " in this program.\n"
90 #if 0
91 "\n"
92 " -wfac fff = Scale the warp by factor 'fff' [default=1.0]\n"
93 " ++ This option doesn't really have much use, except\n"
94 " in making fun movies of how a brain volume deforms.\n"
95 " ++ Also note that '-wfac' applies to the warp AFTER\n"
96 " it is inverted, if '-iwarp' is used.\n"
97 #endif
98 "\n"
99 " -source sss = 'sss' is the name of the source dataset.\n"
100 " ++ That is, the dataset to be warped.\n"
101 " ++ Multiple datasets can be supplied here; they MUST\n"
102 " all be defined over the same 3D grid.\n"
103 " -->>** You can no longer simply supply the source\n"
104 " dataset as the last argument on the command line.\n"
105 "\n"
106 " -master mmm = 'mmm is the name of the master dataset.\n"
107 " ++ Which defines the output grid.\n"
108 " ++ If '-master' is not used, then output\n"
109 " grid is the same as the source dataset grid.\n"
110 " ++ It is often the case that it makes more sense to\n"
111 " use the '-nwarp' dataset as the master, since\n"
112 " that is the grid on which the transformation is\n"
113 " defined, and is (usually) the grid to which the\n"
114 " transformation 'pulls' the source data.\n"
115 " ++ You can use '-master WARP' or '-master NWARP'\n"
116 " for this purpose -- but ONLY if all the warps listed\n"
117 " in the '-nwarp' option have the same 3D grid structure.\n"
118 " ++ In particular, if the transformation includes a\n"
119 " long-distance translation, then the source dataset\n"
120 " grid may not have a lot of overlap with the source\n"
121 " dataset after it is transformed -- in this case, you\n"
122 " really want to use this '-master' option -- or you\n"
123 " will end up cutting of a lot of the output dataset\n"
124 " since it will not overlap with the source dataset.\n"
125 "\n"
126 " -newgrid dd = 'dd' is the new grid spacing (cubical voxels, in mm)\n"
127 " *OR = ++ This lets you resize the master dataset grid spacing.\n"
128 " -dxyz dd = for example, to bring EPI data to a 1 mm template, but at\n"
129 " a coarser resolution, use '-dxyz 2'.\n"
130 " ++ The same grid orientation as the source is used if\n"
131 " the '-master' option is not given.\n"
132 "\n"
133 " -interp iii = 'iii' is the interpolation mode\n"
134 " ++ Default interpolation mode is 'wsinc5' (slowest, bestest)\n"
135 " ++ Available modes are the same as in 3dAllineate:\n"
136 " NN linear cubic quintic wsinc5\n"
137 " ++ The same interpolation mode is used for the warp\n"
138 " itself (if needed) and then for the data being warped.\n"
139 " ++ The warp will be interpolated if the output dataset is\n"
140 " not on the same 3D grid as the warp itself, or if a warp\n"
141 " expression is used in the '-nwarp' option. Otherwise,\n"
142 " it won't need to be interpolated.\n"
143 "\n"
144 " -ainterp jjj = This option lets you specify a different interpolation mode\n"
145 " for the data than might be used for the warp.\n"
146 " ++ In particular, '-ainterp NN' would be most logical for\n"
147 " atlas datasets, where the data values being mapped are\n"
148 " integer labels.\n"
149 #if 0
150 "\n"
151 " -expad EE = Add 'EE' voxels to the warp on input, to help avoid any problems\n"
152 " that might arise if you are catenating multiple warps / matrices\n"
153 " where the nonlinear transform might be shifted far off its\n"
154 " original grid definition. Normally not needed, since the program\n"
155 " internally estimates how much padding is needed.\n"
156 " ++ This option does not affect the use of '-master WARP', which\n"
157 " still refers to the original grid on which the nonlinear\n"
158 " warp is defined.\n"
159 #endif
160 "\n"
161 " -prefix ppp = 'ppp' is the name of the new output dataset\n"
162 " ++ If more than 1 source dataset is supplied, then you\n"
163 " should supply more than one prefix. Otherwise, the\n"
164 " program will invent prefixes for each output, by\n"
165 " attaching the suffix '_Nwarp' to each source\n"
166 " dataset's prefix.\n"
167 "\n"
168 " -suffix sss = If the program generates prefixes, you can change the\n"
169 " default '_Nwarp' suffix to whatever you want (within\n"
170 " reason) by this option.\n"
171 " ++ His Holiness Emperor Zhark defines 'within reason', of course.\n"
172 " ++ By using '-suffix' and NOT using '-prefix', the program\n"
173 " will generate prefix names for all output datasets in\n"
174 " a systematic way -- this might be useful for some people.\n"
175 " ++ Note that only ONE suffix can be supplied even if many source\n"
176 " datasets are input -- unlike the case with '-prefix'.\n"
177 "\n"
178 " -short = Write output dataset using 16-bit short integers, rather than\n"
179 " the usual 32-bit floats.\n"
180 " ++ Intermediate values are rounded to the nearest integer.\n"
181 " No scaling is performed.\n"
182 " ++ This option is intended for use with '-ainterp' and for\n"
183 " source datasets that contain integral values.\n"
184 " ++ If the source dataset is complex-valued, this option will\n"
185 " be ignored.\n"
186 "\n"
187 " -wprefix wp = If this option is used, then every warp generated in the process\n"
188 " of application will be saved to a 3D dataset with prefix 'wp_XXXX',\n"
189 " where XXXX is the index of the sub-brick being created.\n"
190 " For example, '-wprefix Zork.nii' will create datasets with names\n"
191 " 'Zork_0000.nii', et cetera.\n"
192 "\n"
193 " -quiet = Don't be verbose :-(\n"
194 " -verb = Be extra verbose :-)\n"
195 "\n"
196 "SPECIFYING THE NONLINEAR WARP IN '-nwarp'\n"
197 "[If you are catenating warps, read this carefully!]\n"
198 "---------------------------------------------------\n"
199 "A single nonlinear warp (usually created by 3dQwarp) is an AFNI or NIfTI-1\n"
200 "dataset with 3 sub-bricks, holding the 3D displacements of each voxel.\n"
201 "(All coordinates and displacements are expressed in DICOM order.)\n"
202 "\n"
203 "The '-nwarp' option is used to specify the nonlinear transformation used\n"
204 "to create the output dataset from the source dataset. For many purposes,\n"
205 "the only input needed here is the name of a single dataset holding the\n"
206 "warp to be used.\n"
207 "\n"
208 "However, the '-nwarp' option also allows the catenation of a sequence of\n"
209 "spatial transformations (in short, 'warps') that will be combined before\n"
210 "being applied to the source dataset. Each warp is either a nonlinear\n"
211 "warp dataset or a matrix warp (a linear transformation of space).\n"
212 "\n"
213 "A single affine (or linear) warp is a set of 12 numbers, defining a 3x4 matrix\n"
214 " a11 a12 a13 a14\n"
215 " a21 a22 a23 a24\n"
216 " a31 a32 a33 a34\n"
217 "A matrix is stored on a single line, in a file with the extension\n"
218 "'.1D' or '.txt', in this order\n"
219 " a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34\n"
220 "For example, the identity matrix is given by\n"
221 " 1 0 0 0 0 1 0 0 0 0 1 0\n"
222 "This format is output by the '-1Dmatrix_save' options in 3dvolreg and\n"
223 "3dAllineate, for example.\n"
224 "\n"
225 "If the argument 'www' following '-nwarp' is made up of more than one warp\n"
226 "filename, separated by blanks, then the nonlinear warp to be used is\n"
227 "composed on the fly as needed to transform the source dataset. For\n"
228 "example,\n"
229 " -nwarp 'AA_WARP.nii BB.aff12.1D CC_WARP.nii'\n"
230 "specifies 3 spatial transformations, call them A(x), B(x), and C(x) --\n"
231 "where B(x) is just the 3-vector x multipled into the matrix in the\n"
232 "BB.aff12.1D file. The resulting nonlinear warp function N(x) is\n"
233 "obtained by applying these transformations in the order given, A(x) first:\n"
234 " N(x) = C( B( A(x) ) )\n"
235 "That is, the first warp A is applied to the output grid coordinate x,\n"
236 "then the second warp B to that results, then the third warp C. The output\n"
237 "coordinate y = C(B(A(x))) is the coordinate in the source dataset at which\n"
238 "the output value will be interpolated (for the voxel at coordinate x).\n"
239 "\n"
240 "The Proper Order of Catenated Warps:\n"
241 "....................................\n"
242 "To determine the correct order in which to input the warps, it is necessary\n"
243 "to understand what a warp of the source dataset actually computes. Call the\n"
244 "source image S(x) = (scalar) value of source image at voxel location x.\n"
245 "For each x in the output grid, the warped result is S(N(x)) -- that is,\n"
246 "N(x) tells where each output location x must be warped to in order to\n"
247 "find the corresponding value of the source S.\n"
248 "\n"
249 "N(x) does *NOT* tell to where an x in the source image must be moved to in\n"
250 "the output space -- which is what you might think if you mentally prioritize\n"
251 "the idea of 'warping the source image' or 'pushing the source image' -- DO NOT\n"
252 "THINK THIS WAY! It is better to think of N(x) as reaching out from x in the\n"
253 "output space to a location in the source space and then the program will\n"
254 "interpolate from the discrete source space grid at that location -- which\n"
255 "is unlikely to be exactly on a grid node. Another way to think of this is\n"
256 "that the warp 'pulls' the source image back to the coordinate system on which\n"
257 "the warp is defined.\n"
258 "\n"
259 "Now suppose the sequence of operations on an EPI dataset is\n"
260 " (1) Nonlinearly unwarp the dataset via warp AA_WARP.nii (perhaps\n"
261 " from 3dQwarp -plusminus).\n"
262 " (2) Perform linear volume registration on the result from (1) (with\n"
263 " program 3dvolreg) to get affine matrix file BB.aff12.1D -- which\n"
264 " will have 1 line per time point in the EPI dataset.\n"
265 " (3) Linearly register the structural volume to the EPI dataset\n"
266 " (via script align_epi_anat.py). Note that this step transforms\n"
267 " the structural volume to match the EPI, not the EPI to match the\n"
268 " structural volume, so this step does not affect the chain of\n"
269 " transformations being applied to the EPI dataset.\n"
270 " (4) Nonlinearly warp the structural image from (3) to MNI space via\n"
271 " warp CC_WARP.nii (generated by 3dQwarp).\n"
272 "Finally, the goal is to take the original EPI time series dataset, and\n"
273 "warp it directly to MNI space, including the time series registration for\n"
274 "each sub-brick in the dataset, with only one interplation being used --\n"
275 "rather than the 3 interpolations that would come by serially implementing\n"
276 "steps (1), (2), and (4). This one-big-step transformation can be done\n"
277 "with 3dNwarpApply using the '-nwarp' option:\n"
278 " -nwarp 'CC_WARP.nii BB.aff12.1D AA_WARP.nii'\n"
279 "that is, N(x) = A( B( C(x) ) ) -- the opposite order to the sample above,\n"
280 "and with the transformations occuring in the opposite order to the sequence\n"
281 "in which they were calculated. The reason for this apparent backwardness\n"
282 "is that the 'x' being transformed is on the output grid -- in this case, in\n"
283 "MNI-template space. So the warp C(x) transforms such an output grid 'x' to\n"
284 "the EPI-aligned structural space. The warp B(x) then transforms THAT\n"
285 "coordinate from aligned spaced back to the rotated head position of the subject.\n"
286 "And the warp A(x) transforms THAT coordinate back to the original grid that had\n"
287 "to be unwarped (e.g., from susceptibility and/or eddy current artifacts).\n"
288 "\n"
289 "Also note that in step (2), the matrix file BB.aff12.1D has one line for\n"
290 "each time point. When transforming a source dataset, the i-th time point\n"
291 "will be transformed by the warp computed using the i-th line from any\n"
292 "multi-line matrix file in the '-nwarp' specification. (If there are more\n"
293 "dataset time points than matrix lines, then the last line will be re-used.)\n"
294 "\n"
295 "In this way, 3dNwarpApply can be used to carry out time-dependent warping\n"
296 "of time-dependent datasets, provided that the time-dependence in the warp\n"
297 "only occurs in the affine (matrix) parts of the transformation.\n"
298 "\n"
299 "Note that the now-obsolete option '-affter' is subsumed into the new way\n"
300 "that '-nwarp' works. Formerly, the only time-dependent matrix had to\n"
301 "be specified as being at the end of the warp chain, and was given via\n"
302 "the '-affter' option. Now, a time-dependent matrix (or more than one)\n"
303 "can appear anywhere in the warp chain, so there is no need for a special\n"
304 "option. If you DID use '-affter', you will have to alter your script\n"
305 "simply by putting the final matrix filename at the end of the '-nwarp'\n"
306 "chain. (If this seems too hard, please consider another line of work.)\n"
307 "\n"
308 "The other 3dNwarp* programs that take the '-nwarp' option operate similarly,\n"
309 "but do NOT allow time-dependent matrix files. Those programs are built to\n"
310 "operate with one nonlinear warp, so allowing a time-dependent warp doesn't\n"
311 "make sense for them.\n"
312 "\n"
313 "NOTE: If a matrix is NOT time-dependent (just a single set of 12 numbers),\n"
314 " it can be input in the .Xat.1D format of 3 rows, each with 4 values:\n"
315 " a11 a12 a13 a14 } 1 0 0 0\n"
316 " a21 a22 a23 a24 } e.g, identity matrix = 0 1 0 0\n"
317 " a31 a32 a33 a34 } 0 0 1 0\n"
318 " This option is just for convenience. Remember that the coordinates\n"
319 " are DICOM order, and if your matrix comes from Some other PrograM\n"
320 " or from a Fine Software Library, you probably have to change some\n"
321 " signs in the matrix to get things to work correctly.\n"
322 "\n"
323 "RANDOM NOTES:\n"
324 "-------------\n"
325 "* At present, this program doesn't work with 2D warps, only with 3D.\n"
326 " (That is, each warp dataset must have 3 sub-bricks.)\n"
327 "\n"
328 "* At present, the output dataset is stored in float format, no matter what\n"
329 " absurd data format the input dataset uses (but cf. the '-short' option).\n"
330 "\n"
331 "* As described above, 3dNwarpApply allows you to catenate warps directly on\n"
332 " the command line, as if you used 3dNwarpCat before running 3dNwarpApply.\n"
333 " For example:\n"
334 "\n"
335 " ++ You have aligned dataset Fred+orig to MNI-affine space using @auto_tlrc,\n"
336 " giving matrix file Fred.Xaff12.1D\n"
337 "\n"
338 " ++ Then you further aligned from MNI-affine to MNI-qwarp via 3dQwarp,\n"
339 " giving warp dataset Fred_WARP+tlrc\n"
340 "\n"
341 " ++ You can combine the transformations and interpolate Fred+orig directly\n"
342 " to MNI-qwarp space using a command like\n"
343 " 3dNwarpApply -prefix Fred_final \\\n"
344 " -source Fred+orig \\\n"
345 " -master NWARP \\\n"
346 " -nwarp 'Fred_WARP+tlrc Fred.Xaff12.1D'\n"
347 " Note the warps to be catenated are enclosed in quotes to make a single\n"
348 " input argument passed to the program. The processing used for this\n"
349 " purpose is the same as in 3dNwarpCat -- see the help output for that\n"
350 " program for a little more information.\n"
351 "\n"
352 " ++ When you specify a nonlinear warp dataset, you can use the 'SQRT()' and\n"
353 " 'INV()' and 'INVSQRT()' operators, as well as the various 1D-to-3D\n"
354 " displacement prefixes ('AP:' 'RL:' 'IS:' 'VEC:', as well as 'FAC:') -- \n"
355 " for example, the following is a legal (and even useful) definition of a\n"
356 " warp herein:\n"
357 " 'SQRT(AP:epi_BU_yWARP+orig)'\n"
358 " where the 'AP:' transforms the y-displacements in epi_BU_ywarp+orig to a\n"
359 " full 3D warp (with x- and z-displacments set to zero), then calculates the\n"
360 " square root of that warp, then applies the result to some input dataset.\n"
361 " + This is a real example, where the y-displacement-only warp is computed between\n"
362 " blip-up and blip-down EPI datasets, and then the SQRT warp is applied to\n"
363 " warp them into the 'intermediate location' which should be better aligned\n"
364 " with the subject's anatomical datasets.\n"
365 " -->+ However: see also the '-plusminus' option for 3dQwarp for another way to\n"
366 " reach the same goal, as well as the unWarpEPI.py script.\n"
367 " + See the output of 3dNwarpCat -help for a little more information on the\n"
368 " 1D-to-3D warp prefixes ('AP:' 'RL:' 'IS:' 'VEC:').\n"
369 "\n"
370 " ++ You can scale the displacements in a 3D warp file via the 'FAC:' prefix, as in\n"
371 " FAC:0.6,0.4,-0.2:fred_WARP.nii\n"
372 " which will scale the x-displacements by 0.6, the y-displacements by 0.4, and\n"
373 " the z-displacments by -0.2.\n"
374 " + So if you need to reverse the sign of x- and y-displacments, since in AFNI\n"
375 " +x=Left and +y=Posterior while another package uses +x=Right and +y=Anterior,\n"
376 " you could use 'FAC:-1,-1,1:Warpdatasetname' to apply a warp from that\n"
377 " other software package.\n"
378 "\n"
379 " ++ You can also use 'IDENT(dataset)' to define a \"nonlinear\" 3D warp whose\n"
380 " grid is defined by the dataset header -- nothing else from the dataset will\n"
381 " be used. This warp will be filled with all zero displacements, which represents\n"
382 " the identity warp. The purpose of such an object is to let you apply a pure\n"
383 " affine warp -- since this program requires a '-nwarp' option, you can use\n"
384 " -nwarp 'IDENT(dataset)' to define the 3D grid for the 'nonlinear' 3D warp and\n"
385 " then catenate the affine warp.\n"
386 "\n"
387 "* PLEASE note that if you use the '-allineate' option in 3dQwarp, then the affine\n"
388 " warp is already included in the output nonlinear warp from 3dQwarp, and so it\n"
389 " does NOT need to be applied again in 3dNwarpApply! This mistake has been made\n"
390 " in the past, and the results were not good.\n"
391 "\n"
392 "* When using '-allineate' in 3dQwarp, and when there is a large coordinate shift\n"
393 " between the base and source datasets, then the _WARP dataset output by 3dQwarp\n"
394 " will cover a huge grid to encompass both the base and source. In turn, this\n"
395 " can cause 3dNwarpApply to need a lot of memory when it applies that warp.\n"
396 " ++ Some changes were made [Jan 2019] to reduce the size of this problem,\n"
397 " but it still exists.\n"
398 " ++ We have seen this most often in source datasets which have the (0,0,0)\n"
399 " point not in the middle of the volume, but at a corner of the volume.\n"
400 " Since template datasets (such as MNI152_2009_template_SSW.nii.gz) have\n"
401 " (0,0,0) inside the brain, a dataset with (0,0,0) at a corner of the 3D\n"
402 " volume will need a giant coordinate shift to match the template dataset.\n"
403 " And in turn, the encompassing grid that overlaps the source and template\n"
404 " (base) datasets will be huge.\n"
405 " ++ The simplest way to fix this problem is to do something like\n"
406 " @Align_Centers -base MNI152_2009_template_SSW.nii.gz -dset Fred.nii\n"
407 " which will produce dataset Fred_shft.nii, that will have its grid\n"
408 " center approximately lined up with the template (base) dataset.\n"
409 " And from then on, use Fred_shft.nii as your input dataset.\n"
410 ) ;
411
412 PRINT_AFNI_OMP_USAGE("3dNwarpApply",NULL) ; PRINT_COMPILE_DATE ;
413 exit(0) ;
414 }
415
416 /*----------------------------------------------------------------------------*/
417 /* This program is basically a wrapper for function THD_nwarp_dataset() */
418 /*----------------------------------------------------------------------------*/
419
main(int argc,char * argv[])420 int main( int argc , char *argv[] )
421 {
422 THD_3dim_dataset_array *dset_srcar=NULL , *dset_outar=NULL ;
423 THD_3dim_dataset *dset_mast=NULL , *dset_out=NULL , *dset_sss,*dset_ooo ;
424 Nwarp_catlist *nwc=NULL ; char *nwc_string=NULL ;
425 MRI_IMARR *imar_nwarp=NULL , *im_src ;
426 char **prefix = NULL ; int nprefix=0 ; char *suffix = "_Nwarp" ;
427 float dxyz_mast = 0.0f ;
428 int interp_code = MRI_WSINC5 ;
429 int ainter_code = -666 ;
430 int iarg , kk , verb=1 , iv , do_wmast=0 , do_iwarp=0 ;
431 mat44 src_cmat, nwarp_cmat, mast_cmat ;
432 MRI_IMAGE *fim , *wim ; float *ip,*jp,*kp ;
433 int nx,ny,nz,nxyz , toshort=0 ;
434 float wfac=1.0f ;
435 #if 0
436 MRI_IMAGE *awim=NULL ; THD_3dim_dataset *dset_nwarp=NULL ;
437 #endif
438 int expad=0 ; /* 26 Aug 2014 */
439
440 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
441
442 /**----------------------------------------------------------------------*/
443 /**----------------- Help the pitifully ignorant user? -----------------**/
444
445 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ) NWA_help() ;
446
447 /**-------- bookkeeping and marketing --------**/
448
449 mainENTRY("3dNwarpApply"); machdep();
450 AFNI_logger("3dNwarpApply",argc,argv);
451 PRINT_VERSION("3dNwarpApply"); AUTHOR("Zhark the Warped");
452 (void)COX_clock_time() ;
453 putenv("AFNI_WSINC5_SILENT=YES") ;
454
455 /**--- process command line options ---**/
456
457 iarg = 1 ;
458 while( iarg < argc && argv[iarg][0] == '-' ){
459
460 /*---------------*/
461
462 if( strcasecmp(argv[iarg],"-quiet") == 0 ){
463 verb = 0 ; iarg++ ; continue ;
464 }
465 if( strcasecmp(argv[iarg],"-verb") == 0 ){
466 verb++ ; iarg++ ; continue ;
467 }
468 if( strcasecmp(argv[iarg],"-short") == 0 ){
469 toshort = 1 ; iarg++ ; continue ;
470 }
471 if( strcmp(argv[iarg],"-") == 0 ){ /* just skip */
472 iarg++ ; continue ;
473 }
474
475 /*---------------*/
476
477 if( strcasecmp(argv[iarg],"-expad") == 0 ){ /* 26 Aug 2014 */
478 if( ++iarg >= argc ) ERROR_exit("need arg after %s",argv[iarg-1]) ;
479 expad = (int)strtod(argv[iarg],NULL) ;
480 if( expad < 0 ) expad = 0 ;
481 iarg++ ; continue ;
482 }
483
484 /*---------------*/
485
486 #if 0
487 if( strcasecmp(argv[iarg],"-affter") == 0 || strcasecmp(argv[iarg],"-after") == 0 ){
488 MRI_IMAGE *qim ;
489 if( awim != NULL ) ERROR_exit("Can't have multiple %s options :-(" ,argv[iarg] );
490 if( ++iarg >= argc ) ERROR_exit("No argument after '%s' :-(" ,argv[iarg-1]);
491 awim = mri_read_1D(argv[iarg]) ;
492 if( awim == NULL ) ERROR_exit("Can't read affine warp data from file '%s'",argv[iarg] );
493 if( awim->nx == 3 && awim->ny == 4 ){ /* in 3x4 'Xat.1D' format? */
494 MRI_IMAGE *qim = mri_rowmajorize_1D(awim); mri_free(awim); awim = qim; /* make it 1x12 */
495 awim->nx = 12 ; awim->ny = 1 ;
496 } else { /* should be in N x 12 format */
497 if( awim->ny < 12 )
498 ERROR_exit("Affine warp file '%s': fewer than 12 values per row" ,argv[iarg] );
499 else if( awim->ny > 12 )
500 WARNING_message("Affine warp file '%s': more than 12 values per row" ,argv[iarg] );
501 qim = mri_transpose(awim) ; mri_free(awim) ; awim = qim ; /* flip to column major order */
502 }
503 iarg++ ; continue ;
504 }
505 #endif
506
507 /*---------------*/
508
509 if( strcasecmp(argv[iarg],"-nwarp") == 0 ||
510 strcasecmp(argv[iarg],"-warp" ) == 0 ||
511 strcasecmp(argv[iarg],"-qwarp") == 0 ){
512
513 int nwclen=0 ; /* length of string to create */
514
515 if( nwc_string != NULL )
516 ERROR_exit("Can't have multiple %s options :-(",argv[iarg]) ;
517 if( ++iarg >= argc )
518 ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
519
520 for( kk=iarg ; kk < argc && argv[kk][0] != '-' ; kk++ ){ /* to '-' or end */
521 nwclen += strlen(argv[kk])+4 ; /* iarg->kk 27 Jul, 2017 [rickr] */
522 }
523 if( nwclen <= 0 )
524 ERROR_exit("Invalid argument after '%s'",argv[iarg-1]) ;
525 nwc_string = (char *)malloc(sizeof(char)*nwclen) ; nwc_string[0] = '\0' ;
526 for( kk=iarg ; kk < argc && argv[kk][0] != '-' ; kk++,iarg++ ){
527 strcat(nwc_string,argv[kk]) ; strcat(nwc_string," ") ;
528 }
529 continue ; /* don't do iarg++ since iarg already at next option (or end) */
530 }
531
532 /*---------------*/
533
534 if( strcasecmp(argv[iarg],"-source") == 0 ){
535 int nbad=0 ; char *gs=NULL , *hs=NULL ;
536 if( dset_srcar != NULL )
537 ERROR_exit("Can't have multiple %s options :-(",argv[iarg]) ;
538 if( ++iarg >= argc )
539 ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
540 INIT_3DARR(dset_srcar) ;
541 for( ; iarg < argc && argv[iarg][0] != '-' ; iarg++ ){ /* to '-' or end */
542 dset_sss = THD_open_dataset( argv[iarg] ) ;
543 if( dset_sss == NULL ){
544 ERROR_message("Can't open -source dataset '%s' :-(",argv[iarg]) ;
545 nbad++ ; continue ;
546 }
547 if( verb ) INFO_message("opened source dataset '%s'",argv[iarg]) ;
548 DSET_COPYOVER_REAL(dset_sss) ;
549 ADDTO_3DARR(dset_srcar,dset_sss) ;
550 if( dset_srcar->num == 1 ){ /* check geometries of datasets match 1st */
551 gs = EDIT_get_geometry_string(dset_sss) ;
552 } else {
553 hs = EDIT_get_geometry_string(dset_sss) ;
554 if( EDIT_geometry_string_diff(gs,hs) > 0.01f ){
555 ERROR_message(
556 "Geometry mismatch for -source dataset '%s' (compared to first)",
557 argv[iarg]);
558 nbad++ ;
559 }
560 free(hs) ; hs = NULL ;
561 }
562 }
563 if( dset_srcar->num <= 0 )
564 ERROR_exit("No good source datasets found!") ;
565 else if( nbad > 0 )
566 ERROR_exit("Can't continue after above error%s",(nbad > 0)?"s":"\0") ;
567 if( gs != NULL ) free(gs) ;
568 continue ; /* don't do iarg++ ; iarg already at next option (or end) */
569 }
570
571 /*---------------*/
572
573 if( strncasecmp(argv[iarg],"-prefix",5) == 0 ){
574 int nbad=0 , pp,qq ;
575 if( prefix != NULL )
576 ERROR_exit("Can't have multiple %s options :-(",argv[iarg]) ;
577 if( ++iarg >= argc )
578 ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
579
580 for( kk=iarg ; kk < argc && argv[kk][0] != '-' ; kk++ ){ /* to '-' or end */
581 if( !THD_filename_ok(argv[kk]) || strcasecmp(argv[kk],"NULL")==0 ){
582 ERROR_message("Badly formed prefix '%s' :-(",argv[kk]);
583 nbad++ ;
584 }
585 nprefix++ ;
586 }
587 if( nprefix <= 0 )
588 ERROR_exit("No prefix name given after '%s' ???",argv[iarg-1]) ;
589 else if( nbad > 0 )
590 ERROR_exit("Can't continue after above error%s",(nbad > 0)?"s":"\0") ;
591
592 prefix = (char **)malloc(sizeof(char *)*nprefix) ;
593 for( pp=0,kk=iarg ; kk < argc && argv[kk][0] != '-' ; kk++,pp++ ){
594 prefix[pp] = strdup(argv[kk]) ;
595 for( qq=0 ; qq < pp ; qq++ ){
596 if( strcmp(prefix[pp],prefix[qq]) == 0 ){
597 ERROR_message("Prefix name '%s' is duplicated -- this is FORBIDDEN!",
598 prefix[pp] ) ;
599 nbad++ ;
600 }
601 }
602 }
603 if( nbad > 0 )
604 ERROR_exit("Can't continue after above error%s",(nbad > 0)?"s":"\0") ;
605 iarg = kk ; continue ;
606 }
607
608 /*---------------*/
609
610 if( strncasecmp(argv[iarg],"-wprefix",5) == 0 ){ /* 15 Mar 2021 */
611 if( ++iarg >= argc )
612 ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
613 if( strcmp(argv[iarg],"NULL") != 0 )
614 THD_set_nwarp_apply_prefix( argv[iarg] ) ;
615 else
616 THD_set_nwarp_apply_prefix( NULL ) ;
617 iarg++ ; continue ;
618 }
619
620 /*---------------*/
621
622 if( strncasecmp(argv[iarg],"-suffix",5) == 0 ){
623 if( ++iarg >= argc )
624 ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
625 if( !THD_filename_pure(argv[iarg]) )
626 ERROR_exit("Illegal name '%s' after option '%s'",
627 argv[iarg],argv[iarg-1]) ;
628 if( strlen(argv[iarg]) > 128 )
629 ERROR_exit("Too long name '%s' after option '%s'",
630 argv[iarg],argv[iarg-1]) ;
631 suffix = strdup(argv[iarg]) ;
632 iarg++ ; continue ;
633 }
634
635 /*---------------*/
636
637 if( strcasecmp(argv[iarg],"-master") == 0 ){
638 if( dset_mast != NULL || do_wmast )
639 ERROR_exit("Can't have multiple %s options :-(",argv[iarg]) ;
640 if( ++iarg >= argc )
641 ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
642 if( strcmp(argv[iarg],"WARP") == 0 || strcmp(argv[iarg],"NWARP") == 0 ){
643 do_wmast = 1 ;
644 if( verb ) INFO_message("-master dataset is the input warp") ;
645 } else {
646 dset_mast = THD_open_dataset( argv[iarg] ) ;
647 if( dset_mast == NULL )
648 ERROR_exit("can't open -master dataset '%s' :-(",argv[iarg]);
649 DSET_COPYOVER_REAL(dset_mast) ;
650 if( verb ) INFO_message("-master dataset is '%s'",argv[iarg]) ;
651 }
652 iarg++ ; continue ;
653 }
654
655 /*---------------*/
656
657 if( strcasecmp(argv[iarg],"-mast_dxyz") == 0 ||
658 strcasecmp(argv[iarg],"-dxyz_mast") == 0 ||
659 strcasecmp(argv[iarg],"-dxyz") == 0 ||
660 strcasecmp(argv[iarg],"-newgrid" ) == 0 ){
661
662 if( ++iarg >= argc )
663 ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
664 dxyz_mast = (float)strtod(argv[iarg],NULL) ;
665 if( dxyz_mast <= 0.01f )
666 ERROR_exit("Illegal value '%s' after -mast_dxyz :-(",argv[iarg]) ;
667 if( dxyz_mast <= 0.5f )
668 WARNING_message("Small value %g after -mast_dxyz :-?",dxyz_mast) ;
669 INFO_message("output grid size = %.3g mm",dxyz_mast) ;
670 iarg++ ; continue ;
671 }
672
673 /*---------------*/
674
675 if( strcasecmp(argv[iarg],"-NN") == 0 ||
676 strncasecmp(argv[iarg],"-nearest",6) == 0 ){
677 CW_interp = interp_code = MRI_NN ; iarg++ ; continue ;
678 }
679 if( strncasecmp(argv[iarg],"-linear",4) ==0 ||
680 strncasecmp(argv[iarg],"-trilinear",6)==0 ){
681 CW_interp = interp_code = MRI_LINEAR ; iarg++ ; continue ;
682 }
683 if( strncasecmp(argv[iarg],"-cubic",4)==0 ||
684 strncasecmp(argv[iarg],"-tricubic",6)==0 ){
685 CW_interp = interp_code = MRI_CUBIC ; iarg++ ; continue ;
686 }
687 if( strncasecmp(argv[iarg],"-quintic",4)==0 ||
688 strncasecmp(argv[iarg],"-triquintic",6)==0 ){
689 CW_interp = interp_code = MRI_QUINTIC ; iarg++ ; continue ;
690 }
691 if( strncasecmp(argv[iarg],"-wsinc",5) == 0 ){
692 CW_interp = interp_code = MRI_WSINC5 ; iarg++ ; continue ;
693 }
694 if( strncasecmp(argv[iarg],"-interp",5)==0 ){
695 char *inam ;
696 if( ++iarg >= argc ) ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
697 inam = argv[iarg] ; if( *inam == '-' ) inam++ ;
698 if( strcasecmp(inam,"NN")==0 || strncasecmp(inam,"nearest",5)==0 )
699 CW_interp = interp_code = MRI_NN ;
700 else
701 if( strncasecmp(inam,"linear",3)==0 ||
702 strncasecmp(inam,"trilinear",5)==0 )
703 CW_interp = interp_code = MRI_LINEAR ;
704 else
705 if( strncasecmp(inam,"cubic",3)==0 ||
706 strncasecmp(inam,"tricubic",5)==0 )
707 CW_interp = interp_code = MRI_CUBIC ;
708 else
709 if( strncasecmp(inam,"quintic",3)==0 ||
710 strncasecmp(inam,"triquintic",5)==0 )
711 CW_interp = interp_code = MRI_QUINTIC ;
712 else
713 if( strncasecmp(inam,"WSINC",5)==0 )
714 CW_interp = interp_code = MRI_WSINC5 ;
715 else
716 ERROR_exit("Unknown code '%s' after '%s' :-(",argv[iarg],argv[iarg-1]) ;
717 iarg++ ; continue ;
718 }
719 if( strncasecmp(argv[iarg],"-ainterp",5)==0 ){
720 char *inam ;
721 if( ++iarg >= argc ) ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
722 inam = argv[iarg] ; if( *inam == '-' ) inam++ ;
723 if( strcasecmp(inam,"NN")==0 || strncasecmp(inam,"nearest",5)==0 )
724 ainter_code = MRI_NN ;
725 else
726 if( strncasecmp(inam,"linear",3)==0 ||
727 strncasecmp(inam,"trilinear",5)==0 )
728 ainter_code = MRI_LINEAR ;
729 else
730 if( strncasecmp(inam,"cubic",3)==0 ||
731 strncasecmp(inam,"tricubic",5)==0 )
732 ainter_code = MRI_CUBIC ;
733 else
734 if( strncasecmp(inam,"quintic",3)==0 ||
735 strncasecmp(inam,"triquintic",5)==0 )
736 ainter_code = MRI_QUINTIC ;
737 else
738 if( strncasecmp(inam,"WSINC",5)==0 )
739 ainter_code = MRI_WSINC5 ;
740 else
741 ERROR_exit("Unknown code '%s' after '%s' :-(",argv[iarg],argv[iarg-1]) ;
742 iarg++ ; continue ;
743 }
744
745 /*---------------*/
746
747 if( strcasecmp(argv[iarg],"-iwarp") == 0 ){
748 do_iwarp = 1 ; iarg++ ; continue ;
749 }
750
751 /*---------------*/
752
753 ERROR_message("Mysteriously Unknown option '%s' :-( :-( :-(",argv[iarg]) ;
754 suggest_best_prog_option(argv[0],argv[iarg]) ;
755 exit(1) ;
756
757 }
758
759 /*-------- check inputs to see if the user is completely demented ---------*/
760
761 if( nwc_string == NULL )
762 ERROR_exit("No -nwarp option? How do you want to warp? :-(") ;
763
764 if( dset_srcar == NULL )
765 ERROR_exit("No -source option? What do you want to do? :-(") ;
766
767 #ifdef USING_MCW_MALLOC
768 if( verb > 1 ){ enable_mcw_malloc() ; }
769 #endif
770
771 /*-- deal with extra or insufficient prefixes --*/
772
773 if( nprefix > dset_srcar->num ){ /* extra */
774
775 INFO_message(
776 "-prefix had %d names, but only %d source datasets ==> ignoring any extra prefixes",
777 nprefix , dset_srcar->num ) ;
778
779 } else if( nprefix < dset_srcar->num ){ /* insufficient */
780 char *cpt ;
781
782 INFO_message("Number of -prefix names (%d) is fewer than number of source datasets (%d)",
783 nprefix , dset_srcar->num ) ;
784 INFO_message("Making up the missing prefix names from Everest-thin air:") ;
785
786 prefix = (char **)realloc(prefix,sizeof(char *)*dset_srcar->num) ;
787 for( kk=nprefix ; kk < dset_srcar->num ; kk++ ){
788 prefix[kk] = (char *)malloc(sizeof(char)*THD_MAX_NAME) ;
789 dset_sss = DSET_IN_3DARR(dset_srcar,kk) ;
790 strcpy( prefix[kk] , DSET_PREFIX(dset_sss) ) ;
791 cpt = strcasestr(prefix[kk],".nii") ; if( cpt != NULL ) *cpt = '\0' ;
792 strcat(prefix[kk],suffix) ; if( cpt != NULL ) strcat(prefix[kk],".nii") ;
793 ININFO_message("Source %s ==> Prefix %s",DSET_PREFIX(dset_sss),prefix[kk]) ;
794 }
795
796 }
797
798 /*-- default interp codes --*/
799
800 if( ainter_code < 0 ) ainter_code = interp_code ;
801
802 /* verb_nww = global variable used in warping functions in mri_nwarp.c */
803
804 verb_nww = AFNI_yesenv("AFNI_DEBUG_WARP") ? 2 : verb ;
805
806 /*--- read and setup the list of nonlinear warps ---*/
807
808 if( verb || verb_nww ) fprintf(stderr,"++ Processing -nwarp ") ;
809
810 /* read the list */
811
812 nwc = IW3D_read_nwarp_catlist( nwc_string ) ;
813 if( verb && verb_nww < 2 ) fprintf(stderr,"\n") ;
814
815 if( nwc == NULL )
816 ERROR_exit("Cannot process warp string '%s'",nwc_string) ;
817
818 if( do_iwarp ) nwc->flags |= NWC_INVERT_MASK ; /* set inversion flag */
819
820 /* set up the geometry of the constituent warps */
821
822 CW_extra_pad = expad ; /* and ask for some addition padding */
823
824 /* create the master dataset */
825
826 if( do_wmast ){
827 if( nwc->actual_geomstring == NULL )
828 ERROR_exit("Cannot use '-master NWARP' when '-nwarp' contains datasets with mis-matched grids!") ;
829 dset_mast = EDIT_geometry_constructor( nwc->actual_geomstring , "Zhark_Likes_Elephants" ) ;
830 } else if( dset_mast == NULL ){
831 dset_mast = DSET_IN_3DARR(dset_srcar,0) ;
832 }
833
834 if( dxyz_mast > 0.0f ){
835 THD_3dim_dataset *qset ; double dxyz = (double)dxyz_mast ;
836 qset = r_new_resam_dset( dset_mast , NULL ,
837 dxyz,dxyz,dxyz ,
838 NULL , RESAM_NN_TYPE , NULL , 0 , 0) ;
839 if( qset != NULL ){
840 dset_mast = qset ;
841 THD_daxes_to_mat44(dset_mast->daxes) ;
842 }
843 }
844
845 /* here is where geometry of nonlinear warps is harmonized (if needed) */
846
847 if( nwc->actual_geomstring == NULL )
848 IW3D_set_geometry_nwarp_catlist( nwc, EDIT_get_geometry_string(dset_mast) ) ;
849
850 /* combine warps (matrices and datasets) to the extent possible */
851
852 IW3D_reduce_nwarp_catlist( nwc ) ; /* IW3D_read_nwarp_catlist may do this */
853 NI_sleep(1) ;
854
855 /*--------- the actual work of warping ---------*/
856
857 if( verb > 1 || verb_nww > 1 )
858 INFO_message(".......... Starting dataset warping ..........") ;
859
860 dset_outar = THD_nwarp_dataset_array( nwc , dset_srcar , dset_mast , prefix ,
861 interp_code,ainter_code , 0.0f , wfac , 0 ) ;
862
863 if( dset_outar == NULL ) ERROR_exit("Warping failed for some reason :-(((") ;
864
865 /*-- conversion to shorts? --*/
866
867 if( toshort ){
868 int iv,nxyz ; short *sar ; float *far ;
869 for( kk=0 ; kk < dset_outar->num ; kk++ ){
870 dset_ooo = DSET_IN_3DARR(dset_outar,kk) ;
871 nxyz = DSET_NVOX(dset_ooo) ;
872 for( iv=0 ; iv < DSET_NVALS(dset_ooo) ; iv++ ){
873 if( DSET_BRICK_TYPE(dset_ooo,iv) == MRI_float ){
874 far = (float *)DSET_ARRAY(dset_ooo,iv) ;
875 sar = (short *)malloc(sizeof(short)*nxyz) ;
876 EDIT_coerce_type( nxyz , MRI_float,far , MRI_short,sar ) ;
877 EDIT_substitute_brick( dset_ooo , iv , MRI_short,sar ) ;
878 if( iv == 0 )
879 INFO_message("Converting %s to shorts",DSET_HEADNAME(dset_ooo)) ;
880 } else {
881 if( iv == 0 )
882 INFO_message("Cannot convert non-float %s to shorts",DSET_HEADNAME(dset_ooo)) ;
883 }
884 }
885 }
886 }
887
888 /*-- atlas space? --*/
889
890 if( dset_mast != NULL && dset_mast->atlas_space != NULL ){
891 for( kk=0 ; kk < dset_outar->num ; kk++ ){
892 dset_ooo = DSET_IN_3DARR(dset_outar,kk) ;
893 MCW_strncpy( dset_ooo->atlas_space , dset_mast->atlas_space , THD_MAX_NAME ) ;
894 }
895 }
896
897 /*-- output! --*/
898
899 for( kk=0 ; kk < dset_outar->num ; kk++ ){
900 dset_sss = DSET_IN_3DARR(dset_srcar,kk) ;
901 dset_ooo = DSET_IN_3DARR(dset_outar,kk) ;
902 tross_Copy_History( dset_sss , dset_ooo ) ; /* hysterical records */
903 tross_Make_History( "3dNwarpApply" , argc,argv , dset_ooo ) ;
904 DSET_write(dset_ooo) ; if( verb ) WROTE_DSET(dset_ooo) ;
905 DSET_delete(dset_ooo) ; DSET_delete(dset_sss) ;
906 }
907
908 /*-- exeunt omnes, pursued by a bear --*/
909
910 if( verb ){
911 double cput=COX_cpu_time() , clkt=COX_clock_time() ;
912 if( cput >= 0.1 && clkt >= 0.1 )
913 INFO_message("total CPU time = %.1f sec Elapsed = %.1f\n",cput,clkt) ;
914 else if( clkt >= 0.1 )
915 INFO_message("total Elapsed time = %.1f sec\n",clkt) ;
916 }
917
918 exit(0) ;
919 }
920