00001 #include "mrilib.h"
00002 #include "thd.h"
00003 #define MAIN
00004 
00005 /*****************************************************************************
00006    Major portions of this software are copyrighted by the Medical College
00007    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00008    License, Version 2.  See the file README.Copyright for details.
00009 ******************************************************************************/
00010 
00011 /*--- This program will read in a dataset and write it out in axial order ---*/
00012 
00013 /*** 06 Mar 2000: allow sagittal and coronal as well ***/
00014 
00015 int main( int argc , char * argv[] )
00016 {
00017    THD_3dim_dataset * old_dset , * new_dset ;
00018    FD_brick ** fbr , * brax ;
00019    int iarg , a1,a2,a3 , aa1,aa2,aa3 , ival,kk,cmode,npix,dsiz,code ;
00020    THD_ivec3 iv_nxyz   , iv_xyzorient ;
00021    THD_fvec3 fv_xyzorg , fv_xyzdel ;
00022    float xyz_org[4] , xyz_del[4] , brfac_save ;
00023    MRI_IMAGE * im ;
00024    void * imar ;
00025    FILE * far ;
00026    THD_datablock * old_dblk , * new_dblk ;
00027    char new_prefix[THD_MAX_PREFIX] = "axialize" ;
00028    int verbose = 0 , nim , pim=2 ;
00029    int native_order , save_order ;  /* 23 Nov 1999 */
00030 
00031    int axord=0 ;          /* 06 Mar 2000 */
00032    char orients[4]="\0" ; /* 07 Dec 2001 */
00033 
00034    /*- sanity check -*/
00035 
00036    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00037       printf("Usage: 3daxialize [options] dataset\n"
00038              "Purpose: Read in a dataset and write it out as a new dataset\n"
00039              "         with the data brick oriented as axial slices.\n"
00040              "         The input dataset must have a .BRIK file.\n"
00041              "         One application is to create a dataset that can\n"
00042              "         be used with the AFNI volume rendering plugin.\n"
00043              "\n"
00044              "Options:\n"
00045              " -prefix ppp  = Use 'ppp' as the prefix for the new dataset.\n"
00046              "               [default = 'axialize']\n"
00047              " -verb        = Print out a progress report.\n"
00048              "\n"
00049              "The following options determine the order/orientation\n"
00050              "in which the slices will be written to the dataset:\n"
00051              " -sagittal    = Do sagittal slice order [-orient ASL]\n"
00052              " -coronal     = Do coronal slice order  [-orient RSA]\n"
00053              " -axial       = Do axial slice order    [-orient RAI]\n"
00054              "                 This is the default AFNI axial order, and\n"
00055              "                 is the one currently required by the\n"
00056              "                 volume rendering plugin; this is also\n"
00057              "                 the default orientation output by this\n"
00058              "                 program (hence the program's name).\n"
00059              "\n"
00060              " -orient code = Orientation code for output.\n"
00061              "                The code must be 3 letters, one each from the\n"
00062              "                pairs {R,L} {A,P} {I,S}.  The first letter gives\n"
00063              "                the orientation of the x-axis, the second the\n"
00064              "                orientation of the y-axis, the third the z-axis:\n"
00065              "                 R = Right-to-left         L = Left-to-right\n"
00066              "                 A = Anterior-to-posterior P = Posterior-to-anterior\n"
00067              "                 I = Inferior-to-superior  S = Superior-to-inferior\n"
00068              "                If you give an illegal code (e.g., 'LPR'), then\n"
00069              "                the program will print a message and stop.\n"
00070              "          N.B.: 'Neurological order' is -orient LPI\n"
00071             ) ;
00072 
00073       printf("\n" MASTER_SHORTHELP_STRING ) ;
00074       exit(0) ;
00075    }
00076 
00077    mainENTRY("3daxialize main"); machdep(); AFNI_logger("3daxialize",argc,argv);
00078    PRINT_VERSION("3daxialize") ;
00079 
00080    /*- scan options -*/
00081 
00082    iarg = 1 ;
00083    while( argv[iarg][0] == '-' ){
00084 
00085       if( strcmp(argv[iarg],"-orient") == 0 ){    /* 07 Dec 2001 */
00086          int xx,yy,zz ;
00087          MCW_strncpy(orients,argv[++iarg],4) ;
00088          if( strlen(orients) != 3 ){
00089            fprintf(stderr,"** Bad code after -orient: not 3 characters long\n");
00090            exit(1);
00091          }
00092          xx = ORCODE(orients[0]) ;
00093          yy = ORCODE(orients[1]) ; zz = ORCODE(orients[2]) ;
00094          if( xx < 0 || yy < 0 || zz < 0 ){
00095            fprintf(stderr,"** Bad code after -orient: illegal characters\n");
00096            exit(1);
00097          }
00098          if( !OR3OK(xx,yy,zz) ){
00099            fprintf(stderr,"** Bad code after -orient: dependent axes\n");
00100            exit(1);
00101          }
00102          axord = -1 ; iarg++ ; continue ;
00103       }
00104 
00105       if( strcmp(argv[iarg],"-axial") == 0 ){     /* 06 Mar 2000 */
00106          axord = 0 ; iarg++ ; continue ;
00107       }
00108 
00109       if( strcmp(argv[iarg],"-sagittal") == 0 ){  /* 06 Mar 2000 */
00110          axord = 1 ; iarg++ ; continue ;
00111       }
00112 
00113       if( strcmp(argv[iarg],"-coronal") == 0 ){   /* 06 Mar 2000 */
00114          axord = 2 ; iarg++ ; continue ;
00115       }
00116 
00117       if( strcmp(argv[iarg],"-prefix") == 0 ){
00118          strcpy( new_prefix , argv[++iarg] ) ;
00119          if( !THD_filename_ok(new_prefix) ){
00120             fprintf(stderr,"** illegal new prefix: %s\n",new_prefix); exit(1);
00121          }
00122          iarg++ ; continue ;
00123       }
00124 
00125       if( strncmp(argv[iarg],"-verbose",5) == 0 ){
00126          verbose++ ; iarg++ ; continue ;
00127       }
00128 
00129       fprintf(stderr,"** Unknown option: %s\n",argv[iarg]); exit(1);
00130    }
00131 
00132    /*- get input dataset -*/
00133 
00134    old_dset = THD_open_dataset( argv[iarg] ) ;
00135    if( old_dset == NULL ){
00136       fprintf(stderr,"** can't open input dataset: %s\n",argv[iarg]) ; exit(1) ;
00137    }
00138 
00139    if( verbose ) printf("++ Loading input dataset %s\n",argv[iarg]) ;
00140 
00141    DSET_load(old_dset) ;
00142    if( !DSET_LOADED(old_dset) ){
00143       fprintf(stderr,"** can't load input .BRIK: %s\n",argv[iarg]); exit(1);
00144    }
00145 
00146    /*- setup output dataset -*/
00147 
00148    /* use FD bricks for axial, sagittal, coronal displays as basis */
00149 
00150    if( axord >= 0 ){
00151       fbr = THD_setup_bricks( old_dset ) ; brax = fbr[axord] ;
00152    } else {
00153       brax = THD_oriented_brick( old_dset , orients ) ;
00154       if( brax == NULL ){
00155          fprintf(stderr,"** Can't use -orient code: %s\n",orients); exit(1);
00156       }
00157    }
00158 
00159    new_dset = EDIT_empty_copy( old_dset ) ;
00160 
00161    tross_Copy_History( old_dset , new_dset ) ;
00162    tross_Make_History( "3daxialize" , argc,argv , new_dset ) ;
00163 
00164    /* number of points along each axis */
00165 
00166    LOAD_IVEC3( iv_nxyz , brax->n1 , brax->n2 , brax->n3 ) ;
00167 
00168    /* orientation codes for each axis */
00169 
00170    switch( axord ){  /* 06 Mar 2000 */
00171 
00172       case -1:{
00173          int xx=ORCODE(orients[0]) ,
00174              yy=ORCODE(orients[1]) , zz=ORCODE(orients[2]) ;
00175 
00176          LOAD_IVEC3( iv_xyzorient , xx,yy,zz ) ;
00177       }
00178       break ;
00179 
00180       case 0:
00181        LOAD_IVEC3( iv_xyzorient, ORI_R2L_TYPE, ORI_A2P_TYPE, ORI_I2S_TYPE ) ;
00182       break ;
00183 
00184       case 1:
00185        LOAD_IVEC3( iv_xyzorient, ORI_A2P_TYPE, ORI_S2I_TYPE, ORI_L2R_TYPE ) ;
00186       break ;
00187 
00188       case 2:
00189        LOAD_IVEC3( iv_xyzorient, ORI_R2L_TYPE, ORI_S2I_TYPE, ORI_A2P_TYPE ) ;
00190       break ;
00191    }
00192 
00193    /* grid spacing for each axis */
00194 
00195    LOAD_FVEC3( fv_xyzdel ,
00196                ORIENT_sign[iv_xyzorient.ijk[0]]=='+' ? brax->del1 : -brax->del1,
00197                ORIENT_sign[iv_xyzorient.ijk[1]]=='+' ? brax->del2 : -brax->del2,
00198                ORIENT_sign[iv_xyzorient.ijk[2]]=='+' ? brax->del3 : -brax->del3);
00199 
00200    UNLOAD_IVEC3( brax->a123 , a1,a2,a3 ) ;
00201    aa1 = abs(a1) ; aa2 = abs(a2) ; aa3 = abs(a3) ;
00202    xyz_org[1] = new_dset->daxes->xxorg ; xyz_del[1] = new_dset->daxes->xxdel ;
00203    xyz_org[2] = new_dset->daxes->yyorg ; xyz_del[2] = new_dset->daxes->yydel ;
00204    xyz_org[3] = new_dset->daxes->zzorg ; xyz_del[3] = new_dset->daxes->zzdel ;
00205    LOAD_FVEC3( fv_xyzorg ,
00206                (a1 > 0) ? xyz_org[aa1] : xyz_org[aa1]+(brax->n1-1)*xyz_del[aa1],
00207                (a2 > 0) ? xyz_org[aa2] : xyz_org[aa2]+(brax->n2-1)*xyz_del[aa2],
00208                (a3 > 0) ? xyz_org[aa3] : xyz_org[aa3]+(brax->n3-1)*xyz_del[aa3] );
00209 
00210    EDIT_dset_items( new_dset ,
00211                        ADN_nxyz      , iv_nxyz ,
00212                        ADN_xyzdel    , fv_xyzdel ,
00213                        ADN_xyzorg    , fv_xyzorg ,
00214                        ADN_xyzorient , iv_xyzorient ,
00215                        ADN_prefix    , new_prefix ,
00216                     ADN_none ) ;
00217 
00218    if( DSET_NUM_TTOFF(new_dset) > 0 )
00219       EDIT_dset_items( new_dset , ADN_nsl , 0 , ADN_none ) ;  /* 28 Apr 1999 */
00220 
00221    /*- prepare to write new dataset -*/
00222 
00223    new_dblk = new_dset->dblk ;
00224    old_dblk = old_dset->dblk ;
00225 
00226    cmode = THD_get_write_compression() ;
00227    far   = COMPRESS_fopen_write( new_dblk->diskptr->brick_name, cmode ) ;
00228    npix  = brax->n1 * brax->n2 ;
00229 
00230    /*- get slices from input, write to disk -*/
00231 
00232    if( verbose ){
00233       printf("++ Writing new dataset .BRIK"); fflush(stdout);
00234       pim = brax->n3 / 5 ; if( pim <= 1 ) pim = 2 ;
00235    }
00236 
00237    native_order = mri_short_order() ;                           /* 23 Nov 1999 */
00238    save_order   = (DSET_BYTEORDER(new_dset) > 0) ? DSET_BYTEORDER(new_dset)
00239                                                  : THD_get_write_order() ;
00240 
00241    for( nim=ival=0 ; ival < DSET_NVALS(old_dset) ; ival++ ){
00242       brfac_save = DBLK_BRICK_FACTOR(new_dblk,ival) ;
00243       DBLK_BRICK_FACTOR(new_dblk,ival) = 0.0 ;
00244       DBLK_BRICK_FACTOR(old_dblk,ival) = 0.0 ;
00245       dsiz = mri_datum_size( DSET_BRICK_TYPE(new_dset,ival) ) ;
00246       for( kk=0 ; kk < brax->n3 ; kk++ ){
00247          im   = FD_warp_to_mri( kk , ival , brax ) ;
00248          imar = mri_data_pointer(im) ;
00249          if( save_order != native_order ){                   /* 23 Nov 1999 */
00250             switch( im->kind ){
00251                case MRI_short:   mri_swap2(  npix,imar) ; break ;
00252                case MRI_float:
00253                case MRI_int:     mri_swap4(  npix,imar) ; break ;
00254                case MRI_complex: mri_swap4(2*npix,imar) ; break ;
00255             }
00256          }
00257          code = fwrite( imar , dsiz , npix , far ) ;
00258          mri_free(im) ;
00259 
00260          if( verbose && (++nim)%pim == 1 ){ printf("."); fflush(stdout); }
00261       }
00262       DBLK_BRICK_FACTOR(new_dblk,ival) = brfac_save ;
00263       DSET_unload_one(old_dset,ival) ;
00264 
00265       if( verbose ){ printf("!"); fflush(stdout); }
00266    }
00267    COMPRESS_fclose(far) ;
00268    if( verbose ){ printf("\n") ; fflush(stdout); }
00269 
00270    DSET_unload( old_dset ) ;
00271 
00272    /*- do the output header -*/
00273 
00274    DSET_load( new_dset ) ; THD_load_statistics( new_dset ) ;
00275    THD_write_3dim_dataset( NULL,NULL , new_dset , False ) ;
00276    if( verbose ) fprintf(stderr,"++ Wrote new dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
00277 
00278    exit(0) ;
00279 }