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   July 1997: moved warp manipulations routines from afni_warp.c to here.
00012 *************************************************************************/
00013 
00014 /*------------------------------------------------------------------------
00015   Inputs: DICOM coordinate warp from old_dset to new_dset
00016   Output: voxel index coordinates warp that can actually be used
00017 --------------------------------------------------------------------------*/
00018 
00019 THD_warp * AFNI_make_voxwarp( THD_warp * inwarp ,
00020                               THD_3dim_dataset * old_dset ,
00021                               THD_3dim_dataset * new_dset  )
00022 {
00023    THD_warp * newwarp ;
00024    THD_linear_mapping * map ;
00025    THD_dataxes * new_daxes ;
00026 
00027    newwarp       = myXtNew( THD_warp ) ;
00028    newwarp->type = inwarp->type ;
00029    new_daxes     = CURRENT_DAXES(new_dset) ;
00030 
00031    switch( inwarp->type ){
00032 
00033       default:{
00034          fprintf(stderr,"\a\n*** ILLEGAL warp code!!! %d\n",inwarp->type) ;
00035          sleep(1) ; EXIT(1) ;
00036       }
00037       break ;
00038 
00039       case WARP_AFFINE_TYPE:{
00040 
00041          map = AFNI_make_voxmap( &(inwarp->rig_bod.warp),
00042                                  old_dset->daxes , new_daxes ) ;
00043 
00044          /* load the (inclusive) voxel index ranges into the affine map */
00045 
00046          LOAD_FVEC3( map->bot, 0,0,0 ) ;
00047          LOAD_FVEC3( map->top, new_daxes->nxx-1,
00048                                new_daxes->nyy-1, new_daxes->nzz-1 ) ;
00049 
00050 
00051          newwarp->rig_bod.warp = *map ;
00052 
00053          myXtFree( map ) ;
00054       }
00055       break ;
00056 
00057       case WARP_TALAIRACH_12_TYPE:{
00058          int iw ;
00059          for( iw=0 ; iw < 12 ; iw++ ){
00060             map = AFNI_make_voxmap( &(inwarp->tal_12.warp[iw]) ,
00061                                     old_dset->daxes , new_daxes ) ;
00062 
00063             map->bot = THD_dicomm_to_3dmm(new_dset,inwarp->tal_12.warp[iw].bot);
00064             map->top = THD_dicomm_to_3dmm(new_dset,inwarp->tal_12.warp[iw].top);
00065 
00066             map->bot = THD_3dmm_to_3dfind( new_dset , map->bot ) ;
00067             map->top = THD_3dmm_to_3dfind( new_dset , map->top ) ;
00068 
00069             newwarp->tal_12.warp[iw] = *map ;
00070 
00071             myXtFree( map ) ;
00072          }
00073       }
00074       break ;
00075    }
00076 
00077    return newwarp ;
00078 }
00079 
00080 /*----------------------------------------------------------------*/
00081 
00082 #ifdef SOLARIS
00083 void do_nothing(int iii){return ;}  /* for compiler optimization error on Sun */
00084 #else
00085 #define do_nothing(iii)
00086 #endif
00087 
00088 /*-----------------------------------------------------------------
00089   This routine takes as input a linear mapping from Dicom to Dicom
00090   coordinates, and computes a linear mapping from voxel-index
00091   to voxel-index coordinates as output.
00092   Pointers to the old and new dataset axes structures are also input.
00093 
00094   x3dmm_new = [new_dicomm_to_3dmm]
00095               ( [dd_trans] [old_3dmm_to_dicomm] x3dmm_old - dd_base )
00096 
00097   is the conversion between true (3dmm) coordinates.  In turn, we also
00098   have
00099 
00100   x3dfind_new = [new_scale] ( x3dmm_new - new_origin )
00101   x3dmm_old   = [old_scale] x3dfind_old + old_origin
00102 
00103   as the transformations between voxel-index (3dfind) coordinates
00104   and true (3dmm) coordinates.
00105 --------------------------------------------------------------------*/
00106 
00107 THD_linear_mapping * AFNI_make_voxmap( THD_linear_mapping * inmap ,
00108                                        THD_dataxes * old_daxes ,
00109                                        THD_dataxes * new_daxes  )
00110 {
00111    THD_mat33 old_scale , old_3dmm_to_dicomm ,
00112              dd_trans  , new_scale , new_dicomm_to_3dmm ,
00113              mt ; /* temp matrix */
00114 
00115    THD_fvec3 dd_base , new_origin , old_origin , vt0,vt1,vt2 ;
00116 
00117    THD_linear_mapping * newmap ;
00118 
00119    /*--- set up the elements of the transformation ---*/
00120 
00121    dd_trans = inmap->mfor ;
00122 
00123    LOAD_DIAG_MAT( old_scale , old_daxes->xxdel ,
00124                               old_daxes->yydel ,
00125                               old_daxes->zzdel  ) ;
00126 
00127    LOAD_DIAG_MAT( new_scale , 1.0/new_daxes->xxdel ,     /* notice */
00128                               1.0/new_daxes->yydel ,     /* inversion */
00129                               1.0/new_daxes->zzdel  ) ;
00130 
00131    old_3dmm_to_dicomm = old_daxes->to_dicomm ;
00132    new_dicomm_to_3dmm = TRANSPOSE_MAT(new_daxes->to_dicomm) ; /* inversion */
00133 
00134    /* vector elements */
00135 
00136    dd_base = inmap->bvec ;                        /* in new dicomm */
00137 
00138    LOAD_FVEC3( new_origin , new_daxes->xxorg ,    /* in old 3dmm */
00139                             new_daxes->yyorg ,
00140                             new_daxes->zzorg  ) ;
00141 
00142    LOAD_FVEC3( old_origin , old_daxes->xxorg ,    /* in new 3dmm */
00143                             old_daxes->yyorg ,
00144                             old_daxes->zzorg  ) ;
00145 
00146    /* multiply the matrices together */
00147 
00148    mt = MAT_MUL( old_3dmm_to_dicomm , old_scale ) ;   do_nothing(0) ;
00149    mt = MAT_MUL( dd_trans , mt ) ;                    do_nothing(0) ;
00150    mt = MAT_MUL( new_dicomm_to_3dmm , mt ) ;          do_nothing(0) ;
00151    mt = MAT_MUL( new_scale , mt ) ;                   do_nothing(0) ;
00152 
00153    /* compute the new bvec */
00154 
00155    vt0 = MATVEC( old_3dmm_to_dicomm , old_origin ) ;
00156    vt0 = MATVEC( dd_trans , vt0 ) ;
00157    vt0 = MATVEC( new_dicomm_to_3dmm , vt0 ) ;
00158    vt0 = MATVEC( new_scale , vt0 ) ;
00159 
00160    vt1 = MATVEC( new_dicomm_to_3dmm , dd_base ) ;
00161    vt1 = MATVEC( new_scale , vt1 ) ;
00162 
00163    vt2 = MATVEC( new_scale , new_origin ) ;  /* want vt1 + vt2 - vt0 */
00164 
00165    vt2 = ADD_FVEC3( vt1 , vt2 ) ;
00166    vt2 = SUB_FVEC3( vt2 , vt0 ) ;
00167 
00168    /* make the output map */
00169 
00170    newmap = myXtNew( THD_linear_mapping ) ;
00171 
00172    newmap->type = MAPPING_LINEAR_TYPE ;
00173    newmap->mfor = mt ;
00174    newmap->mbac = MAT_INV(mt) ;
00175    newmap->bvec = vt2 ;
00176 
00177    newmap->svec = MATVEC(newmap->mbac,newmap->bvec) ;
00178    NEGATE_FVEC3(newmap->svec) ;
00179 
00180    return newmap ;
00181 }
00182 
00183 /*--------------------------------------------------------------------
00184   this routine takes as input two warps, A and B, and produces the
00185   warp A*B;  B is the warp that occurs before A.
00186   Restriction: A and B cannot both be Talairach_12 warps!
00187   The output is placed into A (*warp_in).
00188 
00189   Note that a NULL warp pointer for warp_prior will have the same
00190   effect as an identity warp, since the routine will return immediately.
00191 ----------------------------------------------------------------------*/
00192 
00193 void AFNI_concatenate_warp( THD_warp * warp_in , THD_warp * warp_prior )
00194 {
00195    THD_linear_mapping * prior_map , * new_map ;
00196 
00197    if( warp_in == NULL || warp_prior == NULL ) return ;
00198 
00199    switch( warp_in->type + 100*warp_prior->type ){
00200 
00201       default:
00202          warp_in->type = -1 ;  /* set error flag! */
00203          return ;
00204 
00205       /*-- 2 affine warps ==> a new affine warp --*/
00206 
00207       case WARP_AFFINE_TYPE + 100*WARP_AFFINE_TYPE:{
00208          prior_map = &(warp_prior->rig_bod.warp) ;
00209          new_map = AFNI_concatenate_lmap(
00210                       &(warp_in->rig_bod.warp) , prior_map ) ;
00211 
00212          warp_in->rig_bod.warp = *new_map ;  /* write over input warp */
00213          myXtFree( new_map ) ;
00214       }
00215       break ;
00216 
00217       /*--- Talairach preceeded by affine ==> new Talairach --*/
00218 
00219       case WARP_TALAIRACH_12_TYPE + 100*WARP_AFFINE_TYPE:{
00220          int iw ;
00221          prior_map = &(warp_prior->rig_bod.warp) ;
00222          for( iw=0 ; iw < 12 ; iw++ ){
00223 
00224             new_map = AFNI_concatenate_lmap(
00225                          &(warp_in->tal_12.warp[iw]) , prior_map ) ;
00226 
00227             warp_in->tal_12.warp[iw] = *new_map ;  /* write over input warp */
00228             myXtFree( new_map ) ;
00229          }
00230       }
00231       break ;
00232 
00233       /*-- affine preceeded by Talairach ==> new Talairach
00234            [this case is not currently used, since there are no warps
00235             AFTER a Talairach warp, but it may be useful in the future]
00236                                         -- RWCox, November 1994 A.D. --*/
00237 
00238       case WARP_AFFINE_TYPE + 100*WARP_TALAIRACH_12_TYPE:{
00239          int iw ;
00240          THD_talairach_12_warp * new_warp = myXtNew( THD_talairach_12_warp ) ;
00241 
00242          new_warp->type = WARP_TALAIRACH_12_TYPE ;
00243          for( iw=0 ; iw < 12 ; iw++ ){
00244             prior_map = &(warp_prior->tal_12.warp[iw]) ;
00245             new_map   = AFNI_concatenate_lmap(
00246                           &(warp_in->rig_bod.warp) , prior_map ) ;
00247             new_warp->warp[iw] = *new_map ;
00248             myXtFree( new_map ) ;
00249          }
00250 
00251          warp_in->tal_12 = *new_warp ;  /* write over input warp */
00252          myXtFree( new_warp ) ;
00253       }
00254       break ;
00255 
00256    }  /* end of switch on warp types */
00257 
00258    return ;
00259 }
00260 
00261 /*-----------------------------------------------------------------------
00262    Concatenate 2 linear maps (including allowing for shift of origin)
00263 -------------------------------------------------------------------------*/
00264 
00265 THD_linear_mapping * AFNI_concatenate_lmap( THD_linear_mapping * map_2 ,
00266                                             THD_linear_mapping * map_1  )
00267 {
00268    THD_linear_mapping * map_out ;
00269    THD_fvec3 tvec ;
00270 
00271    /* make a new linear mapping */
00272 
00273    map_out = myXtNew(THD_linear_mapping) ;
00274    map_out->type = MAPPING_LINEAR_TYPE ;
00275 
00276    /* matrix */
00277 
00278    map_out->mfor = MAT_MUL( map_2->mfor , map_1->mfor ) ;
00279    map_out->mbac = MAT_INV( map_out->mfor ) ;
00280 
00281    /* vector */
00282 
00283    tvec          = MATVEC( map_2->mfor , map_1->bvec ) ;
00284    map_out->bvec = ADD_FVEC3( tvec , map_2->bvec ) ;
00285    map_out->svec = MATVEC( map_out->mbac , map_out->bvec ) ;
00286    NEGATE_FVEC3(map_out->svec) ;
00287 
00288    map_out->bot  = map_2->bot ;
00289    map_out->top  = map_2->top ;
00290 
00291    return map_out ;
00292 }
00293 
00294 /*--------------------------------------------------------------------------*/
00295 /*! Make an affine warp from 12 input numbers:
00296      -         [ a11 a12 a13 ]        [ s1 ]
00297      - x_map = [ a21 a22 a23 ] x_in + [ s2 ]
00298      -         [ a31 a32 a33 ]        [ s3 ]
00299 
00300     27 Aug 2002 - RWCox.
00301 ----------------------------------------------------------------------------*/
00302 
00303 THD_warp * AFNI_make_affwarp_12( float a11, float a12, float a13,  float s1 ,
00304                                  float a21, float a22, float a23,  float s2 ,
00305                                  float a31, float a32, float a33,  float s3  )
00306 {
00307    THD_warp *warp ;
00308    THD_linear_mapping map ;
00309    float dd , nn ;
00310 
00311    warp       = myXtNew( THD_warp ) ;
00312    warp->type = WARP_AFFINE_TYPE ;
00313 
00314    map.type = MAPPING_LINEAR_TYPE ;
00315 
00316    LOAD_MAT(map.mfor,a11,a12,a13,a21,a22,a23,a31,a32,a33) ;
00317    dd = MAT_DET(map.mfor) ; nn = MAT_FNORM(map.mfor) ;
00318    if( fabs(dd) < 1.e-5*nn*nn*nn ) return NULL ;  /* bad input */
00319    LOAD_FVEC3(map.bvec,-s1,-s2,-s3) ;
00320    LOAD_INVERSE_LMAP(map) ;
00321 
00322    warp->rig_bod.warp = map ;
00323 
00324    return warp ;
00325 }
00326 
00327 /*-------------------------------------------------------------------------*/
00328 
00329 THD_warp * AFNI_make_affwarp_mat( THD_mat33 mmm )
00330 {
00331    return AFNI_make_affwarp_12( mmm.mat[0][0], mmm.mat[0][1], mmm.mat[0][2], 0.0 ,
00332                                 mmm.mat[1][0], mmm.mat[1][1], mmm.mat[1][2], 0.0 ,
00333                                 mmm.mat[2][0], mmm.mat[2][1], mmm.mat[2][2], 0.0  ) ;
00334 }
00335 
00336 /*-------------------------------------------------------------------------*/
00337 
00338 THD_warp * AFNI_make_affwarp_matvec( THD_mat33 mmm , THD_fvec3 vvv )
00339 {
00340    return AFNI_make_affwarp_12( mmm.mat[0][0], mmm.mat[0][1], mmm.mat[0][2], vvv.xyz[0] ,
00341                                 mmm.mat[1][0], mmm.mat[1][1], mmm.mat[1][2], vvv.xyz[1] ,
00342                                 mmm.mat[2][0], mmm.mat[2][1], mmm.mat[2][2], vvv.xyz[2]  ) ;
00343 }