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_shear3d.h"
00009 #include <string.h>
00010 #include <stdlib.h>
00011 #include <ctype.h>
00012 
00013 /*--------- Internal Prototype -----------*/
00014 
00015 static void mangle_angle( THD_3dim_dataset *, float, char, float *, int *) ;
00016 
00017 /*----------------------------------------------------------------------------
00018    Routine to take user-input angles (about I,S,R,L,A,P) and return
00019    dataset-oriented angles (about axes 0,1,2).  Adapted from 3drotate.c.
00020 
00021    dset     = dataset defining axes orientations
00022    thI_in   = input angle I (I=1,2,3)
00023    axI_code = axis code for angle I, which is one of
00024                'I' 'S' 'R' 'L' 'A' 'P' 'x' 'y' 'z'
00025 
00026   *thI_out  = output angle I (will be thI_in or -thI_in)
00027   *axI_out  = 0,1,2 indicating rotation about dataset x,y,z axis
00028 
00029   If the input dataset is NULL, then it is treated as in Dicom axes order.
00030 ------------------------------------------------------------------------------*/
00031 
00032 void THD_rotangle_user_to_dset( THD_3dim_dataset *dset ,
00033                                 float th1_in   , char ax1_code ,
00034                                 float th2_in   , char ax2_code ,
00035                                 float th3_in   , char ax3_code ,
00036                                 float *th1_out , int *ax1_out  ,
00037                                 float *th2_out , int *ax2_out  ,
00038                                 float *th3_out , int *ax3_out   )
00039 {
00040 ENTRY("THD_rotangle_user_to_dset") ;
00041 
00042    mangle_angle( dset, th1_in, ax1_code, th1_out, ax1_out ) ;
00043    mangle_angle( dset, th2_in, ax2_code, th2_out, ax2_out ) ;
00044    mangle_angle( dset, th3_in, ax3_code, th3_out, ax3_out ) ;
00045 
00046    if( THD_handedness(dset) < 0 ){
00047       *th1_out = -(*th1_out) ;
00048       *th2_out = -(*th2_out) ;
00049       *th3_out = -(*th3_out) ;
00050    }
00051 
00052    EXRETURN ;
00053 }
00054 
00055 /*----------------------------------------------------------------------------*/
00056 
00057 static void mangle_angle( THD_3dim_dataset *dset ,
00058                           float thin, char ax , float *thout, int *axout )
00059 {
00060    int neg=0 , ax1 ;
00061    float th1=thin ;
00062 
00063    switch( ax ){
00064       default:
00065          if( th1 == 0.0 ){ *thout = 0.0 ; *axout = 0 ; }
00066          else { fprintf(stderr,"*** Illegal ax in mangle_angle\n") ; }
00067       return ;
00068 
00069       case '\0': case 'x': case 'X': ax1 = 0 ; break ;
00070                  case 'y': case 'Y': ax1 = 1 ; break ;
00071                  case 'z': case 'Z': ax1 = 2 ; break ;
00072 
00073       case 'A': case 'P':
00074       case 'R': case 'L':
00075       case 'I': case 'S': ax1 = THD_axcode(dset,ax) ;
00076                           neg = (ax1 < 0) ;
00077                           ax1 = abs(ax1) - 1 ; break ;
00078    }
00079    if( neg ) th1 = -th1 ;
00080 
00081    *thout = th1 ; *axout = ax1 ; return ;
00082 }
00083 
00084 /*---------------------------------------------------------------------------
00085    Input:  ori = orientation code from ISAPLR
00086    Output: +k if ori is along the positive k-axis of the dataset (k=1,2,3)
00087            -k if ori is along the negative k-axis of the dataset
00088 
00089    06 Feb 2001: Promoted to externally visible.
00090    13 Feb 2001: bad dataset -> use Dicom axes order
00091 -----------------------------------------------------------------------------*/
00092 
00093 int THD_axcode( THD_3dim_dataset *dset , char ori )
00094 {
00095    int xxor , yyor , zzor ;
00096 
00097 ENTRY("THD_axcode") ;
00098 
00099    if( ISVALID_DSET(dset) ){
00100       xxor = dset->daxes->xxorient ;
00101       yyor = dset->daxes->yyorient ;
00102       zzor = dset->daxes->zzorient ;
00103    } else {
00104       xxor = ORI_R2L_TYPE ;   /* 13 Feb 2001: Dicom order */
00105       yyor = ORI_A2P_TYPE ;
00106       zzor = ORI_I2S_TYPE ;
00107    }
00108    ori = toupper(ori) ;
00109    if( ori == ORIENT_tinystr[xxor][0] ) RETURN( 1) ;
00110    if( ori == ORIENT_tinystr[xxor][1] ) RETURN(-1) ;
00111    if( ori == ORIENT_tinystr[yyor][0] ) RETURN( 2) ;
00112    if( ori == ORIENT_tinystr[yyor][1] ) RETURN(-2) ;
00113    if( ori == ORIENT_tinystr[zzor][0] ) RETURN( 3) ;
00114    if( ori == ORIENT_tinystr[zzor][1] ) RETURN(-3) ;
00115    RETURN(-99) ;
00116 }
00117 
00118 /*---------------------------------------------------------------------------
00119    Output: +1 if dataset xyz axes are right handed
00120            -1 if dataset xyz axes are left handed
00121    06 Feb 2001: Promoted to externally visible.
00122 -----------------------------------------------------------------------------*/
00123 
00124 int THD_handedness( THD_3dim_dataset *dset )
00125 {
00126    THD_dataxes *dax ;
00127    THD_mat33 q ;
00128    int col ;
00129    float val ;
00130 
00131 ENTRY("THD_handedness") ;
00132 
00133    if( !ISVALID_DSET(dset) ) RETURN(1) ;
00134 
00135    LOAD_ZERO_MAT(q) ;
00136    dax = dset->daxes ;
00137 
00138    col = 0 ;
00139    switch( dax->xxorient ){
00140       case 0: q.mat[0][col] =  1.0 ; break ;
00141       case 1: q.mat[0][col] = -1.0 ; break ;
00142       case 2: q.mat[1][col] = -1.0 ; break ;
00143       case 3: q.mat[1][col] =  1.0 ; break ;
00144       case 4: q.mat[2][col] =  1.0 ; break ;
00145       case 5: q.mat[2][col] = -1.0 ; break ;
00146    }
00147 
00148    col = 1 ;
00149    switch( dax->yyorient ){
00150       case 0: q.mat[0][col] =  1.0 ; break ;
00151       case 1: q.mat[0][col] = -1.0 ; break ;
00152       case 2: q.mat[1][col] = -1.0 ; break ;
00153       case 3: q.mat[1][col] =  1.0 ; break ;
00154       case 4: q.mat[2][col] =  1.0 ; break ;
00155       case 5: q.mat[2][col] = -1.0 ; break ;
00156    }
00157 
00158    col = 2 ;
00159    switch( dax->zzorient ){
00160       case 0: q.mat[0][col] =  1.0 ; break ;
00161       case 1: q.mat[0][col] = -1.0 ; break ;
00162       case 2: q.mat[1][col] = -1.0 ; break ;
00163       case 3: q.mat[1][col] =  1.0 ; break ;
00164       case 4: q.mat[2][col] =  1.0 ; break ;
00165       case 5: q.mat[2][col] = -1.0 ; break ;
00166    }
00167 
00168    val = MAT_DET(q) ;
00169    if( val > 0.0 ) RETURN( 1) ;  /* right handed */
00170    else            RETURN(-1) ;  /* left handed */
00171 }
00172 
00173 /*-------------------------------------------------------------------------
00174   13 Feb 2001: convert a command of the form
00175                  "-rotate 10A 20I 30R -ashift 5I 10A 15L"
00176                to a matrix/vector pair for the given dataset.
00177 ---------------------------------------------------------------------------*/
00178 
00179 THD_dvecmat THD_rotcom_to_matvec( THD_3dim_dataset *dset , char *rotcom )
00180 {
00181    THD_dvecmat out ;
00182    int nn , rpos=0 , nrc ;
00183    char *buf ;
00184 
00185 ENTRY("THD_rotcom_to_matvec") ;
00186 
00187    LOAD_DIAG_DMAT(out.mm,1,1,1) ;  /* identity matrix */
00188    LOAD_DFVEC3(out.vv,0,0,0) ;     /* and zero shift */
00189 
00190    if( rotcom == NULL || rotcom[0] == '\0' ) RETURN(out) ;
00191 
00192    /*-- compute rotation matrix --*/
00193 
00194    nrc = strlen(rotcom) ;
00195    buf = strstr(rotcom,"-rotate") ;
00196    if( buf != NULL && (buf-rotcom)+10 < nrc ){
00197       float th1,th2,th3 , dh1,dh2,dh3 ;
00198       char  cp1,cp2,cp3 ;
00199       int   ax1,ax2,ax3 ;
00200 
00201       nn = sscanf( buf+7, "%f%c %f%c %f%c",
00202                    &th1,&cp1, &th2,&cp2, &th3,&cp3 );
00203       if( nn < 6 ) RETURN(out) ;
00204       if( isspace(cp1) ) cp1 = 'x' ;  /* should not happen */
00205       if( isspace(cp2) ) cp2 = 'y' ;
00206       if( isspace(cp3) ) cp3 = 'z' ;
00207 
00208       THD_rotangle_user_to_dset( dset , th1,cp1 , th2,cp2 , th3,cp3 ,
00209                                  &dh1,&ax1 , &dh2,&ax2 , &dh3,&ax3   ) ;
00210 
00211       out.mm = rot_to_matrix( ax1 , -(PI/180.0)*dh1,
00212                               ax2 , -(PI/180.0)*dh2,
00213                               ax3 , -(PI/180.0)*dh3 ) ;
00214    }
00215 
00216    /*-- compute shift --*/
00217 
00218                      buf = strstr(rotcom,"-ashift") ;
00219    if( buf == NULL ) buf = strstr(rotcom,"-bshift") ;
00220 
00221    if( buf != NULL && (buf-rotcom)+10 < nrc ){
00222       int bs = (buf[1] == 'b') ;  /* save the BS for later */
00223       float dx,dy,dz , qdx=0,qdy=0,qdz=0 ;
00224       char  cdx,cdy,cdz ;
00225       int   adx,ady,adz ;
00226 
00227       nn = sscanf( buf+7, "%f%c %f%c %f%c",
00228                    &dx,&cdx, &dy,&cdy, &dz,&cdz );
00229       if( nn < 6 ) RETURN(out) ;
00230 
00231       adx = THD_axcode(dset,cdx) ;
00232       if( adx > -99 || dx != 0.0 ){
00233          switch( adx ){
00234             case  1: qdx = -dx ; break ;
00235             default:
00236             case -1: qdx =  dx ; break ;
00237             case  2: qdy = -dx ; break ;
00238             case -2: qdy =  dx ; break ;
00239             case  3: qdz = -dx ; break ;
00240             case -3: qdz =  dx ; break ;
00241          }
00242       }
00243 
00244       ady = THD_axcode(dset,cdy) ;
00245       if( ady > -99 || dy != 0.0 ){
00246          switch( ady ){
00247             case  1: qdx = -dy ; break ;
00248             case -1: qdx =  dy ; break ;
00249             case  2: qdy = -dy ; break ;
00250             default:
00251             case -2: qdy =  dy ; break ;
00252             case  3: qdz = -dy ; break ;
00253             case -3: qdz =  dy ; break ;
00254          }
00255       }
00256 
00257       adz = THD_axcode(dset,cdz) ;
00258       if( adz > -99 || dz != 0.0 ){
00259          switch( adz ){
00260             case  1: qdx = -dz ; break ;
00261             case -1: qdx =  dz ; break ;
00262             case  2: qdy = -dz ; break ;
00263             case -2: qdy =  dz ; break ;
00264             case  3: qdz = -dz ; break ;
00265             default:
00266             case -3: qdz =  dz ; break ;
00267          }
00268       }
00269 
00270       LOAD_DFVEC3( out.vv , qdx,qdy,qdz ) ;
00271       if( bs ){
00272          THD_dfvec3 qv = DMATVEC( out.mm , out.vv ) ;
00273          out.vv = qv ;
00274       }
00275    }
00276 
00277    RETURN(out) ;
00278 }
00279 
00280 #if 0
00281 /******************************************
00282    Sample program for above routine
00283 *******************************************/
00284 #include "mrilib.h"
00285 int main( int argc , char *argv[] )
00286 {
00287    THD_3dim_dataset *qset ;
00288    THD_dvecmat vm ;
00289    if( argc < 3 || strcmp(argv[1],"-help")==0 ){
00290      printf("Usage: %s dset \"-rotate a b c -ashift a b c\"\n",argv[0]);exit(0);
00291    }
00292    qset = THD_open_dataset(argv[1]) ;
00293    vm   = THD_rotcom_to_matvec( qset , argv[2] ) ;
00294    DUMP_DVECMAT("r2m output",vm) ; exit(0) ;
00295 }
00296 #endif
00297 
00298 /**************************************************************************/
00299 
00300 #include "thd.h"
00301 
00302 /*---------------------------------------------------------------------
00303   This produces a permutation-like matrix that transforms from
00304   brick axis coordinates to Dicom order coordinates.
00305   [14 Feb 2001 - moved here from 3drotate.c]
00306 -----------------------------------------------------------------------*/
00307 
00308 THD_dmat33 DBLE_mat_to_dicomm( THD_3dim_dataset *dset )
00309 {
00310    THD_dmat33 tod ;
00311 
00312    LOAD_ZERO_DMAT(tod) ;
00313 
00314    switch( dset->daxes->xxorient ){
00315       case ORI_R2L_TYPE: tod.mat[0][0] =  1.0 ; break ;
00316       case ORI_L2R_TYPE: tod.mat[0][0] = -1.0 ; break ;
00317       case ORI_P2A_TYPE: tod.mat[1][0] = -1.0 ; break ;
00318       case ORI_A2P_TYPE: tod.mat[1][0] =  1.0 ; break ;
00319       case ORI_I2S_TYPE: tod.mat[2][0] =  1.0 ; break ;
00320       case ORI_S2I_TYPE: tod.mat[2][0] = -1.0 ; break ;
00321 
00322       default: THD_FATAL_ERROR("illegal xxorient code") ;
00323    }
00324 
00325    switch( dset->daxes->yyorient ){
00326       case ORI_R2L_TYPE: tod.mat[0][1] =  1.0 ; break ;
00327       case ORI_L2R_TYPE: tod.mat[0][1] = -1.0 ; break ;
00328       case ORI_P2A_TYPE: tod.mat[1][1] = -1.0 ; break ;
00329       case ORI_A2P_TYPE: tod.mat[1][1] =  1.0 ; break ;
00330       case ORI_I2S_TYPE: tod.mat[2][1] =  1.0 ; break ;
00331       case ORI_S2I_TYPE: tod.mat[2][1] = -1.0 ; break ;
00332 
00333       default: THD_FATAL_ERROR("illegal yyorient code") ;
00334    }
00335 
00336    switch( dset->daxes->zzorient ){
00337       case ORI_R2L_TYPE: tod.mat[0][2] =  1.0 ; break ;
00338       case ORI_L2R_TYPE: tod.mat[0][2] = -1.0 ; break ;
00339       case ORI_P2A_TYPE: tod.mat[1][2] = -1.0 ; break ;
00340       case ORI_A2P_TYPE: tod.mat[1][2] =  1.0 ; break ;
00341       case ORI_I2S_TYPE: tod.mat[2][2] =  1.0 ; break ;
00342       case ORI_S2I_TYPE: tod.mat[2][2] = -1.0 ; break ;
00343 
00344       default: THD_FATAL_ERROR("illegal zxorient code") ;
00345    }
00346 
00347    return tod ;
00348 }
00349 
00350 /*---------------------------------------------------------------------
00351   This produces a permutation-like matrix that transforms from
00352   brick axis coordinates to Dicom order coordinates.
00353 -----------------------------------------------------------------------*/
00354 
00355 THD_mat33 SNGL_mat_to_dicomm( THD_3dim_dataset *dset )
00356 {
00357    THD_mat33 tod ;
00358 
00359    LOAD_ZERO_MAT(tod) ;
00360 
00361    switch( dset->daxes->xxorient ){
00362       case ORI_R2L_TYPE: tod.mat[0][0] =  1.0 ; break ;
00363       case ORI_L2R_TYPE: tod.mat[0][0] = -1.0 ; break ;
00364       case ORI_P2A_TYPE: tod.mat[1][0] = -1.0 ; break ;
00365       case ORI_A2P_TYPE: tod.mat[1][0] =  1.0 ; break ;
00366       case ORI_I2S_TYPE: tod.mat[2][0] =  1.0 ; break ;
00367       case ORI_S2I_TYPE: tod.mat[2][0] = -1.0 ; break ;
00368 
00369       default: THD_FATAL_ERROR("illegal xxorient code") ;
00370    }
00371 
00372    switch( dset->daxes->yyorient ){
00373       case ORI_R2L_TYPE: tod.mat[0][1] =  1.0 ; break ;
00374       case ORI_L2R_TYPE: tod.mat[0][1] = -1.0 ; break ;
00375       case ORI_P2A_TYPE: tod.mat[1][1] = -1.0 ; break ;
00376       case ORI_A2P_TYPE: tod.mat[1][1] =  1.0 ; break ;
00377       case ORI_I2S_TYPE: tod.mat[2][1] =  1.0 ; break ;
00378       case ORI_S2I_TYPE: tod.mat[2][1] = -1.0 ; break ;
00379 
00380       default: THD_FATAL_ERROR("illegal yyorient code") ;
00381    }
00382 
00383    switch( dset->daxes->zzorient ){
00384       case ORI_R2L_TYPE: tod.mat[0][2] =  1.0 ; break ;
00385       case ORI_L2R_TYPE: tod.mat[0][2] = -1.0 ; break ;
00386       case ORI_P2A_TYPE: tod.mat[1][2] = -1.0 ; break ;
00387       case ORI_A2P_TYPE: tod.mat[1][2] =  1.0 ; break ;
00388       case ORI_I2S_TYPE: tod.mat[2][2] =  1.0 ; break ;
00389       case ORI_S2I_TYPE: tod.mat[2][2] = -1.0 ; break ;
00390 
00391       default: THD_FATAL_ERROR("illegal zxorient code") ;
00392    }
00393 
00394    return tod ;
00395 }