|Home | About | Journals | Submit | Contact Us | Français|
A method for removing arteries that appear bright with intensities similar to white matter in Magnetized Prepared Rapid Gradient Echo images of the ventral medial prefrontal cortex is described. The fast marching method is used to generate a curve within the artery. Then, the largest connected component is selected to segment the artery which is used to mask the image. The surface reconstructed from the masked image yielded cortical thickness maps similar to those generated by manually pruning the arteries from surfaces reconstructed from the original image. The method may be useful in masking vasculature in other cortical regions.
Accurate cortical analysis of gray matter (GM), white matter (WM) volumes and gray/white surfaces can be impeded by the presence of blood vessels. These appear as high intensity structures in gradient echo magnetic resonance (MR) images due to flow related enhancement. As tissues in the desired volume are scanned, their magnetizations become partially saturated, causing longitudinal magnetizations in this region to not recover fully between pulses. Blood outside this region, and therefore not saturated in this manner, then flows into the region. It is imaged and produces a very high signal intensity, causing blood vessels to appear bright in these images [e.g. 1]. A simple thresholding to isolate white matter would result in blood vessels being reconstructed as well (and vice versa). While manual segmentation of the vessels is possible, it is time-consuming and tedious.
Cortical characteristics including cortical thickness metrics, are derived with reference to gray matter/white matter (GM/WM) surfaces, therefore if vessels are not removed from a region of interest they can be misidentified as a GM/WM surface and may artificially influence derived cortical metrics. There is variability in the extent of overlying vasculature in different cortical regions. Certain cortical regions have a significant number of overlying large vessels, which have the potential to induce errors or bias in the neighboring cortical metrics. In the ventral medial prefrontal cortex (VMPFC) there are several major vessels that transverse the region and can cause issues related to cortical characterization. The two main arteries in the region are the callosomarginal and pericallosal arteries . Figure 1 shows the location of the VMPFC relative to the whole brain and these arteries in a sagittal slice of an MR subvolume of the VMPFC. The VMPFC is the lower, central portion of the prefrontal cortex, which is anterior to the precentral sulcus. Studies have found that in clinically depressed subjects there is a significant reduction of gray matter volume in the subgenual prefrontal cortex, [3–8]. There is also data suggesting disruption of the adjacent cortex, the VMPFC (Botteron et al., in preparation and ).
We have previously developed a procedure for cortical thickness analysis of the VMPFC that required manual removal of the arteries by pruning the blood vessels from the reconstructed gray/white surfaces of the VMPFC  but this method is tedious and time-consuming especially in current longitudinal and large-scale neuroimaging studies. Thus a simple, automated method to remove these blood vessels would make cortical analyses of the VMPFC more efficient.
The problem of vessel segmentation in general has been addressed before and multiple methods for solving it exist . Many of these solutions are single scale [12–18]. Several multiscale methods exist as well [19–22]. Many vessel segmentation algorithms involve finding centerlines, or skeletons, of the vessels [23–26]. These algorithms are faster than nonskeleton algorithms, which work without finding a centerline [27–29]. None of these methods have been applied and rigorously tested for specific cortical regions, such as the VMPFC. Thus we develop a skeleton-type method based on a novel application of the fast marching method which is commonly used in level set methods in image analysis [30, 31] and apply it to remove arteries from MR subvolumes of the VMPFC.
The proposed approach finds a path completely contained by the artery and segments it from that path. The problem of finding the skeleton itself was modeled as a minimum path estimation problem. The method to solve the minimal path problem was proposed for 2D problems by Cohen and Kimmel  and later extended to 3D by Deschamps and Cohen . A cost function is defined inside the image such that the desired path is the minimum of the integral of the cost between the two end points. This minimal path problem has been very well studied and a number of solutions exist, such as the Dijkstra method or the Bellman-Ford method . Dijkstra  solved the problem using dynamic programming and graph theory. Cohen and Kimmel  solved the problem by propagating a front between the two end points. This method has the advantage of being geometric, unlike the Dijkstra method, which is purely topologic. As a consequence, it is more precise. The disadvantage of this method is that the cost function used must meet the requirements of the Eikonal equation.
The Eikonal equation can be solved using the level set method  which considers the problem of an evolving interface by increasing the dimensionality of the problem. The initial front is considered the zero level set, and the front at each subsequent time is a higher level set. Thus, a point growing into a circle at uniform speed would be portrayed as a cone with the tip being the zero level set. In the case of monotonically advancing fronts, the surface satisfies the Eikonal equation.
In the blood vessel segmentation method, first a curve is traced inside the blood vessel and is verified by the user. Then the vessel is segmented. Tracing the curve occurs in two steps: calculating a distance map (the front propagation) using a distance function and extracting a minimal path based on this map.
We define the distance in the image between two voxels u and v in two steps. First, for any regular parameterized curve γ(s):[0,1] → 3, such that γ(0) = u and γ(1) = v, define the functional
where x(w) is the image value at location (voxel) w, is the gradient of the curve γ, α and ω are positive constants. Then, the distance between u and v is defined as the minimum value of the functional C over all regular curves γ(s):[0,1] → 3, such that γ(0) = u and γ(1) = v. In the above integrand, |x(γ(s)) − μ|α and ω, are described as data term and regularization term respectively. If there were no data term, C(γ, u, v) would be ω times the length of the curve γ which is minimal when γ is a straight line joining u and v. The data term is small when the image intensity values along γ are close to a constant μ. Deviations from μ are penalized more with large positive α. We experimented with several values of μ, α and ω. Since voxels u and v are selected manually and both x(u) and x(v) are known, we use (x(u) + x(v))/2 for μ which crudely estimates the average intensity along the blood vessel. ω = 1 is chosen and being relatively small means that there is essentially no regularization used. Finally, α = 1 was found to give better results than α = 2.
The user selects a start point and an end point, that is u and v. The Fast Marching algorithm is then used to perform the front propagation. To reduce computation time, the distance is calculated from the starting point, u, outward only until it reached the end point, v, as opposed to over the entire volume. That is, only a partial front propagation was performed.
Once the distance map has been calculated, the actual path must be extracted. A gradient descent on the distance map from the end point is used to trace a path to the start point, which is a global minimum of the distance function. There are several gradient descent methods that can be used. Following Deschamps and Cohen , a steepest gradient descent was chosen because it is easily implemented, consistent and accurate.
After the path has been extracted, the user verifies that all points of the path lie within the blood vessel. If the user is not satisfied, new of start and/or end points are chosen and a new path is delineated. For illustration, Figure 2 shows the curve tracing procedure with initial and start points on the same slice. In reality for the same data in 3D, the initial and start points were located 4 slices apart giving a curve of length of 31.7 mm compared with the Euclidean length of 20.5 mm.
Once a path has been selected, the first step in segmenting the blood vessel is generating a tube around it. The largest possible diameter of the vessel is selected as the radius of this tube because it cannot be assumed that the path is centered in the blood vessel. This value can vary from approximately 1 mm to 3 mm and is selected by the user. The tube at this stage includes the entire blood vessel and some surrounding tissue.
After generating the tube, the next step is thresholding it to remove pixels of intensities below a specified value. This removes surrounding gray matter and other objects that might have been picked up as part of the tube in the area around the blood vessel. Once again the user specifies this value. There is only a lower bound to the range of accepted intensity values because blood vessels have very high intensity values in the T1-weighted MR images.
The final segmentation step consists of selecting the largest connected component. Thresholding removes gray matter, cerebrospinal fluid (CSF) and other darker objects, but neighboring blood vessels and surrounding white matter might still be included in the delineated region. To remove any small pieces that might have been included, the largest connected component is selected and all others discarded. Figure 3 shows the steps of segmenting a blood vessel, following figure 2.
The algorithm was implemented in Matlab1. The front propagation was accomplished using modified functions from the Fast Marching Toolbox2. A blood vessel can be removed in approximately 10 – 15 seconds.
The method was first developed and tested on a large subvolume containing the VMPFC and neighboring structures in one proband (diseased) subject in a neuroimaging study of schizophrenia . Validation was then done on a set of VMPFC subvolumes from ten subjects (three control twin pair, two twin pair with one twin at high risk and other with disease) in a separate and independent neuroimaging study of depression . Both automated and manual methods were applied without knowledge of diagnosis.
The test subvolume was taken using the Magnetized Prepared Rapid Gradient Echo (MPRAGE) sequence from a Siemens 1.5T Vision scanner (with a transmit/receive circularly polarized head coil) with repetition time = 9.7 ms, echo time = 4 ms, flip angle = 10°, section thickness = 1.25 mm, field of view = 224 mm × 256 mm resulting in 1.25×1×1 mm3 voxels. The signed 16 bit MR dataset was compressed to unsigned 8 bit MR dataset using Analyze  by linearly rescaling the voxel intensities such that the voxels with intensity levels at two standard deviations above the mean of white matter (corpus callosum) were mapped to 255 and the voxels with intensity levels at two standard deviations below the mean of CSF (lateral ventricle) were mapped to 0. The means and standard deviations were obtained from sampling of voxels from the respective regions. The MR scan was then interpolated into 0.5 mm × 0.5 mm × 0.5 mm isotropic resolution using trilinear interpolation which allows for more accurate hand segmentation and smoother intensity histograms [3, 10, 37, 39–41].
The ten subvolumes were obtained using T1-weighted MPRAGE sequences in the same Siemens 1.5T Vision scanner and processed using the methods described in Ratnanather et al.  and Botteron et al.  as follows. Each scan had an excitation time of 10–11 ms, a relaxation time of 3–4 ms, with a 256×256 field of view using 160 sections resulting in 1×1×1 mm3 voxels. Three scans were acquired sagittally for each subject. As was done with the test subvolume, Analyze was used to register, sum and interpolate the images, resulting in one image per subject with 0.5×0.5×0.5 mm3 voxels. The scans were averaged to improve gray/white contrast. From these scans, the subvolume containing the VMPFC was defined as tissue in all the medial gyri that were inferior to the corpus callosum and posterior to the first coronal plane that intersects the anterior most portion of the corpus callosum; (brain in standardized AC-PC orientation) the posterior boundary was defined by the plane perpendicular to the AC-PC line where the olfactory nerve nests inside the olfactory sulcus on coronal view.
For the ten VMPFC subvolumes, Bayesian segmentation with an expectation-maximization algorithm was used to classify the CSF, GM and WM tissues . Then a regionally specific intensity correction was performed . Based on the computed gray/white thresholds, gray/white surfaces were reconstructed using isosurface generation algorithms [10, 40]. Labeled Cortical Depth Maps (LCDMs) were then generated from the surfaces [39, 40]. These represent the density of gray matter tissue as a function of gray matter distance from the cortical surface and has been validated for the MPFC .
To determine the accuracy of the method, the blood vessels were segmented manually. Then the automated method was used to segment the blood vessels from the same volume. The distance of each automatically segmented voxel from the closest hand-segmented one was calculated to yield a histogram in figure 4. The Haussdorf (i.e. the maximum) distance was 2.4495 mm (4.8990 voxels). This voxel was an outlier in the lateral, ventral temporal lobe. The blood vessel was immediately next to a bright object, and that object was segmented as part of the vessel. The average distance was 0.1205 mm (0.2410 voxels). 94.0 % of the voxels were within 0.5 mm (1 voxel) and 98.2 % were within 1 mm (2 voxels) of the manually segmented ones. Thus, the method provides an accurate approximation of the ground truth.
Figure 5 shows the superposition of the blood vessels obtained from the original VMPFC gray/white surface on the surface reconstructed with the automated method (cf. Figure 1). Also shown are the axial, coronal and sagittal views of the blood vessel relative to the MRI image. Notice how the blood vessel can compromise the calculation of cortical thickness from LCDMs relative to the green gray/white surface. For validation, we compared the LCDMs generated by the method with those generated by manually pruning of the blood vessel from the triangulated gray/white surface via dynamic programming . Figure 6 illustrates the LCDM curves for both methods plotted by individual subjects demonstrating that both methods gave similar LCDMs. This was statistically verified by the following analysis.
Lilliefors’ test  indicated normality in LCDM distances from manual (p=0.7177) and automated (p=0.8132) methods. Likewise, homogeneity of variances (HOV) for both methods yielded p=0.9091 (Brown-Forsythe test, ). Thus, parametric t-tests were used for pair-wise comparison of the volumes. Volumes from the manual method were larger than those from the automated method (p=0.0002, 2-sided, p=0.0001, manual > automated), but these do not compare trend in shape or size (up to scale) of the LCDMs. Further, correlation analysis indicated that profiles for both methods were similar (Pearson’s r coefficient was 0.9986 with p≈0.0). Kolmogorov-Smirnov tests  for cdfs of volumes showed that the volumes from the two methods were not significantly different (p=0.9945, 2-sided; p=0.6703, manual < automated; manual > automated, p≈1.00) and hence a similar trend in profiles.
LCDM distances generated from manual method were larger than those from the automated method due to gray matter voxels being assigned correct distances relative to the gray/white surface. Pooled LCDM distances for the two methods were found to be significantly non-normal (p≈0.0) and the Brown-Forsythe test for HOV did not hold either (p≈0.0). Kolmogorov-Smirnov tests for cdfs indicated differences in density profiles for LCDMs (p≈0.0, 2-sided; p≈0.0, manual < automated; p≈1.0, manual > automated).
For each subject, LCDM distances for both methods were found to be significantly non-normal using Lilliefors’ tests (p≈0.0). However, Brown-Forsythe tests for HOV were accepted for some subjects and rejected for others. Thus, non-parametric tests such as Wilcoxon (2-sided, > and <) tests  were used for LCDM comparison for each subject. The two methods differed significantly because distances from the manual method were larger than those from the automated method which is not an appropriate test for density profiles. Table 1 shows that Kolmogorov-Smirnov tests comparing density profiles for each subject confirmed distances from the two methods were not significantly different for any subject, indicating that the two methods yielded the same density profile. Similar agreement was obtained for an additional group of ten scans of subjects by an independent analyst.
We have validated the application of a Fast Marching algorithm to segment and remove callosomarginal and pericallosal arteries from MPRAGE subvolumes of the VMPFC, which is the focus of on-going neuroimaging studies of major depressive-disorders. The method is automated and less time consuming than manual pruning of arteries from gray/white surfaces of the VMPFC reconstructed from the original image. Further it is reliable across control and diseased subjects.
Since the blood vessels are relatively large (approximately 3 to 4 voxels in diameter) as well as being bright, the method is not sensitive to the choice of the initial and start points. The cost increased with path length, thus voxels at very large distances from the start point had extremely high costs. Setting α to a higher value caused costs to increase faster as the propagation moved outward. This is a problem for longer blood vessels as the cost function increased to arbitrarily large values at the end point making the path difficult or even impossible to delineate using gradient descent. However, in the case of the VMPFC, long paths were rare in which case the paths were broken up into relatively short ones for bifurcating or multiple vessels.
Other methods such as dynamic programming tracking using a quadratic intensity-based cost function  were initially considered. But dynamic programming was found to be computationally expensive and failed with long paths since it did not attempt to follow the geometric structure of the blood vessel. A region growing method based on the generated path could have been used instead of thresholding and could have worked with high contrast images but would be subject to constraints and would likely not work with multiple or bifurcating blood vessels.
Finally, the method is likely to be useful for isolating vasculature in other cortical regions as well such as the subgenual end of the anterior cingulate cortex  which is incidentally close to the VMPFC and could be adapted and extended to delineating 2D geometric structures in the brain.
This work was supported by the National Institutes of Health (NIH P41-RR15241, P50-MH071616, R01-MH56584, R01-MH626266, P01-AG003991-21, P01-AG026276-01 and R01-EB000975). We thank Suraj Kabadi for technical assistance.
Neeraja Penumetcha received her BS and MSE degrees in Electrical Engineering and Biomedical Engineering at University of Texas Austin and Johns Hopkins University in 2004 and 2006 respectively.
Bruno Jedynak is Assistant Research Professor at the Center for Imaging Science and is affiliated with the Department of Applied Mathematics and Statistics. He is on leave of absence from Université des Sciences et Techologies de Lille.
Malini Hosakere received her B.S. at Loyola College of Maryland and is a senior research technician at the Center for Imaging Science, Johns Hopkins University.
Elvan Ceyhan received his PhD in Applied Mathematics and Statistics at the Johns Hopkins University in 2004. After a year as a postdoctoral fellow at the Center for Imaging Science, Johns Hopkins University, he joined the Department of Mathematics at Koc University, Turkey. His research focus are on statistical methods.
Kelly N Botteron, M.D. is Associate Professor of Psychiatry and Radiology at Washington University at St Louis. Her research interests focus on cortical and subcortical structures such as the prefrontal cortex in major depression disorders.
J. Tilak Ratnanather received his DPhil in mathematics at the University of Oxford. Postdoctoral research at City University and Johns Hopkins University School of Medicine focused on thermal boundary layer separation and mathematical models of cochlear physiology and micromechanics. Since 1998 he has been on the staff and now faculty at the Center for Imaging Science where he works on developing and applying mathematical methods in shape analysis of cortical and sub-cortical structures in a variety of neurodegeneration and neurodevelopmental disorders.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.