00001 /*********************** 3dDWItoDT.c **********************************************/
00002 /* Author: Daniel Glen, 17 Nov 2004 */
00003 /* compute 6 principle direction tensors from multiple gradient vectors*/
00004 /* and corresponding image data */
00005 /* This version includes a more complex algorithm that takes into account*/
00006 /* noise.*/
00007 
00008 #include "thd_shear3d.h"
00009 /*#ifndef FLOATIZE*/
00010 # include "matrix.h"
00011 # include "matrix.c"
00012 /*#endif*/
00013 #include "afni.h"
00014 
00015 #define TINYNUMBER 1E-10
00016 #define SMALLNUMBER 1E-4
00017 #define MAX_CONVERGE_STEPS 10           /* default maximum steps */
00018 #define MAX_RWCONVERGE_STEPS 5
00019 
00020 static char prefix[THD_MAX_PREFIX] = "DT";
00021 static int datum = MRI_float;
00022 static matrix Rtmat;
00023 static double *Rvector;         /* residuals at each gradient */
00024 static double *tempRvector;     /* copy of residuals at each gradient */
00025 static matrix Fmatrix;
00026 static matrix Dmatrix;
00027 static matrix OldD;
00028 static matrix Hplusmatrix, Hminusmatrix;
00029 static vector Dvector;
00030 static matrix tempFmatrix[2];
00031 static matrix tempDmatrix[2];
00032 static matrix tempHplusmatrix[2], tempHminusmatrix[2];
00033 /* static vector tempDvector; */
00034 
00035 static byte *maskptr = NULL;
00036 static double eigs[12];
00037 static double deltatau;
00038 static double *wtfactor;        /* weight factors for time points at each voxel */
00039 static double *bmatrix;         /* b matrix = GiGj for gradient intensities */
00040 static double *cumulativewt;    /* overall wt. factor for each gradient */
00041 static long rewtvoxels;         /* how many voxels were reweighted */ 
00042 static double sigma;            /* std.deviation */
00043 static double ED;               /* error for each iteration - cost function result */
00044 static int automask = 0;        /* automasking flag - user option */
00045 static int reweight_flag = 0;   /* reweight computation flag - user option */
00046 static int method = -1;         /* linear or non-linear method - user option */
00047 static int max_iter = -2;       /* maximum #convergence iteration steps - user option */
00048 static int max_iter_rw = -2;    /* maximum #convergence iteration steps - user option */
00049 static int eigs_flag = 0;       /* eigenvalue calculation in output - user option */
00050 static int cumulative_flag = 0; /* calculate, display cumulative wts for gradients - user option */ 
00051 static int debug_briks = 0;     /* put Ed, Ed0 and Converge step sub-briks in output - user option */
00052 static int verbose = 0;         /* print out info every verbose number of voxels - user option */
00053 static int afnitalk_flag = 0;   /* show convergence in AFNI graph - user option */
00054 
00055 static NI_stream_type * DWIstreamid = 0;     /* NIML stream ID */
00056 
00057 static void Form_R_Matrix (MRI_IMAGE * grad1Dptr);
00058 static void DWItoDT_tsfunc (double tzero, double tdelta, int npts, float ts[], double ts_mean, double ts_slope, void *ud, int nbriks, float *val);
00059 static void EIG_func (void);
00060 static float Calc_FA(float *val);
00061 static float Calc_MD(float *val);
00062 static void ComputeD0 (void);
00063 static double ComputeJ (float ts[], int npts);
00064 static void ComputeDeltaTau (void);
00065 static void Computebmatrix (MRI_IMAGE * grad1Dptr);
00066 static void InitGlobals (int npts);
00067 static void FreeGlobals (void);
00068 static void Store_Computations (int i, int npts, int converge_step);
00069 static void Restore_Computations (int i, int npts, int converge_step);
00070 static void InitWtfactors (int npts);
00071 static void ComputeWtfactors (int npts);
00072 static void ComputeHpHm (double deltatau);
00073 static void ComputeNewD (void);
00074 static int TestConvergence (matrix NewD, matrix OldD);
00075 static void udmatrix_to_vector (matrix m, vector * v);
00076 static void udmatrix_copy (double *udptr, matrix * m);
00077 static double matrix_sumabs (matrix m);
00078 static double *InvertSym3 (double a, double b, double c, double e, double f,
00079                            double i);
00080 static void matrix_copy (matrix a, matrix * b);
00081 static int DWI_Open_NIML_stream(void);
00082 static int DWI_NIML_create_graph(void);
00083 static void DWI_AFNI_update_graph(double *Edgraph, double *dtau, int npts);
00084 
00085 int
00086 main (int argc, char *argv[])
00087 {
00088   THD_3dim_dataset *old_dset, *new_dset;        /* input and output datasets */
00089   int nopt, nbriks;
00090   int i, eigs_brik;
00091   MRI_IMAGE *grad1Dptr = NULL;
00092   MRI_IMAGE *anat_im = NULL;
00093 #if 0
00094   int nvox;
00095   short *sar = NULL;
00096   short *tempsptr = NULL;
00097   byte *tempbptr = NULL;
00098   short tempval;
00099 #endif
00100 
00101   double *cumulativewtptr;
00102   int mmvox=0 ;
00103   int nxyz;
00104 
00105 
00106    /*----- Read command line -----*/
00107   if (argc < 2 || strcmp (argv[1], "-help") == 0)
00108     {
00109       printf ("Usage: 3dDWItoDT [options] gradient-file dataset\n"
00110               "Computes 6 principle direction tensors from multiple gradient vectors\n"
00111               " and corresponding DTI image volumes.\n"
00112               " The program takes two parameters as input :  \n"
00113               "    a 1D file of the gradient vectors with lines of ASCII floats Gxi,Gyi,Gzi.\n"
00114               "    Only the non-zero gradient vectors are included in this file (no G0 line).\n"
00115               "    a 3D bucket dataset with Np+1 sub-briks where the first sub-brik is the\n"
00116               "    volume acquired with no diffusion weighting.\n"
00117               " Options:\n"
00118               "   -prefix pname = Use 'pname' for the output dataset prefix name.\n"
00119               "    [default='DT']\n\n"
00120               "   -automask =  mask dataset so that the tensors are computed only for\n"
00121               "    high-intensity (presumably brain) voxels.  The intensity level is\n"
00122               "    determined the same way that 3dClipLevel works.\n\n"
00123               "   -mask dset = use dset as mask to include/exclude voxels\n\n"
00124               "   -nonlinear = compute iterative solution to avoid negative eigenvalues.\n"
00125               "    This is the default method.\n\n"
00126               "   -linear = compute simple linear solution.\n\n"
00127               "   -reweight = recompute weight factors at end of iterations and restart\n\n"
00128               "   -max_iter n = maximum number of iterations for convergence (Default=10).\n"
00129               "    Values can range from -1 to any positive integer less than 101.\n"
00130               "    A value of -1 is equivalent to the linear solution.\n"
00131               "    A value of 0 results in only the initial estimate of the diffusion tensor\n"
00132               "    solution adjusted to avoid negative eigenvalues.\n\n"
00133               "   -max_iter_rw n = max number of iterations after reweighting (Default=5)\n"
00134               "    values can range from 1 to any positive integer less than 101.\n\n"
00135               "   -eigs = compute eigenvalues, eigenvectors, fractional anisotropy and mean\n"
00136               "    diffusivity in sub-briks 6-19. Computed as in 3dDTeig\n\n"
00137               "   -debug_briks = add sub-briks with Ed (error functional), Ed0 (orig. error),\n"
00138               "     number of steps to convergence and I0 (modeled B0 volume)\n\n"
00139               "   -cumulative_wts = show overall weight factors for each gradient level\n"
00140               "    May be useful as a quality control\n\n"
00141               "   -verbose nnnnn = print convergence steps every nnnnn voxels that survive to\n"
00142               "    convergence loops (can be quite lengthy).\n\n"
00143               "   -drive_afni nnnnn = show convergence graphs every nnnnn voxels that survive\n"
00144               "    to convergence loops. AFNI must have NIML communications on (afni -niml).\n\n"
00145               " Example:\n"
00146               "  3dDWItoDT -prefix rw01 -automask -reweight -max_iter 10 \\\n"
00147               "            -max_iter_rw 10 tensor25.1D grad02+orig.\n\n"
00148               " The output is a 6 sub-brick bucket dataset containing Dxx,Dxy,Dxz,Dyy,Dyz,Dzz.\n"
00149               " Additional sub-briks may be appended with the -eigs and -debug_briks options.\n"
00150               " These results are appropriate as the input to the 3dDTeig program.\n"
00151               "\n");
00152       printf ("\n" MASTER_SHORTHELP_STRING);
00153       exit (0);
00154     }
00155 
00156   mainENTRY ("3dDWItoDT main");
00157   machdep ();
00158   AFNI_logger ("3dDWItoDT", argc, argv);
00159   PRINT_VERSION("3dDWItoDT") ;
00160 
00161   nopt = 1;
00162   nbriks = 6;           /* output contains 6 sub-briks by default */
00163   method = -1;
00164   reweight_flag = 0;
00165 
00166   datum = MRI_float;
00167   while (nopt < argc && argv[nopt][0] == '-')
00168     {
00169 
00170       /*-- prefix --*/
00171 
00172       if (strcmp (argv[nopt], "-prefix") == 0)
00173         {
00174           if (++nopt >= argc)
00175             {
00176               fprintf (stderr, "*** Error - prefix needs an argument!\n");
00177               exit (1);
00178             }
00179           MCW_strncpy (prefix, argv[nopt], THD_MAX_PREFIX);     /* change name from default prefix */
00180           /* check file name to be sure not to overwrite - mod drg 12/9/2004 */
00181           if (!THD_filename_ok (prefix))
00182             {
00183               fprintf (stderr, "*** Error - %s is not a valid prefix!\n", prefix);
00184               exit (1);
00185             }
00186           nopt++;
00187           continue;
00188         }
00189 
00190       /*-- datum --*/
00191 
00192       if (strcmp (argv[nopt], "-datum") == 0)
00193         {
00194           if (++nopt >= argc)
00195             {
00196               fprintf (stderr, "*** Error - datum needs an argument!\n");
00197               exit (1);
00198             }
00199           if (strcmp (argv[nopt], "short") == 0)
00200             {
00201               datum = MRI_short;
00202             }
00203           else if (strcmp (argv[nopt], "float") == 0)
00204             {
00205               datum = MRI_float;
00206             }
00207           else if (strcmp (argv[nopt], "byte") == 0)
00208             {
00209               datum = MRI_byte;
00210             }
00211           else
00212             {
00213               fprintf (stderr, "-datum of type '%s' is not supported!\n",
00214                        argv[nopt]);
00215               exit (1);
00216             }
00217           nopt++;
00218           continue;
00219         }
00220       if (strcmp (argv[nopt], "-automask") == 0)
00221         {
00222          if(maskptr != NULL){
00223            fprintf(stderr,"** ERROR: can't use -mask with -automask!\n");
00224            exit(1) ;
00225          }
00226           automask = 1;
00227           nopt++;
00228           continue;
00229         }
00230 
00231       if( strcmp(argv[nopt],"-mask") == 0 ){
00232          THD_3dim_dataset * mask_dset ;
00233          if( automask ){
00234            fprintf(stderr,"** ERROR: can't use -mask with -automask!\n");
00235            exit(1) ;
00236          }
00237          mask_dset = THD_open_dataset(argv[++nopt]) ;
00238          if( mask_dset == NULL ){
00239             fprintf(stderr,"** ERROR: can't open -mask dataset!\n"); exit(1);
00240          }
00241          if( maskptr != NULL ){
00242             fprintf(stderr,"** ERROR: can't have 2 -mask options!\n"); exit(1);
00243          }
00244          maskptr = THD_makemask( mask_dset , 0 , 1.0,-1.0 ) ;
00245          mmvox = DSET_NVOX( mask_dset ) ;
00246 
00247          DSET_delete(mask_dset) ; nopt++ ; continue ;
00248       }
00249 
00250       if (strcmp (argv[nopt], "-linear") == 0)
00251         {
00252           if(method==1)
00253             {
00254               fprintf(stderr, "*** Error - can not select both linear and non-linear methods at the same time\n");
00255               exit(1);
00256             }
00257           method = 0;
00258           nopt++;
00259           continue;
00260         }
00261 
00262       if ((strcmp (argv[nopt], "-nonlinear") == 0) || (strcmp (argv[nopt], "-non-linear") == 0))
00263         {
00264           if(method==0)
00265             {
00266               fprintf(stderr, "*** Error - can not select both linear and non-linear methods at the same time\n");
00267               exit(1);
00268             }
00269           method = 1;
00270           nopt++;
00271           continue;
00272         }
00273       if (strcmp (argv[nopt], "-reweight") == 0)
00274         {
00275           reweight_flag = 1;
00276           nopt++;
00277           continue;
00278         }
00279 
00280       if (strcmp (argv[nopt], "-max_iter") == 0)
00281         {
00282            if(++nopt >=argc ){
00283               fprintf(stderr,"*** Error - need an argument after -max_iter!\n");
00284               exit(1);
00285            }
00286            max_iter = strtol(argv[nopt], NULL, 10);
00287            if ((max_iter <-1)||(max_iter>100)) {
00288               fprintf(stderr, "Error - max_iter must be between -1 and 100\n");
00289               exit(1);
00290            }
00291           nopt++;
00292           continue;
00293         }
00294         
00295      if (strcmp (argv[nopt], "-max_iter_rw") == 0)
00296         {
00297            if(++nopt >=argc ){
00298               fprintf(stderr,"*** Error - need an argument after -max_iter_rw!\n");
00299               exit(1);
00300            }
00301            max_iter_rw = strtol(argv[nopt], NULL, 10);
00302            if ((max_iter_rw <=0)||(max_iter_rw>100)) {
00303               fprintf(stderr, "*** Error - max_iter_rw must be between 1 and 100\n");
00304               exit(1);
00305            }
00306           nopt++;
00307           continue;
00308         }
00309 
00310      if (strcmp (argv[nopt], "-eigs") == 0)
00311         {
00312           eigs_flag = 1;
00313           nopt++;
00314           continue;
00315         }
00316 
00317      if (strcmp (argv[nopt], "-cumulative_wts") == 0)
00318         {
00319           cumulative_flag = 1;
00320           nopt++;
00321           continue;
00322         }
00323 
00324      if (strcmp (argv[nopt], "-debug_briks") == 0)
00325         {
00326           debug_briks = 1;
00327           nopt++;
00328           continue;
00329         }
00330 
00331      if (strcmp (argv[nopt], "-verbose") == 0)
00332         {
00333            if(++nopt >=argc ){
00334               fprintf(stderr,"*** Error - need an argument after -verbose!\n");
00335               exit(1);
00336            }
00337            verbose = strtol(argv[nopt], NULL, 10);
00338            if (verbose<=0) {
00339               fprintf(stderr, "*** Error - verbose steps must be a positive number !\n");
00340               exit(1);
00341            }
00342           nopt++;
00343           continue;
00344         }
00345 
00346      if (strcmp (argv[nopt], "-drive_afni") == 0)
00347         {
00348            if(++nopt >=argc ){
00349               fprintf(stderr,"*** Error - need an argument after -drive_afni!\n");
00350               exit(1);
00351            }
00352            afnitalk_flag = strtol(argv[nopt], NULL, 10);
00353            if (afnitalk_flag<=0) {
00354               fprintf(stderr, "*** Error - drive_afni steps must be a positive number !\n");
00355               exit(1);
00356            }
00357           nopt++;
00358           continue;
00359         }
00360 
00361         fprintf(stderr, "*** Error - unknown option %s\n", argv[nopt]);
00362         exit(1);
00363     }
00364   
00365   if(method==-1)
00366       method = 1;        /* if not selected, choose non-linear method for now */
00367 
00368   if(max_iter>=-1){
00369      if(method==0)
00370               fprintf(stderr, "+++ Warning - max_iter will be ignored for linear methods.\n");
00371   }
00372   else
00373      max_iter = MAX_CONVERGE_STEPS;
00374 
00375   if(max_iter_rw>=0) {
00376      if(method==0)
00377               fprintf(stderr, "+++ Warning - max_iter_rw will be ignored for linear methods.\n");
00378      if(reweight_flag==0)
00379               fprintf(stderr, "+++ Warning - max_iter_rw will be ignored when not reweighting.\n");
00380   }
00381   else
00382      max_iter_rw = MAX_RWCONVERGE_STEPS;
00383      
00384   if((method==0)&&(reweight_flag==1)) {
00385       fprintf(stderr, "+++ Warning - can not reweight voxels for linear method\n");
00386       reweight_flag = 0;
00387   }
00388 
00389   if(cumulative_flag==1) {
00390      if(method==0) {
00391         fprintf(stderr, "+++ Warning - can not compute cumulative weights for linear method\n");
00392         cumulative_flag = 0;
00393      }
00394      if(reweight_flag == 0) {
00395         fprintf(stderr, "+++ Warning - can not compute cumulative weights if not reweighting\n");
00396         cumulative_flag = 0;
00397      }
00398   }
00399 
00400   if((method==0)&&(debug_briks==1)) {
00401       fprintf(stderr, "+++ Warning - can not compute debug sub-briks for linear method\n");
00402       debug_briks = 0;
00403   }
00404 
00405   if((method==0)&&(afnitalk_flag>0)) {
00406       fprintf(stderr, "+++ Warning - can not graph convergence in AFNI for linear method\n");
00407       afnitalk_flag = 0;
00408   }
00409 
00410   if(eigs_flag)
00411      nbriks += 14;
00412 
00413   if(debug_briks)
00414      nbriks += 4;
00415      
00416 
00417    /*----- read input datasets -----*/
00418 
00419   if (nopt >= argc)
00420     {
00421       fprintf (stderr, "*** Error - No input dataset!?\n");
00422       exit (1);
00423     }
00424 
00425   /* first input dataset - should be gradient vector file of ascii floats Gx,Gy,Gz */
00426 
00427   /* read gradient vector 1D file */
00428   grad1Dptr = mri_read_1D (argv[nopt]);
00429   if (grad1Dptr == NULL)
00430     {
00431       fprintf (stderr, "*** Error reading gradient vector file\n");
00432       exit (1);
00433     }
00434 
00435   if (grad1Dptr->ny != 3)
00436     {
00437       fprintf (stderr, "*** Error - Only 3 columns of gradient vectors allowed\n");
00438       fprintf (stderr, "%d columns found\n", grad1Dptr->nx);
00439       mri_free (grad1Dptr);
00440       exit (1);
00441     }
00442 
00443   if (grad1Dptr->nx < 6)
00444     {
00445       fprintf (stderr, "*** Error - Must have at least 6 gradient vectors\n");
00446       fprintf (stderr, "%d columns found\n", grad1Dptr->nx);
00447       mri_free (grad1Dptr);
00448       exit (1);
00449     }
00450 
00451 
00452   Form_R_Matrix (grad1Dptr);    /* use grad1Dptr to compute R matrix */
00453 
00454   nopt++;
00455 
00456   /* Now read in all the MRI volumes for each gradient vector */
00457   /* assumes first one is no gradient */
00458   old_dset = THD_open_dataset (argv[nopt]);
00459 
00460   if (!ISVALID_DSET (old_dset))
00461     {
00462       fprintf (stderr, "*** Error - Can not open dataset %s\n", argv[nopt]);
00463       mri_free (grad1Dptr);
00464       exit (1);
00465     }
00466 
00467   /* expect at least 7 values per voxel - 7 sub-briks as input dataset */
00468   if (DSET_NVALS (old_dset) != (grad1Dptr->nx + 1))
00469     {
00470       fprintf (stderr,
00471                "*** Error - Dataset must have number of sub-briks equal to one more than number\n");
00472       fprintf (stderr, "  of gradient vectors (B0+Bi)!\n");
00473       mri_free (grad1Dptr);
00474       exit (1);
00475     }
00476    nxyz = DSET_NVOX(old_dset) ;
00477    if( maskptr != NULL && mmvox != nxyz ){
00478       fprintf(stderr,"** Mask and input datasets not the same size!\n") ;
00479       exit(1) ;
00480    }
00481 
00482 
00483   InitGlobals (grad1Dptr->nx + 1);      /* initialize all the matrices and vectors */
00484   Computebmatrix (grad1Dptr);   /* compute bij=GiGj */
00485 
00486   if (automask)
00487     {
00488       DSET_mallocize (old_dset);
00489       DSET_load (old_dset);     /* get B0 (anatomical image) from dataset */
00490       /*anat_im = THD_extract_float_brick( 0, old_dset ); */
00491       anat_im = DSET_BRICK (old_dset, 0);       /* set pointer to the 0th sub-brik of the dataset */
00492       maskptr = mri_automask_image (anat_im);   /* maskptr is a byte pointer for volume */
00493 #if 0
00494       /* convert byte mask to same format type as dataset */
00495       nvox = DSET_NVOX (old_dset);
00496       sar = (short *) calloc (nvox, sizeof (short));
00497       /* copy maskptr values to far ptr */
00498       tempsptr = sar;
00499       tempbptr = maskptr;
00500       for (i = 0; i < nvox; i++)
00501         {
00502           *tempsptr++ = (short) *tempbptr++;
00503           tempval = *(tempsptr - 1);
00504         }
00505 
00506       free (maskptr);
00507 
00508       /*old_dset->dblk->malloc_type = DATABLOCK_MEM_MALLOC; *//* had to set this? */
00509       EDIT_add_brick (old_dset, MRI_short, 0.0, sar);   /* add sub-brik to end */
00510 
00511 #endif
00512     }
00513 
00514   /* temporarily set artificial timing to 1 second interval */
00515   EDIT_dset_items (old_dset,
00516                    ADN_ntt, DSET_NVALS (old_dset),
00517                    ADN_ttorg, 0.0,
00518                    ADN_ttdel, 1.0, ADN_tunits, UNITS_SEC_TYPE, NULL);
00519 
00520   if(afnitalk_flag) {
00521      if(DWI_Open_NIML_stream()!=0) {   /* Open NIML stream */
00522        afnitalk_flag = 0;
00523        fprintf(stderr,"+++could not open NIML communications with AFNI\n");
00524      }
00525      else
00526        if(DWI_NIML_create_graph()!=0) {
00527           afnitalk_flag = 0;
00528           fprintf(stderr,"+++could not create graph within AFNI\n");
00529           /* Close NIML stream */
00530           NI_stream_close(DWIstreamid);
00531           DWIstreamid = 0;
00532        }
00533   }
00534 
00535    /*------------- ready to compute new dataset -----------*/
00536 
00537   new_dset = MAKER_4D_to_typed_fbuc (old_dset,  /* input dataset */
00538                                      prefix,    /* output prefix */
00539                                      datum,     /* output datum  */
00540                                      0, /* ignore count  */
00541                                      0, /* can't detrend in maker function  KRH 12/02 */
00542                                      nbriks,    /* number of briks */
00543                                      DWItoDT_tsfunc,    /* timeseries processor */
00544                                      NULL       /* data for tsfunc */
00545     );
00546 
00547 
00548   if(afnitalk_flag && (DWIstreamid!=0)) {
00549 /* Close NIML stream */
00550     NI_stream_close(DWIstreamid);
00551   }
00552 
00553   if(cumulative_flag && reweight_flag) {
00554     cumulativewtptr = cumulativewt;
00555     printf("Cumulative Wt. factors: ");
00556     for(i=0;i<(grad1Dptr->nx + 1);i++){
00557        *cumulativewtptr = *cumulativewtptr / rewtvoxels;
00558        printf("%5.3f ", *cumulativewtptr++);
00559     }
00560      printf("\n");
00561   }
00562 
00563   FreeGlobals ();
00564   mri_free (grad1Dptr);
00565   matrix_destroy (&Rtmat);      /* clean up */
00566 
00567   if (maskptr)
00568     {
00569       free (maskptr);
00570       if(anat_im)
00571          mri_free (anat_im);
00572 #if 0
00573       DSET_unload_one (old_dset, 0);
00574       sar = NULL;
00575 #endif
00576     }
00577 
00578   if (new_dset != NULL)
00579     {
00580       tross_Copy_History (old_dset, new_dset);
00581       EDIT_dset_items (new_dset, ADN_brick_label_one + 0, "Dxx", ADN_none);
00582       EDIT_dset_items (new_dset, ADN_brick_label_one + 1, "Dxy", ADN_none);
00583       EDIT_dset_items (new_dset, ADN_brick_label_one + 2, "Dxz", ADN_none);
00584       EDIT_dset_items (new_dset, ADN_brick_label_one + 3, "Dyy", ADN_none);
00585       EDIT_dset_items (new_dset, ADN_brick_label_one + 4, "Dyz", ADN_none);
00586       EDIT_dset_items (new_dset, ADN_brick_label_one + 5, "Dzz", ADN_none);
00587       if(eigs_flag) {
00588         eigs_brik = ADN_brick_label_one + 6;   /* 1st eigenvalue brik */
00589         EDIT_dset_items(new_dset, eigs_brik+0, "lambda_1", ADN_none);
00590         EDIT_dset_items(new_dset, eigs_brik+1, "lambda_2",ADN_none);
00591         EDIT_dset_items(new_dset, eigs_brik+2, "lambda_3",ADN_none);
00592         EDIT_dset_items(new_dset, eigs_brik+3, "eigvec_1[1]",ADN_none);
00593         EDIT_dset_items(new_dset, eigs_brik+4, "eigvec_1[2]",ADN_none);
00594         EDIT_dset_items(new_dset, eigs_brik+5, "eigvec_1[3]",ADN_none);
00595         EDIT_dset_items(new_dset, eigs_brik+6, "eigvec_2[1]",ADN_none);
00596         EDIT_dset_items(new_dset, eigs_brik+7, "eigvec_2[2]",ADN_none);
00597         EDIT_dset_items(new_dset, eigs_brik+8, "eigvec_2[3]",ADN_none);
00598         EDIT_dset_items(new_dset, eigs_brik+9, "eigvec_3[1]",ADN_none);
00599         EDIT_dset_items(new_dset, eigs_brik+10,"eigvec_3[2]",ADN_none);
00600         EDIT_dset_items(new_dset, eigs_brik+11,"eigvec_3[3]",ADN_none);
00601         EDIT_dset_items(new_dset, eigs_brik+12,"FA",ADN_none);
00602         EDIT_dset_items(new_dset, eigs_brik+13,"MD",ADN_none);
00603       }
00604 
00605       if(debug_briks) {
00606         EDIT_dset_items (new_dset, ADN_brick_label_one + nbriks - 4, "Converge Step", ADN_none);
00607         EDIT_dset_items (new_dset, ADN_brick_label_one + nbriks - 3, "ED", ADN_none);
00608         EDIT_dset_items (new_dset, ADN_brick_label_one + nbriks - 2, "EDorig", ADN_none);
00609         EDIT_dset_items (new_dset, ADN_brick_label_one + nbriks - 1, "I0", ADN_none);
00610       }
00611 
00612       tross_Make_History ("3dDWItoDT", argc, argv, new_dset);
00613       DSET_write (new_dset);
00614       fprintf(stderr,"--- Output dataset %s\n", DSET_BRIKNAME(new_dset));
00615     }
00616   else
00617     {
00618       fprintf (stderr, "*** Error - Unable to compute output dataset!\n");
00619       exit (1);
00620     }
00621 
00622   exit (0);
00623 }
00624 
00625 /* Form R matrix as matrix of [bxx 2bxy 2bxz byy 2byz bzz] for Np rows */
00626 static void
00627 Form_R_Matrix (MRI_IMAGE * grad1Dptr)
00628 {
00629   matrix Rmat;
00630   register double sf = 1.0;     /* scale factor = 1.0 for now until we know DELTA, delta (gamma = 267.5 rad/ms-mT) */
00631   register double sf2;          /* just scale factor * 2 */
00632   int i, nrows;
00633   register float *imptr, *Gxptr, *Gyptr, *Gzptr;
00634   matrix *nullptr = NULL;
00635   register double Gx, Gy, Gz;
00636 
00637   ENTRY ("Form_R_Matrix");
00638   nrows = grad1Dptr->nx;
00639   matrix_initialize (&Rmat);
00640   matrix_create (nrows, 6, &Rmat);      /* Rmat = Np x 6 matrix */
00641   if (Rmat.elts == NULL)
00642     {                           /* memory allocation error */
00643       fprintf (stderr, "could not allocate memory for Rmat \n");
00644       EXRETURN;
00645     }
00646   sf2 = sf + sf;                /* 2 * scale factor for minor speed improvement */
00647   Gxptr = imptr = MRI_FLOAT_PTR (grad1Dptr);    /* use simple floating point pointers to get values */
00648   Gyptr = imptr + nrows;
00649   Gzptr = Gyptr + nrows;
00650 
00651   for (i = 0; i < nrows; i++)
00652     {
00653       Gx = *Gxptr++;
00654       Gy = *Gyptr++;
00655       Gz = *Gzptr++;
00656       Rmat.elts[i][0] = sf * Gx * Gx;   /* bxx = Gx*Gx*scalefactor */
00657       Rmat.elts[i][1] = sf2 * Gx * Gy;  /* 2bxy = 2GxGy*scalefactor */
00658       Rmat.elts[i][2] = sf2 * Gx * Gz;  /* 2bxz = 2GxGz*scalefactor */
00659       Rmat.elts[i][3] = sf * Gy * Gy;   /* byy = Gy*Gy*scalefactor */
00660       Rmat.elts[i][4] = sf2 * Gy * Gz;  /* 2byz = 2GyGz*scalefactor */
00661       Rmat.elts[i][5] = sf * Gz * Gz;   /* bzz = Gz*Gz*scalefactor */
00662     }
00663 
00664   matrix_initialize (&Rtmat);
00665   matrix_psinv (Rmat, nullptr, &Rtmat); /* compute pseudo-inverse of Rmat=Rtmat */
00666   matrix_destroy (&Rmat);       /*  from the other two matrices */
00667   EXRETURN;
00668 }
00669 
00670 
00671 /**********************************************************************
00672    Function that does the real work
00673 ***********************************************************************/
00674 
00675 static void
00676 DWItoDT_tsfunc (double tzero, double tdelta,
00677                 int npts, float ts[],
00678                 double ts_mean, double ts_slope,
00679                 void *ud, int nbriks, float *val)
00680 {
00681   int i, converge_step, converge, trialstep, ntrial, adjuststep, recordflag;
00682   double orig_deltatau, EDold, J, dt;
00683   static int nvox, ncall, noisecall;
00684   register double i0;
00685   register double dv, dv0;
00686   vector lnvector;
00687   int wtflag;          /* wtflag for recomputing wtfactor*/
00688   int max_converge_step, graphpoint;
00689   double dtau[50], Edgraph[50];
00690   int graphflag;
00691 
00692   ENTRY ("DWItoDT_tsfunc");
00693   /* ts is input vector data of Np+1 floating point numbers.
00694      For each point in volume brik convert vector data to
00695      symmetric matrix */
00696   /* ts should come from data sub-briks in form of I0,I1,...Ip */
00697   /* val is output vector of form Dxx Dxy Dxz Dyy Dyz Dzz for each voxel in 6 sub-briks */
00698   /* the Dij vector is computed as the product of  Rt times ln(I0/Ip)
00699      where Rt is the pseudo-inverse of the [bxx 2bxy 2bxz byy 2byz bzz] for
00700      each gradient vector b */
00701    /** is this a "notification"? **/
00702    if (val == NULL)
00703     {
00704 
00705       if (npts > 0)
00706         {                       /* the "start notification" */
00707           nvox = npts;          /* keep track of   */
00708           ncall = 0;            /* number of calls */
00709           noisecall = 0;
00710         }
00711       else
00712         {                       /* the "end notification" */
00713 
00714           /* nothing to do here */
00715         }
00716       EXRETURN;
00717     }
00718 
00719   ncall++;
00720   /* if there is any mask (automask or user mask), use corresponding voxel as a flag */
00721   if (maskptr)
00722     {
00723 #if 0
00724      npts = npts - 1;
00725      if (ts[npts] == 0)
00726 #endif
00727        if(maskptr[ncall-1]==0)
00728         {                       /* don't include this voxel for mask */
00729           for (i = 0; i < nbriks; i++)  /* faster to copy preset vector */
00730             val[i] = 0.0;       /* return 0 for all Dxx,Dxy,... */
00731           if(debug_briks)  /* use -3 as flag for number of converge steps to mean exited for masked voxels */
00732              val[nbriks-4] = -3.0;
00733           EXRETURN;
00734         }
00735     }
00736   /* load the symmetric matrix vector from the "timeseries" subbrik vector values */
00737   vector_initialize (&lnvector);
00738   vector_create_noinit (npts - 1, &lnvector);
00739   dv0 = ts[0];
00740   if (dv0 > 0.0)
00741     i0 = log (dv0);
00742   else
00743     i0 = 0.0;
00744   for (i = 0; i < (npts - 1); i++)
00745     {
00746       dv = ts[i + 1];
00747       if ((dv > 0.0) && (dv0 > 0.0))
00748         lnvector.elts[i] = i0 - log (dv);       /* ln I0/Ip = ln I0 - ln Ip */
00749       else
00750         lnvector.elts[i] = 0.0;
00751     }
00752 
00753   vector_multiply (Rtmat, lnvector, &Dvector);  /* D = Rt * ln(I0/Ip), allocated Dvector here */
00754 
00755   vector_destroy (&lnvector);   /* free vector elements allocated */
00756 
00757   if((method==0)||(max_iter==-1)) {     /* for linear method,stop here and return D values */
00758      vector_to_array(Dvector, val);
00759 
00760      if(debug_briks) {
00761        InitWtfactors (npts);            /* initialize all weight factors to 1 for all gradient intensities */
00762        J = ComputeJ (ts, npts); /* Ed (error) computed here */
00763        val[nbriks-4] = -1.0;  /* use -1 as flag for number of converge steps to mean exited for */
00764                               /* or initial insignificant deltatau value */
00765        val[nbriks-3] = ED;
00766        val[nbriks-2] = ED;                  /* store original error */
00767        val[nbriks-1] = J;
00768      }
00769 
00770      if(eigs_flag) {                            /* if user wants eigenvalues in output dataset */
00771         EIG_func();                              /* calculate eigenvalues, eigenvectors here */
00772         for(i=0;i<12;i++) 
00773           val[i+6] = eigs[i];
00774        /* calc FA */
00775        val[18] = Calc_FA(val+6);                /* calculate fractional anisotropy */
00776        val[19] = Calc_MD(val+6);                /* calculate mean diffusivity */
00777      }
00778      EXRETURN;
00779   }
00780  
00781   /* now more complex part that takes into account noise */
00782 
00783   /* calculate initial estimate of D using standard linear model */
00784   EIG_func ();                  /* compute eigenvalues, eigenvectors standard way */
00785   
00786 
00787   InitWtfactors (npts);         /* initialize all weight factors to 1 for all gradient intensities */
00788 
00789   ComputeD0 ();                 /* recalculate Dmatrix based on limits on eigenvalues */
00790 
00791 
00792   if(matrix_sumabs(Dmatrix)<=TINYNUMBER) {
00793     for(i=0;i<nbriks;i++)
00794       val[i] = 0.0;
00795      if(debug_briks) {
00796        val[nbriks-4] = -2.0; /* use -2 as flag for number of converge steps to mean exited for insignificant D values*/
00797        val[nbriks-3] = 0;
00798        val[nbriks-2] = 0;                  /* store original error */
00799        val[nbriks-1] = 0;
00800      }
00801 
00802     EXRETURN;
00803   }
00804 
00805 
00806   converge_step = 0;    /* allow up to max_iter=MAX_CONVERGE_STEPS (10) deltatau steps */
00807   max_converge_step = max_iter;   /* 1st time through set limit of converge steps to user option */
00808   converge = 0;
00809   wtflag = reweight_flag;
00810 
00811   /* trial step */
00812   J = ComputeJ (ts, npts);      /* Ed (error) computed here */
00813   Store_Computations (0, npts, wtflag); /* Store 1st adjusted computations */
00814   matrix_copy (Dmatrix, &OldD);   /* store first Dmatrix in OldD too */
00815 
00816   if(debug_briks)
00817      val[nbriks-2] = ED;                  /* store original error */
00818 
00819   EDold = ED;
00820   ComputeDeltaTau ();
00821   if(deltatau<=TINYNUMBER) {         /* deltatau too small, exit */
00822     for(i=0;i<nbriks;i++)
00823       val[i] = 0.0;
00824     if(debug_briks) {
00825       val[nbriks-4] = -1.0; /* use -1 as flag for number of converge steps to mean exited for insignificant deltatau value */
00826       val[nbriks-1] = J;
00827     }
00828     EXRETURN;
00829   }
00830 
00831   if(verbose&&(!(noisecall%verbose)))   /* show verbose messages every verbose=n voxels */
00832      recordflag = 1;
00833      else
00834      recordflag = 0;
00835 
00836   if(afnitalk_flag&&(!(noisecall%afnitalk_flag))) {  /* graph in AFNI convergence steps every afnitalk_flag=n voxels */
00837      graphflag = 1;
00838      graphpoint = 0;
00839    }
00840   else
00841      graphflag = 0;
00842 
00843   noisecall++;
00844 
00845       ntrial = 0;
00846 
00847       while ((converge_step < max_converge_step) && (converge!=1) && (ntrial < 10) )
00848         {
00849       /* find trial step */
00850       /* first do trial step to find acceptable delta tau */
00851       /* Take half of previous tau step to see if less error */
00852       /* if there is less error than delta tau is acceptable */
00853       /* try halving up to 10 times, if it does not work, */
00854       /* use first D from initial estimate */
00855         trialstep = 1;
00856         ntrial = 0;
00857         while ((trialstep) && (ntrial < 10))
00858           {                     /* allow up to 10 iterations to find trial step */
00859             ComputeHpHm (deltatau);
00860             ComputeNewD ();
00861             J = ComputeJ (ts, npts);    /* Ed (error) computed here */
00862             if (ED < EDold)          /* is error less than error of trial step or previous step? */
00863                 {
00864                 /* found acceptable step size of DeltaTau */
00865                 if(recordflag==1)
00866                  printf("ncall= %d, converge_step=%d, deltatau=%f, ED=%f, ntrial %d in find dtau\n", ncall, converge_step, deltatau, ED, ntrial);
00867                 if(graphflag==1) {
00868                   dtau[graphpoint] = deltatau;
00869                   Edgraph[graphpoint] = ED;
00870                   graphpoint++;
00871                 }
00872                                      
00873                 EDold = ED;
00874                 trialstep = 0;
00875                 Store_Computations (1, npts, wtflag);   /* Store current computations */
00876                 }
00877             else {
00878                 Restore_Computations (0, npts, wtflag); /* Restore trial 0 computations */
00879                 dt = deltatau / 2;
00880                 deltatau = dt;  /* find first DeltaTau step with less error than 1st guess */
00881                 /* by trying smaller step sizes */
00882                 ntrial++;
00883               }
00884           }
00885 
00886 
00887         /* end of finding trial step size */
00888         /* in trial step stage, already have result of deltatau step and may have
00889            already tried deltatau*2 if halved (ntrial>=1) */
00890         if(ntrial <10) {
00891           orig_deltatau = deltatau;
00892           adjuststep = 1;
00893 
00894           for(i=0;i<2;i++) {
00895             if(i==0)
00896                deltatau = orig_deltatau*2;
00897             else
00898                deltatau = orig_deltatau/2;
00899                
00900             if((adjuststep==1) && ((i!=0) || (ntrial<2))) {   /* if didn't shrink in initial deltatau step above */
00901               Restore_Computations (0, npts, wtflag);   /* Restore previous Tau step computations */
00902               ComputeHpHm (deltatau);
00903               ComputeNewD ();
00904               J = ComputeJ (ts, npts);  /* computes Intensity without noise,*/
00905                                         /*   Ed, Residuals */
00906               if(ED<EDold){
00907                 adjuststep = 0;
00908                 Store_Computations(0, npts, wtflag);    /* Store Tau+dtau step computations */
00909                 EDold = ED;
00910 
00911                 if(recordflag==1)
00912                   printf("ncall= %d, converge_step=%d, deltatau=%f, ED=%f dt*2 best\n", ncall, converge_step, deltatau, ED);
00913 
00914                 if(graphflag==1) {
00915                   dtau[graphpoint] = deltatau;
00916                   Edgraph[graphpoint] = ED;
00917                   graphpoint++;
00918                 }
00919               }
00920             }
00921           }
00922 
00923           if(adjuststep!=0){            /* best choice was first Delta Tau */
00924              ED = EDold;
00925              deltatau = orig_deltatau;
00926              Restore_Computations (1,  npts, wtflag);   /* restore old computed matrices*/
00927                                            /*   D,Hp,Hm,F,R */
00928              Store_Computations(0, npts, wtflag);       /* Store Tau+dtau step computations */
00929           }
00930          if(recordflag==1)
00931               printf("ncall= %d, converge_step=%d, deltatau=%f, ED=%f dt best\n", ncall, converge_step, deltatau, ED);
00932 
00933          if(graphflag==1) {
00934             dtau[graphpoint] = deltatau;
00935             Edgraph[graphpoint] = ED;
00936             graphpoint++;
00937          }
00938 
00939          if (converge_step != 0) {      /* first time through recalculate*/
00940              /* now see if converged yet */
00941              converge = TestConvergence(Dmatrix, OldD);
00942           }
00943 
00944           matrix_copy (Dmatrix, &OldD);
00945           if(recordflag==1)
00946               printf("ncall= %d, converge_step=%d, deltatau=%f, ED=%f\n", ncall, converge_step, deltatau, ED);
00947 
00948          if(graphflag==1) {
00949             dtau[graphpoint] = deltatau;
00950             Edgraph[graphpoint] = ED;
00951             graphpoint++;
00952          }
00953 
00954           converge_step++;
00955         }
00956         else
00957           {
00958             if(recordflag==1)
00959              printf("ncall= %d, converge_step=%d, deltatau=%f, ED=%f Exiting time evolution\n", ncall, converge_step, deltatau, ED);
00960             Restore_Computations(0, npts, wtflag);       /* Exit with original step calculation */
00961             ED = EDold;
00962           }
00963 
00964 
00965         if(((converge) || (converge_step==max_iter)) && wtflag && reweight_flag) {  /* if at end of iterations the first time*/
00966              converge = 0;                  /* through whole group of iterations */
00967              ComputeWtfactors (npts);       /* compute new weight factors */
00968              converge_step = 1;             /* start over now with computed weight factors */
00969              max_converge_step = max_iter_rw+1;   /* reset limit of converge steps to user option */
00970              wtflag = 0;                    /* only do it once - turn off next reweighting */
00971              J=ComputeJ(ts, npts);            /* compute new Ed value */
00972              EDold = ED;                    /* this avoids having to go through converge loop for two loops */
00973           }
00974         }  /* end while converge loop */
00975 
00976           ED = EDold;     
00977 
00978   val[0] = Dmatrix.elts[0][0];
00979   val[1] = Dmatrix.elts[0][1];
00980   val[2] = Dmatrix.elts[0][2];
00981   val[3] = Dmatrix.elts[1][1];
00982   val[4] = Dmatrix.elts[1][2];
00983   val[5] = Dmatrix.elts[2][2];
00984 
00985   if(eigs_flag) {                            /* if user wants eigenvalues in output dataset */
00986     udmatrix_to_vector(Dmatrix, &Dvector);
00987     EIG_func();                              /* calculate eigenvalues, eigenvectors here */
00988     for(i=0;i<12;i++)
00989       val[i+6] = eigs[i];
00990     /* calc FA */
00991     val[18] = Calc_FA(val+6);                /* calculate fractional anisotropy */
00992     val[19] = Calc_MD(val+6);               /* calculate average (mean) diffusivity */
00993   }
00994 
00995   /* testing information only */
00996   if(recordflag)
00997      printf("ncall= %d, converge_step=%d, deltatau=%f, ED=%f\n", ncall, converge_step, deltatau, ED);
00998   if(debug_briks) {
00999     val[nbriks-4] = converge_step;
01000     val[nbriks-3] = ED;
01001     val[nbriks-1] = ComputeJ(ts, npts);            /* compute J value */;
01002   }
01003   if(graphflag==1) {
01004      DWI_AFNI_update_graph(Edgraph, dtau, graphpoint);
01005   }
01006 
01007   EXRETURN;
01008 }
01009 
01010 /* taken from 3dDTeig.c */
01011 /*! given Dij tensor data for principal directions */
01012 /* calculate 3 principal sets of eigenvalues and eigenvectors */
01013 static void
01014 EIG_func ()
01015 {
01016   /*  THD_dmat33 inmat;
01017      THD_dvecmat eigvmat; */
01018   int i, j;
01019   int maxindex, minindex, midindex;
01020   float temp, minvalue, maxvalue;
01021   int sortvector[3];
01022   double a[9], e[3];
01023   int astart, vstart;
01024 
01025 
01026   ENTRY ("EIG_func");
01027   /* Dvector is vector data of 6 floating point numbers.
01028      For each point in volume brik convert vector data to
01029      symmetric matrix */
01030   /* Dvector should come in form of Dxx,Dxy,Dxz,Dyy,Dyz,Dzz */
01031   /* convert to matrix of form 
01032      [ Dxx Dxy Dxz]
01033      [ Dxy Dyy Dyz]
01034      [ Dxz Dyz Dzz]  */
01035 
01036 
01037   /* load the symmetric matrix vector from the "timeseries" subbrik vector values */
01038 
01039   a[0] = Dvector.elts[0];
01040   a[1] = Dvector.elts[1];
01041   a[2] = Dvector.elts[2];
01042   a[3] = Dvector.elts[1];
01043   a[4] = Dvector.elts[3];
01044   a[5] = Dvector.elts[4];
01045   a[6] = Dvector.elts[2];
01046   a[7] = Dvector.elts[4];
01047   a[8] = Dvector.elts[5];
01048 
01049   symeig_double (3, a, e);      /* compute eigenvalues in e, eigenvectors in a */
01050 
01051   maxindex = 2;                 /* find the lowest, middle and highest eigenvalue */
01052   maxvalue = e[2];
01053   minindex = 0;
01054   minvalue = e[0];
01055   midindex = 1;
01056 
01057   for (i = 0; i < 3; i++)
01058     {
01059       temp = e[i];
01060       if (temp > maxvalue)
01061         {                       /* find the maximum */
01062           maxindex = i;
01063           maxvalue = temp;
01064         }
01065       if (temp < minvalue)
01066         {                       /* find the minimum */
01067           minindex = i;
01068           minvalue = temp;
01069         }
01070     }
01071 
01072   for (i = 0; i < 3; i++)
01073     {                           /* find the middle */
01074       if ((i != maxindex) && (i != minindex))
01075         {
01076           midindex = i;
01077           break;
01078         }
01079     }
01080 
01081   sortvector[0] = maxindex;
01082   sortvector[1] = midindex;
01083   sortvector[2] = minindex;
01084 
01085   /* put the eigenvalues at the beginning of the matrix */
01086   for (i = 0; i < 3; i++)
01087     {
01088       eigs[i] = e[sortvector[i]];       /* copy sorted eigenvalues */
01089 
01090       /* start filling in eigenvector values */
01091       astart = sortvector[i] * 3;       /* start index of double eigenvector */
01092       vstart = (i + 1) * 3;     /* start index of float val vector to copy eigenvector */
01093 
01094       for (j = 0; j < 3; j++)
01095         {
01096           eigs[vstart + j] = a[astart + j];
01097         }
01098     }
01099 
01100   EXRETURN;
01101 }
01102 
01103 /*! calculate Fractional Anisotropy */
01104 /* passed float pointer to start of eigenvalues */
01105 static float
01106 Calc_FA(float *val)
01107 {
01108   float FA;
01109   double ssq, dv0, dv1, dv2, dsq;
01110 
01111   ENTRY("Calc_FA");
01112 
01113   /* calculate the Fractional Anisotropy, FA */
01114   /*   reference, Pierpaoli C, Basser PJ. Microstructural and physiological features 
01115        of tissues elucidated by quantitative-diffusion tensor MRI,J Magn Reson B 1996; 111:209-19 */
01116   if((val[0]<=0.0)||(val[1]<=0.0)||(val[2]<=0.0)) {   /* any negative eigenvalues*/
01117     RETURN(0.0);                                      /* should not see any for non-linear method. Set FA to 0 */  
01118   }
01119 
01120   ssq = (val[0]*val[0])+(val[1]*val[1])+(val[2]*val[2]);        /* sum of squares of eigenvalues */
01121   /* 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 */
01122 
01123   dv0 = val[0]-val[1];
01124   dv0 *= dv0;
01125   dv1 = val[1]-val[2];
01126   dv1 *= dv1;
01127   dv2 = val[2]-val[0];
01128   dv2 *= dv2;
01129   dsq = dv0+dv1+dv2;                 /* sum of differences squared */
01130 
01131   if(ssq!=0.0)
01132     FA = sqrt(dsq/(2.0*ssq));   /* FA calculated here */
01133   else
01134     FA = 0.0;
01135 
01136   RETURN(FA);
01137 }
01138 
01139 /*! calculate Mean Diffusivity */
01140 /* passed float pointer to start of eigenvalues */
01141 static float
01142 Calc_MD(float *val)
01143 {
01144   float MD;
01145 
01146   ENTRY("Calc_MD");
01147 
01148   /* calculate the Fractional Anisotropy, FA */
01149   /*   reference, Pierpaoli C, Basser PJ. Microstructural and physiological features 
01150        of tissues elucidated by quantitative-diffusion tensor MRI,J Magn Reson B 1996; 111:209-19 */
01151   if((val[0]<=0.0)||(val[1]<=0.0)||(val[2]<=0.0)) {   /* any negative eigenvalues*/
01152     RETURN(0.0);                                      /* should not see any for non-linear method. Set FA to 0 */  
01153   }
01154   MD = (val[0] + val[1] + val[2]) / 3;
01155 
01156   RETURN(MD);
01157 }
01158 
01159 
01160 /*! compute initial estimate of D0 */
01161 /* D = estimated diffusion tensor matrix 
01162       [Dxx Dxy Dxz, Dxy Dyy Dyz, Dxz Dyz Dzz] */
01163 /* updates Dvector and Dmatrix */
01164 static void
01165 ComputeD0 ()
01166 {
01167   int i, j;
01168   /*   matrix ULmatrix, Ematrix; */
01169   double mu, alpha, sum;
01170   double e10, e11, e12, e20, e21, e22, e30, e31, e32;
01171   double l1, l2, l3;
01172   double t1, t3, t5, t8, t10, t12, t14, t18, t19, t21, t23, t32, t33, t35,
01173     t37;
01174 
01175   ENTRY ("ComputeD0");
01176   /* create and initialize D0 */
01177 
01178   if (eigs[0] < 0)
01179     {                           /* if all eigenvalues are negative - may never happen */
01180       /* D0 = diag(a,a,a) where a=1/3 Sum(Abs(Lambda_i)) */
01181       sum = 0.0;
01182       for (i = 0; i < 3; i++)
01183         sum += fabs (eigs[i]);
01184       alpha = sum / 3;
01185       for (i = 0; i < 3; i++)
01186         {
01187           for (j = 0; j < 3; j++)
01188             {
01189               if (i == j)
01190                 Dmatrix.elts[i][j] = alpha;
01191               else
01192                 Dmatrix.elts[i][j] = 0.0;
01193             }
01194         }
01195 
01196       udmatrix_to_vector (Dmatrix, &Dvector);   /* convert to vector format for D also */
01197       EXRETURN;
01198     }
01199 
01200   mu = 0.2 * eigs[0];
01201   if (eigs[1] < mu)             /* set limit of eigenvalues to prevent */
01202      eigs[1] = mu;              /*     too much anisotropy */
01203   if (eigs[2] < mu)
01204      eigs[2] = mu;
01205 
01206   /* D0 = U L UT */
01207 /*
01208                             [e10 l1    e20 l2    e30 l3]
01209                             [                          ]
01210                       UL := [e11 l1    e21 l2    e31 l3]
01211                             [                          ]
01212                             [e12 l1    e22 l2    e32 l3]
01213 */
01214 
01215 
01216   /* assign variables to match Maple code */
01217   l1 = eigs[0];
01218   l2 = eigs[1];
01219   l3 = eigs[2];
01220 
01221   e10 = eigs[3];
01222   e11 = eigs[4];
01223   e12 = eigs[5];
01224 
01225   e20 = eigs[6];
01226   e21 = eigs[7];
01227   e22 = eigs[8];
01228 
01229   e30 = eigs[9];
01230   e31 = eigs[10];
01231   e32 = eigs[11];
01232 
01233 
01234 #ifdef lkjsaklfj
01235   matrix_initialize (&Ematrix);
01236   matrix_create (3, 3, &Ematrix);
01237   /* fill Ematrix with Eigenvectors */
01238 
01239   matrix_initialize (&ULmatrix);
01240   matrix_create (3, 3, &ULmatrix);
01241   if (ULmatrix.elts == NULL)
01242     {                           /* memory allocation error */
01243       fprintf (stderr, "could not allocate memory for Rmat \n");
01244       EXRETURN;
01245     }
01246 
01247   for (i = 0; i < 3; i++)
01248     {
01249       for (j = 0; j < 3; j++)
01250         {
01251           ULmatrix.elts[i][j] = Ematrix.elts[i][j] * eigs[i];
01252         }
01253     }
01254 
01255   matrix_multiply (ULmatrix, Ematrix, &Dmatrix);        /* new D is based on modified lambdas */
01256 #endif
01257 
01258 
01259 
01260   t1 = e10 * e10;
01261   t3 = e20 * e20;
01262   t5 = e30 * e30;
01263   t8 = e10 * l1;
01264   t10 = e20 * l2;
01265   t12 = e30 * l3;
01266   t14 = t8 * e11 + t10 * e21 + t12 * e31;
01267   t18 = t8 * e12 + t10 * e22 + t12 * e32;
01268   t19 = e11 * e11;
01269   t21 = e21 * e21;
01270   t23 = e31 * e31;
01271   t32 = e11 * l1 * e12 + e21 * l2 * e22 + e31 * l3 * e32;
01272   t33 = e12 * e12;
01273   t35 = e22 * e22;
01274   t37 = e32 * e32;
01275   Dmatrix.elts[0][0] = t1 * l1 + t3 * l2 + t5 * l3;
01276   Dmatrix.elts[0][1] = t14;
01277   Dmatrix.elts[0][2] = t18;
01278   Dmatrix.elts[1][0] = t14;
01279   Dmatrix.elts[1][1] = t19 * l1 + t21 * l2 + t23 * l3;
01280   Dmatrix.elts[1][2] = t32;
01281   Dmatrix.elts[2][0] = t18;
01282   Dmatrix.elts[2][1] = t32;
01283   Dmatrix.elts[2][2] = t33 * l1 + t35 * l2 + t37 * l3;
01284 
01285   udmatrix_to_vector (Dmatrix, &Dvector);       /* convert to vector format for D */
01286 
01287   EXRETURN;
01288 }
01289 
01290 /*! compute the diffusion weighting matrix bmatrix for q number of gradients */
01291 /* only need to calculate once */
01292 /* bq = diffusion weighting matrix of qth gradient */
01293 /*      GxGx GxGy GxGz
01294         GxGy GyGy GyGz
01295         GxGz GyGz GzGz */
01296 /* b0 is 0 for all 9 elements */
01297 /* bmatrix is really stored as 6 x npts array */
01298 static void
01299 Computebmatrix (MRI_IMAGE * grad1Dptr)
01300 {
01301   int i, n;
01302   register double *bptr;
01303   register float *Gxptr, *Gyptr, *Gzptr;
01304   double Gx, Gy, Gz;
01305 
01306   ENTRY ("Computebmatrix");
01307   n = grad1Dptr->nx;            /* number of gradients other than I0 */
01308   Gxptr = MRI_FLOAT_PTR (grad1Dptr);    /* use simple floating point pointers to get values */
01309   Gyptr = Gxptr + n;
01310   Gzptr = Gyptr + n;
01311 
01312   bptr = bmatrix;
01313   for (i = 0; i < 6; i++)
01314     *bptr++ = 0.0;              /* initialize first 6 elements to 0.0 for the I0 gradient */
01315 
01316   for (i = 0; i < n; i++)
01317     {
01318       Gx = *Gxptr++;
01319       Gy = *Gyptr++;
01320       Gz = *Gzptr++;
01321       *bptr++ = Gx * Gx;
01322       *bptr++ = Gx * Gy;
01323       *bptr++ = Gx * Gz;
01324       *bptr++ = Gy * Gy;
01325       *bptr++ = Gy * Gz;
01326       *bptr++ = Gz * Gz;
01327     }
01328   EXRETURN;
01329 }
01330 
01331 /*! compute non-gradient intensity, J, based on current calculated values of 
01332    diffusion tensor matrix, D */
01333 static double
01334 ComputeJ (float ts[], int npts)
01335 {
01336   /* J = Sum(wq Iq exp(-bq D)) / Sum (wq exp(-2bq D)) */
01337   /*     estimate of b0 intensity without noise and applied gradient */
01338   /* Iq = voxel value for qth gradient */
01339   /* bq = diffusion weighting matrix of qth gradient */
01340   /* wq = weighting factor for qth gradient at Iq voxel */
01341   /* D = estimated diffusion tensor matrix 
01342      [Dxx Dxy Dxz, Dxy Dyy Dyz, Dxz Dyz Dzz] */
01343   /* ts = Iq is time series voxel data from original data of intensities */
01344 
01345   register int i, j;
01346   double sum0, sum1, b1D1, b2D2, b4D4, wtexpbD, J, tempcalc, sumbD, Fscalar;
01347   double *expbD, *expbDptr, *wtfactorptr, *Ftempmatrix;
01348   register double *Fptr, *Rptr, *bptr;
01349   double D0,D1,D2,D3,D4,D5;
01350 
01351   ENTRY ("ComputeJ");
01352   sum0 = sum1 = 0.0;
01353   expbD = malloc (npts * sizeof (double));      /* allocate calculations for speed */
01354   expbDptr = expbD;             /* temporary pointers for indexing */
01355   wtfactorptr = wtfactor;
01356   bptr = bmatrix;               /* npts of b vectors (nx6) */
01357 
01358   D0 = Dmatrix.elts[0][0];
01359   D1 = Dmatrix.elts[0][1];
01360   D2 = Dmatrix.elts[0][2];
01361   D3 = Dmatrix.elts[1][1];
01362   D4 = Dmatrix.elts[1][2];
01363   D5 = Dmatrix.elts[2][2];
01364 
01365 
01366   for (i = 0; i < npts; i++)
01367     {
01368       /* compute bq.D */
01369       /* bq.D is large dot product of b and D at qth gradient */
01370       /* large dot product for Hilbert algebra */
01371       /* regular dot product is for Hilbert space (vectors only)- who knew? */
01372       /* calculate explicitly rather than loop to save time */
01373       b1D1 = *(bptr + 1) * D1;
01374       b1D1 += b1D1;
01375       b2D2 = *(bptr + 2) * D2;
01376       b2D2 += b2D2;
01377       b4D4 = *(bptr + 4) * D4;
01378       b4D4 += b4D4;
01379 
01380       sumbD = *bptr * D0 + b1D1 + b2D2 +        /* bxxDxx + 2bxyDxy +  2bxzDxz + */
01381         (*(bptr + 3) * D3) +    /* byyDyy + */
01382         b4D4 +                  /* 2byzDyz + */
01383         (*(bptr + 5) * D5);     /* bzzDzz */
01384 
01385       /*  exp (-bq.D) */
01386       *expbDptr = exp (-sumbD);
01387       wtexpbD = *(wtfactor + i) * *expbDptr;
01388       sum0 += wtexpbD * ts[i];
01389       sum1 += wtexpbD * *expbDptr;
01390       expbDptr++;
01391       wtfactorptr++;
01392       bptr += 6;                /* increment to next vector of bmatrix */
01393     }
01394 
01395   J = sum0 / sum1;
01396   /* Now compute error functional,E(D,J) and gradient of E with respect to D ,Ed or F in notes */
01397   /* E(D,J)= 1/2 Sum[wq (J exp(-bq.D) - Iq)^2] */
01398   /* F = Ed =  - Sum[wq (J exp(-bq.D) - Iq) bq] *//* Ed is a symmetric matrix */
01399   sum0 = 0.0;
01400   sigma = 0.0;                  /* standard deviation of noise for weight factors */
01401   expbDptr = expbD;
01402   wtfactorptr = wtfactor;
01403   /* initialize F matrix */
01404   Ftempmatrix = malloc(6*sizeof(double));
01405   Fptr = Ftempmatrix;
01406   for(i=0;i<6;i++)
01407     *Fptr++ = 0.0;
01408   Fptr = Ftempmatrix;
01409   Rptr = Rvector;               /* residuals calculated here - used in Wt.factor calculations */
01410   bptr = bmatrix;               /* npts of b vectors (nx6) */
01411   for (i = 0; i < npts; i++)
01412     {
01413       *Rptr = tempcalc = (J * *expbDptr) - ts[i];
01414       Fscalar = -*wtfactorptr * tempcalc;
01415       tempcalc = tempcalc * tempcalc;
01416 
01417       for (j = 0; j < 6; j++)
01418         {                       /* for each entry of Fij (Fxx, Fxy,...) */
01419           /* F = - Sum[wq (J exp(-bq.D) - Iq) bq] = Sum[-wq (J exp(-bq.D) - Iq) bq] */
01420           *(Fptr+j) += Fscalar * (*bptr++);     /*  Fij = Fij + (Fscalar * bij)  */
01421         }
01422 
01423       sum0 += *wtfactorptr * tempcalc;  /* E(D,J) = Sum (wq temp^2) */
01424       sigma += tempcalc;        /* standard deviation of noise for weight factors */
01425       expbDptr++;
01426       wtfactorptr++;
01427       Rptr++;
01428     }
01429 
01430   udmatrix_copy (Ftempmatrix, &Fmatrix);        /* copy upper diagonal vector data into full matrix */
01431 
01432   ED = sum0 / 2;                /* this is the error for this iteration */
01433 
01434   free (Ftempmatrix);
01435   free (expbD);
01436   RETURN (J);
01437 }
01438 
01439 /*! compute initial step size for gradient descent */
01440 static void
01441 ComputeDeltaTau ()
01442 {
01443   double sum0, sum1;
01444   matrix Dsqmatrix, FDsqmatrix, DsqFmatrix, Gmatrix;
01445   /* compute estimate of gradient, dD/dtau */
01446   /*G = [F] [D]^2 + [D]^2 [F] - ask Bob about ^2 and negative for this part to be sure */
01447 
01448   ENTRY ("ComputeDeltaTau");
01449   matrix_initialize (&Dsqmatrix);
01450   matrix_initialize (&FDsqmatrix);
01451   matrix_initialize (&DsqFmatrix);
01452   matrix_initialize (&Gmatrix);
01453 
01454   matrix_multiply (Dmatrix, Dmatrix, &Dsqmatrix);       /* compute D^2 */
01455   matrix_multiply (Fmatrix, Dsqmatrix, &FDsqmatrix);    /* FD^2 */
01456   matrix_multiply (Dsqmatrix, Fmatrix, &DsqFmatrix);    /* D^2F */
01457   matrix_add (FDsqmatrix, DsqFmatrix, &Gmatrix);        /* G= FD^2 +D^2F */
01458 
01459 
01460   /* deltatau = 0.01 * Sum(|Dij|) / Sum (|Gij|) */
01461   sum0 = matrix_sumabs (Dmatrix);
01462   sum1 = matrix_sumabs (Gmatrix);
01463   if (sum1 != 0.0)
01464     deltatau = 0.01 * sum0 / sum1;
01465   else
01466     deltatau = 0.0;
01467   matrix_destroy (&Dsqmatrix);
01468   matrix_destroy (&FDsqmatrix);
01469   matrix_destroy (&DsqFmatrix);
01470   matrix_destroy (&Gmatrix);
01471   EXRETURN;
01472 }
01473 
01474 /*! allocate all the global matrices and arrays once */
01475 static void
01476 InitGlobals (int npts)
01477 {
01478   int i;
01479   double *cumulativewtptr;
01480 
01481   ENTRY ("InitGlobals");
01482   matrix_initialize (&Fmatrix);
01483   matrix_create (3, 3, &Fmatrix);
01484   matrix_initialize (&Dmatrix);
01485   matrix_create (3, 3, &Dmatrix);
01486   matrix_initialize (&Hplusmatrix);
01487   matrix_create (3, 3, &Hplusmatrix);
01488   matrix_initialize (&Hminusmatrix);
01489   matrix_create (3, 3, &Hminusmatrix);
01490   matrix_initialize (&OldD);
01491   matrix_create (3, 3, &OldD);
01492   for(i=0;i<2;i++){
01493     matrix_initialize (&tempFmatrix[i]);
01494     matrix_create (3, 3, &tempFmatrix[i]);
01495     matrix_initialize (&tempDmatrix[i]);
01496     matrix_create (3, 3, &tempDmatrix[i]);
01497     matrix_initialize (&tempHplusmatrix[i]);
01498     matrix_create (3, 3, &tempHplusmatrix[i]);
01499     matrix_initialize (&tempHminusmatrix[i]);
01500     matrix_create (3, 3, &tempHminusmatrix[i]);
01501   }
01502   Rvector = malloc (npts * sizeof (double));
01503   tempRvector = malloc (npts * sizeof(double));
01504   wtfactor = malloc (npts * sizeof (double));
01505 
01506   if(cumulative_flag && reweight_flag) {
01507      cumulativewt = malloc (npts * sizeof (double));
01508      cumulativewtptr = cumulativewt;
01509      for(i=0;i<npts;i++)
01510         *cumulativewtptr++ = 0.0;
01511      rewtvoxels = 0;
01512   }
01513 
01514   bmatrix = malloc (npts * 6 * sizeof (double));
01515 
01516   vector_initialize (&Dvector); /* need to initialize vectors before 1st use-mod drg 12/20/2004 */
01517   /*  vector_initialize (&tempDvector);  vector_create(npts, &tempDvector);*/
01518   EXRETURN;
01519 }
01520 
01521 /*! free up all the matrices and arrays */
01522 static void
01523 FreeGlobals ()
01524 {
01525   int i;
01526 
01527   ENTRY ("FreeGlobals");
01528   matrix_destroy (&Fmatrix);
01529   matrix_destroy (&Dmatrix);
01530   matrix_destroy (&Hplusmatrix);
01531   matrix_destroy (&Hminusmatrix);
01532   matrix_destroy (&OldD);
01533   for(i=0;i<2;i++){
01534     matrix_destroy (&tempFmatrix[i]);
01535     matrix_destroy (&tempDmatrix[i]);
01536     matrix_destroy (&tempHplusmatrix[i]);
01537     matrix_destroy (&tempHminusmatrix[i]);
01538   }
01539 
01540 
01541   free (wtfactor);
01542   wtfactor = NULL;
01543   free (bmatrix);
01544   bmatrix = NULL;
01545   free (Rvector);
01546   Rvector = NULL;
01547   free (tempRvector);
01548   tempRvector = NULL;
01549   vector_destroy (&Dvector);    /* need to free elements of Dvector - mod-drg 12/20/2004 */
01550   /*  vector_destroy (&tempDvector);*/
01551   if(cumulative_flag && reweight_flag) {
01552     free (cumulativewt);
01553     cumulativewt = NULL;
01554   }
01555   EXRETURN;
01556 }
01557 
01558 /*! store current computed matrices D,Hp,Hm, R */
01559 static void
01560 Store_Computations (int i, int npts, int wtflag)
01561 {
01562   ENTRY ("Store_Computations");
01563 
01564   matrix_copy (Fmatrix, &tempFmatrix[i]);
01565   matrix_copy (Dmatrix, &tempDmatrix[i]);
01566   matrix_copy (Hplusmatrix, &tempHplusmatrix[i]);
01567   matrix_copy (Hminusmatrix, &tempHminusmatrix[i]);
01568   if(wtflag==1) 
01569     memcpy(tempRvector, Rvector, npts*sizeof(double));
01570   EXRETURN;
01571 }
01572 
01573 /*! restore old computed matrices D,Hp,Hm, R */
01574 static void
01575 Restore_Computations (int i, int npts, int wtflag)
01576 {
01577   ENTRY ("Restore_Computations");
01578 
01579   matrix_copy (tempFmatrix[i], &Fmatrix);
01580   matrix_copy (tempDmatrix[i], &Dmatrix);
01581   matrix_copy (tempHplusmatrix[i], &Hplusmatrix);
01582   matrix_copy (tempHminusmatrix[i], &Hminusmatrix);
01583   if(wtflag==1)
01584     memcpy(Rvector, tempRvector, npts*sizeof(double));
01585   EXRETURN;
01586 }
01587 
01588 /*! set all wt factors for all gradient levels to be 1.0 the first time through */
01589 static void
01590 InitWtfactors (int npts)
01591 {
01592   double *wtfactorptr;
01593   int i;
01594 
01595   ENTRY ("InitWtfactors");
01596   wtfactorptr = wtfactor;
01597   for (i = 0; i < npts; i++)
01598     *wtfactorptr++ = 1.0;
01599   EXRETURN;
01600 }
01601 
01602 /*! compute wt factors for each gradient level */
01603 static void
01604 ComputeWtfactors (int npts)
01605 {
01606   /* Residuals, rq, computed above in ComputeJ, stored in Rmatrix */
01607   /* unnormalized standard deviation, sigma, computed there too */
01608 /*  wq = 1 / sqrt(1 + (rq/sigma)^2)
01609     where sigma = sqrt[1/Nq Sum(rq^2)] */
01610 /*  and rq = J exp(-bq.D) - Iq */
01611 
01612   int i;
01613   double *wtfactorptr, *Rptr;
01614   double *cumulativewtptr;
01615   double tempcalc, sum;
01616 
01617   ENTRY ("ComputeWtfactors");
01618   sigma = sigma / npts;
01619   sigma = sqrt (sigma);         /* sigma = std.dev. */
01620 
01621   wtfactorptr = wtfactor;
01622   Rptr = Rvector;
01623 
01624   sum = 0.0;
01625   for (i = 0; i < npts; i++)
01626     {
01627       tempcalc = *Rptr++ / sigma;
01628       tempcalc = tempcalc * tempcalc;
01629       tempcalc = 1.0 / (sqrt (1 + tempcalc));
01630       *wtfactorptr++ = tempcalc;
01631       sum += tempcalc;
01632     }
01633   /* now renormalize to avoid changing the relative value of E(D) */
01634   tempcalc = npts / sum;     /* normalization factor */
01635   
01636   wtfactorptr = wtfactor;
01637   for (i=0; i<npts; i++) {
01638       *wtfactorptr = *wtfactorptr * tempcalc;
01639       wtfactorptr++;
01640   }
01641 
01642   if(cumulative_flag) {
01643      wtfactorptr = wtfactor;
01644      cumulativewtptr = cumulativewt;
01645      /*  printf("Wt.factors: ");*/
01646      ++rewtvoxels;
01647      for (i=0; i<npts; i++){
01648         *cumulativewtptr++ += *wtfactorptr++;   /* calculate cumulative wt.factor across all voxels*/
01649      }
01650   }
01651 
01652  EXRETURN;
01653 }
01654 
01655 /*! compute Hplus and Hminus as a function of delta tau */
01656 /* H+- = I +/- 1/2 deltatau F D */
01657 static void
01658 ComputeHpHm (double deltatau)
01659 {
01660   matrix FDmatrix;
01661   double dtau;
01662   int i, j;
01663 
01664   ENTRY ("ComputeHpHm");
01665   dtau = 0.5 * deltatau;
01666 
01667   matrix_initialize (&FDmatrix);
01668   matrix_multiply (Fmatrix, Dmatrix, &FDmatrix);
01669   for (i = 0; i < 3; i++)
01670     for (j = 0; j < 3; j++)
01671       FDmatrix.elts[i][j] = dtau * FDmatrix.elts[i][j];
01672 
01673   for (i = 0; i < 3; i++)
01674     {
01675       for (j = 0; j < 3; j++)
01676         {
01677           if (i == j)
01678             {
01679               Hplusmatrix.elts[i][j] = 1 + FDmatrix.elts[i][j]; /* I + dt/2 * FD */
01680               Hminusmatrix.elts[i][j] = 1 - FDmatrix.elts[i][j];        /* I - dt/2 * FD */
01681             }
01682           else
01683             {
01684               Hplusmatrix.elts[i][j] = Hminusmatrix.elts[i][j] =
01685                 FDmatrix.elts[i][j];
01686             }
01687         }
01688     }
01689 
01690   matrix_destroy (&FDmatrix);
01691   EXRETURN;
01692 }
01693 
01694 /*! compute new D matrix */
01695 /* D(tau+deltatau) = H-  H+^-1  D(tau)  H+^-1 H- */
01696 /*                 = A          D(tau)  A^T */
01697 /* where A = H- H+^-1 */
01698 static void
01699 ComputeNewD ()
01700 {
01701   double *Hpinv;
01702   matrix Hpinvmatrix, Amatrix, ATmatrix, ADmatrix;
01703 
01704   ENTRY ("ComputeNewD");
01705 
01706   Hpinv =
01707     InvertSym3 (Hplusmatrix.elts[0][0], Hplusmatrix.elts[0][1],
01708                 Hplusmatrix.elts[0][2], Hplusmatrix.elts[1][1],
01709                 Hplusmatrix.elts[1][2], Hplusmatrix.elts[2][2]);
01710   matrix_initialize (&Hpinvmatrix);
01711   matrix_initialize (&Amatrix);
01712   matrix_initialize (&ATmatrix);
01713   matrix_initialize (&ADmatrix);
01714 
01715   matrix_create (3, 3, &Hpinvmatrix);
01716   udmatrix_copy (Hpinv, &Hpinvmatrix);  /* copy values from Hpinv vector to Hpinvmatrix */
01717 
01718   matrix_multiply (Hminusmatrix, Hpinvmatrix, &Amatrix);
01719   matrix_multiply (Amatrix, Dmatrix, &ADmatrix);
01720   matrix_transpose (Amatrix, &ATmatrix);
01721   matrix_multiply (ADmatrix, ATmatrix, &Dmatrix);
01722 
01723   matrix_destroy (&ADmatrix);
01724   matrix_destroy (&ATmatrix);
01725   matrix_destroy (&Amatrix);
01726   matrix_destroy (&Hpinvmatrix);
01727 
01728   free (Hpinv);
01729   EXRETURN;
01730 }
01731 
01732 /*! test convergence of calculation of D */
01733 /* if sum of differences hasn't changed by more than 1E-4 */
01734 /*  then the calculations have converged */
01735 /* return 1 for convergence, 0 if not converged */
01736 static int
01737 TestConvergence(matrix NewD, matrix OldD)
01738 { 
01739   int converge;
01740   double convergence;
01741 
01742   ENTRY ("TestConvergence");
01743   /* convergence test */
01744   convergence = fabs (NewD.elts[0][0] - OldD.elts[0][0]) +      /* Dxx */
01745     fabs (NewD.elts[0][1] - OldD.elts[0][1]) +  /* Dxy */
01746     fabs (NewD.elts[0][2] - OldD.elts[0][2]) +  /* Dxz */
01747     fabs (NewD.elts[1][1] - OldD.elts[1][1]) +  /* Dyy */
01748     fabs (NewD.elts[1][2] - OldD.elts[1][2]) +  /* Dyz */
01749     fabs (NewD.elts[2][2] - OldD.elts[2][2]);   /* Dzz */
01750 
01751   if (convergence < SMALLNUMBER)
01752     converge = 1;
01753   else
01754     converge = 0;
01755 
01756   RETURN (converge);
01757 }
01758 
01759 /*! copy an upper diagonal matrix (6 point vector really) into a standard double
01760    array type matrix for n timepoints */
01761 /* ud0 ud1 ud2         m0 m1 m2
01762        ud3 ud4   -->   m3 m4 m5
01763            ud5         m6 m7 m8 */
01764 static void
01765 udmatrix_copy (double *udptr, matrix * m)
01766 {
01767   ENTRY ("udmatrix_copy");
01768 
01769   m->elts[0][0] = *udptr;
01770   m->elts[0][1] = *(udptr + 1);
01771   m->elts[0][2] = *(udptr + 2);
01772   m->elts[1][0] = *(udptr + 1);
01773   m->elts[1][1] = *(udptr + 3);
01774   m->elts[1][2] = *(udptr + 4);
01775   m->elts[2][0] = *(udptr + 2);
01776   m->elts[2][1] = *(udptr + 4);
01777   m->elts[2][2] = *(udptr + 5);
01778   EXRETURN;
01779 }
01780 
01781 /*! copy upper part of 3x3 matrix elements to 6-element vector elements */
01782 /* m1 m2 m3
01783       m4 m5   ->  v = [m1 m2 m3 m4 m5 m6]
01784          m6
01785 */
01786 static void
01787 udmatrix_to_vector (matrix m, vector * v)
01788 {
01789   ENTRY ("udmatrix_to_vector");
01790   v->elts[0] = m.elts[0][0];
01791   v->elts[1] = m.elts[0][1];
01792   v->elts[2] = m.elts[0][2];
01793   v->elts[3] = m.elts[1][1];
01794   v->elts[4] = m.elts[1][2];
01795   v->elts[5] = m.elts[2][2];
01796   EXRETURN;
01797 }
01798 
01799 /*! sum the absolute value of all  elements of a matrix */
01800 static double
01801 matrix_sumabs (matrix m)
01802 {
01803   register int i, j;
01804   register double sum;
01805 
01806   ENTRY ("matrix_sumabs");
01807   sum = 0.0;
01808   for (i = 0; i < 3; i++)
01809     {
01810       for (j = 0; j < 3; j++)
01811         sum += fabs (m.elts[i][j]);
01812     }
01813   RETURN (sum);
01814 }
01815 
01816 /*! calculate inverse of a symmetric 3x3 matrix */
01817 /* returns pointer to 9 element vector corresponding to 3x3 matrix */
01818 /*  a b c */
01819 /*  b e f */
01820 /*  c f i */
01821 /* Maple generated code */
01822 static double *
01823 InvertSym3 (double a, double b, double c, double e, double f, double i)
01824 {
01825   double *symmat, *symmatptr;   /* invert matrix - actually 6 values in a vector form */
01826   double t2, t4, t7, t9, t12, t15, t20, t24, t30;
01827 
01828   ENTRY ("InvertSym3");
01829   symmat = malloc (6 * sizeof (double));
01830   symmatptr = symmat;
01831   t2 = f * f;
01832   t4 = a * e;
01833   t7 = b * b;
01834   t9 = c * b;
01835   t12 = c * c;
01836   t15 = 1 / (t4 * i - a * t2 - t7 * i + 2.0 * t9 * f - t12 * e);
01837   t20 = (b * i - c * f) * t15;
01838   t24 = (b * f - c * e) * t15;
01839   t30 = (a * f - t9) * t15;
01840 
01841   *symmatptr++ = (e * i - t2) * t15;    /*B[0][0] */
01842   *symmatptr++ = -t20;          /*B[0][1] */
01843   *symmatptr++ = t24;           /*B[0][2] */
01844   /* B[1][0] = -t20; */
01845   *symmatptr++ = (a * i - t12) * t15;   /* B[1][1] */
01846   *symmatptr++ = -t30;          /* B [1][2] */
01847   /* B[2][0] = t24; */
01848   /* B[2][1] = -t30; */
01849   *symmatptr = (t4 - t7) * t15; /* B[2][2] */
01850 
01851   RETURN (symmat);
01852 }
01853 
01854 /*! copy elements from matrix a to matrix b */
01855 /*  matrix_equate already exists but creates and initializes new matrix */
01856 /*  steps we don't need to do here */
01857 /* assumes both a and b already exist with equal dimensions */
01858 static void
01859 matrix_copy (matrix a, matrix * b)
01860 {
01861   register int i;
01862   register int rows, cols;
01863 
01864   ENTRY ("matrix_copy");
01865 
01866   rows = a.rows;
01867   cols = a.cols;
01868 
01869   for (i = 0; i < rows; i++)
01870     {
01871 #if 0
01872       register int j;
01873       for (j = 0; j < cols; j++)
01874         b->elts[i][j] = a.elts[i][j];
01875 #else
01876       if (cols > 0)
01877         memcpy (b->elts[i], a.elts[i], sizeof (double) * cols); /* RWCox */
01878 #endif
01879     }
01880   EXRETURN;
01881 }
01882 
01883 
01884 #define DWI_WriteCheckWaitMax 2000
01885 #define DWI_WriteCheckWait 400
01886 /*-----------------------------------------------------*/
01887 /* Stuff for an extra NIML port for non-SUMA programs. */
01888 
01889 #ifndef NIML_TCP_FIRST_PORT
01890 #define NIML_TCP_FIRST_PORT 53212
01891 #endif
01892 
01893 /*! open NIML stream */
01894 static int DWI_Open_NIML_stream()
01895 {
01896    int nn, Wait_tot, tempport;
01897    char streamname[256];
01898 
01899    ENTRY("DWI_Open_NIML_stream");
01900 
01901    /* contact afni */
01902    tempport = NIML_TCP_FIRST_PORT;
01903    sprintf(streamname, "tcp:localhost:%d",tempport);
01904    fprintf(stderr,"Contacting AFNI\n");
01905  
01906    DWIstreamid =  NI_stream_open( streamname, "w" ) ;
01907    if (DWIstreamid==0) {
01908        fprintf(stderr,"+++ Warning - NI_stream_open failed\n");
01909       DWIstreamid = NULL;
01910       RETURN(1);
01911    }
01912 
01913    fprintf (stderr,"\nTrying shared memory...\n");
01914    if (!NI_stream_reopen( DWIstreamid, "shm:DWIDT1M:1M" ))
01915        fprintf (stderr, "Warning: Shared memory communcation failed.\n");
01916    else
01917      fprintf(stderr, "Shared memory connection OK.\n");
01918    Wait_tot = 0;
01919    while(Wait_tot < DWI_WriteCheckWaitMax){
01920       nn = NI_stream_writecheck( DWIstreamid , DWI_WriteCheckWait) ;
01921       if( nn == 1 ){ 
01922          fprintf(stderr,"\n") ; 
01923          RETURN(0) ; 
01924       }
01925       if( nn <  0 ){ 
01926          fprintf(stderr,"Bad connection to AFNI\n"); 
01927          DWIstreamid = NULL;
01928          RETURN(1);
01929       }
01930       Wait_tot += DWI_WriteCheckWait;
01931       fprintf(stderr,".") ;
01932    }
01933 
01934    fprintf(stderr,"+++ WriteCheck timed out (> %d ms).\n",DWI_WriteCheckWaitMax);
01935    RETURN(1);
01936 }
01937 
01938 /*! create the initial graph in AFNI - no points yet*/
01939 static int DWI_NIML_create_graph()
01940 {
01941    NI_element *nel;
01942 
01943    ENTRY("DWI_NIML_create_graph");
01944    nel = NI_new_data_element("ni_do", 0);
01945    NI_set_attribute ( nel, "ni_verb", "DRIVE_AFNI");
01946    NI_set_attribute ( nel, "ni_object", "OPEN_GRAPH_1D DWIConvEd 'DWI Convergence' 25 1 'converge step' 1 0 300000 Ed");
01947    if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
01948       fprintf(stderr,"Failed to send data to AFNI\n");
01949       NI_free_element(nel) ; nel = NULL;
01950       RETURN(1);
01951    }
01952    NI_free_element(nel) ; 
01953    nel = NULL;
01954    RETURN(0);
01955 }
01956 
01957 /*! create new graph with left and right y axes scaled from 0 to max1, max2*/
01958 static int DWI_NIML_create_newgraph(npts, max1, max2)
01959 int npts;
01960 double max1, max2;
01961 {
01962    NI_element *nel;
01963    char stmp[256];
01964    static int nx = -1;
01965 
01966    ENTRY("DWI_NIML_create_newgraph");
01967    nel = NI_new_data_element("ni_do", 0);
01968    NI_set_attribute ( nel, "ni_verb", "DRIVE_AFNI");
01969    if((nx==-1) || (nx<npts))                 /* 1st time through close any existing graph by that name*/
01970       NI_set_attribute ( nel, "ni_object","CLOSE_GRAPH_1D DWIConvEd\n"); /* have to close graph to change axes */
01971    else
01972       NI_set_attribute ( nel, "ni_object","CLEAR_GRAPH_1D DWIConvEd\n");
01973 
01974    if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
01975       fprintf(stderr,"Failed to send data to AFNI\n");
01976       NI_free_element(nel) ; nel = NULL;
01977       RETURN(1);
01978    }
01979 
01980    if((nx==-1) || (nx<npts)) {             /* update the graph only first time or if x-axis not big enough */
01981       nx = max_iter * 4  + 10;
01982       if(reweight_flag==1)
01983         nx += max_iter_rw * 4 + 10;
01984       if(nx<npts)                          /* fix graph to include largest number of steps */
01985         nx = npts;
01986       sprintf(stmp,"OPEN_GRAPH_1D DWIConvEd 'DWI Convergence' %d 1 'converge step' 2 0 100 %%Maximum Ed \\Delta\\tau\n",nx  );
01987       NI_set_attribute ( nel, "ni_object", stmp);
01988       if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
01989         fprintf(stderr,"Failed to send data to AFNI\n");
01990         NI_free_element(nel) ; nel = NULL;
01991         RETURN(1);
01992       }
01993       NI_set_attribute ( nel, "ni_object", "SET_GRAPH_GEOM DWIConvEd geom=700x400+100+400\n");
01994       if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
01995         fprintf(stderr,"Failed to send data to AFNI\n");
01996         NI_free_element(nel) ; nel = NULL;
01997         RETURN(1);
01998       }
01999    }
02000    NI_free_element(nel) ; 
02001    nel = NULL;
02002    RETURN(0);
02003 }
02004 
02005 /*! tell AFNI to graph data for convergence steps */
02006 static void DWI_AFNI_update_graph(double *Edgraph, double *dtau, int npts)
02007 {
02008    NI_element *nel;
02009    char stmp[256];
02010    int i;
02011    double Edmax, dtaumax;
02012    double *Edptr, *dtauptr;
02013    double tempEd, temptau;
02014 
02015    ENTRY("DWI_AFNI_update_graph");
02016 
02017    Edmax = 0.0; dtaumax = 0.0;
02018    Edptr = Edgraph;
02019    dtauptr = dtau;
02020    for(i=0;i<npts;i++) {
02021      if(*Edptr>Edmax)
02022        Edmax = *Edptr;
02023      if(*dtauptr>dtaumax)
02024        dtaumax = *dtauptr;
02025      ++Edptr; ++dtauptr;
02026    }
02027 
02028    NI_write_procins(DWIstreamid, "keep reading");
02029    DWI_NIML_create_newgraph(npts, Edmax, dtaumax);
02030    /* NI_sleep(250);*/
02031 
02032    nel = NI_new_data_element("ni_do", 0);
02033    NI_set_attribute ( nel, "ni_verb", "DRIVE_AFNI");
02034    NI_set_attribute ( nel, "ni_object", "CLEAR_GRAPH_1D DWIConvEd\n");
02035    /*      NI_sleep(25);*/
02036    if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
02037       fprintf(stderr,"Failed to send data to AFNI\n");
02038    }
02039 
02040 
02041    for(i=0;i<npts;i++){
02042       if(Edmax!=0.0)
02043          tempEd = 100*Edgraph[i] / Edmax;
02044       else
02045         tempEd = 0.0;
02046       if(dtaumax!=0.0)
02047         temptau = 100*dtau[i] / dtaumax;
02048       else
02049         temptau = 0.0;
02050 
02051       sprintf(stmp,"ADDTO_GRAPH_1D DWIConvEd %4.2f %4.2f\n", tempEd, temptau);  /* put rel.error, Ed, and deltatau for all the convergence steps */
02052       NI_set_attribute ( nel, "ni_object", stmp);  /* put command and data in stmp */
02053       NI_sleep(25);    /* for dramatic effect */
02054       if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
02055         fprintf(stderr,"Failed to send data to AFNI\n");
02056       }
02057    }
02058    NI_free_element(nel) ; nel = NULL;
02059    EXRETURN;
02060 }
02061