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    3D coordinate conversion routines;
00012      tags for coordinate systems:
00013        fdind  = FD_brick voxel indices (ints)
00014        fdfind = FD_brick voxel indices (floats - added 30 Aug 2001)
00015        3dind  = THD_3dim_dataset voxel indices (ints)
00016        3dfind = THD_3dim_dataset floating point voxel indices
00017                   (used for subvoxel resolution)
00018        3dmm   = THD_3dim_dataset millimetric coordinates (floats)
00019        dicomm = DICOM 3.0 millimetric coordinates (floats)
00020 
00021      The 3dmm coordinate measurements are oriented the same way
00022      as the dicomm coordinates (which means that the ??del values
00023      can be negative if the voxel axes are reflected), but the 3dmm
00024      coordinates are in general a permutation of the DICOM coordinates.
00025 
00026      These routines all take as input and produce as output
00027      THD_?vec3 (i=int,f=float) structs.
00028 ======================================================================*/
00029 
00030 THD_fvec3 THD_3dfind_to_3dmm( THD_3dim_dataset * dset ,
00031                               THD_fvec3 iv )
00032 {
00033    THD_dataxes * daxes ;
00034    THD_fvec3     fv ;
00035 
00036    daxes = CURRENT_DAXES(dset) ;
00037 
00038    fv.xyz[0] = daxes->xxorg + iv.xyz[0] * daxes->xxdel ;
00039    fv.xyz[1] = daxes->yyorg + iv.xyz[1] * daxes->yydel ;
00040    fv.xyz[2] = daxes->zzorg + iv.xyz[2] * daxes->zzdel ;
00041    return fv ;
00042 }
00043 
00044 /*------------------------------------------------------------------*/
00045 
00046 THD_fvec3 THD_3dind_to_3dmm( THD_3dim_dataset * dset ,
00047                              THD_ivec3 iv )
00048 {
00049    THD_dataxes * daxes ;
00050    THD_fvec3     fv ;
00051 
00052    daxes = CURRENT_DAXES(dset) ;
00053 
00054    fv.xyz[0] = daxes->xxorg + iv.ijk[0] * daxes->xxdel ;
00055    fv.xyz[1] = daxes->yyorg + iv.ijk[1] * daxes->yydel ;
00056    fv.xyz[2] = daxes->zzorg + iv.ijk[2] * daxes->zzdel ;
00057    return fv ;
00058 }
00059 
00060 /*--------------------------------------------------------------------*/
00061 /* this version is without using wod dataxes       7 Mar 2005 [rickr] */
00062 THD_fvec3 THD_3dind_to_3dmm_no_wod( THD_3dim_dataset * dset ,
00063                              THD_ivec3 iv )
00064 {
00065    THD_dataxes * daxes ;
00066    THD_fvec3     fv ;
00067 
00068    daxes = dset->daxes ;
00069 
00070    fv.xyz[0] = daxes->xxorg + iv.ijk[0] * daxes->xxdel ;
00071    fv.xyz[1] = daxes->yyorg + iv.ijk[1] * daxes->yydel ;
00072    fv.xyz[2] = daxes->zzorg + iv.ijk[2] * daxes->zzdel ;
00073    return fv ;
00074 }
00075 
00076 /*--------------------------------------------------------------------*/
00077 
00078 THD_fvec3 THD_3dmm_to_3dfind( THD_3dim_dataset * dset ,
00079                               THD_fvec3 fv )
00080 {
00081    THD_dataxes * daxes ;
00082    THD_fvec3     iv ;
00083 
00084    daxes = CURRENT_DAXES(dset) ;
00085 
00086    iv.xyz[0] = (fv.xyz[0] - daxes->xxorg) / daxes->xxdel ;
00087    iv.xyz[1] = (fv.xyz[1] - daxes->yyorg) / daxes->yydel ;
00088    iv.xyz[2] = (fv.xyz[2] - daxes->zzorg) / daxes->zzdel ;
00089 
00090         if( iv.xyz[0] < 0            ) iv.xyz[0] = 0 ;
00091    else if( iv.xyz[0] > daxes->nxx-1 ) iv.xyz[0] = daxes->nxx-1 ;
00092 
00093         if( iv.xyz[1] <  0           ) iv.xyz[1] = 0 ;
00094    else if( iv.xyz[1] > daxes->nyy-1 ) iv.xyz[1] = daxes->nyy-1 ;
00095 
00096         if( iv.xyz[2] < 0            ) iv.xyz[2] = 0 ;
00097    else if( iv.xyz[2] > daxes->nzz-1 ) iv.xyz[2] = daxes->nzz-1 ;
00098 
00099    return iv ;
00100 }
00101 
00102 /*--------------------------------------------------------------------*/
00103 
00104 THD_ivec3 THD_3dmm_to_3dind( THD_3dim_dataset * dset ,
00105                              THD_fvec3 fv )
00106 {
00107    THD_dataxes * daxes ;
00108    THD_ivec3     iv ;
00109 
00110    daxes = CURRENT_DAXES(dset) ;
00111 
00112    iv.ijk[0] = (fv.xyz[0] - daxes->xxorg) / daxes->xxdel + 0.499 ;
00113    iv.ijk[1] = (fv.xyz[1] - daxes->yyorg) / daxes->yydel + 0.499 ;
00114    iv.ijk[2] = (fv.xyz[2] - daxes->zzorg) / daxes->zzdel + 0.499 ;
00115 
00116         if( iv.ijk[0] < 0            ) iv.ijk[0] = 0 ;
00117    else if( iv.ijk[0] > daxes->nxx-1 ) iv.ijk[0] = daxes->nxx-1 ;
00118 
00119         if( iv.ijk[1] < 0            ) iv.ijk[1] = 0 ;
00120    else if( iv.ijk[1] > daxes->nyy-1 ) iv.ijk[1] = daxes->nyy-1 ;
00121 
00122         if( iv.ijk[2] < 0            ) iv.ijk[2] = 0 ;
00123    else if( iv.ijk[2] > daxes->nzz-1 ) iv.ijk[2] = daxes->nzz-1 ;
00124 
00125    return iv ;
00126 }
00127 
00128 /*--------------------------------------------------------------------*/
00129 
00130 /* this version is without using wod dataxes     28 Sep 2004 [rickr] */
00131 THD_ivec3 THD_3dmm_to_3dind_no_wod( THD_3dim_dataset * dset ,
00132                                     THD_fvec3 fv )
00133 {
00134    THD_dataxes * daxes ;
00135    THD_ivec3     iv ;
00136 
00137    daxes = dset->daxes ;
00138 
00139    iv.ijk[0] = (fv.xyz[0] - daxes->xxorg) / daxes->xxdel + 0.499 ;
00140    iv.ijk[1] = (fv.xyz[1] - daxes->yyorg) / daxes->yydel + 0.499 ;
00141    iv.ijk[2] = (fv.xyz[2] - daxes->zzorg) / daxes->zzdel + 0.499 ;
00142 
00143         if( iv.ijk[0] < 0            ) iv.ijk[0] = 0 ;
00144    else if( iv.ijk[0] > daxes->nxx-1 ) iv.ijk[0] = daxes->nxx-1 ;
00145 
00146         if( iv.ijk[1] < 0            ) iv.ijk[1] = 0 ;
00147    else if( iv.ijk[1] > daxes->nyy-1 ) iv.ijk[1] = daxes->nyy-1 ;
00148 
00149         if( iv.ijk[2] < 0            ) iv.ijk[2] = 0 ;
00150    else if( iv.ijk[2] > daxes->nzz-1 ) iv.ijk[2] = daxes->nzz-1 ;
00151 
00152    return iv ;
00153 }
00154 
00155 /*---------------------------------------------------------------------
00156    convert from input image oriented x,y,z to Dicom x,y,z
00157      (x axis = R->L , y axis = A->P , z axis = I->S)
00158 
00159    N.B.: image distances are oriented the same as Dicom,
00160          just in a permuted order.
00161 -----------------------------------------------------------------------*/
00162 
00163 THD_fvec3 THD_3dmm_to_dicomm( THD_3dim_dataset * dset ,
00164                               THD_fvec3 imv )
00165 {
00166    THD_fvec3 dicv ;
00167    float xim,yim,zim , xdic,ydic,zdic ;
00168 
00169    xim = imv.xyz[0] ; yim = imv.xyz[1] ; zim = imv.xyz[2] ;
00170 
00171    switch( dset->daxes->xxorient ){
00172       case ORI_R2L_TYPE:
00173       case ORI_L2R_TYPE: xdic = xim ; break ;
00174       case ORI_P2A_TYPE:
00175       case ORI_A2P_TYPE: ydic = xim ; break ;
00176       case ORI_I2S_TYPE:
00177       case ORI_S2I_TYPE: zdic = xim ; break ;
00178 
00179       default: THD_FATAL_ERROR("illegal xxorient code") ;
00180    }
00181 
00182    switch( dset->daxes->yyorient ){
00183       case ORI_R2L_TYPE:
00184       case ORI_L2R_TYPE: xdic = yim ; break ;
00185       case ORI_P2A_TYPE:
00186       case ORI_A2P_TYPE: ydic = yim ; break ;
00187       case ORI_I2S_TYPE:
00188       case ORI_S2I_TYPE: zdic = yim ; break ;
00189 
00190       default: THD_FATAL_ERROR("illegal yyorient code") ;
00191    }
00192 
00193    switch( dset->daxes->zzorient ){
00194       case ORI_R2L_TYPE:
00195       case ORI_L2R_TYPE: xdic = zim ; break ;
00196       case ORI_P2A_TYPE:
00197       case ORI_A2P_TYPE: ydic = zim ; break ;
00198       case ORI_I2S_TYPE:
00199       case ORI_S2I_TYPE: zdic = zim ; break ;
00200 
00201       default: THD_FATAL_ERROR("illegal zzorient code") ;
00202    }
00203 
00204    dicv.xyz[0] = xdic ; dicv.xyz[1] = ydic ; dicv.xyz[2] = zdic ;
00205    return dicv ;
00206 }
00207 
00208 /*---------------------------------------------------------------------
00209    convert to input image oriented x,y,z from Dicom x,y,z
00210 -----------------------------------------------------------------------*/
00211 
00212 THD_fvec3 THD_dicomm_to_3dmm( THD_3dim_dataset * dset ,
00213                               THD_fvec3 dicv )
00214 {
00215    THD_fvec3 imv ;
00216    float xim,yim,zim , xdic,ydic,zdic ;
00217 
00218    xdic = dicv.xyz[0] ; ydic = dicv.xyz[1] ; zdic = dicv.xyz[2] ;
00219 
00220    switch( dset->daxes->xxorient ){
00221       case ORI_R2L_TYPE:
00222       case ORI_L2R_TYPE: xim = xdic ; break ;
00223       case ORI_P2A_TYPE:
00224       case ORI_A2P_TYPE: xim = ydic ; break ;
00225       case ORI_I2S_TYPE:
00226       case ORI_S2I_TYPE: xim = zdic ; break ;
00227 
00228       default: THD_FATAL_ERROR("illegal xxorient code") ;
00229    }
00230 
00231    switch( dset->daxes->yyorient ){
00232       case ORI_R2L_TYPE:
00233       case ORI_L2R_TYPE: yim = xdic ; break ;
00234       case ORI_P2A_TYPE:
00235       case ORI_A2P_TYPE: yim = ydic ; break ;
00236       case ORI_I2S_TYPE:
00237       case ORI_S2I_TYPE: yim = zdic ; break ;
00238 
00239       default: THD_FATAL_ERROR("illegal yyorient code") ;
00240    }
00241 
00242    switch( dset->daxes->zzorient ){
00243       case ORI_R2L_TYPE:
00244       case ORI_L2R_TYPE: zim = xdic ; break ;
00245       case ORI_P2A_TYPE:
00246       case ORI_A2P_TYPE: zim = ydic ; break ;
00247       case ORI_I2S_TYPE:
00248       case ORI_S2I_TYPE: zim = zdic ; break ;
00249 
00250       default: THD_FATAL_ERROR("illegal zzorient code") ;
00251    }
00252 
00253    imv.xyz[0] = xim ; imv.xyz[1] = yim ; imv.xyz[2] = zim ;
00254    return imv ;
00255 }
00256 
00257 /*---------------------------------------------------------------------*/
00258 
00259 THD_ivec3 THD_fdind_to_3dind( FD_brick * br , THD_ivec3 ib )
00260 {
00261    THD_ivec3 id ;
00262    int qq , ax ;
00263 
00264    for( qq=0 ; qq < 3 ; qq++ ){
00265       ax = abs( br->a123.ijk[qq] ) - 1 ;   /* 0,1,2, for x,y,z */
00266 
00267       if( br->a123.ijk[qq] > 0 ) id.ijk[ax] = ib.ijk[qq] ;
00268       else                       id.ijk[ax] = br->sxyz.ijk[ax] - ib.ijk[qq];
00269    }
00270 
00271    return id ;
00272 }
00273 
00274 THD_ivec3 THD_3dind_to_fdind( FD_brick * br , THD_ivec3 id )
00275 {
00276    THD_ivec3 ib ;
00277    int qq , ax ;
00278 
00279    for( qq=0 ; qq < 3 ; qq++ ){
00280       ax = abs( br->a123.ijk[qq] ) - 1 ;
00281 
00282       if( br->a123.ijk[qq] > 0 ) ib.ijk[qq] = id.ijk[ax] ;
00283       else                       ib.ijk[qq] = br->sxyz.ijk[ax] - id.ijk[ax];
00284    }
00285 
00286    return ib ;
00287 }
00288 
00289 /*---------------------------------------------------------------------*/
00290 
00291 THD_fvec3 THD_fdfind_to_3dfind( FD_brick * br , THD_fvec3 ib ) /* 30 Aug 2001 */
00292 {
00293    THD_fvec3 id ;
00294    int qq , ax ;
00295 
00296    for( qq=0 ; qq < 3 ; qq++ ){
00297       ax = abs( br->a123.ijk[qq] ) - 1 ;   /* 0,1,2, for x,y,z */
00298 
00299       if( br->a123.ijk[qq] > 0 ) id.xyz[ax] = ib.xyz[qq] ;
00300       else                       id.xyz[ax] = br->sxyz.ijk[ax] - ib.xyz[qq];
00301    }
00302 
00303    return id ;
00304 }
00305 
00306 THD_fvec3 THD_3dfind_to_fdfind( FD_brick * br , THD_fvec3 id ) /* 30 Aug 2001 */
00307 {
00308    THD_fvec3 ib ;
00309    int qq , ax ;
00310 
00311    for( qq=0 ; qq < 3 ; qq++ ){
00312       ax = abs( br->a123.ijk[qq] ) - 1 ;
00313 
00314       if( br->a123.ijk[qq] > 0 ) ib.xyz[qq] = id.xyz[ax] ;
00315       else                       ib.xyz[qq] = br->sxyz.ijk[ax] - id.xyz[ax];
00316    }
00317 
00318    return ib ;
00319 }
00320 
00321 /*-------------------------------------------------------------------*/
00322 
00323 void THD_coorder_fill( char * in_orcode , THD_coorder * cord )
00324 {
00325    char acod , orcode[4] ;
00326    int xx,yy,zz , ss1,ss2,ss3 , ii,ll ;
00327 
00328    if( cord == NULL ) return ;
00329 
00330    /* default values */
00331 
00332    cord->xxsign = cord->yysign = cord->zzsign = 1 ;
00333    cord->first  = 0 ;
00334    cord->second = 1 ;
00335    cord->third  = 2 ;
00336    cord->xxor   = ORI_R2L_TYPE ;
00337    cord->yyor   = ORI_A2P_TYPE ;
00338    cord->zzor   = ORI_I2S_TYPE ;
00339    strcpy(cord->orcode,"RAI") ;
00340 
00341    /* check string for OKness */
00342 
00343    if( in_orcode == NULL ) return ;
00344    strncpy(orcode,in_orcode,3) ; orcode[3] = '\0' ;
00345    ll = strlen(orcode) ; if( ll != 3 ) return ;
00346    for( ii=0 ; ii < 3 ; ii++ ) orcode[ii] = toupper(orcode[ii]) ;
00347    if( strncmp(orcode,"FLI",3) == 0 ) strcpy(orcode,"LPI") ;
00348 
00349    /* extract direction codes */
00350 
00351    acod = orcode[0] ; xx = ORCODE(acod) ;
00352    acod = orcode[1] ; yy = ORCODE(acod) ;
00353    acod = orcode[2] ; zz = ORCODE(acod) ;
00354 
00355    /* check direction codes for OKness */
00356 
00357    if( xx<0 || yy<0 || zz<0 || ! OR3OK(xx,yy,zz) ) return ;
00358 
00359    /* all is OK.  get signs of orientations */
00360 
00361    ss1 = (ORIENT_sign[xx] == '-') ? -1 : 1 ;
00362    ss2 = (ORIENT_sign[yy] == '-') ? -1 : 1 ;
00363    ss3 = (ORIENT_sign[zz] == '-') ? -1 : 1 ;
00364 
00365    /* whose on first? */
00366 
00367    cord->first  = xx / 2 ;
00368    cord->second = yy / 2 ;
00369    cord->third  = zz / 2 ;
00370 
00371    cord->xxsign = (cord->first ==0) ? ss1
00372                  :(cord->second==0) ? ss2 : ss3 ;
00373 
00374    cord->yysign = (cord->first ==1) ? ss1
00375                  :(cord->second==1) ? ss2 : ss3 ;
00376 
00377    cord->zzsign = (cord->first ==2) ? ss1
00378                  :(cord->second==2) ? ss2 : ss3 ;
00379 
00380    cord->xxor = xx ;
00381    cord->yyor = yy ;
00382    cord->zzor = zz ;
00383 
00384    strcpy(cord->orcode,orcode) ;
00385    return ;
00386 }
00387 
00388 /*---------------------------------------------------------------------
00389    convert to output order x,y,z from Dicom x,y,z
00390 -----------------------------------------------------------------------*/
00391 
00392 void THD_dicom_to_coorder( THD_coorder * cord ,
00393                            float * xx , float * yy , float * zz )
00394 {
00395    float xval , yval , zval ;
00396 
00397    if( cord == NULL ) return ;
00398 
00399    /* changes signs first */
00400 
00401    xval = cord->xxsign * (*xx) ;
00402    yval = cord->yysign * (*yy) ;
00403    zval = cord->zzsign * (*zz) ;
00404 
00405    /* scramble order */
00406 
00407    *xx = (cord->first == 0) ? xval
00408         :(cord->first == 1) ? yval : zval ;
00409 
00410    *yy = (cord->second == 0) ? xval
00411         :(cord->second == 1) ? yval : zval ;
00412 
00413    *zz = (cord->third == 0) ? xval
00414         :(cord->third == 1) ? yval : zval ;
00415 
00416    return ;
00417 }
00418 
00419 /*-------------------------------------------------------------------
00420    convert to Dicom x,y,z from output order x,y,z
00421 ---------------------------------------------------------------------*/
00422 
00423 void THD_coorder_to_dicom( THD_coorder * cord ,
00424                            float * xx , float * yy , float * zz )
00425 {
00426    float xval , yval , zval ;
00427 
00428    if( cord == NULL ) return ;
00429 
00430    /* unscramble order */
00431 
00432    xval = (cord->first  == 0) ? (*xx)
00433          :(cord->second == 0) ? (*yy) : (*zz) ;
00434 
00435    yval = (cord->first  == 1) ? (*xx)
00436          :(cord->second == 1) ? (*yy) : (*zz) ;
00437 
00438    zval = (cord->first  == 2) ? (*xx)
00439          :(cord->second == 2) ? (*yy) : (*zz) ;
00440 
00441    /* change signs */
00442 
00443    *xx = cord->xxsign * xval ;
00444    *yy = cord->yysign * yval ;
00445    *zz = cord->zzsign * zval ;
00446 
00447    return ;
00448 }