Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation

Doxygen Source Code Documentation


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

rcov.c

Go to the documentation of this file.
00001 #include "mrilib.h"
00002 
00003 typedef struct {
00004    int ndim ;
00005    float * cmat , * cfac , * mvec ;
00006 } covmat ;
00007 
00008 #define CM(i,j) cmat[(i)+(j)*ndim]
00009 #define CH(i,j) cfac[(i)+(j)*ndim]
00010 
00011 /*-----------------------------------------------------------------*/
00012 
00013 void forward_solve_inplace( covmat * cv , float * vec )
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 }
00025 
00026 /*-----------------------------------------------------------------*/
00027 
00028 void backward_solve_inplace( covmat * cv , float * vec )
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 }
00040 
00041 /*-----------------------------------------------------------------*/
00042 
00043 void compute_choleski( covmat * cv )
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 }
00069 
00070 /*-----------------------------------------------------------------*/
00071 
00072 #define CCUT 3.5
00073 #define EPS  1.e-4
00074 
00075 covmat * robust_covar( int ndim , int nvec , float ** vec )
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 }
00225 
00226 /*----------------------------------------------------------------------------*/
00227 
00228 int main( int argc , char * argv[] )
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 }
 

Powered by Plone

This site conforms to the following standards: