1 #include "mrilib.h"
2 #include "thd_brainormalize.h"
3
usage_3dSpatNorm(int detail)4 void usage_3dSpatNorm(int detail) {
5 printf(
6 "Usage: 3dSpatNorm [options] dataset\n"
7 "\n"
8 "Options:\n"
9 " -prefix ppp = Write output dataset using 'ppp' for the prefix.\n"
10 " -orig_space = Write output dataset using the same grid as dataset.\n"
11 " -verb = Write out progress reports\n"
12 " -monkey : Monkey business\n"
13 " -marmost: Marmoset head\n"
14 " -rat: Rat head\n"
15 " -human: Bone head (default)\n"
16 " -bottom_cuts CUTFLAGS = Make approximate cuts at the bottom to shave\n"
17 " non brain areas. CUTFLAGS is a string of \n"
18 " characters indicating which sides to cut.\n"
19 " An 'A' cuts along the anterior side\n"
20 " 'P' for posterior, and 'R', 'L' for right, and\n"
21 " left, respectively.\n"
22 " To cut all four, use: -bottom_cuts APLR\n"
23 " Note: -bottom_cuts only works for Human heads.\n"
24 "--------------------------------------------------------------------\n"
25 "* This program is obsolete, and should not be used by most people. *\n"
26 "--------------------------------------------------------------------\n"
27 ) ;
28 PRINT_COMPILE_DATE ;
29 }
30
main(int argc,char * argv[])31 int main( int argc , char *argv[] )
32 {
33 MRI_IMAGE *imin, *imout , *imout_orig;
34 THD_3dim_dataset *iset, *oset , *ooset;
35 char *prefix = "SpatNorm", *bottom_cuts = NULL;
36 int iarg , verb=0, OrigSpace = 0 , specie = HUMAN;
37 float SpatNormDxyz= 0.0, iset_scaled=1.0;
38 THD_ivec3 orixyz , nxyz ;
39 THD_fvec3 dxyz , orgxyz, originRAIfv, fv2;
40
41
42 mainENTRY("3dSpatNorm main") ; machdep() ;
43 if (argc == 1) { usage_3dSpatNorm(1); exit(0); }
44
45 /*--- options ---*/
46
47 iarg = 1 ;
48 OrigSpace = 0;
49 while( iarg < argc && argv[iarg][0] == '-' ){
50 if (strcmp(argv[iarg],"-h") == 0 || strcmp(argv[iarg],"-help") == 0 ) {
51 usage_3dSpatNorm(strlen(argv[iarg]) > 3 ? 2:1);
52 exit(0);
53 }
54
55 /* -prefix */
56
57 if( strcmp(argv[iarg],"-prefix") == 0 ){
58 if( ++iarg >= argc ){
59 fprintf(stderr,"**ERROR: -prefix requires another argument!\n") ;
60 exit(1) ;
61 }
62 prefix = strdup(argv[iarg]) ;
63 if( !THD_filename_ok(prefix) ){
64 fprintf(stderr,"**ERROR: -prefix value contains forbidden characters!\n") ;
65 exit(1) ;
66 }
67 iarg++ ; continue ;
68 }
69
70 if( strcmp(argv[iarg],"-dxyz") == 0 ){
71 if( ++iarg >= argc ){
72 fprintf(stderr,"**ERROR: -dxyz requires another argument!\n") ;
73 exit(1) ;
74 }
75 SpatNormDxyz = atof(argv[iarg]) ;
76
77 iarg++ ; continue ;
78 }
79 if( strcmp(argv[iarg],"-bottom_cuts") == 0 ){
80 if( ++iarg >= argc ){
81 fprintf(stderr,"**ERROR: -bottom_cuts requires another argument!\n") ;
82 exit(1) ;
83 }
84 bottom_cuts = argv[iarg] ;
85
86 iarg++ ; continue ;
87 }
88 if( strncmp(argv[iarg],"-verb",5) == 0 ){
89 verb++ ; iarg++ ; continue ;
90 }
91 if( strncmp(argv[iarg],"-human",5) == 0 ){
92 specie = HUMAN ; iarg++ ; continue ;
93 }
94 if( strncmp(argv[iarg],"-monkey",5) == 0 ){
95 specie = MONKEY ; iarg++ ; continue ;
96 }
97 if( strncmp(argv[iarg],"-marmoset",5) == 0 ){
98 specie = MARMOSET ; iarg++ ; continue ;
99 }
100 if( strncmp(argv[iarg],"-rat",5) == 0 ){
101 specie = RAT ; iarg++ ; continue ;
102 }
103 if( strncmp(argv[iarg],"-orig_space",10) == 0 ){
104 OrigSpace = 1 ; iarg++ ; continue ;
105 }
106
107 fprintf(stderr,"**ERROR: %s is unknown option!\n",argv[iarg]) ;
108 suggest_best_prog_option(argv[0], argv[iarg]);
109 exit(1) ;
110 }
111
112 if( iarg >= argc ){
113 fprintf(stderr,"**ERROR: no input dataset name on command line?!\n") ;
114 exit(1) ;
115 }
116
117 /*--- read dataset ---*/
118
119 iset = THD_open_dataset( argv[iarg] ) ;
120 if( !ISVALID_DSET(iset) ){
121 fprintf(stderr,"**ERROR: can't open dataset %s\n",argv[iarg]) ;
122 exit(1) ;
123 }
124
125 /*--- get median brick --*/
126
127 if( verb ) fprintf(stderr,"++3dSpatNorm: loading dataset\n") ;
128
129 if (specie == MARMOSET) {
130 iset_scaled = 2.5;
131 THD_volDXYZscale(iset->daxes, iset_scaled, 0);
132 specie = MONKEY;
133 }
134 imin = THD_median_brick( iset ) ;
135 if( imin == NULL ){
136 fprintf(stderr,"**ERROR: can't load dataset %s\n",argv[iarg]) ;
137 exit(1) ;
138 }
139 imin->dx = fabs(iset->daxes->xxdel) ;
140 imin->dy = fabs(iset->daxes->yydel) ;
141 imin->dz = fabs(iset->daxes->zzdel) ;
142
143
144 mri_speciebusiness(specie);
145 mri_brain_normalize_cuts(bottom_cuts);
146
147 if (SpatNormDxyz) {
148 if (verb) fprintf(stderr,"Overriding default resampling\n");
149 mri_brainormalize_initialize(SpatNormDxyz, SpatNormDxyz, SpatNormDxyz);
150 } else {
151 float xxdel, yydel, zzdel, minres;
152 if (specie == MONKEY) minres = 0.5;
153 else if (specie == MARMOSET) minres = 0.2;
154 else if (specie == RAT) minres = 0.1;
155 else minres = 0.5;
156 /* don't allow for too low a resolution, please */
157 if (imin->dx < minres) xxdel = minres;
158 else xxdel = imin->dx;
159 if (imin->dy < minres) yydel = minres;
160 else yydel = imin->dy;
161 if (imin->dz < minres) zzdel = minres;
162 else zzdel = imin->dz;
163 if (verb) {
164 fprintf(stderr,
165 "%s:\n"
166 " Original resolution %f, %f, %f\n"
167 " SpatNorm resolution %f, %f, %f\n",
168 "3dSpatnorm", imin->dx, imin->dy, imin->dz,
169 xxdel, yydel, zzdel);
170 }
171 mri_brainormalize_initialize(xxdel, yydel, zzdel);
172 }
173
174 /* To get around the #define for voxel counts and dimensions */
175 mri_brainormalize_initialize(imin->dz, imin->dy, imin->dz);
176
177 /* me needs the origin of this dset in RAI world */
178 LOAD_FVEC3( originRAIfv ,
179 iset->daxes->xxorg , iset->daxes->yyorg , iset->daxes->zzorg) ;
180 originRAIfv = THD_3dmm_to_dicomm( iset , originRAIfv ) ;
181
182 LOAD_FVEC3(fv2, iset->daxes->xxorg + (iset->daxes->nxx-1)*iset->daxes->xxdel ,
183 iset->daxes->yyorg + (iset->daxes->nyy-1)*iset->daxes->yydel ,
184 iset->daxes->zzorg + (iset->daxes->nzz-1)*iset->daxes->zzdel);
185 fv2 = THD_3dmm_to_dicomm( iset , fv2 ) ;
186
187 if( originRAIfv.xyz[0] > fv2.xyz[0] ) {
188 float tf; tf = originRAIfv.xyz[0];
189 originRAIfv.xyz[0] = fv2.xyz[0]; fv2.xyz[0] = tf; }
190 if( originRAIfv.xyz[1] > fv2.xyz[1] ) {
191 float tf; tf = originRAIfv.xyz[1];
192 originRAIfv.xyz[1] = fv2.xyz[1]; fv2.xyz[1] = tf; }
193 if( originRAIfv.xyz[2] > fv2.xyz[2] ) {
194 float tf; tf = originRAIfv.xyz[2];
195 originRAIfv.xyz[2] = fv2.xyz[2]; fv2.xyz[2] = tf; }
196
197 if (verb) {
198 fprintf(stderr,"++3dSpatNorm (ZSS): RAI origin info: %f %f %f\n",
199 originRAIfv.xyz[0], originRAIfv.xyz[1], originRAIfv.xyz[2]);
200 }
201
202
203 DSET_unload( iset ) ; /* don't need this data no more */
204
205 /*-- convert image to shorts, if appropriate --*/
206
207 if( DSET_BRICK_TYPE(iset,0) == MRI_short ||
208 DSET_BRICK_TYPE(iset,0) == MRI_byte ){
209
210 imout = mri_to_short(0.0,imin) ; /* ZSS Oct 2012: Let function set scaling*/
211 mri_free(imin) ; imin = imout ;
212 }
213
214 /*--- normalize image spatially ---*/
215
216 mri_brainormalize_verbose( verb ) ;
217 if (OrigSpace) {
218 imout = mri_brainormalize( imin , iset->daxes->xxorient,
219 iset->daxes->yyorient,
220 iset->daxes->zzorient , &imout_orig, NULL) ;
221 } else {
222 imout = mri_brainormalize( imin , iset->daxes->xxorient,
223 iset->daxes->yyorient,
224 iset->daxes->zzorient , NULL, NULL) ;
225 }
226 mri_free( imin ) ;
227
228 if( imout == NULL ){
229 fprintf(stderr,"**ERROR: normalization fails!?\n"); exit(1);
230 }
231
232 if (OrigSpace) {
233 if( verb ) fprintf(stderr,"++3dSpatNorm: Output in Orignal space\n") ;
234 mri_free( imout ) ;
235 imout = imout_orig;
236 imout->xo = originRAIfv.xyz[0];
237 imout->yo = originRAIfv.xyz[1];
238 imout->zo = originRAIfv.xyz[2];
239 imout_orig = NULL;
240 } else {
241 if( verb ) fprintf(stderr,"++3dSpatNorm: Output in SpatNorm space\n") ;
242 }
243
244 #if 0
245 if( AFNI_yesenv("WATERSHED") ){
246 imin = mri_watershedize( imout , 0.10 ) ;
247 if( imin != NULL ){ mri_free(imout); imout = imin; }
248 }
249 #endif
250
251 /*--- create output dataset ---*/
252 if( verb )
253 fprintf(stderr,"++3dSpatNorm: Creating output dset\n") ;
254
255 oset = EDIT_empty_copy( NULL ) ;
256
257 tross_Copy_History( iset , oset ) ;
258 tross_Make_History( "3dSpatNorm" , argc,argv , oset ) ;
259
260 LOAD_IVEC3( nxyz , imout->nx , imout->ny , imout->nz ) ;
261 LOAD_FVEC3( dxyz , imout->dx , imout->dy , imout->dz ) ;
262 LOAD_FVEC3( orgxyz , imout->xo , imout->yo , imout->zo ) ;
263 LOAD_IVEC3( orixyz , ORI_R2L_TYPE , ORI_A2P_TYPE , ORI_I2S_TYPE ) ;
264
265 if( verb )
266 fprintf(stderr,"++3dSpatNorm: EDIT_dset_items\n") ;
267 EDIT_dset_items( oset ,
268 ADN_prefix , prefix ,
269 ADN_datum_all , imout->kind ,
270 ADN_nxyz , nxyz ,
271 ADN_xyzdel , dxyz ,
272 ADN_xyzorg , orgxyz ,
273 ADN_xyzorient , orixyz ,
274 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
275 ADN_view_type , VIEW_ORIGINAL_TYPE ,
276 ADN_type , HEAD_ANAT_TYPE ,
277 ADN_func_type , ANAT_BUCK_TYPE ,
278 ADN_none ) ;
279
280 if( verb )
281 fprintf(stderr,"++3dSpatNorm: EDIT_substitute_brick\n") ;
282 EDIT_substitute_brick( oset , 0 , imout->kind , mri_data_pointer(imout) ) ;
283
284 if (OrigSpace) {
285 if( verb )
286 fprintf(stderr,"++3dSpatNorm: Changing orientation from RAI\n") ;
287 ooset = r_new_resam_dset ( oset, iset, 0, 0, 0, NULL, MRI_NN, NULL, 1, 0);
288 if (!ooset) {
289 fprintf(stderr,"**ERROR: Failed to reslice!?\n"); exit(1);
290 }
291 /* put prefix back, r_new_resam_dset puts dummy prefix */
292 EDIT_dset_items( ooset ,
293 ADN_prefix , prefix,
294 ADN_none ) ;
295
296 DSET_delete(oset); oset = ooset; ooset = NULL;
297 }
298
299 if (iset_scaled != 1.0f)
300 THD_volDXYZscale(oset->daxes,
301 1/iset_scaled, 0);
302
303 DSET_write(oset) ;
304 if( verb )
305 fprintf(stderr,"++3dSpatNorm: wrote dataset %s\n",DSET_BRIKNAME(oset)) ;
306
307 exit(0) ;
308 }
309