00001 #include "mrilib.h"
00002 
00003 /*-----------------------------------------------------------------*/
00004 /*! Compute median brick of a dataset - 12 Aug 2001
00005     - 05 Nov 2001: Modified to use THD_extract_array() instead
00006                    of THD_extract_series()
00007 -------------------------------------------------------------------*/
00008 
00009 MRI_IMAGE * THD_median_brick( THD_3dim_dataset *dset )
00010 {
00011    int nvox , nvals , ii ;
00012    MRI_IMAGE *tsim , *medim ;
00013    float *medar ;
00014    float *tsar ;  /* 05 Nov 2001 */
00015 
00016 ENTRY("THD_median_brick") ;
00017 
00018    if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
00019    DSET_load(dset) ;
00020    if( !DSET_LOADED(dset) ) RETURN(NULL) ;
00021 
00022    nvals = DSET_NVALS(dset) ;
00023    tsim  = DSET_BRICK(dset,0) ;
00024 
00025    if( nvals == 1 ){
00026      medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
00027      RETURN(medim) ;
00028    }
00029 
00030    medim = mri_new_conforming( tsim , MRI_float ) ;
00031    medar = MRI_FLOAT_PTR(medim) ;
00032    nvox  = DSET_NVOX(dset) ;
00033 
00034    tsar = (float *) calloc( sizeof(float),nvals+1 ) ; /* 05 Nov 2001 */
00035    for( ii=0 ; ii < nvox ; ii++ ){
00036       THD_extract_array( ii , dset , 0 , tsar ) ;     /* 05 Nov 2001 */
00037       medar[ii] = qmed_float( nvals , tsar ) ;
00038    }
00039 
00040    free(tsar) ; RETURN(medim) ;
00041 }
00042 
00043 /*-----------------------------------------------------------------*/
00044 /*! Compute mean brick of a dataset.  [15 Apr 2005 - RWCox]
00045 -------------------------------------------------------------------*/
00046 
00047 MRI_IMAGE * THD_mean_brick( THD_3dim_dataset *dset )
00048 {
00049    int nvox , nvals , ii , jj ;
00050    MRI_IMAGE *tsim , *medim ;
00051    float *medar , sum,fac ;
00052    float *tsar ;
00053 
00054 ENTRY("THD_mean_brick") ;
00055 
00056    if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
00057    DSET_load(dset) ;
00058    if( !DSET_LOADED(dset) ) RETURN(NULL) ;
00059 
00060    nvals = DSET_NVALS(dset)   ; fac = 1.0 / nvals ;
00061    tsim  = DSET_BRICK(dset,0) ;
00062 
00063    if( nvals == 1 ){
00064      medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
00065      RETURN(medim) ;
00066    }
00067 
00068    medim = mri_new_conforming( tsim , MRI_float ) ;
00069    medar = MRI_FLOAT_PTR(medim) ;
00070    nvox  = DSET_NVOX(dset) ;
00071 
00072    tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
00073    for( ii=0 ; ii < nvox ; ii++ ){
00074       THD_extract_array( ii , dset , 0 , tsar ) ;
00075       for( sum=0.0,jj=0 ; jj < nvals ; jj++ ) sum += tsar[jj] ;
00076       medar[ii] = fac * sum ;
00077    }
00078 
00079    free(tsar) ; RETURN(medim) ;
00080 }
00081 
00082 /*-----------------------------------------------------------------*/
00083 /*! Compute RMS brick of a dataset.  [15 Apr 2005 - RWCox]
00084 -------------------------------------------------------------------*/
00085 
00086 MRI_IMAGE * THD_rms_brick( THD_3dim_dataset *dset )
00087 {
00088    int nvox , nvals , ii , jj ;
00089    MRI_IMAGE *tsim , *medim ;
00090    float *medar , sum,fac ;
00091    float *tsar ;
00092 
00093 ENTRY("THD_rms_brick") ;
00094 
00095    if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
00096    DSET_load(dset) ;
00097    if( !DSET_LOADED(dset) ) RETURN(NULL) ;
00098 
00099    nvals = DSET_NVALS(dset)   ; fac = 1.0 / nvals ;
00100    tsim  = DSET_BRICK(dset,0) ;
00101 
00102    if( nvals == 1 ){
00103      medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
00104      RETURN(medim) ;
00105    }
00106 
00107    medim = mri_new_conforming( tsim , MRI_float ) ;
00108    medar = MRI_FLOAT_PTR(medim) ;
00109    nvox  = DSET_NVOX(dset) ;
00110 
00111    tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
00112    for( ii=0 ; ii < nvox ; ii++ ){
00113       THD_extract_array( ii , dset , 0 , tsar ) ;
00114       for( sum=0.0,jj=0 ; jj < nvals ; jj++ ) sum += tsar[jj]*tsar[jj] ;
00115       medar[ii] = sqrt(fac * sum) ;
00116    }
00117 
00118    free(tsar) ; RETURN(medim) ;
00119 }