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) ;
00093 tv = (float *) malloc(sizeof(float)*ndim) ;
00094 mv = (float *) malloc(sizeof(float)*ndim) ;
00095 wv = (float *) malloc(sizeof(float)*nvec) ;
00096
00097 fnvec = 1.0/nvec ; fndim = 1.0/ndim ;
00098 bcut = 1.0 + CCUT*sqrt(fndim) ;
00099
00100
00101
00102 for( jj=0 ; jj < ndim ; jj++ ) mv[jj] = 0.0 ;
00103
00104 for( kk=0 ; kk < nvec ; kk++ ){
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 ;
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++ ){
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++ ){
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
00127
00128 nite = 0 ;
00129
00130 while(1){
00131
00132 nite++ ;
00133
00134 cmat = cv->cmat = nmat ;
00135 compute_choleski(cv) ;
00136
00137 nmat = (float *) malloc(sizeof(float)*ndim*ndim) ;
00138 cv->mvec = mv ;
00139 mv = (float *) malloc(sizeof(float)*ndim) ;
00140
00141 for( jj=0 ; jj < ndim ; jj++ ){
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
00149
00150 csum = 0.0 ;
00151 for( kk=0 ; kk < nvec ; kk++ ){
00152 vv = vec[kk] ;
00153
00154
00155
00156
00157 for( jj=0 ; jj < ndim ; jj++ ) tv[jj] = vv[jj] - cv->mvec[jj] ;
00158 forward_solve_inplace(cv,tv) ;
00159
00160
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
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 ;
00183
00184
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
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; }
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 }