Forest representation of vessels in cone-beam computed tomographic angiography

Computerized Medical Imaging and Graphics 29 (2005) 1–14 www.elsevier.com/locate/compmedimag

Forest representation of vessels in cone-beam computed tomographic angiography Zikuan Chen*, Ruola Ning Department of Radiology, University of Rochester, Box 648, 601 Elmwood Avenue, Rochester, NY 14642, USA Received 4 November 2003; revised 8 October 2004; accepted 8 October 2004

Abstract Cone-beam computed tomographic angiography (CBCTA) provides a fast three-dimensional (3D) vascular imaging modality, aiming at digitally representing the spatial vascular structure in an angiographic volume. Due to the finite coverage of cone-beam scan, as well as the volume cropping in volumetric image processing, an angiographic volume may fail to contain a whole vascular tree, but rather consist of a multitude of vessel segments or subtrees. As such, it is convenient to represent multitudinal components by a forest. The vessel tracking issue then becomes component characterization/identification in the forest. The forest representation brings several conveniences for vessel tracking: (1) to sort and count the vessels in an angiographic volume, for example, according to spatial occupancy and skeleton pathlength; (2) to single out a vessel and perform in situ 3D measurement and 3D visualization in the support space; (3) to delineate individual vessels from the original angiographic volume; and (4) to cull the forest by getting rid of non-vessels and small vessels. A 3D skeletonization is used to generate component skeletons. For tree construction from skeletons, we suggest a pathlength-based procedure, which lifts the restrictions of unit-width skeleton and root determination. We experimentally demonstrate the forest representation of a dog’s carotid arteries in a CBCTA system. In principle, the forest representation is useful for managing vessels in both 2D angiographic images and 3D angiographic volumes. q 2005 Published by Elsevier Ltd. Keywords: Cone-beam CT angiography; 3D vessel tracking; 3D skeletonization; Tree construction; Forest representation; Forest culling

1. Introduction Through the injection of radiopaque dye (contrast agent), blood vessels become visible in X-ray angiography. Twodimensional (2D) X-ray angiography can snapshot the instant bolus distribution inside vessels by a projection image. Therefore it is mostly used for fast dynamic vascular imaging such as coronary angiography. Due to X-ray projection, the 2D angiography lacks the spatial information. Accordingly, biplane and multiple projections (including tomosynthesis) are designed to recover the spatial information. In terms of volumetric representation of a three-dimensional (3D) configuration, these efforts, however, are not as competent as the computed tomographic angiography (CTA) [1–3]. A CTA provides a digital volume

* Corresponding author. Tel.: C1 585 275 5977. E-mail address: [email protected] (Z. Chen). 0895-6111/$ - see front matter q 2005 Published by Elsevier Ltd. doi:10.1016/j.compmedimag.2004.10.001

for an angiographic entity. However, the CTA application is hindered by its intolerably long data acquisition time (O10 s) for dynamic vascular imaging. The cone-beam CTA (CBCTA) [4–6] considerably reduces this timing barrier (typically by a factor 1/10), therefore may find its wide applications in quasi-static volumetric imaging like carotid angiography [7]. Specifically, a CBCTA acquires a set of projection images in one rotation of gantry, in a time period of a few seconds (!10 s). In carotid angiography, it is feasible to maintain the bolus distribution (in the local carotid region) in a few seconds, therefore CBCTA is suitable for volumetric carotid angiography [7]. With an angiographic volume, the task is to extract the artery information (including geometrical and topological features) by vessel tracking. Volumetric segmentation is an essential step for vessel tracking in an angiographic volume [3]. Among various volumetric segmentation techniques, thresholding is the most common one, as it can quickly produce a result on

2

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

demand with a threshold. Compared with 2D X-ray projection, where the spatial superimposition complicates the vessel segmentation in the projection image, the CT reconstruction facilitates the thresholding-based volumetric segmentation because of its faithful digital reproduction of an entity’s interior (in terms of X-ray attenuation coefficient) [8]. Due to diverse factors, in the binary volume are present not only the vessels (desired) but also the nonvessels (undesired). Spatial truncation may multiply the number of vessels in a finite volume or subvolume, in terms of counting connected components. One truncation case is associated with the finite coverage by the cone-beam angle (usually !108). Another case may happen in the postprocessing stage, where the study is usually confined to a part of the angiographic volume (by volume cropping). Consequently, the digital volume of interest (VOI) may consist of a multitude of vessel segments or subtrees, together with non-vessel objects, rather than a whole tree. To efficiently organize the collection of connected components [9], we suggest the use of forest representation [10]. The forest can be culled by getting rid of non-vessels and small vessels, thus resulting in a neat collection of conspicuous vessels. Each vessel in the culled forest can be numerically characterized in terms of its geometric and topological features [3,11], thereby accomplishing vessel tracking. We have reported a skeleton-based algorithm for automatic vessel tree construction/pruning for conventional CTA [3]. In this paper, we also use 3D skeleton to manage and identify foreground components in a binary volume (containing vessels and non-vessels). Specifically, we obtain the longest pathlength from skeleton, and use it in combination with the component’s occupancy to classify a component into vessel and non-vessel. To construct an undirected tree from connected skeletons, we adopt the pathlength-based iterative algorithm, which establishes a tree by progressively providing paths in a decreasing order of pathlength. As will be shown in this paper, this pathlength-based algorithm offers two advantages: (1) it allows the non-unit-width skeletons [12,13], thereby speeding up the 3D skeletonization by terminating the procedure early; (2) it relaxes the root determination requirement for vessel tree construction [3]; and (3) it accomplishes tree construction and tree pruning [3,14] at the same time. This paper is organized as follows. In Section 2, we briefly introduce the angiographic volume reconstruction in CBCTA. In Section 3 we present the forest representation of vessels in CBCTA, which involves volumetric segmentation, component identification, 3D skeletonization, vessel identification, vessel delineation, tree construction and forest culling. In Section 4, we report an experiment on a dog’s carotid angiography, which is followed by a discussion. Finally, in Section 5 is the conclusive summary.

2. Cone-beam CT angiography (CBCTA) Briefly speaking, CBCTA is one kind of angiography representing blood vessels in a digital volume through the use of cone-beam X-ray scanning and volumetric tomographic reconstruction. Fig. 1 illustrates a CBCTA system using a flat-panel detector. The system consists of X-ray source generation, gantry rotation, projection image detection, angiographic entity (patient, animal or phantom), dye injection, and data acquisition/processing. The detailed architecture and specifications have been reported elsewhere [5,6]. In terms of functional modules, a CBCTA consists of (1) projection image acquisition of the entity with and without radiopaque dye injection, (2) digital subtraction, (3) volumetric tomographic reconstruction, and (4) vessel tracking in the angiographic volume. In this paper, we are concerned with vessel tracking in the angiographic volume, especially with the forest representation of the multitudinal vessels and non-vessels therein. Fig. 2 displays the modules of the volume reconstruction in CBCTA, followed by the vessel tracking. We will concentrate on the implementation of the vessel tracking modules in Fig. 2(b). Like conventional 2D projection angiography that eliminates background or non-vessel stationary portions by using digital subtraction angiography (DSA) [3], the CBCTA performs digital subtraction (between the postinjection and pre-injection projection images) before 3D reconstruction. Let P0q ½p; q represent pre-injection image acquired by cone-beam X-ray projection, where the subscript q denotes the gantry rotation angle and [p,q] addresses the pixels of the flat panel detector. While maintaining all the configurations, post-injection projection images are acquired in the second rotation of the gantry (with the vessel filled with radiopaque dye), as denoted by Pq[p,q]. Therefore, two rotations of the gantry produce two sets of projection images: Pq[p,q] and P0q ½p; q, corresponding to with and without radiopaque dyes, respectively.

Fig. 1. Schematic diagram of cone-beam CT angiography (CBCTA).

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

Fig. 2. (a) Flowchart of angiographic volume reconstruction in CBCTA; (b) forest-representation-based vessel tracking modules in an angiographic volume.

Then the digital subtraction angiographic (DSA) image is given by DSAq ½p; q Z jln Pq ½p; q K ln P0q ½p; qj;

q 2½0; 360+ Þ (1)

where the absolute value is adopted to eliminate negative values, the logarithm operation is used to convert projection image intensity to the line integral of X-ray attenuation coefficient along the ray path. With the DSA images, we construct a CBCTA volume by a cone-beam reconstruction (using Feldkamp algorithm [15]) [5,6,8]. Let V represent the 3D CBCTA volume, then we can symbolically express the volume reconstruction process by [8] V Z CBCTfDSAq ½p; q; q 2½0; 360+ Þg

(2)

where the gantry rotation angles and detector array are explicitly designated by q and [p,q], while other setup specifications, such as source-to-detector distance and gantry orbit diameter, are implicitly absorbed in the operator CBCT. In our experiments, we adopted Feldkamp (FDK) algorithm [15] for the implementation of CBCT in Eq. (2). The voxel value, V[m,n,k], is interpreted as the X-ray linear attenuation coefficient. Concerning the CBCTA volume, we should address its finiteness in size. Since a CBCT system acquires the X-ray projection images in one rotation of the circular gantry, its scan coverage is limited by the detector-source geometry. This finite coverage is also required by a meaningful conebeam reconstruction (for example, the FDK algorithm [15] is valid for a small cone-angle, typically !108). Consequently, a CBCTA volume is always confined in a small

3

Fig. 3. Vessel tree and subtrees. The overall image contains a single vessel tree, which is decomposed into subtrees as the image region is cropped. The dotted subregion contains a forest of three subtrees, and the hatched subregion contains a forest of four subtrees.

part, for example, at the neck in carotid angiography. Instead of containing a whole artery tree, a carotid angiographic volume consists of a bunch of vessels (mutually disconnected) due to spatial truncation. Another spatial truncation encountered in tree-object representation is volume cropping that is widely used in volumetric image analysis. Since the handling of a high-resolution, for example a 512!512!512, demands huge memory and tremendous computations, we usually divide the whole volume into subvolumes, or confine the study to a small subvolume (or regions of interest), for example to an array of 128!128!256. Such volume reduction may fail to contain a whole tree, and rather intercept subtrees of that tree, without awareness of the parent tree. In other words, a smaller volume always sees more trees (or subtrees), albeit small and partial, than a larger volume. Fig. 3 illustrates this phenomenon. Specifically, the overall region contains a single tree. The cropped regions, however, contain more than one tree, in their own sovereigns, which are actually the subtrees cut off from a big tree. These facts indicate that an angiographic volume may fail to contain a single artery tree. Instead, it usually contains a multitude of vessel trees (or subtrees), without the knowledge about their parents or the inter-vessel relationships. Therefore, the vessel tracking of CBCTA should not only characterize individual vessels (in terms of geometry, topology, and densitometry) but also manage the collection of vessels and non-vessels.

3. Forest representation of vessels In this section, we perform basic segmentation (globally thresholding) on an angiographic volume and produce a binary volume, in which the foreground components are

4

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

organized by a forest representation. By component separation in the binary volume and 3D skeletonization, we can identify a component by its spatial occupancy and skeleton pathlength. Based on component identification, we perform forest culling and vessel tracking. For convenient reference, the notations are listed in Appendix. 3.1. Volumetric segmentation The digital volume V is represented in a 3D grid support space, as denoted by supp(V), with [m,n,k] addressing the voxel. For a M!N!K grid space, a 8-bit digital volume is expressed by VðvÞ 2½0; 255; v Z ½m; n; k 2suppðVÞ with suppðVÞ

(3)

Z f½m; n; kj1% m% M; 1% n% N; 1% k% Kg For conventional CT, the volume V is obtained by stacking slices, which usually results in K/M and MzN, indicating the inherent anisotropy of the grid space. In such cases, a preprocessing step of inter-slice interpolation or intra-slice decimation operation is often needed for balancing the resolutions in the three principal orthogonal directions. In cone-beam CT, the volumetric reconstruction in an isotropic resolution is readily available, thus leading to MZNZK. Before volumetric segmentation, the digital volume V is often subject to some helpful preprocessing steps, such as denoising, defining region of interest (volume cropping), or sampling rate conversion. Given an appropriate threshold th, the basic volumetric segmentation is a global thresholding operation, as expressed by, ( 1 VðvÞO th VB ðvÞ Z with suppðVB Þ Z suppðVÞ: (4) 0 Otherwise This operation classifies voxels into foreground (1-voxel) and background (0-voxel). It can be easily implemented at a high speed, thus serving as a basic 3D image processing/ visualization tool (suitable for interactively manipulation). In the binary volume VB, we can estimate the sparseness of vessels in the support space by P v2suppðVÞ VB ðvÞ (5) sparse Z MNK where the numerator corresponds to the spatial occupancy of foreground objects (dominated by vessels), and the denominator MNK represents the size of 3D support grid space. Meanwhile, we can visualize the binary volume by isosurface rendering with an isovalue 0. If vessels are spatially occluded (or overlapped) each other from a viewpoint, especially for small vessels staying closely and mutually intermingling in space, 3D visualization is evasive and illusive in rendering. The spatial proximity also

hampers the individual vessel measurement. For example, a dilation operation may bring a touch or merge of two close vessels. However, suppose vessels can be separated (extracted out) according to connectedness, then we can focus on a single vessel and perform vessel measurement in an ample space (as large as the whole support space), without the interference from other vessels. To describe the inter-vessel disconnectedness and intra-vessel connectedness, also including the agglomerates and noise, we suggest a forest representation for the collection of vessels. 3.2. Forest representation In terms of binary image processing, the 1-pixels in a binary image are usually classified into connected components [9]. This treatment is also applicable to 1-voxels in a binary volume like VB. In doing so, the foreground (the collection of all 1-voxels) in VB is considered as a collection of components, which can be conveniently described by a component forest, or a forest for short. In the forest representation, each component corresponds to a collection of connected 1-voxels (mathematically a connected graph), mutually exclusively occupying a small portion of the support space. In the forest representation of an angiographic volume, what we are interested in are vessel trees, which may get entangled each other in space, but never touch (otherwise they merge into one component). The forest may also contain non-vessel segments and noise. By definitions, an isolated point or a small stick is also treated as a component, and a forest consisting of a single component is a special forest, i.e. a tree. Therefore, the forest representation of a binary volume VB is a good structure for component management. In a 3D binary volume, the 26-adjacency is usually adopted for 1-voxel connection, 6-adjacency is for 0-voxel connection, and 26-adjacency is for the connection between 1- and 0-voxels. In a 3!3!3 neighborhood, its central voxel has 26 26-adjacent voxels and six 6-adjacent voxels. Based on this convention, we can partition the foreground (1-voxels) in VB into connected components. With the forest representation, the partition scheme can be considered as identifying components in the forest. Let FC denote the component forest consisting of I components (as denoted by Ci, iZ1,2,.,I), then the forest representation can be expressed by FC Z fC1 ; C2 ; .; CI g; with Ci h Cj Z :; i sj:

(6)

The components are mutually exclusive in the support space. For visualization purpose, each component can be extracted from the forest, but represented in situ in the support space. This consideration can be expressed by VB Z

I X iZ1

V Ci

(7a)

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

( with VC Z

1;

v 2Ci

0;

otherwise

(7b)

VCi !VCj Z 0; i sj

(7c)

suppðVCi Þ Z suppðVB Þ Z suppðVÞ:

(7d)

It is noted that Eqs. (6) and (7a)–(7d) represent the same foreground in the binary volume by two different formats. Specifically, Eq. (6) describes the forest of the foreground components in terms of set theory, which has the advantage in storage due to the omission of vast 0-voxels (background), while Eq. (7a)–(7d) represents the foreground components in form of 3D array, which is useful for in situ visualization in the whole support space. The spatial occupancy of a binary component can be easily determined by summing up the 1-voxels in VCi, or counting the elements in Ci. For the forest representation, it is observed that X X X cardðFC Þ Z cardðCi Þ Z VCi ðvÞ!!MNK i

i

v2suppðVÞ

5

manipulated in an ample space, i.e. supp(V), for 3D visualization. This advantage is illustrated in Fig. 4, where (a) represents a binary volume resulting from a volumetric segmentation, which consists of four components: two vessel trees, one isolated stick and one cluster. Obviously, its visualization suffers from spatial occlusion, and the vessel analysis in the binary volume is spatially jostled each other due to proximity. With the component separation, vessel analysis/visualization can be carried out individually in an ample environment in the absence of others as illustrated in (b1)–(b4), where each component assumes its spatial configuration in situ in the whole support space. Besides manipulating the individual components in an ample space, we can also allow the presences of selected components, while hiding others. Fig. 4(c) displays the first two longest components, in a neat exhibition. In this sense, we say that the forest is culled. So far, vessels, as well as non-vessels, are represented in components in a forest. For vessel tracking purpose, vessel information (geometrical and topological features) can be efficiently grasped from skeletons. Therefore, we come to 3D skeletonization for component characterization.

(8) where card() denotes ‘cardinality’, meaning counting the elements in a set [3]. At this stage, we can sort the components according to its spatial occupancy, i.e. by card(Ci). Generally speaking, a large component corresponds to a large vessel tree (or vessel segment), and a small component usually results from an isolated stick or noise. Without interference or occlusion illusion from other components, each component is

3.3. 3D skeletonization and tree construction Skeletons or midlines provide the essential characteristics of a 3D tree-like object, which can be generated by any existent 3D thinning algorithms [12,13]. The skeletonization algorithms involve designing a set of deleting templates. The spatial arrangement of the voxels in the neighborhood of the object voxel under consideration is

Fig. 4. Illustrations of component forest and forest culling. (a) A binary volume consisting of a forest of four components; (b1)–(b4) the individual components; (c) a forest culling by getting rid of the small components.

6

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

checked against these templates, and a voxel is marked for deletion if this arrangement is compliant to one of them. The deleting templates are required to preserve the topology and geometry. As aforementioned and will be seen later, it is not necessary to generate a unit-width skeleton, which is a strong condition pursued by 3D skeletonization algorithm. This finding offers two advantages. First, is allows us to speed up the 3D skeletonization procedure. Second, we are no longer particular with implementation algorithms. In our program, we adopt the deleting templates as reported in [12], and modify the algorithm reported in [3]. The revised algorithm is stated as follows. Algorithm 1. (3D skeletonization allowing non-unit-width voxels) Repeat for each 1-voxel in VB 1. Find the border voxels 2. For each border voxel, check its 3!3!3 neighborhood against the templates in [12] 3. Delete the matched 1-voxels, record the number of deleted voxels. 4. If the number of deleted voxels is less than a specific value, say 10, then halt. Otherwise, continue the loop. A border voxel is defined as a 1-voxel that is 6-adjacent to a 0-voxel. During 3D thinning process, border pixels are gradually peeled off. As soon as a matching case occurs, the border 1-voxel is deleted by simply marking it as a nonforeground voxel. Each loop involves matching comparisons with six basic templates (in six directions: UP, DOWN, NORTH, SOUTH, EAST, and WEST), and their rotations of 90, 180 and 2708 around the directions [3,12]. As the border voxels are peeled off, the interior voxels become border voxels. In the first few epochs, the number of deleted border voxels is considerable. As the loop proceeds, this number decreases drastically. Since each loop should match all templates, the elapse time in each epoch is almost invariant to the number of deletable voxels. In the pursuit of unit-width skeleton, even a few deletable voxels demands another loop. Consequently, the skeletonization progression becomes less efficient in its last few epochs. By allowing the non-unit-width voxels (as seen soon), we can terminate the program early, thereby speeding up the 3D skeletonization procedure. The time gain is implemented by inserting a comparison step 4 in Algorithm 1, where the number of allowed non-unit-width voxels can be set to a default value, say 10. In comparison, the default number associated with the algorithm in [3] is zero, meaning that the looping continues until no deletable voxel is found. Upon completion of 3D skeletonization, each component Ci produces a connected skeleton, as denoted by Gi. Mathematically, the skeleton is a subset of the component, as expressed by Gi Z skelðCi Þ 4Ci

(9)

where skel(Ci) denotes extracting skeletons of component Ci. Since the thinning process removes majority of the 1voxels in a component, it is usually observed that card(Gi)/card(Ci). With the skeleton, we can estimate the thickness of the stick-like vessels by a ratio card(Gi) /card(Ci). Specifically, a thicker vessel produces a smaller ratio. In the extremity, a large ratio may indicate a nonvessel component, such as a solid clump or a ball. Since the skeletonization preserves the topological connections, Gi is a connected graph. Now, our task is to span a tree from the connected graph. Since the vessels are very sparsely populated in an angiographic volume (as characterized by Eq. (5)), the most appropriate way of representing the tree structure in computer is using a linked list. Usually, spanning a tree from a connected graph needs determining a starting point as the tree root, especially for constructing a directed tree. In 3D grid space, we refer to a line-end point of skeleton as an endvoxel. In most cases of vessel tracking, it is not necessary to construct a directed tree from a connected graph. Rather, an undirected tree structure is needed to vessel delineation. We herein propose an undirected-tree-construction algorithm, which does not require a root. Algorithm 2. (Tree construction) Find all the endvoxels in the skeleton;

1. 2. 3. 4.

Start from any one of endvoxels and trace the skeleton voxels; Repeat: Find the maximum pathlength (in the first epoch, record it by Pi1 for late use); Delete voxels on the maximum pathlength (including endvoxels on the path); Update the skeleton; Trace skeleton voxels from endvoxels until reaching any deleted voxel; Until all endvoxels are deleted.

During the process, the deletion of a path is done by assigning its voxels to non-skeleton labels (different from background voxels). The algorithm updates the skeleton and iterates the deletion and tracing process. In each epoch of the iteration, the tracing always starts from an endvoxel and stops at a deleted voxel (generated in previous epochs), except the first epoch where it stops at another endvoxel. The number of loops in Algorithm 2 depends on that of endvoxels. To reduce this number, we impose stronger conditions on the endvoxel definition. Let N3!3!3(v) 1 represent the 3!3!3 neighborhoods of v, and N3!3!3 ðvÞ the collection of 1-voxels therein, analogously for N5!5! 1 5(v) and N 5!5!5(v), i.e. 1 ðvÞ Z fv 0 jv 0 2N3!3!3 ðvÞ; Cðv 0 Þ Z 1g N3!3!3

(10a)

N3!3!3 ðvÞ Z fv 0 jmaxðjv 0 K vjÞ Z 1g

(10b)

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14 1 N5!5!5 ðvÞ Z fv 0 jv 0 2N5!5!5 ðvÞ; Cðv 0 Þ Z 1g

(10c)

N5!5!5 ðvÞ Z fv 0 jmaxðjv 0 vjÞ Z 2g

(10d)

then, we define a reinforced endvoxel by 1 1 1 cardðN3!3!3 ðvÞnfvgÞ Z cardðN5!5!5 ðvÞnN3!3!3 ðvÞÞ Z 1

(11a) !vv1 v2 O is a 26 K adjacent path; for v1 1 1 1 2N3!3!3 ðvÞnfvg;v2 2N5!5!5 ðvÞnN3!3!3 ðvÞ

(11b)

where N5!5!5(v)\N3!3!3(v) stands for the collection of 1voxels in N5!5!5(v) but not in N5!5!5(v), and N3!3!3(v)\{v} stands for the collection of 1-voxels in N3!3!3(v) minus the voxel v under consideration. The conditions (11a) and (11b) mean that the 1-voxel in 26-adjacent to exact one 1-voxel which itself is adjacent to another different 1-voxel, hence making a 26-adjacent path of length 3 in N5!5!5(v). The reinforcement condition requires at least two-voxels clearance from the voxel in question to other unconnected skeleton voxels. Due to a small number of reinforced endvoxels, the computation in the tree construction is reduced considerably.

7

It is not important for tree construction to begin with a root in our algorithm. This conclusion is explained in Fig. 5 by arbitrarily choosing a starting point from the identified endvoxels. For the skeletons in Fig. 5(a), one procedure (Fig. 5(b1)–(b3)) starts from an endvoxel at the top, while another procedure (Fig. 5(c1)–(c3)) starts from an endvoxel at the bottom; and in the result, both produce the same tree structure after three epochs of iteration, in terms of the visited paths. The allowance of non-unit-width skeleton is also explained in the illustration of Fig. 5, where the presence of a non-unit-width skeleton portion (marked by a cluster of grayed voxels) does not influence the tree construction. (In other cases, influence from the non-unit-width skeleton may exist, but ignorable). The insignificance of non-unit-width skeleton can also be seen from the tree pruning ability of the tree construction algorithm [3]. Since the algorithm sorts the paths in a descending order of pathlengths, the tree structure can be pruned by simply discarding those paths whose pathlengths are shorter than a specific length. For example, in Fig. 5(b3) and (c3), a neat tree consisting of conspicuous branches is constructed by three epochs of the iteration. As the tree pruning comes into effect, it brings an early termination of the iterative process, equivalently speeding

Fig. 5. Illustrations for constructing undirected trees from non-unit-width skeletons, by two different root selections. (a) Skeletons with nine endvoxels; (b1)– (b3) the first three epochs during tree construction that start from an endvoxel at the top; (c1)–(c3) the first three epochs that start from an endvoxel at the bottom. The non-unit-width skeletons have no effect on the pathlength-based tree construction.

8

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

Fig. 5 (continued)

up tree construction. In the result, each skeleton Gi gives rise to a tree, denoted by Ti, as expressed by Ti Z treeðGi Þ Z g Pij ðvÞ 4Gi v2G

(12)

where Pij represent the jth path (generated in the jth epoch) of skeleton Gi. The tree construction from skeletons gives rise to card(Ti)%card(Gi), where the strict inequality implies that the pruning operation has come into action. Moreover, the noises and spurs removed by the algorithm can be quantitatively measured by card(Gi)-card(Ti). 3.4. Forest culling and vessel delineation For each component, we extract the skeleton and further construct the tree structure. As a result, the component forest evolves in a general sense as a tree forest, in which both vessels and non-vessels are considered as trees. A stick or deserted vessel segment is considered as a special vessel tree, which has only two endvoxels, i.e. a branchless tree. However, a bubble-like component also results in a sticklike skeleton, its discrimination from a stick vessel should

resort to additional parameters, for example, the spatial occupancy. To characterize the shape of tubular vessels, we define a shape parameter by longest pathlength cardðPi1 Þ ci Z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi spatial occupancy cardðCi Þ

(13)

where Pi1 represents the longest pathlength in the skeleton, which is obtained in the first epoch in the execution of Algorithm 2. The Pi1 collects skeleton voxels on the longest path, therefore, card(Pi1) is also interpreted as the height of a tree. A tree-like component produces a large ci, and a bloblike component produces a small one. The trick of using a square root in Eq. (13) is to compress the large variation among the component occupancies. Now, with two parameters, ci and card(Pi1), we can classify the components into vessels and non-vessels by ( vessel; if cardðPi1 ÞO t1 and ci O t2 Ci is (14) non vessel; otherwise where t1 and t2 are default settings. The two-class recognition issue is diagramed in Fig. 6. The condition

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

Fig. 6. Diagram for vessel/non-vessel classification using two parameters: skeleton pathlength and shape factor.

card(Pi1)!t1 excludes the small blobs, isolated twigs, and random noise (in region a). The condition ci!t2 implies a long stick-like pattern, excluding the large round agglomerates (in region g). The region b is dominated by noise. In principle, the total skeleton pathlength card(Gi) can be used to substitute card(Pi1) in Eq. (13). In practice, this consideration may fail to discriminate a porous but connected agglomerate (resulting in a connected complex graph full of cycles) from a tree structure. In the forest representation, we can remove the nonvessel skeleton trees and generate a vessel forest. In such a way, a forest is said to be culled. The forest and its culled version are expressed by FT Z fTi ; i Z 1; 2; .; Ig

(15a)

FT0 Z fTi ; i Z 1; 2; .; I 0 ! Ig

(15b)

FC0 Z fCi ; i Z 1; 2; .; I 0 ! Ig

(15c)

where Eq. (15a) represents the forest representation of constructed skeleton trees, (15b) represents a culled forest of skeleton trees, and (15c) a culled forest of components. The forest culling needs some rules. For example, we suggest a culling rule by the vessel/non-vessel separation scheme in Fig. 6. More reliable and meaningful rules may utilize more tree attributes, such as summed pathlength, summed spatial occupancy, number of branching points. In program implementation, the forest culling can be numerically characterized by the following relations X X cardðFT0 Þ! cardðFT Þ Z cardðTi Þ% cardðGi Þ /

X

i

cardðCi Þ Z cardðFC Þ

9

the procedure; while card(FT)/Sicard(Gi) implies that many small vessels have been removed in tree construction, which can be quantitatively measured by Sicard(Gi)K card(FT). Other relations are interpreted and characterized in similar fashion. Besides forest culling, we can extract individual vessels from the binary volume or the original grey-level angiographic volume by component occupancy. We call this procedure the vessel delineation. In this sense, Eq. (7a)–(7d) is also interpreted as the vessel delineation from the binary volume. Each vessel is measured and visualized with VCi in its original location in the support space, where it is manipulated in an ample space with empty background. The vessel delineation retains the grey-level representation of vessels in a CBCTA volume. Due to the binarization in Eq. (4), all voxels with density values greater than threshold value are unambiguously assigned to 1. This degeneracy is not desired when the CBCTA is used for tracing dye bolus in blood vessels. To restore the densitometry information along the vessel or over the cross section, we can use the solid tree as 3D masks to extract its spatial occupancy in the original volume data. This action lead to what we call vessel delineation in the original volume, as given by ( VCo i ðvÞ

Z

VðvÞ;

v 2Ci

0;

otherwise

o with suppðVCi Þ Z suppðVÞ

(17) 0 Since VCi assume grey-level values, its visualization needs volume rendering, maximum intensity projection 0 (MIP), or slice display. In VCi , the subtle bolus distributions can be scrutinized. As a result, an ample space is provided for true 3D vessel measurement (in terms of both geometry and 0 densitometry) with VCi . A vessel tree can be endowed with vessel features by 3D measurement, thus producing a valued tree [3]. Specifically, the voxel value in the original 3D volume provides the bolus density, the vessel cross-section in the 3D binary component provides lumen diameter; and connectivity number in the 3D skeleton provides the branching feature and ending feature. Most of these features can be extracted from the skeletons. By consulting the binary components and

i

ð16Þ

i

Again, computations are governed by these relationships, which themselves are in turn meaningfully interpreted by the computational values. For example, card(F)ZSicard(Gi) implies neither forest culling nor tree pruning happened in

Fig. 7. MIP displays of an angiographic volume of a dog’s carotid artery, with the projections along (a) x-axis, (b) y-axis, and (c) z-axis.

10

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

the delineated densitometrical subvolumes, a vessel can be numerically characterized.

4. Experimental results

Fig. 8. Isosurface rendering of a binary volume (a component forest) resulting from volumetric segmentation.

In this section, we demonstrate the vessel forest representation scheme and its application to 3D vessel tracking in CBCTA through the study of a dog’s carotid arteries. The experimental settings and angiographic volume reconstruction of our cone-beam CT system have been reported elsewhere [5,6]. Fig. 7 displays three MIPs of the angiographic volume, with the projections along the principal grid axes. The volume is a 180!180!80 8-bit 3D gray-level array, which was obtained by preprocessing the original digital volume (a 512!512!512 array). Specifically, the preprocessing procedure includes (1) grid space decimation by a factor 2, (2) volume cropping, (3) 3D smoothing with a Gaussian kernel in a 3!3!3 neighborhood, and (4) data conversion from 12-bit to 8-bit. We begin the vessel tracking with the volume in Fig. 7, therefore, it is considered as the original grey-level digital volume in this demonstration. Fig. 8 displays the binary volume resulting from volumetric segmentation of Fig. 7, using Eq. (4) with thZ 230. The vessel sparseness in the volume is 0.0168 (calculated by Eq. (5)). In this display, we can see

Fig. 9. The first 18 components in the forest in Fig. 8, sorted in decreasing order of occupancy. Above each figure is the component label (in parenthesis) followed by the spatial occupancy.

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

11

Fig. 9 (continued)

the presence of stick-like vessels as well as the clumpy nonvessels, which obstruct the vessel visualization and measurement. Fig. 9 displays the 18 individual components in Fig. 8. These components are sorted in descending order of spatial occupancy, i.e. card(Ci)Ocard(C iC1), iZ {1,2,.,17}, where i is provided in parenthesis in Fig. 9, together with card(Ci). Usually, the components with large spatial occupancies represent primary vessels. This is not always true. For example, the two subfigures corresponding to iZ{3,4} are two clumps instead of vessels. Such nonvessel components can be characterized (recognized) through the use of the shape parameter defined in Eq. (13). Using the 3D skeletonization, we generate the skeletons for each component in the forest. By constructing trees using Algorithm 2, we obtain the longest pathlengths for the components in Fig. 9, as given by {card(Pil), iZ 1,2,.,18}Z{290, 260, 132, 83, 186, 162, 175, 179, 81, 99, 114, 9, 28, 6, 10, 13, 15, 11}, in units of voxel spacing. Incorporating the spatial occupancy of a component, we calculate the shape parameter ci by Eq. (13). The components are then identified in a feature space spanned by ci and card(Pi1), as shown in Fig. 10. In this diagram, the vessels and non-vessels are separated by two discrimination lines, corresponding to settings t1Z150, and t2Z1.5. The forest culling is achieved by getting rid of the non-vessels and small vessels (pathlength!t1), hence producing a neat

display. Fig. 11 shows the vessel forest resulting from forest culling, it consists of six components iZ{1,2,5,6,7,8}. The vessel tracking needs to extract the geometrical and topological features of vessels. The forest representation provides a convenient arena to manipulate the vessels, individually as well as collectively. The most important features on vessel structure include endpoints, branches,

Fig. 10. Component classification diagram. The number by the marker corresponds to the label (in parenthesis) in Fig. 9.

12

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

Fig. 11. A display of the culled forest, which consists of six vessels, FT0 Z fCi ; iZ 1; 2; 5; 6; 7; 8g.

and lumen diameters, which are usually determined through the skeletons. Fig. 12 displays a typical result of 3D skeletonization, tree construction/pruning, and vessel feature extractions for the two significant vessels (dog’s carotid

common arteries). Specifically, Fig. 12(a) displays the 3D skeletons embedded in the transparent displays of the vessel solids. The result generated from tree construction/pruning is displayed Fig. 12(b), where the tree pruning effect is easily seen. The vessel endvoxels and bifurcations (detected by a procedure reported in [3]) are displayed in Fig. 12(c), serving as attributes in valued trees [3]. Fig. 13 displays three MIPs of the 3D vessel solids delineated from the original volume. Compared with Fig. 7, Fig. 13 displays the vessels of interest in their original densitometry. The elimination of background interference is attributed to the forest culling. To show the usefulness of forest representation for vessel tracking in volumetric angiography, Fig. 14 provides two more experimental results on the dog’s carotid angiography. The angiographic volumes were generated on the same CBCT system, and the forest culling was done using the same procedure. Concerning the forest representation for vessel tracking in angiographic volume, we discuss the following aspects. The forest representation provides a good structure for managing the objects present in a volume, with which it is convenient to describe inter-component relationships and to study individual components. We can examine an individual vessel in an ample space, without interference from other vessels. Either for visualization or for quantitative measurement, the component isolation scheme and forest culling scheme are very useful. More work should extend to

Fig. 12. Skeleton-based vessel tracking. (a) Skeleton generation; (b) skeleton pruning; (c) skeleton features: endvoxels (marked by ‘D’) and branching points (marked by ‘†’).

Fig. 13. MIP displays of the vessels delineated from the original angiographic volume. The arrangement is the same as that used in Fig. 7.

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

13

Fig. 14. Two more examples that show the effectiveness of forest culling. (a1) and (b1) Two anigographic volumes of dog’s carotid arteries, and (a2) and (b2) the results from the forest culling.

grafting vessels (or subtrees) from subvolumes and to assemble into large trees in a large support volume. It is possible to automate the vessel tracking in both CTA and CBCTA. One crucial step in the automatic procedure is the volumetric segmentation, which is commonly based on a thresholding operation. Another key procedure is to automatically construct a tree from a connected skeleton graph. We have suggested to construct an undirected tree in Algorithm 2, where the root determination is not necessary, therefore facilitating the automatic implementation. Due to high-dimension complexity, there exits no perfect solution to 3D skeletonization. Presently, great endeavors are being exerted on the pursuit of unit-width skeleton under the condition of maintaining the geometrical and topological structures. Instead of pursuing a unit-width skeleton generation, we relax this condition for vessel tree construction by a pathlength-based iterative algorithm. We have used this technique to prune as well as to construct vessel trees from non-unit-width skeletons in our previous report [12]. In this paper, we emphasize the allowance of non-unit width in skeleton, which can significantly reduce the 3D skeletonization time (by an early termination). The pathlength-based technique can also lift the requirement of root selection in undirected tree construction, thereby assuring the automatic implementation. Skeletons are very helpful for vessel characterization (including measurement and identification). With skeletons, we can determine vessel segment length, branching, and spatial location and orientation. Vessel lumen diameter should be determined for the binary volumetric component, and the vessel bolus densitometry should consult the original grey-level angiographic volume by a procedure as

called vessel delineation. For skeleton-based forest representation, we suggest a vessel identification scheme using the skeleton pathlength and a vessel shape factor (defined as variation of the ratio length/volume). The purpose of vessel identification is for forest culling. In most cases, the vessels can be identified by the shape factor defined in Eq. (13). The suggestion of two-parameter recognition scheme is for robustness.

5. Summary The goal of vessel tracking in an angiographic volume is to numerically characterize vessels. In this paper, we meet the goal by collectively representing the vessels, including the non-vessels, in a forest and then performing tree construction and forest culling. This forest representation of the foreground components in a binary volume (resulting form volumetric segmentation on an angiographic volume) offers a bunch of advantages for vessel tracking. First, we can decompose the forest into individual components, with each component being studied separately in the whole support volume, thereby facilitating numerical analysis and 3D visualization. At this stage, the components can be classified into vessels and non-vessels. For each component, we extract its skeleton through 3D skeletonization. The longest pathlength in a skeleton, together with the component occupancy, is used to characterize a component. A shape factor (defined by a variant of the ratio length/ occupancy) is used to identify the tubular shapeliness of vessels. A more reliable vessel identification scheme is based on a pattern recognition scheme by two parameters:

14

Z. Chen, R. Ning / Computerized Medical Imaging and Graphics 29 (2005) 1–14

the longest pathlength and the shape factor. By getting rid of non-vessels and small vessels, we cull the forest and thereby provide a neat collection of conspicuous vessels. Subsequently, vessel delineation can be readily done by extracting the spatial occupancy of a component from the original angiographic volume, where the densitometry of the bolus distribution inside the vessel can be examined. We have justified vessel forest representation in cone-beam CT carotid angiography. This technique is generally applicable to vessel tracking in 2D angiographic images and 3D angiographic volumes. Future work includes grafting the subtrees (distributed over partitioned subvolumes) and thereby restoring larger parental vessel trees (in a larger support volume or a merged volume), as an inverse to the forest representation of subtrees.

Acknowledgements This project was supported in part by NIH grant 2R01HL4803-07 and 1R01CA85904-01.

[4] Ning R, Kruger RA. Image intensifier-based computed tomography volume scanner for angiography. Acad Radiol 1996;3:334–50. [5] X. Wang, Volume tomographic angiography. PhD Dissertation, University of Rochester; 1997. [6] Ning R, Chen B, Yu R, Conover D, Tang X, Ning Y. Flat panel detector-based cone-beam volume CT angiography imaging: system evaluation. IEEE Trans Med Imaging 2000;19:949–63. [7] Leclerc X, Godefroy O, Lucas C, Benhaim JF, Michel TS, Leys D, et al. Internal carotid arterial stenosis: CT angiography with volume rendering. Radiology 1999;210:673–82. [8] Chen Z, Ning R. Why should breast tumor detection go threedimensional? Phys Med Biol 2003;48:2217–28. [9] Ronse C, Devijver RA. Connected components in binary image: the detection problem. New York: Wiley; 1984. [10] Chen Z, Karim M. Forest representation of wavelet transform and feature detection. Opt Eng 2000;39:1194–202. [11] Wink O, Niessen WJ, Viergever MA. Fast delineation and visualization of vessels in 3-D angiographic images. IEEE Trans Med Imaging 2000;19:337–46. [12] Palagyi K, Kuba A. A 3D 6-subiteration thinning algorithm for extracting medial lines. Pattern Recogn Lett 1998;19:613–27. [13] Borgefors G, Nystro¨m I, Baja GD. Computing skeletons in three dimensions. Pattern Recogn 1999;32:1225–36. [14] Shaked D, Bruckstein AM. Pruning medial axes. Comput Vision Image Understand 1998;69:156–69. [15] Feldkamp LA, Davis LC, Kress JW. Practical cone-beam algorithm. J Opt Soc Am A 1984;1(6):612–9.

Appendix. Notations and connotations (some of them were defined in [3]) Notation

Connotation

V, or V(v) vZ[m,n,k] supp(V) VB(v) or VB C FC G VC N3!3!3(v) 1 N3!3!3 (v) N5!5!5(v) 1 N5!5!5 (v) card skel tree T FT c

3D volume array (grey-level values) Voxel coordinates in 3D grid Support space of the volume V Binary volume resulting from V Binary connected component Component forest Collection of skeleton voxels Binary representation of C in supp(V) 3!3!3 neighborhood of voxel v Collection of 1-voxels in N3!3!3(v) 5!5!5 neighborhood of voxel v Collection of 1-voxels in N5!5!5(v) Counting the elements in a set Skeletonization operation Tree construction from skeleton Mathematical tree Tree forest Vessel shape parameter

References [1] Kim HC, Min BG, Lee TS, Lee SJ, Lee CW, Park JH, et al. Threedimensional digital subtraction angiography. IEEE Trans Med Imaging 1982;MI-1(2):152–8. [2] Higgins WE, Karwoski RA, Spyra WJT, Ritman EL. System for analyzing high-resolution three-dimensional coronary angiograms. IEEE Trans Med Imaging 1996;15:377–85. [3] Chen Z, Molloi S. Automatic 3D vascular tree construction in CT angiography. Comp Med Imag Graphics 2003;27:469–79.

Zikuan Chen received his BS (in 1985) and MS (in 1998) in physics from Yunnan University and his DrEng (in 1993) in optical information processing from Nankai University. From 1993 to 1997 he was an assistant professor and an associate professor with the Institute of Modern Optics, Nankai University, working on fast fingerprint identification and biometrics applications in security access control. From 1998 to 2000, he was a postdoctoral fellow at the University of Tennessee and the University of Arkansas, working on multiresolution image processing and fast pattern recognition techniques. From 2001 to 2002, he joined the University of California-Irvine as a research associate, working on medical image analysis (blood vessel tracking and measurement in X-ray angiographic images). Since 2003, Dr Chen has been working on cone-beam CT breast imaging and volumetric angiography in the University of Rochester. He is a member of SPIE and IEEE.

Ruola Ning received a BS degree in Electronic-Physics from Zhongshan University in Guangzhou, China in 1982, a MS degree in Electrical Engineering in 1986 and a PhD degree in Electrical Engineering in 1989 from the University of Utah. He has been involved in the development of an image-intensifier based volume CT system since 1985 at the Medical Imaging Research Laboratory, University of Utah School of Medicine. He is a member of IEEE, SPIE and AAPM. From 1989, he has been on the faculty of University of Rochester. He is an ABR certified medical physicist and is the author of more than 60 referred journal articles and proceedings. He is currently an associate professor of Radiology, Biomedical Engineering, Electrical and computer engineering, Radiation Oncology and Oncology at the University of Rochester School of Medicine, where his research interests include the development of cone-beam reconstruction algorithms, cone-beam CT imaging, cone-beam CT breast imaging, cone-beam CT angiography imaging, 3D and 4D medical imaging, and computer aided detection and diagnosis.