1 #include "mrilib.h"
2 #include "thd.h"
3 #define MAIN
4
5 /*****************************************************************************
6 Major portions of this software are copyrighted by the Medical College
7 of Wisconsin, 1994-2000, and are released under the Gnu General Public
8 License, Version 2. See the file README.Copyright for details.
9 ******************************************************************************/
10
11 /*--- This program will read in a dataset and write it out in axial order ---*/
12
13 /*** 06 Mar 2000: allow sagittal and coronal as well ***/
14
15 static RwcBoolean write_output = False; /* 08 Aug 2006 [dg] -force rewrite as in 3drefit by rickr */
16 static RwcBoolean NIFTI_mode = False; /* saving NIFTI output */
17 static int cmode = COMPRESS_NOFILE; /* check compression mode for NIFTI separately */
18 static int AXIAL_frugal = 0 ; /* Saving sub-bricks one at a time is optional 18 Mar 2009 [dg] */
19
20
main(int argc,char * argv[])21 int main( int argc , char * argv[] )
22 {
23 THD_3dim_dataset * old_dset , * new_dset ;
24 FD_brick ** fbr , * brax ;
25 int iarg , a1,a2,a3 , aa1,aa2,aa3 , ival,kk,npix,dsiz,code ;
26 THD_ivec3 iv_nxyz , iv_xyzorient ;
27 THD_fvec3 fv_xyzorg , fv_xyzdel ;
28 float xyz_org[4] , xyz_del[4] , brfac_save ;
29 MRI_IMAGE * im ;
30 void * imar ;
31 FILE * data_file=NULL ;
32 THD_datablock * old_dblk , * new_dblk ;
33 char new_prefix[THD_MAX_PREFIX] = "axialize" ;
34 int verbose = 0 , nim , pim=2 ;
35 int native_order , save_order ; /* 23 Nov 1999 */
36
37 int axord=0 ; /* 06 Mar 2000 */
38 char orients[4]="\0" ; /* 07 Dec 2001 */
39
40 MRI_IMARR *im_array=NULL; /* 12 Mar 2009 */
41 MRI_IMAGE *svol_im=NULL;
42 void * svol=NULL ;
43
44 /*- sanity check -*/
45 memset(&iv_xyzorient, 0, sizeof(THD_ivec3));
46
47 WARNING_message("This program (3daxialize) is old, not maintained, and probably useless!") ;
48
49 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
50 printf("Usage: 3daxialize [options] dataset\n"
51 "Purpose: Read in a dataset and write it out as a new dataset\n"
52 " with the data brick oriented as axial slices.\n"
53 " The input dataset must have a .BRIK file.\n"
54 " One application is to create a dataset that can\n"
55 " be used with the AFNI volume rendering plugin.\n"
56 "\n"
57 "Options:\n"
58 " -prefix ppp = Use 'ppp' as the prefix for the new dataset.\n"
59 " [default = 'axialize']\n"
60 " -verb = Print out a progress report.\n"
61 "\n"
62 "The following options determine the order/orientation\n"
63 "in which the slices will be written to the dataset:\n"
64 " -sagittal = Do sagittal slice order [-orient ASL]\n"
65 " -coronal = Do coronal slice order [-orient RSA]\n"
66 " -axial = Do axial slice order [-orient RAI]\n"
67 " This is the default AFNI axial order, and\n"
68 " is the one currently required by the\n"
69 " volume rendering plugin; this is also\n"
70 " the default orientation output by this\n"
71 " program (hence the program's name).\n"
72 "\n"
73 " -orient code = Orientation code for output.\n"
74 " The code must be 3 letters, one each from the\n"
75 " pairs {R,L} {A,P} {I,S}. The first letter gives\n"
76 " the orientation of the x-axis, the second the\n"
77 " orientation of the y-axis, the third the z-axis:\n"
78 " R = Right-to-left L = Left-to-right\n"
79 " A = Anterior-to-posterior P = Posterior-to-anterior\n"
80 " I = Inferior-to-superior S = Superior-to-inferior\n"
81 " If you give an illegal code (e.g., 'LPR'), then\n"
82 " the program will print a message and stop.\n"
83 " N.B.: 'Neurological order' is -orient LPI\n"
84 " -frugal = Write out data as it is rotated, a sub-brick at\n"
85 " a time. This saves a little memory and was the\n"
86 " previous behavior.\n"
87 " Note the frugal option is not available with NIFTI\n"
88 " datasets\n"
89 ) ;
90
91 printf("\n" MASTER_SHORTHELP_STRING ) ;
92 PRINT_COMPILE_DATE ; exit(0) ;
93 }
94
95 mainENTRY("3daxialize main"); machdep(); AFNI_logger("3daxialize",argc,argv);
96 PRINT_VERSION("3daxialize") ;
97
98 /*- scan options -*/
99
100 iarg = 1 ;
101 while( argv[iarg][0] == '-' ){
102
103 if( strcmp(argv[iarg],"-orient") == 0 ){ /* 07 Dec 2001 */
104 int xx,yy,zz ;
105 MCW_strncpy(orients,argv[++iarg],4) ;
106 if( strlen(orients) != 3 ){
107 ERROR_exit("Bad code after -orient: not 3 characters long");
108 }
109 xx = ORCODE(orients[0]) ;
110 yy = ORCODE(orients[1]) ; zz = ORCODE(orients[2]) ;
111 if( xx < 0 || yy < 0 || zz < 0 ){
112 ERROR_exit("Bad code after -orient: illegal characters");
113 }
114 if( !OR3OK(xx,yy,zz) ){
115 ERROR_exit("Bad code after -orient: dependent axes");
116 }
117 axord = -1 ; iarg++ ; continue ;
118 }
119
120 if( strcmp(argv[iarg],"-axial") == 0 ){ /* 06 Mar 2000 */
121 axord = 0 ; iarg++ ; continue ;
122 }
123
124 if( strcmp(argv[iarg],"-sagittal") == 0 ){ /* 06 Mar 2000 */
125 axord = 1 ; iarg++ ; continue ;
126 }
127
128 if( strcmp(argv[iarg],"-coronal") == 0 ){ /* 06 Mar 2000 */
129 axord = 2 ; iarg++ ; continue ;
130 }
131
132 if( strcmp(argv[iarg],"-prefix") == 0 ){
133 iarg++ ;
134 if( iarg >= argc ){
135 ERROR_exit("Need argument after -prefix!") ;
136 }
137 MCW_strncpy( new_prefix , argv[iarg++] , THD_MAX_PREFIX ) ;
138
139 if( strstr(new_prefix,".nii") != NULL ) { /* check for NIFTI mod - drg 08 Aug 2006 */
140 write_output = True;
141 NIFTI_mode = True;
142 if( strstr(new_prefix,".nii.gz") != NULL ) {
143 cmode = 0; /* force gzip compression (actually zlib from nifti library)*/
144 }
145 }
146 else if( !THD_filename_ok(new_prefix) ){
147 ERROR_exit("Illegal new prefix: %s",new_prefix);
148 }
149 continue ;
150 }
151
152 if( strncmp(argv[iarg],"-verbose",5) == 0 ){
153 verbose++ ; iarg++ ; continue ;
154 }
155
156 if( strncmp(argv[iarg],"-frugal",4) == 0 ){
157 AXIAL_frugal = 1; iarg++ ; continue ;
158 }
159
160 ERROR_exit("Unknown option: %s",argv[iarg]);
161 }
162
163 if(NIFTI_mode && AXIAL_frugal){
164 ERROR_message("Frugality and NIFTI output do not mix.\n");
165 ERROR_exit("Try without -frugal or with different output type.\n");
166 }
167
168 /*- get input dataset -*/
169
170 old_dset = THD_open_dataset( argv[iarg] ) ;
171 if( old_dset == NULL ){
172 ERROR_exit("Can't open input dataset: %s",argv[iarg]) ;
173 }
174
175 if( verbose ) INFO_message("Loading input dataset %s",argv[iarg]) ;
176
177 DSET_load(old_dset) ; CHECK_LOAD_ERROR(old_dset) ;
178
179 /*- setup output dataset -*/
180
181 /* use FD bricks for axial, sagittal, coronal displays as basis */
182
183 if( axord >= 0 ){
184 fbr = THD_setup_bricks( old_dset ) ; brax = fbr[axord] ;
185 } else {
186 brax = THD_oriented_brick( old_dset , orients ) ;
187 if( brax == NULL ){
188 ERROR_exit("Can't use -orient code: %s",orients);
189 }
190 }
191
192 new_dset = EDIT_empty_copy( old_dset ) ;
193
194 /* lose obliquity */
195 /* recompute Tc (Cardinal transformation matrix for new grid output */
196 THD_make_cardinal(new_dset);
197
198 tross_Copy_History( old_dset , new_dset ) ;
199 tross_Make_History( "3daxialize" , argc,argv , new_dset ) ;
200
201 /* number of points along each axis */
202
203 LOAD_IVEC3( iv_nxyz , brax->n1 , brax->n2 , brax->n3 ) ;
204
205 /* orientation codes for each axis */
206
207 switch( axord ){ /* 06 Mar 2000 */
208
209 case -1:{
210 int xx=ORCODE(orients[0]) ,
211 yy=ORCODE(orients[1]) , zz=ORCODE(orients[2]) ;
212
213 LOAD_IVEC3( iv_xyzorient , xx,yy,zz ) ;
214 }
215 break ;
216
217 case 0:
218 LOAD_IVEC3( iv_xyzorient, ORI_R2L_TYPE, ORI_A2P_TYPE, ORI_I2S_TYPE ) ;
219 break ;
220
221 case 1:
222 LOAD_IVEC3( iv_xyzorient, ORI_A2P_TYPE, ORI_S2I_TYPE, ORI_L2R_TYPE ) ;
223 break ;
224
225 case 2:
226 LOAD_IVEC3( iv_xyzorient, ORI_R2L_TYPE, ORI_S2I_TYPE, ORI_A2P_TYPE ) ;
227 break ;
228 }
229
230 /* grid spacing for each axis */
231
232 LOAD_FVEC3( fv_xyzdel ,
233 ORIENT_sign[iv_xyzorient.ijk[0]]=='+' ? brax->del1 : -brax->del1,
234 ORIENT_sign[iv_xyzorient.ijk[1]]=='+' ? brax->del2 : -brax->del2,
235 ORIENT_sign[iv_xyzorient.ijk[2]]=='+' ? brax->del3 : -brax->del3);
236
237 UNLOAD_IVEC3( brax->a123 , a1,a2,a3 ) ;
238 aa1 = abs(a1) ; aa2 = abs(a2) ; aa3 = abs(a3) ;
239 xyz_org[1] = new_dset->daxes->xxorg ; xyz_del[1] = new_dset->daxes->xxdel ;
240 xyz_org[2] = new_dset->daxes->yyorg ; xyz_del[2] = new_dset->daxes->yydel ;
241 xyz_org[3] = new_dset->daxes->zzorg ; xyz_del[3] = new_dset->daxes->zzdel ;
242 LOAD_FVEC3( fv_xyzorg ,
243 (a1 > 0) ? xyz_org[aa1] : xyz_org[aa1]+(brax->n1-1)*xyz_del[aa1],
244 (a2 > 0) ? xyz_org[aa2] : xyz_org[aa2]+(brax->n2-1)*xyz_del[aa2],
245 (a3 > 0) ? xyz_org[aa3] : xyz_org[aa3]+(brax->n3-1)*xyz_del[aa3] );
246
247 /*-- open output BRIK file --*/
248 if (cmode == COMPRESS_NOFILE) { /* ignore compression for NIFTI - do in write
249 automatically later */
250 cmode = THD_get_write_compression() ; /* check env. variable for compression*/
251 }
252
253 EDIT_dset_items( new_dset ,
254 ADN_nxyz , iv_nxyz ,
255 ADN_xyzdel , fv_xyzdel ,
256 ADN_xyzorg , fv_xyzorg ,
257 ADN_xyzorient , iv_xyzorient ,
258 ADN_prefix , new_prefix ,
259 ADN_none ) ;
260
261 if( DSET_NUM_TTOFF(new_dset) > 0 )
262 EDIT_dset_items( new_dset , ADN_nsl , 0 , ADN_none ) ; /* 28 Apr 1999 */
263
264 /*- prepare to write new dataset -*/
265
266 new_dblk = new_dset->dblk ;
267 old_dblk = old_dset->dblk ;
268
269 if (AXIAL_frugal)
270 data_file = COMPRESS_fopen_write( new_dblk->diskptr->brick_name, cmode ) ;
271 else
272 INIT_IMARR(im_array);
273
274 npix = brax->n1 * brax->n2 ;
275
276 /*- get slices from input, write to disk -*/
277
278 if( verbose ){
279 INFO_message("Writing new dataset .BRIK"); fflush(stdout);
280 pim = brax->n3 / 5 ; if( pim <= 1 ) pim = 2 ;
281 }
282
283 native_order = mri_short_order() ; /* 23 Nov 1999 */
284 save_order = (DSET_BYTEORDER(new_dset) > 0) ? DSET_BYTEORDER(new_dset)
285 : THD_get_write_order() ;
286
287 for( nim=ival=0 ; ival < DSET_NVALS(old_dset) ; ival++ ){
288 brfac_save = DBLK_BRICK_FACTOR(new_dblk,ival) ;
289 DBLK_BRICK_FACTOR(new_dblk,ival) = 0.0 ;
290 DBLK_BRICK_FACTOR(old_dblk,ival) = 0.0 ;
291 dsiz = mri_datum_size( DSET_BRICK_TYPE(new_dset,ival) ) ;
292
293 if(!AXIAL_frugal) {
294 /* copy the sub-brick volume, svol, into an MRI_IMAGE structure and
295 then append that to the image array */
296 svol_im = mri_new_vol( brax->n1 , brax->n2 , brax->n3, DSET_BRICK_TYPE(new_dset,ival)) ;
297 svol = mri_data_pointer(svol_im);
298 }
299
300
301 for( kk=0 ; kk < brax->n3 ; kk++ ){
302 im = FD_warp_to_mri( kk , ival , brax ) ;
303 imar = mri_data_pointer(im) ;
304 if( save_order != native_order ){ /* 23 Nov 1999 */
305 switch( im->kind ){
306 case MRI_short: mri_swap2( npix,imar) ; break ;
307 case MRI_float:
308 case MRI_int: mri_swap4( npix,imar) ; break ;
309 case MRI_complex: mri_swap4(2*npix,imar) ; break ;
310
311 default:
312 ERROR_exit(
313 "** Illegal input brick type=%s\n",MRI_TYPE_name[im->kind]) ;
314 }
315 }
316
317 if(AXIAL_frugal) {
318 code = fwrite( imar , dsiz , npix , data_file ) ;
319 mri_free(im) ;
320 }
321 else {
322 memcpy((char*)svol+(kk*npix*dsiz), imar, npix*dsiz );
323 }
324
325 if( verbose && (++nim)%pim == 1 ){ printf("."); fflush(stdout); }
326 }
327
328
329 DBLK_BRICK_FACTOR(new_dblk,ival) = brfac_save ;
330 DSET_unload_one(old_dset,ival) ;
331
332 if(!AXIAL_frugal) {
333 /* copy the sub-brick volume, svol, into an MRI_IMAGE structure and
334 then append that to the image array */
335 ADDTO_IMARR(im_array, svol_im);
336 }
337
338 if( verbose ){ printf("!"); fflush(stdout); }
339 }
340 if(AXIAL_frugal)
341 COMPRESS_fclose(data_file) ;
342 else
343 new_dset->dblk->brick = im_array; /* update pointer to data */
344
345 if( verbose ){ printf("\nUnloading old dataset\n") ; fflush(stdout); }
346
347 DSET_unload( old_dset ) ;
348
349 /*- do the output header -*/
350 if(AXIAL_frugal) {
351 DSET_load( new_dset ) ;
352 }
353 else write_output = 1; /* have to write everything out now*/
354 THD_load_statistics( new_dset ) ;
355
356 if( verbose ){ printf("\nWriting new dataset\n") ; fflush(stdout); }
357
358 code = THD_write_3dim_dataset(NULL,NULL,new_dset,write_output); /* (re)write output file */
359 if( !code )
360 ERROR_exit("Could not write dataset");
361
362 if( verbose )
363 INFO_message("Wrote new dataset: %s",DSET_BRIKNAME(new_dset)) ;
364 DSET_unload( new_dset ) ;
365
366 exit(0) ;
367 }
368