00001
00002
00003
00004
00005
00006
00007
00008 #include "thd_shear3d.h"
00009
00010 # include "matrix.h"
00011 # include "matrix.c"
00012
00013 #include "afni.h"
00014
00015 #define TINYNUMBER 1E-10
00016 #define SMALLNUMBER 1E-4
00017 #define MAX_CONVERGE_STEPS 10
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;
00024 static double *tempRvector;
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
00034
00035 static byte *maskptr = NULL;
00036 static double eigs[12];
00037 static double deltatau;
00038 static double *wtfactor;
00039 static double *bmatrix;
00040 static double *cumulativewt;
00041 static long rewtvoxels;
00042 static double sigma;
00043 static double ED;
00044 static int automask = 0;
00045 static int reweight_flag = 0;
00046 static int method = -1;
00047 static int max_iter = -2;
00048 static int max_iter_rw = -2;
00049 static int eigs_flag = 0;
00050 static int cumulative_flag = 0;
00051 static int debug_briks = 0;
00052 static int verbose = 0;
00053 static int afnitalk_flag = 0;
00054
00055 static NI_stream_type * DWIstreamid = 0;
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;
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
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;
00163 method = -1;
00164 reweight_flag = 0;
00165
00166 datum = MRI_float;
00167 while (nopt < argc && argv[nopt][0] == '-')
00168 {
00169
00170
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);
00180
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
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;
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
00418
00419 if (nopt >= argc)
00420 {
00421 fprintf (stderr, "*** Error - No input dataset!?\n");
00422 exit (1);
00423 }
00424
00425
00426
00427
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);
00453
00454 nopt++;
00455
00456
00457
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
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);
00484 Computebmatrix (grad1Dptr);
00485
00486 if (automask)
00487 {
00488 DSET_mallocize (old_dset);
00489 DSET_load (old_dset);
00490
00491 anat_im = DSET_BRICK (old_dset, 0);
00492 maskptr = mri_automask_image (anat_im);
00493 #if 0
00494
00495 nvox = DSET_NVOX (old_dset);
00496 sar = (short *) calloc (nvox, sizeof (short));
00497
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
00509 EDIT_add_brick (old_dset, MRI_short, 0.0, sar);
00510
00511 #endif
00512 }
00513
00514
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) {
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
00530 NI_stream_close(DWIstreamid);
00531 DWIstreamid = 0;
00532 }
00533 }
00534
00535
00536
00537 new_dset = MAKER_4D_to_typed_fbuc (old_dset,
00538 prefix,
00539 datum,
00540 0,
00541 0,
00542 nbriks,
00543 DWItoDT_tsfunc,
00544 NULL
00545 );
00546
00547
00548 if(afnitalk_flag && (DWIstreamid!=0)) {
00549
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);
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;
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
00626 static void
00627 Form_R_Matrix (MRI_IMAGE * grad1Dptr)
00628 {
00629 matrix Rmat;
00630 register double sf = 1.0;
00631 register double sf2;
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);
00641 if (Rmat.elts == NULL)
00642 {
00643 fprintf (stderr, "could not allocate memory for Rmat \n");
00644 EXRETURN;
00645 }
00646 sf2 = sf + sf;
00647 Gxptr = imptr = MRI_FLOAT_PTR (grad1Dptr);
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;
00657 Rmat.elts[i][1] = sf2 * Gx * Gy;
00658 Rmat.elts[i][2] = sf2 * Gx * Gz;
00659 Rmat.elts[i][3] = sf * Gy * Gy;
00660 Rmat.elts[i][4] = sf2 * Gy * Gz;
00661 Rmat.elts[i][5] = sf * Gz * Gz;
00662 }
00663
00664 matrix_initialize (&Rtmat);
00665 matrix_psinv (Rmat, nullptr, &Rtmat);
00666 matrix_destroy (&Rmat);
00667 EXRETURN;
00668 }
00669
00670
00671
00672
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;
00688 int max_converge_step, graphpoint;
00689 double dtau[50], Edgraph[50];
00690 int graphflag;
00691
00692 ENTRY ("DWItoDT_tsfunc");
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702 if (val == NULL)
00703 {
00704
00705 if (npts > 0)
00706 {
00707 nvox = npts;
00708 ncall = 0;
00709 noisecall = 0;
00710 }
00711 else
00712 {
00713
00714
00715 }
00716 EXRETURN;
00717 }
00718
00719 ncall++;
00720
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 {
00729 for (i = 0; i < nbriks; i++)
00730 val[i] = 0.0;
00731 if(debug_briks)
00732 val[nbriks-4] = -3.0;
00733 EXRETURN;
00734 }
00735 }
00736
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);
00749 else
00750 lnvector.elts[i] = 0.0;
00751 }
00752
00753 vector_multiply (Rtmat, lnvector, &Dvector);
00754
00755 vector_destroy (&lnvector);
00756
00757 if((method==0)||(max_iter==-1)) {
00758 vector_to_array(Dvector, val);
00759
00760 if(debug_briks) {
00761 InitWtfactors (npts);
00762 J = ComputeJ (ts, npts);
00763 val[nbriks-4] = -1.0;
00764
00765 val[nbriks-3] = ED;
00766 val[nbriks-2] = ED;
00767 val[nbriks-1] = J;
00768 }
00769
00770 if(eigs_flag) {
00771 EIG_func();
00772 for(i=0;i<12;i++)
00773 val[i+6] = eigs[i];
00774
00775 val[18] = Calc_FA(val+6);
00776 val[19] = Calc_MD(val+6);
00777 }
00778 EXRETURN;
00779 }
00780
00781
00782
00783
00784 EIG_func ();
00785
00786
00787 InitWtfactors (npts);
00788
00789 ComputeD0 ();
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;
00797 val[nbriks-3] = 0;
00798 val[nbriks-2] = 0;
00799 val[nbriks-1] = 0;
00800 }
00801
00802 EXRETURN;
00803 }
00804
00805
00806 converge_step = 0;
00807 max_converge_step = max_iter;
00808 converge = 0;
00809 wtflag = reweight_flag;
00810
00811
00812 J = ComputeJ (ts, npts);
00813 Store_Computations (0, npts, wtflag);
00814 matrix_copy (Dmatrix, &OldD);
00815
00816 if(debug_briks)
00817 val[nbriks-2] = ED;
00818
00819 EDold = ED;
00820 ComputeDeltaTau ();
00821 if(deltatau<=TINYNUMBER) {
00822 for(i=0;i<nbriks;i++)
00823 val[i] = 0.0;
00824 if(debug_briks) {
00825 val[nbriks-4] = -1.0;
00826 val[nbriks-1] = J;
00827 }
00828 EXRETURN;
00829 }
00830
00831 if(verbose&&(!(noisecall%verbose)))
00832 recordflag = 1;
00833 else
00834 recordflag = 0;
00835
00836 if(afnitalk_flag&&(!(noisecall%afnitalk_flag))) {
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
00850
00851
00852
00853
00854
00855 trialstep = 1;
00856 ntrial = 0;
00857 while ((trialstep) && (ntrial < 10))
00858 {
00859 ComputeHpHm (deltatau);
00860 ComputeNewD ();
00861 J = ComputeJ (ts, npts);
00862 if (ED < EDold)
00863 {
00864
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);
00876 }
00877 else {
00878 Restore_Computations (0, npts, wtflag);
00879 dt = deltatau / 2;
00880 deltatau = dt;
00881
00882 ntrial++;
00883 }
00884 }
00885
00886
00887
00888
00889
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))) {
00901 Restore_Computations (0, npts, wtflag);
00902 ComputeHpHm (deltatau);
00903 ComputeNewD ();
00904 J = ComputeJ (ts, npts);
00905
00906 if(ED<EDold){
00907 adjuststep = 0;
00908 Store_Computations(0, npts, wtflag);
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){
00924 ED = EDold;
00925 deltatau = orig_deltatau;
00926 Restore_Computations (1, npts, wtflag);
00927
00928 Store_Computations(0, npts, wtflag);
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) {
00940
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);
00961 ED = EDold;
00962 }
00963
00964
00965 if(((converge) || (converge_step==max_iter)) && wtflag && reweight_flag) {
00966 converge = 0;
00967 ComputeWtfactors (npts);
00968 converge_step = 1;
00969 max_converge_step = max_iter_rw+1;
00970 wtflag = 0;
00971 J=ComputeJ(ts, npts);
00972 EDold = ED;
00973 }
00974 }
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) {
00986 udmatrix_to_vector(Dmatrix, &Dvector);
00987 EIG_func();
00988 for(i=0;i<12;i++)
00989 val[i+6] = eigs[i];
00990
00991 val[18] = Calc_FA(val+6);
00992 val[19] = Calc_MD(val+6);
00993 }
00994
00995
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); ;
01002 }
01003 if(graphflag==1) {
01004 DWI_AFNI_update_graph(Edgraph, dtau, graphpoint);
01005 }
01006
01007 EXRETURN;
01008 }
01009
01010
01011
01012
01013 static void
01014 EIG_func ()
01015 {
01016
01017
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
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
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);
01050
01051 maxindex = 2;
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 {
01062 maxindex = i;
01063 maxvalue = temp;
01064 }
01065 if (temp < minvalue)
01066 {
01067 minindex = i;
01068 minvalue = temp;
01069 }
01070 }
01071
01072 for (i = 0; i < 3; i++)
01073 {
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
01086 for (i = 0; i < 3; i++)
01087 {
01088 eigs[i] = e[sortvector[i]];
01089
01090
01091 astart = sortvector[i] * 3;
01092 vstart = (i + 1) * 3;
01093
01094 for (j = 0; j < 3; j++)
01095 {
01096 eigs[vstart + j] = a[astart + j];
01097 }
01098 }
01099
01100 EXRETURN;
01101 }
01102
01103
01104
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
01114
01115
01116 if((val[0]<=0.0)||(val[1]<=0.0)||(val[2]<=0.0)) {
01117 RETURN(0.0);
01118 }
01119
01120 ssq = (val[0]*val[0])+(val[1]*val[1])+(val[2]*val[2]);
01121
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;
01130
01131 if(ssq!=0.0)
01132 FA = sqrt(dsq/(2.0*ssq));
01133 else
01134 FA = 0.0;
01135
01136 RETURN(FA);
01137 }
01138
01139
01140
01141 static float
01142 Calc_MD(float *val)
01143 {
01144 float MD;
01145
01146 ENTRY("Calc_MD");
01147
01148
01149
01150
01151 if((val[0]<=0.0)||(val[1]<=0.0)||(val[2]<=0.0)) {
01152 RETURN(0.0);
01153 }
01154 MD = (val[0] + val[1] + val[2]) / 3;
01155
01156 RETURN(MD);
01157 }
01158
01159
01160
01161
01162
01163
01164 static void
01165 ComputeD0 ()
01166 {
01167 int i, j;
01168
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
01177
01178 if (eigs[0] < 0)
01179 {
01180
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);
01197 EXRETURN;
01198 }
01199
01200 mu = 0.2 * eigs[0];
01201 if (eigs[1] < mu)
01202 eigs[1] = mu;
01203 if (eigs[2] < mu)
01204 eigs[2] = mu;
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
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
01238
01239 matrix_initialize (&ULmatrix);
01240 matrix_create (3, 3, &ULmatrix);
01241 if (ULmatrix.elts == NULL)
01242 {
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);
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);
01286
01287 EXRETURN;
01288 }
01289
01290
01291
01292
01293
01294
01295
01296
01297
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;
01308 Gxptr = MRI_FLOAT_PTR (grad1Dptr);
01309 Gyptr = Gxptr + n;
01310 Gzptr = Gyptr + n;
01311
01312 bptr = bmatrix;
01313 for (i = 0; i < 6; i++)
01314 *bptr++ = 0.0;
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
01332
01333 static double
01334 ComputeJ (float ts[], int npts)
01335 {
01336
01337
01338
01339
01340
01341
01342
01343
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));
01354 expbDptr = expbD;
01355 wtfactorptr = wtfactor;
01356 bptr = bmatrix;
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
01369
01370
01371
01372
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 +
01381 (*(bptr + 3) * D3) +
01382 b4D4 +
01383 (*(bptr + 5) * D5);
01384
01385
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;
01393 }
01394
01395 J = sum0 / sum1;
01396
01397
01398
01399 sum0 = 0.0;
01400 sigma = 0.0;
01401 expbDptr = expbD;
01402 wtfactorptr = wtfactor;
01403
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;
01410 bptr = bmatrix;
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 {
01419
01420 *(Fptr+j) += Fscalar * (*bptr++);
01421 }
01422
01423 sum0 += *wtfactorptr * tempcalc;
01424 sigma += tempcalc;
01425 expbDptr++;
01426 wtfactorptr++;
01427 Rptr++;
01428 }
01429
01430 udmatrix_copy (Ftempmatrix, &Fmatrix);
01431
01432 ED = sum0 / 2;
01433
01434 free (Ftempmatrix);
01435 free (expbD);
01436 RETURN (J);
01437 }
01438
01439
01440 static void
01441 ComputeDeltaTau ()
01442 {
01443 double sum0, sum1;
01444 matrix Dsqmatrix, FDsqmatrix, DsqFmatrix, Gmatrix;
01445
01446
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);
01455 matrix_multiply (Fmatrix, Dsqmatrix, &FDsqmatrix);
01456 matrix_multiply (Dsqmatrix, Fmatrix, &DsqFmatrix);
01457 matrix_add (FDsqmatrix, DsqFmatrix, &Gmatrix);
01458
01459
01460
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
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);
01517
01518 EXRETURN;
01519 }
01520
01521
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);
01550
01551 if(cumulative_flag && reweight_flag) {
01552 free (cumulativewt);
01553 cumulativewt = NULL;
01554 }
01555 EXRETURN;
01556 }
01557
01558
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
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
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
01603 static void
01604 ComputeWtfactors (int npts)
01605 {
01606
01607
01608
01609
01610
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);
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
01634 tempcalc = npts / sum;
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
01646 ++rewtvoxels;
01647 for (i=0; i<npts; i++){
01648 *cumulativewtptr++ += *wtfactorptr++;
01649 }
01650 }
01651
01652 EXRETURN;
01653 }
01654
01655
01656
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];
01680 Hminusmatrix.elts[i][j] = 1 - FDmatrix.elts[i][j];
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
01695
01696
01697
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);
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
01733
01734
01735
01736 static int
01737 TestConvergence(matrix NewD, matrix OldD)
01738 {
01739 int converge;
01740 double convergence;
01741
01742 ENTRY ("TestConvergence");
01743
01744 convergence = fabs (NewD.elts[0][0] - OldD.elts[0][0]) +
01745 fabs (NewD.elts[0][1] - OldD.elts[0][1]) +
01746 fabs (NewD.elts[0][2] - OldD.elts[0][2]) +
01747 fabs (NewD.elts[1][1] - OldD.elts[1][1]) +
01748 fabs (NewD.elts[1][2] - OldD.elts[1][2]) +
01749 fabs (NewD.elts[2][2] - OldD.elts[2][2]);
01750
01751 if (convergence < SMALLNUMBER)
01752 converge = 1;
01753 else
01754 converge = 0;
01755
01756 RETURN (converge);
01757 }
01758
01759
01760
01761
01762
01763
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
01782
01783
01784
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
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
01817
01818
01819
01820
01821
01822 static double *
01823 InvertSym3 (double a, double b, double c, double e, double f, double i)
01824 {
01825 double *symmat, *symmatptr;
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;
01842 *symmatptr++ = -t20;
01843 *symmatptr++ = t24;
01844
01845 *symmatptr++ = (a * i - t12) * t15;
01846 *symmatptr++ = -t30;
01847
01848
01849 *symmatptr = (t4 - t7) * t15;
01850
01851 RETURN (symmat);
01852 }
01853
01854
01855
01856
01857
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);
01878 #endif
01879 }
01880 EXRETURN;
01881 }
01882
01883
01884 #define DWI_WriteCheckWaitMax 2000
01885 #define DWI_WriteCheckWait 400
01886
01887
01888
01889 #ifndef NIML_TCP_FIRST_PORT
01890 #define NIML_TCP_FIRST_PORT 53212
01891 #endif
01892
01893
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
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
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
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))
01970 NI_set_attribute ( nel, "ni_object","CLOSE_GRAPH_1D DWIConvEd\n");
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)) {
01981 nx = max_iter * 4 + 10;
01982 if(reweight_flag==1)
01983 nx += max_iter_rw * 4 + 10;
01984 if(nx<npts)
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
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
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
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);
02052 NI_set_attribute ( nel, "ni_object", stmp);
02053 NI_sleep(25);
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