PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Magn Reson Med. Author manuscript; available in PMC 2018 January 1.
Published in final edited form as:
PMCID: PMC4947568
NIHMSID: NIHMS741390

Automatic recognition of subject-specific cerebrovascular trees

Abstract

Purpose

An image filter designed for reconstructing cerebrovascular trees from MR images is described. Current imaging techniques capture major cerebral vessels reliably, but often fail to detect small vessels, whose contrast is suppressed due to limited resolution, slow blood flow rate, and distortions around bifurcations or non-vascular structures. An incomplete view of angioarchitecture limits the information available to physicians.

Methods

A novel Hessian-based filter for contrast-enhancement in MR angiography and venography for blood vessel reconstruction without introducing dangling segments is presented. We quantify filter performance with receiver-operating-characteristic and dice-similarity-coefficient analysis. Total extracted vascular length, number-of-segments, volume, surface-to-distance, and positional error are calculated for validation.

Results

Reconstruction of cerebrovascular trees from MR images of six volunteers show that the new filter renders more complete representations of subject-specific cerebrovascular networks. Validation with phantom models shows the filter correctly detects blood vessels across all length scales without failing at bifurcations or distorting diameters.

Conclusion

The novel filter can potentially improve the diagnosis of cerebrovascular diseases by delivering metrics and anatomy of the vasculature. It also facilitates the automated analysis of large datasets by computing biometrics free of operator subjectivity. The high quality reconstruction enables computational mesh generation for subject-specific hemodynamic simulations.

Keywords: vessels, segmentation, image processing

Introduction

Medical imaging is a common technique for non-invasive diagnosis of cerebrovascular diseases (CVD), which constitute the fourth leading cause of death in the United States (1). Magnetic resonance imaging (MRI) enables physicians to diagnose CVDs which often have no apparent symptoms. Crisp images also provide data for treatment planning (2). In raw magnetic resonance angiography (MRA) and venography (MRV), anatomical structures including gray and white matter overlap with the patient's cerebral vasculature hindering blood vessel detection. Small vessels like communicating arteries in the Circle of Willis are often invisible in raw MRA images due to their relatively weak contrast compared to stronger signals from larger arteries, such as the middle cerebral artery (MCA). Therefore, a method of separating different components of the MR signal without missing portions of the vessel network is needed. A vesselness filter is a digital image processing procedure for vessel contrast enhancement that aims at accentuating blood vessel contrast while eliminating other structures from the image.

In addition to vessel filtering for diagnosis, there is also a growing need for patient-specific vasculature reconstructions for computational fluid dynamics (CFD) studies. Researchers and physicians are beginning to demonstrate the benefits of cerebral hemodynamic simulations to acquire a more profound understanding of why CVDs occur, better plan endovascular interventions such as bypass surgery (3,4), or assess aneurysm (5) and arteriovenous malformation rupture risk (6). All CFD studies require subject-specific computational meshes delineating local features such as the carotid bifurcation (710), or extending to larger portions of the arterial tree, such as the Circle of Willis (1113).

Numerous prior studies presented excellent Hessian-based filters for vessel segmentation (14-20). Eigenvalues of the Hessian matrices of the image intensity map can guide automatic algorithms to track network centerlines. Sato (14) proposed an algorithm based on experimental analysis of ideal tubes. The Frangi filter (15,16) has three adjustable parameters to select the best tradeoff between noise elimination and shape sensitivity. Erdt (17,18) developed a purely analytical function for fully automatic reconstruction.

All Hessian-based filters discussed so far amplify the contrast between tubular objects and their surroundings, thus enhancing elongated blood vessels while suppressing other anatomical features and noise. However, the eigenvalues of the Hessian experience sign changes in addition to wide variations in magnitude at bifurcations. This situation often causes existing filters to erroneously eliminate bifurcations, introducing dangling segments in the reconstructed network. The Shikata filter (19,20) is able to preserve network connectivity by also enhancing the contrast from plate or blob-like shapes, which favors bifurcations. Unfortunately, Shikata's improvement for bifurcations has the unwanted side effect of suppressing neighboring segments in the vicinity as pointed out by Drechsler (21). This filter behavior is undesirable because closeness of bifurcations to other branches frequently occurs in the pial network of the cerebral cortex. It would be ideal to design a novel filter, which accentuates tubular shapes, while maintaining network connectivity in bifurcations without losing neighboring segments. Accordingly, our proposed algorithm aims to achieve the following three objectives:

  • Enhance the contrast of blood vessels across the entire spectrum of length scales from the internal carotid artery or superior sagittal sinus down to the pial network.
  • Suppress the intensity of non-vascular tissues and reduce noise in the image.
  • Preserve the connectivity of the reconstructed cerebrovascular network especially at bifurcations without creating dangling segments or suppressing neighboring branches.

This paper is organized as follows. First, MRA/MRV image acquisition parameters for vessel detection are presented. We then introduce a novel multi-scale composite filter (MCF) algorithm and its implementation. We compare its performance against previous algorithms in two phantom studies as well as in six human subject data sets. Finally, we discuss findings and indicate future research directions.

Methods

MR images were acquired from healthy volunteers to demonstrate the ability to create subject-specific reconstructions of their cerebral vasculature automatically. Although the methodology can be applied to any number of images, data and results for only six volunteers are presented to discuss performance metrics in detail. Moreover, filter performance was analyzed statistically in two additional phantom studies.

Image Acquisition

Six healthy human subjects with no known cerebral vascular disease were recruited and underwent MR imaging on a General Electric 3T MR750 scanner using a 32 channel phased array coil (Nova Medical, Inc., Wilmington, MA). MR angiograms and venograms were acquired under Institutional Review Board approval at our institution. Before the scan, stability pads were provided to reduce motion artifacts. No motion artifacts affecting the scan were observed in all six cases. For MR angiography, a 3D TOF pulse sequence was used with the following parameters: TR=26ms, TE=3.4ms, NEX=1, Flip Angle=18°, acceleration factor=2, number of slab=4, magnetization transfer=on, matrix size=512×512×408, voxel size=0.39×0.39×0.3mm3, scan time=30min.

MR venography was performed using a 2D INHANCE pulse sequence with the following imaging parameters: TR=18.5ms, TE=5.65ms, NEX=1, Flip Angle=8°, matrix size=512×248×512, voxel size=0.47×0.8×0.47mm3, scan time=15min. Automatic rigid coregistration with a one-plus-one evolutionary optimizer of the MRA and MRV was performed with the final voxel size of 0.39×0.39×0.3mm3. The MRA and MRV captured major branches of the cerebral arterial and venous systems, respectively.

Vesselness Filters

We first processed all angiograms and venograms with our MATLAB implementation of existing Hessian-based filters including Sato (14), Frangi (15,16), Erdt (17,18), and Shikata (19, 20). The mathematical background for each filter is given in Supporting Material. All filter functions use the raw 3D image intensity map I(x, y, z) as input, where x,y,z represent the spatial coordinates. The Hessian matrix, H(x, y, z, σ) of the image measures the curvature of the intensity at each voxel. Construction of Hessian matrices involved computations of gradients and second order derivatives, which are functions of voxel spacing in each direction, [partial differential]x, [partial differential]y, [partial differential]z. To reduce noise in derivatives, Gaussian filtering was performed on the raw image as given in Eq. 1. Noise in directional derivatives around arteries shown in Supporting Fig. S1 is best suppressed by Gaussian filtering with its variance (σ-parameter) equal or close to the vessel diameter. Therefore, vesselness filters sensitive to different diameter ranges should use corresponding σ-values. The human cerebral vasculature covers multiple length scales from the 5mm range in the carotid to about 200μm in pial arteries. To address the need for multi-scale smoothing, multiple runs of Gaussian filtering with different σ-values was performed. Different Gaussian σ-settings produced distinct intensity maps, Ĩ(x, y, z, σ). Hessian matrices, H(x, y, z, σ), for each scale were computed as in Eq. 2. Furthermore, the three eigenvalues [lambda with tilde]1, [lambda with tilde]2, [lambda with tilde]3, of each Hessian matrix (x, y, z, σ), were calculated for each voxel. Supporting table S1 shows the range of eigenvalues characterizing different geometrical shapes. For example, an ideal cylinder has [lambda with tilde]1 = 0, [lambda with tilde]2 < 0, [lambda with tilde]3 < 0. The eigenvalues permit the calculation of the image enhancement response f(x, y, z, σ, [lambda with tilde]) for each of the five Hessian based filter algorithms. Each image filtered with a different σ-value leads to a different f-map. The final enhanced image intensity If(x, y, z) is assembled by selecting the highest voxel values among the different f-maps in Eq. 3. By electing locally the scale σ which achieves the maximum f-value, the filtered image benefits from the optimal contrast enhancement in each voxel. For example, a map of the corresponding σ-values that optimally maximized intensity in an axial image projection is shown in Supporting Fig. S2.

Gσ(x,y,z)=1(σ2π)3e(x2+y2+z2)2σ2I(x,y,z,σ)=Gσ(x,y,z)I(x,y,z)
Eq.1
H(x,y,z,σ)=[2Ix22Ixy2Ixz2Iyx2Iy22Iyz2Izx2Izy2Iz2]
Eq.2
λ1,λ2,λ3eig{H(x,y,z,σ)}If(x,y,z)=maxσf(x,y,z,σ,λ)
Eq.3

Specifically, for MRA we chose fixed parameter settings, σ = {1, 2, 3, 4, 5mm}. This range spans the relevant length scale of the human arterial tree. It is worth noting that Gaussian σ-settings are not adjustable parameters, because physiological diameter ranges do not vary significantly between subjects. The recommended range gives optimal results in all human MRA we performed. In MRV, a fixed range of σ = {1, 2, 3, 4, 5, 6, 7, 8mm} was used. A wider range than for arteries was needed to cover the larger dimensions of veins.

The final scalar field If(x, y, z) reflects the intensity map of the enhanced image with multi-scale Gaussian filtering for each of the five filter algorithms. Some filters also have adjustable parameters as summarized in Supporting table S2. The mathematical background of a novel multi-scale composite filter (MCF) is presented next.

Multi-scale Composite Filter Function (MCF)

We discovered that gap creation in vesselness filters can be avoided by suitable combination of contrast enhancement of the Frangi algorithm with the bifurcation treatment of the Shikata algorithm. The novel multi-scale composite filter (MCF) has two stages as given in Eqs. 4 and Eqs. 5. Hessian matrices for each length scale, H(x, y, z, σ) are generated from smoothed intensity maps Ĩ(x, y, z, σ). For each voxel, its three eigenvalues are calculated and ordered by magnitude |[lambda with tilde]1| ≤ |[lambda with tilde]2| ≤ |[lambda with tilde]3|.

Step 1 uses the auxiliary scalar function, C1, which is computed from the eigenvalues of the Hessian H(x, y, z, σ) according to a modification of the Frangi algorithm (15,16). For each scale, the C1 value at each voxel is normalized with respect to the maximum intensity in the entire image as given by the maximum norm in Eqs. 4. Again, the maximum intensity among the different length scales is selected for each voxel. Thus, C1 filtering produces a new normalized intensity map IC1 (x, y, z).

To preserve bifurcations without cutting off segments, step 2 applies a modified Shikata procedure onto the normalized intensity map IC1 (x, y, z). After a second pass of multiscale smoothing identical to step 1, the second largest eigenvalue exhibits consistent magnitude along a connected network, while the smallest eigenvalue often experiences sign changes as discussed in Shikata et al (19,20). Accordingly, the auxiliary scalar function, C2, uses only the second largest eigenvalue λ͌2 of Hessians, (x, y, z, σ), derived from the C1-filtered normalized intensity map C1 (x, y, z, σ). Moreover, positive second eigenvalues, λ͌2 > 0, mark dark objects which cannot be blood vessels in MRA. Therefore, we set the final intensity to zero whenever λ͌2 > 0 as indicated in Eqs. 5. Furthermore, we add an additional inequality λ͌1/λ͌ < 0.25. If the inequality holds, the filter response is set to zero to suppress plate-like structures. We found in all case studies that the two inequality constraints do not induce speckle artifacts.

Step 1:

C1(H)=(1era22a2)erb22b2(1erS22c2)(1era22a2)erb22b2(1erS22c2)rA=|λ2||λ3|;rB=|λ1||λ2λ3|;rS=λ12+λ22+λ32IC1(x,y,z)maxσC1(H(x,y,z,σ))
Eq.4

Step 2:

IC1(x,y,z,σ)=Gσ(x,y,z)IC1(x,y,z)λ1,λ2,λ3,eig{H(x,y,z,σ)}C2(H(x,y,z,σ))={0λ2>00λ1/λ2<0.25σ2λ2σ2λ2else;IMCF=maxσC2(H(x,y,z,σ))
Eq.5

For each length scale, the sequential application of C1 and C2 filter components produces a final scalar value at each voxel. The length scale leading to the highest intensity is selected for the final intensity IMCF for each voxel. The three dimensional matrix of maxima, obtained perhaps with different σ for each voxel, gives the vessel enhanced intensity image. The novel filter also has three adjustable parameters (a, b, c) to allow for filter tuning.

Note that both the local C1 and C2 functions are normalized with highest intensity found in the entire image. Normalization also ensures that the relative contributions from the two auxiliary functions are commensurate.

The value of the scalar IMCF at each voxel signifies the probability that it belongs to the vascular network. At vessel centerlines, the IMCF values assume a local maximum. The centerline of the vascular network is the locus of extrema in the IMCF map. The centerline trajectories form three dimensional space curves that trace the subject's vascular trees.

Computational Mesh Generation and Validation

Filtered images were processed to create logically connected networks for vectorized display or computational purposes. We applied the filters and the segmentation process to the human data for reconstructing computational meshes representing subject-specific cerebral vascular trees. Subject-specific computational meshes are needed for individualized high quality visualization of the angioarchitecture (27), solid mechanic injury simulations (28) or hemodynamic CFD simulations (29,30). The workflow for subject-specific automatic mesh generation is depicted in Fig. 1.

Figure 1
Workflow for automatic computational mesh generation of subject-specific vascular trees. First, raw MRA/MRV images are processed with vesselness filters for contrast enhancement of the cerebral vessel network. The filtered image crisply outlines the vasculature ...

First, we create a distance map between the vessel centerlines and the vessel surface, where a zero level set corresponds to the vessel walls. The fast marching algorithm with cutoff threshold (intensity threshold=0.01) was used to create a binary mask of the connected vascular tree (22). The binary mask delineates a connected domain. An active contour procedure is used to calculate the geodesic level of each interior point, in which the computed levels measure the distances of a point to the vessel wall, whose outer edges are assigned a zero level (23). With the entire image domain now parameterized in terms of wall distance, the marching cubes algorithm retrieves smooth physical coordinates of the vessel walls (40). Note that physical coordinates for the surface obtained by marching cubes is not limited to the original discrete voxel grid spacing. Optionally, isolated dangling surfaces that are not connected to the main arterial tree can be discarded. This situation typically affects external surface arteries. Because of the special attention to preserving connectivity with the C2 filter, dangling segments do not occur in the territory of the main cerebral arteries.

The surface enclosing the arterial tree serves as input for unstructured mesh generation (25). Alternatively, parametric volumetric meshes (26) can be created with the help of centerline and diameter information. Coordinates of the centerline trajectory as well as corresponding vessel radii were tracked with the maximal inscribed spheres method (24). By using the fast marching algorithm with cutoff threshold, our method could eliminate the tendency to overestimate diameters in the maximal inscribed sphere method (39). With precise centerlines and radii, body-fitted smooth hexahedral computational meshes of the entire arterial trees were generated (33). Venous trees were created similarly from MRV data.

Validation

Volume, number of segments, as well as total length of the reconstructed arterial and venous trees are used as metrics for quantitative comparison of the filter performance. A segment is defined as a cylindrical connection between two adjacent centerline points. For further validation, we created two phantoms and rigorously tested filter performance in detecting blood vessel, connectivity at bifurcations, and diameter reconstruction against a ground truth reference.

Results

Image enhancement of MRA and MRV data from healthy volunteers

All filters were custom-coded in MATLAB on an Intel Xeon CPU E5645(x00040)2.40 GHz computer. The CPU time for existing filters is approximately 60s. The sequential implementation of the MCF consumes less than twice the CPU time. The top row of Fig. 2 shows an axial, sagittal, and coronal view of the original MRA image acquired from the first of the six volunteers. Each view is a maximum intensity projection displayed with its maximum and minimum filter response as the contrast window. The middle cerebral arteries (MCA) and external cerebral arteries appear the brightest. The sagittal view also depicts the major branches of the anterior (ACA) and posterior cerebral artery (PCA) territories. The Circle of Willis (COW) can only be seen partially in the coronal view. The thin anterior communicating (ACOM) and posterior communicating arteries (PCOM) are barely visible in the original MRA as highlighted in the magnified view in Fig. 2 with arrows indicating their usual anatomical position. In general, smaller arteries, particularly in the ACA and PCA territory, are hard to discern in the raw image.

Figure 2
Maximum intensity projections of the raw MRA and five filter responses (A-E) in axial, sagittal, and coronal views as well as zoomed in detail around the circle of Willis. Top row. Raw MRA images. Filtered images from: A. Sato, B. Frangi, C. Shikata, ...

The images in Fig. 2A to 2D display axial, sagittal and coronal views produced by Sato, Frangi, Shikata and Erdt filters, respectively. The Sato filter in Fig. 2A suppresses noise and soft tissues, while enhancing blood vessels as seen most clearly in the coronal view. The Frangi algorithm offers tunable control. Its application with parameters, a=0.5, b=0.5, c=0.05, in Fig. 2B sharpens vessels better than the Sato filtered image. Structures not belonging to the cerebrovascular network such as gray and white matter are successfully suppressed in the Shikata filter as shown in Fig. 2C. The Erdt filter, which has previously been successful in computed tomography images, did not improve the vascular network contrast substantially as depicted in Fig. 2D.

Results for the novel MCF with parameters, a=0.5, b=0.5, c=0.05, are shown in Fig. 2E. These parameters can easily be set in the MATLAB implementation and gave the best results in all six case studies. The MCF clearly produces the brightest outline of arterial trees. The magnified view of the COW also clearly depicts the PCOMs, which are entirely missed in the raw images, and are hardly discernable in prior vesselness filters. This result demonstrates that MCF not merely sharpens vessel contrast, but more importantly crisply delineates small vessels whose signals are too faint to be seen with prior filtered or raw images. MCF is also able to reconstruct a good portion of the cortical surface network, which consists of small bifurcating pial arteries that are entirely missed in all existing filter techniques. The bright and crisp delineation of arterial trees show that the vessel intensity is enhanced for all diameters ranging from the large carotid arteries down to small pial vessels of about 400μm in diameter. The novel MCF also preserves the contrast around bifurcations, which is an essential feature for maintaining network connectivity without introducing dangling segments.

Statistical analysis of filter performance

To systematically evaluate the performance of the filters and to quantify the advantage of the novel MCF, we conducted receiver operating characteristic (ROC) analysis, dice similarity coefficient analysis, surface distance errors, and positional offsets with two artificially generated image phantoms. For the ROC analysis, thresholds, Ithreshold, covering the entire range between zero and unity were chosen, and the number of correctly detected voxel belonging to vessels were counted. We also counted false positives, which is equal to voxels erroneously assigned to vessels. High threshold levels decreased the number of false positives, but miss many vessels. In general, there is a trade-off between specificity (avoidance of false positives) and recognition rate (fraction of detected vessels). The phantoms offer exact vessel coordinates and diameters as ground truth reference for objectively evaluating filter performance.

3D artificial tube network (lattice phantom)

The first phantom consists of a rectangular lattice of tubular objects embedded in a spherical domain to emulate a configuration of blood vessels embedded inside a brain image. Fig. 3 depicts the original phantom lattice, the filter responses and the ROC curves obtained for each filter. Fig. 3A shows that the Sato response loses signals in bifurcations and that it was unable to suppress the surrounding sphere. The Frangi filter in Fig. 3B suffers from artifacts at the intersection of the sphere with tubes. The quality of intensity enhancement achieved by the Shikata procedure is poor as shown in Fig. 3C. The lattice structure detected by the Erdt filter in Fig. 3D is severely distorted by the surrounding sphere. The image filtered by MCF in Fig. 3E renders the lattice with bright intensity. There is slight blurring where tubes and background intersect. The ROC curve confirms that the MCF filter has the best statistical performance with the analysis of area under the curve of 0.94. Moreover, for every threshold, the MCF curve is pareto-optimal, meaning that it has higher or equal true positives, for a given false positive level. The MCF also has the highest dice coefficient of 0.89 among all filters, including Frangi with 0.85, Sato with 0.82, Shikata with 0.75 and Erdt with 0.63.

Figure 3
Statistical analysis of filter performance in a 3D lattice phantom. Top left, raw phantom image representing a tubular network inside a spherical background. A. Sato filter response. B. Frangi filter response. C. Shikata filter response. D. Erdt filter ...

We further analyzed the performance of all filters with respect to connectivity in this case study with ten bifurcations. The Sato and Frangi filters give faint signals for all bifurcations as seen in Fig. 3A and 3B. Moreover, the two filter algorithms introduce dark stripe artifacts close to bifurcations. The Shikata, Erdt and MCF filters have no stripe artifacts. We counted the number of correctly detected bifurcations as well as false positives with respect to different intensity thresholds that mark a cutoff above which a voxel is considered to belong to a blood vessel. We count a bifurcation as connected when all six neighboring voxels are categorized as members of blood vessels. This analysis summarized in Table 1 demonstrates that the MCF correctly labeled all ten bifurcations without any false positives for the selected intensity thresholds.

Table 1
Number of bifurcations detected from filter responses in phantom studies. Intensity thresholds mark the cutoff point in intensity above which a voxel is considered to belong to a blood vessel.

We further used the ground truth diameter information to measure the fidelity of diameter reconstruction. All tubes in the artificial tube network had a diameter of 5 voxels. Diameter reconstruction with the Sato filter results in average of 5.2 with standard deviation of 0.3, the Frangi filter with 5.1±0.2, the Shikata filter with 5.3±0.5, and the Erdt filter with 5.4±0.2. The MCF with 5.1±0.1 was closest to the correct value.

3D microcirculation (network phantom)

The second phantom establishes a more physiologically realistic ground truth scenario against which filter performances were evaluated. It consists of a 3D microcirculatory network overlapped with an oval shaped backdrop to resemble arterioles within a cerebral sulcus. Bifurcating trees and branches were generated by a constructive growth algorithm to emulate the topology of the cerebral microcirculation (34, 41). The MCF filter eliminates the background as shown in Fig. 4a. The other filters erroneously characterize the sulcus as blood vessels. This artifact is especially severe in the Shikata filter. The Frangi filter suppresses the background, but fails to amplify vessel intensity. The MCF filter correctly identifies all arterioles, and eliminates the backdrop as desired. The ROC curve in Fig. 4b also indicates that the MCF has the best statistical performance with an area under the curve of 0.81. The MCF excels over all other filter by detecting more or equal vessels at any false positive level, making the MCF pareto-optimal. The dice coefficient similarity analysis of the MCF gives 0.73, Frangi with 0.68, Erdt with 0.64, Shikata with 0.62 and Sato with 0.56.

Figure 4
Visual and quantitative comparison of filter performance in a 3D network phantom study. (Panel a. top-left). The phantom consists of a microvascular tree embedded in an oval object mimicking a cerebral sulcus. A. Sato filter response. B. Frangi filter ...

We further measured how many of the 114 bifurcations in the microcirculatory phantom were correctly detected by each filter. Using the same criteria introduced above, the number of correctly detected bifurcation connections are given in Table 1 for each filter. All filters lose some bifurcations in the smallest diameter range (d~1 voxel). However, the MCF always surpasses other algorithms for all threshold levels without creating false positives.

The microcirculatory network phantom also featured branches with thin diameters between one to three voxels, a dimension close to the theoretical image feature threshold. Results of the statistical analysis for three diameter ranges (1-3 voxels) are summarized in Fig. 4c. It demonstrates that all filters tend to overestimate the size of the thinnest branches (d=1 voxel). For three voxels thick vessels (d=3 voxels), reconstruction with Sato filter gives a mean value of 2.5 voxels with standard deviation of 0.5 voxels, the Frangi filter 2.2±0.2, Shikata 3.1±0.3, and Erdt 3.3±0.5. The MCF gave 3.0±0.1, which is again the closest estimate. The bar graph in Fig. 4c also indicates that diameter reconstructions with Sato and Erdt filters have very large variance. The Frangi underestimates vessel diameters in the two and three voxel range. The diameter reconstruction of the Shikata response was limited to vessels outside the sulcus, since the filter wrongly recognizes the sulcus as a blood vessel, preventing diameter reconstruction of vessels in the interior of the domain. We also found that MCF overestimates vessels of one voxel in diameter by 10%, which is inevitable when using Gaussian smoothing.

Surface distance error and positional offset

Surface distance errors and positional offset were further analyzed in the 3D microcirculation phantom. The phantom serves as an ideal reference, because the true diameter and coordinates of the original object were exactly known and could be compared objectively to filtered artifacts. We rigorously computed for each vascular segment the surface distance error (SDE) between the original network coordinates and the coordinates of the MCF segmented model. Fig. 4a also shows a color coded SDE map, according to which, surface distance error with MCF is less than 0.5% everywhere, which is an acceptable deviation. Fig. 4d provides the statistical analysis of the SDE, indicating that relative surface distance errors are largest in smallest vessels. Since one voxel thick vessels are equal to the imaging threshold, this error cannot be eliminated entirely. We also computed the probability density function as depicted in Fig. 4d.

Moreover, results of the analysis of positional offset due to imaging, filtering and segmentation are summarized in a color coded map in Fig. 4a. Here, the exact position of the centerlines in the ground truth network was subtracted from the locations of the segmented network. The norm of the difference was plotted for all segments, indicating that the largest offset occurs also at smaller vessels as expected. The positional error is always less than 10μm as shown in Fig. 4d.

These two results demonstrate that the MCF and segmentation do not induce significant position or diameter errors. Deviations from the ground truth at the resolution level are minor and unavoidable.

Subject-Specific Computational Meshes from MRA and MRV

We also quantified filter performance using subject-specific image data in a sample group of six subjects. The workflow of Fig. 1 illustrates steps for generating three dimensional computational meshes from filtered images of six volunteers. The networks reconstructed with prior filters encompass only the COW and the initial segments of the ACA, MCA, and PCA (data not shown). The ACA territory is poorly reconstructed in all previous filters. By contrast, the MCF recognizes the arterial tree including a large portion of the pial network while maintaining physiological network connectivity without dangling segments. Fig. 5 shows the raw image data and the vessel enhanced intensity map, as well as reconstructed arterial and venous trees for six subjects. The reconstructed trees clearly delineate the subject-specific vascular topology. The cortical surfaces shown in Fig. 5 were reconstructed using the method introduced by Dale et al (32).

Figure 5
Reconstruction of the arterial and venous angioarchitecture from six subjects. First column. Superposition of raw MRA and MRV images. Second column. Superposition of filtered MRA and MRV images showing both the arterial and venous signal enhancement. ...

Quantification of filter performance with in vivo image data

Number of segments, length, and volume of the segmented vascular trees were used as metrics to assess filter performance. The size of the reconstructed vasculature was estimated by computing the number of segments, total length of segments and the volume enclosed by the surface mesh. Larger and longer trees indicate that a larger portion of the subjects' vascular network was reconstructed by the filter. All statistics are listed in Table 2, including total number of segments, arterial length, venous length, arterial volume, and venous volume.

Table 2
Total segment numbers, length and volume of the segmented arterial and venous tree.

On average, the MCF detects arterial trees with more than 6000 segments and a total length exceeding five meters. Statistics for each of the six individuals are summarized in the bar charts of Fig. 6 and the tabulation of Table 2. The number of segments and length of the arterial networks identified with MCF are at least three times more complete than the best of the prior filters. A similar advantage of three times the number of segments and three times of length is achieved in the venous trees. Moreover, the novel MCF successfully segmented on average three times larger arterial volumes with 6.09 mL. By comparison, Frangi filter only recovered a tree with 1.59 mL, Shikata filter with 1.32 mL, Sato filter with 1.20 mL, and Erdt filter with 0.62 mL. We also measured the volumes of the venous trees. The MRV captured major branches of the cerebral venous system, including the superficial and deep structures. The novel MCF detects 20.28 mL of the venous volume for subject 1. The Frangi filter gives a venous tree with only 4.83 mL, Erdt filter with 3.57 mL, Sato filter with 4.05 mL, and Shikata filter with 3.52 mL. The MCF typically achieved three times higher volume recognition rate compared to other filters.

Figure 6
Quantification of recognition rate of arterial and venous trees in numbers of segments, length, and total volume for six human subjects. The top bar graphs show the number of segments in the arterial and venous trees. MCF gives on average 6000 segments ...

These metrics show that MCF surpasses all prior vesselness filters. The MCF also preserves network connectivity which greatly simplifies the labeling of vascular territories. Color coded ACA, MCA, and PCA territories for the six subjects are depicted in Fig. 5. The labeling shown here was done manually, but can readily be automated. The statistics of volume, surface area, mean diameter and segment number for each territory are summarized in Table 3. The data also show that the MCF diameter reconstruction of main arteries and pial vessels agree with values reported in the literature (35-37). The agreement supports the notion that the MCF segmentation procedure does not induce diameter artifacts.

Table 3
Statistics of vascular territories for model construction of six subjects.

Discussion

Automatic vessel segmentation aims at reconstructing the entire cerebral vascular tree from MR images by enhancing blood vessel contrast across multiple length scales. Prior vesselness filters often fail to detect small tortuous vessels close to bifurcations resulting in incomplete segmentation of vascular trees. The novel MCF preserves the connectivity by combining the advantages of two separate contrast enhancement functions C1 and C2. C1 identifies tubular objects, while C2 magnifies vessel contrast without creating dangling segments. The multi-scale composite filter optimizes the contrast of the cerebral angioarchitecture across all length scales. Normalization along with thresholded levelsets in the fast marching achieved reasonably accurate diameter reconstructions.

The computational time for the MCF is roughly twice than that of the previous filters with total CPU time under 99s. The computational time could be further accelerated with parallel processing and GPU computation (31). Filter performance in terms of contrast enhancement, connectivity, surface distance error and positional offset was quantified in phantom studies. In all cases, MCF achieved the best score outperforming all prior algorithms.

We also performed expert manual delineation of the Circle of Willis. The ground truth was established by a neurosurgeon with more than 20 years of experience. The manually traced trajectory of major arteries in the raw image was compared to automatically reconstructed centerlines of the MCF image. No significant positional offset could be detected between manual and automatic segmentation (data not shown). However, the surgeon failed to identify significant structures such as the left and right posterior communicating arteries (PCOM) due to extensive blur in the raw image. This can be seen in Fig. 2. On the other hand, the MCF image clearly depicted both PCOMs as shown with arrows in Fig 2E. A similar limitation of the expert manual delineation was observed when trying to trace pial arteries (data not shown). Weak contrast impeded manual delineation of pial vessels, which only came to light in the MCF enhanced image. Given the inability even to trace prominent arteries in the COW by hand, we chose not to proceed with further manual delineation as a test for filter performance. Moreover, the large number of more than 6000 segments in each subject poses a practical limit on manual segmentation by expert surgeons or radiologists.

We successfully segmented subject-specific arterial and venous trees from in vivo data acquired in healthy human volunteers. The novel filter segmented a three times larger portion of the arterial and venous trees than any previous filter. The MCF enables the fully automatic reconstruction of the cerebrovascular angioarchitecture without requiring user supervision or input. The reconstructed subject-specific computational meshes can be used for CFD simulations of blood flow in the entire arterial or venous tree (39).

Limitations

The quality of the reconstructed vessel networks depend on the resolution of the acquired images. In the imaging study, the cerebral microcirculation could not be visualized since it is well below the image resolution of current MR scanners. The reconstruction of the venous network also suffers from the limited image quality because of slower venous flow. Because a 2D INHANCE sequence was used to acquire MRV, inter-slice discontinuity may diminish the reformatted 3D image quality. This limitation could be eliminated by using partially overlapping slices (i.e., negative slice gap) at the expense of increased number of slices to cover the brain, or by a 3D pulse sequence at the cost of longer acquisition time. The original MRA image contained some signals coming from the veins; these were eliminated by subtracting segments from the arterial image that are also registered in the MRV.

Conclusion

A novel multi-scale composite filter utilizes the eigenvalues of the Hessian matrix to achieve blood vessel enhancement over all relevant length scales. Our results confirm the feasibility of generating subject-specific computational meshes, with sufficient detail to track the ACA, MCA, and PCA territories in individual subjects almost down to the level of the pial network. Enhanced images provide more complete information about the pial blood supply for physicians to assess the status of the patient's cerebrovascular angioarchitecture.

Moreover, our filtering process enhances both large and small vessels without compromising network connectivity. The MCF allows the fully automatic reconstruction of major portions of the cerebrovascular network for the generation of high quality subject-specific computational meshes. The proposed automatic image segmentation could also be useful for computer-aided reconstruction of vascular trees to extract biometrics associated with pathological conditions in large data sets without variations caused by operator subjectivity (38). The novel filter will enable future studies on hemodynamic simulations of blood flow through the subject-specific angioarchitecture. This technology will be useful in predicting the outcome of planned intervention such as bypass surgery, stenting, or coiling of aneurysms.

Supplementary Material

Supp Material & Fig S1-S2 & Table S1-S2

Supporting Figure S1. Gaussian convolution eliminates noise in the signals for smoother representation. Intensity is extracted from a vessel in the original image along the red line with bar representation. Gaussian convolution with three different σ is applied with results shown on the right column. The convolution provides the most accurate smoothing as the σ approaches the real radius of the vessel. Smaller σ tends to suppress larger arteries while larger σ might overestimate the vessel diameter or merge with adjacent vessels.

Supporting Figure S2. Multi-scale MRA map showing the σ-value out of the arterial range (5 mm, 4 mm, 3 mm, 2 mm, 1 mm) that gave optimal vessel enhancement. Areas where σ=5mm are drawn in red. Note that these areas correspond to the large arteries in the CoW. Areas where σ=3mm is optimal, shown in green, map the initial segments of the MCA, and ACA as well as larger pial arteries. σ=1mm in blue was the optimal setting for the smallest pial arteries in the diameter range of 400-1000μm. Black indicate areas where the intensity is zero, therefore no σ value is chosen.

Supporting Table S1. Eigenvalue behaviors due to shapes of the feature. L indicates low value, H indicates high value, + indicates positive value, and – indicates negative value.

Supporting Table S2. This summarizes the number of adjustable parameters for each filter. It is important to note that the eigenvalue ordering differs between methods as indicated in row entitled eigenvalue ordering.

Acknowledgments

The authors would like to gratefully acknowledge partial support of this project by NSF grants CBET-0756154 and CBET-1301198. X.J.Z. acknowledges support by the following NIH grants: 1S10RR028898 and UL1RR029879. No conflicts of interest were posed in the conduct of this research.

References

1. Go AS, Mozaffarian D, Roger VL, et al. Heart disease and stroke statistics--2013 update: a report from the American Heart Association. Circulation. 2013;127:e6–e245. doi: 10.1161/CIR.0b013e31828124ad. [PMC free article] [PubMed] [Cross Ref]
2. Neumann-Haefelin T, Moseley ME, Albers GW. New magnetic resonance imaging methods for cerebrovascular disease: Emerging clinical applications. Ann Neurol. 2000;47:559–570. [PubMed]
3. Kodoma T, Suzuki Y, Yano T, Watanabe K, Ueda T, Asada K. Phase-contrast MRA in the evaluation of EC-IC bypass patency. Clin Radiol. 1995;50:459–465. [PubMed]
4. Taylor CA, Draney MT, Ku JP, Parker D, Steele BN, Wang K, Zarins CK. Predictive medicine: computational techniques in therapeutic decision-making. Comput Aided Surg Off J Int Soc Comput Aided Surg. 1999;4:231–247. [PubMed]
5. Steinman DA, Milner JS, Norley CJ, Lownie SP, Holdsworth DW. Image-Based Computational Simulation of Flow Dynamics in a Giant Intracranial Aneurysm. Am J Neuroradiol. 2003;24:559–566. [PubMed]
6. Hademenos GJ, Massoud TF. Risk of Intracranial Arteriovenous Malformation Rupture Due to Venous Drainage Impairment A Theoretical Analysis. Stroke. 1996;27:1072–1083. [PubMed]
7. Cebral JR, Yim PJ, Löhner R, Soto O, Choyke PL. Blood flow modeling in carotid arteries with computational fluid dynamics and MR imaging. Acad Radiol. 2002;9:1286–1299. [PubMed]
8. Milner JS, Moore JA, Rutt BK, Steinman DA. Hemodynamics of human carotid artery bifurcations: Computational studies with models reconstructed from magnetic resonance imaging of normal subjects. J Vasc Surg. 1998;28:143–156. [PubMed]
9. Zhao SZ, Xu XY, Hughes AD, Thom SA, Stanton AV, Ariff B, Long Q. Blood flow and vessel mechanics in a physiologically realistic model of a human carotid arterial bifurcation. J Biomech. 2000;33:975–984. [PubMed]
10. Stroud JS, Berger SA, Saloner D. Numerical analysis of flow through a severely stenotic carotid artery bifurcation. J Biomech Eng. 2002;124:9–20. [PubMed]
11. Grinberg L, Anor T, Cheever E, Madsen JR, Karniadakis GE. Simulation of the human intracranial arterial tree. Philos Trans R Soc Lond Math Phys Eng Sci. 2009;367:2371–2386. [PubMed]
12. Grinberg L, Anor T, Madsen J, Yakhot A, Karniadakis G. Large-Scale Simulation of the Human Arterial Tree. Clin Exp Pharmacol Physiol. 2009;36:194–205. [PubMed]
13. Zagzoule M, Marc-Vergnes JP. A global mathematical model of the cerebral circulation in man. J Biomech. 1986;19:1015–1022. [PubMed]
14. Sato Y, Nakajima S, Shiraga N, Atsumi H, Yoshida S, Koller T, Gerig G, Kikinis R. Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images. Med Image Anal. 1998;2:143–168. [PubMed]
15. Frangi AF, Niessen WJ, Vincken KL, Viergever MA. Multiscale vessel enhancement filtering. In: Wells WM, Colchester A, Delp S, editors. Medical Image Computing and Computer-Assisted Interventation — MICCAI'98. Lecture Notes in Computer Science; Springer; Berlin Heidelberg: 1998. pp. 130–137.
16. Frangi AF, Niessen WJ, Hoogeveen RM, van Walsum T, Viergever MA. Model-based quantitation of 3-D magnetic resonance angiographic images. IEEE Trans Med Imaging. 1999;18:946–956. [PubMed]
17. Erdt M, Raspe M, Suehling M. Automatic Hepatic Vessel Segmentation Using Graphics Hardware. In: Dohi T, Sakuma I, Liao H, editors. Medical Imaging and Augmented Reality. Springer; Berlin Heidelberg: 2008. pp. 403–412. Lecture Notes in Computer Science.
18. Kaftan JN, Kiraly AP, Bakai A, Das M, Novak CL, Aach T. Fuzzy pulmonary vessel segmentation in contrast enhanced CT data. 2008;6914:69141Q–69141Q–12.
19. Shikata H, McLennan G, Hoffman EA, Sonka M. Segmentation of Pulmonary Vascular Trees from Thoracic 3D CT Images. Int J Biomed Imaging. 2009;2009:e636240. [PMC free article] [PubMed]
20. Shikata H, Hoffman EA, Sonka M. Automated segmentation of pulmonary vascular tree from 3D CT images. 2004;5369:107–116.
21. Drechsler K, Laura CO. Comparison of vesselness functions for multiscale analysis of the liver vasculature. 2010 10th IEEE International Conference on Information Technology and Applications in Biomedicine (ITAB) 2010:1–5.
22. Sethian JA. A fast marching level set method for monotonically advancing fronts. Proc Natl Acad Sci. 1996;93:1591–1595. [PubMed]
23. Caselles V, Kimmel R, Sapiro G. Geodesic Active Contours. Int J Comput Vis. 1997;22:61–79.
24. Antiga L, Piccinelli M, Botti L, Ene-Iordache B, Remuzzi A, Steinman DA. An image-based modeling framework for patient-specific computational hemodynamics. Med Biol Eng Comput. 2008;46:1097–1112. [PubMed]
25. Fang Q, Boas DA. Tetrahedral mesh generation from volumetric binary and grayscale images. IEEE International Symposium on Biomedical Imaging: From Nano to Macro, 2009. 2009:1142–1145. ISBI '09.
26. Houman Borouchaki LP. Parametric surface meshing using a combined advancing-front generalized Delaunay approach. Int J Numer Methods Eng - Int J Numer Method Eng. 2000;49:233–259.
27. Nowinski WL, Chua BC, Marchenko Y, Puspitsari F, Volkau I, Knopp MV. Three-dimensional reference and stereotactic atlas of human cerebrovasculature from 7 Tesla. NeuroImage. 2011;55:986–998. [PubMed]
28. Ho J, Kleiven S. Can sulci protect the brain from traumatic injury? J Biomech. 2009;42:2074–2080. [PubMed]
29. Linninger AA, Tsakiris C, Zhu DC, Xenos M, Roycewicz P, Danziger Z, Penn R. Pulsatile cerebrospinal fluid dynamics in the human brain. IEEE Trans Biomed Eng. 2005;52:557–565. [PubMed]
30. Linninger AA, Xenos M, Zhu DC, Somayaji MR, Kondapalli S, Penn RD. Cerebrospinal fluid flow in the normal and hydrocephalic human brain. IEEE Trans Biomed Eng. 2007;54:291–302. [PubMed]
31. Asano S, Maruyama T, Yamaguchi Y. Performance comparison of FPGA, GPU and CPU in image processing. International Conference on Field Programmable Logic and Applications, 2009. 2009;FPL 2009:126–131.
32. Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis. I. Segmentation and surface reconstruction. NeuroImage. 1999;9:179–194. [PubMed]
33. Ghaffari M, Hsu CY, Linninger AA. Automatic Reconstruction and Generation of Structured Hexahedral Mesh for Non-planar Bifurcations in Vascular Networks. PSE2015 ESCAPE25. 75
34. Linninger AA, Gould IG, Marinnan T, Hsu CY, Chojecki M, Alaraj A. Cerebral microcirculation and oxygen tension in the human secondary cortex. Ann Biomed Eng. 2013;41:2264–2284. [PMC free article] [PubMed]
35. Kamath S. Observations on the length and diameter of vessels forming the circle of Willis. J Anat. 1981;133:419–423. [PubMed]
36. Bevan JA, Dodge J, Walters CL, Wellman T, Bevan RD. As human pial arteries get smaller, their wall thickness and capacity to develop tension relative to their diameter increase. Life Sci. 1999;65:1153–1161. [PubMed]
37. Jain KK. Some observations on the the anatomy of the middle cerebral artery. Can J Surg J Can Chir. 1964;7:134–139. [PubMed]
38. Doi K. Computer-aided diagnosis in medical imaging: Historical review, current status and future potential. Comput Med Imaging Graph. 2007;31:198–211. [PMC free article] [PubMed]
39. Hsu CY, Schneller B, Ghaffari M, Alaraj A, Linninger A. Medical Image Processing for Fully Integrated Subject Specific Whole Brain Mesh Generation. Technologies. 2015;3:126–141.
40. Lorensen WE, Cline HE. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques; SIGGRAPH '87; ACM: New York, NY, USA. 1987. pp. 163–169.
41. Karch R, Neumann F, Neumann M, Schreiner W. Staged growth of optimized arterial model trees. Annals of Biomedical Engineering. 2000;28:495–511. [PubMed]