Doxygen Source Code Documentation
3dDTeig.c File Reference
#include "mrilib.h"#include "thd_shear3d.h"Go to the source code of this file.
Functions | |
| void | EIG_tsfunc (double tzero, double tdelta, int npts, float ts[], double ts_mean, double ts_slope, void *ud, int nbriks, float *val) |
| int | main (int argc, char *argv[]) |
Variables | |
| char | prefix [THD_MAX_PREFIX] = "eig" |
| int | datum = MRI_float |
Function Documentation
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 150 of file 3dDTeig.c. References a, ENTRY, i, and symeig_double(). Referenced by main().
00154 {
00155 /* THD_dmat33 inmat;
00156 THD_dvecmat eigvmat;*/
00157 int i,j;
00158 static int nvox , ncall ;
00159 int maxindex, minindex, midindex;
00160 float temp, minvalue, maxvalue;
00161 int sortvector[3];
00162 double a[9], e[3];
00163 int astart, vstart;
00164 double ssq, dsq;
00165 double dv0, dv1, dv2;
00166
00167 ENTRY("EIG_tsfunc");
00168 /* ts is input vector data of 6 floating point numbers.
00169 For each point in volume brik convert vector data to
00170 symmetric matrix */
00171 /* ts should come from data sub-briks in form of Dxx,Dxy,Dxz,Dyy,Dyz,Dzz */
00172 /* convert to matrix of form
00173 [ Dxx Dxy Dxz]
00174 [ Dxy Dyy Dyz]
00175 [ Dxz Dyz Dzz] */
00176
00177 /** is this a "notification"? **/
00178 if( val == NULL ){
00179
00180 if( npts > 0 ){ /* the "start notification" */
00181
00182 nvox = npts ; /* keep track of */
00183 ncall = 0 ; /* number of calls */
00184
00185 } else { /* the "end notification" */
00186
00187 /* nothing to do here */
00188 }
00189 return ;
00190 }
00191
00192 /* load the symmetric matrix vector from the "timeseries" subbrik vector values */
00193
00194 a[0]=ts[0]; a[1]=ts[1]; a[2]=ts[2];
00195 a[3]=ts[1]; a[4]=ts[3]; a[5]=ts[4];
00196 a[6]=ts[2]; a[5]=ts[4]; a[8]=ts[5];
00197
00198 symeig_double(3, a, e); /* compute eigenvalues in e, eigenvectors in a */
00199
00200 maxindex=2; /* find the lowest, middle and highest eigenvalue */
00201 maxvalue=e[2];
00202 minindex=0;
00203 minvalue=e[0];
00204 midindex = 1;
00205
00206 for(i=0;i<3;i++) {
00207 temp = e[i];
00208 if(temp>maxvalue) { /* find the maximum */
00209 maxindex = i;
00210 maxvalue = temp;
00211 }
00212 if(temp<minvalue) { /* find the minimum */
00213 minindex = i;
00214 minvalue = temp;
00215 }
00216 }
00217
00218 for(i=0;i<3;i++){ /* find the middle */
00219 if((i!=maxindex) && (i!=minindex))
00220 {
00221 midindex = i;
00222 break;
00223 }
00224 }
00225
00226 sortvector[0] = maxindex;
00227 sortvector[1] = midindex;
00228 sortvector[2] = minindex;
00229
00230 /* put the eigenvalues at the beginning of the matrix */
00231 for(i=0;i<3;i++) {
00232 val[i] = e[sortvector[i]]; /* copy sorted eigenvalues */
00233
00234 /* start filling in eigenvector values */
00235 astart=sortvector[i]*3; /* start index of double eigenvector */
00236 vstart=(i+1)*3; /* start index of float val vector to copy eigenvector */
00237
00238 for(j=0;j<3;j++){
00239 val[vstart+j] = a[astart+j];
00240 }
00241 }
00242
00243 /* calculate the Fractional Anisotropy, FA */
00244 /* reference, Pierpaoli C, Basser PJ. Microstructural and physiological features
00245 of tissues elucidated by quantitative-diffusion tensor MRI,J Magn Reson B 1996; 111:209-19 */
00246 if((val[0]<=0.0)||(val[1]<=0.0)||(val[2]<=0.0)) { /* any negative eigenvalues?*/
00247 val[12]=0.0; /* set FA to 0 */
00248 val[13]=0.0; /* set MD to 0 */
00249 EXRETURN;
00250 }
00251
00252 ssq = (val[0]*val[0])+(val[1]*val[1])+(val[2]*val[2]); /* sum of squares of eigenvalues */
00253 /* dsq = pow((val[0]-val[1]),2.0) + pow((val[1]-val[2]),2.0) + pow((val[2]-val[0]),2.0);*/ /* sum of differences squared */
00254
00255 dv0 = val[0]-val[1];
00256 dv0 *= dv0;
00257 dv1 = val[1]-val[2];
00258 dv1 *= dv1;
00259 dv2 = val[2]-val[0];
00260 dv2 *= dv2;
00261 dsq = dv0+dv1+dv2; /* sum of differences squared */
00262
00263 if(ssq!=0.0)
00264 val[12] = sqrt(dsq/(2.0*ssq)); /* FA calculated here */
00265 else
00266 val[12] = 0.0;
00267 val[13] = (val[0]+val[1]+val[2]) / 3; /* MD - mean diffusivity=average of eigenvalues */
00268
00269 EXRETURN;
00270 }
|
|
||||||||||||
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 12 of file 3dDTeig.c. References ADN_brick_label_one, ADN_none, ADN_ntt, ADN_ttdel, ADN_ttorg, ADN_tunits, AFNI_logger(), argc, datum, DSET_BRIKNAME, DSET_NVALS, DSET_write, EDIT_dset_items(), EIG_tsfunc(), ISVALID_DSET, machdep(), mainENTRY, MAKER_4D_to_typed_fbuc(), MASTER_SHORTHELP_STRING, MCW_strncpy, prefix, PRINT_VERSION, THD_filename_ok(), THD_MAX_PREFIX, THD_open_dataset(), tross_Copy_History(), tross_Make_History(), and UNITS_SEC_TYPE.
00013 {
00014 THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */
00015 int nopt, nbriks ;
00016
00017 /*----- Read command line -----*/
00018 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00019 printf("Usage: 3dDTeig [options] dataset\n"
00020 "Computes eigenvalues and eigenvectors for an input dataset of\n"
00021 " 6 sub-bricks Dxx,Dxy,Dxz,Dyy,Dyz,Dzz.\n"
00022 " The results are stored in a 14-subbrick bucket dataset.\n"
00023 " The resulting 14-subbricks are\n"
00024 " lambda_1,lambda_2,lambda_3,\n"
00025 " eigvec_1[1-3],eigvec_2[1-3],eigvec_3[1-3],\n"
00026 " FA,MD.\n\n"
00027 "The output is a bucket dataset. The input dataset\n"
00028 "may use a sub-brick selection list, as in program 3dcalc.\n"
00029 " Mean diffusivity (MD) calculated as simple average of eigenvalues.\n"
00030 " Fractional Anisotropy (FA) calculated according to Pierpaoli C, Basser PJ.\n"
00031 " Microstructural and physiological features of tissues elucidated by\n"
00032 " quantitative-diffusion tensor MRI, J Magn Reson B 1996; 111:209-19\n"
00033 ) ;
00034 printf("\n" MASTER_SHORTHELP_STRING ) ;
00035 exit(0) ;
00036 }
00037
00038 mainENTRY("3dDTeig main"); machdep(); AFNI_logger("3dDTeig",argc,argv);
00039 PRINT_VERSION("3dDTeig") ;
00040
00041 nopt = 1 ;
00042 nbriks = 14 ;
00043 datum = MRI_float;
00044 while( nopt < argc && argv[nopt][0] == '-' ){
00045
00046 /*-- prefix --*/
00047
00048 if( strcmp(argv[nopt],"-prefix") == 0 ){
00049 if( ++nopt >= argc ){
00050 fprintf(stderr,"*** -prefix needs an argument!\n"); exit(1);
00051 }
00052 MCW_strncpy(prefix,argv[nopt],THD_MAX_PREFIX) ;
00053 if( !THD_filename_ok(prefix) ){
00054 fprintf(stderr,"*** %s is not a valid prefix!\n",prefix); exit(1);
00055 }
00056 nopt++ ; continue ;
00057 }
00058
00059 /*-- datum --*/
00060
00061 if( strcmp(argv[nopt],"-datum") == 0 ){
00062 if( ++nopt >= argc ){
00063 fprintf(stderr,"*** -datum needs an argument!\n"); exit(1);
00064 }
00065 if( strcmp(argv[nopt],"short") == 0 ){
00066 datum = MRI_short ;
00067 } else if( strcmp(argv[nopt],"float") == 0 ){
00068 datum = MRI_float ;
00069 } else if( strcmp(argv[nopt],"byte") == 0 ){
00070 datum = MRI_byte ;
00071 } else {
00072 fprintf(stderr,"-datum of type '%s' is not supported!\n",
00073 argv[nopt] ) ;
00074 exit(1) ;
00075 }
00076 nopt++ ; continue ;
00077 }
00078 }
00079
00080 /*----- read input dataset -----*/
00081
00082 if( nopt >= argc ){
00083 fprintf(stderr,"*** No input dataset!?\n"); exit(1);
00084 }
00085
00086 old_dset = THD_open_dataset( argv[nopt] ) ;
00087 if( !ISVALID_DSET(old_dset) ){
00088 fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]); exit(1);
00089 }
00090
00091 /* expect 6 values per voxel - 6 sub-briks as input dataset */
00092 if( DSET_NVALS(old_dset) < 6 ){ /* allows 6 or greater sub-briks */
00093 fprintf(stderr,"*** Can't use dataset that is not at least 6 values per voxel!\n") ;
00094 exit(1) ;
00095 }
00096
00097 EDIT_dset_items( old_dset ,
00098 ADN_ntt , DSET_NVALS(old_dset) ,
00099 ADN_ttorg , 0.0 ,
00100 ADN_ttdel , 1.0 ,
00101 ADN_tunits , UNITS_SEC_TYPE ,
00102 NULL ) ;
00103
00104
00105 /*------------- ready to compute new dataset -----------*/
00106
00107 new_dset = MAKER_4D_to_typed_fbuc(
00108 old_dset , /* input dataset */
00109 prefix , /* output prefix */
00110 datum , /* output datum */
00111 0 , /* ignore count */
00112 0 , /* can't detrend in maker function KRH 12/02*/
00113 nbriks , /* number of briks */
00114 EIG_tsfunc , /* timeseries processor */
00115 NULL /* data for tsfunc */
00116 ) ;
00117
00118 if( new_dset != NULL ){
00119 tross_Copy_History( old_dset , new_dset ) ;
00120 EDIT_dset_items(new_dset, ADN_brick_label_one+0, "lambda_1", ADN_none);
00121 EDIT_dset_items(new_dset, ADN_brick_label_one+1, "lambda_2",ADN_none);
00122 EDIT_dset_items(new_dset, ADN_brick_label_one+2, "lambda_3",ADN_none);
00123 EDIT_dset_items(new_dset, ADN_brick_label_one+3, "eigvec_1[1]",ADN_none);
00124 EDIT_dset_items(new_dset, ADN_brick_label_one+4, "eigvec_1[2]",ADN_none);
00125 EDIT_dset_items(new_dset, ADN_brick_label_one+5, "eigvec_1[3]",ADN_none);
00126 EDIT_dset_items(new_dset, ADN_brick_label_one+6, "eigvec_2[1]",ADN_none);
00127 EDIT_dset_items(new_dset, ADN_brick_label_one+7, "eigvec_2[2]",ADN_none);
00128 EDIT_dset_items(new_dset, ADN_brick_label_one+8, "eigvec_2[3]",ADN_none);
00129 EDIT_dset_items(new_dset, ADN_brick_label_one+9, "eigvec_3[1]",ADN_none);
00130 EDIT_dset_items(new_dset, ADN_brick_label_one+10,"eigvec_3[2]",ADN_none);
00131 EDIT_dset_items(new_dset, ADN_brick_label_one+11,"eigvec_3[3]",ADN_none);
00132 EDIT_dset_items(new_dset, ADN_brick_label_one+12,"FA",ADN_none);
00133 EDIT_dset_items(new_dset, ADN_brick_label_one+13,"MD",ADN_none);
00134
00135 tross_Make_History( "3dDTeig" , argc,argv , new_dset ) ;
00136 DSET_write( new_dset ) ;
00137 fprintf(stderr,"--- Output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
00138 } else {
00139 fprintf(stderr,"*** Unable to compute output dataset!\n") ;
00140 exit(1) ;
00141 }
00142
00143 exit(0) ;
00144 }
|
Variable Documentation
|
|
Definition at line 7 of file 3dDTeig.c. Referenced by main(). |
|
|
Definition at line 6 of file 3dDTeig.c. Referenced by main(). |
