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 
00009 /*** NOT 7D SAFE ***/
00010 
00011 /* Collect stats on sequence of images.
00012    Usage: call repeatedly, once for each image;  will return void.
00013           after last image, call with NULL as argument, will return
00014           pointer to array of two pointers to (float) images, the
00015           first the mean, the second the stdev.
00016 */
00017 
00018 MRI_IMAGE ** mri_stat_seq( MRI_IMAGE * imin )
00019 {
00020 
00021 /* static variables (exist between calls) */
00022 
00023    static MRI_IMAGE * imsum , * imsumq ;
00024    static int nim = 0 , npix , nx , ny ;
00025    static double * dbsum , * dbsumq ;
00026 
00027 /* local variables */
00028 
00029    register int ii ;
00030    MRI_IMAGE * imfl , * imsd , ** retval ;
00031    float     * flar , * sdar ;
00032    double    scl , vscl , vvv ;
00033 
00034 /*** case: set up new problem ***/
00035 
00036    if( nim == 0 ){
00037 
00038       if( imin == NULL ){
00039          fprintf(stderr,"mri_stat_seq:  NULL argument for initial call!\n") ;
00040          EXIT(1) ;
00041       }
00042 
00043       if( ! MRI_IS_2D(imin) ){
00044          fprintf(stderr,"\n*** mri_stat_seq: only works on 2D images!\n") ;
00045          EXIT(1) ;
00046       }
00047 
00048       nx   = imin->nx ;
00049       ny   = imin->ny ;
00050       npix = nx * ny ;
00051 
00052       imsum  = mri_new( nx , ny , MRI_double ) ;
00053       imsumq = mri_new( nx , ny , MRI_double ) ;
00054       dbsum  = mri_data_pointer( imsum ) ;
00055       dbsumq = mri_data_pointer( imsumq ) ;
00056 
00057       MRI_COPY_AUX(imsum,imin) ;
00058 
00059       for( ii=0 ; ii < npix ; ii++ ) dbsum[ii] = dbsumq[ii] = 0.0 ;
00060 
00061    }
00062 
00063 /*** case: add new image into problem ***/
00064 
00065    if( imin != NULL ){
00066 
00067       if( imin->nx != nx || imin->ny != ny ){
00068          fprintf(stderr,"mri_stat_seq: input image size mismatch!\n") ;
00069          EXIT(1) ;
00070       }
00071 
00072       if( ! MRI_IS_2D(imin) ){
00073          fprintf(stderr,"\n*** mri_stat_seq: only works on 2D images!\n") ;
00074          EXIT(1) ;
00075       }
00076 
00077       if( imin->kind == MRI_float ){
00078          imfl = imin ;
00079       } else {
00080          imfl = mri_to_float( imin ) ;
00081       }
00082       flar = mri_data_pointer( imfl ) ;
00083 
00084       for( ii=0 ; ii < npix ; ii++ ){
00085          dbsum[ii]  += flar[ii] ;
00086          dbsumq[ii] += flar[ii] * flar[ii] ;
00087       }
00088 
00089       if( imin != imfl ) mri_free( imfl ) ;
00090       nim++ ;
00091       return NULL ;
00092    }
00093 
00094 /*** case: make report on results, and reset static data ***/
00095 
00096    if( nim < 1 ){
00097       fprintf(stderr,"mri_stat_seq: # images input=%d; too small!\n",nim) ;
00098       EXIT(1) ;
00099    }
00100 
00101    imfl = mri_new( nx , ny , MRI_float ) ;
00102    imsd = mri_new( nx , ny , MRI_float ) ;
00103    flar = mri_data_pointer( imfl ) ;
00104    sdar = mri_data_pointer( imsd ) ;
00105 
00106    MRI_COPY_AUX(imfl,imsum) ;
00107 
00108    scl  = 1.0 / nim ;
00109    vscl = (nim==1) ? 1.0 : sqrt( ((double)(nim))/(nim-1.0) ) ;
00110 
00111    for( ii=0 ; ii < npix ; ii++ ){
00112       flar[ii] = scl * dbsum[ii] ;
00113       vvv      = scl * dbsumq[ii] - flar[ii] * flar[ii] ;
00114       sdar[ii] = (vvv > 0.0) ? (sqrt(vvv)*vscl) : 0.0 ;
00115    }
00116 
00117    mri_free( imsum ) ; mri_free( imsumq ) ;
00118    nim = 0 ;
00119 
00120    retval = (MRI_IMAGE **) malloc( 2 * sizeof(MRI_IMAGE *) ) ;
00121    if( retval == NULL ){
00122       fprintf(stderr,"mri_stat_seq: malloc error for retval!\n") ;
00123       EXIT(1) ;
00124    }
00125 
00126    retval[0] = imfl ;
00127    retval[1] = imsd ;
00128 
00129    return retval ;
00130 }