1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 #include "mri_genalign_util.c"
8 #include "mri_genalign.c"
9 #include "mri_nwarp.c"
10
11 /*----------------------------------------------------------------------------*/
12 /* This program is mostly a wrapper for function IW3D_load_bsv() */
13
main(int argc,char * argv[])14 int main( int argc , char *argv[] )
15 {
16 THD_3dim_dataset *dset_nwarp=NULL , *dset_out=NULL ;
17 char *prefix = "NFunc" ;
18 int iarg , verb=1 , do_bulk=0, do_shear=0, do_vort=0 , nxyz ;
19 float *bim=NULL , *sim=NULL , *vim=NULL ;
20 IndexWarp3D *AA ;
21
22 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
23
24 /**----------------------------------------------------------------------*/
25 /**----------------- Help the pitifully ignorant user? -----------------**/
26
27 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ){
28 printf(
29 "Usage: 3dNwarpFuncs [options]\n"
30 "\n"
31 "This program reads in a nonlinear 3D warp (from 3dQwarp, etc.) and\n"
32 "computes some functions of the displacements. See the OPTIONS below\n"
33 "for information on what can be computed. The NOTES sections describes\n"
34 "the formulae of the functions that are available.\n"
35 "\n"
36 "--------\n"
37 "OPTIONS:\n"
38 "--------\n"
39 " -nwarp www = 'www' is the name of the 3D warp dataset\n"
40 " (this is a mandatory option!)\n"
41 " ++ This can be computed on the fly, as in 3dNwarpApply.\n"
42 "\n"
43 " -prefix ppp = 'ppp' is the name of the new output dataset\n"
44 "\n"
45 " -bulk = Compute the (fractional) bulk volume change.\n"
46 " ++ e.g., Jacobian determinant minus 1.\n"
47 " ++ see 'MORE...' (below) for interpreting the sign of '-bulk'.\n"
48 " -shear = Compute the shear energy.\n"
49 " -vorticity = Compute the vorticity enerty.\n"
50 " -all = Compute all 3 of these fun fun functions.\n"
51 "\n"
52 "If none of '-bulk', '-shear', or '-vorticity' are given, then '-bulk'\n"
53 "will be assumed.\n"
54 "\n"
55 "------\n"
56 "NOTES:\n"
57 "------\n"
58 "Denote the displacement vector field (warp) by\n"
59 " [ p(x,y,z) , q(x,y,z) , r(x,y,z) ]\n"
60 "Define the Jacobian matrix by\n"
61 "\n"
62 " [ 1+dp/dx dp/dy dp/dz ] [ Jxx Jxy Jxz ]\n"
63 " J = [ dq/dx 1+dq/dy dq/dz ] = [ Jyx Jyy Jyz ]\n"
64 " [ dr/dx dr/dy 1+dr/dz ] [ Jzx Jzy Jzz ]\n"
65 "\n"
66 "* The '-bulk' output is the determinant of this matrix (det[J]), minus 1.\n"
67 "* It measures the fractional amount of volume distortion.\n"
68 "* Negative means the warped coordinates are shrunken (closer together)\n"
69 " than the input coordinates. Also see the 'MORE...' section below.\n"
70 "\n"
71 "* The '-shear' output is the sum of squares of the J matrix elements --\n"
72 " which equals the sum of squares of its eigenvalues -- divided by\n"
73 " det[J]^(2/3), then minus 3.\n"
74 "* It measures the amount of shearing distortion (normalized by the amount\n"
75 " of volume distortion).\n"
76 "\n"
77 "* The '-vorticity' output is the sum of squares of the skew part of\n"
78 " the J matrix = [ Jxy-Jyx , Jxz-Jzx , Jyz-Jzy ], divided by det[J]^(2/3).\n"
79 "* It measures the amount of twisting distortion (also normalized).\n"
80 "\n"
81 "* All 3 of these functions are dimensionless.\n"
82 "\n"
83 "* The penalty used in 3dQwarp is a combination of the bulk, shear,\n"
84 " and vorticity functions.\n"
85 "\n"
86 "------------------------------\n"
87 "MORE about interpreting -bulk:\n"
88 "------------------------------\n"
89 "If the warp N(x,y,z) is the '_WARP' output from 3dQwarp, then N(x,y,z)\n"
90 "maps the base dataset (x,y,z) coordinates to the source dataset (x,y,z)\n"
91 "coordinates. If the source dataset has to expand in size to match\n"
92 "the base dataset, then going from base coordinates to source must\n"
93 "be a shrinkage. Thus, negative '-bulk' in this '_WARP' dataset\n"
94 "corresponds to expansion going from source to base. Conversely,\n"
95 "in this situation, positive '-bulk' will show up in the '_WARPINV'\n"
96 "dataset from 3dQwarp as that is the map from source (x,y,z) to\n"
97 "base (x,y,z).\n"
98 "\n"
99 "The situation above happens a lot when using one of the MNI152 human\n"
100 "brain templates as the base dataset. This family of datasets is larger\n"
101 "than the average human brain, due to the averaging process used to\n"
102 "define the first MNI152 template back in the 1990s.\n"
103 "\n"
104 "I have no easy interpretation handy for the '-shear' and '-vorticity'\n"
105 "outputs, alas. They are computed as part of the penalty function used\n"
106 "to control weirdness in the 3dQwarp optimization process.\n"
107 "\n"
108 "---------------------------\n"
109 "AUTHOR -- RWCox == @AFNIman\n"
110 "---------------------------\n"
111 ) ;
112
113 PRINT_COMPILE_DATE ;
114 exit(0) ;
115 }
116
117 /**--- bookkeeping and marketing ---**/
118
119 mainENTRY("3dNwarpFuncs"); machdep();
120 AFNI_logger("3dNwarpFuncs",argc,argv);
121 PRINT_VERSION("3dNwarpFuncs"); AUTHOR("Zhark the Warped");
122 (void)COX_clock_time() ;
123 putenv("AFNI_WSINC5_SILENT=YES") ;
124
125 /**--- process command line options ---**/
126
127 iarg = 1 ;
128 while( iarg < argc && argv[iarg][0] == '-' ){
129
130 /*---------------*/
131
132 if( strcasecmp(argv[iarg],"-quiet") == 0 ){
133 verb = 0 ; iarg++ ; continue ;
134 }
135 if( strcasecmp(argv[iarg],"-verb") == 0 ){
136 verb++ ; iarg++ ; continue ;
137 }
138
139 /*---------------*/
140
141 if( strncasecmp(argv[iarg],"-prefix",5) == 0 ){
142 if( ++iarg >= argc ) ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
143 if( !THD_filename_ok(argv[iarg]) )
144 ERROR_exit("badly formed filename: '%s' '%s' :-(",argv[iarg-1],argv[iarg]) ;
145 if( strcasecmp(argv[iarg],"NULL") == 0 )
146 ERROR_exit("can't use filename: '%s' '%s' :-(",argv[iarg-1],argv[iarg]) ;
147 else
148 prefix = strdup(argv[iarg]) ;
149 iarg++ ; continue ;
150 }
151
152 /*---------------*/
153
154 if( strcasecmp(argv[iarg],"-nwarp") == 0 || strcasecmp(argv[iarg],"-warp") == 0 ){
155 if( dset_nwarp != NULL ) ERROR_exit("Can't have multiple %s options :-(",argv[iarg]) ;
156 if( ++iarg >= argc ) ERROR_exit("No argument after '%s' :-(",argv[iarg-1]) ;
157 dset_nwarp = IW3D_read_catenated_warp( argv[iarg] ) ;
158 if( dset_nwarp == NULL ) ERROR_exit("can't open warp dataset '%s' :-(",argv[iarg]);
159 if( DSET_NVALS(dset_nwarp) < 3 ) ERROR_exit("dataset '%s' isn't a 3D warp",argv[iarg]);
160 iarg++ ; continue ;
161 }
162
163 /*---------------*/
164
165 if( strncasecmp(argv[iarg],"-bulk",4) == 0 ){
166 do_bulk = 1 ; iarg++ ; continue ;
167 }
168 if( strncasecmp(argv[iarg],"-shear",4) == 0 ){
169 do_shear = 1 ; iarg++ ; continue ;
170 }
171 if( strncasecmp(argv[iarg],"-vort",4) == 0 ){
172 do_vort = 1 ; iarg++ ; continue ;
173 }
174 if( strcasecmp(argv[iarg],"-all") == 0 ) {
175 do_bulk = do_shear = do_vort = 1 ; iarg++ ; continue ;
176 }
177
178 /*---------------*/
179
180 ERROR_message("Bogusly Unknown option '%s' :-( :-( :-(",argv[iarg]) ;
181 suggest_best_prog_option(argv[0],argv[iarg]) ;
182 exit(1) ;
183
184 }
185
186 /*-------- check inputs to see if the user is completely demented ---------*/
187
188 if( dset_nwarp == NULL && iarg >= argc )
189 ERROR_exit("No -nwarp option? How do you want to warp? :-(") ;
190
191 if( dset_nwarp == NULL ){
192 dset_nwarp = IW3D_read_catenated_warp( argv[iarg] ) ;
193 if( dset_nwarp == NULL ) ERROR_exit("can't open warp dataset '%s' :-(",argv[iarg]);
194 if( DSET_NVALS(dset_nwarp) < 3 ) ERROR_exit("dataset '%s' isn't a 3D warp",argv[iarg]);
195 }
196
197 if( do_bulk == 0 && do_shear == 0 && do_vort == 0 ){
198 do_bulk = 1 ;
199 INFO_message("No outputs ordered ==> '-bulk' is what you'll get") ;
200 }
201
202 nxyz = DSET_NVOX(dset_nwarp) ;
203
204 if( do_bulk ) bim = (float *)calloc(sizeof(float),nxyz) ;
205 if( do_shear ) sim = (float *)calloc(sizeof(float),nxyz) ;
206 if( do_vort ) vim = (float *)calloc(sizeof(float),nxyz) ;
207
208 /*--- the actual work (bow your head in reverence) ---*/
209
210 verb_nww = 0 ;
211
212 AA = IW3D_from_dataset( dset_nwarp , 0,0 ) ;
213
214 IW3D_load_bsv( AA , fabsf(DSET_DX(dset_nwarp)) ,
215 fabsf(DSET_DY(dset_nwarp)) ,
216 fabsf(DSET_DZ(dset_nwarp)) , bim,sim,vim ) ;
217
218 DSET_unload(dset_nwarp) ; IW3D_destroy(AA) ;
219
220 dset_out = EDIT_empty_copy(dset_nwarp) ;
221 EDIT_dset_items( dset_out ,
222 ADN_prefix , prefix ,
223 ADN_nvals , do_bulk+do_shear+do_vort ,
224 ADN_ntt , 0 ,
225 ADN_datum_all , MRI_float ,
226 ADN_none ) ;
227 tross_Copy_History( dset_nwarp , dset_out ) ; /* hysterical records */
228 tross_Make_History( "3dNwarpFuncs" , argc,argv , dset_out ) ;
229
230 iarg = 0 ;
231 if( do_bulk ){
232 EDIT_substitute_brick( dset_out , iarg , MRI_float , bim ) ;
233 EDIT_BRICK_LABEL( dset_out , iarg , "bulk" ) ;
234 iarg++ ;
235 }
236 if( do_shear ){
237 EDIT_substitute_brick( dset_out , iarg , MRI_float , sim ) ;
238 EDIT_BRICK_LABEL( dset_out , iarg , "shear" ) ;
239 iarg++ ;
240 }
241 if( do_vort ){
242 EDIT_substitute_brick( dset_out , iarg , MRI_float , vim ) ;
243 EDIT_BRICK_LABEL( dset_out , iarg , "vorticity" ) ;
244 iarg++ ;
245 }
246
247 DSET_write(dset_out) ; if( verb ) WROTE_DSET(dset_out) ;
248
249 exit(0) ;
250 }
251