Image Processing for Arterial Tree Morphometry

Early studies of vascular structures relied solely on qualitative visual interpretation of the acquired imagery. Severe stenoses and vascular abnormalities such as arteriovenous shunts could be observed and documented, for example as an aid to surgical planning. Interpretation of histological structure, electron micrographic ultrastructure and early "3D" methods such as scanning microscopy of freeze-dried tissue or corrosion casts was also purely qualitative. However, a great deal of physiological knowledge was gained, in particular about the cellular constituents and ultrastructure of the vascular wall, the connectedness of the microvascular network and the spatial relationships, for example, between alveolar sacs and the flat capillary sheets that bathe them in blood separated from the air by a two-cell-thick endothelial-epithelial wall [73]. The state of early research in the field is well summarized by a collection of papers covering all aspects of small vessel angiography from the effects of focal spot size, X-ray spectrum, magnification, film-screen types and combinations, and electronic imaging chain imperfections to clinically relevant applications in the cerebral, pulmonary, renal, pancreatic, and coronary circulations [74]. However, for clinical methods to be reliably applied in routine practice, the need for quantitative measures of vascular morphometry was clearly evident. In the basic science arena, it is equally clear that quantitative statistical methods are required to distinguish significant differences between normal and diseased structures, and between diseased structures of varying pathology or pathogenesis.

An important early application of quantitative methods was in the area of assessing the degree of stenosis, particularly of coronary vessels [75]. Given a gray-level image of a vascular tree such as those available from radiography, fluoroscopy, or MRA, loosely called "angiograms" hereafter, there are two basic approaches for extraction of quantitative metrics. (Either may involve judicious preprocessing of the gray-level image [76-78].) The first is to create a binary image from which measurements can be made in a straightforward manner, and the other is to operate on the gray-level image directly. To make a binary image is to decide, based upon some criterion or set of criteria, which pixels in the image belong to the vessels, usually the intraluminal space filled with contrast agent, and which do not. The vessel pixels are then turned on and all the other pixels turned off, creating the two-phase or binary image. The process of creating a binary image in this case is called segmentation of the arterial tree. Possible criteria for selecting vessel pixels range from thresholding procedures, wherein (for the case of a global, or simple threshold) every pixel below a specified gray level is called vessel and all pixels above the threshold background, to more complex methods involving, for example, vessel tracking with decisions based upon the location of the maximum gradient magnitude at the vessel edge. Adaptive threshold values may be based upon neighborhood gray-level statistics such as the local histogram or more complex parameters such as entropy [79]. Recently, rather sophisticated techniques for segmentation of MRA data have been developed. One utilizes assumed distributions for tissue-type gray levels and a modification of the expectation maximization (EM) algorithm to distinguish vessel from nonvessel voxels and has been applied to cerebral time-of-flight MRA data [70]. Viergever et al. developed a model-based approach wherein the vessel medial axis was modeled using a B-spline curve and the vessel wall with a deformable tensor product B-spline surface to obtain a segmentation from which diameter assessment could be performed with subvoxel precision [67].

Mathematical morphology-based methods including watershed analysis show promise for segmentation of complex structures [80-82]. In gray-level mathematical morphology, an image is treated as a topographical surface in which the gray level corresponds to elevation [83,84]. The watershed transform works by piercing holes in local minima in the terrain and flooding the surface from below [85]. Each minimum starts its own "catchment basin." As the flood rises, when the water from adjacent catchment basins is about to meet, dams are built to keep them from joining, pixel by pixel. When the process is complete, the water has risen to the level of the highest peak, and the transformed image consists of the crest lines of the dams separating catchment basins. In the case of a vascular image, these crest lines would ideally correspond to vessel medial axes. Vessel contours can be segmented by watershed transformation of the gradient image. Brute force computation of the watershed transform results in severe oversegmentation: Too many crest lines are extracted if all local minima are used. Instead, a subset of minima has to selected by use of a marker function whose catchment basins correspond to the objects to be segmented. The marker function is generally defined using a priori knowledge about the nature of the objects to be segmented. Region growing methods are also quite popular [86]. These involve specification of a pixel or pixels, called "seeds," that are definitely contained within the vessel lumen, then growing the regions around the pixels so designated either to include all candidates between an upper and a lower threshold bound ("hysteresis thresholding" [87]) or until a boundary satisfying some criterion, such as gradient magnitude, is reached. These boundaries are then the vessel edges. The analogy to over-segmentation by watershed transformation in the case of region growing would mean that portions of arterial tree remained disconnected. In this case, boundaries may be dissolved according to a chosen weakness criterion. One reason for the popularity of seeded region growing methods for segmenting structures such as arterial trees is that vascular networks are known to be connected: Every pixel in the lumen has to be connected to the arterial inlet. Incorporating a connectivity constraint into a segmentation algorithm can help to eliminate spurious background pixels or regions from being included as vessel.

A method for segmenting tubular objects from gray-level medical images has been developed by Aylward and Pizer et al. [88,89]. It is a variant of their multiscale medial analysis algorithms, such as marching cores [90,91] and exploits the facts that vessels are tubular and the medial axes of tubular objects imaged with imperfect imaging systems are well approximated by their intensity ridges. The vessel intensity ridges are tracked, and convolution with a Laplacian of a Gaussian is used to calculate the local maximum in multiscale medialness at each point along the ridge as an estimate ofvessel width. This method for vessel diameter estimation benefits from the smoothly varying, nearly circular shape of tubular cross-sections. The tubular object segmentation method makes it possible, without imposing severe computational demands, to extract representations that preserve the location, size, and topology of objects in 3D images with minimal user interaction, and is stable in the presence of noise. Figure 2a is a volume rendering of 48 MRA slices depicting the cerebral vasculature. The challenges for automated extraction of vessel boundaries, including variations in vessel intensity and width, and the complexity of the vascular structures in the noisy background are clearly evident. Figure 2b shows the result of the tubular segmentation algorithm. The surface rendering is based upon 105 ridges and their widths extracted from the data set of

Fig. 2a. The ventricular contours and brain surface are shown in green and red, respectively, for reference.

Creating a binary image is a required prerequisite for producing surface shaded renderings and can facilitate some feature extraction tasks such as obtaining the tree skeleton. Binary skeletonization, also called thinning, is a well-established technique, at least in two dimensions, whereas gray-level thinning is not [92]. Skeletonization, in turn, achieves dramatic data reduction and makes the specification of bifurcation locations (branch points) and branching angles more straightforward, particularly in 2D images. The skeleton or approximation of the medial axis for each segment can be fitted to a line or a curve and the three midlines of a parent and two daughter segments at a bifurcation forced to meet at a common point as shown in Fig. 1. This midline juncture then defines the branch point unambiguously. Skeletons properly extracted preserve the topology of the branching tree, as well as segment lengths and branching angles, though diameters, unless retained as separate items in a data structure linked either to medial axis points or segments, are lost. However, although segmentation followed by skeletonization may be appropriate in some circumstances, in general ( particularly where the accuracy of vessel diameter measurements is of paramount importance, as is so often the case) it may not be advantageous to segment the vascular image. Instead it is more exact and precise to operate on the gray-level image data directly. Gray-level image-based tree extraction schemes

(a)

capture the vessels' medial axes and borders. If the whole tree or large portions thereof are to be identified, they usually involve tracking the vessel centerlines, starting at the arterial inlet or other user specified point, and searching a rectangle or sector distal to the current centerline point, symmetric about the current centerline direction [93-96]. The dimensions of the sector are dynamically adjusted to take into account the local noise level [93,94] and vessel tortuosity, or curvature [97]: If either or both of these are high, the search-ahead distance is shortened in the interest of more reliable tracking. The direction of search is updated by use of a matched filter, for example a rect or triangular function, based on a model incorporating the current vessel width and the expected shape of the vessel's cross-sectional intensity profile. Branch points are detected by comparing the image brightness level outside the current borders to the brightness at the vessel center. In recursive tracking, once a major trunk and its branch points have been identified, they are deleted from the image to avoid repeated tracking. The new branch centerline directions at each bifurcation are then determined, again by the use of a matched filter, and smaller and smaller subtrees are tracked and deleted in succession. Several parameters, including a minimal acceptable vessel length, need to be adjusted depending on image quality to avoid capture of nonvessel structures or noise in the image.

Some of the most rigorous processing methods have been developed for computerized systems designed to improve the

(b)

FIGURE 2 (a) Avolume rendering of 48 MRA slices of the cerebral vasculature, (b) The result of the tubular segmentation algorithm. Images courtesy of Stephen Aylward, University of North Carolina at Chapel Hill.

accuracy and objectivity of clinical assessment of stenosis severity [75,98]. These methods are based on the classic work of Canny who developed a computational approach to edge detection for antisymmetric edges (f(x) = —/(—x), where the true edge point is at x = 0) such as ridge, roof, and step edges [99]. The requirements for reliable performance are good detection and good localization. Good detection means that a true edge is highly likely to be a detected and that the likelihood of detecting an edge where there is none is low. Good localization means that the error in the location of a detected edge relative to that of the true edge should be small. Canny developed mathematical forms for the criteria for low error rate and good localization and showed that there is a trade-off between the two that dictates a unique shape for the impulse response of the optimal detector. Further, the detection-localization trade-off varies with the spatial extent of the filter and the noise level in the image, mandating a multiscale approach followed by synthesis of the features detected by operators at the various scales. Avariable-scale approach makes sense in light of the fact that intensity changes occur over a range of scale in real images [100]. An approximate realization of the ideal filter for antisymmetric edges can be realized in the spatial domain by designating edge points as the maxima in the gradient magnitude of a Gaussian smoothed image, also called the Gaussian-weighted gradient operator. Intuitively, this filter is a derivative operator with a smoothing effect — in fact, the first derivative of a Gaussian — so edges are located at the maxima in the smoothed derivative of a line scan perpendicular to the edge. In fact, Marr and Hildreth had pointed out much earlier in the derivation of their Laplacian of Gaussian operator that the optimal filter for the image averaging at different scales had to be Gaussian, since this provides the best compromise for the conflicting requirements for band limit-edness (localization) in both spatial (because edge information is highly localized) and frequency (to reduce the range of scales over which intensity changes take place) domains [ 100]. In one dimension, the zero crossings [101] produced by the Marr-Hildreth operator correspond exactly to the maxima in Canny's Gaussian-weighted gradient operator, but Canny claimed two advantages to his approach: The directional properties of his operator enhance its detection and localization properties in 2D implementations; and the amplitude of the response at maxima provides a good estimation of edge strength, which can be exploited by a subsequent thresholding operation on the gradient image. An important contribution of Reiber's group was to extend Canny's computational approach to the analysis of asymmetric edges like those of density profiles through arterial cross-sections. The consequence of asymmetric edges is that two additional constraints on the operator must be added: the precise shape or phase shift of the filter has to be optimized for a particular arterial diameter and, when a filter of the same spatial extent is applied to arteries of different sizes, correction factors to the diameter estimate have to be invoked. A suite of angiographic image analysis software incorporating the optimal edge detection algorithms was developed at the Leiden University Medical Center to assist in clinical decision making regarding severity of coronary artery stenosis [102]. Figure 3 shows one display mode of the image and results panel. An enlarged view of the portion of the angiogram containing the lesion is shown, with the detected arterial borders and highlighted stenosis displayed as a graphical overlay on the image. A plot of artery diameter versus length in the diseased region and a panel of quantitative results are also displayed on the same screen. Phantom and clinical experiments showed that the method performed accurately and reproducibly for vessels between 1.2 and 6 mm in diameter even in the presence of realistic noise levels, and for lower noise levels performance remained acceptable for vessels with diameters greater than 0.7mm [98].

Focal stenoses can be readily appreciated and quantified with reasonable accuracy from single angiograms by the methods just described, but subtle disease may elude detection on single-projection angiograms. Three-dimensional reconstruction techniques may facilitate accurate assessment of coronary arterial disease from routinely acquired biplane angiograms, even in cases of complex morphology. For example, at the German Heart Institute of Berlin these methods have been successfully applied to the quantification of diffuse coronary artery diseases, conditions that may remain occult if only single

FIGURE 3 Portion of a clinical coronary angiogram processed with the software developed at the Leiden University Medical Center. The stenotic lesion is highlighted, a plot of vessel diameter versus length in the diseased region is displayed, along with a panel of the important quantitative results including percent stenosis. Image courtesy of Johan Reiber, Leiden University Medical Center. See also Plate 26.

FIGURE 3 Portion of a clinical coronary angiogram processed with the software developed at the Leiden University Medical Center. The stenotic lesion is highlighted, a plot of vessel diameter versus length in the diseased region is displayed, along with a panel of the important quantitative results including percent stenosis. Image courtesy of Johan Reiber, Leiden University Medical Center. See also Plate 26.

FIGURE 4 (a) Source biplane angiograms: the upper pair depicts the left coronary artery system (LCA) in time-equivalent projections from right and left anterior oblique (RAO/LAO) views; the lower pair shows the right coronary artery system of the same subject in a second acquisition. (b) Visualization of the 3D model reconstructed from the angiograms of Fig. 4a. The right coronary artery (RCA) is in the front, at the same level behind it to the left circumflex artery (LCX), and to the right the left anterior descending artery (LAD). Images courtesy of Andreas Wahle, University of Iowa.

angiograms are evaluated in two dimensions [103-104]. Figure 4a shows four angiographic views of a normal subject: near-orthogonal views of the left coronary arterial system in the upper panels and of the right coronary arterial system (RCA) of the same patient in the bottom panels. The left heart pair was acquired simultaneously, as was the right at a later point in time. A surface shaded rendering of the 3D reconstruction from the views of Fig. 4a is shown in Fig. 4b. The RCA is positioned anteriorly with the left circumflex artery behind it and the left anterior descending (LAD) to the right. The reconstructed 3D model allows accurate determination of morphometric parameters of the vessel, such as spatial lengths and volumes and mean segment or subtree diameter. The length/volume evaluation can be performed on single vessel segments, on a set of segments, or on subtrees. Avolume model based on generalized elliptical conic sections is created for the selected segments. Volumes and lengths (measured along the vessel course) of those elements are summed and used to derive the mean diameter. In this way, the morphological parameters of a vessel subsystem can be set in relation to the parameters of the proximal segment supplying it. These relations allow objective assessments of diffuse coronary artery diseases.

Given this summary of various clinical techniques, the remainder of this chapter describes the work in our laboratory to develop methods for morphometric and mechanical characterization of pulmonary arterial trees in a basic physiology research laboratory setting.

0 0

Post a comment