Doxygen Source Code Documentation
1ddot.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Functions | |
int | main (int argc, char *argv[]) |
Function Documentation
|
convert three sub-briks to a raw dataset with consecutive triplets Definition at line 3 of file 1ddot.c. References ADDTO_IMARR, argc, calloc, DESTROY_IMARR, far, free, IMARR_COUNT, IMARR_SUBIM, INIT_IMARR, malloc, MRI_FLOAT_PTR, mri_read_1D(), MRI_IMAGE::nx, MRI_IMAGE::ny, and symeig_double().
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 /* help? */ 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 /* options */ 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 /* input 1D files */ 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 /* create normalized vectors from 1D files */ 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 /* create matrix from dot product of vectors */ 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 /* print matrix out */ 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 /* compute eigendecomposition */ 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 /* compute matrix inverse */ 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 /* normalize matrix inverse */ 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 /* done */ 00201 00202 exit(0) ; 00203 } |