00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 #include "thd.h"
00009 
00010 /*-----------------------------------------------------------------------
00011   Create FD_bricks for viewing purposes.
00012 
00013   FD_bricks are a relic of the very earliest code that evolved into
00014   AFNI.  By the time I invented THD_3dim_datasets, I had enough code
00015   in place to keep this data structure alive.  Basically, it exists
00016   only to provide a structure that enables fast copying of data out
00017   of a 3D array into 2D arrays.
00018 -------------------------------------------------------------------------*/
00019 
00020 FD_brick ** THD_setup_bricks( THD_3dim_dataset * dset )
00021 {
00022    int r2l=0 , a2p=0 , i2s=0 ;
00023    THD_dataxes * daxes ;
00024    FD_brick ** br ;
00025 
00026    if( ! ISVALID_3DIM_DATASET(dset) ) return NULL ;
00027 
00028    daxes = CURRENT_DAXES(dset) ;
00029    if( ! ISVALID_DATAXES(daxes) ) return NULL ;
00030 
00031    /*----- create FD_bricks for viewing purposes -----*/
00032 
00033    switch( daxes->xxorient ){
00034       case ORI_R2L_TYPE: r2l =  1 ; break ;
00035       case ORI_L2R_TYPE: r2l = -1 ; break ;
00036       case ORI_P2A_TYPE: a2p = -1 ; break ;
00037       case ORI_A2P_TYPE: a2p =  1 ; break ;
00038       case ORI_I2S_TYPE: i2s =  1 ; break ;
00039       case ORI_S2I_TYPE: i2s = -1 ; break ;
00040    }
00041 
00042    switch( daxes->yyorient ){
00043       case ORI_R2L_TYPE: r2l =  2 ; break ;
00044       case ORI_L2R_TYPE: r2l = -2 ; break ;
00045       case ORI_P2A_TYPE: a2p = -2 ; break ;
00046       case ORI_A2P_TYPE: a2p =  2 ; break ;
00047       case ORI_I2S_TYPE: i2s =  2 ; break ;
00048       case ORI_S2I_TYPE: i2s = -2 ; break ;
00049    }
00050 
00051    switch( daxes->zzorient ){
00052       case ORI_R2L_TYPE: r2l =  3 ; break ;
00053       case ORI_L2R_TYPE: r2l = -3 ; break ;
00054       case ORI_P2A_TYPE: a2p = -3 ; break ;
00055       case ORI_A2P_TYPE: a2p =  3 ; break ;
00056       case ORI_I2S_TYPE: i2s =  3 ; break ;
00057       case ORI_S2I_TYPE: i2s = -3 ; break ;
00058    }
00059 
00060    if( r2l==0 || a2p==0 || i2s==0 ){
00061       char buf[256] ;
00062       sprintf(buf,"Illegal orientation codes: %d %d %d",
00063                   daxes->xxorient,daxes->yyorient,daxes->zzorient ) ;
00064       THD_FATAL_ERROR(buf) ;
00065    }
00066 
00067    /* now we can set up views: axial, sagittal, coronal;
00068       the top option is the way I think it ought to be, so that
00069       in the axial and coronal views, left is left and right is right;
00070       the bottom options are the radiologists conventions, so sue me! */
00071 
00072    br = (FD_brick **) XtMalloc( sizeof(FD_brick *) * 3 ) ;
00073 
00074 #undef LEFT_IS_LEFT
00075 #ifdef LEFT_IS_LEFT
00076       br[0] = THD_3dim_dataset_to_brick(dset,-r2l, a2p,-i2s); /* axi */
00077       br[1] = THD_3dim_dataset_to_brick(dset, a2p,-i2s,-r2l); /* sag */
00078       br[2] = THD_3dim_dataset_to_brick(dset,-r2l,-i2s,-a2p); /* cor */
00079 #else
00080       br[0] = THD_3dim_dataset_to_brick(dset, r2l, a2p, i2s); /* axi */
00081       br[1] = THD_3dim_dataset_to_brick(dset, a2p,-i2s,-r2l); /* sag */
00082       br[2] = THD_3dim_dataset_to_brick(dset, r2l,-i2s, a2p); /* cor */
00083 #endif
00084 
00085    strcpy( br[0]->namecode , "Axial" ) ;
00086    strcpy( br[1]->namecode , "Sagittal" ) ;
00087    strcpy( br[2]->namecode , "Coronal" ) ;
00088 
00089    return br ;
00090 }
00091 
00092 /*----------------------------------------------------------------------
00093   07 Dec 2001 - orient an FD brick in any legal way
00094                 (e.g., orients = "RAI" for standard axial)
00095 ------------------------------------------------------------------------*/
00096 
00097 FD_brick * THD_oriented_brick( THD_3dim_dataset *dset , char *orients )
00098 {
00099    int r2l=0 , a2p=0 , i2s=0 , xx,yy,zz , pp=0,qq=0,rr=0 ;
00100    THD_dataxes *daxes ;
00101    FD_brick *br ;
00102 
00103 ENTRY("THD_oriented_brick") ;
00104 
00105    if( !ISVALID_DSET(dset) ||
00106        orients == NULL     ||
00107        strlen(orients) < 3   ) RETURN(NULL) ;
00108 
00109    daxes = CURRENT_DAXES(dset) ;
00110    if( !ISVALID_DATAXES(daxes) ) RETURN(NULL) ;
00111 
00112    xx = ORCODE( orients[0] ) ;
00113    yy = ORCODE( orients[1] ) ;
00114    zz = ORCODE( orients[2] ) ;
00115    if( !OR3OK(xx,yy,zz) ) RETURN(NULL) ;
00116 
00117    /*----- create FD_bricks for viewing purposes -----*/
00118 
00119    switch( daxes->xxorient ){
00120       case ORI_R2L_TYPE: r2l =  1 ; break ;
00121       case ORI_L2R_TYPE: r2l = -1 ; break ;
00122       case ORI_P2A_TYPE: a2p = -1 ; break ;
00123       case ORI_A2P_TYPE: a2p =  1 ; break ;
00124       case ORI_I2S_TYPE: i2s =  1 ; break ;
00125       case ORI_S2I_TYPE: i2s = -1 ; break ;
00126    }
00127 
00128    switch( daxes->yyorient ){
00129       case ORI_R2L_TYPE: r2l =  2 ; break ;
00130       case ORI_L2R_TYPE: r2l = -2 ; break ;
00131       case ORI_P2A_TYPE: a2p = -2 ; break ;
00132       case ORI_A2P_TYPE: a2p =  2 ; break ;
00133       case ORI_I2S_TYPE: i2s =  2 ; break ;
00134       case ORI_S2I_TYPE: i2s = -2 ; break ;
00135    }
00136 
00137    switch( daxes->zzorient ){
00138       case ORI_R2L_TYPE: r2l =  3 ; break ;
00139       case ORI_L2R_TYPE: r2l = -3 ; break ;
00140       case ORI_P2A_TYPE: a2p = -3 ; break ;
00141       case ORI_A2P_TYPE: a2p =  3 ; break ;
00142       case ORI_I2S_TYPE: i2s =  3 ; break ;
00143       case ORI_S2I_TYPE: i2s = -3 ; break ;
00144    }
00145 
00146    if( r2l==0 || a2p==0 || i2s==0 ) RETURN(NULL) ;
00147 
00148    switch( xx ){
00149       case ORI_R2L_TYPE: pp =  r2l ; break ;
00150       case ORI_L2R_TYPE: pp = -r2l ; break ;
00151       case ORI_P2A_TYPE: pp = -a2p ; break ;
00152       case ORI_A2P_TYPE: pp =  a2p ; break ;
00153       case ORI_I2S_TYPE: pp =  i2s ; break ;
00154       case ORI_S2I_TYPE: pp = -i2s ; break ;
00155    }
00156 
00157    switch( yy ){
00158       case ORI_R2L_TYPE: qq =  r2l ; break ;
00159       case ORI_L2R_TYPE: qq = -r2l ; break ;
00160       case ORI_P2A_TYPE: qq = -a2p ; break ;
00161       case ORI_A2P_TYPE: qq =  a2p ; break ;
00162       case ORI_I2S_TYPE: qq =  i2s ; break ;
00163       case ORI_S2I_TYPE: qq = -i2s ; break ;
00164    }
00165 
00166    switch( zz ){
00167       case ORI_R2L_TYPE: rr =  r2l ; break ;
00168       case ORI_L2R_TYPE: rr = -r2l ; break ;
00169       case ORI_P2A_TYPE: rr = -a2p ; break ;
00170       case ORI_A2P_TYPE: rr =  a2p ; break ;
00171       case ORI_I2S_TYPE: rr =  i2s ; break ;
00172       case ORI_S2I_TYPE: rr = -i2s ; break ;
00173    }
00174 
00175    if( pp==0 || qq==0 || rr==0 ) RETURN(NULL) ;
00176 
00177    br = THD_3dim_dataset_to_brick(dset,pp,qq,rr) ;
00178    RETURN(br) ;
00179 }
00180 
00181 /*======================================================================
00182     routines adapted from original afni code for display of a 3D dataset
00183 ========================================================================*/
00184 
00185 /*----------------------------------------------------------------------
00186   prepare a "display brick" for a 3D dataset;
00187     ax_n = 1,2,3 for displaying the x,y,z axis along the n-th display
00188               direction (n=1 --> across, n=2 --> down, n=3 --> slices )
00189            negative value means reverse
00190 
00191    a return of NULL means something bad happened
00192 ------------------------------------------------------------------------*/
00193 
00194 FD_brick * THD_3dim_dataset_to_brick( THD_3dim_dataset * dset ,
00195                                       int ax_1, int ax_2, int ax_3 )
00196 {
00197    FD_brick      * br ;    /* will be output */
00198    THD_dataxes   * daxes ; /* connection to actual axes */
00199 
00200    int   xyz_dim[4] , xyz_stp[4] , xyz_dir[4] ;
00201    float xyz_del[4] , xyz_org[4] ;
00202 
00203    int x_dir,y_dir,z_dir , sx,sy,sz , aax_1,aax_2,aax_3 , nx,ny,nz ;
00204 
00205    /*-- sanity check --*/
00206 
00207    if( ! ISVALID_3DIM_DATASET(dset) ) return NULL ;
00208 
00209    daxes = CURRENT_DAXES(dset) ;
00210 
00211    aax_1 = abs(ax_1) ;  /* which axes, regardless of + or - */
00212    aax_2 = abs(ax_2) ;
00213    aax_3 = abs(ax_3) ;
00214 
00215    if( aax_1 < 1 || aax_1 > 3 ||   /* range checks */
00216        aax_2 < 1 || aax_2 > 3 ||
00217        aax_3 < 1 || aax_3 > 3   ) return NULL ;
00218 
00219    xyz_dir[1] = xyz_dir[2] = xyz_dir[3] = 0 ;
00220 
00221    xyz_dir[aax_1] = ax_1 ;  /* assign to original directions */
00222    xyz_dir[aax_2] = ax_2 ;
00223    xyz_dir[aax_3] = ax_3 ;
00224 
00225    x_dir = xyz_dir[1] ;  /* if any |ax_n| is duplicated */
00226    y_dir = xyz_dir[2] ;  /* then one of these will end  */
00227    z_dir = xyz_dir[3] ;  /* up as zero --> bad inputs!  */
00228 
00229    if( x_dir == 0 || y_dir == 0 || z_dir == 0 ) return NULL ;
00230 
00231    /*-- the inputs are good, so create a brick: --*/
00232 
00233    br             = myXtNew(FD_brick) ;  /* new brick */
00234    br->dset       = dset ;               /* dataset */
00235    br->resam_code = RESAM_NN_TYPE ;      /* crudest type */
00236    br->parent     = NULL ;
00237 
00238    br->thr_resam_code = RESAM_NN_TYPE ;  /* 09 Dec 1997 */
00239 
00240    /*-- at this point, x_dir is +1 or -1, y_dir is +2 or -2, etc. --*/
00241 
00242    nx = daxes->nxx ; ny = daxes->nyy ; nz = daxes->nzz ;
00243 
00244    sx = (x_dir > 0) ? (0) : (nx-1) ;  /* starting voxel indices */
00245    sy = (y_dir > 0) ? (0) : (ny-1) ;  /* for each original dimension */
00246    sz = (z_dir > 0) ? (0) : (nz-1) ;
00247 
00248    br->start = sx + sy*nx + sz*nx*ny ; /* overall starting voxel index */
00249 
00250    /*-- assign original dimensions to arrays,
00251         then pick out the permuted dimensions --*/
00252 
00253    xyz_dim[1] = nx ;  /* dimensions */
00254    xyz_dim[2] = ny ;
00255    xyz_dim[3] = nz ;
00256 
00257    LOAD_IVEC3( br->nxyz , nx,ny,nz ) ;       /* save stuff in br */
00258    LOAD_IVEC3( br->sxyz , sx,sy,sz ) ;
00259    LOAD_IVEC3( br->a123 , ax_1,ax_2,ax_3 ) ;
00260 
00261    xyz_stp[1] = 1 ;            /* index step sizes */
00262    xyz_stp[2] = nx ;
00263    xyz_stp[3] = nx * ny ;
00264 
00265    xyz_del[1] = daxes->xxdel ; /* voxel physical step sizes (mm) */
00266    xyz_del[2] = daxes->yydel ;
00267    xyz_del[3] = daxes->zzdel ;
00268 
00269    xyz_org[1] = daxes->xxorg ; /* voxel origins (mm) */
00270    xyz_org[2] = daxes->yyorg ;
00271    xyz_org[3] = daxes->zzorg ;
00272 
00273    br->n1 = xyz_dim[aax_1] ;   /* permute dimensions, etc. */
00274    br->n2 = xyz_dim[aax_2] ;
00275    br->n3 = xyz_dim[aax_3] ;
00276 
00277    br->d1 = (ax_1 > 0) ? (xyz_stp[aax_1]) : (-xyz_stp[aax_1]) ;
00278    br->d2 = (ax_2 > 0) ? (xyz_stp[aax_2]) : (-xyz_stp[aax_2]) ;
00279    br->d3 = (ax_3 > 0) ? (xyz_stp[aax_3]) : (-xyz_stp[aax_3]) ;
00280 
00281    br->e1 = br->n1 * br->d1 ;  /* last indices for readout */
00282    br->e2 = br->n2 * br->d2 ;
00283 
00284    br->del1 = fabs(xyz_del[aax_1]) ;  /* dimensions */
00285    br->del2 = fabs(xyz_del[aax_2]) ;
00286    br->del3 = fabs(xyz_del[aax_3]) ;
00287 
00288    br->namecode[0] = '\0' ;
00289 
00290    return br ;
00291 }