Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
afni_warp.c

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 #undef MAIN
00008 
00009 /****** This file contains routines related to warping datasets ******/
00010 /****** and getting pieces of them into MRI_IMAGE structures.   ******/
00011 
00012 #include "afni_warp.h"
00013 
00014 /*------------------------------------------------------------------------
00015    Return a slice from a dataset, possibly to be warped on-the-fly
00016    from its parent:
00017      dset        = input dataset
00018      fixed_axis  = 1, 2, or 3 (for x,y,z)
00019      fixed_index = value of frozen index
00020      ival        = index of value (or sub-brick) to extract
00021      resam_mode  = used if "warp-on-the-fly" is used
00022 
00023    Output is in the MRI_IMAGE format, and will be
00024       nyy * nzz for fixed_axis = 1
00025       nzz * nxx for fixed_axis = 2
00026       nxx * nyy for fixed_axis = 3
00027 --------------------------------------------------------------------------*/
00028 
00029 MRI_IMAGE * AFNI_dataset_slice( THD_3dim_dataset * dset ,
00030                                 int fixed_axis , int fixed_index ,
00031                                 int ival , int resam_mode )
00032 {
00033    MRI_IMAGE * newim ;
00034    void * sar , * bar ;
00035    int nxx,nyy,nzz ;
00036    MRI_TYPE typ ;
00037    THD_dataxes * daxes ;
00038    THD_3dim_dataset * parent_dset ;
00039    THD_warp parent_to_child_warp ;
00040 
00041 ENTRY("AFNI_dataset_slice") ;
00042 
00043    /*--------- sanity checks ---------*/
00044 
00045    if( ! ISVALID_3DIM_DATASET(dset) || fixed_index < 0 ||
00046        ival >= dset->dblk->nvals                         ) RETURN(NULL) ;
00047 
00048    /** must allow for ival < 0 --> return an empty (not NULL) image **/
00049 
00050    typ = (ival < 0 ) ? MRI_short : DSET_BRICK_TYPE(dset,ival) ;
00051    if( ! AFNI_GOOD_DTYPE(typ) ) RETURN(NULL) ;
00052 
00053    /*--------- brick dimensions ---------*/
00054 
00055    daxes = CURRENT_DAXES(dset) ;
00056    nxx   = daxes->nxx ;
00057    nyy   = daxes->nyy ;
00058    nzz   = daxes->nzz ;
00059 
00060 #if 0
00061    /*** June 1996: if resampling a threshold, always use NN ***/
00062 
00063    if( ival == DSET_THRESH_VALUE(dset) ) resam_mode = RESAM_NN_TYPE ;
00064 #endif
00065 
00066 if(PRINT_TRACING)
00067 { char str[256] ;
00068   sprintf(str,"Input dataset = %s",DSET_FILECODE(dset)) ; STATUS(str) ;
00069   sprintf(str,"nxx=%d nyy=%d nzz=%d",nxx,nyy,nzz) ;  STATUS(str) ;
00070   sprintf(str,"fixed_axis=%d fixed_index=%d ival=%d resam=%d",
00071           fixed_axis,fixed_index,ival,resam_mode ) ; STATUS(str) ; }
00072 
00073    /*--------- setup output image ---------*/
00074 
00075    switch( fixed_axis ){
00076       default: RETURN(NULL) ;  /* should not happen */
00077 
00078       /* fixed x --> image is y-z */
00079 
00080       case 1:
00081          if( fixed_index >= nxx ) RETURN(NULL) ;
00082          newim     = mri_new( nyy , nzz , typ ) ;
00083          newim->dx = fabs(daxes->yydel) ;
00084          newim->dy = fabs(daxes->zzdel) ;
00085          newim->dz = fabs(daxes->xxdel) ;
00086          sar       = mri_data_pointer( newim ) ;
00087       break ;
00088 
00089       /* fixed y --> image is z-x */
00090 
00091       case 2:
00092          if( fixed_index >= nyy ) RETURN(NULL) ;
00093          newim     = mri_new( nzz , nxx , typ ) ;
00094          newim->dx = fabs(daxes->zzdel) ;
00095          newim->dy = fabs(daxes->xxdel) ;
00096          newim->dz = fabs(daxes->yydel) ;
00097          sar       = mri_data_pointer( newim ) ;
00098       break ;
00099 
00100       /* fixed z --> image is x-y */
00101 
00102       case 3:
00103          if( fixed_index >= nzz ) RETURN(NULL) ;
00104          newim     = mri_new( nxx , nyy , typ ) ;
00105          newim->dx = fabs(daxes->xxdel) ;
00106          newim->dy = fabs(daxes->yydel) ;
00107          newim->dz = fabs(daxes->zzdel) ;
00108          sar       = mri_data_pointer( newim ) ;
00109       break ;
00110    }
00111 
00112    /** return empty image? **/
00113 
00114    if( ival < 0 ) RETURN(newim) ;
00115 
00116 if(PRINT_TRACING)
00117 { char str[256] ;
00118   sprintf(str,"newim nx=%d ny=%d dx=%f dy=%f",
00119           newim->nx , newim->ny , newim->dx , newim->dy ) ;
00120   STATUS(str) ; }
00121 
00122    /*----- if datablock exists and not forcing warp-on-demand, use it -----*/
00123 
00124    if( !dset->wod_flag && DSET_INMEMORY(dset) ){
00125 
00126       bar = DSET_ARRAY(dset,ival) ;  /* pointer to data brick array */
00127 
00128       if( bar == NULL ){  /* if data needs to be loaded from disk */
00129          (void) THD_load_datablock( dset->dblk ) ;
00130          bar = DSET_ARRAY(dset,ival) ;
00131          if( bar == NULL ){
00132             STATUS("couldn't load dataset!") ;
00133             mri_free(newim) ;
00134             RETURN(NULL) ;  /* couldn't load data --> return nothing */
00135          }
00136       }
00137 
00138 STATUS("reading from memory") ;
00139 
00140       switch( typ ){  /* copy data from brick (bar) to newim (sar) */
00141 
00142          default:              /* should not happen! */
00143             mri_free(newim) ;
00144             RETURN(NULL) ;
00145 
00146          case MRI_short:{
00147             AFNI_br2sl_short( nxx,nyy,nzz , fixed_axis,fixed_index ,
00148                               (short *) bar , (short *) sar         ) ;
00149 
00150             if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00151                MRI_IMAGE * qim ;
00152                qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00153                mri_free(newim) ; newim = qim ;
00154 STATUS("scaling slice to floats") ;
00155             }
00156             RETURN(newim) ;
00157          }
00158 
00159          case MRI_float:{
00160             AFNI_br2sl_float( nxx,nyy,nzz , fixed_axis,fixed_index ,
00161                               (float *) bar , (float *) sar         ) ;
00162 
00163             if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00164                MRI_IMAGE * qim ;
00165                qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00166                mri_free(newim) ; newim = qim ;
00167 STATUS("scaling slice to floats") ;
00168             }
00169             RETURN(newim) ;
00170          }
00171 
00172          case MRI_byte:{
00173             AFNI_br2sl_byte( nxx,nyy,nzz , fixed_axis,fixed_index ,
00174                              (byte *) bar , (byte *) sar           ) ;
00175 
00176             if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00177                MRI_IMAGE * qim ;
00178                qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00179                mri_free(newim) ; newim = qim ;
00180 STATUS("scaling slice to floats") ;
00181             }
00182             RETURN(newim) ;
00183          }
00184 
00185          case MRI_complex:{
00186             AFNI_br2sl_complex( nxx,nyy,nzz , fixed_axis,fixed_index ,
00187                                 (complex *) bar , (complex *) sar         ) ;
00188             RETURN(newim) ;
00189          }
00190 
00191          case MRI_rgb:{
00192             AFNI_br2sl_rgbyte( nxx,nyy,nzz , fixed_axis,fixed_index ,
00193                                 (rgbyte *) bar , (rgbyte *) sar         ) ;
00194             RETURN(newim) ;
00195          }
00196       }
00197    } /* end of if on dataset brick existing */
00198 
00199    /*------------------ must warp from parent dataset -------------------*/
00200 
00201    if( dset->warp != NULL ){
00202 STATUS("setting parent_to_child_warp to stored warp") ;
00203       parent_to_child_warp = *(dset->warp) ;
00204    } else {
00205 STATUS("setting parent_to_child_warp to identity") ;
00206       parent_to_child_warp = IDENTITY_WARP ;  /* no warp ==> use identity */
00207    }
00208 
00209    if( dset->warp_parent != NULL ){
00210 if(PRINT_TRACING)
00211 { char str[256] ;
00212   sprintf(str,"setting parent_dset to stored warp_parent=%p  this dset=%p",
00213           (void *) dset->warp_parent , (void *) dset ) ; STATUS(str) ; }
00214 
00215       parent_dset = dset->warp_parent ;
00216    } else {
00217 STATUS("setting parent_dset to self, and parent_to_child_warp to identity") ;
00218       parent_dset = dset ;                    /* self-parenting */
00219 
00220       if( dset->self_warp != NULL )
00221         parent_to_child_warp = *(dset->self_warp) ;  /* 26 Aug 2002 */
00222       else
00223         parent_to_child_warp = IDENTITY_WARP ;  /* use identity warp */
00224    }
00225 
00226    /*----- make the voxel-to-voxel warp, if needed -----*/
00227 
00228    if( ! ISVALID_WARP(dset->vox_warp) ){
00229       THD_warp * qwarp ;
00230       qwarp = AFNI_make_voxwarp( &parent_to_child_warp, parent_dset, dset ) ;
00231 
00232       if( dset->vox_warp == NULL ){    /* totally new */
00233          dset->vox_warp = qwarp ;
00234          ADDTO_KILL(dset->kl,dset->vox_warp) ;
00235       } else {
00236          *(dset->vox_warp) = *qwarp ;  /* just copy insides */
00237          myXtFree( qwarp ) ;
00238       }
00239    }
00240 
00241    if( DSET_ARRAY(parent_dset,ival) == NULL ){  /* reload from disk */
00242       Boolean good ;
00243       good = THD_load_datablock( parent_dset->dblk ) ;
00244       if( ! good ){
00245          STATUS("couldn't load parent dataset!") ;
00246          mri_free(newim) ;
00247          RETURN(NULL) ;  /* couldn't load data --> return nothing */
00248       }
00249    }
00250 
00251    bar = DSET_ARRAY(parent_dset,ival) ;
00252 
00253 STATUS("warp-on-demand") ;
00254 
00255    /******************************************************************/
00256    /*** Select warp routine based on data type and slice direction ***/
00257    /**********   [See template routines in afni_slice.c]    **********/
00258    /******************************************************************/
00259 
00260 #undef USE_CLIP  /* 10 Dec 1997 -- CUBIC_CLIP is broken */
00261 
00262    switch( typ ){
00263 
00264       default:               /** Illegal type: should not happen! **/
00265          mri_free(newim) ;
00266          RETURN(NULL) ;
00267 
00268       /**************************** short ****************************/
00269 #undef  DTYPE
00270 #undef  LMAP_XNAME
00271 #undef  LMAP_YNAME
00272 #undef  LMAP_ZNAME
00273 #undef  CUBIC_CLIP
00274 #define DTYPE      short
00275 #define LMAP_XNAME TWO_TWO(AFNI_lmap_to_xslice_,DTYPE)
00276 #define LMAP_YNAME TWO_TWO(AFNI_lmap_to_yslice_,DTYPE)
00277 #define LMAP_ZNAME TWO_TWO(AFNI_lmap_to_zslice_,DTYPE)
00278 #ifdef USE_CLIP
00279 #define CUBIC_CLIP                                                         \
00280    if( resam_mode == RESAM_CUBIC_TYPE && ISVALID_STATISTIC(dset->stats) ){ \
00281      int ii , npix = newim->nx * newim->ny ;                               \
00282      DTYPE * ar = mri_data_pointer(newim) ;                                \
00283      DTYPE bot = dset->stats->bstat[ival].min ,                            \
00284            top = dset->stats->bstat[ival].max ;                            \
00285      if( bot < top ) for( ii=0 ; ii < npix ; ii++ )                        \
00286                         if( ar[ii] < bot ) ar[ii] = bot ;                  \
00287                    else if( ar[ii] > top ) ar[ii] = top ; }
00288 #else
00289 #define CUBIC_CLIP /* nada */
00290 #endif
00291 
00292       case TWO_TWO(MRI_,DTYPE):
00293       switch( fixed_axis ){
00294 
00295          case 1:{
00296 
00297             switch( dset->vox_warp->type ){
00298 
00299                case WARP_AFFINE_TYPE:{
00300                   LMAP_XNAME( &(dset->vox_warp->rig_bod.warp) ,
00301                               resam_mode ,
00302                               parent_dset->daxes ,
00303                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00304                }
00305                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00306                   MRI_IMAGE * qim ;
00307                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00308                   mri_free(newim) ; newim = qim ;
00309 STATUS("scaling slice to floats") ;
00310                }
00311                CUBIC_CLIP ;
00312                RETURN(newim) ;
00313 
00314                case WARP_TALAIRACH_12_TYPE:{
00315                   int iw ;
00316                   for( iw=0 ; iw < 12 ; iw++ )
00317                      LMAP_XNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00318                                  resam_mode ,
00319                                  parent_dset->daxes ,
00320                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00321                   }
00322                }
00323                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00324                   MRI_IMAGE * qim ;
00325                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00326                   mri_free(newim) ; newim = qim ;
00327 STATUS("scaling slice to floats") ;
00328                }
00329                CUBIC_CLIP ;
00330                RETURN(newim) ;
00331          }
00332          break ;
00333 
00334          case 2:{
00335 
00336             switch( dset->vox_warp->type ){
00337 
00338                case WARP_AFFINE_TYPE:{
00339                   LMAP_YNAME( &(dset->vox_warp->rig_bod.warp) ,
00340                               resam_mode ,
00341                               parent_dset->daxes ,
00342                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00343                }
00344                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00345                   MRI_IMAGE * qim ;
00346                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00347                   mri_free(newim) ; newim = qim ;
00348 STATUS("scaling slice to floats") ;
00349                }
00350                CUBIC_CLIP ;
00351                RETURN(newim) ;
00352 
00353                case WARP_TALAIRACH_12_TYPE:{
00354                   int iw ;
00355                   for( iw=0 ; iw < 12 ; iw++ )
00356                      LMAP_YNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00357                                  resam_mode ,
00358                                  parent_dset->daxes ,
00359                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00360                   }
00361                }
00362                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00363                   MRI_IMAGE * qim ;
00364                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00365                   mri_free(newim) ; newim = qim ;
00366 STATUS("scaling slice to floats") ;
00367                }
00368                CUBIC_CLIP ;
00369                RETURN(newim) ;
00370          }
00371          break ;
00372 
00373          case 3:{
00374 
00375             switch( dset->vox_warp->type ){
00376 
00377                case WARP_AFFINE_TYPE:{
00378                   LMAP_ZNAME( &(dset->vox_warp->rig_bod.warp) ,
00379                               resam_mode ,
00380                               parent_dset->daxes ,
00381                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00382                }
00383                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00384                   MRI_IMAGE * qim ;
00385                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00386                   mri_free(newim) ; newim = qim ;
00387 STATUS("scaling slice to floats") ;
00388                }
00389                CUBIC_CLIP ;
00390                RETURN(newim) ;
00391 
00392                case WARP_TALAIRACH_12_TYPE:{
00393                   int iw ;
00394                   for( iw=0 ; iw < 12 ; iw++ )
00395                      LMAP_ZNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00396                                  resam_mode ,
00397                                  parent_dset->daxes ,
00398                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00399                   }
00400                }
00401                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00402                   MRI_IMAGE * qim ;
00403                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00404                   mri_free(newim) ; newim = qim ;
00405 STATUS("scaling slice to floats") ;
00406                }
00407                CUBIC_CLIP ;
00408                RETURN(newim) ;
00409          }
00410          break ;
00411       }
00412 
00413       /**************************** float ****************************/
00414 #undef  DTYPE
00415 #undef  LMAP_XNAME
00416 #undef  LMAP_YNAME
00417 #undef  LMAP_ZNAME
00418 #undef  CUBIC_CLIP
00419 #define DTYPE      float
00420 #define LMAP_XNAME TWO_TWO(AFNI_lmap_to_xslice_,DTYPE)
00421 #define LMAP_YNAME TWO_TWO(AFNI_lmap_to_yslice_,DTYPE)
00422 #define LMAP_ZNAME TWO_TWO(AFNI_lmap_to_zslice_,DTYPE)
00423 #ifdef USE_CLIP
00424 #define CUBIC_CLIP                                                         \
00425    if( resam_mode == RESAM_CUBIC_TYPE && ISVALID_STATISTIC(dset->stats) ){ \
00426      int ii , npix = newim->nx * newim->ny ;                               \
00427      DTYPE * ar = mri_data_pointer(newim) ;                                \
00428      DTYPE bot = dset->stats->bstat[ival].min ,                            \
00429            top = dset->stats->bstat[ival].max ;                            \
00430      if( bot < top ) for( ii=0 ; ii < npix ; ii++ )                        \
00431                         if( ar[ii] < bot ) ar[ii] = bot ;                  \
00432                    else if( ar[ii] > top ) ar[ii] = top ; }
00433 #else
00434 #define CUBIC_CLIP /* nada */
00435 #endif
00436 
00437       case TWO_TWO(MRI_,DTYPE):
00438       switch( fixed_axis ){
00439 
00440          case 1:{
00441 
00442             switch( dset->vox_warp->type ){
00443 
00444                case WARP_AFFINE_TYPE:{
00445                   LMAP_XNAME( &(dset->vox_warp->rig_bod.warp) ,
00446                               resam_mode ,
00447                               parent_dset->daxes ,
00448                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00449                }
00450                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00451                   MRI_IMAGE * qim ;
00452                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00453                   mri_free(newim) ; newim = qim ;
00454 STATUS("scaling slice to floats") ;
00455                }
00456                CUBIC_CLIP ;
00457                RETURN(newim) ;
00458 
00459                case WARP_TALAIRACH_12_TYPE:{
00460                   int iw ;
00461                   for( iw=0 ; iw < 12 ; iw++ )
00462                      LMAP_XNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00463                                  resam_mode ,
00464                                  parent_dset->daxes ,
00465                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00466                   }
00467                }
00468                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00469                   MRI_IMAGE * qim ;
00470                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00471                   mri_free(newim) ; newim = qim ;
00472 STATUS("scaling slice to floats") ;
00473                }
00474                CUBIC_CLIP ;
00475                RETURN(newim) ;
00476          }
00477          break ;
00478 
00479          case 2:{
00480 
00481             switch( dset->vox_warp->type ){
00482 
00483                case WARP_AFFINE_TYPE:{
00484                   LMAP_YNAME( &(dset->vox_warp->rig_bod.warp) ,
00485                               resam_mode ,
00486                               parent_dset->daxes ,
00487                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00488                }
00489                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00490                   MRI_IMAGE * qim ;
00491                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00492                   mri_free(newim) ; newim = qim ;
00493 STATUS("scaling slice to floats") ;
00494                }
00495                CUBIC_CLIP ;
00496                RETURN(newim) ;
00497 
00498                case WARP_TALAIRACH_12_TYPE:{
00499                   int iw ;
00500                   for( iw=0 ; iw < 12 ; iw++ )
00501                      LMAP_YNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00502                                  resam_mode ,
00503                                  parent_dset->daxes ,
00504                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00505                   }
00506                }
00507                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00508                   MRI_IMAGE * qim ;
00509                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00510                   mri_free(newim) ; newim = qim ;
00511 STATUS("scaling slice to floats") ;
00512                }
00513                CUBIC_CLIP ;
00514                RETURN(newim) ;
00515          }
00516          break ;
00517 
00518          case 3:{
00519 
00520             switch( dset->vox_warp->type ){
00521 
00522                case WARP_AFFINE_TYPE:{
00523                   LMAP_ZNAME( &(dset->vox_warp->rig_bod.warp) ,
00524                               resam_mode ,
00525                               parent_dset->daxes ,
00526                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00527                }
00528                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00529                   MRI_IMAGE * qim ;
00530                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00531                   mri_free(newim) ; newim = qim ;
00532 STATUS("scaling slice to floats") ;
00533                }
00534                CUBIC_CLIP ;
00535                RETURN(newim) ;
00536 
00537                case WARP_TALAIRACH_12_TYPE:{
00538                   int iw ;
00539                   for( iw=0 ; iw < 12 ; iw++ )
00540                      LMAP_ZNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00541                                  resam_mode ,
00542                                  parent_dset->daxes ,
00543                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00544                   }
00545                }
00546                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00547                   MRI_IMAGE * qim ;
00548                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00549                   mri_free(newim) ; newim = qim ;
00550 STATUS("scaling slice to floats") ;
00551                }
00552                CUBIC_CLIP ;
00553                RETURN(newim) ;
00554          }
00555          break ;
00556       }
00557 
00558       /**************************** byte ****************************/
00559 #undef  DTYPE
00560 #undef  LMAP_XNAME
00561 #undef  LMAP_YNAME
00562 #undef  LMAP_ZNAME
00563 #undef  CUBIC_CLIP
00564 #define DTYPE      byte
00565 #define LMAP_XNAME TWO_TWO(AFNI_lmap_to_xslice_,DTYPE)
00566 #define LMAP_YNAME TWO_TWO(AFNI_lmap_to_yslice_,DTYPE)
00567 #define LMAP_ZNAME TWO_TWO(AFNI_lmap_to_zslice_,DTYPE)
00568 #ifdef USE_CLIP
00569 #define CUBIC_CLIP                                                         \
00570    if( resam_mode == RESAM_CUBIC_TYPE && ISVALID_STATISTIC(dset->stats) ){ \
00571      int ii , npix = newim->nx * newim->ny ;                               \
00572      DTYPE * ar = mri_data_pointer(newim) ;                                \
00573      DTYPE bot = dset->stats->bstat[ival].min ,                            \
00574            top = dset->stats->bstat[ival].max ;                            \
00575      if( bot < top ) for( ii=0 ; ii < npix ; ii++ )                        \
00576                         if( ar[ii] < bot ) ar[ii] = bot ;                  \
00577                    else if( ar[ii] > top ) ar[ii] = top ; }
00578 #else
00579 #define CUBIC_CLIP /* nada */
00580 #endif
00581 
00582       case TWO_TWO(MRI_,DTYPE):
00583       switch( fixed_axis ){
00584 
00585          case 1:{
00586 
00587             switch( dset->vox_warp->type ){
00588 
00589                case WARP_AFFINE_TYPE:{
00590                   LMAP_XNAME( &(dset->vox_warp->rig_bod.warp) ,
00591                               resam_mode ,
00592                               parent_dset->daxes ,
00593                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00594                }
00595                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00596                   MRI_IMAGE * qim ;
00597                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00598                   mri_free(newim) ; newim = qim ;
00599 STATUS("scaling slice to floats") ;
00600                }
00601                CUBIC_CLIP ;
00602                RETURN(newim) ;
00603 
00604                case WARP_TALAIRACH_12_TYPE:{
00605                   int iw ;
00606                   for( iw=0 ; iw < 12 ; iw++ )
00607                      LMAP_XNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00608                                  resam_mode ,
00609                                  parent_dset->daxes ,
00610                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00611                   }
00612                }
00613                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00614                   MRI_IMAGE * qim ;
00615                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00616                   mri_free(newim) ; newim = qim ;
00617 STATUS("scaling slice to floats") ;
00618                }
00619                CUBIC_CLIP ;
00620                RETURN(newim) ;
00621          }
00622          break ;
00623 
00624          case 2:{
00625 
00626             switch( dset->vox_warp->type ){
00627 
00628                case WARP_AFFINE_TYPE:{
00629                   LMAP_YNAME( &(dset->vox_warp->rig_bod.warp) ,
00630                               resam_mode ,
00631                               parent_dset->daxes ,
00632                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00633                }
00634                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00635                   MRI_IMAGE * qim ;
00636                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00637                   mri_free(newim) ; newim = qim ;
00638 STATUS("scaling slice to floats") ;
00639                }
00640                CUBIC_CLIP ;
00641                RETURN(newim) ;
00642 
00643                case WARP_TALAIRACH_12_TYPE:{
00644                   int iw ;
00645                   for( iw=0 ; iw < 12 ; iw++ )
00646                      LMAP_YNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00647                                  resam_mode ,
00648                                  parent_dset->daxes ,
00649                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00650                   }
00651                }
00652                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00653                   MRI_IMAGE * qim ;
00654                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00655                   mri_free(newim) ; newim = qim ;
00656 STATUS("scaling slice to floats") ;
00657                }
00658                CUBIC_CLIP ;
00659                RETURN(newim) ;
00660          }
00661          break ;
00662 
00663          case 3:{
00664 
00665             switch( dset->vox_warp->type ){
00666 
00667                case WARP_AFFINE_TYPE:{
00668                   LMAP_ZNAME( &(dset->vox_warp->rig_bod.warp) ,
00669                               resam_mode ,
00670                               parent_dset->daxes ,
00671                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00672                }
00673                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00674                   MRI_IMAGE * qim ;
00675                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00676                   mri_free(newim) ; newim = qim ;
00677 STATUS("scaling slice to floats") ;
00678                }
00679                CUBIC_CLIP ;
00680                RETURN(newim) ;
00681 
00682                case WARP_TALAIRACH_12_TYPE:{
00683                   int iw ;
00684                   for( iw=0 ; iw < 12 ; iw++ )
00685                      LMAP_ZNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00686                                  resam_mode ,
00687                                  parent_dset->daxes ,
00688                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00689                   }
00690                }
00691                if( DSET_BRICK_FACTOR(dset,ival) != 0.0 ){
00692                   MRI_IMAGE * qim ;
00693                   qim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,ival) , newim ) ;
00694                   mri_free(newim) ; newim = qim ;
00695 STATUS("scaling slice to floats") ;
00696                }
00697                CUBIC_CLIP ;
00698                RETURN(newim) ;
00699          }
00700          break ;
00701       }
00702 
00703       /**************************** complex ****************************/
00704 #undef  DTYPE
00705 #undef  LMAP_XNAME
00706 #undef  LMAP_YNAME
00707 #undef  LMAP_ZNAME
00708 #undef  CUBIC_CLIP
00709 #define DTYPE      complex
00710 #define LMAP_XNAME TWO_TWO(AFNI_lmap_to_xslice_,DTYPE)
00711 #define LMAP_YNAME TWO_TWO(AFNI_lmap_to_yslice_,DTYPE)
00712 #define LMAP_ZNAME TWO_TWO(AFNI_lmap_to_zslice_,DTYPE)
00713 
00714       case TWO_TWO(MRI_,DTYPE):
00715       switch( fixed_axis ){
00716 
00717          case 1:{
00718 
00719             switch( dset->vox_warp->type ){
00720 
00721                case WARP_AFFINE_TYPE:{
00722                   LMAP_XNAME( &(dset->vox_warp->rig_bod.warp) ,
00723                               resam_mode ,
00724                               parent_dset->daxes ,
00725                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00726                }
00727                RETURN(newim) ;
00728 
00729                case WARP_TALAIRACH_12_TYPE:{
00730                   int iw ;
00731                   for( iw=0 ; iw < 12 ; iw++ )
00732                      LMAP_XNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00733                                  resam_mode ,
00734                                  parent_dset->daxes ,
00735                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00736                   }
00737                }
00738                RETURN(newim) ;
00739          }
00740          break ;
00741 
00742          case 2:{
00743 
00744             switch( dset->vox_warp->type ){
00745 
00746                case WARP_AFFINE_TYPE:{
00747                   LMAP_YNAME( &(dset->vox_warp->rig_bod.warp) ,
00748                               resam_mode ,
00749                               parent_dset->daxes ,
00750                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00751                }
00752                RETURN(newim) ;
00753 
00754                case WARP_TALAIRACH_12_TYPE:{
00755                   int iw ;
00756                   for( iw=0 ; iw < 12 ; iw++ )
00757                      LMAP_YNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00758                                  resam_mode ,
00759                                  parent_dset->daxes ,
00760                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00761                   }
00762                }
00763                RETURN(newim) ;
00764          }
00765          break ;
00766 
00767          case 3:{
00768 
00769             switch( dset->vox_warp->type ){
00770 
00771                case WARP_AFFINE_TYPE:{
00772                   LMAP_ZNAME( &(dset->vox_warp->rig_bod.warp) ,
00773                               resam_mode ,
00774                               parent_dset->daxes ,
00775                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00776                }
00777                RETURN(newim) ;
00778 
00779                case WARP_TALAIRACH_12_TYPE:{
00780                   int iw ;
00781                   for( iw=0 ; iw < 12 ; iw++ )
00782                      LMAP_ZNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00783                                  resam_mode ,
00784                                  parent_dset->daxes ,
00785                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00786                   }
00787                }
00788                RETURN(newim) ;
00789          }
00790          break ;
00791       }
00792 
00793       /**************************** rgb ****************************/
00794 #undef  DTYPE
00795 #undef  LMAP_XNAME
00796 #undef  LMAP_YNAME
00797 #undef  LMAP_ZNAME
00798 #undef  CUBIC_CLIP
00799 #define DTYPE      rgbyte
00800 #define LMAP_XNAME TWO_TWO(AFNI_lmap_to_xslice_,DTYPE)
00801 #define LMAP_YNAME TWO_TWO(AFNI_lmap_to_yslice_,DTYPE)
00802 #define LMAP_ZNAME TWO_TWO(AFNI_lmap_to_zslice_,DTYPE)
00803 
00804       case TWO_TWO(MRI_,DTYPE):
00805       switch( fixed_axis ){
00806 
00807          case 1:{
00808 
00809             switch( dset->vox_warp->type ){
00810 
00811                case WARP_AFFINE_TYPE:{
00812                   LMAP_XNAME( &(dset->vox_warp->rig_bod.warp) ,
00813                               resam_mode ,
00814                               parent_dset->daxes ,
00815                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00816                }
00817                RETURN(newim) ;
00818 
00819                case WARP_TALAIRACH_12_TYPE:{
00820                   int iw ;
00821                   for( iw=0 ; iw < 12 ; iw++ )
00822                      LMAP_XNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00823                                  resam_mode ,
00824                                  parent_dset->daxes ,
00825                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00826                   }
00827                }
00828                RETURN(newim) ;
00829          }
00830          break ;
00831 
00832          case 2:{
00833 
00834             switch( dset->vox_warp->type ){
00835 
00836                case WARP_AFFINE_TYPE:{
00837                   LMAP_YNAME( &(dset->vox_warp->rig_bod.warp) ,
00838                               resam_mode ,
00839                               parent_dset->daxes ,
00840                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00841                }
00842                RETURN(newim) ;
00843 
00844                case WARP_TALAIRACH_12_TYPE:{
00845                   int iw ;
00846                   for( iw=0 ; iw < 12 ; iw++ )
00847                      LMAP_YNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00848                                  resam_mode ,
00849                                  parent_dset->daxes ,
00850                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00851                   }
00852                }
00853                RETURN(newim) ;
00854          }
00855          break ;
00856 
00857          case 3:{
00858 
00859             switch( dset->vox_warp->type ){
00860 
00861                case WARP_AFFINE_TYPE:{
00862                   LMAP_ZNAME( &(dset->vox_warp->rig_bod.warp) ,
00863                               resam_mode ,
00864                               parent_dset->daxes ,
00865                               (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00866                }
00867                RETURN(newim) ;
00868 
00869                case WARP_TALAIRACH_12_TYPE:{
00870                   int iw ;
00871                   for( iw=0 ; iw < 12 ; iw++ )
00872                      LMAP_ZNAME( &(dset->vox_warp->tal_12.warp[iw]) ,
00873                                  resam_mode ,
00874                                  parent_dset->daxes ,
00875                                  (DTYPE *)bar , daxes , fixed_index , (DTYPE *)sar ) ;
00876                   }
00877                }
00878                RETURN(newim) ;
00879          }
00880          break ;
00881       }
00882 
00883       /**************************** DONE ****************************/
00884 #undef  DTYPE
00885 #undef  LMAP_XNAME
00886 #undef  LMAP_YNAME
00887 #undef  LMAP_ZNAME
00888 #undef  CUBIC_CLIP
00889    }  /** end of switch over legal types for warping **/
00890 
00891    /*--- should never reach this code! ---*/
00892 
00893    mri_free(newim) ; RETURN(NULL) ;
00894 }
00895 
00896 /*-------------------------------------------------------------
00897    A version of FD_brick_to_mri that allows on-the-fly warping.
00898    25 Jul 2001 -- moved the flipping stuff to new function
00899    AFNI_slice_flip() -- RWCox
00900 ---------------------------------------------------------------*/
00901 
00902 MRI_IMAGE * FD_warp_to_mri( int kslice , int ival , FD_brick * br )
00903 {
00904    int ax_1 , ax_2 , ax_3 ;
00905    MRI_IMAGE * flim ;
00906    int resam_code ;
00907 
00908 ENTRY("FD_warp_to_mri") ;
00909 
00910    /*--- get the image from the dataset ---*/
00911 
00912    if( br == NULL || kslice < 0 ) RETURN(NULL) ;  /* should not happen! */
00913 
00914    ax_1 = br->a123.ijk[0] ;  /* axis codes for the desired image */
00915    ax_2 = br->a123.ijk[1] ;
00916    ax_3 = br->a123.ijk[2] ;
00917 
00918    resam_code = ( DSET_BRICK_STATCODE(br->dset,ival) > 0 )
00919                 ? br->thr_resam_code
00920                 : br->resam_code ;
00921 
00922 if(PRINT_TRACING){
00923  char str[256] ;
00924  sprintf(str,"thr_resam_code=%d fim_resam_code=%d resam_code=%d",
00925          br->thr_resam_code,br->resam_code,resam_code) ;
00926  STATUS(str); }
00927 
00928    flim = AFNI_slice_flip( kslice , ival , resam_code ,
00929                            ax_1 , ax_2 , ax_3 , br->dset ) ;
00930 
00931    RETURN(flim) ;
00932 }
00933 
00934 /*-----------------------------------------------------------------
00935    Get the #kslice-th slice from sub-brick #ival with resampling
00936    mode resam, and flip it so that it is oriented along directions
00937    (ax_1,ax_2,ax_3) -- 25 Jul 2001 - RWCox
00938 -------------------------------------------------------------------*/
00939 
00940 MRI_IMAGE * AFNI_slice_flip( int kslice , int ival , int resam ,
00941                              int ax_1, int ax_2, int ax_3 ,
00942                              THD_3dim_dataset * dset       )
00943 {
00944    int fixed_axis , fixed_index , dsl_1 , dsl_2 , rot,mir;
00945    MRI_IMAGE * dsim , * flim ;
00946    THD_dataxes * daxes ;
00947 
00948 ENTRY("AFNI_slice_flip") ;
00949 
00950    /*--- get the image from the dataset ---*/
00951 
00952    if( dset == NULL || kslice < 0 ) RETURN(NULL) ;  /* should not happen! */
00953 
00954    daxes = CURRENT_DAXES(dset) ;
00955 
00956    /* determine which axis kslice refers to,
00957       which determines fixed_index, and also compute
00958       dsl_1 and dsl_2 = axis codes for image returned from AFNI_dataset_slice() */
00959 
00960    switch( ax_3 ){
00961       default: RETURN(NULL) ;
00962 
00963       case  1:  fixed_axis  = 1 ; dsl_1 = 2 ; dsl_2 = 3 ;
00964                 fixed_index = kslice ; break ;
00965 
00966       case -1:  fixed_axis  = 1 ; dsl_1 = 2 ; dsl_2 = 3 ;
00967                 fixed_index = daxes->nxx - kslice - 1 ; break ;
00968 
00969       case  2:  fixed_axis  = 2 ; dsl_1 = 3 ; dsl_2 = 1 ;
00970                 fixed_index = kslice ; break ;
00971 
00972       case -2:  fixed_axis  = 2 ; dsl_1 = 3 ; dsl_2 = 1 ;
00973                 fixed_index = daxes->nyy - kslice - 1 ; break ;
00974 
00975       case  3:  fixed_axis  = 3 ; dsl_1 = 1 ; dsl_2 = 2 ;
00976                 fixed_index = kslice ; break ;
00977 
00978       case -3:  fixed_axis  = 3 ; dsl_1 = 1 ; dsl_2 = 2 ;
00979                 fixed_index = daxes->nzz - kslice - 1 ; break ;
00980    }
00981 
00982    dsim = AFNI_dataset_slice( dset ,
00983                               fixed_axis , fixed_index , ival , resam ) ;
00984 
00985    if( dsim == NULL ) RETURN(NULL) ;
00986 
00987    /*--- the image may need to be flippo-ed to be
00988          in the orientation that the FD_brick wants ---*/
00989 
00990           if( dsl_1 ==  ax_1 && dsl_2 ==  ax_2 ){  /* the easy case */
00991       RETURN(dsim) ;
00992 
00993    } else if( dsl_1 ==  ax_1 && dsl_2 == -ax_2 ){
00994       rot = MRI_ROT_180 ;
00995       mir = True ;
00996 
00997    } else if( dsl_1 == -ax_1 && dsl_2 ==  ax_2 ){
00998       rot = MRI_ROT_0 ;
00999       mir = True ;
01000 
01001    } else if( dsl_1 == -ax_1 && dsl_2 == -ax_2 ){
01002       rot = MRI_ROT_180 ;
01003       mir = False ;
01004 
01005    } else if( dsl_1 ==  ax_2 && dsl_2 == -ax_1 ){
01006       rot = MRI_ROT_270 ;
01007       mir = False ;
01008 
01009    } else if( dsl_1 ==  ax_2 && dsl_2 ==  ax_1 ){
01010       rot = MRI_ROT_270 ;
01011       mir = True ;
01012 
01013    } else if( dsl_1 == -ax_2 && dsl_2 ==  ax_1 ){
01014       rot = MRI_ROT_90 ;
01015       mir = False ;
01016 
01017    } else if( dsl_1 == -ax_2 && dsl_2 == -ax_1 ){
01018       rot = MRI_ROT_90 ;
01019       mir = True ;
01020 
01021    } else {  /* should never happen! */
01022 
01023       fprintf(stderr,"\a\n*** dsl_1=%d dsl_2=%d ax_1=%d ax_2=%d\n",
01024               dsl_1,dsl_2,ax_1,ax_2 ) ;
01025       mri_free( dsim ) ;
01026       RETURN(NULL) ;
01027    }
01028 
01029    flim = mri_flippo( rot , mir , dsim ) ;  /* will set dx,dy OK */
01030 
01031    if( dsim != flim ) mri_free( dsim ) ;
01032 
01033    RETURN(flim) ;
01034 }
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

afni_warp.c
