1 /*********************** 3dDWItoDT.c **********************************************/
2 /* Author: Daniel Glen, 17 Nov 2004 */
3 /* compute 6 principle direction tensors from multiple gradient vectors*/
4 /* and corresponding image data */
5 /* This version includes a more complex algorithm that takes into account*/
6 /* noise.*/
7
8
9 /* Gregory Baxter (gbaxter@ucsd.edu), 2 Nov 2012 */
10 /* UCSD Center for Scientific Computation in Imaging (csci.ucsd.edu) */
11 /* */
12 /* Added modifications to allow for reading a b-matrix file directly rather */
13 /* constructing it from gradient directions. */
14 /* */
15 /* When -bmatrix is specified on the command line, the gradient file is */
16 /* read as a bmatrix file instead and expected to have 6 columns (rather */
17 /* than 3), with elements ordered Bxx, Byy, Bzz, Bxy, Bxz, Byz. */
18 /* */
19 /* For calculations with the -bmatrix option, the diffusion weighted dataset */
20 /* is assumed to have a single b0 image first, followed by all the diffusion */
21 /* weighted images. The bmatrix file should have the same number of rows as */
22 /* the total number of diffusion weighted images (including the b0 image). */
23
24
25 // Dec, 2013 (PT):
26 // + (probable) bug fix and new matrix format allowed for input (no
27 // zeros at top). search for '$$' to see the new lines for allowing
28 // bmatrix input with no zeros in top row, and for '@@' for potential
29 // bug fix for original `-bmatrix' options.
30
31 // Apr, 2016 (PT):
32 // + allow for non-unit gradient magnitudes, so that gradients can be
33 // weighted by b-values- affects Form_R_Matrix() and computebmatrix().
34 // search for 'apt,2016' to see changes.
35
36 // May, 2016 (PT):
37 // + extra badness condition for when fits go really wrong, for example if
38 // data is missing. By default, will censor when values are >100*(MD_0),
39 // where MD_0 is the default 'CSF value' based on the b_max.
40
41 // (end of) May, 2016 (PT):
42 // + allow user to scaled DT/MD/L? values by 1000 internally. If
43 // bweighted matrix/gradient values have been input, this may be
44 // nicer than having 0.0007 (mm^2/s), etc., values, changing them to
45 // ~0.7 (x10^{-3} mm^2/s), etc., which is more the style propounded in
46 // reported literature.
47
48 // Sep, 2016 (PT):
49 // + introduce user opt BMAX_REF into program, so we can identify if
50 // bvalue is really > 0, which can happen
51
52 // Oct, 2016 (PT):
53 // + now calcs/outputs RD as part of '-eigs'
54
55 // May, 2017 (PT):
56 // + can now calc goodness of fit chi-squared values (signal-based
57 // ones, at the moment), if user wants; now happens when
58 // '-debug_briks' is used.
59
60 // Aug, 2017 (PT):
61 // + can now have cumulative weights dumped into a file, as well as
62 // shown in terminal -> CORRECTLY now...
63
64
65 #include "thd_shear3d.h"
66 /*#ifndef FLOATIZE*/
67 # include "matrix.h"
68 /*#endif*/
69 #include "mrilib.h"
70
71 #define SIG_EPS 1E-8
72 #define TINYNUMBER 1E-10
73 #define SMALLNUMBER 1E-4
74 #define HUGENUMBER 1E38
75 #define MAX_CONVERGE_STEPS 10 /* default maximum steps */
76 #define MAX_RWCONVERGE_STEPS 5
77
78 static char prefix[THD_MAX_PREFIX] = "DT";
79 static int datum = MRI_float;
80 static matrix Rtmat;
81 static double *Rvector; /* residuals at each gradient */
82 static double *tempRvector; /* copy of residuals at each gradient */
83 static int *B0list = NULL;
84 static matrix Fmatrix;
85 static matrix Dmatrix;
86 static matrix OldD;
87 static matrix Hplusmatrix, Hminusmatrix;
88 static vector Dvector;
89 static matrix tempFmatrix[2];
90 static matrix tempDmatrix[2];
91 static matrix tempHplusmatrix[2], tempHminusmatrix[2];
92 /* static vector tempDvector; */
93
94 static byte *maskptr = NULL;
95 static double eigs[12];
96 static double deltatau;
97 static double *wtfactor; /* weight factors for time points at
98 each voxel */
99 static double *bmatrix; /* b matrix = GiGj for gradient
100 intensities */
101 static double *cumulativewt; /* overall wt. factor for each
102 gradient */
103 static long rewtvoxels; /* how many voxels were reweighted */
104 static double sigma; /* std.deviation */
105 static double ED; /* error for each iteration - cost
106 function result */
107 static int automask = 0; /* automasking flag - user option */
108 static int reweight_flag = 0; /* reweight computation flag - user
109 option */
110 static int method = -1; /* linear or non-linear method - user
111 option */
112 static int max_iter = -2; /* maximum #convergence iteration
113 steps - user option */
114 static int max_iter_rw = -2; /* maximum #convergence iteration
115 steps - user option */
116 static int eigs_flag = 0; /* eigenvalue calculation in output -
117 user option */
118 static int cumulative_flag = 0; /* calculate, display cumulative wts
119 for gradients - user option */
120 static int debug_briks = 0; /* put Ed, Ed0 and Converge step
121 sub-briks in output - user
122 option */
123 static int verbose = 0; /* print out info every verbose number
124 of voxels - user option */
125 static int afnitalk_flag = 0; /* show convergence in AFNI graph -
126 user option */
127 static int bmatrix_given = 0; /* user input file is b matrix
128 (instead of gradient direction
129 matrix) with 6 columns: Bxx, Byy,
130 Bzz, Bxy, Bxz, Byz */
131 static int BMAT_NZ = 0; /* non-zero bmatrix supplied */
132 static int use_mean_b0 = 0; /* compute initial linear estimate
133 with first b=0 volume rather than
134 mean of all b=0 volumes */
135 static int opt_method = 2; /* use gradient descent instead of
136 Powell's new optimize method*/
137 static int voxel_opt_method = 0;/* hybridize optimization between
138 Powell and gradient descent */
139 static int Powell_npts = 1; /* number of points in input dataset
140 for Powell optimization function */
141 static float *Powell_ts; /* pointer to time-wise voxel data for
142 Powell optimization function */
143 static double Powell_J;
144 static double backoff_factor = 0.2; /* minimum allowable factor for
145 lambda2,3 relative to lambda1
146 eigenvalues*/
147 static float csf_val = 1.0; /* a default value for diffusivity for
148 CSF at 37C; apr,2016: now set to
149 '1' */
150 static float rd_val = 1.0; /* a default value for diffusivity for
151 CSF at 37C; apr,2016: now set to
152 '1' */
153
154 static float csf_fa = 0.012345678; /* default FA value where there is
155 CSF */
156 static float MAX_BVAL = -1.; /* use in scaling csf_val for MD and
157 such, *if* user doesn't enter
158 one */
159 int USER_CSF_VAL = 0; // flag so we know if the user entered
160 // his/her own scaling
161 static float MIN_BAD_SCALE_MD = 100.; // May,2016: mult this number
162 // times csf_val to detect some
163 // baaadness
164 float SCALE_VAL_OUT = -1.0 ; // allow users to scaled physical
165 // values by 1000 easily; will be
166 // set to 1 if negative after
167 // reading in inputs
168 float BMAX_REF = 0.01; // for identifying reference bvalues
169
170 // [PT: May, 2017]
171 static float *I0_ptr; /* I0 matrix */
172 static int nvols = -1;
173 // used in copying DT info: NOTE THE ORDER!
174 static int sublist[7]={6,0,1,2,3,4,5};
175
176 static NI_stream_type * DWIstreamid = 0; /* NIML stream ID */
177
178
179 int ChiSq_GOF( THD_3dim_dataset *DWmeas,
180 THD_3dim_dataset *DWfit,
181 int *b0list,
182 float **ccc,
183 byte *maskp);
184
185 // *Shamelessly* copied+pasted from 3dDTtoDWI.c for use here, because
186 // *it was a pain to navigate around static variables and linking.
187 static void DTtoDWI_tsfunc (double tzero, double tdelta, int npts, float ts[],
188 double ts_mean, double ts_slope, void *ud,
189 int nbriks, float *val);
190 static void Form_R_Matrix (MRI_IMAGE * grad1Dptr);
191 static void DWItoDT_tsfunc ( double tzero, double tdelta, int npts,
192 float ts[], double ts_mean, double ts_slope,
193 void *ud, int nbriks, float *val);
194 static void EIG_func (void);
195 static float Calc_FA(float *val);
196 static float Calc_MD(float *val);
197 static float Calc_RD(float *val); // pt,Oct,2016: for RD
198 static void ComputeD0 (void);
199 static double ComputeJ (float ts[], int npts);
200 static void ComputeDeltaTau (void);
201 static void Computebmatrix (MRI_IMAGE * grad1Dptr, int NO_ZERO_ROW1);
202 static void InitGlobals (int npts);
203 static void FreeGlobals (void);
204 static void Store_Computations (int i, int npts, int converge_step);
205 static void Restore_Computations (int i, int npts, int converge_step);
206 static void InitWtfactors (int npts);
207 static void ComputeWtfactors (int npts);
208 static void ComputeHpHm (double deltatau);
209 static void ComputeNewD (void);
210 static int TestConvergence (matrix NewD, matrix OldD);
211 static void udmatrix_to_vector (matrix m, vector * v);
212 static void udmatrix_copy (double *udptr, matrix * m);
213 static double matrix_sumabs (matrix m);
214 static double *InvertSym3 (double a, double b, double c, double e, double f,
215 double i);
216 static void matrix_copy (matrix a, matrix * b);
217 static int DWI_Open_NIML_stream(void);
218 static int DWI_NIML_create_graph(void);
219 static void DWI_AFNI_update_graph(double *Edgraph, double *dtau, int npts);
220 static void vals_to_NIFTI(float *val);
221 static void Save_Sep_DTdata(THD_3dim_dataset *, char *, int);
222 static void Copy_dset_array(THD_3dim_dataset *, int,int,char *, int);
223 static int ComputeDwithPowell(float *ts, float *val, int npts, int nbriks);
224 static int bad_DWI_data(int npts, float *ts);
225 static int bad_DWI_data_2_MD(float v1, float v2, float v3); // May,2016: new badness
226 // check, based on MD val
227 static double compute_mean_B0(int npts, float *ts, char B0flag);
228 static void Assign_CSF_values(float *val);
229 static int all_dwi_zero(int npts, float *ts);
230
231 int
main(int argc,char * argv[])232 main (int argc, char *argv[])
233 {
234 THD_3dim_dataset *old_dset=NULL, *new_dset=NULL; // input and output datasets
235 THD_3dim_dataset *fit_dset=NULL; // intermed dset: DWIs on tensor
236 THD_3dim_dataset *new_dsetDT=NULL, *old_dset0=NULL;
237 int nopt, nbriks;
238 int i, eigs_brik;
239 MRI_IMAGE *grad1Dptr = NULL;
240 MRI_IMAGE *anat_im = NULL;
241 #if 0
242 int nvox;
243 short *sar = NULL;
244 short *tempsptr = NULL;
245 byte *tempbptr = NULL;
246 short tempval;
247 #endif
248
249 double *cumulativewtptr=NULL;
250 int mmvox=0 ;
251 int nxyz;
252 int sep_dsets = 0;
253
254 // [PT: May, 2017]
255 MRI_IMAGE *data_im = NULL;
256 double fac;
257 int n;
258
259 // [PT: Aug, 2017]
260 float tmp_cwvalue=0.0; // for outputting cumulative_wt values
261 int ncwts=0;
262 FILE *fout1=NULL;
263
264
265 /*----- Read command line -----*/
266 if (argc < 2 || strcmp (argv[1], "-help") == 0)
267 {
268 printf (
269 "\n"
270 " Usage: 3dDWItoDT [options] gradient-file dataset\n"
271 "\n"
272 " Computes 6 principle direction tensors from multiple gradient vectors\n"
273 " and corresponding DTI image volumes.\n"
274 " The program takes two parameters as input : \n"
275 " a 1D file of the gradient vectors with lines of ASCII floats:\n"
276 " Gxi, Gyi, Gzi.\n"
277 " Only the non-zero gradient vectors are included in this file (no G0 \n"
278 " line). \n"
279 " ** Now, a '1D' file of b-matrix elements can alternatively be input,\n"
280 " and *all* the gradient values are included!**\n"
281 " A 3D bucket dataset with Np+1 sub-briks where the first sub-brik is the\n"
282 " volume acquired with no diffusion weighting.\n"
283 "\n"
284 " OUTPUTS: \n"
285 " + you can output all 6 of the independent tensor values (Dxx, Dyy, \n"
286 " etc.), as well as all three eigenvalues (L1, L2, L3) and \n"
287 " eigenvectors (V1, V2, V3), and useful DTI parameters FA, MD and\n"
288 " RD.\n"
289 " + 'Debugging bricks' can also be output, see below.\n"
290 "\n"
291 " Options:\n"
292 " -prefix pname = Use 'pname' for the output dataset prefix name.\n"
293 " [default='DT']\n\n"
294 " -automask = mask dataset so that the tensors are computed only for\n"
295 " high-intensity (presumably brain) voxels. The intensity \n"
296 " level is determined the same way that 3dClipLevel works.\n\n"
297 " -mask dset = use dset as mask to include/exclude voxels\n\n"
298 " -bmatrix_NZ FF = switch to note that the input dataset is b-matrix, \n"
299 " not gradient directions, and there is *no* row of zeros \n"
300 " at the top of the file, similar to the format for the grad\n"
301 " input: N-1 rows in this file for N vols in matched data set.\n"
302 " There must be 6 columns of data, representing either elements\n"
303 " of G_{ij} = g_i*g_j (i.e., dyad of gradients, without b-value\n"
304 " included) or of the DW scaled version, B_{ij} = b*g_i*g_j.\n"
305 " The order of components is: G_xx G_yy G_zz G_xy G_xz G_yz.\n"
306 " -bmatrix_Z FF = similar to '-bmatrix_NZ' above, but assumes that first\n"
307 " row of the file is all zeros (or whatever the b-value for\n"
308 " the reference volume was!), i.e. there are N rows to the\n"
309 " text file and N volumes in the matched data set.\n\n"
310 " -bmatrix_FULL FF = exact same as '-bmatrix_Z FF' above (i.e. there are N\n"
311 " rows to the text file and N volumes in the matched data set)\n"
312 " with just a lot more commonsensical name. Definitely would\n"
313 " be preferred way to go, for ease of usage!\n"
314 "\n"
315 " -scale_out_1000 = increase output parameters that have physical units\n"
316 " (DT, MD, RD, L1, L2 and L3) by multiplying them by 1000. This\n"
317 " might be convenient, as the input bmatrix/gradient values \n"
318 " can have their physical magnitudes of ~1000 s/mm^2, for\n"
319 " which typical adult WM has diffusion values of MD~0.0007\n"
320 " (in physical units of mm^2/s), and people might not like so\n"
321 " many decimal points output; using this option rescales the\n"
322 " input b-values and would lead to having a typical MD~0.7\n"
323 " (now in units of x10^{-3} mm^2/s). If you are not using\n"
324 " bmatrix/gradient values that have their physical scalings,\n"
325 " then using this switch probably wouldn't make much sense.\n"
326 " FA, V1, V2 and V3 are unchanged.\n\n"
327 " -bmax_ref THR = if the 'reference' bvalue is actually >0, you can flag\n"
328 " that here. Otherwise, it is assumed to be zero.\n"
329 " At present, this is probably only useful/meaningful if\n"
330 " using the '-bmatrix_Z ...' or '-bmatrix_FULL ...' \n"
331 " option, where the reference bvalue must be found and \n"
332 " identified from the input info alone.\n"
333 "\n"
334 " -nonlinear = compute iterative solution to avoid negative eigenvalues.\n"
335 " This is the default method.\n\n"
336 " -linear = compute simple linear solution.\n\n"
337 " -reweight = recompute weight factors at end of iterations and restart\n\n"
338 " -max_iter n = maximum number of iterations for convergence (Default=10).\n"
339 " Values can range from -1 to any positive integer less than\n"
340 " 101. A value of -1 is equivalent to the linear solution.\n"
341 " A value of 0 results in only the initial estimate of the\n"
342 " diffusion tensor solution adjusted to avoid negative\n"
343 " eigenvalues.\n\n"
344 " -max_iter_rw n = max number of iterations after reweighting (Default=5)\n"
345 " values can range from 1 to any positive integer less\n"
346 " than 101.\n\n"
347 " -eigs = compute eigenvalues, eigenvectors, fractional anisotropy and mean\n"
348 " diffusivity in sub-briks 6-19. Computed as in 3dDTeig\n\n"
349 " -debug_briks = add sub-briks with Ed (error functional), Ed0 (orig.\n"
350 " error), number of steps to convergence and I0 (modeled B0\n"
351 " volume).\n"
352 " [May, 2017] This also now calculates two goodness-of-fit\n"
353 " measures and outputs a new PREFIX_CHI* dset that has two\n"
354 " briks:\n"
355 " brik [0]: chi^2_p,\n"
356 " brik [1]: chi^2_c.\n"
357 " These values are essentially calculated according to\n"
358 " Papadakis et al. (2003, JMRI), Eqs. 4 and 3,\n"
359 " respectively (in chi^2_c, the sigma value is the\n"
360 " variance of measured DWIs *per voxel*). Note for both\n"
361 " chi* values, only DWI signal values are used in the\n"
362 " calculation (i.e., where b>THR; by default,\n"
363 " THR=0.01, which can be changed using '-bmax_ref ...').\n"
364 " In general, chi^2_p values seem to be <<1, consistent\n"
365 " with Papadakis et al.'s Fig. 4; the chi^2_c values are\n"
366 " are also pretty consistent with the same fig and seem to\n"
367 " be best viewed with the upper limit being roughly =Ndwi\n"
368 " or =Ndwi-7 (with the latter being the given degrees\n"
369 " of freedom value by Papadakis et al.)\n"
370 " -cumulative_wts = show overall weight factors for each gradient level\n"
371 " May be useful as a quality control\n\n"
372 " -verbose nnnnn = print convergence steps every nnnnn voxels that survive\n"
373 " to convergence loops (can be quite lengthy).\n\n"
374 " -drive_afni nnnnn = show convergence graphs every nnnnn voxels that\n"
375 " survive to convergence loops. AFNI must have NIML\n"
376 " communications on (afni -niml)\n\n"
377 " -sep_dsets = save tensor, eigenvalues, vectors, FA, MD in separate\n"
378 " datasets\n\n"
379 " -csf_val n.nnn = assign diffusivity value to DWI data where the mean\n"
380 " values for b=0 volumes is less than the mean of the\n"
381 " remaining volumes at each voxel. The default value is\n"
382 " '1.0 divided by the max bvalue in the grads/bmatrices'.\n"
383 " The assumption is that there are flow artifacts in CSF\n"
384 " and blood vessels that give rise to lower b=0 voxels.\n"
385 " NB: MD, RD L1, L2, L3, Dxx, Dyy, etc. values are all\n"
386 " scaled in the same way.\n\n"
387 " -min_bad_md N = change the min MD value used as a 'badness check' for\n"
388 " tensor fits that have veeery (-> unreasonably) large MD\n"
389 " values. Voxels where MD > N*(csf_val) will be treated\n"
390 " like CSF and turned into spheres with radius csf_val \n"
391 " (default N=100).\n"
392 " -csf_fa n.nnn = assign a specific FA value to those voxels described\n"
393 " above The default is 0.012345678 for use in tractography\n"
394 " programs that may make special use of these voxels\n\n"
395 " -opt mname = if mname is 'powell', use Powell's 2004 method for \n"
396 " optimization. If mname is 'gradient' use gradient descent\n"
397 " method. If mname is 'hybrid', use combination of methods.\n"
398 " MJD Powell, \"The NEWUOA software for unconstrained \n"
399 " optimization without derivatives\", Technical report DAMTP\n"
400 " 2004/NA08, Cambridge University Numerical Analysis Group:\n"
401 " See: http://www.ii.uib.no/~lennart/drgrad/Powell2004.pdf\n\n"
402 " -mean_b0 = use mean of all b=0 volumes for linear computation and initial\n"
403 " linear for nonlinear method\n\n"
404 " Example:\n"
405 " 3dDWItoDT -prefix rw01 -automask -reweight -max_iter 10 \\\n"
406 " -max_iter_rw 10 tensor25.1D grad02+orig.\n\n"
407 " The output is a 6 sub-brick bucket dataset containing \n"
408 " Dxx, Dxy, Dyy, Dxz, Dyz, Dzz\n"
409 " (the lower triangular, row-wise elements of the tensor in symmetric matrix\n"
410 " form). Additional sub-briks may be appended with the -eigs and -debug_briks\n"
411 " options. These results are appropriate as the input to 3dDTeig.\n"
412 "\n");
413 printf ("\n" MASTER_SHORTHELP_STRING);
414 exit (0);
415 }
416
417 mainENTRY ("3dDWItoDT main");
418 machdep ();
419 AFNI_logger ("3dDWItoDT", argc, argv);
420 PRINT_VERSION("3dDWItoDT") ; AUTHOR("Daniel Glen") ;
421 THD_check_AFNI_version("3dDWItoDT") ;
422
423 nopt = 1;
424 nbriks = 6; /* output contains 6 sub-briks by default */
425 method = -1;
426 reweight_flag = 0;
427
428 datum = MRI_float;
429 while (nopt < argc && argv[nopt][0] == '-')
430 {
431
432 /*-- prefix --*/
433
434 if (strcmp (argv[nopt], "-prefix") == 0)
435 {
436 if (++nopt >= argc)
437 {
438 ERROR_exit("Error - prefix needs an argument!");
439 }
440 /* change name from default prefix */
441 MCW_strncpy (prefix, argv[nopt], THD_MAX_PREFIX);
442 /* check file name to be sure not to overwrite - mod
443 drg 12/9/2004 */
444 if (!THD_filename_ok (prefix))
445 {
446 ERROR_exit("Error - %s is not a valid prefix!", prefix);
447 }
448 nopt++;
449 continue;
450 }
451
452 /*-- datum --*/
453
454 if (strcmp (argv[nopt], "-datum") == 0)
455 {
456 if (++nopt >= argc)
457 {
458 ERROR_exit("Error - datum needs an argument!");
459 }
460 if (strcmp (argv[nopt], "short") == 0)
461 {
462 datum = MRI_short;
463 }
464 else if (strcmp (argv[nopt], "float") == 0)
465 {
466 datum = MRI_float;
467 }
468 else if (strcmp (argv[nopt], "byte") == 0)
469 {
470 datum = MRI_byte;
471 }
472 else
473 {
474 ERROR_exit("-datum of type '%s' is not supported!",
475 argv[nopt]);
476 }
477 nopt++;
478 continue;
479 }
480 if (strcmp (argv[nopt], "-automask") == 0)
481 {
482 if(maskptr != NULL){
483 ERROR_exit("ERROR: can't use -mask with -automask!");
484 }
485 automask = 1;
486 nopt++;
487 continue;
488 }
489
490 if( strcmp(argv[nopt],"-mask") == 0 ){
491 THD_3dim_dataset * mask_dset ;
492 if( automask ){
493 ERROR_exit("ERROR: can't use -mask with -automask!");
494 }
495 mask_dset = THD_open_dataset(argv[++nopt]) ;
496 CHECK_OPEN_ERROR(mask_dset,argv[nopt]) ;
497 if( maskptr != NULL ){
498 ERROR_exit("** ERROR: can't have 2 -mask options!");
499 }
500 maskptr = THD_makemask( mask_dset , 0 , 1.0,-1.0 ) ;
501 mmvox = DSET_NVOX( mask_dset ) ;
502
503 DSET_delete(mask_dset) ; nopt++ ; continue ;
504 }
505
506 // input bmatrix with initial B=0 line
507 if (strcmp(argv[nopt], "-bmatrix_Z") == 0) {
508 bmatrix_given = 1;
509 nopt++;
510 continue;
511 }
512 // duplicate behavior of above, but with a LOT easier name
513 if (strcmp(argv[nopt], "-bmatrix_FULL") == 0) {
514 bmatrix_given = 1;
515 nopt++;
516 continue;
517 }
518
519 // input bmatrix without first B=0 line
520 if (strcmp(argv[nopt], "-bmatrix_NZ") == 0) {
521 BMAT_NZ = 1;
522 nopt++;
523 continue;
524 }
525
526 // compute mean across b=0 vols for init. linear
527 if (strcmp(argv[nopt], "-mean_b0") == 0) {
528 use_mean_b0 = 1;
529 nopt++;
530 continue;
531 }
532
533 // can mult output diff values by 1000
534 if (strcmp(argv[nopt], "-scale_out_1000") == 0) {
535 SCALE_VAL_OUT = 1000.;
536 nopt++;
537 continue;
538 }
539
540 if (strcmp (argv[nopt], "-bmax_ref") == 0){
541 if(++nopt >=argc )
542 ERROR_exit("Error: need an argument after -bmax_ref!");
543 BMAX_REF = (float) strtod(argv[nopt], NULL);
544 nopt++;
545 continue;
546 }
547
548 // May,2016: essentially, turn off badness criterion of
549 // superlarge MD by setting the value SOOO high
550 if (strcmp (argv[nopt], "-min_bad_md") == 0){
551 if(++nopt >=argc )
552 ERROR_exit("Error - need an argument after -min_bad_md!");
553 MIN_BAD_SCALE_MD = (float) strtod(argv[nopt], NULL);
554 if(MIN_BAD_SCALE_MD < 0.00001)
555 ERROR_exit("Bad MD scale! Needs to be >0.00001-- "
556 " and probably much more so!");
557 nopt++;
558 continue;
559 }
560
561
562 if (strcmp (argv[nopt], "-linear") == 0)
563 {
564 if(method==1)
565 {
566 ERROR_exit("Error - cannot select both linear and "
567 "non-linear methods at the same time");
568 }
569 method = 0;
570 nopt++;
571 continue;
572 }
573
574 if ((strcmp (argv[nopt], "-nonlinear") == 0) ||
575 (strcmp (argv[nopt], "-non-linear") == 0))
576 {
577 if(method==0)
578 {
579 ERROR_exit("Error - cannot select both linear "
580 "and non-linear methods at the same time");
581 exit(1);
582 }
583 method = 1;
584 nopt++;
585 continue;
586 }
587 if (strcmp (argv[nopt], "-reweight") == 0) {
588 reweight_flag = 1;
589 nopt++;
590 continue;
591 }
592
593 if (strcmp (argv[nopt], "-max_iter") == 0)
594 {
595 if(++nopt >=argc ){
596 ERROR_exit("Error - need an argument after -max_iter!");
597 }
598 max_iter = strtol(argv[nopt], NULL, 10);
599 if ((max_iter <-1)||(max_iter>100)) {
600 ERROR_exit("Error - max_iter must be between -1 and 100");
601 }
602 nopt++;
603 continue;
604 }
605
606 if (strcmp (argv[nopt], "-max_iter_rw") == 0)
607 {
608 if(++nopt >=argc ){
609 ERROR_exit("Error - need an argument after -max_iter_rw!");
610 }
611 max_iter_rw = strtol(argv[nopt], NULL, 10);
612 if ((max_iter_rw <=0)||(max_iter_rw>100)) {
613 ERROR_exit("Error - max_iter_rw must be between 1 and 100");
614 }
615 nopt++;
616 continue;
617 }
618
619 if (strcmp (argv[nopt], "-eigs") == 0)
620 {
621 eigs_flag = 1;
622 nopt++;
623 continue;
624 }
625
626 if (strcmp (argv[nopt], "-cumulative_wts") == 0)
627 {
628 cumulative_flag = 1;
629 nopt++;
630 continue;
631 }
632
633 if ((strcmp (argv[nopt], "-debug_briks") == 0) ||
634 (strcmp (argv[nopt], "-debug_bricks") == 0))
635 {
636 debug_briks = 1;
637 nopt++;
638 continue;
639 }
640
641 if (strcmp (argv[nopt], "-verbose") == 0)
642 {
643 if(++nopt >=argc ){
644 ERROR_exit("*** Error - need an argument after -verbose!");
645 }
646 verbose = strtol(argv[nopt], NULL, 10);
647 if (verbose<=0) {
648 ERROR_exit("Error- verbose steps must be a positive number!");
649 }
650 nopt++;
651 continue;
652 }
653
654 if (strcmp (argv[nopt], "-drive_afni") == 0)
655 {
656 if(++nopt >=argc ){
657 ERROR_exit("Error - need an argument after -drive_afni!");
658 }
659 afnitalk_flag = strtol(argv[nopt], NULL, 10);
660 if (afnitalk_flag<=0) {
661 ERROR_exit("Error - drive_afni steps must be "
662 "a positive number!");
663 }
664 nopt++;
665 continue;
666 }
667
668 if (strcmp (argv[nopt], "-sep_dsets") == 0)
669 {
670 sep_dsets = 1; /* save data in separate datasets */
671 nopt++;
672 continue;
673 }
674
675 if (strcmp (argv[nopt], "-opt") == 0)
676 {
677 if (++nopt >= argc)
678 {
679 ERROR_exit("Error - opt should be followed by "
680 "gradient or powell!");
681 }
682 if (strcmp(argv[nopt], "gradient") == 0)
683 {
684 opt_method = 0;
685 }
686 else if (strcmp(argv[nopt], "powell") == 0)
687 {
688 opt_method = 1; /* use Powell's new optimize
689 method instead of gradient
690 descent*/
691 }
692 else if (strcmp(argv[nopt], "hybrid") == 0)
693 {
694 opt_method = 2; /* use combination of Powell and
695 gradient descent*/
696 }
697 else
698 {
699 ERROR_exit("-opt method '%s' is not supported!",
700 argv[nopt]);
701 }
702 nopt++;
703 continue;
704 }
705
706 if (strcmp (argv[nopt], "-backoff") == 0)
707 {
708 backoff_factor = strtod(argv[++nopt],NULL);
709 if(backoff_factor<=0.0 || backoff_factor>1.0)
710 ERROR_exit("-backoff factor must be > 0 and <= 1");
711 nopt++;
712 continue;
713 }
714
715 if (strcmp (argv[nopt], "-csf_val") == 0){
716 if(++nopt >=argc )
717 ERROR_exit("Error - need an argument after -csf_val!");
718 csf_val = (float) strtod(argv[nopt], NULL);
719 rd_val = csf_val; // pt,Oct,2016: for RD
720 USER_CSF_VAL = 1;
721 nopt++;
722 continue;
723 }
724
725 if (strcmp (argv[nopt], "-csf_fa") == 0){
726 if(++nopt >=argc )
727 ERROR_exit("Error - need an argument after -csf_fa!");
728 csf_fa = (float) strtod(argv[nopt], NULL);
729 nopt++;
730 continue;
731 }
732
733 ERROR_exit("Error - unknown option %s", argv[nopt]);
734 }
735
736 if(method==-1)
737 method = 1; /* if not selected, choose non-linear method for now */
738
739 if(max_iter>=-1){
740 if(method==0)
741 WARNING_message("Warning - max_iter will be ignored "
742 "for linear methods");
743 }
744 else
745 max_iter = MAX_CONVERGE_STEPS;
746
747 if(max_iter_rw>=0) {
748 if(method==0)
749 WARNING_message("Warning - max_iter_rw will be ignored "
750 "for linear methods");
751 if(reweight_flag==0)
752 WARNING_message("Warning - max_iter_rw will be ignored "
753 "when not reweighting");
754 }
755 else
756 max_iter_rw = MAX_RWCONVERGE_STEPS;
757
758 if((method==0)&&(reweight_flag==1)) {
759 WARNING_message("Warning - can not reweight voxels for linear method");
760 reweight_flag = 0;
761 }
762
763 if(cumulative_flag==1) {
764 if(method==0) {
765 WARNING_message("Warning - can not compute cumulative weights "
766 "for linear method");
767 cumulative_flag = 0;
768 }
769 if(reweight_flag == 0) {
770 WARNING_message("Warning - can not compute cumulative weights "
771 "if not reweighting");
772 cumulative_flag = 0;
773 }
774 }
775
776 if((method==0)&&(debug_briks==1)) {
777 WARNING_message("Warning - can not compute debug sub-briks "
778 "for linear method");
779 debug_briks = 0;
780 }
781
782 if((method==0)&&(afnitalk_flag>0)) {
783 WARNING_message("Warning - can not graph convergence in AFNI "
784 "for linear method");
785 afnitalk_flag = 0;
786 }
787
788 if(eigs_flag)
789 nbriks += 15; // pt,Oct,2016: for RD
790
791 if(debug_briks)
792 nbriks += 4;
793
794
795 /*----- read input datasets -----*/
796
797 if (nopt >= argc)
798 {
799 ERROR_exit("Error - No input dataset!?");
800 }
801
802 /* first input dataset - should be gradient vector file of ascii
803 floats Gx,Gy,Gz */
804
805 /* read gradient vector 1D file */
806 grad1Dptr = mri_read_1D (argv[nopt]);
807 if (grad1Dptr == NULL)
808 {
809 ERROR_exit("Error reading gradient vector file");
810 }
811
812 // if ((grad1Dptr->ny != 3 && !bmatrix_given) || (grad1Dptr->ny !=
813 // 6 && bmatrix_given)) // $$old
814 if ((grad1Dptr->ny != 3 && !(bmatrix_given || BMAT_NZ) )
815 || (grad1Dptr->ny != 6 && (bmatrix_given || BMAT_NZ ) )) // BMAT_NZ opts
816 {
817 ERROR_message("Error - Only 3 columns of gradient vectors "
818 "(or 6 columns for b matrices) "
819 "allowed: %d columns found", grad1Dptr->ny);
820 mri_free (grad1Dptr);
821 exit (1);
822 }
823
824 if (grad1Dptr->nx < 6)
825 {
826 mri_free (grad1Dptr);
827 ERROR_message("Error - Must have at least 6 gradient vectors");
828 ERROR_exit("%d columns found", grad1Dptr->nx);
829 }
830
831 // -----------------------------------------------------------------
832 // apply 1/SCALE_VAL_OUT to the bvalue as a way to scale;
833 // shamelessly pillaging format of mri_histog.c; later might allow
834 // users to choose other values
835 if( SCALE_VAL_OUT < 0 )
836 SCALE_VAL_OUT = 1.;
837 else
838 INFO_message("Implementing user-selected scaling of diffusion values by %.0f",
839 SCALE_VAL_OUT);
840
841 if( grad1Dptr->kind != MRI_float )
842 ERROR_exit("How are the gradients not floats?");
843
844 int npix;
845 register float *flar=mri_data_pointer(grad1Dptr);
846 npix = grad1Dptr->nvox;
847 for( i=0 ; i<npix ; i++ ){
848 flar[i]/= SCALE_VAL_OUT;
849 }
850 BMAX_REF/= SCALE_VAL_OUT;
851
852 //fprintf(stderr,"\n\n nyx=%d; npix=%d\n",grad1Dptr->nxy, npix);
853 //for( i=0 ; i<npix ; i++)
854 // fprintf(stderr,"\n %f", flar[i]);
855 // -----------------------------------------------------------------
856
857 Form_R_Matrix (grad1Dptr); /* use grad1Dptr to compute R matrix */
858
859 nopt++;
860
861 /* Now read in all the MRI volumes for each gradient vector */
862 /* assumes first one is no gradient */
863 old_dset = THD_open_dataset (argv[nopt]);
864 CHECK_OPEN_ERROR(old_dset,argv[nopt]) ;
865
866 /* expect at least 7 values per voxel - 7 sub-briks as input dataset */
867 if (!(!bmatrix_given && DSET_NVALS (old_dset) == (grad1Dptr->nx + 1))
868 && !((bmatrix_given && DSET_NVALS (old_dset) == (grad1Dptr->nx)))
869 && !((BMAT_NZ && DSET_NVALS (old_dset) == (grad1Dptr->nx))) )
870 {
871 mri_free (grad1Dptr);
872 ERROR_message("Error - Dataset must have number of sub-briks "
873 "equal to one more than number");
874 ERROR_exit(" of gradient vectors (B0+Bi)!");
875 }
876 nxyz = DSET_NVOX(old_dset) ;
877 if( maskptr != NULL && mmvox != nxyz ){
878 ERROR_exit("Mask and input datasets not the same size!") ;
879 }
880
881 nvols = DSET_NVALS (old_dset);
882
883 if (bmatrix_given) {
884 InitGlobals (grad1Dptr->nx); /* initialize all the matrs and vecs */
885 }
886 else { // $$same: can leave same, because BMAT_NZ will have `nx+1'
887 InitGlobals (grad1Dptr->nx + 1); /* initialize all the matrs and vecs */
888 }
889
890 Computebmatrix (grad1Dptr, BMAT_NZ); /* compute bij=GiGj */
891 // after this, the MAX_BVAL is known, as well; apply it to the CSF
892 // val *if* the user didn't set his/her own scale. apr,2016
893
894 if(MAX_BVAL<=0)
895 ERROR_exit("Whooa! How did the max bval get to look like %f??", MAX_BVAL);
896
897 if (!bmatrix_given)
898 INFO_message("The maximum magnitude of the gradients "
899 "appears to be: %.5f", MAX_BVAL);
900 else
901 INFO_message("The maximum magnitude of the bmatrix "
902 "appears to be: %.2f", MAX_BVAL);
903 if(!USER_CSF_VAL) {
904 csf_val/= MAX_BVAL;
905 rd_val/= MAX_BVAL;
906 INFO_message("-> and, by scale, the 'CSF value' (e.g., MD)"
907 " for any unfit voxels will be: %.8f", csf_val);
908 }
909 INFO_message("Voxels with MD>%f will be deemed unfit for "
910 "fitting and given 'CSF' status",
911 MIN_BAD_SCALE_MD*csf_val);
912
913 if (automask)
914 {
915 DSET_mallocize (old_dset);
916 DSET_load (old_dset); /* get B0 (anatomical image) from
917 dataset */
918 /*anat_im = THD_extract_float_brick( 0, old_dset ); */
919 anat_im = DSET_BRICK (old_dset, 0); /* set pointer to the 0th
920 sub-brik of the
921 dataset */
922 maskptr = mri_automask_image (anat_im); /* maskptr is a
923 byte pointer for
924 volume */
925 }
926
927 /* temporarily set artificial timing to 1 second interval */
928 EDIT_dset_items (old_dset,
929 ADN_ntt, DSET_NVALS (old_dset),
930 ADN_ttorg, 0.0,
931 ADN_ttdel, 1.0, ADN_tunits,
932 UNITS_SEC_TYPE, NULL);
933
934 if(afnitalk_flag) {
935 if(DWI_Open_NIML_stream()!=0) { /* Open NIML stream */
936 afnitalk_flag = 0;
937 WARNING_message("Could not open NIML communications with AFNI");
938 }
939 else
940 if(DWI_NIML_create_graph()!=0) {
941 afnitalk_flag = 0;
942 WARNING_message("Could not create graph within AFNI");
943 /* Close NIML stream */
944 NI_stream_close(DWIstreamid);
945 DWIstreamid = 0;
946 }
947 }
948
949 /*------------- ready to compute new dataset -----------*/
950
951 new_dset = MAKER_4D_to_typed_fbuc (old_dset, /* input dataset */
952 prefix, /* output prefix */
953 datum, /* output datum */
954 0, /* ignore count */
955 0, /* can't detrend in maker
956 function KRH 12/02 */
957 nbriks, /* number of briks */
958 DWItoDT_tsfunc, /* timeseries processor */
959 NULL, /* data for tsfunc */
960 NULL, /* mask */
961 0 /* Allow auto scaling of output */
962 );
963
964 if(afnitalk_flag && (DWIstreamid!=0)) {
965 /* Close NIML stream */
966 NI_stream_close(DWIstreamid);
967 }
968
969 if(cumulative_flag && reweight_flag) {
970
971 // [PT: Aug, 2017] Functionality to output cumulative weights to
972 // text file
973 char cwprefix[THD_MAX_PREFIX], cprefix[THD_MAX_PREFIX];
974 char *ext, nullch;
975 sprintf(cprefix,"%s", prefix);
976 if(has_known_non_afni_extension(prefix)){ /* for NIFTI, 3D, Niml,
977 Analyze,...*/
978 ext = find_filename_extension(prefix);
979 cprefix[strlen(prefix) - strlen(ext)] = '\0'; /* remove
980 non-afni-extension
981 for now*/
982 }
983 else {
984 nullch = '\0';
985 ext = &nullch;
986 }
987 sprintf(cwprefix,"%s_cwts.1D", cprefix);
988
989 if( (fout1 = fopen(cwprefix, "w")) == NULL) {
990 fprintf(stderr, "Error opening file %s.",cwprefix);
991 exit(11);
992 }
993
994 cumulativewtptr = cumulativewt;
995 INFO_message("Cumulative Wt. factors (and output to %s:", cwprefix);
996 // need to get the number based on whether full mat was input
997 if(bmatrix_given)
998 ncwts = grad1Dptr->nx;
999 else
1000 ncwts = grad1Dptr->nx + 1;
1001 for(i=0; i<ncwts; i++){
1002 *cumulativewtptr = *cumulativewtptr / rewtvoxels;
1003 tmp_cwvalue = *cumulativewtptr++;
1004 INFO_message("%5.3f ", tmp_cwvalue);
1005 fprintf(fout1,"%5.3f\n",tmp_cwvalue);
1006 }
1007 fclose(fout1);
1008
1009 /* printf("\n");*/
1010 }
1011
1012
1013 if (new_dset != NULL)
1014 {
1015 // copy the tensor over. Ugh.
1016 if( debug_briks) {
1017 new_dsetDT = THD_copy_dset_subs(new_dset, sublist);
1018 old_dset0 = THD_copy_one_sub( new_dset, nbriks-1);
1019 }
1020
1021 tross_Copy_History (old_dset, new_dset);
1022 // Warning for future generations: carefully note the order
1023 // of elements here for the DT!!!
1024 EDIT_dset_items (new_dset, ADN_brick_label_one + 0, "Dxx", ADN_none);
1025 EDIT_dset_items (new_dset, ADN_brick_label_one + 1, "Dxy", ADN_none);
1026 EDIT_dset_items (new_dset, ADN_brick_label_one + 2, "Dyy", ADN_none);
1027 EDIT_dset_items (new_dset, ADN_brick_label_one + 3, "Dxz", ADN_none);
1028 EDIT_dset_items (new_dset, ADN_brick_label_one + 4, "Dyz", ADN_none);
1029 EDIT_dset_items (new_dset, ADN_brick_label_one + 5, "Dzz", ADN_none);
1030 if(eigs_flag) {
1031 eigs_brik = ADN_brick_label_one + 6; /* 1st eigenvalue brik */
1032 EDIT_dset_items(new_dset, eigs_brik+0, "lambda_1", ADN_none);
1033 EDIT_dset_items(new_dset, eigs_brik+1, "lambda_2", ADN_none);
1034 EDIT_dset_items(new_dset, eigs_brik+2, "lambda_3", ADN_none);
1035 EDIT_dset_items(new_dset, eigs_brik+3, "eigvec_1[1]", ADN_none);
1036 EDIT_dset_items(new_dset, eigs_brik+4, "eigvec_1[2]", ADN_none);
1037 EDIT_dset_items(new_dset, eigs_brik+5, "eigvec_1[3]", ADN_none);
1038 EDIT_dset_items(new_dset, eigs_brik+6, "eigvec_2[1]", ADN_none);
1039 EDIT_dset_items(new_dset, eigs_brik+7, "eigvec_2[2]", ADN_none);
1040 EDIT_dset_items(new_dset, eigs_brik+8, "eigvec_2[3]", ADN_none);
1041 EDIT_dset_items(new_dset, eigs_brik+9, "eigvec_3[1]", ADN_none);
1042 EDIT_dset_items(new_dset, eigs_brik+10,"eigvec_3[2]", ADN_none);
1043 EDIT_dset_items(new_dset, eigs_brik+11,"eigvec_3[3]", ADN_none);
1044 EDIT_dset_items(new_dset, eigs_brik+12,"FA", ADN_none);
1045 EDIT_dset_items(new_dset, eigs_brik+13,"MD", ADN_none);
1046 // pt,Oct,2016: for RD
1047 EDIT_dset_items(new_dset, eigs_brik+14,"RD", ADN_none);
1048 }
1049
1050 if(debug_briks) {
1051 EDIT_dset_items (new_dset, ADN_brick_label_one + nbriks - 4,
1052 "Converge Step", ADN_none);
1053 EDIT_dset_items (new_dset, ADN_brick_label_one + nbriks - 3,
1054 "ED", ADN_none);
1055 EDIT_dset_items (new_dset, ADN_brick_label_one + nbriks - 2,
1056 "EDorig", ADN_none);
1057 EDIT_dset_items (new_dset, ADN_brick_label_one + nbriks - 1,
1058 "I0", ADN_none);
1059 }
1060
1061 tross_Make_History ("3dDWItoDT", argc, argv, new_dset);
1062 if(sep_dsets)
1063 Save_Sep_DTdata(new_dset, prefix, datum);
1064 else {
1065 DSET_write (new_dset);
1066 INFO_message("--- Output dataset %s", DSET_BRIKNAME(new_dset));
1067 }
1068 }
1069 else
1070 {
1071 ERROR_exit("*** Error - Unable to compute output dataset!");
1072 }
1073
1074 INFO_message("Finished main tensor+parameter calcs.");
1075
1076 if ( debug_briks && (new_dset != NULL) ) {
1077 // pointer new_dset should just be to the AFNI-style DT, so this
1078 // should be fine...
1079
1080 INFO_message("Start creating stats dset");
1081
1082 char ffprefix[THD_MAX_PREFIX], fitprefix[THD_MAX_PREFIX];
1083 char *ext=NULL, nullch;
1084
1085 THD_3dim_dataset *CHIS=NULL;
1086 float **ccc=NULL;
1087
1088 // ------------------ allocate ----------------------------
1089 // make output arrays
1090 ccc = calloc( nxyz, sizeof(ccc)); // major output set
1091 for(i=0 ; i<2 ; i++)
1092 ccc[i] = calloc( nxyz, sizeof(float));
1093
1094 if( ccc == NULL ) {
1095 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
1096 exit(4);
1097 }
1098
1099 // -------------------- calc ---------------------------
1100
1101 EDIT_dset_items ( new_dsetDT,
1102 ADN_func_type, old_dset->func_type,
1103 ADN_type , old_dset->type,
1104 ADN_ntt , DSET_NVALS (new_dsetDT),
1105 ADN_ttorg , 0.0,
1106 ADN_ttdel , 1.0,
1107 ADN_ttdur , old_dset->taxis->ttdur,
1108 ADN_nsl , old_dset->taxis->nsl,
1109 ADN_zorg_sl , old_dset->taxis->zorg_sl,
1110 ADN_toff_sl , old_dset->taxis->toff_sl,
1111 ADN_tunits , UNITS_SEC_TYPE ,
1112 ADN_none);
1113
1114 // make the reference volume for the modeled signals be the I0
1115 // from the Powell_J value
1116 data_im = DSET_BRICK (old_dset0, 0);
1117 if( data_im == NULL ) {
1118 ERROR_exit("Null pointer when trying to make data_im set.");
1119 }
1120
1121 fac = DSET_BRICK_FACTOR(old_dset0, 0); // get scale factor for
1122 // each sub-brik
1123 if(fac==0.0) fac=1.0;
1124 if((data_im->kind != MRI_float)) {
1125 fprintf (stderr, "*** Error - Can only open float datasets. "
1126 "Use 3dcalc to convert.\n");
1127 mri_free (grad1Dptr);
1128 mri_free (data_im);
1129 exit (1);
1130 }
1131
1132 I0_ptr = mri_data_pointer(data_im) ; /* pointer to I0 data */
1133
1134 if( I0_ptr == NULL ) {
1135 ERROR_exit("Null pointer when trying to make I0 set.");
1136 }
1137
1138 sprintf(ffprefix, "%s", prefix);
1139 if(has_known_non_afni_extension(prefix)){ // for NIFTI, 3D, Niml,
1140 // Analyze,...
1141 ext = find_filename_extension(prefix);
1142 // remove non-afni-extension for now
1143 ffprefix[strlen(prefix) - strlen(ext)] = '\0';
1144 }
1145 else {
1146 nullch = '\0';
1147 ext = &nullch;
1148 }
1149 sprintf(fitprefix,"%s_CHI%s", ffprefix, ext);
1150
1151 if ( new_dset == NULL)
1152 ERROR_exit("New dset got reset to null?");
1153
1154 // calculates FIT dset: !!! prob not use fitprefix *here*
1155 fit_dset = MAKER_4D_to_typed_fbuc (new_dsetDT, /* input dataset */
1156 fitprefix, /* output prefix */
1157 datum, /* output datum */
1158 0, /* ignore count */
1159 0, /* can't detrend in maker
1160 function KRH 12/02 */
1161 nvols, /* number of briks */
1162 DTtoDWI_tsfunc, /* timeseries processor */
1163 NULL, /* data for tsfunc */
1164 NULL, /* mask */
1165 0 /* Allow auto scaling of output */
1166 );
1167
1168 if ( fit_dset == NULL)
1169 ERROR_message("'fit' dset was somehow null!");
1170
1171 // calc the goodness of fit!
1172 i = ChiSq_GOF( old_dset,
1173 fit_dset,
1174 B0list,
1175 ccc,
1176 maskptr);
1177
1178 CHIS = EDIT_empty_copy( fit_dset );
1179
1180 EDIT_dset_items( CHIS,
1181 ADN_datum_all , MRI_float,
1182 ADN_ntt , 2,
1183 ADN_nvals , 2,
1184 ADN_prefix , fitprefix,
1185 ADN_none );
1186
1187 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(CHIS)) )
1188 ERROR_exit("Can't overwrite existing dataset '%s'",
1189 DSET_HEADNAME(CHIS));
1190
1191 for( n=0; n<2 ; n++) {
1192 EDIT_substitute_brick(CHIS, n, MRI_float, ccc[n]);
1193 ccc[n]=NULL;
1194 }
1195
1196 EDIT_BRICK_LABEL(CHIS, 0, "chi_p");
1197 EDIT_BRICK_LABEL(CHIS, 1, "chi_c");
1198
1199 THD_load_statistics(CHIS);
1200 tross_Make_History("3dDWItoDT", argc, argv, CHIS);
1201 THD_write_3dim_dataset(NULL, NULL, CHIS, True);
1202 INFO_message("--- Output dataset %s", DSET_BRIKNAME(CHIS));
1203
1204 // Count number of non-b0s
1205 n = 0;
1206 for(i=0 ; i<nvols ; i++)
1207 if( !B0list[i] )
1208 n++;
1209 INFO_message("There were %d non-reference dsets (b>%.4f), "
1210 "and %d reference ones.",
1211 n, BMAX_REF, nvols-n);
1212
1213 if(ccc) {
1214 for(i=0 ; i<2 ; i++)
1215 free(ccc[i]);
1216 free(ccc);
1217 }
1218 if (CHIS) {
1219 DSET_delete(CHIS);
1220 free(CHIS);
1221 }
1222
1223 /*
1224 tross_Copy_History (old_dset, fit_dset);
1225 tross_Make_History ("3dDWItoDT", argc, argv, fit_dset);
1226 //for(i=0 ; i<nvols ; i++)
1227 // EDIT_dset_items (fit_dset, ADN_brick_label_one + i, "aa", ADN_none);
1228
1229 DSET_write (fit_dset);
1230 INFO_message("--- Output dataset %s", DSET_BRIKNAME(fit_dset));
1231 */
1232 }
1233 INFO_message("Finished with all calcs: cleanup.");
1234
1235
1236
1237 // ---------------- older cleanup ----------------------
1238
1239 if (maskptr)
1240 {
1241 free (maskptr);
1242 if(anat_im)
1243 mri_free (anat_im);
1244 #if 0
1245 DSET_unload_one (old_dset, 0);
1246 sar = NULL;
1247 #endif
1248 }
1249
1250
1251 // !!!!! ADD MORE !!!!
1252 if (fit_dset) {
1253 DSET_delete(fit_dset);
1254 free(fit_dset);
1255 }
1256 if (old_dset0) {
1257 DSET_delete(old_dset0);
1258 free(old_dset0);
1259 }
1260 if (new_dsetDT) {
1261 DSET_delete(new_dsetDT);
1262 free(new_dsetDT);
1263 }
1264
1265 FreeGlobals ();
1266 mri_free (grad1Dptr);
1267 matrix_destroy (&Rtmat);
1268
1269 exit (0);
1270 }
1271
1272 /*! save separate datasets for each kind of output */
1273 static void
Save_Sep_DTdata(whole_dset,prefix,output_datum)1274 Save_Sep_DTdata(whole_dset, prefix, output_datum)
1275 THD_3dim_dataset *whole_dset; /* whole dataset */
1276 char *prefix;
1277 int output_datum;
1278 {
1279 /* takes base prefix and appends to it for DT, eigvalues,
1280 eigvectors, FA, MD, debug bricks */
1281 int nbriks;
1282 char nprefix[THD_MAX_PREFIX], tprefix[THD_MAX_PREFIX];
1283 char *ext, nullch;
1284
1285 ENTRY("Save_Sep_DTdata");
1286 sprintf(tprefix,"%s",prefix);
1287 if(has_known_non_afni_extension(prefix)){ /* for NIFTI, 3D, Niml,
1288 Analyze,...*/
1289 ext = find_filename_extension(prefix);
1290 tprefix[strlen(prefix) - strlen(ext)] = '\0'; /* remove
1291 non-afni-extension
1292 for now*/
1293 }
1294 else {
1295 nullch = '\0';
1296 ext = &nullch;
1297 }
1298
1299 sprintf(nprefix,"%s_DT%s", tprefix,ext);
1300 Copy_dset_array(whole_dset,0,6, nprefix, output_datum);
1301 if(eigs_flag) {
1302 sprintf(nprefix,"%s_L1%s", tprefix,ext);
1303 Copy_dset_array(whole_dset,6,1, nprefix, output_datum);
1304 sprintf(nprefix,"%s_L2%s", tprefix,ext);
1305 Copy_dset_array(whole_dset,7,1, nprefix, output_datum);
1306 sprintf(nprefix,"%s_L3%s", tprefix,ext);
1307 Copy_dset_array(whole_dset,8,1, nprefix, output_datum);
1308 sprintf(nprefix,"%s_V1%s", tprefix,ext);
1309 Copy_dset_array(whole_dset,9,3, nprefix, output_datum);
1310 sprintf(nprefix,"%s_V2%s", tprefix,ext);
1311 Copy_dset_array(whole_dset,12,3, nprefix, output_datum);
1312 sprintf(nprefix,"%s_V3%s", tprefix,ext);
1313 Copy_dset_array(whole_dset,15,3, nprefix, output_datum);
1314 sprintf(nprefix,"%s_FA%s", tprefix,ext);
1315 Copy_dset_array(whole_dset,18,1, nprefix, output_datum);
1316 sprintf(nprefix,"%s_MD%s", tprefix,ext);
1317 Copy_dset_array(whole_dset,19,1, nprefix, output_datum);
1318 sprintf(nprefix,"%s_RD%s", tprefix,ext); // pt,Oct,2016: for RD
1319 Copy_dset_array(whole_dset,20,1, nprefix, output_datum);
1320 }
1321 if(debug_briks) {
1322 sprintf(nprefix,"%s_debugbriks%s", tprefix,ext);
1323 nbriks = whole_dset->dblk->nvals;
1324 Copy_dset_array(whole_dset,nbriks-4,4, nprefix, output_datum);
1325 }
1326
1327 EXRETURN;
1328 }
1329
1330 /*! create new dataset from part of existing dataset in memory */
1331 static void
Copy_dset_array(whole_dset,startbrick,nbriks,prefix,output_datum)1332 Copy_dset_array(whole_dset,startbrick,nbriks,prefix,output_datum)
1333 THD_3dim_dataset *whole_dset;
1334 int startbrick, nbriks;
1335 char *prefix;
1336 int output_datum;
1337 {
1338 THD_3dim_dataset *out_dset;
1339
1340 int i, ierror;
1341 MRI_IMAGE *fim;
1342 void *dataptr;
1343 float *fbuf;
1344
1345 ENTRY("Copy_dset_array");
1346
1347 out_dset = EDIT_empty_copy(whole_dset) ;
1348 fbuf = (float *) malloc (sizeof(float) * nbriks);
1349
1350 tross_Copy_History (whole_dset, out_dset);
1351 ierror = EDIT_dset_items( out_dset ,
1352 ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
1353 ADN_prefix , prefix ,
1354 ADN_datum_all, output_datum,
1355 ADN_nvals, nbriks,
1356 ADN_ntt, 0,
1357 ADN_type , ISHEAD(whole_dset) /* dataset type */
1358 ? HEAD_FUNC_TYPE
1359 : GEN_FUNC_TYPE ,
1360 ADN_func_type , FUNC_BUCK_TYPE , /* function type */
1361 ADN_none ) ;
1362
1363 if(ierror>0)
1364 ERROR_exit("*** Error - Unable to edit dataset!");
1365
1366 THD_init_datablock_keywords( out_dset->dblk ) ;
1367 THD_init_datablock_stataux( out_dset->dblk ) ; /* for some reason,
1368 need to do this
1369 for single brick
1370 NIFTI files */
1371
1372 /* attach brick, factors and labels to new dataset using existing
1373 brick pointers */
1374 for(i=0;i<nbriks;i++) {
1375 fim = DSET_BRICK(whole_dset,startbrick+i);
1376 dataptr = mri_data_pointer(fim);
1377 fbuf[i] = whole_dset->dblk->brick_fac[startbrick+i];
1378 /* copy labels here too.....*/
1379 EDIT_dset_items (out_dset, ADN_brick_label_one + i,
1380 whole_dset->dblk->brick_lab[startbrick+i], ADN_none);
1381 /*----- attach mri_image pointer to to be sub-brick #i -----*/
1382 EDIT_substitute_brick(out_dset, i, output_datum, dataptr);
1383 }
1384
1385 (void) EDIT_dset_items( out_dset , ADN_brick_fac , fbuf , ADN_none ) ;
1386 DSET_write (out_dset);
1387 INFO_message("--- Output dataset %s", DSET_BRIKNAME(out_dset));
1388 /*----- deallocate memory -----*/
1389 THD_delete_3dim_dataset (out_dset, False); out_dset = NULL ;
1390 free (fbuf); fbuf = NULL;
1391 EXRETURN;
1392 }
1393
1394
1395 /* Form R matrix as matrix of [bxx 2bxy 2bxz byy 2byz bzz] for Np rows */
1396 static void
Form_R_Matrix(MRI_IMAGE * grad1Dptr)1397 Form_R_Matrix (MRI_IMAGE * grad1Dptr)
1398 {
1399 matrix Rmat;
1400 register double sf = 1.0; /* scale factor = 1.0 for now until we
1401 know DELTA, delta (gamma = 267.5
1402 rad/ms-mT) */
1403 register double sf2; /* just scale factor * 2 */
1404 int i, nrows, noff; // @@new: offset counter, 'noff'= number for offset
1405 register float *imptr, *Gxptr, *Gyptr, *Gzptr;
1406 register float *Bxxptr, *Byyptr, *Bzzptr, *Bxyptr, *Bxzptr, *Byzptr;
1407 matrix *nullptr = NULL;
1408 register double Gx, Gy, Gz, Bxx, Byy, Bzz, Bxy, Bxz, Byz;
1409 double gscale;
1410
1411 ENTRY ("Form_R_Matrix");
1412 sf2 = sf + sf; /* 2 * scale factor for minor speed improvement */
1413
1414 if (grad1Dptr->ny == 6) {
1415 /* just read from b-matrix and scale as necessary */
1416 noff = grad1Dptr->nx; // offset is always whole dim of array
1417 if(BMAT_NZ) // for non-zero first row of BMATs
1418 nrows = grad1Dptr->nx;
1419 else
1420 nrows = grad1Dptr->nx - 1; //$$old: same for 'bmatrix_given' option
1421 matrix_initialize (&Rmat);
1422 matrix_create (nrows, 6, &Rmat); /* Rmat = Np x 6 matrix */
1423 if (Rmat.elts == NULL)
1424 { /* memory allocation error */
1425 ERROR_message("Could not allocate memory for Rmat");
1426 EXRETURN;
1427 }
1428
1429 /* use simple floating point pointers to get values */
1430 Bxxptr = imptr = MRI_FLOAT_PTR (grad1Dptr);
1431 if( !BMAT_NZ ) {// @@new switch: start at the nonzero grad
1432 Bxxptr+= 1;
1433 imptr+= 1;
1434 }
1435
1436 Byyptr = imptr + noff; // @@new: needed to iterate by dim of array
1437 Bzzptr = Byyptr + noff;
1438 Bxyptr = Bzzptr + noff;
1439 Bxzptr = Bxyptr + noff;
1440 Byzptr = Bxzptr + noff;
1441
1442 for (i = 0; i < nrows; i++)
1443 {
1444 Bxx = *Bxxptr++;
1445 Byy = *Byyptr++;
1446 Bzz = *Bzzptr++;
1447 Bxy = *Bxyptr++;
1448 Bxz = *Bxzptr++;
1449 Byz = *Byzptr++;
1450 Rmat.elts[i][0] = sf * Bxx; /* bxx = Gx*Gx*scalefactor */
1451 Rmat.elts[i][1] = sf2 * Bxy; /* 2bxy = 2GxGy*scalefactor */
1452 Rmat.elts[i][2] = sf2 * Bxz; /* 2bxz = 2GxGz*scalefactor */
1453 Rmat.elts[i][3] = sf * Byy; /* byy = Gy*Gy*scalefactor */
1454 Rmat.elts[i][4] = sf2 * Byz; /* 2byz = 2GyGz*scalefactor */
1455 Rmat.elts[i][5] = sf * Bzz; /* bzz = Gz*Gz*scalefactor */
1456 }
1457 }
1458 else {
1459 nrows = grad1Dptr->nx;
1460 matrix_initialize (&Rmat);
1461 matrix_create (nrows, 6, &Rmat); /* Rmat = Np x 6 matrix */
1462 if (Rmat.elts == NULL)
1463 { /* memory allocation error */
1464 ERROR_message("Could not allocate memory for Rmat");
1465 EXRETURN;
1466 }
1467
1468 /* use simple floating point pointers to get values */
1469 Gxptr = imptr = MRI_FLOAT_PTR (grad1Dptr);
1470 Gyptr = imptr + nrows;
1471 Gzptr = Gyptr + nrows;
1472
1473 for (i = 0; i < nrows; i++)
1474 {
1475 Gx = *Gxptr++;
1476 Gy = *Gyptr++;
1477 Gz = *Gzptr++;
1478 gscale = sqrt(Gx*Gx + Gy*Gy + Gz*Gz); // apr,2016: scale by this
1479 if ( gscale < TINYNUMBER) // for b=0
1480 gscale = 1.;
1481 /* bxx = Gx*Gx*scalefactor */
1482 Rmat.elts[i][0] = sf * Gx * Gx / gscale;
1483 /* 2bxy = 2GxGy*scalefactor */
1484 Rmat.elts[i][1] = sf2 * Gx * Gy / gscale;
1485 /* 2bxz = 2GxGz*scalefactor */
1486 Rmat.elts[i][2] = sf2 * Gx * Gz / gscale;
1487 /* byy = Gy*Gy*scalefactor */
1488 Rmat.elts[i][3] = sf * Gy * Gy / gscale;
1489 /* 2byz = 2GyGz*scalefactor */
1490 Rmat.elts[i][4] = sf2 * Gy * Gz / gscale;
1491 /* bzz = Gz*Gz*scalefactor */
1492 Rmat.elts[i][5] = sf * Gz * Gz / gscale;
1493 }
1494 }
1495
1496 /*
1497 fprintf(stderr,"\n\n TESTING: Bmatrix values...\n");
1498 for (i = 0; i < nrows; i++)
1499 INFO_message("%f %f %f %f %f %f",
1500 Rmat.elts[i][0],
1501 Rmat.elts[i][1],
1502 Rmat.elts[i][2],
1503 Rmat.elts[i][3],
1504 Rmat.elts[i][4],
1505 Rmat.elts[i][5]);
1506 */
1507
1508 matrix_initialize (&Rtmat);
1509 /* compute pseudo-inverse of Rmat=Rtmat */
1510 matrix_psinv (Rmat, nullptr, &Rtmat);
1511 matrix_destroy (&Rmat); /* from the other two matrices */
1512 EXRETURN;
1513 }
1514
1515
1516
1517
1518
1519 /*
1520 [PT: May, 2017]
1521 !! NB !! This has been copy/pasted from 3dDTtoDWI.c, so if changes
1522 happen there/here, likely make them here/there as well.
1523 This inglorious route was chosen rather than deal with
1524 compile-time issues with static variables. I'm sure I'll
1525 adjust this more properly later. Definitely.
1526 */
1527 static void
DTtoDWI_tsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)1528 DTtoDWI_tsfunc (double tzero, double tdelta,
1529 int npts, float ts[],
1530 double ts_mean, double ts_slope,
1531 void *ud, int nbriks, float *val)
1532 {
1533 int i;
1534 static int nvox2, ncall2;
1535 double I0, bq_d;
1536 double *bptr;
1537 float *tempptr;
1538
1539 ENTRY ("DTtoDWI_tsfunc");
1540 /* ts is input vector data of 6 floating point numbers.*/
1541 /* ts should come from data sub-briks in form of Dxx, Dxy, Dyy, Dxz, Dyz, Dzz */
1542 /* if automask is turned on, ts has 7 floating point numbers */
1543 /* val is output vector of form DWI0, DWI1, DWIn sub-briks */
1544 /* where n = number of gradients = nbriks-1 */
1545
1546 /** is this a "notification"? **/
1547 if (val == NULL) {
1548 if (npts > 0) { /* the "start notification" */
1549 nvox2 = npts; /* keep track of */
1550 ncall2 = 0; /* number of calls */
1551 }
1552 else { /* the "end notification" */
1553
1554 /* nothing to do here */
1555 }
1556 EXRETURN;
1557 }
1558
1559 ncall2++;
1560
1561 if (automask)
1562 npts = npts - 1;
1563
1564 tempptr = I0_ptr+ncall2-1;
1565 I0 = *tempptr;
1566
1567 val[0] = I0; /* the first sub-brik is the I0 sub-brik */
1568 bptr = bmatrix+6; /* start at the first gradient */
1569
1570 // NOTE THE ORDER HERE! Need to match bmatrix and DT orderings,
1571 // which are different. Yay.
1572 for(i=1;i<nbriks;i++) {
1573 bptr = bmatrix+(6*i); /* start at the first gradient */
1574 bq_d = *bptr++ * ts[0]; /* GxGxDxx */
1575 bq_d += *bptr++ * ts[1] * 2; /* 2GxGyDxy */
1576 bq_d += *bptr++ * ts[3] * 2; /* 2GxGzDxz */
1577 bq_d += *bptr++ * ts[2]; /* GyGyDyy */
1578 bq_d += *bptr++ * ts[4] * 2; /* 2GyGzDyz */
1579 bq_d += *bptr++ * ts[5]; /* GzGzDzz */
1580
1581 // for each gradient,q, Iq = J e -(bq.D)
1582 val[i] = I0 * exp(-bq_d);
1583 // where bq.D is the large dot product of bq and D
1584
1585 }
1586
1587 EXRETURN;
1588 }
1589
1590
1591 /**********************************************************************
1592 Function that does the real work
1593 ***********************************************************************/
1594
1595 static void
DWItoDT_tsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)1596 DWItoDT_tsfunc (double tzero, double tdelta,
1597 int npts, float ts[],
1598 double ts_mean, double ts_slope,
1599 void *ud, int nbriks, float *val)
1600 {
1601 int i, converge_step, converge, trialstep, ntrial, adjuststep, recordflag;
1602 double orig_deltatau, best_deltatau, EDold, J;
1603 static int ncall, noisecall;
1604 register double i0;
1605 register double dv, dv0;
1606 vector lnvector;
1607 int wtflag; /* wtflag for recomputing wtfactor*/
1608 int max_converge_step, graphpoint;
1609 double dtau[50], Edgraph[50];
1610 int graphflag;
1611 int IS_BAD_MD = 0;
1612
1613 ENTRY ("DWItoDT_tsfunc");
1614 /* ts is input vector data of Np+1 floating point numbers.
1615 For each point in volume brik convert vector data to
1616 symmetric matrix */
1617 /* ts should come from data sub-briks in form of I0,I1,...Ip */
1618 /* val is output vector of form
1619 Dxx Dxy Dxz Dyy Dyz Dzz
1620 for each voxel in 6 sub-briks */
1621 /* the Dij vector is computed as the product of Rt times ln(I0/Ip)
1622 where Rt is the pseudo-inverse of the [bxx 2bxy 2bxz byy 2byz
1623 bzz] for each gradient vector b */
1624 /** is this a "notification"? **/
1625 if (val == NULL){
1626
1627 if (npts > 0)
1628 { /* the "start notification" */
1629 /*nvox = npts;*/ /* keep track of */
1630 ncall = 0; /* number of calls */
1631 noisecall = 0;
1632 }
1633 else
1634 { /* the "end notification" */
1635
1636 /* nothing to do here */
1637 }
1638 EXRETURN;
1639 }
1640
1641 ncall++;
1642 /* if there is any mask (automask or user mask), use corresponding
1643 voxel as a flag */
1644 if((maskptr && (maskptr[ncall-1]==0)) ||
1645 all_dwi_zero(npts, ts)){
1646 /* don't include this voxel for mask or if all zeros*/
1647 for (i = 0; i < nbriks; i++) /* faster to copy preset vector */
1648 val[i] = 0.0; /* return 0 for all Dxx,Dxy,... */
1649 if(debug_briks) /* use -3 as flag for number of converge steps
1650 to mean exited for masked voxels */
1651 val[nbriks-4] = -3.0;
1652 EXRETURN;
1653 }
1654
1655 /* check for valid data at this series */
1656 if(bad_DWI_data(npts, ts)) {
1657 for (i = 0; i < nbriks; i++) /* faster to copy preset vector */
1658 val[i] = 0.0; /* return 0 for all Dxx,Dxy,... */
1659 Assign_CSF_values(val);
1660 EXRETURN;
1661 }
1662
1663
1664 /* load the symmetric matrix vector from the "timeseries" subbrik
1665 vector values */
1666 vector_initialize (&lnvector);
1667 vector_create_noinit (npts - 1, &lnvector);
1668 if(use_mean_b0)
1669 dv0 = compute_mean_B0(npts, ts, 1);
1670 else
1671 dv0 = ts[0];
1672
1673
1674 if (dv0 > 0.0)
1675 i0 = log (dv0);
1676 else
1677 i0 = 0.0;
1678 for (i = 0; i < (npts - 1); i++)
1679 {
1680 dv = ts[i + 1];
1681 if ((dv > 0.0) && (dv0 > 0.0))
1682 lnvector.elts[i] = i0 - log (dv); /* ln I0/Ip = ln I0 - ln Ip */
1683 else
1684 lnvector.elts[i] = 0.0;
1685 }
1686
1687 vector_multiply (Rtmat, lnvector, &Dvector); /* D = Rt * ln(I0/Ip),
1688 allocated Dvector
1689 here */
1690
1691 vector_destroy (&lnvector); /* free vector elements allocated */
1692
1693 if((method==0)||(max_iter==-1)) { /* for linear method,stop here
1694 and return D values */
1695 vector_to_array(Dvector, val);
1696
1697 if(debug_briks) {
1698 InitWtfactors (npts); /* initialize all weight factors to 1
1699 for all gradient intensities */
1700 J = ComputeJ (ts, npts); /* Ed (error) computed here */
1701 val[nbriks-4] = -1.0; /* use -1 as flag for number of
1702 converge steps to mean exited
1703 for */
1704 /* or initial insignificant deltatau
1705 value */
1706 val[nbriks-3] = ED;
1707 val[nbriks-2] = ED; /* store original error */
1708 val[nbriks-1] = J;
1709 }
1710
1711 if(eigs_flag) { /* if user wants eigenvalues in
1712 output dataset */
1713 EIG_func(); /* calculate eigenvalues,
1714 eigenvectors here */
1715 for(i=0;i<12;i++)
1716 val[i+6] = eigs[i];
1717 /* calc FA */
1718 val[18] = Calc_FA(val+6); /* calculate fractional anisotropy */
1719 val[19] = Calc_MD(val+6); /* calculate mean diffusivity */
1720 }
1721
1722 vals_to_NIFTI(val);
1723
1724 EXRETURN;
1725 }
1726
1727
1728 /* now more complex part that takes into account noise */
1729
1730 /* calculate initial estimate of D using standard linear model */
1731 EIG_func (); /* compute eigenvalues, eigenvectors standard
1732 way */
1733
1734
1735 InitWtfactors (npts); /* initialize all weight factors to 1 for all
1736 gradient intensities */
1737
1738 ComputeD0 (); /* recalculate Dmatrix based on limits on
1739 eigenvalues */
1740
1741
1742 if(matrix_sumabs(Dmatrix)<=TINYNUMBER) {
1743 for(i=0;i<nbriks;i++)
1744 val[i] = 0.0;
1745 if(debug_briks) {
1746 val[nbriks-4] = -2.0; /* use -2 as flag for number of converge
1747 steps to mean exited for insignificant
1748 D values*/
1749 val[nbriks-3] = 0;
1750 val[nbriks-2] = 0; /* store original error */
1751 val[nbriks-1] = 0;
1752 }
1753 vals_to_NIFTI(val); /* swap D tensor values for NIFTI
1754 standard */
1755 EXRETURN;
1756 }
1757
1758 /* show verbose messages every verbose=n voxels */
1759 if(verbose&&(!(noisecall%verbose)))
1760 recordflag = 1;
1761 else
1762 recordflag = 0;
1763 if(ncall==202293)
1764 recordflag = 1;
1765
1766 /* graph in AFNI convergence steps every afnitalk_flag=n voxels */
1767 if(afnitalk_flag&&(!(noisecall%afnitalk_flag))) {
1768 graphflag = 1;
1769 graphpoint = 0;
1770 }
1771 else
1772 graphflag = 0;
1773
1774 noisecall++;
1775
1776 /* allow up to max_iter=MAX_CONVERGE_STEPS (10) deltatau steps */
1777 converge_step = 0;
1778 /* 1st time through set limit of converge steps to user option */
1779 max_converge_step = max_iter;
1780
1781 /* need to use Powell optimize method instead */
1782 if( (opt_method==1) || ((opt_method==2) && (voxel_opt_method==1))) {
1783 if(recordflag)
1784 printf("using powell method\n");
1785 converge_step = ComputeDwithPowell(ts, val, npts, nbriks); /*compute D */
1786 Dmatrix.elts[0][0] = val[0];
1787 Dmatrix.elts[0][1] = val[1];
1788 Dmatrix.elts[0][2] = val[2];
1789 Dmatrix.elts[1][1] = val[3];
1790 Dmatrix.elts[1][2] = val[4];
1791 Dmatrix.elts[2][2] = val[5];
1792
1793 goto Other_Bricks; /* compute the other bricks for eigenvalues
1794 and debugging */
1795 }
1796
1797 converge = 0;
1798 wtflag = reweight_flag;
1799
1800 /* trial step */
1801 J = ComputeJ (ts, npts); /* Ed (error) computed here */
1802 Store_Computations (0, npts, wtflag); /* Store 1st adjusted
1803 computations */
1804 matrix_copy (Dmatrix, &OldD); /* store first Dmatrix in OldD
1805 too */
1806
1807 if(debug_briks)
1808 val[nbriks-2] = ED; /* store original error */
1809
1810 EDold = ED;
1811 ComputeDeltaTau ();
1812 if(deltatau<=TINYNUMBER) { /* deltatau too small, exit */
1813 for(i=0;i<nbriks;i++)
1814 val[i] = 0.0;
1815 if(debug_briks) {
1816 val[nbriks-4] = -1.0; /* use -1 as flag for number
1817 of converge steps to mean
1818 exited for insignificant
1819 deltatau value */
1820 val[nbriks-1] = J;
1821 }
1822 vals_to_NIFTI(val); /* swap D tensor values for
1823 NIFTI standard */
1824 EXRETURN;
1825 }
1826
1827 ntrial = 0;
1828
1829 while ( (converge_step < max_converge_step) &&
1830 (converge!=1) && (ntrial < 10) )
1831 {
1832 /* find trial step */
1833 /* first do trial step to find acceptable delta tau */
1834 /* first trial step is same size as previous time step */
1835 /* Then take half of previous tau step to see if less error */
1836 /* if there is less error, then delta tau is acceptable */
1837 /* Stop at first time step giving less error */
1838 /* try halving up to 10 times, if it does not work, */
1839 /* use first D from initial estimate (previous iteration) */
1840 trialstep = 1;
1841 ntrial = 0;
1842 orig_deltatau = best_deltatau = deltatau;
1843 while ((trialstep) && (ntrial < 10))
1844 { /* allow up to 10 iterations to find trial step */
1845 ComputeHpHm (deltatau);
1846 ComputeNewD ();
1847 J = ComputeJ (ts, npts); /* Ed (error) computed here */
1848 if (ED < EDold) /* is error less than error
1849 of trial step or previous
1850 step? */
1851 {
1852 /* found acceptable step size of DeltaTau */
1853 if(recordflag==1)
1854 INFO_message("ncall= %d, converge_step=%d, deltatau=%f, ED=%f, ntrial %d in find dtau",
1855 ncall, converge_step, deltatau, ED, ntrial);
1856 if(graphflag==1) {
1857 dtau[graphpoint] = deltatau;
1858 Edgraph[graphpoint] = ED;
1859 graphpoint++;
1860 }
1861
1862 EDold = ED;
1863 best_deltatau = deltatau;
1864 trialstep = 0; /* leave trial stepping */
1865 Store_Computations (1, npts, wtflag); /* Store
1866 current
1867 computations */
1868 }
1869 else {
1870 Restore_Computations (0, npts, wtflag); /* Restore
1871 trial 0
1872 computations */
1873 deltatau = deltatau / 2; /* find first DeltaTau
1874 step with less error
1875 than 1st guess */
1876 /* by trying smaller step sizes */
1877 ntrial++;
1878 }
1879 }
1880
1881 deltatau = best_deltatau;
1882
1883 /* end of finding trial step size */
1884 /* in trial step stage, already have result of deltatau step
1885 and may have already tried deltatau*2 if halved
1886 (ntrial>=1) */
1887 if(ntrial <10) {
1888 orig_deltatau = best_deltatau = deltatau;
1889 adjuststep = 1;
1890
1891 for(i=0;i<2;i++) {
1892 if(i==0)
1893 deltatau = orig_deltatau*2.0;
1894 else
1895 deltatau = orig_deltatau/2.0;
1896
1897 if((adjuststep==1) && ((i!=0) || (ntrial<2))) {
1898 /* if didn't shrink in initial deltatau step above */
1899
1900 Restore_Computations (1, npts, wtflag); /* Restore
1901 previous
1902 Tau
1903 step
1904 computations */
1905 ComputeHpHm (deltatau);
1906 ComputeNewD ();
1907 J = ComputeJ (ts, npts); /* computes Intensity
1908 without noise,*/
1909 /* Ed, Residuals */
1910 if(ED<EDold){
1911 best_deltatau = deltatau;
1912 adjuststep = 0;
1913 Store_Computations(0, npts, wtflag); /* Store
1914 Tau+dtau
1915 step
1916 computations */
1917 EDold = ED;
1918
1919 if(recordflag==1) {
1920 if(i==0)
1921 INFO_message("ncall= %d, converge_step=%d, deltatau=%f, ED=%f dt*2 best",
1922 ncall, converge_step, deltatau, ED);
1923 else
1924 INFO_message("ncall= %d, converge_step=%d, deltatau=%f, ED=%f dt/2 best",
1925 ncall, converge_step, deltatau, ED);
1926 }
1927 if(graphflag==1) {
1928 dtau[graphpoint] = deltatau;
1929 Edgraph[graphpoint] = ED;
1930 graphpoint++;
1931 }
1932 }
1933 }
1934 }
1935
1936 deltatau = best_deltatau;
1937
1938 if(adjuststep!=0){ /* best choice was first
1939 Delta Tau */
1940 ED = EDold;
1941 Restore_Computations (1, npts, wtflag); /* restore
1942 old
1943 computed
1944 matrices*/
1945 /* D,Hp,Hm,F,R */
1946 Store_Computations(0, npts, wtflag); /* Store
1947 Tau+dtau step
1948 computations */
1949 if(recordflag==1)
1950 INFO_message("ncall= %d, converge_step=%d, deltatau=%f, ED=%f dt best",
1951 ncall, converge_step, deltatau, ED);
1952 }
1953
1954 if(graphflag==1) {
1955 dtau[graphpoint] = deltatau;
1956 Edgraph[graphpoint] = ED;
1957 graphpoint++;
1958 }
1959
1960 if (converge_step != 0) { /* first time through
1961 recalculate*/
1962 /* now see if converged yet */
1963 converge = TestConvergence(Dmatrix, OldD);
1964 }
1965
1966 matrix_copy (Dmatrix, &OldD);
1967
1968 if(graphflag==1) {
1969 dtau[graphpoint] = deltatau;
1970 Edgraph[graphpoint] = ED;
1971 graphpoint++;
1972 }
1973
1974 converge_step++;
1975 }
1976 else
1977 {
1978 if(recordflag==1)
1979 INFO_message("ncall= %d, converge_step=%d, deltatau=%f, ED=%f Exiting time evolution",
1980 ncall, converge_step, deltatau, ED);
1981 Restore_Computations(0, npts, wtflag); /* Exit
1982 with
1983 original
1984 step
1985 calculation */
1986 ED = EDold;
1987 }
1988
1989
1990 if(((converge) || (converge_step==max_iter)) && wtflag && reweight_flag) {
1991 /* if at end of iterations the first time*/
1992
1993 converge = 0; /* through whole group of iterations */
1994 ComputeWtfactors (npts); /* compute new weight factors */
1995 converge_step = 1; /* start over now with computed weight factors */
1996 max_converge_step = max_iter_rw+1; /* reset limit of converge steps to user option */
1997 wtflag = 0; /* only do it once - turn off next reweighting */
1998 J=ComputeJ(ts, npts); /* compute new Ed value */
1999 EDold = ED; /* this avoids having to go through converge loop for two loops */
2000 }
2001 } /* end while converge loop */
2002
2003 ED = EDold;
2004
2005 val[0] = Dmatrix.elts[0][0];
2006 val[1] = Dmatrix.elts[0][1];
2007 val[2] = Dmatrix.elts[0][2];
2008 val[3] = Dmatrix.elts[1][1];
2009 val[4] = Dmatrix.elts[1][2];
2010 val[5] = Dmatrix.elts[2][2];
2011
2012
2013
2014
2015
2016 Other_Bricks:
2017
2018 // May,2016: final tensor has been estimated-- check on badness 2!
2019 if(bad_DWI_data_2_MD(val[0], val[3], val[5])) {
2020 for (i = 0; i < nbriks; i++) /* faster to copy preset vector */
2021 val[i] = 0.0; /* return 0 for all Dxx,Dxy,... */
2022 Assign_CSF_values(val);
2023 EXRETURN;
2024 }
2025
2026 if(eigs_flag) { /* if user wants eigenvalues in output dataset */
2027 double length, dl;
2028 udmatrix_to_vector(Dmatrix, &Dvector);
2029 EIG_func(); /* calculate eigenvalues, eigenvectors here */
2030 for(i=0;i<3;i++) {
2031 if(fabs(eigs[i])<SMALLNUMBER)
2032 eigs[i] = 0.0;
2033 }
2034 for(i=0;i<12;i++)
2035 val[i+6] = eigs[i];
2036 /* calc FA */
2037 val[18] = Calc_FA(val+6); /* calculate fractional anisotropy */
2038 val[19] = Calc_MD(val+6); /* calculate average (mean) diffusivity */
2039 val[20] = Calc_RD(val+6); // calculate radial diffusivity; pt,Oct,2016
2040 length = sqrt(eigs[3]*eigs[3]+eigs[4]*eigs[4]+eigs[5]*eigs[5]);
2041 dl = fabs(1.0 - length);
2042 if (dl>SMALLNUMBER) {
2043 if(recordflag) {
2044 printf("V1 not a unit vector, length = %f\n", length);
2045 printf(" V1 eigs %f %f %f\n", eigs[3], eigs[4], eigs[5]);
2046 printf("D tensor %f %f %f %f %f %f\n",
2047 val[0],val[1],val[2],val[3],val[4],val[5]);
2048 printf("Time series:\n");
2049 for(i=0;i<npts;i++) printf("%f ", ts[i]);
2050 printf("\n");
2051 }
2052 }
2053 }
2054
2055 /* testing information only */
2056 if(recordflag)
2057 INFO_message("ncall= %d, converge_step=%d, deltatau=%f, ED=%f",
2058 ncall, converge_step, deltatau, ED);
2059 if(debug_briks && ((opt_method==0) || ((opt_method==2) && (voxel_opt_method==0)) )){
2060 val[nbriks-4] = converge_step;
2061 val[nbriks-3] = ED;
2062 val[nbriks-1] = ComputeJ(ts, npts); /* compute J value */;
2063 }
2064 if(graphflag==1) {
2065 DWI_AFNI_update_graph(Edgraph, dtau, graphpoint);
2066 }
2067
2068 vals_to_NIFTI(val); /* swap D tensor values for NIFTI standard */
2069
2070 EXRETURN;
2071 }
2072
2073 /* some simple checks for bad DWI data */
2074 static int
bad_DWI_data(int npts,float * ts)2075 bad_DWI_data(int npts, float *ts)
2076 {
2077 double m0, m1;
2078
2079 ENTRY("bad_DWI_data");
2080 /* check if mean of B0 is less than mean of rest of data */
2081 m0 = compute_mean_B0(npts, ts, 1);
2082 m1 = compute_mean_B0(npts, ts, 0);
2083 if(m0<=m1) RETURN(1); /* bad data */
2084 RETURN(0);
2085 }
2086
2087 // May,2016-- added; extra badness check
2088 static int
bad_DWI_data_2_MD(float v1,float v2,float v3)2089 bad_DWI_data_2_MD(float v1, float v2, float v3)
2090 {
2091 double m0;
2092 float v[3] = {0., 0., 0.};
2093
2094 ENTRY("bad_DWI_data_2_MD");
2095
2096 v[0] = v1;
2097 v[1] = v2;
2098 v[2] = v3;
2099
2100 m0 = Calc_MD(v);
2101 if(m0 > MIN_BAD_SCALE_MD*csf_val) {
2102 //WARNING_message("Whoa! MD is reeeeally big: %f", m0);
2103 RETURN(1); /* bad data */
2104 }
2105
2106 RETURN(0);
2107 }
2108
2109 static int
all_dwi_zero(int npts,float * ts)2110 all_dwi_zero(int npts, float *ts)
2111 {
2112 int i;
2113
2114 ENTRY("all_DWI_zero");
2115 for(i=0;i<npts;i++) {
2116 if(ts[i]!=0.0) RETURN(0);
2117 }
2118 RETURN(1);
2119 }
2120
2121 static double
compute_mean_B0(int npts,float * ts,char B0flag)2122 compute_mean_B0(int npts, float *ts, char B0flag)
2123 {
2124 int i, nb=0;
2125 double m0=0.0, sum=0.0;
2126
2127 ENTRY("compute_mean_B0");
2128 for(i=0;i<npts;i++) {
2129 if(B0list[i] == B0flag) {
2130 sum += (double) ts[i];
2131 nb++;
2132 }
2133 }
2134 m0 = sum / nb;
2135 RETURN(m0);
2136 }
2137
2138 static void
Assign_CSF_values(float * val)2139 Assign_CSF_values(float *val)
2140 {
2141 /* assign default CSF values to tensor */
2142 // apr,2016: lots of '1.0' values changed to 'csf_val' for
2143 // consistency, and to not cause visualization issues anywhere.
2144 // Hopefully. Maybe. It still doesn't match technically with the
2145 // csf_fa, but at least that is *almost* zero.
2146 val[0] = val[3] = val[5] = csf_val;
2147 val[1] = val[2] = val[4] = 0.0;
2148 if(eigs_flag) { /* if user wants eigenvalues in
2149 output dataset */
2150 val[6] = csf_val;
2151 val[7] = csf_val;
2152 val[8] = csf_val;
2153 val[9] = 1.0; val[10] = 0.0; val[11] = 0.0; // eigenvectors
2154 val[12] = 0.0; val[13] = 1.0; val[14] = 0.0;
2155 val[15] = 0.0; val[16] = 0.0; val[17] = 1.0;
2156 val[18] = csf_fa; /* calculate fractional
2157 anisotropy */
2158 val[19] = csf_val; /* calculate average (mean)
2159 diffusivity */
2160 val[20] = rd_val; /* calculate average (mean)
2161 diffusivity */
2162 }
2163 }
2164
2165 /* taken from 3dDTeig.c */
2166 /*! given Dij tensor data for principal directions */
2167 /* calculate 3 principal sets of eigenvalues and eigenvectors */
2168 static void
EIG_func()2169 EIG_func ()
2170 {
2171 /* THD_dmat33 inmat;
2172 THD_dvecmat eigvmat; */
2173 int i, j;
2174 int maxindex, minindex, midindex;
2175 float temp, minvalue, maxvalue;
2176 int sortvector[3];
2177 double a[9], e[3];
2178 int astart, vstart;
2179
2180
2181 ENTRY ("EIG_func");
2182 /* Dvector is vector data of 6 floating point numbers.
2183 For each point in volume brik convert vector data to
2184 symmetric matrix */
2185 /* Dvector should come in form of Dxx,Dxy,Dxz,Dyy,Dyz,Dzz */
2186 /* convert to matrix of form
2187 [ Dxx Dxy Dxz]
2188 [ Dxy Dyy Dyz]
2189 [ Dxz Dyz Dzz] */
2190
2191
2192 /* load the symmetric matrix vector from the "timeseries" subbrik
2193 vector values */
2194
2195 a[0] = Dvector.elts[0];
2196 a[1] = Dvector.elts[1];
2197 a[2] = Dvector.elts[2];
2198 a[3] = Dvector.elts[1];
2199 a[4] = Dvector.elts[3];
2200 a[5] = Dvector.elts[4];
2201 a[6] = Dvector.elts[2];
2202 a[7] = Dvector.elts[4];
2203 a[8] = Dvector.elts[5];
2204
2205 symeig_double (3, a, e); /* compute eigenvalues in e,
2206 eigenvectors in a */
2207
2208 maxindex = 2; /* find the lowest, middle and highest
2209 eigenvalue */
2210 maxvalue = e[2];
2211 minindex = 0;
2212 minvalue = e[0];
2213 midindex = 1;
2214
2215 for (i = 0; i < 3; i++)
2216 {
2217 temp = e[i];
2218 if (temp > maxvalue)
2219 { /* find the maximum */
2220 maxindex = i;
2221 maxvalue = temp;
2222 }
2223 if (temp < minvalue)
2224 { /* find the minimum */
2225 minindex = i;
2226 minvalue = temp;
2227 }
2228 }
2229
2230 for (i = 0; i < 3; i++)
2231 { /* find the middle */
2232 if ((i != maxindex) && (i != minindex))
2233 {
2234 midindex = i;
2235 break;
2236 }
2237 }
2238
2239 sortvector[0] = maxindex;
2240 sortvector[1] = midindex;
2241 sortvector[2] = minindex;
2242
2243 /* put the eigenvalues at the beginning of the matrix */
2244 for (i = 0; i < 3; i++)
2245 {
2246 eigs[i] = e[sortvector[i]]; /* copy sorted eigenvalues */
2247 if(fabs (eigs[i]) < TINYNUMBER)
2248 eigs[i] = 0.0;
2249 /* start filling in eigenvector values */
2250 astart = sortvector[i] * 3; /* start index of double
2251 eigenvector */
2252 vstart = (i + 1) * 3; /* start index of float val
2253 vector to copy eigenvector */
2254
2255 for (j = 0; j < 3; j++)
2256 {
2257 eigs[vstart + j] = a[astart + j];
2258 }
2259 }
2260
2261 EXRETURN;
2262 }
2263
2264 /*! calculate Fractional Anisotropy */
2265 /* passed float pointer to start of eigenvalues */
2266 static float
Calc_FA(float * val)2267 Calc_FA(float *val)
2268 {
2269 float FA;
2270 double ssq, dv0, dv1, dv2, dsq;
2271
2272 ENTRY("Calc_FA");
2273
2274 /* calculate the Fractional Anisotropy, FA */
2275 /* reference, Pierpaoli C, Basser PJ. Microstructural and
2276 physiological features of tissues elucidated by
2277 quantitative-diffusion tensor MRI,J Magn Reson B 1996;
2278 111:209-19 */
2279 if((val[0]<=0.0)||(val[1]<0.0)||(val[2]<0.0)) {
2280
2281 /* any negative eigenvalues-- should not see any for non-linear
2282 method. Set FA to 0*/
2283 RETURN(0.0);
2284 }
2285
2286 /* sum of squares of eigenvalues */
2287 ssq = (val[0]*val[0])+(val[1]*val[1])+(val[2]*val[2]);
2288 /* dsq = pow((val[0]-val[1]),2.0) + pow((val[1]-val[2]),2.0) +
2289 pow((val[2]-val[0]),2.0);*/ /* sum of differences squared */
2290
2291 dv0 = val[0]-val[1];
2292 dv0 *= dv0;
2293 dv1 = val[1]-val[2];
2294 dv1 *= dv1;
2295 dv2 = val[2]-val[0];
2296 dv2 *= dv2;
2297 dsq = dv0+dv1+dv2; /* sum of differences squared */
2298
2299 if(ssq!=0.0)
2300 FA = sqrt(dsq/(2.0*ssq)); /* FA calculated here */
2301 else
2302 FA = 0.0;
2303
2304 RETURN(FA);
2305 }
2306
2307 /*! calculate Mean Diffusivity */
2308 /* passed float pointer to start of eigenvalues */
2309 static float
Calc_MD(float * val)2310 Calc_MD(float *val)
2311 {
2312 float MD;
2313
2314 ENTRY("Calc_MD");
2315
2316 /* calculate the Fractional Anisotropy, FA */
2317 /* reference, Pierpaoli C, Basser PJ. Microstructural and
2318 physiological features of tissues elucidated by
2319 quantitative-diffusion tensor MRI,J Magn Reson B 1996;
2320 111:209-19 */
2321 if((val[0]<=0.0)||(val[1]<0.0)||(val[2]<0.0)) {
2322 /* any negative eigenvalues-- should not see any for non-linear
2323 method. Set FA to 0*/
2324 RETURN(0.0);
2325 }
2326 MD = (val[0] + val[1] + val[2]) / 3;
2327
2328 RETURN(MD);
2329 }
2330
2331 //! calculate radial Diffusivity; pt,Oct,2016
2332 /* passed float pointer to start of eigenvalues */
2333 static float
Calc_RD(float * val)2334 Calc_RD(float *val)
2335 {
2336 float RD;
2337
2338 ENTRY("Calc_RD");
2339
2340 /* calculate the Fractional Anisotropy, RD */
2341 if((val[0]<=0.0)||(val[1]<0.0)||(val[2]<0.0)) {
2342 /* any negative eigenvalues-- should not see any for non-linear
2343 method. Set FA to 0*/
2344 RETURN(0.0);
2345 }
2346 RD = (val[1] + val[2]) / 2.;
2347
2348 RETURN(RD);
2349 }
2350
2351
2352
2353
2354
2355
2356 /*! compute initial estimate of D0 */
2357 /* D = estimated diffusion tensor matrix
2358 [Dxx Dxy Dxz, Dxy Dyy Dyz, Dxz Dyz Dzz] */
2359 /* updates Dvector and Dmatrix */
2360 static void
ComputeD0()2361 ComputeD0 ()
2362 {
2363 int i, j;
2364 /* matrix ULmatrix, Ematrix; */
2365 double mu, alpha, sum;
2366 double e10, e11, e12, e20, e21, e22, e30, e31, e32;
2367 double l1, l2, l3;
2368 double t1, t3, t5, t8, t10, t12, t14, t18, t19, t21, t23, t32, t33, t35,
2369 t37;
2370
2371 ENTRY ("ComputeD0");
2372 /* create and initialize D0 */
2373
2374 if (eigs[0] < 0)
2375 { /* if all eigenvalues are negative - may never happen */
2376 /* D0 = diag(a,a,a) where a=1/3 Sum(Abs(Lambda_i)) */
2377 sum = 0.0;
2378 for (i = 0; i < 3; i++)
2379 sum += fabs (eigs[i]);
2380 alpha = sum / 3;
2381 for (i = 0; i < 3; i++)
2382 {
2383 for (j = 0; j < 3; j++)
2384 {
2385 if (i == j)
2386 Dmatrix.elts[i][j] = alpha;
2387 else
2388 Dmatrix.elts[i][j] = 0.0;
2389 }
2390 }
2391
2392 /* convert to vector format for D also */
2393 udmatrix_to_vector (Dmatrix, &Dvector);
2394 EXRETURN;
2395 }
2396
2397 mu = backoff_factor * eigs[0];
2398 voxel_opt_method = 0;
2399 /* mu = SMALLNUMBER;*/
2400 if (eigs[1] < mu) { /* set limit of eigenvalues to prevent */
2401 eigs[1] = mu; /* too much anisotropy */
2402 voxel_opt_method = 1; /* switch to Powell optimization for
2403 this voxel */
2404 }
2405 if (eigs[2] < mu) {
2406 eigs[2] = mu;
2407 voxel_opt_method = 1;
2408 }
2409 /* D0 = U L UT */
2410 /*
2411 [e10 l1 e20 l2 e30 l3]
2412 [ ]
2413 UL := [e11 l1 e21 l2 e31 l3]
2414 [ ]
2415 [e12 l1 e22 l2 e32 l3]
2416 */
2417
2418
2419 /* assign variables to match Maple code */
2420 l1 = eigs[0];
2421 l2 = eigs[1];
2422 l3 = eigs[2];
2423
2424 e10 = eigs[3];
2425 e11 = eigs[4];
2426 e12 = eigs[5];
2427
2428 e20 = eigs[6];
2429 e21 = eigs[7];
2430 e22 = eigs[8];
2431
2432 e30 = eigs[9];
2433 e31 = eigs[10];
2434 e32 = eigs[11];
2435
2436
2437 #ifdef lkjsaklfj
2438 matrix_initialize (&Ematrix);
2439 matrix_create (3, 3, &Ematrix);
2440 /* fill Ematrix with Eigenvectors */
2441
2442 matrix_initialize (&ULmatrix);
2443 matrix_create (3, 3, &ULmatrix);
2444 if (ULmatrix.elts == NULL)
2445 { /* memory allocation error */
2446 ERROR_message("Could not allocate memory for Rmat");
2447 EXRETURN;
2448 }
2449
2450 for (i = 0; i < 3; i++)
2451 {
2452 for (j = 0; j < 3; j++)
2453 {
2454 ULmatrix.elts[i][j] = Ematrix.elts[i][j] * eigs[i];
2455 }
2456 }
2457
2458 /* new D is based on modified lambdas */
2459 matrix_multiply (ULmatrix, Ematrix, &Dmatrix);
2460 #endif
2461
2462
2463
2464 t1 = e10 * e10;
2465 t3 = e20 * e20;
2466 t5 = e30 * e30;
2467 t8 = e10 * l1;
2468 t10 = e20 * l2;
2469 t12 = e30 * l3;
2470 t14 = t8 * e11 + t10 * e21 + t12 * e31;
2471 t18 = t8 * e12 + t10 * e22 + t12 * e32;
2472 t19 = e11 * e11;
2473 t21 = e21 * e21;
2474 t23 = e31 * e31;
2475 t32 = e11 * l1 * e12 + e21 * l2 * e22 + e31 * l3 * e32;
2476 t33 = e12 * e12;
2477 t35 = e22 * e22;
2478 t37 = e32 * e32;
2479 Dmatrix.elts[0][0] = t1 * l1 + t3 * l2 + t5 * l3;
2480 Dmatrix.elts[0][1] = t14;
2481 Dmatrix.elts[0][2] = t18;
2482 Dmatrix.elts[1][0] = t14;
2483 Dmatrix.elts[1][1] = t19 * l1 + t21 * l2 + t23 * l3;
2484 Dmatrix.elts[1][2] = t32;
2485 Dmatrix.elts[2][0] = t18;
2486 Dmatrix.elts[2][1] = t32;
2487 Dmatrix.elts[2][2] = t33 * l1 + t35 * l2 + t37 * l3;
2488
2489 /* convert to vector format for D */
2490 udmatrix_to_vector (Dmatrix, &Dvector);
2491
2492 EXRETURN;
2493 }
2494
2495 /*! compute the diffusion weighting matrix bmatrix for q number of gradients */
2496 /* only need to calculate once */
2497 /* bq = diffusion weighting matrix of qth gradient */
2498 /* GxGx GxGy GxGz
2499 GxGy GyGy GyGz
2500 GxGz GyGz GzGz */
2501 /* b0 is 0 for all 9 elements */
2502 /* bmatrix is really stored as 6 x npts array */
2503 // -----> !! if you change anything here, probably should copy/paste
2504 // -----> !! the function into 3dDTtoDWI, for consistency across
2505 // -----> !! functions (as of apr,2016)!
2506 static void
Computebmatrix(MRI_IMAGE * grad1Dptr,int NO_ZERO_ROW1)2507 Computebmatrix (MRI_IMAGE * grad1Dptr, int NO_ZERO_ROW1)
2508 // flag to differentiate bmatrix cases that include a zero row or not
2509 {
2510 int i, n;
2511 register double *bptr;
2512 register float *Gxptr, *Gyptr, *Gzptr;
2513 register float *Bxxptr, *Byyptr, *Bzzptr, *Bxyptr, *Bxzptr, *Byzptr;
2514 double Gx, Gy, Gz, Bxx, Byy, Bzz, Bxy, Bxz, Byz;
2515 double gscale;
2516
2517 ENTRY ("Computebmatrix");
2518 n = grad1Dptr->nx; /* number of gradients other than I0 */
2519
2520 if ( (grad1Dptr->ny == 6) && !NO_ZERO_ROW1 ) { // extra switch to
2521 // keep OLD,
2522 // zero-row version
2523 /* just read in b matrix */
2524 Bxxptr = MRI_FLOAT_PTR (grad1Dptr); /* use simple floating point
2525 pointers to get values */
2526 Byyptr = Bxxptr + n;
2527 Bzzptr = Byyptr + n;
2528 Bxyptr = Bzzptr + n;
2529 Bxzptr = Bxyptr + n;
2530 Byzptr = Bxzptr + n;
2531
2532 bptr = bmatrix;
2533
2534 /* B0list[0]= 1;*/ /* keep a record of which volumes have no
2535 gradient: first one always assumed */
2536
2537 for (i = 0; i < n; i++){
2538 Bxx = *Bxxptr++;
2539 Byy = *Byyptr++;
2540 Bzz = *Bzzptr++;
2541 Bxy = *Bxyptr++;
2542 Bxz = *Bxzptr++;
2543 Byz = *Byzptr++;
2544 *bptr++ = Bxx;
2545 *bptr++ = Bxy;
2546 *bptr++ = Bxz;
2547 *bptr++ = Byy;
2548 *bptr++ = Byz;
2549 *bptr++ = Bzz;
2550
2551 // if(Bxx==0.0 && Byy==0.0 && Bzz==0.0)
2552
2553 // is this a zero gradient volume also? -> user can input a
2554 // larger value, if necessary.
2555 if( (Bxx+Byy+Bzz)<BMAX_REF )
2556 B0list[i] = 1;
2557 else{
2558 B0list[i] = 0;
2559 // apr,2016: need the MAX_BVAL for scaling
2560 gscale = Bxx + Byy + Bzz;
2561 if(gscale > MAX_BVAL)
2562 MAX_BVAL = gscale;
2563 }
2564
2565 }
2566 }
2567 else if( (grad1Dptr->ny == 6) && NO_ZERO_ROW1 ) { // very similar to old bmatrix option
2568 /* just read in b matrix */
2569 Bxxptr = MRI_FLOAT_PTR (grad1Dptr); /* use simple floating point pointers to get values */
2570 Byyptr = Bxxptr + n;
2571 Bzzptr = Byyptr + n;
2572 Bxyptr = Bzzptr + n;
2573 Bxzptr = Bxyptr + n;
2574 Byzptr = Bxzptr + n;
2575
2576 bptr = bmatrix;
2577
2578 // do as grads below
2579 for (i = 0; i < 6; i++)
2580 *bptr++ = 0.0; /* initialize first 6 elements to 0.0 for the I0 gradient */
2581 B0list[0]= 1; /* keep a record of which volumes have no gradient */
2582
2583
2584 for (i = 0; i < n; i++){
2585 Bxx = *Bxxptr++;
2586 Byy = *Byyptr++;
2587 Bzz = *Bzzptr++;
2588 Bxy = *Bxyptr++;
2589 Bxz = *Bxzptr++;
2590 Byz = *Byzptr++;
2591 *bptr++ = Bxx;
2592 *bptr++ = Bxy;
2593 *bptr++ = Bxz;
2594 *bptr++ = Byy;
2595 *bptr++ = Byz;
2596 *bptr++ = Bzz;
2597
2598 // is this a zero gradient volume also?
2599 //if(Bxx==0.0 && Byy==0.0 && Bzz==0.0)
2600 if( (Bxx+Byy+Bzz)<BMAX_REF )
2601 B0list[i+1] = 1;
2602 else{
2603 B0list[i+1] = 0;
2604 // apr,2016: need the MAX_BVAL for scaling
2605 gscale = Bxx + Byy + Bzz;
2606 if(gscale > MAX_BVAL)
2607 MAX_BVAL = gscale;
2608 }
2609
2610 }
2611 }
2612 else {
2613 Gxptr = MRI_FLOAT_PTR (grad1Dptr); /* use simple floating point pointers to get values */
2614 Gyptr = Gxptr + n;
2615 Gzptr = Gyptr + n;
2616
2617 bptr = bmatrix;
2618 for (i = 0; i < 6; i++)
2619 *bptr++ = 0.0; /* initialize first 6 elements to 0.0 for the I0 gradient */
2620
2621 B0list[0]= 1; /* keep a record of which volumes have no gradient */
2622
2623 for (i = 0; i < n; i++)
2624 {
2625 gscale = 1.; // apr,2016: allow for non-unit gradient magnitudes
2626 Gx = *Gxptr++;
2627 Gy = *Gyptr++;
2628 Gz = *Gzptr++;
2629
2630 //if((Gx==0.0) && (Gy==0.0) && (Gz==0.0))
2631 gscale = sqrt(Gx*Gx + Gy*Gy + Gz*Gz);
2632 if( gscale<BMAX_REF )
2633 B0list[i+1] = 1; /* no gradient applied*/
2634 else{
2635 B0list[i+1] = 0;
2636 if(gscale > MAX_BVAL)
2637 MAX_BVAL = gscale; // apr,2016
2638 }
2639
2640 *bptr++ = Gx * Gx / gscale; // apr,2016: allow for non-unit gradient magnitudes
2641 *bptr++ = Gx * Gy / gscale;
2642 *bptr++ = Gx * Gz / gscale;
2643 *bptr++ = Gy * Gy / gscale;
2644 *bptr++ = Gy * Gz / gscale;
2645 *bptr++ = Gz * Gz / gscale;
2646 }
2647 }
2648 EXRETURN;
2649 }
2650
2651
2652 /*! compute non-gradient intensity, J, based on current calculated values of
2653 diffusion tensor matrix, D */
2654 static double
ComputeJ(float ts[],int npts)2655 ComputeJ (float ts[], int npts)
2656 {
2657 /* J = Sum(wq Iq exp(-bq D)) / Sum (wq exp(-2bq D)) */
2658 /* estimate of b0 intensity without noise and applied gradient */
2659 /* Iq = voxel value for qth gradient */
2660 /* bq = diffusion weighting matrix of qth gradient */
2661 /* wq = weighting factor for qth gradient at Iq voxel */
2662 /* D = estimated diffusion tensor matrix
2663 [Dxx Dxy Dxz, Dxy Dyy Dyz, Dxz Dyz Dzz] */
2664 /* ts = Iq is time series voxel data from original data of intensities */
2665
2666 register int i, j;
2667 double sum0, sum1, b1D1, b2D2, b4D4, wtexpbD, J, tempcalc, sumbD, Fscalar;
2668 double *expbD, *expbDptr, *wtfactorptr, *Ftempmatrix;
2669 register double *Fptr, *Rptr, *bptr;
2670 double D0,D1,D2,D3,D4,D5;
2671
2672 ENTRY ("ComputeJ");
2673 sum0 = sum1 = 0.0;
2674 expbD = malloc (npts * sizeof (double)); /* allocate calculations for speed */
2675 expbDptr = expbD; /* temporary pointers for indexing */
2676 wtfactorptr = wtfactor;
2677 bptr = bmatrix; /* npts of b vectors (nx6) */
2678
2679 D0 = Dmatrix.elts[0][0];
2680 D1 = Dmatrix.elts[0][1];
2681 D2 = Dmatrix.elts[0][2];
2682 D3 = Dmatrix.elts[1][1];
2683 D4 = Dmatrix.elts[1][2];
2684 D5 = Dmatrix.elts[2][2];
2685
2686
2687 for (i = 0; i < npts; i++)
2688 {
2689 /* compute bq.D */
2690 /* bq.D is large dot product of b and D at qth gradient */
2691 /* large dot product for Hilbert algebra */
2692 /* regular dot product is for Hilbert space (vectors only)- who knew? */
2693 /* calculate explicitly rather than loop to save time */
2694 b1D1 = *(bptr + 1) * D1;
2695 b1D1 += b1D1;
2696 b2D2 = *(bptr + 2) * D2;
2697 b2D2 += b2D2;
2698 b4D4 = *(bptr + 4) * D4;
2699 b4D4 += b4D4;
2700
2701 sumbD = *bptr * D0 + b1D1 + b2D2 + /* bxxDxx + 2bxyDxy + 2bxzDxz + */
2702 (*(bptr + 3) * D3) + /* byyDyy + */
2703 b4D4 + /* 2byzDyz + */
2704 (*(bptr + 5) * D5); /* bzzDzz */
2705
2706 /* exp (-bq.D) */
2707 *expbDptr = exp (-sumbD);
2708 wtexpbD = *(wtfactor + i) * *expbDptr;
2709 sum0 += wtexpbD * ts[i];
2710 sum1 += wtexpbD * *expbDptr;
2711 expbDptr++;
2712 wtfactorptr++;
2713 bptr += 6; /* increment to next vector of bmatrix */
2714 }
2715
2716 J = sum0 / sum1;
2717 /* Now compute error functional,E(D,J) and gradient of E with respect to D ,Ed or F in notes */
2718 /* E(D,J)= 1/2 Sum[wq (J exp(-bq.D) - Iq)^2] */
2719 /* F = Ed = - Sum[wq (J exp(-bq.D) - Iq) bq] *//* Ed is a symmetric matrix */
2720 sum0 = 0.0;
2721 sigma = 0.0; /* standard deviation of noise for weight factors */
2722 expbDptr = expbD;
2723 wtfactorptr = wtfactor;
2724 /* initialize F matrix */
2725 Ftempmatrix = malloc(6*sizeof(double));
2726 Fptr = Ftempmatrix;
2727 for(i=0;i<6;i++)
2728 *Fptr++ = 0.0;
2729 Fptr = Ftempmatrix;
2730 Rptr = Rvector; /* residuals calculated here - used in Wt.factor calculations */
2731 bptr = bmatrix; /* npts of b vectors (nx6) */
2732 for (i = 0; i < npts; i++)
2733 {
2734 *Rptr = tempcalc = (J * *expbDptr) - ts[i];
2735 Fscalar = -*wtfactorptr * tempcalc;
2736 tempcalc = tempcalc * tempcalc;
2737
2738 for (j = 0; j < 6; j++)
2739 { /* for each entry of Fij (Fxx, Fxy,...) */
2740 /* F = - Sum[wq (J exp(-bq.D) - Iq) bq] = Sum[-wq (J exp(-bq.D) - Iq) bq] */
2741 *(Fptr+j) += Fscalar * (*bptr++); /* Fij = Fij + (Fscalar * bij) */
2742 }
2743
2744 sum0 += *wtfactorptr * tempcalc; /* E(D,J) = Sum (wq temp^2) */
2745 sigma += tempcalc; /* standard deviation of noise for weight factors */
2746 expbDptr++;
2747 wtfactorptr++;
2748 Rptr++;
2749 }
2750
2751 udmatrix_copy (Ftempmatrix, &Fmatrix); /* copy upper diagonal vector data into full matrix */
2752
2753 ED = sum0 / 2; /* this is the error for this iteration */
2754
2755 free (Ftempmatrix);
2756 free (expbD);
2757 RETURN (J);
2758 }
2759
2760 /*! compute initial step size for gradient descent */
2761 static void
ComputeDeltaTau()2762 ComputeDeltaTau ()
2763 {
2764 double sum0, sum1;
2765 matrix Dsqmatrix, FDsqmatrix, DsqFmatrix, Gmatrix;
2766 /* compute estimate of gradient, dD/dtau */
2767 /*G = [F] [D]^2 + [D]^2 [F] - ask Bob about ^2 and negative for this part to be sure */
2768
2769 ENTRY ("ComputeDeltaTau");
2770 matrix_initialize (&Dsqmatrix);
2771 matrix_initialize (&FDsqmatrix);
2772 matrix_initialize (&DsqFmatrix);
2773 matrix_initialize (&Gmatrix);
2774
2775 matrix_multiply (Dmatrix, Dmatrix, &Dsqmatrix); /* compute D^2 */
2776 matrix_multiply (Fmatrix, Dsqmatrix, &FDsqmatrix); /* FD^2 */
2777 matrix_multiply (Dsqmatrix, Fmatrix, &DsqFmatrix); /* D^2F */
2778 matrix_add (FDsqmatrix, DsqFmatrix, &Gmatrix); /* G= FD^2 +D^2F */
2779
2780
2781 /* deltatau = 0.01 * Sum(|Dij|) / Sum (|Gij|) */
2782 sum0 = matrix_sumabs (Dmatrix);
2783 sum1 = matrix_sumabs (Gmatrix);
2784 if (sum1 != 0.0)
2785 deltatau = 0.01 * sum0 / sum1;
2786 else
2787 deltatau = 0.0;
2788 matrix_destroy (&Dsqmatrix);
2789 matrix_destroy (&FDsqmatrix);
2790 matrix_destroy (&DsqFmatrix);
2791 matrix_destroy (&Gmatrix);
2792 EXRETURN;
2793 }
2794
2795 /*! allocate all the global matrices and arrays once */
2796 static void
InitGlobals(int npts)2797 InitGlobals (int npts)
2798 {
2799 int i;
2800 double *cumulativewtptr;
2801
2802 ENTRY ("InitGlobals");
2803 matrix_initialize (&Fmatrix);
2804 matrix_create (3, 3, &Fmatrix);
2805 matrix_initialize (&Dmatrix);
2806 matrix_create (3, 3, &Dmatrix);
2807 matrix_initialize (&Hplusmatrix);
2808 matrix_create (3, 3, &Hplusmatrix);
2809 matrix_initialize (&Hminusmatrix);
2810 matrix_create (3, 3, &Hminusmatrix);
2811 matrix_initialize (&OldD);
2812 matrix_create (3, 3, &OldD);
2813 for(i=0;i<2;i++){
2814 matrix_initialize (&tempFmatrix[i]);
2815 matrix_create (3, 3, &tempFmatrix[i]);
2816 matrix_initialize (&tempDmatrix[i]);
2817 matrix_create (3, 3, &tempDmatrix[i]);
2818 matrix_initialize (&tempHplusmatrix[i]);
2819 matrix_create (3, 3, &tempHplusmatrix[i]);
2820 matrix_initialize (&tempHminusmatrix[i]);
2821 matrix_create (3, 3, &tempHminusmatrix[i]);
2822 }
2823 Rvector = malloc (npts * sizeof (double));
2824 tempRvector = malloc (npts * sizeof(double));
2825 wtfactor = malloc (npts * sizeof (double));
2826 B0list = malloc (npts * sizeof (int));
2827
2828 if(cumulative_flag && reweight_flag) {
2829 cumulativewt = malloc (npts * sizeof (double));
2830 cumulativewtptr = cumulativewt;
2831 for(i=0;i<npts;i++)
2832 *cumulativewtptr++ = 0.0;
2833 rewtvoxels = 0;
2834 }
2835
2836 bmatrix = malloc (npts * 6 * sizeof (double));
2837
2838 vector_initialize (&Dvector); /* need to initialize vectors before 1st use-mod drg 12/20/2004 */
2839 /* vector_initialize (&tempDvector); vector_create(npts, &tempDvector);*/
2840 EXRETURN;
2841 }
2842
2843 /*! free up all the matrices and arrays */
2844 static void
FreeGlobals()2845 FreeGlobals ()
2846 {
2847 int i;
2848
2849 ENTRY ("FreeGlobals");
2850 matrix_destroy (&Fmatrix);
2851 matrix_destroy (&Dmatrix);
2852 matrix_destroy (&Hplusmatrix);
2853 matrix_destroy (&Hminusmatrix);
2854 matrix_destroy (&OldD);
2855 for(i=0;i<2;i++){
2856 matrix_destroy (&tempFmatrix[i]);
2857 matrix_destroy (&tempDmatrix[i]);
2858 matrix_destroy (&tempHplusmatrix[i]);
2859 matrix_destroy (&tempHminusmatrix[i]);
2860 }
2861
2862
2863 free (wtfactor);
2864 wtfactor = NULL;
2865 free (bmatrix);
2866 bmatrix = NULL;
2867 free (Rvector);
2868 Rvector = NULL;
2869 free (tempRvector);
2870 tempRvector = NULL;
2871 free(B0list);
2872 B0list = NULL;
2873
2874 vector_destroy (&Dvector); /* need to free elements of Dvector - mod-drg 12/20/2004 */
2875 /* vector_destroy (&tempDvector);*/
2876 if(cumulative_flag && reweight_flag) {
2877 free (cumulativewt);
2878 cumulativewt = NULL;
2879 }
2880 EXRETURN;
2881 }
2882
2883 /*! store current computed matrices D,Hp,Hm, R */
2884 static void
Store_Computations(int i,int npts,int wtflag)2885 Store_Computations (int i, int npts, int wtflag)
2886 {
2887 ENTRY ("Store_Computations");
2888
2889 matrix_copy (Fmatrix, &tempFmatrix[i]);
2890 matrix_copy (Dmatrix, &tempDmatrix[i]);
2891 matrix_copy (Hplusmatrix, &tempHplusmatrix[i]);
2892 matrix_copy (Hminusmatrix, &tempHminusmatrix[i]);
2893 if(wtflag==1)
2894 memcpy(tempRvector, Rvector, npts*sizeof(double));
2895 EXRETURN;
2896 }
2897
2898 /*! restore old computed matrices D,Hp,Hm, R */
2899 static void
Restore_Computations(int i,int npts,int wtflag)2900 Restore_Computations (int i, int npts, int wtflag)
2901 {
2902 ENTRY ("Restore_Computations");
2903
2904 matrix_copy (tempFmatrix[i], &Fmatrix);
2905 matrix_copy (tempDmatrix[i], &Dmatrix);
2906 matrix_copy (tempHplusmatrix[i], &Hplusmatrix);
2907 matrix_copy (tempHminusmatrix[i], &Hminusmatrix);
2908 if(wtflag==1)
2909 memcpy(Rvector, tempRvector, npts*sizeof(double));
2910 EXRETURN;
2911 }
2912
2913 /*! set all wt factors for all gradient levels to be 1.0 the first time through */
2914 static void
InitWtfactors(int npts)2915 InitWtfactors (int npts)
2916 {
2917 double *wtfactorptr;
2918 int i;
2919
2920 ENTRY ("InitWtfactors");
2921 wtfactorptr = wtfactor;
2922 for (i = 0; i < npts; i++)
2923 *wtfactorptr++ = 1.0;
2924 EXRETURN;
2925 }
2926
2927 /*! compute wt factors for each gradient level */
2928 static void
ComputeWtfactors(int npts)2929 ComputeWtfactors (int npts)
2930 {
2931 /* Residuals, rq, computed above in ComputeJ, stored in Rmatrix */
2932 /* unnormalized standard deviation, sigma, computed there too */
2933 /* wq = 1 / sqrt(1 + (rq/sigma)^2)
2934 where sigma = sqrt[1/Nq Sum(rq^2)] */
2935 /* and rq = J exp(-bq.D) - Iq */
2936
2937 int i;
2938 double *wtfactorptr, *Rptr;
2939 double *cumulativewtptr;
2940 double tempcalc, sum;
2941
2942 ENTRY ("ComputeWtfactors");
2943 sigma = sigma / npts;
2944 sigma = sqrt (sigma); /* sigma = std.dev. */
2945
2946 wtfactorptr = wtfactor;
2947 Rptr = Rvector;
2948
2949 sum = 0.0;
2950 for (i = 0; i < npts; i++)
2951 {
2952 tempcalc = *Rptr++ / sigma;
2953 tempcalc = tempcalc * tempcalc;
2954 tempcalc = 1.0 / (sqrt (1 + tempcalc));
2955 *wtfactorptr++ = tempcalc;
2956 sum += tempcalc;
2957 }
2958 /* now renormalize to avoid changing the relative value of E(D) */
2959 tempcalc = npts / sum; /* normalization factor */
2960
2961 wtfactorptr = wtfactor;
2962 for (i=0; i<npts; i++) {
2963 *wtfactorptr = *wtfactorptr * tempcalc;
2964 wtfactorptr++;
2965 }
2966
2967 if(cumulative_flag) {
2968 wtfactorptr = wtfactor;
2969 cumulativewtptr = cumulativewt;
2970 /* printf("Wt.factors: ");*/
2971 ++rewtvoxels;
2972 for (i=0; i<npts; i++){
2973 *cumulativewtptr++ += *wtfactorptr++; /* calculate cumulative wt.factor across all voxels*/
2974 }
2975 }
2976
2977 EXRETURN;
2978 }
2979
2980 /*! compute Hplus and Hminus as a function of delta tau */
2981 /* H+- = I +/- 1/2 deltatau F D */
2982 static void
ComputeHpHm(double deltatau)2983 ComputeHpHm (double deltatau)
2984 {
2985 matrix FDmatrix;
2986 double dtau;
2987 int i, j;
2988
2989 ENTRY ("ComputeHpHm");
2990 dtau = 0.5 * deltatau;
2991
2992 matrix_initialize (&FDmatrix);
2993 matrix_multiply (Fmatrix, Dmatrix, &FDmatrix);
2994 for (i = 0; i < 3; i++)
2995 for (j = 0; j < 3; j++)
2996 FDmatrix.elts[i][j] = dtau * FDmatrix.elts[i][j];
2997
2998 for (i = 0; i < 3; i++)
2999 {
3000 for (j = 0; j < 3; j++)
3001 {
3002 if (i == j)
3003 {
3004 Hplusmatrix.elts[i][j] = 1 + FDmatrix.elts[i][j]; /* I + dt/2 * FD */
3005 Hminusmatrix.elts[i][j] = 1 - FDmatrix.elts[i][j]; /* I - dt/2 * FD */
3006 }
3007 else
3008 {
3009 Hplusmatrix.elts[i][j] = Hminusmatrix.elts[i][j] =
3010 FDmatrix.elts[i][j];
3011 }
3012 }
3013 }
3014
3015 matrix_destroy (&FDmatrix);
3016 EXRETURN;
3017 }
3018
3019 /*! compute new D matrix */
3020 /* D(tau+deltatau) = H- H+^-1 D(tau) H+^-1 H- */
3021 /* = A D(tau) A^T */
3022 /* where A = H- H+^-1 */
3023 static void
ComputeNewD()3024 ComputeNewD ()
3025 {
3026 double *Hpinv;
3027 matrix Hpinvmatrix, Amatrix, ATmatrix, ADmatrix;
3028
3029 ENTRY ("ComputeNewD");
3030
3031 Hpinv =
3032 InvertSym3 (Hplusmatrix.elts[0][0], Hplusmatrix.elts[0][1],
3033 Hplusmatrix.elts[0][2], Hplusmatrix.elts[1][1],
3034 Hplusmatrix.elts[1][2], Hplusmatrix.elts[2][2]);
3035 matrix_initialize (&Hpinvmatrix);
3036 matrix_initialize (&Amatrix);
3037 matrix_initialize (&ATmatrix);
3038 matrix_initialize (&ADmatrix);
3039
3040 matrix_create (3, 3, &Hpinvmatrix);
3041 udmatrix_copy (Hpinv, &Hpinvmatrix); /* copy values from Hpinv vector to Hpinvmatrix */
3042
3043 matrix_multiply (Hminusmatrix, Hpinvmatrix, &Amatrix);
3044 matrix_multiply (Amatrix, Dmatrix, &ADmatrix);
3045 matrix_transpose (Amatrix, &ATmatrix);
3046 matrix_multiply (ADmatrix, ATmatrix, &Dmatrix);
3047
3048 matrix_destroy (&ADmatrix);
3049 matrix_destroy (&ATmatrix);
3050 matrix_destroy (&Amatrix);
3051 matrix_destroy (&Hpinvmatrix);
3052
3053 free (Hpinv);
3054 EXRETURN;
3055 }
3056
3057 /*! test convergence of calculation of D */
3058 /* if sum of differences hasn't changed by more than 1E-4 */
3059 /* then the calculations have converged */
3060 /* return 1 for convergence, 0 if not converged */
3061 static int
TestConvergence(matrix NewD,matrix OldD)3062 TestConvergence(matrix NewD, matrix OldD)
3063 {
3064 int converge;
3065 double convergence;
3066
3067 ENTRY ("TestConvergence");
3068 /* convergence test */
3069 convergence = fabs (NewD.elts[0][0] - OldD.elts[0][0]) + /* Dxx */
3070 fabs (NewD.elts[0][1] - OldD.elts[0][1]) + /* Dxy */
3071 fabs (NewD.elts[0][2] - OldD.elts[0][2]) + /* Dxz */
3072 fabs (NewD.elts[1][1] - OldD.elts[1][1]) + /* Dyy */
3073 fabs (NewD.elts[1][2] - OldD.elts[1][2]) + /* Dyz */
3074 fabs (NewD.elts[2][2] - OldD.elts[2][2]); /* Dzz */
3075
3076 if (convergence < SMALLNUMBER)
3077 converge = 1;
3078 else
3079 converge = 0;
3080
3081 RETURN (converge);
3082 }
3083
3084 /*! copy an upper diagonal matrix (6 point vector really) into a standard double
3085 array type matrix for n timepoints */
3086 /* ud0 ud1 ud2 m0 m1 m2
3087 ud3 ud4 --> m3 m4 m5
3088 ud5 m6 m7 m8 */
3089 static void
udmatrix_copy(double * udptr,matrix * m)3090 udmatrix_copy (double *udptr, matrix * m)
3091 {
3092 ENTRY ("udmatrix_copy");
3093
3094 m->elts[0][0] = *udptr;
3095 m->elts[0][1] = *(udptr + 1);
3096 m->elts[0][2] = *(udptr + 2);
3097 m->elts[1][0] = *(udptr + 1);
3098 m->elts[1][1] = *(udptr + 3);
3099 m->elts[1][2] = *(udptr + 4);
3100 m->elts[2][0] = *(udptr + 2);
3101 m->elts[2][1] = *(udptr + 4);
3102 m->elts[2][2] = *(udptr + 5);
3103 EXRETURN;
3104 }
3105
3106 /*! copy upper part of 3x3 matrix elements to 6-element vector elements */
3107 /* m1 m2 m3
3108 m4 m5 -> v = [m1 m2 m3 m4 m5 m6]
3109 m6
3110 */
3111 static void
udmatrix_to_vector(matrix m,vector * v)3112 udmatrix_to_vector (matrix m, vector * v)
3113 {
3114 ENTRY ("udmatrix_to_vector");
3115 v->elts[0] = m.elts[0][0];
3116 v->elts[1] = m.elts[0][1];
3117 v->elts[2] = m.elts[0][2];
3118 v->elts[3] = m.elts[1][1];
3119 v->elts[4] = m.elts[1][2];
3120 v->elts[5] = m.elts[2][2];
3121 EXRETURN;
3122 }
3123
3124 /*! sum the absolute value of all elements of a matrix */
3125 static double
matrix_sumabs(matrix m)3126 matrix_sumabs (matrix m)
3127 {
3128 register int i, j;
3129 register double sum;
3130
3131 ENTRY ("matrix_sumabs");
3132 sum = 0.0;
3133 for (i = 0; i < 3; i++)
3134 {
3135 for (j = 0; j < 3; j++)
3136 sum += fabs (m.elts[i][j]);
3137 }
3138 RETURN (sum);
3139 }
3140
3141 /*! calculate inverse of a symmetric 3x3 matrix */
3142 /* returns pointer to 9 element vector corresponding to 3x3 matrix */
3143 /* a b c */
3144 /* b e f */
3145 /* c f i */
3146 /* Maple generated code */
3147 static double *
InvertSym3(double a,double b,double c,double e,double f,double i)3148 InvertSym3 (double a, double b, double c, double e, double f, double i)
3149 {
3150 double *symmat, *symmatptr; /* invert matrix - actually 6 values in a vector form */
3151 double t2, t4, t7, t9, t12, t15, t20, t24, t30;
3152
3153 ENTRY ("InvertSym3");
3154 symmat = malloc (6 * sizeof (double));
3155 symmatptr = symmat;
3156 t2 = f * f;
3157 t4 = a * e;
3158 t7 = b * b;
3159 t9 = c * b;
3160 t12 = c * c;
3161 t15 = 1 / (t4 * i - a * t2 - t7 * i + 2.0 * t9 * f - t12 * e);
3162 t20 = (b * i - c * f) * t15;
3163 t24 = (b * f - c * e) * t15;
3164 t30 = (a * f - t9) * t15;
3165
3166 *symmatptr++ = (e * i - t2) * t15; /*B[0][0] */
3167 *symmatptr++ = -t20; /*B[0][1] */
3168 *symmatptr++ = t24; /*B[0][2] */
3169 /* B[1][0] = -t20; */
3170 *symmatptr++ = (a * i - t12) * t15; /* B[1][1] */
3171 *symmatptr++ = -t30; /* B [1][2] */
3172 /* B[2][0] = t24; */
3173 /* B[2][1] = -t30; */
3174 *symmatptr = (t4 - t7) * t15; /* B[2][2] */
3175
3176 RETURN (symmat);
3177 }
3178
3179 /*! copy elements from matrix a to matrix b */
3180 /* matrix_equate already exists but creates and initializes new matrix */
3181 /* steps we don't need to do here */
3182 /* assumes both a and b already exist with equal dimensions */
3183 static void
matrix_copy(matrix a,matrix * b)3184 matrix_copy (matrix a, matrix * b)
3185 {
3186 register int i;
3187 register int rows, cols;
3188
3189 ENTRY ("matrix_copy");
3190
3191 rows = a.rows;
3192 cols = a.cols;
3193
3194 for (i = 0; i < rows; i++)
3195 {
3196 #if 0
3197 register int j;
3198 for (j = 0; j < cols; j++)
3199 b->elts[i][j] = a.elts[i][j];
3200 #else
3201 if (cols > 0)
3202 memcpy (b->elts[i], a.elts[i], sizeof (double) * cols); /* RWCox */
3203 #endif
3204 }
3205 EXRETURN;
3206 }
3207
3208
3209 #define DWI_WriteCheckWaitMax 2000
3210 #define DWI_WriteCheckWait 400
3211 /*-----------------------------------------------------*/
3212 /* Stuff for an extra NIML port for non-SUMA programs. */
3213
3214 /* ZSS June 2011
3215 #ifndef NIML_TCP_FIRST_PORT
3216 #define NIML_TCP_FIRST_PORT 53212
3217 #endif
3218 Replace with:
3219 get_port_named("AFNI_DEFAULT_LISTEN_NIML");
3220 */
3221
3222 /*! open NIML stream */
DWI_Open_NIML_stream()3223 static int DWI_Open_NIML_stream()
3224 {
3225 int nn, Wait_tot, tempport;
3226 char streamname[256];
3227
3228 ENTRY("DWI_Open_NIML_stream");
3229
3230 /* contact afni */
3231 tempport = get_port_named("AFNI_DEFAULT_LISTEN_NIML");
3232 sprintf(streamname, "tcp:localhost:%d",tempport);
3233 INFO_message("Contacting AFNI");
3234
3235 DWIstreamid = NI_stream_open( streamname, "w" ) ;
3236 if (DWIstreamid==0) {
3237 WARNING_message("Warning - NI_stream_open failed");
3238 DWIstreamid = NULL;
3239 RETURN(1);
3240 }
3241
3242 INFO_message("Trying shared memory...");
3243 if (!NI_stream_reopen( DWIstreamid, "shm:DWIDT1M:1M" ))
3244 INFO_message("Warning: Shared memory communcation failed.");
3245 else
3246 INFO_message("Shared memory connection OK.");
3247 Wait_tot = 0;
3248 while(Wait_tot < DWI_WriteCheckWaitMax){
3249 nn = NI_stream_writecheck( DWIstreamid , DWI_WriteCheckWait) ;
3250 if( nn == 1 ){
3251 fprintf(stderr,"\n") ;
3252 RETURN(0) ;
3253 }
3254 if( nn < 0 ){
3255 WARNING_message("Bad connection to AFNI");
3256 DWIstreamid = NULL;
3257 RETURN(1);
3258 }
3259 Wait_tot += DWI_WriteCheckWait;
3260 fprintf(stderr,".") ;
3261 }
3262
3263 WARNING_message("WriteCheck timed out (> %d ms).",DWI_WriteCheckWaitMax);
3264 RETURN(1);
3265 }
3266
3267 /*! create the initial graph in AFNI - no points yet*/
DWI_NIML_create_graph()3268 static int DWI_NIML_create_graph()
3269 {
3270 NI_element *nel;
3271
3272 ENTRY("DWI_NIML_create_graph");
3273 nel = NI_new_data_element("ni_do", 0);
3274 NI_set_attribute ( nel, "ni_verb", "DRIVE_AFNI");
3275 NI_set_attribute ( nel, "ni_object",
3276 "OPEN_GRAPH_1D DWIConvEd 'DWI Convergence' "
3277 "25 1 'converge step' 1 0 300000 Ed");
3278 if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
3279 WARNING_message("Failed to send data to AFNI");
3280 NI_free_element(nel) ; nel = NULL;
3281 RETURN(1);
3282 }
3283 NI_free_element(nel) ;
3284 nel = NULL;
3285 RETURN(0);
3286 }
3287
3288 /*! create new graph with left and right y axes scaled from 0 to max1, max2*/
DWI_NIML_create_newgraph(npts,max1,max2)3289 static int DWI_NIML_create_newgraph(npts, max1, max2)
3290 int npts;
3291 double max1, max2;
3292 {
3293 NI_element *nel;
3294 char stmp[256];
3295 static int nx = -1;
3296
3297 ENTRY("DWI_NIML_create_newgraph");
3298 nel = NI_new_data_element("ni_do", 0);
3299 NI_set_attribute ( nel, "ni_verb", "DRIVE_AFNI");
3300 if((nx==-1) || (nx<npts)) /* 1st time through close any existing graph by that name*/
3301 NI_set_attribute ( nel, "ni_object","CLOSE_GRAPH_1D DWIConvEd\n"); /* have to close graph to change axes */
3302 else
3303 NI_set_attribute ( nel, "ni_object","CLEAR_GRAPH_1D DWIConvEd\n");
3304
3305 if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
3306 WARNING_message("Failed to send data to AFNI");
3307 NI_free_element(nel) ; nel = NULL;
3308 RETURN(1);
3309 }
3310
3311 if((nx==-1) || (nx<npts)) { /* update the graph only
3312 first time or if x-axis
3313 not big enough */
3314 nx = max_iter * 4 + 10;
3315 if(reweight_flag==1)
3316 nx += max_iter_rw * 4 + 10;
3317 if(nx<npts) /* fix graph to include
3318 largest number of
3319 steps */
3320 nx = npts;
3321 sprintf(stmp,"OPEN_GRAPH_1D DWIConvEd 'DWI Convergence' %d 1 'converge step' 2 0 100 %%Maximum Ed \\Delta\\tau\n",nx );
3322 NI_set_attribute ( nel, "ni_object", stmp);
3323 if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
3324 WARNING_message("Failed to send data to AFNI");
3325 NI_free_element(nel) ; nel = NULL;
3326 RETURN(1);
3327 }
3328 NI_set_attribute ( nel, "ni_object",
3329 "SET_GRAPH_GEOM DWIConvEd geom=700x400+100+400\n");
3330 if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
3331 WARNING_message("Failed to send data to AFNI");
3332 NI_free_element(nel) ; nel = NULL;
3333 RETURN(1);
3334 }
3335 }
3336 NI_free_element(nel) ;
3337 nel = NULL;
3338 RETURN(0);
3339 }
3340
3341 /*! tell AFNI to graph data for convergence steps */
DWI_AFNI_update_graph(double * Edgraph,double * dtau,int npts)3342 static void DWI_AFNI_update_graph(double *Edgraph, double *dtau, int npts)
3343 {
3344 NI_element *nel;
3345 char stmp[256];
3346 int i;
3347 double Edmax, dtaumax;
3348 double *Edptr, *dtauptr;
3349 double tempEd, temptau;
3350
3351 ENTRY("DWI_AFNI_update_graph");
3352
3353 Edmax = 0.0; dtaumax = 0.0;
3354 Edptr = Edgraph;
3355 dtauptr = dtau;
3356 for(i=0;i<npts;i++) {
3357 if(*Edptr>Edmax)
3358 Edmax = *Edptr;
3359 if(*dtauptr>dtaumax)
3360 dtaumax = *dtauptr;
3361 ++Edptr; ++dtauptr;
3362 }
3363
3364 NI_write_procins(DWIstreamid, "keep reading");
3365 DWI_NIML_create_newgraph(npts, Edmax, dtaumax);
3366 /* NI_sleep(250);*/
3367
3368 nel = NI_new_data_element("ni_do", 0);
3369 NI_set_attribute ( nel, "ni_verb", "DRIVE_AFNI");
3370 NI_set_attribute ( nel, "ni_object", "CLEAR_GRAPH_1D DWIConvEd\n");
3371 /* NI_sleep(25);*/
3372 if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
3373 WARNING_message("Failed to send data to AFNI");
3374 }
3375
3376
3377 for(i=0;i<npts;i++){
3378 if(Edmax!=0.0)
3379 tempEd = 100*Edgraph[i] / Edmax;
3380 else
3381 tempEd = 0.0;
3382 if(dtaumax!=0.0)
3383 temptau = 100*dtau[i] / dtaumax;
3384 else
3385 temptau = 0.0;
3386
3387 /* put rel.error, Ed, and deltatau for all the convergence steps */
3388 sprintf(stmp,"ADDTO_GRAPH_1D DWIConvEd %4.2f %4.2f\n", tempEd, temptau);
3389 NI_set_attribute ( nel, "ni_object", stmp); /* put command and data in stmp */
3390 NI_sleep(25); /* for dramatic effect */
3391 if (NI_write_element( DWIstreamid, nel, NI_BINARY_MODE ) < 0) {
3392 WARNING_message("Failed to send data to AFNI");
3393 }
3394 }
3395 NI_free_element(nel) ; nel = NULL;
3396 EXRETURN;
3397 }
3398
3399 /* need to put D tensor in NIFTI format standard */
vals_to_NIFTI(float * val)3400 static void vals_to_NIFTI(float *val)
3401 {
3402 float temp;
3403
3404 /* D tensor as lower triangular for NIFTI standard */
3405 temp = val[2];
3406 val[2] = val[3];
3407 val[3] = temp;
3408 }
3409
3410 /* function called at each optimization step */
DT_Powell_optimize_fun(int n,double * x)3411 double DT_Powell_optimize_fun(int n, double *x)
3412 {
3413 /* n should always be 6 here */
3414 /* returns J*exp(-b/a) - data values */
3415 register int i;
3416 double sum0, sum1, b1D1, b2D2, b4D4, wtexpbD, tempcalc, sumbD;
3417 double *expbD, *expbDptr, *wtfactorptr;
3418 register double *Rptr,*bptr;
3419 double D0,D1,D2,D3,D4,D5;
3420
3421 sum0 = sum1 = 0.0;
3422 expbD = malloc (Powell_npts * sizeof (double));
3423 expbDptr = expbD; /* temporary pointers for indexing */
3424 wtfactorptr = wtfactor;
3425 bptr = bmatrix; /* npts of b vectors (nx6) */
3426
3427 D0 = x[0]*x[0]; /* D is used as upper triangular */
3428 D1 = x[0]*x[1];
3429 D2 = x[0]*x[3];
3430 D3 = (x[1]*x[1])+(x[2]*x[2]);
3431 D4 = (x[1]*x[3])+(x[2]*x[4]);
3432 D5 = (x[3]*x[3]) + (x[4]*x[4]) + (x[5]*x[5]);
3433
3434
3435 for (i = 0; i < Powell_npts; i++)
3436 {
3437 /* compute bq.D */
3438 /* bq.D is large dot product of b and D at qth gradient */
3439 /* large dot product for Hilbert algebra */
3440 /* regular dot product is for Hilbert space (vectors only)-
3441 who knew? */
3442 /* calculate explicitly rather than loop to save time */
3443 b1D1 = *(bptr + 1) * D1;
3444 b1D1 += b1D1;
3445 b2D2 = *(bptr + 2) * D2;
3446 b2D2 += b2D2;
3447 b4D4 = *(bptr + 4) * D4;
3448 b4D4 += b4D4;
3449
3450 sumbD = *bptr * D0 + b1D1 + b2D2 + /* bxxDxx + 2bxyDxy + 2bxzDxz + */
3451 (*(bptr + 3) * D3) + /* byyDyy + */
3452 b4D4 + /* 2byzDyz + */
3453 (*(bptr + 5) * D5); /* bzzDzz */
3454
3455 /* exp (-bq.D) */
3456 *expbDptr = exp (-sumbD);
3457 wtexpbD = *(wtfactor + i) * *expbDptr;
3458 sum0 += wtexpbD * Powell_ts[i];
3459 sum1 += wtexpbD * *expbDptr;
3460 expbDptr++;
3461 wtfactorptr++;
3462 bptr += 6; /* increment to next vector of bmatrix */
3463 }
3464
3465 Powell_J = sum0 / sum1;
3466 /* Now compute error functional,E(D,J) and gradient of E with
3467 respect to D ,Ed or F in notes */
3468 /* E(D,J)= 1/2 Sum[wq (J exp(-bq.D) - Iq)^2] */
3469 sum0 = 0.0;
3470 sigma = 0.0; /* standard deviation of noise for weight factors */
3471 expbDptr = expbD;
3472 wtfactorptr = wtfactor;
3473 Rptr = Rvector; /* residuals calculated here - used in
3474 Wt.factor calculations */
3475 for (i = 0; i < Powell_npts; i++)
3476 {
3477 *Rptr = tempcalc = (Powell_J * *expbDptr) - Powell_ts[i];
3478 tempcalc = tempcalc * tempcalc;
3479 sum0 += *wtfactorptr * tempcalc; /* E(D,J) = Sum (wq temp^2) */
3480 sigma += tempcalc; /* standard deviation of noise for
3481 weight factors */
3482 expbDptr++;
3483 wtfactorptr++;
3484 Rptr++;
3485 }
3486
3487 /* sum0 is the error for this iteration */
3488 ED = sum0 / 2; /* this is the error for this iteration */
3489
3490 free (expbD);
3491 return(sum0);
3492 }
3493
3494
3495 /*! compute using optimization method by Powell, 2004*/
ComputeDwithPowell(float * ts,float * val,int npts,int nbriks)3496 static int ComputeDwithPowell(float *ts, float *val, int npts, int nbriks)
3497 /* compute D tensor */
3498 /* ts is input time-wise voxel data, val is output tensor data, npts
3499 is number of time points */
3500 {
3501 /* assumes initial estimate for Dtensor already store in Dvector
3502 and Dmatrix above*/
3503 double *x, tx;
3504 int i, icalls;
3505
3506 ENTRY("ComputeDwithPowell");
3507
3508 Powell_npts = npts;
3509 Powell_ts = ts;
3510
3511 x = (double *)malloc(sizeof(double)*6) ;
3512
3513 /* move data into lower triangular format */
3514 x[0] = sqrt(Dvector.elts[0]);
3515 x[1] = Dvector.elts[1] / x[0];
3516 x[2] = sqrt(Dvector.elts[3] - (x[1]*x[1]));
3517 x[3] = Dvector.elts[2] / x[0];
3518 x[4] = (Dvector.elts[4] - (x[1]*x[3]))/x[2];
3519 x[5] = sqrt(Dvector.elts[5] - (x[3]*x[3])-(x[4]*x[4]));
3520
3521 /*printf("Dvector.elts[] %f %f %f %f %f %f\n",
3522 Dvector.elts[0],Dvector.elts[1],Dvector.elts[2],
3523 Dvector.elts[3],Dvector.elts[4],Dvector.elts[5]);*/
3524 if(debug_briks) {
3525 DT_Powell_optimize_fun(6, x); /* calculate original error */
3526 val[nbriks-2] = ED; /* store original error */
3527 }
3528
3529 tx = TINYNUMBER;
3530 for(i=0;i<6;i++) { /* find the largest element of the initial D tensor */
3531 if(x[i]>tx) tx = x[i];
3532 }
3533
3534 icalls = powell_newuoa( 6 , x , 0.1*tx , 0.000001 * tx , 99999 ,
3535 DT_Powell_optimize_fun ) ;
3536
3537
3538 if(reweight_flag) {
3539 ComputeWtfactors (npts); /* compute new weight factors */
3540 tx = TINYNUMBER;
3541 for(i=0;i<6;i++) { /* find the largest element of the initial D tensor */
3542 if(x[i]>tx) tx = x[i];
3543 }
3544 /* parameters to powell_newuoa (not constrained)s
3545 ndim = 6 Solving for D tensor with 6 elements
3546 x variable for input and output (elements of D tensor)
3547 rstart = 0.1*tx size of search region aoround initial value of x
3548 rend = 0.001*tx size of final search region (desired accuracy)
3549 maxcall = 99999 maximum number times to call cost functin
3550 ufunc = DT_Powell_optimize_fun cost function
3551 */
3552 i = powell_newuoa( 6 , x , 0.1*tx , 0.001 * tx ,
3553 99999 , DT_Powell_optimize_fun ) ;
3554 }
3555
3556 val[0] = x[0]*x[0]; /* D is used as upper triangular */
3557 val[1] = x[0]*x[1];
3558 val[2] = x[0]*x[3];
3559 val[3] = (x[1]*x[1])+(x[2]*x[2]);
3560 val[4] = (x[1]*x[3])+(x[2]*x[4]);
3561 val[5] = (x[3]*x[3]) + (x[4]*x[4]) + (x[5]*x[5]);
3562 /*
3563 printf("D tensor %f %f %f %f %f %f\n",
3564 val[0],val[1],val[2],val[3],val[4],val[5]);
3565 */
3566 if(debug_briks) {
3567 val[nbriks-4] = (float) icalls;
3568 if(icalls<1) {
3569 printf("x values %12.9g %12.9g %12.9g %12.9g %12.9g %12.9g tx %g\n", \
3570 x[0],x[1],x[2],x[3],x[4],x[5],tx );
3571 DT_Powell_optimize_fun(6, x); /* compute J value if not
3572 already computed */
3573 }
3574 val[nbriks-3] = ED;
3575 val[nbriks-1] = Powell_J; /* compute J value */;
3576 }
3577 free(x);
3578
3579 RETURN(icalls);
3580 }
3581
3582 // ========================================================================
3583
3584 // [PT: May, 2017] Calculate goodness-of-fit in two ways, following
3585 // Papadakis et al. (2003, JMRI): Eqs. 1-4.
3586
ChiSq_GOF(THD_3dim_dataset * DWmeas,THD_3dim_dataset * DWfit,int * b0list,float ** ccc,byte * maskp)3587 int ChiSq_GOF( THD_3dim_dataset *DWmeas,
3588 THD_3dim_dataset *DWfit,
3589 int *b0list,
3590 float **ccc,
3591 byte *maskp)
3592 {
3593 int i,j,k;
3594 int Nvox = -1, Ndwi = -1;
3595 int Nwei = 0;
3596 double sumx=0., sumxx=0., diff=0.;
3597 float xm, xf;
3598
3599 ENTRY("ChiSq_GOF");
3600
3601 Nvox = DSET_NVOX(DWmeas);
3602 Ndwi = DSET_NVALS(DWmeas); // len of b0list
3603
3604 // Count number of non-b0s
3605 for(j=0 ; j<Ndwi ; j++)
3606 if( !b0list[j] )
3607 Nwei++;
3608
3609 if( Nwei < 2 )
3610 ERROR_exit("Somehow only %d non-b0 values? "
3611 "How can that be?", Nwei);
3612
3613 // calc the chis
3614 for(i=0 ; i<Nvox ; i++)
3615 if( *(maskp+i) ) {
3616 sumx = 0.;
3617 sumxx = 0.;
3618 diff = 0.;
3619 for(j=0 ; j<Ndwi ; j++) {
3620 if( !b0list[j] ) {
3621 xm = THD_get_voxel(DWmeas, i, j);
3622 xf = THD_get_voxel(DWfit, i, j) - xm;
3623 sumx+= xm;
3624 sumxx+= xm*xm;
3625 diff+= xf * xf;
3626 }
3627 }
3628
3629 if( sumxx > SIG_EPS )
3630 ccc[0][i] = (float) (diff / sumxx); // fine chi_p
3631 else
3632 ccc[0][i] = -1; // bad chi_p
3633
3634 // Calc var; have already guarded against badness here.
3635 sumxx-= sumx*sumx/Nwei;
3636 sumxx/= (Nwei - 1);
3637
3638 if( sumxx > SIG_EPS )
3639 ccc[1][i] = (float) (diff / sumxx); // fine chi_c
3640 else
3641 ccc[1][i] = -1; // bad chi_c
3642 }
3643
3644 INFO_message("Calc'ed chi values. Writing out now.");
3645
3646
3647 RETURN(0);
3648 }
3649