Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
1ddot.c
Go to the documentation of this file.00001 #include "mrilib.h"
00002
00003 int main( int argc , char *argv[] )
00004 {
00005 int iarg , ii,jj,kk,mm , nvec , do_one=0 , nx=0,ny , ff ;
00006 MRI_IMAGE *tim ;
00007 MRI_IMARR *tar ;
00008 double sum , *eval , *amat , **tvec , *bmat ;
00009 float *far ;
00010
00011
00012
00013 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00014 printf("Usage: 1ddot [options] 1Dfile 1Dfile ...\n"
00015 "- Prints out correlation matrix of the 1D files and\n"
00016 " their inverse correlation matrix.\n"
00017 "- Output appears on stdout.\n"
00018 "\n"
00019 "Options:\n"
00020 " -one = Make 1st vector be all 1's.\n"
00021 ) ;
00022 exit(0) ;
00023 }
00024
00025
00026
00027 iarg = 1 ; nvec = 0 ;
00028 while( iarg < argc && argv[iarg][0] == '-' ){
00029
00030 if( strcmp(argv[iarg],"-one") == 0 ){
00031 do_one = 1 ; iarg++ ; continue ;
00032 }
00033
00034 fprintf(stderr,"** Unknown option: %s\n",argv[iarg]); exit(1);
00035 }
00036
00037 if( iarg == argc ){
00038 fprintf(stderr,"** No 1D files on command line!?\n"); exit(1);
00039 }
00040
00041
00042
00043 ff = iarg ;
00044 INIT_IMARR(tar) ; if( do_one ) nvec = 1 ;
00045 for( ; iarg < argc ; iarg++ ){
00046 tim = mri_read_1D( argv[iarg] ) ;
00047 if( tim == NULL ){
00048 fprintf(stderr,"** Can't read 1D file %s\n",argv[iarg]); exit(1);
00049 }
00050 if( nx == 0 ){
00051 nx = tim->nx ;
00052 } else if( tim->nx != nx ){
00053 fprintf(stderr,"** 1D file %s doesn't match first file in length!\n",
00054 argv[iarg]); exit(1);
00055 }
00056 nvec += tim->ny ;
00057 ADDTO_IMARR(tar,tim) ;
00058 }
00059
00060 printf("\n") ;
00061 printf("++ 1ddot input vectors:\n") ;
00062 jj = 0 ;
00063 if( do_one ){
00064 printf("00..00: all ones\n") ; jj = 1 ;
00065 }
00066 for( mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
00067 tim = IMARR_SUBIM(tar,mm) ;
00068 printf("%02d..%02d: %s\n", jj,jj+tim->ny-1, argv[ff+mm] ) ;
00069 jj += tim->ny ;
00070 }
00071
00072
00073
00074 tvec = (double **) malloc( sizeof(double *)*nvec ) ;
00075 for( jj=0 ; jj < nvec ; jj++ )
00076 tvec[jj] = (double *) malloc( sizeof(double)*nx ) ;
00077
00078 kk = 0 ;
00079 if( do_one ){
00080 sum = 1.0 / sqrt((double)nx) ;
00081 for( ii=0 ; ii < nx ; ii++ ) tvec[0][ii] = sum ;
00082 kk = 1 ;
00083 }
00084
00085 for( mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
00086 tim = IMARR_SUBIM(tar,mm) ;
00087 far = MRI_FLOAT_PTR(tim) ;
00088 for( jj=0 ; jj < tim->ny ; jj++,kk++ ){
00089 sum = 0.0 ;
00090 for( ii=0 ; ii < nx ; ii++ ){
00091 tvec[kk][ii] = far[ii+jj*nx] ;
00092 sum += tvec[kk][ii] * tvec[kk][ii] ;
00093 }
00094 if( sum == 0.0 ){
00095 fprintf(stderr,"** Input column is all zero!\n"); exit(1);
00096 }
00097 sum = 1.0 / sqrt(sum) ;
00098 for( ii=0 ; ii < nx ; ii++ ) tvec[kk][ii] *= sum ;
00099 }
00100 }
00101 DESTROY_IMARR(tar) ;
00102
00103
00104
00105 amat = (double *) calloc( sizeof(double) , nvec*nvec ) ;
00106
00107 for( kk=0 ; kk < nvec ; kk++ ){
00108 for( jj=0 ; jj <= kk ; jj++ ){
00109 sum = 0.0 ;
00110 for( ii=0 ; ii < nx ; ii++ ) sum += tvec[jj][ii] * tvec[kk][ii] ;
00111 amat[jj+nvec*kk] = sum ;
00112 if( jj < kk ) amat[kk+nvec*jj] = sum ;
00113 }
00114 }
00115
00116 for( jj=0 ; jj < nvec ; jj++ ) free(tvec[jj]) ;
00117 free(tvec) ;
00118
00119
00120
00121 printf("++ Correlation Matrix:\n ") ;
00122 for( jj=0 ; jj < nvec ; jj++ ) printf(" %02d ",jj) ;
00123 printf("\n ") ;
00124 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
00125 printf("\n") ;
00126 for( kk=0 ; kk < nvec ; kk++ ){
00127 printf("%02d:",kk) ;
00128 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",amat[jj+kk*nvec]) ;
00129 printf("\n") ;
00130 }
00131
00132
00133
00134 eval = (double *) malloc( sizeof(double)*nvec ) ;
00135 symeig_double( nvec , amat , eval ) ;
00136
00137 printf("\n"
00138 "++ Eigensolution:\n " ) ;
00139 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",eval[jj]) ;
00140 printf("\n ") ;
00141 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
00142 printf("\n") ;
00143 for( kk=0 ; kk < nvec ; kk++ ){
00144 printf("%02d:",kk) ;
00145 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",amat[kk+jj*nvec]) ;
00146 printf("\n") ;
00147 }
00148
00149
00150
00151 for( jj=0 ; jj < nvec ; jj++ ) eval[jj] = 1.0 / eval[jj] ;
00152
00153 bmat = (double *) calloc( sizeof(double) , nvec*nvec ) ;
00154
00155 for( ii=0 ; ii < nvec ; ii++ ){
00156 for( jj=0 ; jj < nvec ; jj++ ){
00157 sum = 0.0 ;
00158 for( kk=0 ; kk < nvec ; kk++ )
00159 sum += amat[ii+kk*nvec] * amat[jj+kk*nvec] * eval[kk] ;
00160 bmat[ii+jj*nvec] = sum ;
00161 }
00162 }
00163
00164 printf("\n") ;
00165 printf("++ Correlation Matrix Inverse:\n ") ;
00166 for( jj=0 ; jj < nvec ; jj++ ) printf(" %02d ",jj) ;
00167 printf("\n ") ;
00168 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
00169 printf("\n") ;
00170 for( kk=0 ; kk < nvec ; kk++ ){
00171 printf("%02d:",kk) ;
00172 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",bmat[jj+kk*nvec]) ;
00173 printf("\n") ;
00174 }
00175
00176
00177
00178 for( ii=0 ; ii < nvec ; ii++ ){
00179 for( jj=0 ; jj < nvec ; jj++ ){
00180 sum = bmat[ii+ii*nvec] * bmat[jj+jj*nvec] ;
00181 if( sum > 0.0 )
00182 amat[ii+jj*nvec] = bmat[ii+jj*nvec] / sqrt(sum) ;
00183 else
00184 amat[ii+jj*nvec] = 0.0 ;
00185 }
00186 }
00187
00188 printf("\n") ;
00189 printf("++ Correlation Matrix Inverse Normalized:\n ") ;
00190 for( jj=0 ; jj < nvec ; jj++ ) printf(" %02d ",jj) ;
00191 printf("\n ") ;
00192 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
00193 printf("\n") ;
00194 for( kk=0 ; kk < nvec ; kk++ ){
00195 printf("%02d:",kk) ;
00196 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",amat[jj+kk*nvec]) ;
00197 printf("\n") ;
00198 }
00199
00200
00201
00202 exit(0) ;
00203 }
