Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search

rcov.c File Reference

#include "mrilib.h"


Data Structures
struct	covmat
Defines
#define	CM(i, j) cmat[(i)+(j)*ndim]
#define	CH(i, j) cfac[(i)+(j)*ndim]
#define	CCUT 3.5
#define	EPS 1.e-4
#define	DDD csum
Functions
void	forward_solve_inplace (covmat cv, float vec)
void	backward_solve_inplace (covmat cv, float vec)
void	compute_choleski (covmat *cv)
covmat *	robust_covar (int ndim, int nvec, float **vec)
int	main (int argc, char *argv[])

Define Documentation

#define CCUT 3.5

Definition at line 72 of file rcov.c.

Referenced by robust_covar().

#define CH	(	i,
		j		)	cfac[(i)+(j)*ndim]

Definition at line 9 of file rcov.c.

#define CM	(	i,
		j		)	cmat[(i)+(j)*ndim]

Definition at line 8 of file rcov.c.

#define DDD csum

#define EPS 1.e-4

Definition at line 73 of file rcov.c.

Referenced by robust_covar().

Function Documentation

void backward_solve_inplace	(	covmat *	cv,
		float *	vec
	)

Definition at line 28 of file rcov.c.

References covmat::cfac, CH, covmat::ndim, and vec.

Referenced by evaluate_span().

00029 {
00030    register int     ndim=cv->ndim , ii,jj ;
00031    register float * cfac=cv->cfac , sum ;
00032 
00033    for( ii=ndim-1 ; ii >= 0 ; ii-- ){
00034       sum = vec[ii] ;
00035       for( jj=ii+1 ; jj < ndim ; jj++ ) sum -= CH(jj,ii) * vec[jj] ;
00036       vec[ii] = sum / CH(ii,ii) ;
00037    }
00038    return ;
00039 }

void compute_choleski

(

covmat *

)

Definition at line 43 of file rcov.c.

References covmat::cfac, CH, CM, covmat::cmat, free, malloc, and covmat::ndim.

Referenced by evaluate_span(), and robust_covar().

00044 {
00045    register int     ndim=cv->ndim ,          ii,jj,kk ;
00046    register float * cmat=cv->cmat , * cfac , sum ;
00047 
00048    if( ndim < 2 || cmat == NULL ) return ;
00049 
00050    if( cv->cfac == NULL )
00051       cv->cfac = (float *) malloc(sizeof(float)*ndim*ndim) ;
00052 
00053    cfac = cv->cfac ;
00054 
00055    for( ii=0 ; ii < ndim ; ii++ ){
00056       for( jj=0 ; jj < ii ; jj++ ){
00057          sum = CM(ii,jj) ;
00058          for( kk=0 ; kk < jj ; kk++ ) sum -= CH(ii,kk) * CH(jj,kk) ;
00059          CH(ii,jj) = sum / CH(jj,jj) ;
00060       }
00061       sum = CM(ii,ii) ;
00062       for( kk=0 ; kk < ii ; kk++ ) sum -= CH(ii,kk) * CH(ii,kk) ;
00063       if( sum <= 0.0 ){ free(cv->cfac); cv->cfac = NULL; return; }
00064       CH(ii,ii) = sqrt(sum) ;
00065       for( jj=ii+1 ; jj < ndim ; jj++ ) CH(ii,jj) = 0.0 ;
00066    }
00067    return ;
00068 }

void forward_solve_inplace	(	covmat *	cv,
		float *	vec
	)

Definition at line 13 of file rcov.c.

References covmat::cfac, CH, covmat::ndim, and vec.

Referenced by evaluate_span(), and robust_covar().

00014 {
00015    register int     ndim=cv->ndim , ii,jj ;
00016    register float * cfac=cv->cfac , sum ;
00017 
00018    for( ii=0 ; ii < ndim ; ii++ ){
00019       sum = vec[ii] ;
00020       for( jj=0 ; jj < ii ; jj++ ) sum -= CH(ii,jj) * vec[jj] ;
00021       vec[ii] = sum / CH(ii,ii) ;
00022    }
00023    return ;
00024 }

int main	(	int	argc,
		char *	argv[]
	)

\** File : SUMA.c

Author:: : Ziad Saad Date : Thu Dec 27 16:21:01 EST 2001

Purpose :

Input paramters :

Parameters:

param

Usage : SUMA ( )

Returns :

Returns:: Support :

See also:: OpenGL prog. Guide 3rd edition , varray.c from book's sample code

Side effects :

Definition at line 228 of file rcov.c.

References argc, far, MRI_IMAGE::kind, malloc, MRI_FLOAT_PTR, mri_free(), mri_read_just_one(), mri_to_float(), mri_transpose(), MRI_IMAGE::nx, MRI_IMAGE::ny, robust_covar(), and vec.

00229 {
00230    MRI_IMAGE * im , * qim ;
00231    float * far , ** vec ;
00232    int ii,jj , ndim,nvec ;
00233 
00234    if( argc < 2 ){printf("Usage: rcov imfile\n"); exit(0); }
00235 
00236    im = mri_read_just_one( argv[1] ) ;
00237    if( im == NULL ) exit(1) ;
00238    if( im->kind != MRI_float ){
00239       qim = mri_to_float(im) ; mri_free(im) ; im = qim ;
00240    }
00241    qim = mri_transpose(im) ; mri_free(im) ; im = qim ;
00242 
00243    ndim = im->nx ; nvec = im->ny ;
00244    far = MRI_FLOAT_PTR(im) ;
00245    vec = (float **) malloc(sizeof(float *)*nvec) ;
00246    for( jj=0 ; jj < nvec ; jj++ ) vec[jj] = far + jj*ndim ;
00247 
00248 #if 0
00249    far = (float *) malloc(sizeof(float)*ndim) ;
00250    for( ii=0 ; ii < ndim ; ii++ ) far[ii] = 0.0 ;
00251    for( jj=0 ; jj < nvec ; jj++ )
00252       for( ii=0 ; ii < ndim ; ii++ ) far[ii] += vec[jj][ii] ;
00253    for( ii=0 ; ii < ndim ; ii++ ) far[ii] /= nvec ;
00254    for( jj=0 ; jj < nvec ; jj++ )
00255       for( ii=0 ; ii < ndim ; ii++ ) vec[jj][ii] -= far[ii] ;
00256 #endif
00257 
00258    (void) robust_covar( ndim , nvec , vec ) ;
00259    exit(0) ;
00260 }

covmat* robust_covar	(	int	ndim,
		int	nvec,
		float **	vec
	)

Definition at line 75 of file rcov.c.

References CCUT, covmat::cfac, covmat::cmat, compute_choleski(), EPS, forward_solve_inplace(), free, malloc, covmat::mvec, covmat::ndim, and vec.

Referenced by evaluate_span(), and main().

00076 {
00077    covmat * cv ;
00078    float *nmat, *cmat , fnvec,fndim,cnorm,csum , *tv , *vv , *mv , *wv ;
00079    int ii , jj , kk , nite ;
00080    float bcut , cwt ;
00081 
00082 fprintf(stderr,"Enter robust_covar:  ndim=%d  nvec=%d\n",ndim,nvec) ;
00083 
00084    if( ndim < 2 || nvec < ndim || vec == NULL ) return NULL ;
00085 
00086    cv = (covmat *) malloc(sizeof(covmat)) ;
00087    cv->ndim = ndim ;
00088    cv->cmat = NULL ;
00089    cv->cfac = NULL ;
00090    cv->mvec = NULL ;
00091 
00092    nmat = (float *) malloc(sizeof(float)*ndim*ndim) ;  /* matrix     */
00093    tv   = (float *) malloc(sizeof(float)*ndim) ;       /* temp vector */
00094    mv   = (float *) malloc(sizeof(float)*ndim) ;       /* mean vector  */
00095    wv   = (float *) malloc(sizeof(float)*nvec) ;       /* weight vector */
00096 
00097    fnvec = 1.0/nvec ; fndim = 1.0/ndim ;
00098    bcut  = 1.0 + CCUT*sqrt(fndim) ;
00099 
00100    /* compute initial mean & covariance matrix with all weights = 1 */
00101 
00102    for( jj=0 ; jj < ndim ; jj++ ) mv[jj] = 0.0 ;
00103 
00104    for( kk=0 ; kk < nvec ; kk++ ){   /* mean vector sum */
00105       vv = vec[kk] ;
00106       for( jj=0 ; jj < ndim ; jj++ ) mv[jj] += vv[jj] ;
00107    }
00108    for( jj=0 ; jj < ndim ; jj++ ) mv[jj] *= fnvec ;  /* scale mean vector */
00109 
00110    for( jj=0 ; jj < ndim ; jj++ )
00111       for( ii=0 ; ii < ndim ; ii++ ) nmat[ii+jj*ndim] = 0.0 ;
00112 
00113    for( kk=0 ; kk < nvec ; kk++ ){   /* covariance matrix sum */
00114       vv = vec[kk] ;
00115       for( jj=0 ; jj < ndim ; jj++ ){
00116          for( ii=0 ; ii <= jj ; ii++ )
00117             nmat[ii+jj*ndim] += (vv[ii]-mv[ii])*(vv[jj]-mv[jj]) ;
00118       }
00119    }
00120    for( jj=0 ; jj < ndim ; jj++ ){   /* scale covariance matrix */
00121       for( ii=0 ; ii < jj ; ii++ )
00122          nmat[jj+ii*ndim] = (nmat[ii+jj*ndim] *= fnvec) ;
00123       nmat[jj+jj*ndim] *= fnvec ;
00124    }
00125 
00126    /* now iterate */
00127 
00128    nite = 0 ;
00129 
00130    while(1){
00131 
00132       nite++ ;
00133 
00134       cmat = cv->cmat = nmat ;  /* put old matrix into cv */
00135       compute_choleski(cv) ;    /* decompose it */
00136 
00137       nmat = (float *) malloc(sizeof(float)*ndim*ndim) ; /* new matrix */
00138       cv->mvec = mv ;                                    /* old mean vector */
00139       mv = (float *) malloc(sizeof(float)*ndim) ;        /* new mean vector */
00140 
00141       for( jj=0 ; jj < ndim ; jj++ ){  /* initialize new things to zero */
00142          mv[jj] = 0.0 ;
00143          for( ii=0 ; ii < ndim ; ii++ ) nmat[ii+jj*ndim] = 0.0 ;
00144       }
00145 
00146 fprintf(stderr,"\niteration %2d:\n",nite) ;
00147 
00148       /* update mean */
00149 
00150       csum = 0.0 ;
00151       for( kk=0 ; kk < nvec ; kk++ ){
00152          vv = vec[kk] ;
00153 
00154          /*                    -1/2          */
00155          /* compute tv = [cmat]    (vv-mvec) */
00156 
00157          for( jj=0 ; jj < ndim ; jj++ ) tv[jj] = vv[jj] - cv->mvec[jj] ;
00158          forward_solve_inplace(cv,tv) ;
00159 
00160          /* compute norm of tv, then weighting factor for this vector */
00161 
00162          cnorm = 0.0 ; for( ii=0 ; ii < ndim ; ii++ ) cnorm += tv[ii]*tv[ii] ;
00163          cnorm = cnorm*fndim ;
00164          cnorm = (cnorm <= bcut) ? 1.0 : bcut/cnorm ;
00165          wv[kk] = cnorm ; csum += cnorm ;
00166 #if 0
00167 fprintf(stderr,"  weight %2d = %12.4g\n",kk,cnorm) ;
00168 #endif
00169 #if 0
00170 for( jj=0 ; jj < ndim ; jj++ ) fprintf(stderr," %f:%f:%f",vv[jj],tv[jj],cv->mvec[jj]) ;
00171 fprintf(stderr,"\n") ;
00172 #endif
00173 
00174          /* add vv into accumulating mean, with weight cnorm */
00175 
00176          for( jj=0 ; jj < ndim ; jj++ ) mv[jj] += cnorm*vv[jj] ;
00177       }
00178 #if 0
00179 fprintf(stderr,"  wv:");for(kk=0;kk<nvec;kk++)fprintf(stderr," %11.3g",wv[kk]);fprintf(stderr,"\n csum: %11.3g\n",csum);
00180 #endif
00181       csum = 1.0 / csum ; cwt = nvec*csum ;
00182       for( jj=0 ; jj < ndim ; jj++ ) mv[jj] *= csum ;  /* scale new mean */
00183 
00184       /* update covariance */
00185 
00186       for( kk=0 ; kk < nvec ; kk++ ){
00187          vv = vec[kk] ; cnorm = wv[kk] ;
00188          for( jj=0 ; jj < ndim ; jj++ ){
00189             for( ii=0 ; ii <= jj ; ii++ )
00190                nmat[ii+jj*ndim] +=
00191                   cnorm*(vv[ii]-cv->mvec[ii])*(vv[jj]-cv->mvec[jj]) ;
00192          }
00193       }
00194 #define DDD csum
00195       for( jj=0 ; jj < ndim ; jj++ ){
00196          for( ii=0 ; ii < jj ; ii++ )
00197             nmat[jj+ii*ndim] = (nmat[ii+jj*ndim] *= DDD) ;
00198          nmat[jj+jj*ndim] *= DDD ;
00199       }
00200 
00201       /* check for convergence - L1 norm */
00202 
00203       cnorm = csum = 0.0 ;
00204       for( jj=0 ; jj < ndim ; jj++ ){
00205          for( ii=0 ; ii <= jj ; ii++ ){
00206             cnorm += fabs( nmat[ii+jj*ndim] - cmat[ii+jj*ndim] ) ;
00207             csum  += fabs( nmat[ii+jj*ndim] ) ;
00208          }
00209       }
00210 
00211 fprintf(stderr,"  |dif|=%12.4g  |mat|=%12.4g   cwt=%12.4g\n",cnorm,csum,cwt) ;
00212 fprintf(stderr,"  matrix:\n") ;
00213 for( ii=0 ; ii < ndim ; ii++ ){
00214    fprintf(stderr,"  Row%2d: %12.4g    ",ii,mv[ii]) ;
00215    for( jj=0 ; jj < ndim ; jj++ ) fprintf(stderr," %12.4g",nmat[ii+jj*ndim]) ;
00216    fprintf(stderr,"\n") ;
00217 }
00218 
00219       free(cv->cmat) ; free(cv->mvec) ;
00220       if( cnorm <= EPS*csum ){ cv->cmat = nmat; cv->mvec = mv; break; }  /* exit loop */
00221    }
00222 
00223    free(wv) ; free(tv) ; compute_choleski(cv) ; return cv ;
00224 }

AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

rcov.c File Reference

Data Structures

Defines

Functions

Define Documentation

Function Documentation