Virtual Contrast for Coronary Vessels Based on Level Set Generated Subvoxel Accurate Centerlines

We present a tool for tracking coronary vessels in MRI scans of the human heart to aid in the screening of heart diseases. The vessels are identified through a single click inside each vessel present in a standard orthogonal view. The vessel identification results from a series of computational steps including eigenvalue analysis of the Hessian of the MRI image followed by a level set-based extraction of the vessel centerline. All identified vessels are highlighted using a virtual contrast agent and displayed simultaneously in a spherical curved reformation view. In cases of over segmentation, the vessel traces can be shortened by a click on each vessel end point. Intermediate analysis results of the vessel computation steps can be displayed as well. We successfully validated the tool on 40 MRI scans demonstrating accuracy and significant time savings over manual vessel tracing.


INTRODUCTION
Coronary artery disease (CAD) is one of the leading causes of mortality and morbidity in the USA and other industrialized nations [1]. Although conventional cardiac angiography remains the "gold standard" for the evaluation of CAD, it is an invasive procedure that is associated with morbidity (1.5%) and mortality (0.15%) risks [2]. Coronary CT angiography (CTA) is emerging as a promising noninvasive alternative [3]. However, this technology requires patient exposure to substantial amounts of radiation [4] and potentially nephrotoxic contrast agents as in conventional angiography. As a result, coronary magnetic resonance angiography (CMRA) provides a more patient friendly option for CAD assessment [5] without the use of contrast agents and radiation. Unfortunately, currently the MRI image signal-to-noise ratios and the maximally achievable resolution are not as high as for CTA. This complicates the process of identifying the vessels and MRA targeted vessel segmentation and analysis tools usually fail. Thus, the common solution is still time consuming, manual vessel tracing. Our research presents an MRI coronary identification software tool that is able to track and intuitively display the MRI data along all three main coronary vessels. A typical screenshot of the software is presented in Figure 1.
A core component of the system is the computation of a centerline for each vessel. As the vessels are only a few voxels thick, it is important to compute the vessel centerlines at subvoxel accuracy. They must also be inherently smooth to yield a smooth display of the vessel in the spherical curved reformation view.
The next section describes the technical background of our methods followed by related work and a description of the novel methods used in our system. We conclude with showing and discussing results acquired with our software.

BACKGROUND
Many automatic and semiautomatic skeletonization techniques compute the centerline of an object on the voxel grid with optional subsequent smoothing [6][7][8][9][10][11][12][13][14][15][16]. These discrete centerline solutions are inappropriate for vessels that are only a few voxels thick. Subsequent smoothing may displace the centerline from the vessel center and is thus inappropriate as well. Another method [17] computes the centerline as a minimum cost B-spline. This delivers an inherently smooth centerline, but is computationally expensive due to the explicit global optimization. The same holds for [18], which iteratively computes a globally optimized NURBS curve that locally minimizes the vessel cross sections perpendicular to 2 International Journal of Biomedical Imaging the curve. In contrast, our prior centerline algorithm efficiently and automatically extracts a smooth, continuous centerline directly at subvoxel precision [19]. Our algorithm is based on level sets. Level set methods evolve an isosurface in the direction of the surface normal [20]. In its general form the evolution speed can depend on position, normal direction, curvature, and shape; and the isosurface can cross over the same point multiple times. In our centerline method the evolution speed is always positive and depends only on position. Hence, its boundary front moves only outwards. With these restrictions the isosurface can be represented by an Eikonal equation: where T is the arrival time function, F is the speed of evolution function, and Γ is the initial isosurface at time zero. An efficient method to numerically evaluate the solution to the Eikonal equation is the fast marching method [20]. It processes the voxels in a sorted order based on increasing values of T. The fast marching method calculates a time crossing map, which indicates for each pixel how much time it would take for the level set front to arrive at the pixel location. The evolution only needs to be computed on a rectilinear grid. However, values at nongrid locations can be interpolated from these grid positions properly to simulate the true propagation value.
Several other centerline methods based on level sets have been previously presented [10][11][12]. One approach [10] computes centerlines by first detecting medial axis points at the locations where the level set fronts collide and form sharp discontinuities. However, the level sets are only computed on two dimensional cross-sections of the three dimensional data, which are not identical to the 3D discontinuities. Next, the algorithm performs topological thinning and filling in of gaps with voxels along straight lines which may not result in positions on the skeleton. Other methods [11,12] make use of the full 3D data in its level set propagation and guarantee a minimum cost, discrete solution, but as pointed out before, do not extract the centerline with subvoxel accuracy. In addition, algorithms [6-10] require a segmentation of the vessel as input. In our images it is very difficult to segment the vessels accurately and completely. Hence, we were looking for an algorithm that does not require a segmented vessel as input. An algorithm that is subvoxel accurate and does not require vessel segmentation is presented in [21]. It directly traverses the centerline along ridges in a Hessian medialness measure, but it performs a sequence of local greedy decisions that do not guarantee a globally optimal solution. The methods in [11][12][13][14][15][16][17]21] use Hessian matrix analysis to track a vessel with having to segment its boundaries first. Computing the Hessian at different scales proved to be beneficial for vessels varying greatly in thickness, but was not necessary for our data. Table 1 lists the prior methods and classifies them according to the main algorithm ideas. In [22], 94 vessel extraction algorithms are compared: only 50% of them do compute an explicit centerline, only one uses level set methods, but not for centerline tracking and only one uses Hessian eigenvalues, but not in combination with level sets. Our proposed algorithm combines the benefits of the previous methods without their shortcomings. It can find minimum cost, subvoxel accurate centerlines of thin vessels without the need of first segmenting them.

Subvoxel accurate centerline algorithm
The most closely related prior algorithm is our input segmentation dependent, subvoxel accurate centerline algorithm [19]. It uses a level set segmentation of the vessel to obtain a subvoxel accurate surface and a Euclidean distance transform of the object. This distance transform is then used as a speed image in a fast marching level set method with propagation starting at the global maximum point of the distance transform. The fast marching method propagation is augmented to calculate the geodesic distance in addition to the time crossing map. The furthest geodesic point from the global maximum point is used as the start point of the vessel centerline and the remaining points of the centerline are determined by performing a gradient descent on the time crossing map with a subvoxel step size.
The algorithm presented in this paper is an extension of this previous method, which handles the absence of vessel segmentation and improves upon the computation speed of the level set propagation.

Vessel enhancement
In order to track thin vessels without an explicit representation, we found it necessary to process the MRI images using vessel enhancing image filters. Given the eigenvalues λ 3 ≤ λ 2 ≤ λ 1 of the 3 × 3 Hessian matrix for each 3D image pixel, it is possible to compute a likelihood of the pixel being part of a linear structure [23,24]. This measure, which we Table 1: Comparison of ideas used in various vessel centerline computation methods. "+" stands for the best idea within a group, "0" for average, and "−" for the least desirable idea in a group.

Curved reformation vessel view
A good vessel centerline can be used to create a curved reformation vessel view. One such method is to stretch the vessel and display its surroundings with as little distortion as possible [26]; however this is not appropriate for a view that is supposed to include multiple vessels. The "soap bubble" method [27] allows projection of multiple vessels to a plane, preserving the relationship of the vessels to each other, but for roughly spherically arranged vessels, the projected vessels may overlap or surrounding tissue may be severely distorted. We use the spherical curved reformation method [28], which eliminates the problems of the "soap bubble" method for the specific case of coronary vessels. It achieves this through projecting a spherical approximation of the heart onto the vessels with little distortion, followed by a standard globe unrolling onto a rectangular view as done for any world map. The minimal distortion is the consequence of minimizing the energy of a "thin plate spline" being fit to the vessel points [28].

METHODS
While the system is designed to look at all vessels at once, processing is done one vessel at a time. The user identifies each vessel initially by clicking on one landmark point for each vessel.

Vessel centerline computation
The vessel centerline computation results from a series of computational steps ( Figure 2). Initially, noise removal is performed on the MRI data by using edge preserving anisotropic diffusion filtering. Next, the intensities are normalized through a sigmoid window whose parameters are  determined from the MRI intensity I at the first landmark point (A in Figure 3) which must be inside the vessel and is assumed to be the maximum intensity in the local neighborhood. The width of the sigmoid window is equal to I and the center is at I/2. Subsequently the vesselness, υ, is computed for each pixel in the image by (2) from the Hessian matrix of the image. The partial derivatives that form the Hessian matrix are a result of convolving the smoothed MRI image with the derivatives of a Gaussian with 3σ covering 2 mm, which is the median of the expected vessel diameter. This vesselness map is then normalized using a sigmoid window. Again the window parameters are relative to the landmark point (width =õ, center =õ/2), which maximizes the contrast in the transition region. Values in the normalized vesselness map are clamped to zero if they are less than 33% of the maximum value. The clamp threshold and σ were determined empirically for a single dataset and applied for all others.
This resulting vesselness map (middle image in Figure 4) is used as a speed function for a fast marching level set method that starts at the initial vessel landmark point A. However, instead of computing the fast marching through the complete image or at least until the entire vessel is covered (as in [11]), the computations are stopped when a point 1 cm from the landmark point is reached. Due to the nature of the vesselness computation, the highest speed values are found in the center of the vessel, and thus the point first reached at 1 cm distance (larger than the vessel radius) must be central to the vessel. This point is then defined as the "spearhead point." The fast marching method is now continued, but only newly discovered points that are within 1 cm of the "spearhead point" are allowed to be added to the evolving surface of the fast marching method. Each time a new furthest geodesic distance point is found, the "spearhead point" is updated. Consequently, only a small band of voxels along the vessel is involved in the computation. When the modified fast marching method has processed all of the points in the connected object, the final "spearhead point" is the most distant trackable vessel point. This becomes the second, automatic landmark point (B in Figure 3). The steepest gradient decent from the second to the first landmark point yields a partial vessel centerline. This centerline is not based on local greedy decisions, but is the minimum cost path with respect to the given vesselness speed image. Next, the speed image is intensified along this partial vessel. With the updated speed image, the above described fast marching method is started from the second automatic landmark. Again, only a small band of voxels along the vessel are involved in the computation. When the second fast marching method algorithm completes, the resulting "spearhead point" becomes the third, automatic landmark point (C in Figure 3). Intensifying the speed image along the first partial centerline is necessary to guarantee that the initial path of the second fast marching does not detour into a vessel branch. Finally, the steepest gradient decent from the third (C) to the second (B) landmark point yields the complete vessel centerline ( Figure 5).
Once the centerline is computed, a virtual contrast dataset can be created. The virtual contrast dataset has the original data as its basis, but each pixel is intensified that is within 2 mm (expected vessel radius) of the computed centerline and also has at least 33% of the maximal vesselness intensity in the initial speed image.
The above process can be repeated for each of the desired vessels. After the centerline for each vessel has been computed, a new spherical curved reformation can be generated from all current vessel centerlines. The length of each of the vessels is displayed on the graphical user interface. In order to be able to bridge areas of stenosis or low signal, short (< 1 cm) sections of low vesselness (below 33% of the maximum intensity value) are allowed, as long as the tracking can be continued with more obvious vesselness pixels after the gap. Unfortunately, this sometimes also causes the tracking to go beyond the vessel ends into other nearby vessels or to jump onto the edge of the heart, which may not be totally suppressed in the speed image. In this case, the curved reformation view can be used  to allow the user to manually relocate the second and/or third landmark point to the desired vessel endpoint(s) (D and E in Figure 2). The system then adjusts the centerline to only cover the vessel between these updated vessel end points. System validation was based on visual assessment of completeness of the tracked vessels, partial success was defined as a section of the vessel being visible in the scan, but not tracked. The scanning protocols used for the validation were (A) standard imaging parameter, (B) shortened acquisition window, (C) isotropic acquisition voxel resolution, (D) short-axis plane aligned with the right coronary artery.

RESULTS
The algorithm presented was validated on 40 MRI cardiac scans with volume sizes of 512× 512 × 100 to 512 × 512 × 300 containing 0.5 × 0.5 mm images with 0.5 to 1 mm spacing in the z-direction. The data came from 10 volunteers scanned each with four scanning protocols. The right coronary artery (RCA) was found completely for 90% of the volunteers on two protocols (A, B), for the other protocols 80% (B) and 70% (D) of the MRI scans had completely tracked vessels. For the incomplete scans it was possible to complete them by treating the missed vessel section as a new vessel. Figure 3 shows some intermediate and final results.
Completing the interactive part of the vessel tracking was accomplished within one minute. On a professional medical image analysis workstation a trained radiologist took 2.5 minutes to hand segment a vessel via coarse contours on every third slice that were then interpolated by the system.

DISCUSSION
The results in Figures 1, 3, and 5 show spherical curved reformations of the three main coronary arteries. In this example for two arteries joined near the aorta a single landmark was sufficient, but for the third both endpoints needed to be corrected. In either case, the user interaction time required is minimal when compared to manual vessel tracking.
All numerical parameters listed in this algorithm degrade gracefully. A 10% change of the parameter value has only a small impact on the final result, but doubling or halving it usually significantly shortens the identification of vessel segment.
The novelty of this research is two fold. First, it lies in the creation of a time saving tool that combines the idea of semi-automatic tracking with spherical curved reformation. Second, it improves over prior work on the method of finding the vessel centerline. Due to the low signal-tonoise ratio on the MRI input images, the vesselness map is a network of mutually connected vessels and pseudovessels. Following all branches as done in [10] frequently results in automatically found vessel end landmarks that are very far from the intended vessel end. The restriction of expanding the fast marching only near the "spearhead point" allows for a much more intuitive behavior of the algorithm. The manual clipping of the traced path to only the portion within the vessels is easily performed on the spherical curved reformation view and no scrolling through slices is needed.

CONCLUSIONS
This paper has presented a semi-automatic algorithm for determining centerlines of the main coronary vessels and its application to creating virtual contrast enhanced MRI scans that are displayed in an intuitive spherical curved reformation view. The method can track vessels even in the presence of low signal-to-noise ratios, is subvoxel accurate, and is more computationally efficient than previous methods.