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   Routine to extract a time-series (fixed index, variable ival)
00012     from a dataset.
00013   ixyz = spatial index of desired voxel
00014        = ix + jy * n1 + kz * n1*n2
00015   raw  = 0 if you always want floats
00016        = 1 if you want the truly stored data type
00017 
00018   05 Nov 2001: Split into 2 functions -
00019       THD_extract_series() produces a new image each time
00020       THD_extract_array()  copies data into a user-supplied array
00021 -----------------------------------------------------------------*/
00022 
00023 MRI_IMAGE * THD_extract_series( int ind , THD_3dim_dataset *dset , int raw )
00024 {
00025    int nv , typ , ii ;
00026    MRI_IMAGE *im ;
00027    void *imar ;
00028 
00029 ENTRY("THD_extract_series") ;
00030 
00031    if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
00032 
00033    nv  = dset->dblk->nvals ;
00034    if( raw ) typ = DSET_BRICK_TYPE(dset,0) ;  /* type of output array */
00035    else      typ = MRI_float ;
00036 
00037    im   = mri_new( nv , 1 , typ ) ;           /* output image */
00038    imar = mri_data_pointer(im) ;
00039 
00040    ii = THD_extract_array( ind , dset , raw , imar ) ; /* get data */
00041 
00042    if( ii != 0 ){ mri_free(im) ; RETURN(NULL) ; }      /* bad */
00043 
00044    if( dset->taxis != NULL ){  /* 21 Oct 1996 */
00045       float zz , tt ;
00046       int kz = ind / ( dset->daxes->nxx * dset->daxes->nyy ) ;
00047 
00048       zz = dset->daxes->zzorg + kz * dset->daxes->zzdel ;
00049       tt = THD_timeof( 0 , zz , dset->taxis ) ;
00050 
00051       im->xo = tt ; im->dx = dset->taxis->ttdel ;   /* origin and delta */
00052 
00053       if( dset->taxis->units_type == UNITS_MSEC_TYPE ){ /* convert to sec */
00054          im->xo *= 0.001 ; im->dx *= 0.001 ;
00055       }
00056    } else {
00057       im->xo = 0.0 ; im->dx = 1.0 ;  /* 08 Nov 1996 */
00058    }
00059 
00060    RETURN(im) ;
00061 }
00062 
00063 /*---------------------------------------------------------------------------
00064   Return value is 0 for all is good, -1 for all is bad.
00065   If raw == 0, uar is of type float.
00066   If raw != 0, uar is of the same type as the dataset brick
00067     (the user need to know the type ahead of time if raw != 0).
00068   Data goes into a user-supplied array.
00069 -----------------------------------------------------------------------------*/
00070 
00071 int THD_extract_array( int ind, THD_3dim_dataset *dset, int raw, void *uar )
00072 {
00073    MRI_TYPE typ ;
00074    int nv , ival , nb ;
00075    char  *iar ;      /* brick in the input */
00076    float *far=NULL ; /* non-raw output */
00077    static void *tar=NULL ; static int ntar=0 ;
00078 
00079 ENTRY("THD_extract_array") ;
00080 
00081    if( ind < 0             || uar == NULL           ||
00082        !ISVALID_DSET(dset) || ind >= DSET_NVOX(dset)  ) RETURN(-1) ;
00083 
00084    nv  = dset->dblk->nvals ;
00085    iar = DSET_ARRAY(dset,0) ;
00086    if( iar == NULL ){         /* load data from disk? */
00087      DSET_load(dset) ;
00088      iar = DSET_ARRAY(dset,0); if( iar == NULL ) RETURN(-1) ;
00089    }
00090    typ = DSET_BRICK_TYPE(dset,0) ;  /* raw data type */
00091 
00092    /* will extract nb bytes of raw data into array tar */
00093 
00094    nb = mri_datum_size(typ) * (nv+1) ;
00095    if( nb > ntar ){ tar = AFREALL(tar,void *,nb) ; ntar = nb ; }
00096    memset(tar,0,nb) ;
00097 
00098    if( !raw ) far = (float *) uar ;  /* non-raw output */
00099 
00100    switch( typ ){
00101 
00102       default:           /* don't know what to do --> return nada */
00103          RETURN(-1);
00104       break ;
00105 
00106       case MRI_byte:{
00107          byte *ar = (byte *)tar , *bar ;
00108          for( ival=0 ; ival < nv ; ival++ ){
00109             bar = (byte *) DSET_ARRAY(dset,ival) ;
00110             if( bar != NULL ) ar[ival] = bar[ind] ;
00111          }
00112          if( !raw ){
00113             for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
00114          }
00115       }
00116       break ;
00117 
00118       case MRI_short:{
00119          short *ar = (short *)tar , *bar ;
00120          for( ival=0 ; ival < nv ; ival++ ){
00121             bar = (short *) DSET_ARRAY(dset,ival) ;
00122             if( bar != NULL ) ar[ival] = bar[ind] ;
00123          }
00124          if( !raw ){
00125             for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
00126          }
00127       }
00128       break ;
00129 
00130       case MRI_float:{
00131          float *ar = (float *)tar , *bar ;
00132          for( ival=0 ; ival < nv ; ival++ ){
00133             bar = (float *) DSET_ARRAY(dset,ival) ;
00134             if( bar != NULL ) ar[ival] = bar[ind] ;
00135          }
00136          if( !raw ){
00137             for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
00138          }
00139       }
00140       break ;
00141 
00142       case MRI_int:{
00143          int *ar = (int *)tar , *bar ;
00144          for( ival=0 ; ival < nv ; ival++ ){
00145             bar = (int *) DSET_ARRAY(dset,ival) ;
00146             if( bar != NULL ) ar[ival] = bar[ind] ;
00147          }
00148          if( !raw ){
00149             for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
00150          }
00151       }
00152       break ;
00153 
00154       case MRI_double:{
00155          double *ar = (double *)tar , *bar ;
00156          for( ival=0 ; ival < nv ; ival++ ){
00157             bar = (double *) DSET_ARRAY(dset,ival) ;
00158             if( bar != NULL ) ar[ival] = bar[ind] ;
00159          }
00160          if( !raw ){
00161             for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
00162          }
00163       }
00164       break ;
00165 
00166       case MRI_complex:{
00167          complex *ar = (complex *)tar , *bar ;
00168          for( ival=0 ; ival < nv ; ival++ ){
00169             bar = (complex *) DSET_ARRAY(dset,ival) ;
00170             if( bar != NULL ) ar[ival] = bar[ind] ;
00171          }
00172          if( !raw ){
00173             for( ival=0 ; ival < nv ; ival++ ) far[ival] = CABS(ar[ival]) ;
00174          }
00175       }
00176       break ;
00177 
00178    }
00179 
00180    if( raw ){ memcpy(uar,tar,nb); RETURN(0); }
00181 
00182    if( THD_need_brick_factor(dset) ){
00183       for( ival=0 ; ival < nv ; ival++ )
00184          if( DSET_BRICK_FACTOR(dset,ival) > 0.0 )
00185            far[ival] *= DSET_BRICK_FACTOR(dset,ival) ;
00186    }
00187 
00188    RETURN(0);
00189 }
00190 
00191 /*----------------------------------------------------------------------------
00192    04 Feb 2000: do a bunch of timeseries at once (for efficiency)
00193    27 Feb 2003: rearranged slightly for more efficiency
00194 ------------------------------------------------------------------------------*/
00195 
00196 MRI_IMARR * THD_extract_many_series( int ns, int *ind, THD_3dim_dataset *dset )
00197 {
00198    MRI_IMARR *imar ; /* output */
00199    MRI_IMAGE *im ;
00200    int nv , ival , kk ;
00201    char *iar ;      /* brick in the input */
00202    float **far ;    /* 27 Feb 2003: ptrs to output */
00203 
00204 ENTRY("THD_extract_many_series") ;
00205 
00206    if( ns <= 0 || ind == NULL | dset == NULL ) RETURN( NULL );
00207 
00208    /* try to load dataset */
00209 
00210    nv  = dset->dblk->nvals ;
00211    iar = DSET_ARRAY(dset,0) ;
00212    if( iar == NULL ){  /* if data needs to be loaded from disk */
00213      (void) THD_load_datablock( dset->dblk ) ;
00214      iar = DSET_ARRAY(dset,0) ;
00215      if( iar == NULL ) RETURN( NULL );
00216    }
00217 
00218    /* create output */
00219 
00220    far = (float **) malloc(sizeof(float *)*ns) ;  /* 27 Feb 2003 */
00221    INIT_IMARR(imar) ;
00222    for( kk=0 ; kk < ns ; kk++ ){
00223      im = mri_new( nv , 1 , MRI_float ) ;  /* N.B.: now does 0 fill */
00224      far[kk] = MRI_FLOAT_PTR(im) ;         /* ptr to kk-th output series */
00225      ADDTO_IMARR(imar,im) ;
00226    }
00227 
00228    /* fill the output */
00229 
00230    switch( DSET_BRICK_TYPE(dset,0) ){
00231 
00232       default:             /* don't know what to do --> return nada */
00233         DESTROY_IMARR(imar) ; free(far) ;
00234         RETURN( NULL );
00235 
00236       case MRI_byte:{
00237         byte * bar ;
00238         for( ival=0 ; ival < nv ; ival++ ){
00239           bar = (byte *) DSET_ARRAY(dset,ival) ;
00240           if( bar != NULL ){
00241             for( kk=0 ; kk < ns ; kk++ ){
00242               far[kk][ival] = (float)bar[ind[kk]] ;
00243             }
00244           }
00245         }
00246       }
00247       break ;
00248 
00249       case MRI_short:{
00250         short * bar ;
00251         for( ival=0 ; ival < nv ; ival++ ){
00252           bar = (short *) DSET_ARRAY(dset,ival) ;
00253           if( bar != NULL ){
00254             for( kk=0 ; kk < ns ; kk++ ){
00255               far[kk][ival] = (float)bar[ind[kk]] ;
00256             }
00257           }
00258         }
00259       }
00260       break ;
00261 
00262       case MRI_float:{
00263          float * bar ;
00264          for( ival=0 ; ival < nv ; ival++ ){
00265             bar = (float *) DSET_ARRAY(dset,ival) ;
00266             if( bar != NULL ){
00267               for( kk=0 ; kk < ns ; kk++ ){
00268                  far[kk][ival] = bar[ind[kk]] ;
00269               }
00270             }
00271          }
00272       }
00273       break ;
00274 
00275       case MRI_int:{
00276          int * bar ;
00277          for( ival=0 ; ival < nv ; ival++ ){
00278             bar = (int *) DSET_ARRAY(dset,ival) ;
00279             if( bar != NULL ){
00280               for( kk=0 ; kk < ns ; kk++ ){
00281                  far[kk][ival] = bar[ind[kk]] ;
00282               }
00283             }
00284          }
00285       }
00286       break ;
00287 
00288       case MRI_double:{
00289          double * bar ;
00290          for( ival=0 ; ival < nv ; ival++ ){
00291             bar = (double *) DSET_ARRAY(dset,ival) ;
00292             if( bar != NULL ){
00293               for( kk=0 ; kk < ns ; kk++ ){
00294                  far[kk][ival] = (float)bar[ind[kk]] ;
00295               }
00296             }
00297          }
00298       }
00299       break ;
00300 
00301       case MRI_complex:{
00302          complex * bar ;
00303          for( ival=0 ; ival < nv ; ival++ ){
00304             bar = (complex *) DSET_ARRAY(dset,ival) ;
00305             if( bar != NULL ){
00306               for( kk=0 ; kk < ns ; kk++ ){
00307                  far[kk][ival] = bar[ind[kk]].r ;
00308               }
00309             }
00310          }
00311       }
00312       break ;
00313 
00314    }
00315 
00316    /* scale outputs, if needed */
00317 
00318    if( THD_need_brick_factor(dset) ){
00319       MRI_IMAGE * qim ;
00320       for( kk=0 ; kk < ns ; kk++ ){
00321          im  = IMARR_SUBIMAGE(imar,kk) ;
00322          qim = mri_mult_to_float( dset->dblk->brick_fac , im ) ;
00323          mri_free(im) ;
00324          IMARR_SUBIMAGE(imar,kk) = qim ;
00325       }
00326    }
00327 
00328 #if 0  /* 27 Feb 2003 */
00329    /* convert to floats, if needed */
00330 
00331    if( IMARR_SUBIMAGE(imar,0)->kind != MRI_float ){
00332       MRI_IMAGE * qim ;
00333       for( kk=0 ; kk < ns ; kk++ ){
00334          im  = IMARR_SUBIMAGE(imar,kk) ;
00335          qim = mri_to_float( im ) ;
00336          mri_free(im) ;
00337          IMARR_SUBIMAGE(imar,kk) = qim ;
00338       }
00339    }
00340 #endif
00341 
00342    /* add time axis stuff to output images, if present */
00343 
00344    if( dset->taxis != NULL ){
00345       float zz , tt ;
00346       int kz ;
00347 
00348       for( kk=0 ; kk < ns ; kk++ ){
00349          kz = ind[kk] / ( dset->daxes->nxx * dset->daxes->nyy ) ;
00350          zz = dset->daxes->zzorg + kz * dset->daxes->zzdel ;
00351          tt = THD_timeof( 0 , zz , dset->taxis ) ;
00352          im = IMARR_SUBIMAGE(imar,kk) ;
00353          im->xo = tt ; im->dx = dset->taxis->ttdel ;   /* origin and delta */
00354          if( dset->taxis->units_type == UNITS_MSEC_TYPE ){ /* convert to sec */
00355             im->xo *= 0.001 ; im->dx *= 0.001 ;
00356          }
00357       }
00358    } else {
00359       for( kk=0 ; kk < ns ; kk++ ){
00360          im = IMARR_SUBIMAGE(imar,kk) ;
00361          im->xo = 0.0 ; im->dx = 1.0 ;
00362       }
00363    }
00364 
00365    free(far) ; RETURN(imar);
00366 }