Edited by: John Ashburner, University College London Institute of Neurology, UK
Reviewed by: Andreas Deistung, University of Jena, Germany; Russell A. Poldrack, Stanford University, USA
*Correspondence: Phillip G. D. Ward
This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Important physiological information including venous oxygenation can be quantified from MRI scans in a minimally invasive manner. Venous oxygen extraction fraction (OEF) can be derived from quantitative susceptibility mapping (QSM) (Salomir et al.,
Analysis of OEF and other physiological measures from MRI scans requires delineation of regions of interest (ROIs), in the form of binary masks, to identify which voxels represent the structures of interest. The non-regular geometry of biological structures (e.g., vessels) and the regular grid of MRI are fundamentally mismatched, which manifests as partial volume effects at boundaries. Interestingly, in some cases small biological structures such as veins are only visible in MR gradient echo magnitude images due to partial volume which causes signal cancelation between parenchyma and venous blood (Reichenbach et al.,
An intuitive approach to handling partial volume errors in the analysis of large structures is to morphologically erode an ROI, such that the remaining central voxels have minimal or no partial volume effects. Following this erosion a mean or median operation can be applied, however, for veins smaller than three voxels in diameter the erosion operation may leave few (if any) voxels in the mask. The mean of such a small sample is highly sensitive to noise, and approaches the maximum voxel intensity as the sample size decreases. Without the erosion operation, the partial volume effects are likely to cause systematic underestimation by a mean or median approach.
Other approaches to perfect partial volume correction (PPC) include fitting Gaussian mixture models (Shattuck et al.,
Instead, an established approach is to use the maximum intensity voxel (MIV) to represent the vessel (Fan et al.,
Promising vein-specific approaches based on modeling magnetic moments in gradient echo MRI have been recently proposed to estimate the magnetic susceptibility, orientation and size of veins (Hsieh et al.,
Instead of directly modeling MR gradient echo signals, we used QSM techniques to calculate the magnetic susceptibility in veins. The resulting voxel intensities depict local magnetic susceptibility with negligible orientation or echo time dependencies (Fan et al.,
The objectives of this paper are to introduce a novel cylindrical geometry estimation method called Iterative Cylindrical Fitting (ICF) and to validate its use for correcting partial volume errors in small veins. The ICF method resolves the position and radius of small cylindrical structures, to construct partial volume maps, using the voxel intensities of cross-sectional image slices. The partial volume maps enable improved estimation of the true intensities of tubular structures such as veins by using all available voxels, including boundary voxels. The ICF method reduces the effect of noise in the measurements by increasing the number of data points used. In this paper the ICF method is applied to small veins in QSM maps to improve OEF measurements. However, the technique can also be applied to any linear partial volume problem for cylindrical geometries.
The ICF technique takes two inputs, an image (in this case a QSM map) and an initial binary vein mask, and outputs the vein geometry, vein magnetic susceptibility, and a partial volume map for each cross-sectional slice. The initial vein mask may be obtained through manual delineation or by an automated method (Ward et al.,
Cross-sectional slices of the Cartesian image grid, in which a vein appears circular or elliptical, are processed individually. Slices are cropped to contain the vein and immediate neighborhood (in practice 4–10 voxels are included around the vein in all directions). Analysis of each slice provides an estimate of vein radius,
ICF has the following iterative pattern and uses an initial estimate of χbackground from an ROI outside of the dilated binary vein mask:
The initial partial volume values,
A “vein-only image” is produced by subtracting the background component of the image,
The vein geometry (
A new set of non-binary
Steps 2–4 are then repeated until a convergence criterion is met, e.g., a minimal change in
This process is depicted in Figure
The section of a grid line (
The vein radius (
OEF measurements are made using the magnetic susceptibility difference between venous voxels and CSF. In this study we use a uniform hematocrit (Hct) value of 0.4 (Fan et al.,
In numerical simulations χbackground was the mean of the voxels outside of a dilated vein mask. For
Numerical synthetic QSM maps were used to assess systematic error and the performance characteristics of ICF with different vein geometries, magnetic field orientations, noise characteristics and sequence parameters. A description of all default simulation parameters used can be found in Table
Vein position relative to voxel center (in 2 transverse axes) | [−0.5, 0.5]* voxels |
Vein radius (“apparent vein radius” in final image) | 1.3 voxels |
Vein radius (input vein radius prior to truncation) | 8 voxels |
Vein orientation relative to image plane (Θ) | [−45°, 45°]* |
Vein orientation relative to image x-axis (ϕ) | [0°, 45°]* |
Magnetic susceptibility contrast (χvein−χbackground) | 0.30 ppm |
Field strength | 7 tesla |
Transverse decay constant for vein (Deistung et al., |
7.4 ms |
Transverse decay constant for background (Deistung et al., |
33.2 ms |
Longitudinal decay constant for vein (Deistung et al., |
2,587 ms |
Longitudinal decay constant for background (Deistung et al., |
2,132 ms |
Proton density for vein (Deistung et al., |
0.90 |
Proton density for background (Deistung et al., |
0.77 |
Repetition time (TR) | 25 ms |
Echo time (TE) | 7.65 ms |
Gaussian noise (standard deviation, std. noise) | 0.1 |
For a given set of parameters, a perfect cylinder of infinite length was modeled with a radius of 8 voxels. This “ideal case” model was sampled onto a discrete matrix by taking the complex sum of 200 randomly sampled points per voxel. The complex gradient recalled echo signal at each point was simulated with longitudinal and transverse decay (Deistung et al.,
QSM maps were reconstructed using iterative least-squares regression (iLSQR) from the STI-Suite (Wu et al.,
The three experiments explored the effects of echo time, noise, and vein radius. The ranges explored can be found in Table
1 | Echo time (TE) | [3, 24] ms |
2 | Gaussian noise (standard deviation, std. noise) | [0.05, 0.5] |
3 | Vein radius (“apparent vein radius” due to image resolution) | [0.56, 2.1] voxels |
The phase images used to compute the QSM maps contained noise, truncation artifacts, and partial volume effects, along with signal-to-noise profiles influenced by signal decay and decoherent phase. These effects are all influenced by varying vein orientation, vein position, elliptical cross-sections, and vein radius/image resolution. Note that resolution and radius are coupled as the apparent vein radius was controlled by k-space truncation. Experimentally the degree of downsampling varied from approximately 16:1 to 4:1 depending on the required apparent vein radius in different experimental conditions.
The contrast-to-noise ratio (CNR) of each simulation was approximated as the QSM difference between background and vein divided by noise. The noise was estimated by the standard deviation of the QSM map outside of the dilated vein mask. Fit error (ϵfit) was defined as the mean square of the residuals across voxels (
Three analysis methods were employed alongside the ICF method. The first method attempted to mitigate partial volume by taking the MIV from the centermost slice. The second method incorporated no partial volume correction (NPC) by calculating the mean of all voxels with non-zero partial volume. Finally, an ideal method of PPC with ground truth partial volume maps (i.e.,
OEF values were derived according to Equation (8) using the vein magnetic susceptibility estimates from all four methods (ICF, MIV, NPC, and PPC). The error was calculated (ϵ
Numerical simulations of specific vein center positions, zero added noise and a higher resolution matrix were used to examine error from the QSM reconstruction process and the rasterization of the gradient-recalled echo signal.
The complex gradient-recalled echo signal was simulated on a 512 × 512 × 512 matrix with a vein radius of 32 voxels. The matrix was then downsampled by a factor of 24.6 to yield an apparent vein radius of 1.3 voxels. Default experimental parameters were used for all other values (Table
Three vein center positions were simulated, the corner of a voxel, the center of a voxel, and half way between these two positions. Both parallel and perpendicular magnetic fields were applied, resulting in six simulated images that explored asymmetric partial volume profiles and field orientations.
The PPC, ICF, and MIV methods were performed upon this set of images and compared.
All
The vein segments were manually identified and each cross-sectional slice was analyzed. ICF was performed in two passes, the first assuming the vein was perpendicular to the slice, and the second using a vein orientation calculated from a linear fit of the vein positions calculated in the first pass. The convergence and termination criteria were the same as for the synthetic data. The mean radius across all slices was used to estimate
Geometry estimates were examined across all simulations. The mean partial volume map error (ϵρ) was 12.9%±0.3%, the mean vein position error (ϵ
The mean absolute error of OEF measurements (ϵ
1. Echo time | 0° | 6.2 ± 0.3 | 13.4 ± 0.6 | 11.6 ± 0.2 | 5.0 ± 0.2 |
2. Noise | 0° | 11.2 ± 0.5 | 18.4 ± 0.7 | 15.9 ± 0.3 | 8.0 ± 0.5 |
3. Radius | 0° | 5.5 ± 0.3 | 10.5 ± 0.4 | 13.4 ± 0.2 | 3.6 ± 0.2 |
1. Echo time | 90° | 7.1 ± 0.3 | 6.4 ± 0.3 | 15.0 ± 0.1 | 5.6 ± 0.3 |
2. Noise | 90° | 9.8 ± 0.4 | 17.7 ± 0.9 | 15.6 ± 0.2 | 6.2 ± 0.2 |
3. Radius | 90° | 6.5 ± 0.3 | 7.8 ± 0.4 | 14.6 ± 0.2 | 4.5 ± 0.2 |
The ICF error was lower than the NPC error in all experimental conditions, and lower than the MIV error in 5/6 experimental conditions. The MIV method had lower mean error in the remaining experiment (simulations of varied echo time with a perpendicular magnetic field), particularly at longer echo times (Figure
The median OEF error (and inter-quartile range) as a function of varied parameters is displayed in Figure
The MIV approach had a consistent positive bias, whereas the NPC technique has a severe negative bias. The MIV over-estimation was exacerbated by larger vein radius (radius > 1.25 voxels, Figures
ICF computation time was less than 3 s per slice. Twelve maps failed to converge prior to reaching the termination condition (from the total of 1,800 simulated maps). These 12 cases were included in the results above as the value of ϵfit in these cases was similar to the average value for all simulations.
The results of the systematic error investigation are presented in Table
(0,0) | 0° | 1.2% | 0.32 | 1.8% | 0.32 | 1.2% | 0.32 | 0.6% | 1.31 |
(0.25, 0.25) | 0° | 0.7% | 0.31 | 6.8% | 0.39 | 5.5% | 0.37 | −10.8% | 1.16 |
(0.5, 0.5) | 0° | 3.4% | 0.35 | 7.9% | 0.41 | 6.9% | 0.39 | −12.4% | 1.14 |
(0,0) | 90° | 0.1% | 0.30 | 0.7% | 0.31 | −0.1% | 0.30 | 1.4% | 1.32 |
(0.25, 0.25) | 90° | 1.3% | 0.32 | 7.5% | 0.40 | 5.9% | 0.38 | −9.4% | 1.18 |
(0.5, 0.5) | 90° | 4.2% | 0.36 | 8.9% | 0.42 | 6.9% | 0.39 | −11.5% | 1.15 |
An
Frontal cortical vein | 6.6 | 5° | 0.43 ± 0.05 | 0.71 ± 0.08 | 0.22 ± 0.03 | 0.23 ± 0.03 |
Frontal cortical vein | 5.4 | 15° | 0.50 ± 0.06 | 0.84 ± 0.10 | 0.20 ± 0.03 | 0.22 ± 0.06 |
Frontal cortical vein | 7.2 | 20° | 0.55 ± 0.08 | 0.91 ± 0.13 | 0.29 ± 0.07 | 0.28 ± 0.05 |
Frontal cortical vein | 6.6 | 85° | 0.58 ± 0.07 | 0.97 ± 0.12 | 0.25 ± 0.05 | 0.26 ± 0.05 |
Inferior cerebral vein | 6.6 | 10° | 0.58 ± 0.06 | 0.97 ± 0.09 | 0.23 ± 0.05 | 0.22 ± 0.05 |
Inferior cerebral vein | 4.8 | 15° | 0.63 ± 0.05 | 1.05 ± 0.09 | 0.31 ± 0.05 | 0.29 ± 0.05 |
Inferior frontal vein | 6 | 80° | 0.48 ± 0.05 | 0.80 ± 0.08 | 0.28 ± 0.02 | 0.25 ± 0.02 |
Inferior sagittal sinus | 5.4 | 85° | 0.61 ± 0.04 | 1.01 ± 0.07 | 0.26 ± 0.06 | 0.24 ± 0.06 |
Inferior sagittal sinus | 6 | 85° | 0.62 ± 0.06 | 1.03 ± 0.10 | 0.28 ± 0.06 | 0.24 ± 0.04 |
Inferior sagittal sinus | 4.8 | 90° | 0.70 ± 0.07 | 1.17 ± 0.11 | 0.26 ± 0.06 | 0.24 ± 0.05 |
Inferior sagittal sinus | 6.6 | 90° | 0.70 ± 0.08 | 1.17 ± 0.13 | 0.31 ± 0.04 | 0.28 ± 0.04 |
Inferior temporal vein | 7.8 | 85° | 0.44 ± 0.07 | 0.74 ± 0.12 | 0.29 ± 0.03 | 0.23 ± 0.03 |
Inferior temporal vein | 4.8 | 85° | 0.49 ± 0.08 | 0.81 ± 0.13 | 0.28 ± 0.03 | 0.29 ± 0.04 |
Inferior temporal vein | 4.8 | 90° | 0.51 ± 0.08 | 0.85 ± 0.13 | 0.18 ± 0.03 | 0.21 ± 0.03 |
Inferior temporal vein | 4.2 | 35° | 0.55 ± 0.06 | 0.91 ± 0.10 | 0.19 ± 0.04 | 0.16 ± 0.03 |
Internal cerebral vein | 3.6 | 30° | 0.96 ± 0.09 | 1.59 ± 0.15 | 0.21 ± 0.03 | 0.18 ± 0.02 |
Septal vein | 4.8 | 50° | 0.55 ± 0.08 | 0.92 ± 0.14 | 0.35 ± 0.01 | 0.25 ± 0.01 |
Septal vein | 4.8 | 65° | 0.56 ± 0.08 | 0.94 ± 0.14 | 0.35 ± 0.04 | 0.27 ± 0.02 |
Superior cerebral vein | 7.2 | 5° | 0.38 ± 0.05 | 0.63 ± 0.08 | 0.21 ± 0.05 | 0.18 ± 0.03 |
Superior cerebral vein | 7.8 | 20° | 0.49 ± 0.05 | 0.81 ± 0.09 | 0.20 ± 0.03 | 0.22 ± 0.03 |
Superior cerebral vein | 6 | 85° | 0.49 ± 0.08 | 0.81 ± 0.14 | 0.18 ± 0.04 | 0.18 ± 0.05 |
Superior cerebral vein | 6.6 | 15° | 0.51 ± 0.06 | 0.86 ± 0.10 | 0.27 ± 0.05 | 0.17 ± 0.05 |
Superior cerebral vein | 7.2 | 15° | 0.54 ± 0.09 | 0.89 ± 0.14 | 0.20 ± 0.04 | 0.22 ± 0.04 |
Superior cerebral vein | 6 | 15° | 0.58 ± 0.11 | 0.97 ± 0.18 | 0.22 ± 0.04 | 0.20 ± 0.05 |
Superior cerebral vein | 4.8 | 10° | 0.61 ± 0.12 | 1.01 ± 0.19 | 0.22 ± 0.05 | 0.20 ± 0.04 |
Superior cerebral vein | 6.6 | 10° | 0.64 ± 0.07 | 1.06 ± 0.11 | 0.23 ± 0.05 | 0.23 ± 0.05 |
Superior cerebral vein (trib.) | 4.2 | 80° | 0.52 ± 0.04 | 0.87 ± 0.07 | 0.42 ± 0.06 | 0.39 ± 0.05 |
Superior cerebral vein (trib.) | 4.2 | 85° | 0.54 ± 0.05 | 0.89 ± 0.08 | 0.18 ± 0.03 | 0.18 ± 0.03 |
Tentorial vein | 7.2 | 90° | 0.63 ± 0.03 | 1.05 ± 0.04 | 0.23 ± 0.04 | 0.22 ± 0.04 |
Vein of the cingulate sulcus | 6 | 90° | 0.52 ± 0.05 | 0.86 ± 0.08 | 0.44 ± 0.08 | 0.40 ± 0.06 |
Average | 5.8 | 50° | 0.56 ± 0.10 | 0.94 ± 0.17 | 0.35 ± 0.09 | 0.32 ± 0.08 |
The ICF technique provided OEF estimates with lower mean error compared to the MIV and NPC method on synthetic QSM maps. The improved performance of the ICF method can be attributed to the inclusion of information in boundary voxels using a simple two-component linear partial volume model. Four key factors affecting imaging of small structures in QSM were examined: echo time, noise, image resolution and orientation with respect to the main magnetic field.
The close alignment between ICF behavior and PPC behavior, as a function of parameter values, is indicative of low partial volume errors and the accuracy of the ICF approach to estimating vein geometry. However, for smaller veins (
The sensitivity of the MIV method to noise is expected as it uses a single voxel value (maxima) that may be an outlier. In noise and radius simulations, the MIV error in OEF rapidly increased as noise and radius increased. The NPC method depicted a sustained and significant under-estimation of OEF, and a mean error of similar magnitude to the MIV approach in most cases. This under-estimation is expected when taking the mean of small veins where partial volume effects are substantial.
As discussed in the introduction more sophisticated models, such as a Gaussian mixture model, are unsuited as there are too few purely venous voxels to estimate model parameters from. By considering partial volume the ICF method allows more voxels to be included in the analysis, thereby decreasing sensitivity to noise. Additionally, all intermediate steps of the ICF method, such as column and row summations, employ all available voxels yielding robustness to noise in intermediate steps.
The radius experiments show the ICF method outperforming MIV for
Similar results were obtained with both parallel and perpendicular magnetic field orientations. A sustained under-estimation of magnetic susceptibility (by the PPC method) was found when the magnetic field was perpendicular to the vein orientation. This effect was most pronounced in vein radius experiments and indicates an orientation bias in the QSM reconstruction.
In studies which examined larger draining veins, and brain regions using PET 15O, changes in OEF on the order of 10–20% have been found in many conditions including mild traumatic brain injury, multiple-sclerosis, tumor and stroke (Leenders et al.,
The reliability of ICF in
Echo time experiments explored higher phase accumulation (longer echo times). All four methods examined had similar characteristic curves following the expected behavior of phase accumulation and partial volume. As the phase difference inside and outside the cylinder approached π, increased signal cancelation occurred due to destructive interference. This signal loss is in addition to
The PPC method error as a function of echo time, particularly where the magnetic field was parallel to the vein orientation, depicts the effects of phase accumulation and phase decoherence in QSM reconstruction without any partial volume artifacts in the analysis. Intra- and extravascular phase destructively interfere at an echo time of 11–12 ms. At these echo times, the signal-to-noise ratio is reduced due to signal cancelation, which may impact on QSM reconstruction. Similar echo time dependent errors have been identified in other imaging studies that focused on voxel constituents (Sood et al.,
In this work, the ICF method takes a QSM map as an input and therefore results may be affected by susceptibility errors due to the QSM reconstruction itself. Variance was observed in the results from the iLSQR simulations, which may be attributable to the random number sequence used in the sampling, the integration of geometry, and/or the solvers within the iLSQR algorithm. The choice of QSM algorithm and the reconstruction parameters affect the accuracy of the QSM estimates. The systematic error in this study is indicative of this, with non-zero PPC error and a positive bias in maximum voxel intensities in the absence of noise. The systematic error was sensitive to the partial volume profile, i.e., the vein center position within the voxel. These findings suggest residual errors in the QSM reconstruction from the ill-posed nature of the dipole kernel inversion, and from the phase information lost in the complex sum of voxel constituents.
The dependence of MR gradient echo phase images on vessel (or cylinder) orientation to the main field also affects the accuracy of QSM maps due to the differing intra- and extravascular field contributions, and resultant signal-to-noise properties (Li et al.,
Techniques have been proposed to preserve high frequency components and address orientation-specific biases by incorporating structural priors into QSM reconstruction (Tang et al.,
The ICF technique is applicable to cylindrical geometries in any image where partial volume effects are approximately linear (Equation 1). In this work ICF was applied to QSM images, as such χ was used in many of the equations, however it could be replaced with “flow,” “proton density” or in the general case “signal.” ICF may prove useful for analyzing both arterial and venous vessels in different body regions and modalities (other than QSM), as well as other applications where cylindrical geometries are present. As such, other suitable vessel-based applications include CT angiography, spin-echo imaging and phase contrast in arteries for flow velocity. There are also cylindrical image processing applications outside of medical imaging that may also benefit from the technique.
A limitation of ICF is the requirement that the vein intersects two grid lines in each dimension. This is not guaranteed to occur for veins smaller than one voxel in radius, which restricts the smallest veins suitable for ICF. This may be partially mitigated by vein position and tilt angle, as an elliptical cross-section has a higher in-slice radius than the vein itself. For particularly small veins, the contrast is unlikely to be sufficient to initially identify the veins. For larger veins, where more than 2 intersections are possible, it is expected that the centermost grid lines will provide the best estimates in Equation (4) due to the even distribution of area between the segments used.
The lack of ground truth for
Normal aging and other neurological disease processes affect vessel size and metabolism in the brain, and ICF may aid in characterizing and distinguishing these two effects. Demonstration of its clinical utility requires a larger population dataset and application of ICF
The ICF technique employs PPC to improve the accuracy of vein magnetic susceptibility estimates, and subsequent OEF measurements. Compared to the MIV and No Partial volume Correction (NPC) approaches, results from synthetic data demonstrated that ICF provides more accurate estimates of partial volume and venous OEF in the presence of noise and with varying vessel radius and echo time. Application of ICF in quantitation of vessel structure and physiology, such as OEF, will improve the accuracy and reliability of measurements.
PW originally conceived of the method. AN, AF, PR, and PW contributed to the data acquisition. All authors (AF, AN, DD, DB, GE, PR, and PW) contributed to the analysis and interpretation of the data. Drafting was lead by PW and all authors (AF, AN, DB, DD, GE, PR) contributed to the critical revision and final presentation of the work.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer RP declared a shared affiliation, though no other collaboration, with one of the authors AF to the handling Editor, who ensured that the process nevertheless met the standards of a fair and objective review.
The Alzheimer's Australia Dementia Research Foundation (AADRF), the Victorian Life Sciences Computation Initiative (VLSCI), and the Multi-model Australian ScienceS Imaging and Visualization Environment (MASSIVE) supported this work. The authors acknowledge the facilities, and the scientific and technical assistance of the National Imaging Facility at the Melbourne Brain Centre Imaging Unit, The University of Melbourne. Thanks to Ms. Shawna Farquharson for acquiring the MRI data. Preliminary results of this method were previously presented at the International Society of Magnetic Resonance in Medicine annual meeting (Ward et al.,
The Supplementary Material for this article can be found online at: