1 #include "mrilib.h"
2 #include "r_new_resam_dset.h"
3
4 #ifdef USE_OMP
5 #include <omp.h>
6 #endif
7
8 #include "mri_genalign.c"
9 #include "mri_genalign_util.c"
10 #include "mri_nwarp.c"
11
12 /*----------------------------------------------------------------------------*/
13
NWC_help(void)14 void NWC_help(void)
15 {
16 printf(
17 "Usage: 3dNwarpCat [options] warp1 warp2 ...\n"
18 "------\n"
19 " * This program catenates (composes) 3D warps defined on a grid,\n"
20 " OR via a matrix.\n"
21 " ++ All transformations are from DICOM xyz (in mm) to DICOM xyz.\n"
22 "\n"
23 " * Matrix warps are in files that end in '.1D' or in '.txt'. A matrix\n"
24 " warp file should have 12 numbers in it, as output (for example), by\n"
25 " '3dAllineate -1Dmatrix_save'.\n"
26 " ++ The matrix (affine) warp can have either 12 numbers on one row,\n"
27 " or be in the 3x4 format.\n"
28 " ++ The 12-numbers-on-one-row format is preferred, and is the format\n"
29 " output by the '-1Dmatrix_save' option in 3dvolreg and 3dAllineate.\n"
30 " ++ The matrix warp is a transformation of coordinates, not voxels,\n"
31 " and its use presumes the correctness of the voxel-to-coordinate\n"
32 " transformation stored in the header of the datasets involved.\n"
33 "\n"
34 " * Nonlinear warps are in dataset files (AFNI .HEAD/.BRIK or NIfTI .nii)\n"
35 " with 3 sub-bricks giving the DICOM order xyz grid displacements in mm.\n"
36 " ++ Note that it is not required that the xyz order of voxel storage be in\n"
37 " DICOM order, just that the displacements be in DICOM order (and sign).\n"
38 " ++ However, it is important that the warp dataset coordinate order be\n"
39 " properly specified in the dataset header, since warps are applied\n"
40 " based on coordinates, not on voxels.\n"
41 " ++ Also note again that displacements are in mm, NOT in voxel.\n"
42 " ++ You can 'edit' the warp on the command line by using the 'FAC:'\n"
43 " scaling prefix, described later. This input editing could be used\n"
44 " to change the sign of the xyz displacments, if needed.\n"
45 "\n"
46 " * If all the input warps are matrices, then the output is a matrix\n"
47 " and will be written to the file 'prefix.aff12.1D'.\n"
48 " ++ Unless the prefix already contains the string '.1D', in which case\n"
49 " the filename is just the prefix.\n"
50 " ++ If 'prefix' is just 'stdout', then the output matrix is written\n"
51 " to standard output.\n"
52 " ++ In any of these cases, the output format is 12 numbers in one row.\n"
53 "\n"
54 " * If any of the input warps are datasets, they must all be defined on\n"
55 " the same 3D grid!\n"
56 " ++ And of course, then the output will be a dataset on the same grid.\n"
57 " ++ However, you can expand the grid using the '-expad' option.\n"
58 "\n"
59 " * The order of operations in the final (output) warp is, for the\n"
60 " case of 3 input warps:\n"
61 "\n"
62 " OUTPUT(x) = warp3( warp2( warp1(x) ) )\n"
63 "\n"
64 " That is, warp1 is applied first, then warp2, et cetera.\n"
65 " The 3D x coordinates are taken from each grid location in the\n"
66 " first dataset defined on a grid.\n"
67 "\n"
68 " * For example, if you aligned a dataset to a template with @auto_tlrc,\n"
69 " then further refined the alignment with 3dQwarp, you would do something\n"
70 " like this:\n"
71 " warp1 is the output of 3dQwarp\n"
72 " warp2 is the matrix from @auto_tlrc\n"
73 " This is the proper order, since the desired warp takes template xyz\n"
74 " to original dataset xyz, and we have\n"
75 " 3dQwarp warp: takes template xyz to affinely aligned xyz, and\n"
76 " @auto_tlrc matrix: takes affinely aligned xyz to original xyz\n"
77 "\n"
78 " 3dNwarpCat -prefix Fred_total_WARP -warp1 Fred_WARP+tlrc.HEAD -warp2 Fred.Xat.1D \n"
79 "\n"
80 " The dataset Fred_total_WARP+tlrc.HEAD could then be used to transform original\n"
81 " datasets directly to the final template space, as in\n"
82 "\n"
83 " 3dNwarpApply -prefix Wilma_warped \\\n"
84 " -nwarp Fred_total_WARP+tlrc \\\n"
85 " -source Wilma+orig \\\n"
86 " -master Fred_total_WARP+tlrc\n"
87 "\n"
88 " * If you wish to invert a warp before it is used here, supply its\n"
89 " input name in the form of\n"
90 " INV(warpfilename)\n"
91 " To produce the inverse of the warp in the example above:\n"
92 "\n"
93 " 3dNwarpCat -prefix Fred_total_WARPINV \\\n"
94 " -warp2 'INV(Fred_WARP+tlrc.HEAD)' \\\n"
95 " -warp1 'INV(Fred.Xat.1D)' \n"
96 "\n"
97 " Note the order of the warps is reversed, in addition to the use of 'INV()'.\n"
98 "\n"
99 " * The final warp may also be inverted simply by adding the '-iwarp' option, as in\n"
100 "\n"
101 " 3dNwarpCat -prefix Fred_total_WARPINV -iwarp -warp1 Fred_WARP+tlrc.HEAD -warp2 Fred.Xat.1D \n"
102 "\n"
103 " * Other functions you can apply to modify a 3D dataset warp are:\n"
104 " SQRT(datasetname) to get the square root of a warp\n"
105 " SQRTINV(datasetname) to get the inverse square root of a warp\n"
106 " However, you can't do more complex expressions, such as 'SQRT(SQRT(warp))'.\n"
107 " If you think you need something so rococo, use 3dNwarpCalc. Or think again.\n"
108 "\n"
109 " * You can also manufacture a 3D warp from a 1-brick dataset with displacments\n"
110 " in a single direction. For example:\n"
111 " AP:0.44:disp+tlrc.HEAD (note there are no blanks here!)\n"
112 " means to take the 1-brick dataset disp+tlrc.HEAD, scale the values inside\n"
113 " by 0.44, then load them into the y-direction displacements of a 3-brick 3D\n"
114 " warp, and fill the other 2 directions with zeros. The prefixes you can use\n"
115 " here for the 1-brick to 3-brick displacment trick are\n"
116 " RL: for x-displacements (Right-to-Left)\n"
117 " AP: for y-displacements (Anterior-to-Posterior)\n"
118 " IS: for z-displacements (Inferior-to-Superior)\n"
119 " VEC:a,b,c: for displacements in the vector direction (a,b,c),\n"
120 " which vector will be scaled to be unit length.\n"
121 " Following the prefix's colon, you can put in a scale factor followed\n"
122 " by another colon (as in '0.44:' in the example above). Then the name\n"
123 " of the dataset with the 1D displacments follows.\n"
124 " * You might reasonably ask of what possible value is this peculiar format?\n"
125 " This was implemented to use Bz fieldmaps for correction of EPI datasets,\n"
126 " which are distorted only along the phase-encoding direction. This format\n"
127 " for specifying the input dataset (the fieldmap) is built to make the\n"
128 " scripting a little easier. Its principal use is in the program 3dNwarpApply.\n"
129 "\n"
130 " * You can scale the displacements in a 3D warp file via the 'FAC:' prefix, as in\n"
131 " FAC:0.6,0.4,-0.2:fred_WARP.nii\n"
132 " which will scale the x-displacements by 0.6, the y-displacements by 0.4, and\n"
133 " the z-displacments by -0.2.\n"
134 "\n"
135 " * Finally, you can input a warp catenation string directly as in the '-nwarp'\n"
136 " option of 3dNwarpApply, as in\n"
137 "\n"
138 " 3dNwarpCat -prefix Fred_total_WARP 'Fred_WARP+tlrc.HEAD Fred.Xat.1D' \n"
139 "\n"
140 "\n"
141 "OPTIONS\n"
142 "-------\n"
143 " -interp iii == 'iii' is the interpolation mode:\n"
144 " ++ Modes allowed are a subset of those in 3dAllineate:\n"
145 " linear quintic wsinc5\n"
146 " ++ The default interpolation mode is 'wsinc5'.\n"
147 " ++ 'linear' is much faster but less accurate.\n"
148 " ++ 'quintic' is between 'linear' and 'wsinc5',\n"
149 " in both accuracy and speed.\n"
150 "\n"
151 " -verb == print (to stderr) various fun messages along the road.\n"
152 "\n"
153 " -prefix ppp == prefix name for the output dataset that holds the warp.\n"
154 " -space sss == attach string 'sss' to the output dataset as its atlas\n"
155 " space marker.\n"
156 "\n"
157 " -warp1 ww1 == alternative way to specify warp#1\n"
158 " -warp2 ww2 == alternative way to specify warp#2 (etc.)\n"
159 " ++ If you use any '-warpX' option for X=1..99, then\n"
160 " any addition warps specified after all command\n"
161 " line options appear AFTER these enumerated warps.\n"
162 " That is, '-warp1 A+tlrc -warp2 B+tlrc C+tlrc'\n"
163 " is like using '-warp3 C+tlrc'.\n"
164 " ++ At most 99 warps can be used. If you need more,\n"
165 " PLEASE back away from the computer slowly, and\n"
166 " get professional counseling.\n"
167 "\n"
168 " -iwarp == Invert the final warp before output.\n"
169 "\n"
170 " -expad PP == Pad the nonlinear warps by 'PP' voxels in all directions.\n"
171 " The warp displacements are extended by linear extrapolation\n"
172 " from the faces of the input grid.\n"
173 ) ;
174
175 printf(
176 "\n"
177 "AUTHOR -- RWCox -- March 2013\n"
178 ) ;
179
180 PRINT_AFNI_OMP_USAGE("3dNwarpCat",NULL) ; PRINT_COMPILE_DATE ;
181 exit(0) ;
182 }
183
184 /*----------------------------------------------------------------------------*/
185
186 #define NWMAX 99
187
188 static int nwtop=0 ;
189 static char *cwarp[NWMAX] ;
190 static char *sname=NULL ;
191
192 static int verb=0 , interp_code=MRI_WSINC5 ;
193
194 /*----------------------------------------------------------------------------*/
195 /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
196 /*----------------------------------------------------------------------------*/
197
main(int argc,char * argv[])198 int main( int argc , char *argv[] )
199 {
200 int iarg=1 , ii , do_iwarp=0 ;
201 char *prefix = "NwarpCat" ;
202 mat44 wmat , smat , qmat ;
203 THD_3dim_dataset *oset=NULL ;
204 char *cwarp_all=NULL ; int ntot=0 ;
205
206 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
207
208 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ) NWC_help() ;
209
210 /*-- bureaucracy --*/
211
212 mainENTRY("3dNwarpCat"); machdep();
213 AFNI_logger("3dNwarpCat",argc,argv);
214 PRINT_VERSION("3dNwarpCat"); AUTHOR("Zhark the Warper");
215 (void)COX_clock_time() ;
216 putenv("AFNI_WSINC5_SILENT=YES") ;
217
218 /*-- initialization --*/
219
220 CW_no_expad = 1 ; /* don't allow automatic padding of input warp */
221 Hverb = 0 ; /* don't be verbose inside mri_nwarp.c */
222 for( ii=0 ; ii < NWMAX ; ii++ ) cwarp[ii] = NULL ;
223
224 /*-- scan args --*/
225
226 while( iarg < argc && argv[iarg][0] == '-' ){
227
228 /*---------------*/
229
230 if( strcasecmp(argv[iarg],"-iwarp") == 0 ){
231 do_iwarp = 1 ; iarg++ ; continue ;
232 }
233
234 /*---------------*/
235
236 if( strcasecmp(argv[iarg],"-space") == 0 ){
237 sname = strdup(argv[++iarg]) ; iarg++ ; continue ;
238 }
239
240 /*---------------*/
241
242 if( strcasecmp(argv[iarg],"-NN") == 0 || strncasecmp(argv[iarg],"-nearest",6) == 0 ){
243 WARNING_message("NN interpolation not legal here -- switched to linear") ;
244 interp_code = MRI_LINEAR ; iarg++ ; continue ;
245 }
246 if( strncasecmp(argv[iarg],"-linear",4)==0 || strncasecmp(argv[iarg],"-trilinear",6)==0 ){
247 interp_code = MRI_LINEAR ; iarg++ ; continue ;
248 }
249 if( strncasecmp(argv[iarg],"-cubic",4)==0 || strncasecmp(argv[iarg],"-tricubic",6)==0 ){
250 WARNING_message("cubic interplation not legal here -- switched to quintic") ;
251 interp_code = MRI_QUINTIC ; iarg++ ; continue ;
252 }
253 if( strncasecmp(argv[iarg],"-quintic",4)==0 || strncasecmp(argv[iarg],"-triquintic",6)==0 ){
254 interp_code = MRI_QUINTIC ; iarg++ ; continue ;
255 }
256 if( strncasecmp(argv[iarg],"-wsinc",5) == 0 ){
257 interp_code = MRI_WSINC5 ; iarg++ ; continue ;
258 }
259
260 /*---------------*/
261
262 if( strcasecmp(argv[iarg],"-expad") == 0 ){
263 int expad ;
264 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
265 expad = (int)strtod(argv[iarg],NULL) ;
266 if( expad < 0 ){
267 WARNING_message("-expad %d is illegal and is set to zero",expad) ;
268 expad = 0 ;
269 }
270 CW_extra_pad = expad ; /* this is how we force extra padding */
271 iarg++ ; continue ;
272 }
273
274 /*---------------*/
275
276 if( strncasecmp(argv[iarg],"-interp",5)==0 ){
277 char *inam ;
278 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
279 inam = argv[iarg] ; if( *inam == '-' ) inam++ ;
280 if( strcasecmp(inam,"NN")==0 || strncasecmp(inam,"nearest",5)==0 ){
281 WARNING_message("NN interpolation not legal here -- changed to linear") ;
282 interp_code = MRI_LINEAR ;
283 } else if( strncasecmp(inam,"linear",3)==0 || strncasecmp(inam,"trilinear",5)==0 ){
284 interp_code = MRI_LINEAR ;
285 } else if( strncasecmp(inam,"cubic",3)==0 || strncasecmp(inam,"tricubic",5)==0 ){
286 WARNING_message("cubic interplation not legal here -- changed to quintic") ;
287 interp_code = MRI_QUINTIC ;
288 } else if( strncasecmp(inam,"quintic",3)==0 || strncasecmp(inam,"triquintic",5)==0 ){
289 interp_code = MRI_QUINTIC ;
290 } else if( strncasecmp(inam,"wsinc",4)==0 ){
291 interp_code = MRI_WSINC5 ;
292 } else {
293 ERROR_exit("Unknown code '%s' after '%s' :-(",argv[iarg],argv[iarg-1]) ;
294 }
295 iarg++ ; continue ;
296 }
297
298 /*---------------*/
299
300 if( strcasecmp(argv[iarg],"-verb") == 0 ){
301 verb++ ; NwarpCalcRPN_verb(verb) ; iarg++ ; continue ;
302 }
303
304 /*---------------*/
305
306 if( strcasecmp(argv[iarg],"-prefix") == 0 ){
307 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
308 prefix = argv[iarg] ;
309 if( !THD_filename_ok(prefix) ) ERROR_exit("Illegal name after '%s'",argv[iarg-1]) ;
310 iarg++ ; continue ;
311 }
312
313 /*---------------*/
314
315 if( strncasecmp(argv[iarg],"-warp",5) == 0 ){
316 int nn ;
317 if( iarg >= argc-1 ) ERROR_exit("no argument after '%s' :-(",argv[iarg]) ;
318 if( !isdigit(argv[iarg][5]) ) ERROR_exit("illegal format for '%s' :-(",argv[iarg]) ;
319 nn = (int)strtod(argv[iarg]+5,NULL) ;
320 if( nn <= 0 || nn > NWMAX )
321 ERROR_exit("illegal warp index in '%s' :-(",argv[iarg]) ;
322 if( cwarp[nn-1] != NULL )
323 ERROR_exit("'%s': you can't specify warp #%d more than once :-(",argv[iarg],nn) ;
324 cwarp[nn-1] = strdup(argv[++iarg]) ;
325 if( nn > nwtop ) nwtop = nn ;
326 iarg++ ; continue ;
327 }
328
329 /*---------------*/
330
331 ERROR_message("Confusingly Unknown option '%s' :-(",argv[iarg]) ;
332 suggest_best_prog_option(argv[0],argv[iarg]) ;
333 exit(1) ;
334
335 }
336
337 /*-- load any warps left on the command line, after options --*/
338
339 for( ; iarg < argc && nwtop < NWMAX-1 ; iarg++ )
340 cwarp[nwtop++] = strdup(argv[iarg]) ;
341
342 /*-- check if all warp strings are affine matrices --*/
343
344 #undef AFFINE_WARP_STRING
345 #define AFFINE_WARP_STRING(ss) \
346 ( strstr((ss)," ") == NULL && \
347 ( strcasestr((ss),".1D") != NULL || strcasestr((ss),".txt") != NULL ) )
348
349 for( ntot=ii=0 ; ii < nwtop ; ii++ ){
350 if( cwarp[ii] == NULL ) continue ;
351 ntot += strlen(cwarp[ii]) ;
352 if( ! AFFINE_WARP_STRING(cwarp[ii]) ) break ; /* not affine */
353 }
354 if( ntot == 0 ) ERROR_exit("No warps on command line?!") ;
355
356 if( ii == nwtop ){ /* all are affine (this is for Ziad) */
357 char *fname = malloc(sizeof(char)*(strlen(prefix)+16)) ; FILE *fp ;
358 float a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34 ;
359
360 LOAD_IDENT_MAT44(wmat) ;
361 for( ii=0 ; ii < nwtop ; ii++ ){
362 if( cwarp[ii] == NULL ) continue ;
363 smat = CW_read_affine_warp_OLD(cwarp[ii]) ;
364 qmat = MAT44_MUL(smat,wmat) ; wmat = qmat ;
365 }
366
367 if( strcmp(prefix,"-") == 0 || strncmp(prefix,"stdout",6) == 0 ){
368 fp = stdout ; strcpy(fname,"stdout") ;
369 } else {
370 strcpy(fname,prefix) ;
371 if( strstr(fname,".1D") == NULL ) strcat(fname,".aff12.1D") ;
372 fp = fopen(fname,"w") ;
373 if( fp == NULL ) ERROR_exit("Can't open output matrix file %s",fname) ;
374 }
375 if( do_iwarp ){
376 qmat = MAT44_INV(wmat) ; wmat = qmat ;
377 }
378 UNLOAD_MAT44(wmat,a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34) ;
379 fprintf(fp,
380 " %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
381 a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34 ) ;
382 if( verb && fp != stdout ) INFO_message("Wrote matrix to %s",fname) ;
383 if( fp != stdout ) fclose(fp) ;
384 exit(0) ;
385 }
386
387 /*** at least one nonlinear warp ==> cat all strings, use library function to read ***/
388
389 cwarp_all = (char *)calloc(sizeof(char),(ntot+NWMAX)*2) ;
390 for( ii=0 ; ii < nwtop ; ii++ ){
391 if( cwarp[ii] != NULL ){ strcat(cwarp_all,cwarp[ii]) ; strcat(cwarp_all," ") ; }
392 }
393
394 oset = IW3D_read_catenated_warp( cwarp_all ) ; /* process all of them at once */
395
396 if( do_iwarp ){ /* 18 Jul 2014 */
397 THD_3dim_dataset *qwarp ;
398 if( verb ) fprintf(stderr,"Applying -iwarp option") ;
399 qwarp = THD_nwarp_invert(oset) ;
400 DSET_delete(oset) ;
401 oset = qwarp ;
402 if( verb ) fprintf(stderr,"\n") ;
403 }
404 tross_Make_History( "3dNwarpCat" , argc,argv , oset ) ;
405 if( sname != NULL ) MCW_strncpy( oset->atlas_space , sname , THD_MAX_NAME ) ;
406 EDIT_dset_items( oset , ADN_prefix,prefix , ADN_none ) ;
407 DSET_write(oset) ; WROTE_DSET(oset) ;
408
409 /*--- run away screaming into the night, never to be seen again ---*/
410
411 INFO_message("total CPU time = %.1f sec Elapsed = %.1f\n",
412 COX_cpu_time() , COX_clock_time() ) ;
413
414 exit(0) ;
415 }
416