Content-type: text/html
Interpolation is optional, and is different for the probabilistic and deterministic cases. Deterministic trackers interpolate tensor or vector data from the 8 neighbouring tensors and tracking proceeds with a fixed step size. Probabilistic trackers also use a fixed step size, but the interpolation follows the method of Behrens et al (Magnetic Resonance in Medicine, 50:1077-1088, 2003). Noninterpolated tracking uses a modified FACT (Mori et al, Annals Neurology 45:265-269, 1999) method. The stopping criteria is different to the original FACT algorithm (see below).
The tensor deflection (tend) algorithm may be used for deterministic tracking in tensor images (Lazar et al, Human Brain Mapping 18:306-321, 2003). This algorithm is similar to FACT except that the tracking direction in each voxel is deflected by the local diffusion tensor.
Stopping criteria for streamlines are based on on voxel classification, anisotropy and / or curvature. Streamlines will always terminate upon entry to background voxels. Anisotropy thresholds can be used to specify the minimum anisotropy in a voxel. Streamlines will terminate if the anisotropy is below the threshold. The default curvature threshold will terminate streamlines if they curve by more than 80 degrees across the length of one voxel.
The program procstreamlines (see procstreamlines(1)) contains many tools for post-processing of streamlines.
The seed file is an Analyze / Meta / ITK image. The image need not be of the same resolution as the diffusion data, but they must be aligned. A point in physical space, measured from voxel (0,0,0) of the diffusion data, must be the same point in the seed image.
The seed file contains the regions of interest for tracking. Streamlines are seeded at the centre of all voxels with an intensity greater than zero. Voxels that share the same intensity are considered part of the same ROI, they do not have to be contiguous.
The seed, waypoint and target file (see procstreamlines(1)) can be of any data type, but the maximum value in these files must not exceed 2^15 - 1.
A single seed point may be specified on the command line (see OPTIONS).
The output is streamlines in one of three formats: raw binary format, voxel lists, or OOGL format (as used by Geomview).
If -outputroot is specified, streamline output is contained in a single file for each ROI.
for (all rois numbered r) {
open file outputFile
for (all z) {
for (all y) {
for (all x) {
s = seed voxel (x,y,z);
if (s == r) {
for (all iterations i) {
calculate streamline for seed (x,y,z)
write streamline to outputFile
}
}
}
}
}
close outputFile
}
If no output root is specified, then all streamlines are written to a single stream, in the same order as above. The stream is stdout by default but can be set to a file with the -outputfile option.
The coordinate system of tracts is measured in millimeters from the origin of the image, where voxel (0,0,0) is the first to be read from a file in Camino voxel-order format.
The seed file may contain multiple ROIs, with different intensity values. If an output root is specified with -outputroot <root>, then the output is split into different files for each ROI. The output for each ROI with intensity N in the seed file is rootN.suffix, where the suffix depends on the output type. If no output root is given, then all output is written to stdout.
The raw streamline format is 32 bit float. For each streamline, the program outputs the number of points N in the streamline, the index of the seed point, followed by the (x,y,z) coordinates (in mm) of each point: [<N>, <seed point index>, <x_1>, <y_1>, <z_1>,...,<x_numPoints>, <y_N>, <z_N>, <N>,...,<z_N>], where the <seed point index> is the point on the streamline where tracking began.
The voxels format is 16-bit signed integer, and lists the integer indices of all M voxels that the streamline passes through, in the format [<M>, <seed point index>, <tan_x>, <tan_y>, <tan_z>, <ntan_x>, <ntan_y>, <ntan_z>, <vx_1>, <vy_1>, <vz_1>,...,<vx_M>, <vy_M>, <vz_M>, M,...,<vz_M>]. The vector tan is the tangent to the tract at the seed point, pointing from the seed point to the last point in the streamline. The vector ntan is the tangent pointing from the seed point to point 0. Both tan and ntan are unit vectors that have been scaled by 2^15-1.
The OOGL format outputs a LIST of OOGL VECT objects in ASCII format. The colour at each point on the streamline is an RGB vector that describes the local orientation as a combination of red (x), green (y) and blue (z).
The vtkstreamlines program converts Camino streamlines to VTK polylines.
The number of preprocessing steps between raw data to tractography varies according to the data and what kind of tracking is required.
The most common task is to track streamlines through tensor data, and this is the simplest case to deal with. The steps are:
1. Fit the tensors (eg with dtfit).
2. Define the seed ROIs.
3. Run track.
The handedness of the coordinate system used by the scanner may not agree with that used within Camino. If this happens, the anisotropy will be correct but the orientations of the tensors will be incorrect. This problem is often only noticed during tractography. A simple way to check for this is to fit the diffusion tensors and visualize the principal directions with the dtpdview program. If the principal directions appear to be rotated by 180 degrees about the X, Y or Z directions, then the likely cause is that the gradient directions do not agree. This can be remedied by negating the relevant entries in the scheme file.
For streamlines from PAS or Q-ball data, there is an additional step.
1. Calculate the inversion (eg with mesd).
2. Extract the principal directions with sfpeaks.
3. Define the seed ROIs.
4. Run track.
For PICo, one additional step is required per scanner sequence, and one per subject. A lookup table (LUT) of PDF parameters must be created, for tensor data this is done with dtlutgen. The PDFs describe the noise-based uncertainty in principal directions, and therefore the values are dependent on the acquisition scheme and the type of inversion. Once the LUT has been calculated, it can be used for all subjects scanned using the same acquisition parameters and processed with the same tensor fitting routine. See dtlutgen(1) for details of the LUT generation process for tensor data. Once this is done, the steps for a subject are:
1. Fit the tensors (eg with dtfit).
2. Define ROIs.
3. Run picopdfs to construct PDFs in each voxel in the image.
4. Run track with -inputmodel pico and pass the output of
picopdfs as the input file. The -numpds option tells track
the maximum number of principal directions in each voxel. The default
is 1.
PICo can also be performed using the output from multiple-fibre reconstruction techniques such as PASMRI and Q-Ball. See picocalibdata(1) and sflutgen(1) for details on how to create the lookup tables. The rest of the steps remain the same.
When using PICo or PD data, you may still wish to apply an anisotropy threshold. This can be done by passing an anisotropy image with the -anisfile option. In fact, any scalar image may be passed with the -anisfile option; tracking will be restricted to voxels where the intensity is not less than the threshold specified by the -anisthresh option.
Bootstrap tracking requires the raw DWI data and a reconstruction algorithm. The principal direction or directions in each voxel are determined independently for each bootstrap sample of the data.
Currently, diffusion tensor is the only model supported. Both repetition and wild bootstrapping are available. Please see datasynth(1) for more information on the bootstrap techniques.
Using the repetition bootstrap, one and two-tensor models may be fitted to the bootstrap data. The
reconstruction parameters [see modelfit(1) should be passed to track along with the other
parameters. For example, given 4 repeats of a scan SubjectA_[1,2,3,4].Bfloat, (in voxel order),
we can track using repetition bootstrapping and DTI:
track -inputmodel bootstrap -bsdatafiles SubjectA_1.Bfloat SubjectA_2.Bfloat SubjectA_3.Bfloat
SubjectA_4.Bfloat -schemefile A.scheme -inversion 1 -bgmask A_BrainMask.hdr
-iterations 1000 -seedfile ROI.hdr -bsmodel dt > A_bs.Bfloat
To use a two-tensor model, we must pass the voxel classification from voxelclassify.
track -inputmodel bootstrap -bsdatafiles SubjectA_1.Bfloat SubjectA_2.Bfloat SubjectA_3.Bfloat
SubjectA_4.Bfloat -schemefile A.scheme -inversion 21 -voxclassmap A_vc.Bint
-iterations 1000 -seedfile ROI.hdr -bsmodel multitensor > A_bs.Bfloat
The voxel classifications are fixed; they are not re-determined dynamically.
Note that you may pass either -voxclassmap or -bgmask, but not both. If you are using a voxel classification map, the brain / background mask should be passed to voxelclassify. You may always restrict tracking to any volume of the brain by using the -anisfile and -anisthresh options.
Wild bootstrapping requires a single DWI data set. Multi-tensor reconstruction is not supported. The scheme file is required.
track -inputfile SubjectA_1.Bfloat -inputmodel bootstrap -schemefile A.scheme -bgmask A_BrainMask.hdr
-iterations 1000 -seedfile ROI.hdr -wildbsmodel dt > A_wildbs.Bfloat
This method was presented by Friman et al, IEEE-TMI 25:965-978 (2006). In each voxel, we compute the likelihood of the fibre orientation being the axis X, given the data and the model of the data. We wish to sample from
P(X | data) = P(data | X) P(X) / P(data)
in each voxel. We first fit a model to the data (like the diffusion tensor), the model yields m_i, the predicted measurement i given a principal direction X. The observed data y_i is a noisy estimate of m_i. The noise is modelled on the log data as
ln(y_i) = ln(m_i) + epsilon,
where epsilon is Gaussian distributed as N(0, sigma^2 / m_i^2), where sigma^2 is the variance of the noise in the complex MR data. Therefore,
P(data | X) = P(y_1 | m_1)P(y_2 | m_2)...P(y_N | m_N)
where there are N measurements. The prior distribution for all parameters except X is a dirac delta function, so P(data) is the integral of P(data | X) over the sphere. In the case of the diffusion tensor, for example, the priors of S(0) and the tensor eigenvalues L1, and L2 = L3 are fixed around the maximum-likelihood estimate (MLE). The function P(data | X) is then evaluated by setting the tensor principal direction to X and computing the likelihood of the observed data.
The prior on X, P(X), may be set to favor low tract curvature. With the -curvepriork option, the user may set a Watson concentration parameter k. Given a previous tract orientation T, P(X) = W(X, T, k), where k >= 0. The default is k = 0, which is a uniform distribution. Higher values of k increase the sharpness of P(X) around its peak axis T. Suggested values of k are in the range of 0 to 5. You may also use -curvepriorg to implement Friman's curvature prior. Note that a curvature prior does not directly impose a curvature threshold, which may be imposed separately.
An external prior may also be added, in the form of a PICo PDF O(X) defined for each voxel in the image. The full prior is then W(X, T, k)O(X). Pass a PICo image with -extpriorfile.
Example:
track -inputfile SubjectA_1.Bfloat -inputmodel bayesdirac -schemefile A.scheme -bgmask A_BrainMask.hdr
-iterations 1000 -seedfile ROI.hdr > A_bd.Bfloat
This section deals with the various types input to track, which are specified with the -inputmodel option.
The default data type is "float" for input models "bootstrap" and "bayesdirac", and "double" for everything else.
The available input models are:
dt - diffusion tensor data, as produced by modeltfit.
The -anisthresh option may be specified without supplying a separate anisotropy map; in this case the fractional anisotropy of the DT is used.
multitensor - diffusion tensor data, as produced by multitenfit.
The maximum number of tensors in each voxel is specified by the -maxcomponents option. The number of tensors in individual voxels is encoded in the data, so no voxel classification map is required. As with single DT data, a fractional anisotropy mask can be derived from the data. Voxels containing multiple tensors are always included in the mask; single-tensor voxels are included if their anisotropy is above the specified threshold.
pds - principal directions, as produced by sfpeaks.
The maximum number of PDs in each voxel is specified by the -numpds option. The number of PDs in individual voxels is encoded in the data.
ballstick - Ball and stick partial volume model, as produced by ballstickfit.
The ball and stick model [Behrens et al, MRM 50:1077-1088, (2003)] has one principal direction per voxel.
pico - PICo PDFs, as produced by picopdfs.
The maximum number of PDs in each voxel is specified by the -numpds option. The number of PDs in individual voxels is encoded in the data.
bootstrap - raw data for bootstrapping.
Repetition bootstrap data is passed to the program with -bsdatafiles. Wild bootstrap data is passed on standard input or with -inputfile.
bayesdirac - raw data for Bayesian tractography.
This inputmodel is for probabilistic tracking using the Bayesian method
with dirac priors. The assumption of dirac delta priors on the "nuisance parameters"
(parameters other than the fibre orientation) simplifies computation so that the PDF
may be sampled in real time.
The voxel-ordered diffusion-weighted data is passed on standard input or with -inputfile. Relevant options include -curvepriork, -datamodel, and -extpriorfile.
Track within a region of interest defined in an analyze image subAROI.[hdr, img | img.gz ]. The ROI is defined by a collection of voxels with the intensity value 1.
track -inputmodel dt -seedfile subAROI.hdr -anisthresh 0.1 -outputroot A_oneDT_ < SubjectA.oneDT.Bdouble
This outputs A_oneDT_1.Bfloat, containing all streamlines from the ROI. If there were a total of R separate ROIs in the seed file, there would be another output file for each ROI: A_oneDT_{1...R}.Bfloat. In this example we do not specify voxel or data dimensions, so they are taken from the Analyze header of the seed file.
Track from the same ROI, but get PICo connection probability maps. This first requires us to use picopdfs to get the PDF volume A.picopdfs.Bdouble. We then track with
track -inputmodel pico -seedfile subAROI -outputroot A_oneDTpico_ < A.picopdfs.Bdouble
And pass the results to the procstreamlines program
procstreamlines -outputroot A_oneDTpico_ -seedfile subAROI -outputcp < A_oneDT_1.Bfloat
The output here is A_oneDT_nonInt_1_N_1.[hdr,img], where N is 1 for the first seed point and is incremented for each subsequent seed. The last 1 refers to the principal directions of the seed point. If there are P principal directions at a seed point, then the above command will produce output numbered A_oneDT_nonInt_1_N_{1...P}.[hdr,img]. If there were further ROIs numbered 2...R, then there would be output A_oneDT_nonInt_{2...R}_N_{1...P}.[hdr,img]. Note that the number of seeds and principal directions may be different for each ROI.
In this example we want only streamlines that intersect two or more volumes within the brain, Let subA2ROI.img contain the two ROIs, the first defined by voxels with intensity 1, the second by voxels with intensity 2.
We track as before:
track -inputmodel dt -seedfile subA2ROI -outputroot A_oneDT_ < A.picopdfs.Bdouble
And pass the results to the procstreamlines program
cat A_oneDT_1.Bfloat A_oneDT_2.Bfloat | procstreamlines -outputroot A_twoROI_ -outputtracts raw -waypointfile subA2ROI
The procstreamlines program discards all streamlines that do not enter both ROIs.
Equivalently, we could do the above in one line:
cat A.picopdfs.Bdouble | track -inputmodel dt -seedfile subA2ROI | procstreamlines -outputroot A_twoROI_ -outputtracts raw -waypointfile subA2ROI
If the -outputroot option is not specified, then track writes all output to stdout.
The -seedindex option allows seed points within an ROI to be processed independently. The output of parallel PICo tracking will NOT be the same as that of serial PICo. In the parallel case, track is run many times, and a different random number generator is created for each run. When all seed points are processed in series, the random number generator is created only once. By default, the random number generator is seeded from the system time in milliseconds. To avoid multiple tracking operations using the same sequence of random numbers, include a short delay between successive runs of track.
To track using -seedindex, we must know the total number of seeds within the ROI. The program countseeds does this for us:
countseeds subA2ROI
We can then process each seed in parallel with the bash command, subject to one constraint: because track outputs all streamlines from an ROI to the same file. We therefore have to name the output root from each seed uniquely:
for seed in `seq 1 $(countseeds subA2ROI)`; do
track -inputmodel pico -datadims 128 128 60 -voxeldims 1.7 1.7
2.3 -seedfile subA2ROI -outputroot A_oneDT_${seed}_ -seedindex $seed < A.picopdfs.Bdouble
done
This command produces the files A_oneDT_{1,2,...,N}_R_1.Bfloat for each of the N seeds in the seed file (the program determines the ROI R to which each seed belongs). These can be concatenated to produce the same output we get from tracking without parallel processing.
If we do not wish to process each seed in parallel, we can process each ROI in parallel instead. The -regionindex option is similar to the -seedindex option, but all seeds in each ROI are processed in series. Therefore, if there are four ROIs in the seed file subA2ROI.img, then we use a similar command to process each ROI in parallel:
for roi in `seq 1 4`; do
track -inputmodel pico -datadims 128 128 60 -voxeldims 1.7 1.7
2.3 -seedfile subA2ROI -outputroot A_oneDT_ -regionindex $roi < A.picopdfs.Bdouble
done
To produce connection probability maps directly, the procedure is
for roi in `seq 1 4`; do
track -inputmodel pico -datadims 128 128 60 -voxeldims 1.7 1.7
2.3 -seedfile subA2ROI -regionindex $roi < A.picopdfs.Bdouble | procstreamlines
-seedfile subA2ROI -regionindex $roi -outputcp -outputroot A_oneDT_cp_
done
The advantage of this approach is that there is no need to concatenate streamline output after tracking. In some circumstances, this method may also be more efficient, since the program is only initialized once per region instead of once per seed.
The following list details the options pertaining to the input data, the tractography parameters, the output, and the PICo parameters.
dt - diffusion tensor data, as produced by dtfit.
multitensor - diffusion tensor data, as produced by multitenfit. Presently, track supports a maximum of two tensors per voxel for PICo. Three tensors per voxel is allowed for streamline tracking. The maximum number of tensor components in a voxel is controlled by the maxcomponents option, and the default is 2.
pds - principal directions, as produced by sfpeaks. The maximum number of PDs in a voxel (if different from the default of 3) is specified by the numpds option.
pico - PICo PDFs, as produced by picopdfs.
bootstrap - raw data to be bootstrapped.
ballstick - as produced by ballstickfit.
bayesdirac - raw data for Bayesian tracking.
The maximum number of PDs in a voxel for input models pds and pico. The default is 3 for input model pds and 1 for input model pico. This option determines the size of the voxels in the input file and does not affect tracking. For tensor data, this option is not used. If the input model is "dt", the number of PDs in each voxel is 1. If the input model is "multitensor", the number of PDs is specified by the -maxcomponents option.
The maximum number of tensor components in a voxel. This determines the size of the input file and does not say anything about the voxel classification. The default is 2 if the input model is multitensor and 1 if the input model is dt.
For probabilistic methods, the "interpolated" tracking direction is chosen from the voxels surrounding the current point. The algorithm for choosing a direction is described by Behrens et al Magnetic Resonance in Medicine, 50:1077-1088, 2003.
1.0 2.0 3.0
4.0 5.0 6.0
...
The entire list of points is treated as a single ROI by track.
Output streamlines, in one of the following formats: raw binary, voxels, oogl, or ooglbinary. The default is to output raw tracts.
The raw streamline format is 32 bit float. For each streamline, the program outputs the number of points N in the streamline, the index of the seed point, followed by the (x,y,z) coordinates (in mm) of each point: [<N>, <seed point index>, <x_1>, <y_1>, <z_1>,...,<x_numPoints>, <y_N>, <z_N>, <N>,...,<z_N>], where the <seed point index> is the point on the streamline where tracking began.
The voxels format is 16-bit signed integer, and lists the integer indices of all M voxels that the streamline passes through, in the format [<M>, <seed point index>, <vx_1>, <vy_1>, <vz_1>,...,<vx_M>, <vy_M>, <vz_M>, M,...,<vz_M>]. The voxels are defined in the seed space.
The OOGL format option, oogl, outputs a LIST of OOGL VECT objects in ASCII format. The colour at each point on the streamline is an RGB vector that describes the local orientation as a combination of red (x), green (y) and blue (z). The ooglbinary option outputs the same objects in the OOGL binary format.