modelfit is an interface to various model fitting routines for diffusion MRI data that fit models of the spin-displacement density function. In particular, modelfit will fit the diffusion tensor to a set of measurements as well as various other models including two or three-tensor models. The program can read input data from a file or can generate synthetic data using various test functions for testing and simulations.
The program sends its output to the standard output by default. The input data must be in voxel order.
Fit the diffusion tensor to a data file SubjectA.Bfloat using linear regression to the log measurements (see Options, -inversion):
modelfit -inputfile SubjectA.Bfloat -model ldt -schemefile A.scheme > DiffTenA.Bdouble
The above command is completely equivalent to:
dtfit SubjectA.Bfloat A.scheme > /tmp/DiffTenA.Bdouble
Fit two cylindrically symmetric tensors to data file SubjectA.Bfloat with starting point determined from the single tensor fit to the log measurements (see Options, --model):
modelfit -inputfile SubjectA.Bfloat -model cylcyl ldt -schemefile A.scheme \
> TwoTensorCSymA.Bdouble
The above command is completely equivalent to:
twotenfit SubjectA.Bfloat A.scheme -cylsym > /tmp/DiffTenA.Bdouble
Run a simulation experiment on data synthesized from standard test function 1, which is a Gaussian displacement density (i.e., the single diffusion-tensor model) with diffusion tensor diag(17, 2, 2)*10^{-10} m^2 s^{-1}:
modelfit -testfunc 1 -model ldt -snr 16 -voxels 10000 -schemefile A.scheme
This simulation runs 10000 trials and reconstructs the tensor by linear regression to the log measurements. Independent noise is added in each trial so that the signal to noise ratio with no diffusion weighting is 16. The program uses the acquisition scheme specified in the file N.scheme (see Options, -schemefile). The command above is equivalent to the pipeline:
datasynth -testfunc 1 -snr 16 -voxels 10000 -schemefile A.scheme | modelfit -model ldt 1 \
-schemefile A.scheme
Run a similar experiment fitting two positive definite tensors to data synthesized from a two-tensor model with diffusion tensors diag(16, 7, 4)*10^{-10} m^2 s^{-1} and diag(3, 18, 9)*10^{-10} m^2 s^{-1} mixed in proportion 6:4.
modelfit -gaussmix 2 16E-10 0 0 7E-10 0 4E-10 0.6 3E-10 0 0 18E-10 0 9E-10 0.4 -model \
pospos ldt -snr 16 -voxels 10000 -schemefile A.scheme
Fit multi-compartment models (see tutorial White Matter Analytic Models) to a data file SubjectA.Bfloat using non linear regression to the log measurements:
modelfit -inputfile SubjectA.Bfloat -inputdatatype float -fitmodel ZEPPELINCYLINDERDOT \
-fitalgorithm LM -schemefile A.scheme -voxels 1 -outputfile ZCD.Bdouble
To perform multiple runs of the fitting:
modelfit -inputfile SubjectA.Bfloat -inputdatatype float -fitmodel ZEPPELINCYLINDERDOT \
-fitalgorithm MULTIRUNLM -samples 1000 -schemefile A.scheme -voxels 1 -noisemodel \
offGauss -sigma 0.0333 -outputfile ZCDmulti.Bfloat
To introduce a specific starting point, use the flag "startpoint", which would make the fitting procedure faster. First you define how many compartments the model you want to fit has by adding "compartment" after the startpoint and the number of compartments next to it, then you define the intra-axonal model and its parameters, then the extra-axonal with its parameters and then the third compartment with its parameters. The following example uses cylinders with gamma distributed radii with a cylindrically symetric tensor and a third compartment of astrocyliders:
modelfit -startpoint compartment 3 gammadistribradiicylinders 0.4 1.8 3E-6 6E-10 1.5 1.5 \
ZEPPELIN 0.2 6E-10 1.5 1.5 1E-10 Astrocylinders 6E-10 2E-6 -inputfile SubjectA.Bfloat \
-inputdatatype float -fitmodel ZEPPELINGDRCYLINDERSASTROCYLINDERS -fitalgorithm LM \
-schemefile A.scheme -voxels 1 -outputfile ZGAc.Bfloat
You can generate synthetic data with these models using the parameter estimates from the model fitting. Here is an example using the ZeppelinCylinderDot model:
datasynth -synthmodel compartment 3 CYLINDERGPD 0.345 8.289E-10 1.293 -4.97 4.3E-6 \
zeppelin 0.331 8.289E-9 1.293 -4.97 3.615E-10 Dot -schemefile 59.scheme -voxels 1 \
-outputfile ZCD.Bfloat
modelfit processes options in command line order.
Single-tensor models
inversion model
-2 restore
1 ldt, dt (default)
2 nldt_pos (nonlinear optimization, constrained to be positive semi-definite)
4 nldt (unconstrained nonlinear optimization)
7 ldt_wtd (weighted linear)
Two-tensor models
Two tensor models are specified in two parts: <two tensor model> <single tensor model> in the same manner as you would specify two digits with -inversion, for example "-model pospos ldt" is equivalent to "-inversion 31". The single tensor model is used to initialize the two-tensor solution and as a fallback if the two-tensor fitting fails. The model names refer to the constraints on the diffusion tensors, which are
pos - tensors must be positive, ie they cannot have negative eigenvalues
cyl - tensors must be positive and also cylindrically symmetric, that is the second
and third eigenvalues must be equal
_eq - enforces that the mixing parameter of the two compartments is 0.5.
inversion model
1? cylcyl
2? cylcyl_eq
3? pospos
4? pospos_eq
5? poscyl
6? poscyl_eq
Three-tensor models Three tensor models are specified in two parts: <three tensor model> <single tensor model>.
inversion model
21? cylcylcyl
22? cylcylcyl_eq
23? pospospos
24? pospospos_eq
25? posposcyl
26? posposcyl_eq
27? poscylcyl
28? poscylcyl_eq
Other models
inversion model
-1 adc
-3 ball_stick
inversion model
-1 adc
-3 ball_stick
Single-tensor inversions
1. Compute the least-squares-fit diffusion tensor to the log measurements by linear regression.
2. Compute the least-squares-fit diffusion tensor to the raw measurements by non-linear optimization using a Levenburg--Marquardt algorithm. The diffusion tensor is constrained to be positive definite by fitting its Cholesky decomposition.
7. Compute the weighted least-squares-fit diffusion tensor to the log measurements by linear regression (see wdtfit(1)).
-2. Compute the diffusion tensor using the RESTORE method.
The program outputs the results voxel by voxel in the same order as the input data file. The output data type is big-endian double-precision (8-byte) floating point values by default. For single-tensor inversions, the output for each voxel contains [exitcode, ln(S(0)), D_xx, D_xy, D_xz, D_yy, D_yz, D_zz], where S(0) is the estimate of the measurement at q=0, the fitted diffusion tensor is D = [D_xx, D_xy, D_xz] [D_xy, D_yy, D_yz] [D_xz, D_yz, D_zz]
and the following exit codes can arise:
-100. Bad data. Inversion could not be performed and all output values are zero.
-1. Background voxel. The voxel was classified as background using a threshold on the mean A^tar(0) measurements, see Options, -bgthresh.
0. No problems.
2. An iterative fitting algorithm (linear or non-linear) failed to converge.
6. Bad data, but the inversion could be performed by substituting a different value for one or more measurements. Usually this means that a measurement was zero or negative so it has no logarithm. For a linear model, this measurement is ignored. Non-linear fitters use the negative value.
Two-tensor inversions All two-tensor inversions use a Levenburg--Marquardt algorithm to fit a mixture of two Gaussian densities to the data. The starting point comes from a single-tensor fit to the data. The indices for all two-tensor inversions are two digit integers. The first digit specifies the type of two-tensor model. The second digit comes from the list of single tensor inversions above and specifies which single tensor fit to use to define the starting point for the optimization.
To determine the starting point, suppose the single tensor fit has eigenvectors e_1, e_2 and e_3 with eigenvalues l_1 >= l_2 >= l_3, respectively. We initialize one component of the two-tensor model, D_1, to have eigenvectors e_1, e_2 and e_3 with eigenvalues (2 l_1 - l_3), l_3 and l_3, respectively. The second component D_2 initially has eigenvectors e_2, e_1 and e_3 with eigenvalues (2 l_2 - l_3, l_3, l_3, respectively. This ensures that the initial D_1 and D_2 have cylindrical symmetry, principal eigenvectors along e_1 and e_2, respectively, and that the average eigenvalue along each e_i is l_i. The mixing parameter is initially 0.5.
In what follws, the question mark "?" denotes a wildcard in the output.
1?. Both diffusion tensors are constrained to be cylindrically symmetric, i.e., each tensor has two equal eigenvalues. An index of 11 means that the program fits two cylindrically symmetric diffusion tensors to the data with the starting point determined using the least-squares-fit diffusion tensor to the log measurements (inversion 1). An index of 12 means that the starting point comes from the diffusion tensor fit to the raw measurements by non-linear optimization (inversion 2).
2?. As 1?, but the mixing parameter is fixed at 0.5.
3?. Both diffusion tensors are constrained to be positive definite.
4?. As 3?, but the mixing parameter is fixed at 0.5.
5?. One diffusion tensor is cylindrically symmetric, the other is positive definite.
6?. As 5?, but the mixing parameter is fixed at 0.5.
For two-tensor inversions, the output for each voxel contains [exitcode, ln(S(0)), N, a_1, D_1xx, D_1xy, D_1xz, D_1yy, D_1yz, D_1zz, a_2, D_2xx, D_2xy, D_2xz, D_2yy, D_2yz, D_2zz], where N is the number of components (here always 2), a_1 is the mixing parameter for D_1 and a_2 is that for D2. This output format is consistent with the multiple-tensor format output by multitenfit(1). The exit codes are as in the single-tensor case.
Three-tensor inversions The three-tensor inversion indices all have three digits, the first of which is a 2. The last digit indicates the kind of single tensor inversion used to obtain the starting point for the optimization that fits the three-tensor model.
The starting point is three cylindrically symmetric diffusion tensors with the largest to smallest eigenvalues in the ratio 8:1. The mixing parameters are all equal. The trace of each tensor is chosen so that the average eigenvalue along each e_i is l_i.
21?. All three diffusion tensors are cylindrically symmetric. For example, an index of 211 means that the program fits three cylindrically symmetric diffusion tensors to the data with the starting point determined using the least-squares-fit diffusion tensor to the log measurements (inversion 1).
22?. As 21?, but the mixing parameters are both fixed at 1/3.
23?. All three diffusion tensors are constrained to be positive definite.
24?. As 23?, but the mixing parameter is fixed at 1/3.
25?. One diffusion tensor is cylindrically symmetric, the other two are positive definite.
26?. As 25?, but the mixing parameter is fixed at 1/3.
27?. Two diffusion tensors are cylindrically symmetric, the other one is positive definite.
28?. As 27?, but the mixing parameter is fixed at 1/3.
For three-tensor inversions, the output for each voxel contains [exitcode, ln(S(0)), N, a_1, D_1xx, D_1xy, D_1xz, D_1yy, D_1yz, D_1zz, a_2, D_2xx, D_2xy, D_2xz, D_2yy, D_2yz, D_2zz, a_3, D_3xx, D_3xy, D_3xz, D_3yy, D_3yz, D_3zz], where N is the number of components (here always 3), a_i is the mixing parameter for D_i. This output format is consistent with the multiple-tensor format output by multitenfit(1). The exit codes are as in the single-tensor case.
Other inversions -1. Fits the apparent diffusion coeffient on the assumption of isotropic Gaussian distributed displacements. The output is [exitcode, ln(S(0)), ADC].
-2. Uses the RESTORE algorithm to fit a single diffusion tensor, see restore(1). Output is [exitcode, ln(S(0)), D_xx, D_xy, D_xz, D_yy, D_yz, D_zz] (see restore(1) for details of exit codes). See restore(1).
-3. Fits the ball and stick partial volume model via nonlinear optimization [Behrens et al, MRM 50:1077-1088, (2003)]. Output is [exitcode, ln(S(0)), d, f, x, y, z]. See ballstickfit(1).
7. Computes the weighted least-squares-fit diffusion tensor to the log measurements by linear regression. See wdtfit(1).
-fixedmodq <M> <N> <Q> <tau> Specifies a spherical acquisition scheme with M measurements with q=0 and N measurements with |q|=Q and diffusion time tau. The N measurements with |q|=Q have unique directions. The program reads in the directions from the files in directory PointSets. The value of N must be in the range 3, ..., 150 or 246. The point sets with 3 up to 150 points minimize the electrostatic energy of pairs of equal and opposite points on the sphere. They are computed using the method outlined in Jansons and Alexander, Inverse Problems, Vol 19, pp. 1031--1046, 2003.. The point set with 246 points is an icosahedral tesselation.
This option is deprecated and should be replaced with -fixedbvalue where possible.