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 previously set up FD_brick structure.
00013   ixyz = spatial index of desired voxel (in brick coordinates)
00014        = ix + jy * n1 + kz * n1*n2
00015   Return value is an image of the type of the dataset (assumed
00016   uniform).
00017 -----------------------------------------------------------------*/
00018 
00019 MRI_IMAGE * FD_brick_to_series( int ixyz , FD_brick * br )
00020 {
00021    MRI_IMAGE * im ;  /* output */
00022    int nv , ival ;
00023    char * iar ;      /* brick in the input */
00024    MRI_TYPE typ ;
00025    int ix,jy,kz , ind ;
00026    THD_ivec3 ind_fd , ind_ds ;
00027 
00028    if( ixyz < 0 || ixyz >= br->n1 * br->n2 * br->n3 ) return NULL ;
00029 
00030    /** otherwise, get ready for a real image **/
00031 
00032    ix  = ixyz % br->n1 ;
00033    jy  = ( ixyz % (br->n1 * br->n2) ) / br->n1 ;
00034    kz  = ixyz / (br->n1 * br->n2) ;
00035    LOAD_IVEC3( ind_fd , ix,jy,kz ) ; ind_ds = THD_fdind_to_3dind( br , ind_fd ) ;
00036    ix  = ind_ds.ijk[0] ;
00037    jy  = ind_ds.ijk[1] ;
00038    kz  = ind_ds.ijk[2] ;
00039    ind = (kz * br->dset->daxes->nyy + jy) * br->dset->daxes->nxx + ix ;
00040 
00041    nv = br->dset->dblk->nvals ;
00042 
00043    iar = DSET_ARRAY(br->dset,0) ;
00044    if( iar == NULL ){  /* if data needs to be loaded from disk */
00045       (void) THD_load_datablock( br->dset->dblk ) ;
00046       iar = DSET_ARRAY(br->dset,0) ;
00047       if( iar == NULL ) return NULL ;
00048    }
00049 
00050    /* 15 Sep 2004: allow for nonconstant datum */
00051 
00052    if( !DSET_datum_constant(br->dset) ){  /* only for stupid users */
00053      float *ar ;
00054      im = mri_new( nv , 1 , MRI_float ) ; ar = MRI_FLOAT_PTR(im) ;
00055      for( ival = 0 ; ival < nv ; ival++ )
00056        ar[ival] = THD_get_voxel( br->dset , ind , ival ) ;
00057      goto image_done ;
00058    }
00059 
00060    /* the older (more efficient) way */
00061 
00062    typ = DSET_BRICK_TYPE(br->dset,0) ;
00063    im  = mri_new( nv , 1 , typ ) ;
00064 #if 0
00065    mri_zero_image(im) ;             /* 18 Oct 2001 */
00066 #endif
00067 
00068    switch( typ ){
00069 
00070       default:             /* don't know what to do --> return nada */
00071          mri_free( im ) ;
00072          return NULL ;
00073 
00074       case MRI_byte:{
00075          byte * ar  = MRI_BYTE_PTR(im) , * bar ;
00076          for( ival=0 ; ival < nv ; ival++ ){
00077             bar = (byte *) DSET_ARRAY(br->dset,ival) ;
00078             if( bar != NULL ) ar[ival] = bar[ind] ;
00079          }
00080       }
00081       break ;
00082 
00083       case MRI_short:{
00084          short * ar  = MRI_SHORT_PTR(im) , * bar ;
00085          for( ival=0 ; ival < nv ; ival++ ){
00086             bar = (short *) DSET_ARRAY(br->dset,ival) ;
00087             if( bar != NULL ) ar[ival] = bar[ind] ;
00088          }
00089       }
00090       break ;
00091 
00092       case MRI_float:{
00093          float * ar  = MRI_FLOAT_PTR(im) , * bar ;
00094          for( ival=0 ; ival < nv ; ival++ ){
00095             bar = (float *) DSET_ARRAY(br->dset,ival) ;
00096             if( bar != NULL ) ar[ival] = bar[ind] ;
00097          }
00098       }
00099       break ;
00100 
00101       case MRI_int:{
00102          int * ar  = MRI_INT_PTR(im) , * bar ;
00103          for( ival=0 ; ival < nv ; ival++ ){
00104             bar = (int *) DSET_ARRAY(br->dset,ival) ;
00105             if( bar != NULL ) ar[ival] = bar[ind] ;
00106          }
00107       }
00108       break ;
00109 
00110       case MRI_double:{
00111          double * ar  = MRI_DOUBLE_PTR(im) , * bar ;
00112          for( ival=0 ; ival < nv ; ival++ ){
00113             bar = (double *) DSET_ARRAY(br->dset,ival) ;
00114             if( bar != NULL ) ar[ival] = bar[ind] ;
00115          }
00116       }
00117       break ;
00118 
00119       case MRI_complex:{
00120          complex * ar  = MRI_COMPLEX_PTR(im) , * bar ;
00121          for( ival=0 ; ival < nv ; ival++ ){
00122             bar = (complex *) DSET_ARRAY(br->dset,ival) ;
00123             if( bar != NULL ) ar[ival] = bar[ind] ;
00124          }
00125       }
00126       break ;
00127 
00128       /* 15 Apr 2002: RGB types */
00129 
00130       case MRI_rgb:{
00131          rgbyte *ar  = (rgbyte *) MRI_RGB_PTR(im) , *bar ;
00132          for( ival=0 ; ival < nv ; ival++ ){
00133             bar = (rgbyte *) DSET_ARRAY(br->dset,ival) ;
00134             if( bar != NULL ) ar[ival] = bar[ind] ;
00135          }
00136       }
00137       break ;
00138 
00139       case MRI_rgba:{
00140          rgba *ar  = (rgba *) MRI_RGBA_PTR(im) , *bar ;
00141          for( ival=0 ; ival < nv ; ival++ ){
00142             bar = (rgba *) DSET_ARRAY(br->dset,ival) ;
00143             if( bar != NULL ) ar[ival] = bar[ind] ;
00144          }
00145       }
00146       break ;
00147 
00148    }
00149 
00150    if( THD_need_brick_factor(br->dset) ){
00151       MRI_IMAGE * qim ;
00152       qim = mri_mult_to_float( br->dset->dblk->brick_fac , im ) ;
00153       mri_free(im) ; im = qim ;
00154    }
00155 
00156    /* at this point, the image is ready to ship out;
00157       but first, maybe attach a time origin and spacing */
00158 
00159 image_done:
00160    if( br->dset->taxis != NULL ){  /* 21 Oct 1996 */
00161       float zz , tt ;
00162 
00163       zz = br->dset->daxes->zzorg + kz * br->dset->daxes->zzdel ;
00164       tt = THD_timeof( 0 , zz , br->dset->taxis ) ;
00165 
00166       im->xo = tt ; im->dx = br->dset->taxis->ttdel ;   /* origin and delta */
00167 
00168       if( br->dset->taxis->units_type == UNITS_MSEC_TYPE ){ /* convert to sec */
00169          im->xo *= 0.001 ; im->dx *= 0.001 ;
00170       }
00171    } else {
00172       im->xo = 0.0 ; im->dx = 1.0 ;  /* 08 Nov 1996 */
00173    }
00174 
00175    return im ;
00176 }