Doxygen Source Code Documentation
rcov.c File Reference
#include "mrilib.h"
Go to the source code of this file.
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
|
Definition at line 72 of file rcov.c. Referenced by robust_covar(). |
|
|
|
|
|
|
|
Definition at line 73 of file rcov.c. Referenced by robust_covar(). |
Function Documentation
|
Definition at line 28 of file rcov.c. References covmat::cfac, CH, covmat::ndim, and vec. Referenced by evaluate_span().
|
|
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 } |
|
Definition at line 13 of file rcov.c. References covmat::cfac, CH, covmat::ndim, and vec. Referenced by evaluate_span(), and robust_covar().
|
|
\** File : SUMA.c
Input paramters :
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 } |
|
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 } |