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