00001 /*********************** 3dDTtoDWI.c **********************************************/
00002 /* Author: Daniel Glen, 17 Nov 2004 */
00003 /* compute multiple gradient images from 6 principle direction tensors*/
00004 /* and corresponding gradient vector coordinates */
00005 /* used for modeling diffusion tensor data */
00006 
00007 #include "thd_shear3d.h"
00008 /*#ifndef FLOATIZE*/
00009 # include "matrix.h"
00010 # include "matrix.c"
00011 /*#endif*/
00012 #include "afni.h"
00013 
00014 #define TINYNUMBER 1E-10
00015 #define SMALLNUMBER 1E-4
00016 
00017 static char prefix[THD_MAX_PREFIX] = "DT";
00018 static int datum = MRI_float;
00019 static matrix Dmatrix;
00020 static vector Dvector;
00021 
00022 static double *bmatrix;         /* b matrix = GiGj for gradient intensities */
00023 
00024 static float *I0_ptr;           /* I0 matrix */
00025 static int automask = 0;        /* automasking flag - user option */
00026 
00027 static void DTtoDWI_tsfunc (double tzero, double tdelta, int npts, float ts[], double ts_mean, double ts_slope, void *ud, int nbriks, float *val);
00028 static void Computebmatrix (MRI_IMAGE * grad1Dptr);
00029 static void InitGlobals (int npts);
00030 static void FreeGlobals (void);
00031 
00032 int
00033 main (int argc, char *argv[])
00034 {
00035   THD_3dim_dataset *old_dset, *new_dset, *I0_dset;      /* input and output datasets */
00036   int nopt, nbriks, nvox;
00037   int i, eigs_brik;
00038   MRI_IMAGE *grad1Dptr = NULL;
00039   MRI_IMAGE *anat_im = NULL;
00040   MRI_IMAGE *data_im = NULL;
00041   double fac;
00042   short *sar = NULL, *tempsptr = NULL, tempval;
00043   byte *maskptr = NULL, *tempbptr = NULL;
00044   char tempstr[25];
00045 
00046    /*----- Read command line -----*/
00047   if (argc < 2 || strcmp (argv[1], "-help") == 0)
00048     {
00049       printf ("Usage: 3dDTtoDWI [options] gradient-file I0-dataset DT-dataset\n"
00050               "Computes  multiple gradient images from 6 principle direction tensors and\n"
00051               "    corresponding gradient vector coordinates applied to the I0-dataset.\n"
00052               " The program takes three parameters as input :  \n"
00053               "    a 1D file of the gradient vectors with lines of ASCII floats Gxi,Gyi,Gzi.\n"
00054               "    Only the non-zero gradient vectors are included in this file (no G0 line).\n"
00055               " The I0 dataset is a volume without any gradient applied.\n"
00056               " The DT dataset is the 6-sub-brick dataset containing the diffusion tensor data,\n"
00057               "    Dxx, Dxy, Dxz, Dyy, Dyz, Dzz\n"
00058               " Options:\n"
00059               "   -prefix pname = Use 'pname' for the output dataset prefix name.\n"
00060               "    [default='DWI']\n"
00061               "   -automask =  mask dataset so that the gradient images are computed only for\n"
00062               "    high-intensity (presumably brain) voxels.  The intensity level is\n"
00063               "    determined the same way that 3dClipLevel works.\n\n"
00064               "   -datum type = output dataset type [float/short/byte] (default is float).\n"
00065               "   -help = show this help screen.\n"
00066               " Example:\n"
00067               "  3dDTtoDWI -prefix DWI -automask tensor25.1D 'DT+orig[26]' DT+orig.\n\n"
00068               " The output is a n sub-brick bucket dataset containing computed DWI images.\n"
00069               "    where n is the number of vectors in the gradient file + 1\n"
00070               "\n");
00071       printf ("\n" MASTER_SHORTHELP_STRING);
00072       exit (0);
00073     }
00074 
00075   mainENTRY ("3dDTtoDWI main");
00076   machdep ();
00077   AFNI_logger ("3dDTtoDWI", argc, argv);
00078   PRINT_VERSION("3dDTtoDWI") ;
00079 
00080   nopt = 1;
00081   datum = MRI_float;
00082 
00083 
00084   while (nopt < argc && argv[nopt][0] == '-')
00085     {
00086 
00087       /*-- prefix --*/
00088 
00089       if (strcmp (argv[nopt], "-prefix") == 0)
00090         {
00091           if (++nopt >= argc)
00092             {
00093               fprintf (stderr, "*** Error - prefix needs an argument!\n");
00094               exit (1);
00095             }
00096           MCW_strncpy (prefix, argv[nopt], THD_MAX_PREFIX);     /* change name from default prefix */
00097           /* check file name to be sure not to overwrite - mod drg 12/9/2004 */
00098           if (!THD_filename_ok (prefix))
00099             {
00100               fprintf (stderr, "*** Error - %s is not a valid prefix!\n", prefix);
00101               exit (1);
00102             }
00103           nopt++;
00104           continue;
00105         }
00106 
00107       /*-- datum --*/
00108 
00109       if (strcmp (argv[nopt], "-datum") == 0)
00110         {
00111           if (++nopt >= argc)
00112             {
00113               fprintf (stderr, "*** Error - datum needs an argument!\n");
00114               exit (1);
00115             }
00116           if (strcmp (argv[nopt], "short") == 0)
00117             {
00118               datum = MRI_short;
00119             }
00120           else if (strcmp (argv[nopt], "float") == 0)
00121             {
00122               datum = MRI_float;
00123             }
00124           else if (strcmp (argv[nopt], "byte") == 0)
00125             {
00126               datum = MRI_byte;
00127             }
00128           else
00129             {
00130               fprintf (stderr, "-datum of type '%s' is not supported!\n",
00131                        argv[nopt]);
00132               exit (1);
00133             }
00134           nopt++;
00135           continue;
00136         }
00137       if (strcmp (argv[nopt], "-automask") == 0)
00138         {
00139           automask = 1;
00140           nopt++;
00141           continue;
00142         }
00143 
00144         fprintf(stderr, "*** Error - unknown option %s\n", argv[nopt]);
00145         exit(1);
00146     }
00147   
00148    /*----- read input datasets -----*/
00149 
00150   if (nopt >= argc)
00151     {
00152       fprintf (stderr, "*** Error - No input dataset!?\n");
00153       exit (1);
00154     }
00155 
00156   /* first input dataset - should be gradient vector file of ascii floats Gx,Gy,Gz */
00157 
00158   /* read gradient vector 1D file */
00159   grad1Dptr = mri_read_1D (argv[nopt]);
00160   if (grad1Dptr == NULL)
00161     {
00162       fprintf (stderr, "*** Error reading gradient vector file\n");
00163       exit (1);
00164     }
00165 
00166   if (grad1Dptr->ny != 3)
00167     {
00168       fprintf (stderr, "*** Error - Only 3 columns of gradient vectors allowed\n");
00169       fprintf (stderr, "%d columns found\n", grad1Dptr->nx);
00170       mri_free (grad1Dptr);
00171       exit (1);
00172     }
00173 
00174   if (grad1Dptr->nx < 6)
00175     {
00176       fprintf (stderr, "*** Error - Must have at least 6 gradient vectors\n");
00177       fprintf (stderr, "%d columns found\n", grad1Dptr->nx);
00178       mri_free (grad1Dptr);
00179       exit (1);
00180     }
00181 
00182   nbriks = grad1Dptr->nx + 1;    /* number of gradients specified here from file */     
00183   nopt++;
00184 
00185   /* open I0 dataset - idealized no gradient image */
00186   I0_dset = THD_open_dataset (argv[nopt]);
00187   if (!ISVALID_DSET (I0_dset))
00188     {
00189        fprintf (stderr, "*** Error - Can not open dataset %s\n", argv[nopt]);
00190        mri_free (grad1Dptr);
00191        exit (1);
00192     }
00193 
00194    DSET_mallocize (I0_dset);
00195    DSET_load (I0_dset);                 /* load dataset */
00196    data_im = DSET_BRICK (I0_dset, 0);   /* set pointer to the 0th sub-brik of the dataset */
00197    fac = DSET_BRICK_FACTOR(I0_dset, 0); /* get scale factor for each sub-brik*/
00198    if(fac==0.0) fac=1.0;
00199    if((data_im->kind != MRI_float) || (fac!=1.0)) {
00200        fprintf (stderr, "*** Error - Can only open float datasets with scale factors of 1\n");
00201        mri_free (grad1Dptr);
00202        mri_free (data_im);
00203        exit (1);
00204    }
00205 
00206    I0_ptr = mri_data_pointer(data_im) ; /* pointer to I0 data */
00207 
00208    nopt++;
00209 
00210   /* Now read in all the MRI volumes for each gradient vector */
00211   /* assumes first one is no gradient */
00212   old_dset = THD_open_dataset (argv[nopt]);
00213 
00214   if (!ISVALID_DSET (old_dset))
00215     {
00216       fprintf (stderr, "*** Error - Can not open dataset %s\n", argv[nopt]);
00217       mri_free (grad1Dptr);
00218       exit (1);
00219     }
00220 
00221   /* expect at least 6 values per voxel - 6 sub-briks as input dataset */
00222   if (DSET_NVALS (old_dset) <6)
00223     {
00224       fprintf (stderr,
00225       "*** Error - Dataset must have at least 6 sub-briks to describe the diffusion tensor\n");
00226       mri_free (grad1Dptr);
00227       mri_free (data_im);
00228       exit (1);
00229     }
00230 
00231 
00232   InitGlobals (grad1Dptr->nx + 1);      /* initialize all the matrices and vectors */
00233   Computebmatrix (grad1Dptr);   /* compute bij=GiGj */
00234 
00235   if (automask)
00236     {
00237       DSET_mallocize (old_dset);
00238       DSET_load (old_dset);     /* get B0 (anatomical image) from dataset */
00239       /*anat_im = THD_extract_float_brick( 0, old_dset ); */
00240       anat_im = DSET_BRICK (old_dset, 0);       /* set pointer to the 0th sub-brik of the dataset */
00241       maskptr = mri_automask_image (anat_im);   /* maskptr is a byte pointer for volume */
00242 
00243       /* convert byte mask to same format type as dataset */
00244       nvox = DSET_NVOX (old_dset);
00245       sar = (short *) calloc (nvox, sizeof (short));
00246       /* copy maskptr values to far ptr */
00247       tempsptr = sar;
00248       tempbptr = maskptr;
00249       for (i = 0; i < nvox; i++)
00250         {
00251           *tempsptr++ = (short) *tempbptr++;
00252           tempval = *(tempsptr - 1);
00253         }
00254 
00255       free (maskptr);
00256 
00257       /*old_dset->dblk->malloc_type = DATABLOCK_MEM_MALLOC; *//* had to set this? */
00258       EDIT_add_brick (old_dset, MRI_short, 0.0, sar);   /* add sub-brik to end */
00259 
00260 
00261     }
00262 
00263   /* temporarily set artificial timing to 1 second interval */
00264   EDIT_dset_items (old_dset,
00265                    ADN_ntt, DSET_NVALS (old_dset),
00266                    ADN_ttorg, 0.0,
00267                    ADN_ttdel, 1.0, ADN_tunits, UNITS_SEC_TYPE, NULL);
00268 
00269    /*------------- ready to compute new dataset -----------*/
00270 
00271   new_dset = MAKER_4D_to_typed_fbuc (old_dset,  /* input dataset */
00272                                      prefix,    /* output prefix */
00273                                      datum,     /* output datum  */
00274                                      0, /* ignore count  */
00275                                      0, /* can't detrend in maker function  KRH 12/02 */
00276                                      nbriks,    /* number of briks */
00277                                      DTtoDWI_tsfunc,    /* timeseries processor */
00278                                      NULL       /* data for tsfunc */
00279     );
00280 
00281 
00282 
00283   FreeGlobals ();
00284   mri_free (grad1Dptr);
00285 
00286 
00287   if (automask)
00288     {
00289       mri_free (anat_im);
00290       DSET_unload_one (old_dset, 0);
00291       sar = NULL;
00292     }
00293 
00294   if (new_dset != NULL)
00295     {
00296       tross_Copy_History (old_dset, new_dset);
00297       for(i=0;i<nbriks;i++) {
00298         sprintf(tempstr,"grad%3.3d", i);
00299         EDIT_dset_items (new_dset, ADN_brick_label_one + i, tempstr, ADN_none);
00300       }
00301       tross_Make_History ("3dDTtoDWI", argc, argv, new_dset);
00302       DSET_write (new_dset);
00303       fprintf(stderr,"--- Output dataset %s\n", DSET_BRIKNAME(new_dset));
00304     }
00305   else
00306     {
00307       fprintf (stderr, "*** Error - Unable to compute output dataset!\n");
00308       exit (1);
00309     }
00310 
00311   exit (0);
00312 }
00313 
00314 /**********************************************************************
00315    Function that does the real work
00316 ***********************************************************************/
00317 
00318 static void
00319 DTtoDWI_tsfunc (double tzero, double tdelta,
00320                 int npts, float ts[],
00321                 double ts_mean, double ts_slope,
00322                 void *ud, int nbriks, float *val)
00323 {
00324   int i, j;
00325   static int nvox, ncall;
00326   int allzeros;
00327   double I0, bq_d;
00328   double *bptr;
00329   float *tempptr;
00330  
00331   ENTRY ("DTtoDWI_tsfunc");
00332   /* ts is input vector data of 6 floating point numbers.*/
00333   /* ts should come from data sub-briks in form of Dxx, Dxy, Dxz, Dyy, Dyz, Dzz */
00334   /* if automask is turned on, ts has 7 floating point numbers */
00335   /* val is output vector of form DWI0, DWI1, DWIn sub-briks */
00336   /* where n = number of gradients = nbriks-1 */
00337 
00338    /** is this a "notification"? **/
00339    if (val == NULL)
00340     {
00341 
00342       if (npts > 0)
00343         {                       /* the "start notification" */
00344           nvox = npts;          /* keep track of   */
00345           ncall = 0;            /* number of calls */
00346         }
00347       else
00348         {                       /* the "end notification" */
00349 
00350           /* nothing to do here */
00351         }
00352       EXRETURN;
00353     }
00354 
00355   ncall++;
00356 
00357   if (automask)
00358      npts = npts - 1;
00359 
00360   tempptr = I0_ptr+ncall-1;
00361   I0 = *tempptr;
00362 
00363 #if 0
00364   /* check for reasons not to calculate anything and just return zeros */
00365 if(ncall==1)
00366    printf("checking for all zeros\n");
00367   allzeros = 0;
00368   i=0;
00369 
00370   while ((i<6)&&(I0!=0.0))        /* check if all the DT values are 0 */
00371   { 
00372      if(ts[i++]!=0)
00373         allzeros = 0;
00374   }
00375 
00376 
00377   for(i=0;i<nbriks;i++)
00378     val[i] = 0.0;
00379   EXRETURN;
00380 
00381 
00382 /* return zeros if all the DT values are 0, if the I0 value is 0 or 
00383    the mask is off at this voxel*/
00384   if(allzeros||(I0==0.0) || (automask&&(ts[npts] == 0)))
00385     {
00386           for (i = 0; i < nbriks; i++)
00387             val[i] = 0.0;       /* return 0 for all DWIn points */
00388           EXRETURN;
00389     }
00390 
00391 #endif
00392 
00393   val[0] = I0; /* the first sub-brik is the I0 sub-brik */
00394   bptr = bmatrix+6;   /* start at the first gradient */
00395 
00396   for(i=1;i<nbriks;i++) {
00397      bptr = bmatrix+(6*i);   /* start at the first gradient */
00398      bq_d = *bptr++ * ts[0];          /* GxGxDxx  */
00399      bq_d += *bptr++ * ts[1] * 2;      /* 2GxGyDxy */
00400      bq_d += *bptr++ * ts[2] * 2;      /* 2GxGzDxz */
00401      bq_d += *bptr++ * ts[3];          /* GyGyDyy  */
00402      bq_d += *bptr++ * ts[4] * 2;      /* 2GyGzDyz */
00403      bq_d += *bptr++ * ts[5];          /* GzGzDzz  */
00404 
00405      val[i] = I0 * exp(-bq_d);   /* for each gradient,q, Iq = J e -(bq.D) */
00406                                  /* where bq.D is the large dot product of bq and D */
00407   }
00408 
00409   EXRETURN;
00410 }
00411 
00412 
00413 /*! compute the diffusion weighting matrix bmatrix for q number of gradients */
00414 /* only need to calculate once */
00415 /* bq = diffusion weighting matrix of qth gradient */
00416 /*      GxGx GxGy GxGz
00417         GxGy GyGy GyGz
00418         GxGz GyGz GzGz */
00419 /* b0 is 0 for all 9 elements */
00420 /* bmatrix is really stored as 6 x npts array */
00421 static void
00422 Computebmatrix (MRI_IMAGE * grad1Dptr)
00423 {
00424   int i, j, n;
00425   register double *bptr;
00426   register float *Gxptr, *Gyptr, *Gzptr;
00427   double Gx, Gy, Gz;
00428 
00429   ENTRY ("Computebmatrix");
00430   n = grad1Dptr->nx;            /* number of gradients other than I0 */
00431   Gxptr = MRI_FLOAT_PTR (grad1Dptr);    /* use simple floating point pointers to get values */
00432   Gyptr = Gxptr + n;
00433   Gzptr = Gyptr + n;
00434 
00435   bptr = bmatrix;
00436   for (i = 0; i < 6; i++)
00437     *bptr++ = 0.0;              /* initialize first 6 elements to 0.0 for the I0 gradient */
00438 
00439   for (i = 0; i < n; i++)
00440     {
00441       Gx = *Gxptr++;
00442       Gy = *Gyptr++;
00443       Gz = *Gzptr++;
00444       *bptr++ = Gx * Gx;
00445       *bptr++ = Gx * Gy;
00446       *bptr++ = Gx * Gz;
00447       *bptr++ = Gy * Gy;
00448       *bptr++ = Gy * Gz;
00449       *bptr++ = Gz * Gz;
00450       for(j=0;j<6;j++)
00451          printf("%-13.6g ", *(bmatrix + 6 + (i*6)+ j) );
00452       printf("\n");
00453     }
00454 
00455 
00456   EXRETURN;
00457 }
00458 
00459 /*! allocate all the global matrices and arrays once */
00460 static void
00461 InitGlobals (int npts)
00462 {
00463   int i;
00464   double *cumulativewtptr;
00465 
00466   ENTRY ("InitGlobals");
00467 
00468   bmatrix = malloc (npts * 6 * sizeof (double));
00469 
00470   EXRETURN;
00471 }
00472 
00473 /*! free up all the matrices and arrays */
00474 static void
00475 FreeGlobals ()
00476 {
00477 
00478   ENTRY ("FreeGlobals");
00479   free (bmatrix);
00480   bmatrix = NULL;
00481   EXRETURN;
00482 }