00001 /*********************** 3dDTeig.c **********************************************/
00002 /* Author: Daniel Glen, 15 Nov 2004 */
00003 #include "mrilib.h"
00004 #include "thd_shear3d.h"
00005 
00006 static char prefix[THD_MAX_PREFIX] = "eig" ;
00007 static int datum                   = MRI_float ;
00008 static void EIG_tsfunc( double tzero , double tdelta ,
00009                          int npts , float ts[] , double ts_mean ,
00010                          double ts_slope , void * ud , int nbriks, float * val ) ;
00011 
00012 int main( int argc , char * argv[] )
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 }
00145 
00146 /**********************************************************************
00147    Function that does the real work
00148 ***********************************************************************/
00149 
00150 static void EIG_tsfunc( double tzero, double tdelta ,
00151                           int npts, float ts[],
00152                           double ts_mean, double ts_slope,
00153                           void * ud, int nbriks, float * val          )
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 }
00271