- Email: [email protected]
Automatic navigation path generation based on two-phase adaptive region-growing algorithm for virtual angioscopy
Medical Engineering & Physics 28 (2006) 339–347
Automatic navigation path generation based on two-phase adaptive region-growing algorithm for virtual angioscopy Do-Yeon Kim a,∗ , Sung-Mo Chung b , Jong-Won Park a a
b
Department of Information and Communication Engineering, Chungnam National University, 220 Gung-Dong, Yuseong-Gu, Taejon 305-764, Republic of Korea Department of Plastic & Reconstructive Surgery, Chonbuk National University, Jeonju, Republic of Korea Received 21 October 2004; received in revised form 12 July 2005; accepted 12 July 2005
Abstract In this paper, we propose a fast and automated navigation path generation algorithm to visualize inside of carotid artery using MR angiography images. The carotid artery is one of the body regions not accessible by real optical probe but can be visualized with virtual endoscopy. By applying two-phase adaptive region-growing algorithm, the carotid artery segmentation is started at the initial seed, which is located on the initially thresholded binary image. This segmentation algorithm automatically detects the branch position with stack feature. Combining with a priori knowledge of anatomic structure of carotid artery, the detected branch position is used to separate the carotid artery into internal carotid artery and external carotid artery. A fly-through path is determined to automatically move the virtual camera based on the intersecting coordinates of two bisectors on the circumscribed quadrangle of segmented carotid artery. In consideration of the interactive rendering speed and the usability of standard graphic hardware, endoscopic view of carotid artery is generated by using surface rendering algorithm with perspective projection method. In addition, the endoscopic view is provided with ray casting algorithm for off-line navigation of carotid artery. Experiments have been conducted on both mathematical phantom and clinical data sets. This algorithm is more effective than key-framing and topological thinning method in terms of automated features and computing time. This algorithm is also applicable to generate the centerline of renal artery, coronary artery, and airway tree which has tree-like cylinder shape of organ structures in the medical imagery. © 2005 IPEM. Published by Elsevier Ltd. All rights reserved. Keywords: Carotid artery; Virtual angioscopy; Adaptive region growing; Navigation path; Medical image segmentation
1. Introduction With the increasing size and number of medical images, the use of computers in facilitating their processing and analysis has become necessary. In particular, computer algorithm for the visualization and delineation of anatomical structures and viewing the inner surfaces of organs are a key component in assisting and automating specific radiological tasks. Endoscopy is a diagnostic method using optical, video-assisted technology to view the inner surfaces of hollow organs in a continuous fashion. By changing the position ∗ Corresponding author at: 108-302 CheongGu Apt., 462-4 JeonminDong, Yuseong-Gu, Taejon 305-729, Republic of Korea. Tel.: +82 42 868 4342; fax: +82 42 861 1488. E-mail address: [email protected] (D.-Y. Kim).
of the endoscope, the operator can view the inside of an organ while controlling the viewing position and angle of the probe. During this invasive and interactive exploration, the endoscopist has full control of navigation within the hollow organ. Physical endoscopy applies the optical probe; however, virtual endoscopy is a new method of diagnosis using computer processing of three-dimensional (3D) image data sets. These are computed tomography (CT) and magnetic resonance imaging (MRI) scans to provide simulated visualizations [1,2] of specific organs of a patient. Visualization is the process of exploring, transforming, and viewing data as images to gain understanding and insight into the data. The visual effects of virtual endoscopy are similar or equivalent to those produced by standard endoscopic procedures [3]. Conventional CT and MRI scans produce axial slices of the
1350-4533/$ – see front matter © 2005 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2005.07.011
340
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
body that are viewed sequentially by radiologists who must imagine or extrapolate from these views what the actual 3D anatomy should be. By using sophisticated algorithms and high performance computing, these axial slices may be rendered as direct 3D representations of human anatomy. Specific anatomic data appropriate for realistic endoscopic simulations can be obtained from 3D MRI digital imaging examinations [4] or 3D acquired spiral CT data [5]. In fact, virtual endoscopy provides viewing control and options that are not possible with real endoscopy such as direction, angle of view, scale of view, immediate translocation to new views, lighting, and measurement. In addition, there are many body regions not accessible by optical endoscopy but can be explored with virtual endoscopy. Several important body systems are compatible with invasive probes including the heart, spinal canal, inner ear, biliary, pancreatic ducts, and large blood vessels. These are important anatomic structures ideally suited for virtual endoscopy [6]. By contrast, optical endoscopic procedures can be uncomfortable and sedation or anesthesia may be required [7]. Furthermore, optical endoscopic procedures are invasive and have serious side effects such as perforation, infection and hemorrhage. Therefore, virtual endoscopic visualization avoids the risks associated with real endoscopy, and when used prior to performing an actual endoscopic examination can minimize procedural difficulties and decrease the rate of morbidity, especially for endoscopists in training.
2. Related work The navigation paths for virtual camera play a crucial role in virtual endoscopy application. The navigation technique decides how the physician controls the position and orientation of a virtual camera to examine the inner surface of organ. The planned navigation scheme pre-computes a flythrough path to automatically move the virtual camera from one end of the organ to another end. On the other hand, the free navigation scheme aims to allow a physician to interactively control the virtual camera position and orientation. The physician is required to control the virtual camera at each frame to navigate through the virtual model. Although the physician is specifically well trained, it can be difficult to navigate to target. This may cause the problem that the camera position and orientation are lost in the virtual environment. The guided navigation scheme has combined the benefits of both planned and interactive navigation. By using a guided navigation scheme, the virtual camera moves automatically from a source point to a target point along a predefined flythrough path. When necessary, the physician can take over the control at any time by interactively and easily adjusting the camera position and direction. The switch scheme between these two modes is seamless, resulting in a smooth navigation [8]. For the purpose of guided navigation scheme, the determination of fly-through path is an essential and prerequisite
task. The well-defined flight paths provide a stable viewing direction, avoiding sudden twists and turns. The manually planning flight path such as key framing [2] is a tedious and time-consuming task. The key framing is that the user specifies manually the pose of a small subset of the total number of frames to be rendered. Then, interpolating curves are fitted to generate a continuous path and set of frames. The other method called distance mapping, a technique used in robot path planning, has also been applied to navigation path planning for virtual camera. A goal voxel is selected within the volume and a distance to the goal is computed for each voxel within the structure of interest. The distance assigned is the shortest distance along a path to the goal through the structure. However, the shortest distance approach has a tendency to produce paths that hug the wall of the organ rather than follow the central axis, thus limiting what is visible to the virtual camera. Although the improvement has been conducted by penalized distance [9], the distance mapping still requires substantial computing resources. On the other hand, the medial axis is the locus of the centers of all circles that are tangent to the boundary of the object at two or more disjoint points. Based on the medial axis transform (MAT) [10], continuously generated centerline has been used as a fly-through path for virtual camera [11]. However, MAT is heavily dependent on the circular shape of object on 2D image space. As a result of unstable viewing direction and abrupt twists, operator interaction is required to adjust the only jagged navigation path. In addition, the skeleton of a volumetric data provides a compact description of interesting object. The skeletonization also provides an efficient method for visualization and analysis such as feature extraction, feature tracking, and automatic navigation [12]. The topological thinning or onion peeling [13,14] is one of the skeletonization methods that one layer at a time of voxels is peeled off the object until just the skeleton remains. The skeleton algorithm automatically provides a 3D central axis of structure; however, onion peeling is computationally expensive. In this paper, therefore, we propose a fast and automated navigation path generation algorithm for virtual angioscopy based on two-phase adaptive region-growing and bounding box method. In the following, Section 3 describes the two-phase adaptive region-growing algorithm for carotid artery segmentation and then Section 4 provides a navigation path generation method for virtual angioscopy. Section 5 describes the mathematical phantom for validation of segmentation and navigation path generation. In Section 6, experimental results of segmentation and navigation path generation on both phantom and patient data sets are presented. Finally, Section 7 provides the discussions and concluding remarks.
3. Two-phase adaptive region growing The 3D time-of-flight (TOF) carotid artery (CA) MRA images are used to automatically generate the navigation
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
path for virtual angioscopy based on two-phase adaptive region-growing segmentation algorithm. The carotid artery is segmented as the following sequence: (1) image analysis, (2) initial thresholding for seed selection, (3) 2D region growing, (4) bifurcation detection, and (5) 3D region growing. Fig. 1 shows the overall segmentation process of carotid artery.
341
3.1. Image analysis The pairs of carotid and vertebral artery in the neck mainly deliver the blood to the brain about 80%, 20% of blood flow rate, respectively. Anatomically, the carotid artery is classified into common carotid artery (CCA), internal carotid artery (ICA), and external carotid artery (ECA). The CCA is branched off ICA and ECA at the position of thyroid cartilage and then connected to Circle of Willis through anterior and posterior communication artery. In this paper, the image analysis is defined as preliminary tasks to automatically segment the carotid artery. These consist of two image examination sub-tasks, which are: (1) the gray level distribution between objects and background and (2) the anatomic structure of carotid artery. Fig. 2(a) shows the pairs of CCA and vertebral artery (VA) before branch; Fig. 2(b) shows the ICA and ECA just after bifurcation from CCA. As shown in Fig. 2, CA and VA are brighter than surrounding tissues and background. This is the nature of MRA scanning protocol that focus on the blood flow area and suppresses the signal intensity of surrounding tissues. Fig. 3 shows the histogram of Fig. 2 that is plotted the frequency of intensity occurrence in the image [15–17]. Compared to background area and surrounding tissues in Fig. 2, blood flow area is small portion of the image. As indicated in the histogram, the frequency of gray level occurrence is biased to the left side. Consequently, the dominant background pixels provide uni-modal shape of histogram. While global and local scheme of thresholding is difficult to segment the object or region, region growing is suitable for the delineation of small structures such as tumors [18], lesions, and blood vessels in the medical imagery. Regarding to the anatomic structure of carotid artery, a priori knowledge is used as the following: (1) spatial information between VA and CA for seed selection of region growing, (2) spatial information and branch direction of ICA and ECA for navigation path generation. 3.2. Initial thresholding for seed selection
Fig. 1. Block diagram of carotid artery segmentation process.
Image segmentation is the process that partitions a digital image into disjoint connected sets of pixels, each of which corresponds to an object or region. The various image processing algorithms could be combined to segment the objects on the image space. In the case of pipeline operation, thresholding technique could be used as an initial or intermediate processing step to completely segment the specific objects or regions. Thresholding is a pixel classification process to identify the pixels of a given image into two classes: those pertaining to objects and those pertaining to background. While one class includes pixels with gray values that are below or equal to a certain threshold value, other class includes those with gray values above the threshold. Thresholding is the transformation of an input image f to an output binary image g as follows:
342
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
Fig. 2. Carotid artery MRA axial image: (a) before and (b) after bifurcation.
g(x, y) =
1
for f (x, y) ≥ T
0
for f (x, y) ≺ T
(1)
where T is the threshold, g(x,y) = 1 for image elements of objects, and g(x,y) = 0 for image elements for the background. In order to determine the region-growing seed in the sequence images, the initial segmentation is performed with threshold value T = (2h − 1)/2, where h is the number of bits for representation of each pixel. Fig. 4 shows the initial segmentation result that is the binary version of Fig. 2(a). The CA is isometric and almost symmetrically resides in left and right side of the neck, even though anatomic structures are varied among the subjects. By applying a priori knowledge of CA anatomic structure, the region R is divided into two separate sub-regions to locate two initial seeds on R1 and R2 for pairs of CA segmentation, respectively. The sub-region R1 and R2 are as follows: R1 = {(x, y), 1 ≤ x ≤ xm , 1 ≤ y ≤ yn /2}
(2)
R2 = {(x, y), 1 ≤ x ≤ xm , yn /2 + 1 ≤ y ≤ yn }
(3)
although the (row, column) orientation used in matrices is also often used in digital image processing. These sub-regions R1 and R2 are used to recognize the left and right carotid artery for region-growing seeds, respectively. In consideration of the spatial relationship between VA and CA as shown in Fig. 4, the every pixel, which have value of g1 (x,y) = 1, is examined from the top left corner toward the horizontal direction on each sub-region. By using eight connectivity features, the connected neighboring pixels are labeled as carotid artery starting at the first located pixel that have value of g1 (x,y) = 1. The labeled component of each sub-region is the binary version of carotid artery segmentation result in the first sequence image that can be used as seed region for further region growing. 3.3. 2D region growing The image can be modeled by a continuous function of two or three variables; in the 2D case arguments are coordinates (x,y) in a plane, while if images are expended to 3D a third variable might be added. The image sequence function F is as follows:
where xm , yn represent maximal image coordinates. The customary orientation of coordinates in an image is in the Cartesian fashion (horizontal x-axis, vertical y-axis),
Fig. 3. Histogram of carotid artery MRA axial image (Fig. 2).
Fig. 4. Binary version image (Fig. 2(a)).
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
343
Fig. 5. Candidate requirements of start pixel selection for region growing: (a) segmented image (current slice), (b) unsegmented image that have smaller region (next slice) than segmented current slice image (a), and (c) unsegmented image that have larger region (next slice) than segmented current slice image (a).
F=
n
fi (x, y)
(4)
i=1
where n is image sequence number, x and y are the spatial coordinates. CT and MRI scans produce series of axial slice images. These volumetric data could be used as input image sequence for segmentation of specific object or region. Based on Eqs. (1) and (4), the segmented output image g1 is a binary form of the input image f1 . In order to apply the region growing on the first sequence image, the intensity values are averaged only for CA region in image f1 , which is corresponding region of labeled component in image g1 . This averaged value is used as homogeneity criterion for region growing on the first sequence image f1 starting at arbitrary selected pixel with eight-connectivity feature. After the region growing is finished on the first sequence image f1 , homogeneity criterion is adjusted for segmentation of next sequence image f2 . Based on the connectivity-preserving features of sequence image, region growing is extended to the next slice image with adaptive homogeneity criterion that is re-calculated in the current slice image. In order to select the start pixel of region growing for next slice segmentation, the pixel must meet the following criteria. These conditions are: (1) pixel must be labeled in the current slice, (2) pixel must not be labeled in the next slice, (3) coordinates must be same on the current and next slice, (4) pixel of next slice must be passed the homogeneity test, and (5) next slice must not be last slice. Fig. 5 shows the candidate requirements of start pixel selection for next slice image segmentation. These example images (10 × 10) provide detail explanation of top three conditions that are described above to select the start pixel for region growing. Fig. 5(a) shows the segmentation result on current slice image. Fig. 5(b) shows the un-segmented next slice image that have smaller region; Fig. 5(c) shows the un-segmented image that have larger region than segmented current slice image. After start pixel is selected, pointer for current slice number is incremented and region growing is performed on next slice image. The 2D region-growing process is performed until next slice number is not last slice image. 3.4. Bifurcation detection After segmentation is finished on sequence image f2 , the face connection test is performed to detect bifurcation
between current slice image f2 and previous slice image f1 . In order to find the branch in the current slice, every pixel in the current slice must meet the following criteria. These conditions are: (1) pixel must be labeled in the previous slice, (2) pixel must not be labeled in the current slice, (3) coordinates must be same on the previous and current slice, and (4) pixel of current slice must be passed the homogeneity test. If the bifurcation is detected in the current slice image, the current slice number and homogeneity criterion are saved into the stack S and retrieved from the stack S to segment the other branch (ICA or ECA). In addition, stack S3D is duplicated from the stack S for second-phase 3D region growing. The stack S and S3D are identical that have slice number and homogeneity criterion in the case of bifurcation is detected. The stack S3D is only used to save in first-phase region growing and to retrieve the slice number and homogeneity criterion in second-phase region growing. The branch detection process is performed until next slice number is not last slice image. 3.5. 3D region growing The region growing is a procedure that group pixels or sub-regions into larger regions based on predefined criteria. The basic approach is to start with a set of seed points and from these grow regions by appending to each seed those neighboring pixels that have properties similar to the seed, such as specific ranges of gray level or color [19–21]. In similar, region growing could be expanded to 3D image space with voxel connectivity. The 3D region-growing method perform the homogeneous test from the start voxel (or 3D block) to the neighbor voxel (or 3D block) using graylevel, texture, color as an acceptance criterion. According to the homogeneous test result, the neighbor voxel (or 3D block) is included or excluded until termination condition is satisfied. In addition, the voxel connectivity is defined as the following: (1) six-connectivity with joint faces, (2) 18connectivity with joint edges, and (3) 26-connectivity with joint corners. If the next slice is last slice image and if stack is not empty, slice number and homogeneity criterion are retrieved from the stack S for the other branch segmentation. The 2D region-growing process is performed until next slice is last slice image and if stack is empty. After first-phase 2D adaptive region growing is terminated including branch
344
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
Fig. 7. (a)–(e) Five models of circumscribed quadrangle.
The circumscribed quadrangle is connected to the four tangent lines that run parallel with x- and y-axis in the image, respectively. The four tangent points are located outermost in the segmented carotid artery from left, right, top, and bottom directions, respectively. Fig. 7 shows the five models of circumscribed quadrangle including one circular and four ellipse shapes. The dashed lines represent the bisector of each circumscribe quadrangle. The generated centerline on each axial slice of carotid artery is used as a fly-through path for virtual angioscopy. 4.2. Rendering of endoscopic view
Fig. 6. Two-phase adaptive region-growing algorithm.
segmentation, the second-phase 3D region growing is performed as the following. These process are: (1) the intensity values of labeled pixels, which are from 2D region growing on whole sequence image, are averaged and used as an homogeneity criterion, (2) 26-connectivity feature is used, (3) only apply to the unlabeled voxel for preserving the labeled pixel in the first-phase region growing, and (4) use the duplicated stack S3D for branch handling. To summarize the carotid artery segmentation process, the algorithm for the two-phase adaptive region growing is specified in Fig. 6.
4. Navigation path generation 4.1. Navigation path By using a priori knowledge such as spatial information and branch direction of carotid artery, common carotid artery and internal carotid artery are segmented to determine the navigation path for virtual camera before and after the bifurcation, respectively. In order to visualize the inside of internal carotid artery, a fly-through path is pre-computed to automatically move the virtual camera from common carotid artery to internal carotid artery. Our method for positioning the virtual camera is based upon the intersecting coordinates of two bisectors on the circumscribed quadrangle of segmented carotid artery.
The endoscopic view of the internal carotid artery can be achieved by either surface or volume rendering computer graphic techniques. Surface rendering is an indirect method of obtaining an image from a volume data set. The surfaces are produced by mapping data values onto a set of geometric primitives in a process known as isosurfacing. These isosurfaces can then be rendered into a displayable image using standard computer graphics techniques. On the other hand, volume rendering is a more direct method for reconstruction of 3D structures. Volume rendering represents 3D objects as a collection of cube-like building blocks called voxels. The main advantage of this type of rendering is its ability to preserve the integrity of the original data throughout the visualization process. This technique, however, is generally expensive than conventional surface rendering technique. Although surface rendering lacks rendering quality and information beyond surface in terms of the usage of partial volume data, the surface rendering technique provides interactive rendering speed and usability of standard graphic hardware. Marching cubes is an algorithm for generating isosurfaces from volumetric data known as surface rendering technique. This algorithm uses the information at the corners of a voxel to construct a surface that approximates the original surface. The virtual angioscopy is implemented to provide the endoscopic view of internal carotid artery by using both Marching cubes algorithm [22] and volume rendering algorithm with perspective projection technique.
5. Phantom generation For the purpose of the performance quantification for segmentation method and specific medical application such as proposed automatic navigation path generation algorithm, validation experiments are necessary. Validation is typically
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
Fig. 8. Segmentation result of: (a) common carotid artery before bifurcation and (b) internal and external carotid artery after bifurcation.
Fig. 9. 3D reconstruction view of: (a) phantom and (b) clinical data set using segmentation results.
performed using one of two different types of truth models. The most straightforward approach to validation is by comparing the automated segmentations with manually obtained segmentations [23,24]. This approach, besides suffering from the drawbacks outlined of laborious and time consuming, does not guarantee a perfect truth model since an operator’s performance can also be flawed. The other common approach to validating segmentation methods is by use of physical phantom [25] or computational phantoms which consisting of digital (voxel-based) phantom [26,27] and mathematical (analytical) phantom [28]. Physical phantoms provide an accurate depiction of the image acquisition process but typically do not present a realistic representation of anatomy. Digital phantoms are mainly derived from segmented tomographic images of the human anatomy obtained by either CT or MRI. Mathematical phantoms consist of regularly shaped continuous objects defined by combinations of mathematical geometries such as spheres, ellipsoids, cylinders, and cones. Although anatomically less realistic than phantoms
345
Fig. 11. Projected centerline for: (a) phantom and (b) clinical data set.
derived from CT or MR images of patients, the mathematical phantom has the advantage that it can be easily modified to simulate a wide variety of patient anatomies. By considering the anatomical structure of blood vessel, the cylinder shape of mathematical phantom is devised to simulate the carotid artery bifurcation. A carotid artery phantom is implemented that consist of two cylinders model without branches (simulation of vertebral artery) and two cylinders model with bifurcation for each cylinder (simulation of carotid artery). The ellipse objects are drawn on 2D image space to simulate the blood vessels with a graphic package. To simulate the tomographic image sequence, the ellipse objects are located at the corresponding coordinates that are determined on clinical data set. For the purpose of intensity value assignment to blood vessels and background, the inside and outside of each ellipse objects are recognized on 2D image space. Based on the intensity assessment, the initial threshold value T = (2h − 1)/2, where h is the number of bits for representation of each pixel, is used to assign the intensity value for blood vessel and background. The intensity above the threshold value T is assigned to blood vessels and below the threshold value T is assigned to background based upon the examination result of intensity distribution between objects and background in the clinical data set. A random number generator also used to determine the specific intensity value both blood vessels and background after range are checked with threshold value T. Regarding to the determination of blood vessel size, the diameter of common carotid artery and internal carotid artery are designed to same radius; external carotid artery and vertebral artery are also implemented with same radius. According to the geometrical structure of carotid artery, isometrically simulated carotid and vertebral artery are located in the sub-regions R1 and R2 , respectively. The branch positions determined on clinical data
Fig. 10. (a)–(f) Center coordinates of segmented carotid artery based upon intersecting coordinates of two bisectors on circumscribed quadrangle.
346
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
Fig. 12. Endoscopic view of carotid artery where virtual camera was positioned at: (a)–(c) common carotid artery and (d)–(f) internal carotid artery by using ray casting algorithm.
set are also used to implement the bifurcations of the carotid artery phantom.
6. Experimental results The Intel P4 processor with MS Windows, Visual C++ and Matlab were used to segment the carotid artery, to generate the navigation path, and to implement the virtual angioscopy system. Our algorithm was applied on six clinical data sets of 3D time-of-flight (TOF) MRA images from GE Signa Horizon Echospeed scanner (TR = 25 ms, TE = 6.9 ms, FA = 25◦ , FOV = 220 mm). Each case consisted of 136 axial images and the dimension of image was 512 × 512 pixels, pixel spacing was 0.43 mm and slice thickness was 1.4 mm. The phantom data set was also used to validate our algorithm. Fig. 8(a) and (b) shows the segmentation result of CCA before bifurcation and ECA and ICA after bifurcation by using two-phase adaptive regiongrowing algorithm, respectively. Fig. 9(a) and (b) presents the fully rotatable external view of 3D reconstructed phantom and clinical data set using segmentation results. As a result of navigation path generation, Fig. 10 shows the center coordinates of segmented carotid artery based upon intersecting coordinates of two bisectors on the circumscribed quadrangle. Fig. 11 also shows the projected centerlines of phantom and clinical data set by using maximum intensity projection algorithm. The times to compute the navigation paths for 12 carotid arteries in six clinical data sets ranged from 19 to 21 s
depending on branch location of each carotid artery including segmentation times. For the purpose of visualization of carotid artery, Marching cubes algorithm was used to generate the endoscopic view. This algorithm provided interactive rendering speed at 15 frames per second, so effective navigation of carotid artery was possible. On the other hand, ray casting is one of the volume rendering techniques using backward projection. Rays are cast from each pixel of the image plane into the volume data. At locations along each ray, a sample value and a surface normal approximation are calculated using values of surrounding voxels. Using the sample value and normal, a sample opacity and color are dynamically assigned by a lookup table or in a preprocessing phase. Then a local shading model is applied and the samples along the ray are composed into a pixel value of the final image [29,30]. Fig. 12(a)–(c) and (d)–(f) shows the endoscopic view of carotid artery where the virtual camera was positioned at common carotid artery and internal carotid artery by using ray casting algorithm, respectively. This algorithm required lots of computing time, so only provided as an off-line navigation of internal carotid artery.
7. Discussion and conclusions Based on two-phase adaptive region-growing algorithm, we have presented a fast and automated navigation path generation method to visualize the carotid artery. This method
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
was validated by phantom data set and applied on clinical data sets. This algorithm is more effective than key framing and topological thinning method in terms of automated features and computing time. This algorithm is also applicable to generate the centerline of cylinder shape structures and tree-like blood vessels such as renal artery and coronary artery. The limitation of this algorithm occurs on occlusion or near occlusion cases resulted in segmentation failure of carotid artery. Caused by further problem on navigation path generation, it is impossible to visualize the inside of carotid artery. In addition, this algorithm is applicable to convex shape of blood vessel. The central axis of carotid artery, which was determined by using two-phase region growing and bounding box method, could be used as navigation path for virtual camera. In the case of that blood vessels run horizontally across the scan plane, however, the images provide non-circular shape of blood vessel. With this limitation, the navigation path cannot be easily determined by proposed algorithm. Hence, the further works are required to generate the oblique slices which scanning planes are perpendicular to the central axis of blood vessel. These slices provide the circular shape of blood vessel even though artery runs horizontally across the scan plane. Furthermore, it is beneficial that circular objects are restored to original shape, especially for distorted objects caused by stenosis, scanning artifacts and possible noise. In conclusion, newly developed virtual angioscopy system aims at supporting physicians in determining the precise spatial location and shape of the carotid artery stenosis that are the critical indicators of the transient ischemic stroke. Combining with a computerized quantification method for carotid artery stenosis [31], virtual angioscopy system will provide better diagnostic results of carotid artery disease.
References [1] Geiger B, Kikinis R. Simulation of endoscopy. In: Proceedings of the AAAI spring symposium series: applications of computer vision in medical image processing. 1994. p. 138–40. [2] Rubin GD, Beaulieu CF, Argiro V, Ringl H, Norbash AM, Feller JF, et al. Perspective volume rendering of CT and MR images: applications for endoscopic imaging. Radiology 1996;199:321–30. [3] Wickham JEA. Minimally invasive surgery: feature developments. Br Med J 1994;308:193–6. [4] Rusinek H, Mourino MR, Firooznia H, Weinreb JC, Chase NE. Volumetric rendering of MR images. Radiology 1989;171:269–72. [5] Fishman EK, Jeffrey RB. Spiral CT: principles, techniques, and applications. Philadelphia: Lippincott, Williams & Wilkins; 1998. [6] Robb RA. Virtual (computed) endoscopy: development and evaluation using the visible human datasets. In: Proceedings of the visible human project conference. 1996. p. 19–20. [7] Jolesz FA, Lorensen WE, Shinmoto H, Atsumi H, Nakajima S, kavanaugh P, et al. Interactive virtual endoscopy. Am J Roentgenol 1997;169(5):1229–35.
347
[8] He T, Hong L, Chen D, Liang Z. Reliable path for virtual endoscopy: ensuring complete examination of human organs. IEEE Trans Visual Comput Graph 2001;7(4):333–42. [9] Bitter I, Kaufman AE, Sato M. Penalized distance volumetric skeleton algorithm. IEEE Trans Visual and Comput Graph 2001;7(3):195–206. [10] Jain R, Kasturi R, Schunck BG. Machine vision. New York: McGraw-Hill; 1995. [11] Kim DY, Park JW. Visualization of the internal carotid artery using MRA images. Magn Reson Imag 2005;23(1):27–33. [12] Zhou Y, Toga AW. Efficient skeletonization of volumetric objects. IEEE Trans Visual Comput Graph 1999;5(3):196–209. [13] Hong L, Muraki S, Kaufman AE, Bartz D, He T. Virtual voyage: interactive navigation in the human colon. In: Proceeding of ACM SIGGRAPH. 1997. p. 27–34. [14] Paik DS, Beaulieu CF, Jeffery RB, Rubin GD, Napel S. Automated flight path planning for virtual endoscopy. Med Phys 1998;25(5):629–37. [15] Gonzalez RC, Woods RE. Digital image processing. New York: Addision-Wesley; 2000. [16] Parker JR. Algorithm for image processing and computer vision. New York: John Wiley & Sons; 1997. [17] Calotto MJ. Histogram analysis using a scale-space approach. IEEE Trans Pattern Anal Mach Intell 1987;9:121–9. [18] Kim DY, Park JW. Computer aided detection of kidney tumor on abdominal computed tomography scans. Acta Radiol 2004;45(7):791–5. [19] Sonka M, Hlavac V, Boyle R. Image processing, analysis, and machine vision. Pacific Grove: PWS Publishing; 1999. [20] Adams R, Bischof L. Seeded region growing. IEEE Trans Pattern Anal Mach Intell 1994;16:641–7. [21] Mehnert A, Jackway P. An improved seeded region growing algorithm. Pattern Recogn Lett 1997;18:1065–71. [22] Lorensen W, Cline H. Marching cubes: a high resolution 3D surface construction algorithm. In: Proceeding of ACM SIGGRAPH. 1997. p. 163–9. [23] Vaidyanathan M, Clarke LP, Hall LO, Heidtman C, Velthuizen R, Gosche K, et al. Monitoring brain tumor response to therapy using MRI segmentation. Magn Reson Imag 1997;15:323– 34. [24] Wust P, Gellermann J, Beier J, Wegner S, Tilly W, Troger J, et al. Evaluation of segmentation algorithms for generation of patient models in radiofrequency hyperthermia. Phys Med Biol 1998;43:3295–307. [25] Lei T, Sewchand W. Statistical approach to X-ray CT imaging and its applications in image analysis. IEEE Trans Med Imag 1992;11: 62–9. [26] Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, et al. Design and construction of a realistic digital brain phantom. IEEE Trans Med Imag 1998;17:463–8. [27] Zubal IG, Harrell CR, Smith EO, Rattner Z, Gindi G, Hoffer PB. Computerized 3-dimensional segmented human anatomy. Med Phys 1994;21:299–302. [28] Zaidi H. Relevance of accurate Monte Carlo modeling in nuclear medical imaging. Med Phys 1999;26:574–608. [29] Watt A. 3D computer graphics. Boston: Addison-Wesley; 1993. [30] Noordmans HJ, Voort H, Smeulders A. Spectral volume rendering. IEEE Trans Visual Comput Graph 2000;6(3):196–207. [31] Kim DY, Park JW. Computerized quantification of carotid artery stenosis using MRA axial images. Magn Reson Imag 2004;22(3):353–9.
Automatic navigation path generation based on two-phase adaptive region-growing algorithm for virtual angioscopy Do-Yeon Kim a,∗ , Sung-Mo Chung b , Jong-Won Park a a
b
Department of Information and Communication Engineering, Chungnam National University, 220 Gung-Dong, Yuseong-Gu, Taejon 305-764, Republic of Korea Department of Plastic & Reconstructive Surgery, Chonbuk National University, Jeonju, Republic of Korea Received 21 October 2004; received in revised form 12 July 2005; accepted 12 July 2005
Abstract In this paper, we propose a fast and automated navigation path generation algorithm to visualize inside of carotid artery using MR angiography images. The carotid artery is one of the body regions not accessible by real optical probe but can be visualized with virtual endoscopy. By applying two-phase adaptive region-growing algorithm, the carotid artery segmentation is started at the initial seed, which is located on the initially thresholded binary image. This segmentation algorithm automatically detects the branch position with stack feature. Combining with a priori knowledge of anatomic structure of carotid artery, the detected branch position is used to separate the carotid artery into internal carotid artery and external carotid artery. A fly-through path is determined to automatically move the virtual camera based on the intersecting coordinates of two bisectors on the circumscribed quadrangle of segmented carotid artery. In consideration of the interactive rendering speed and the usability of standard graphic hardware, endoscopic view of carotid artery is generated by using surface rendering algorithm with perspective projection method. In addition, the endoscopic view is provided with ray casting algorithm for off-line navigation of carotid artery. Experiments have been conducted on both mathematical phantom and clinical data sets. This algorithm is more effective than key-framing and topological thinning method in terms of automated features and computing time. This algorithm is also applicable to generate the centerline of renal artery, coronary artery, and airway tree which has tree-like cylinder shape of organ structures in the medical imagery. © 2005 IPEM. Published by Elsevier Ltd. All rights reserved. Keywords: Carotid artery; Virtual angioscopy; Adaptive region growing; Navigation path; Medical image segmentation
1. Introduction With the increasing size and number of medical images, the use of computers in facilitating their processing and analysis has become necessary. In particular, computer algorithm for the visualization and delineation of anatomical structures and viewing the inner surfaces of organs are a key component in assisting and automating specific radiological tasks. Endoscopy is a diagnostic method using optical, video-assisted technology to view the inner surfaces of hollow organs in a continuous fashion. By changing the position ∗ Corresponding author at: 108-302 CheongGu Apt., 462-4 JeonminDong, Yuseong-Gu, Taejon 305-729, Republic of Korea. Tel.: +82 42 868 4342; fax: +82 42 861 1488. E-mail address: [email protected] (D.-Y. Kim).
of the endoscope, the operator can view the inside of an organ while controlling the viewing position and angle of the probe. During this invasive and interactive exploration, the endoscopist has full control of navigation within the hollow organ. Physical endoscopy applies the optical probe; however, virtual endoscopy is a new method of diagnosis using computer processing of three-dimensional (3D) image data sets. These are computed tomography (CT) and magnetic resonance imaging (MRI) scans to provide simulated visualizations [1,2] of specific organs of a patient. Visualization is the process of exploring, transforming, and viewing data as images to gain understanding and insight into the data. The visual effects of virtual endoscopy are similar or equivalent to those produced by standard endoscopic procedures [3]. Conventional CT and MRI scans produce axial slices of the
1350-4533/$ – see front matter © 2005 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2005.07.011
340
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
body that are viewed sequentially by radiologists who must imagine or extrapolate from these views what the actual 3D anatomy should be. By using sophisticated algorithms and high performance computing, these axial slices may be rendered as direct 3D representations of human anatomy. Specific anatomic data appropriate for realistic endoscopic simulations can be obtained from 3D MRI digital imaging examinations [4] or 3D acquired spiral CT data [5]. In fact, virtual endoscopy provides viewing control and options that are not possible with real endoscopy such as direction, angle of view, scale of view, immediate translocation to new views, lighting, and measurement. In addition, there are many body regions not accessible by optical endoscopy but can be explored with virtual endoscopy. Several important body systems are compatible with invasive probes including the heart, spinal canal, inner ear, biliary, pancreatic ducts, and large blood vessels. These are important anatomic structures ideally suited for virtual endoscopy [6]. By contrast, optical endoscopic procedures can be uncomfortable and sedation or anesthesia may be required [7]. Furthermore, optical endoscopic procedures are invasive and have serious side effects such as perforation, infection and hemorrhage. Therefore, virtual endoscopic visualization avoids the risks associated with real endoscopy, and when used prior to performing an actual endoscopic examination can minimize procedural difficulties and decrease the rate of morbidity, especially for endoscopists in training.
2. Related work The navigation paths for virtual camera play a crucial role in virtual endoscopy application. The navigation technique decides how the physician controls the position and orientation of a virtual camera to examine the inner surface of organ. The planned navigation scheme pre-computes a flythrough path to automatically move the virtual camera from one end of the organ to another end. On the other hand, the free navigation scheme aims to allow a physician to interactively control the virtual camera position and orientation. The physician is required to control the virtual camera at each frame to navigate through the virtual model. Although the physician is specifically well trained, it can be difficult to navigate to target. This may cause the problem that the camera position and orientation are lost in the virtual environment. The guided navigation scheme has combined the benefits of both planned and interactive navigation. By using a guided navigation scheme, the virtual camera moves automatically from a source point to a target point along a predefined flythrough path. When necessary, the physician can take over the control at any time by interactively and easily adjusting the camera position and direction. The switch scheme between these two modes is seamless, resulting in a smooth navigation [8]. For the purpose of guided navigation scheme, the determination of fly-through path is an essential and prerequisite
task. The well-defined flight paths provide a stable viewing direction, avoiding sudden twists and turns. The manually planning flight path such as key framing [2] is a tedious and time-consuming task. The key framing is that the user specifies manually the pose of a small subset of the total number of frames to be rendered. Then, interpolating curves are fitted to generate a continuous path and set of frames. The other method called distance mapping, a technique used in robot path planning, has also been applied to navigation path planning for virtual camera. A goal voxel is selected within the volume and a distance to the goal is computed for each voxel within the structure of interest. The distance assigned is the shortest distance along a path to the goal through the structure. However, the shortest distance approach has a tendency to produce paths that hug the wall of the organ rather than follow the central axis, thus limiting what is visible to the virtual camera. Although the improvement has been conducted by penalized distance [9], the distance mapping still requires substantial computing resources. On the other hand, the medial axis is the locus of the centers of all circles that are tangent to the boundary of the object at two or more disjoint points. Based on the medial axis transform (MAT) [10], continuously generated centerline has been used as a fly-through path for virtual camera [11]. However, MAT is heavily dependent on the circular shape of object on 2D image space. As a result of unstable viewing direction and abrupt twists, operator interaction is required to adjust the only jagged navigation path. In addition, the skeleton of a volumetric data provides a compact description of interesting object. The skeletonization also provides an efficient method for visualization and analysis such as feature extraction, feature tracking, and automatic navigation [12]. The topological thinning or onion peeling [13,14] is one of the skeletonization methods that one layer at a time of voxels is peeled off the object until just the skeleton remains. The skeleton algorithm automatically provides a 3D central axis of structure; however, onion peeling is computationally expensive. In this paper, therefore, we propose a fast and automated navigation path generation algorithm for virtual angioscopy based on two-phase adaptive region-growing and bounding box method. In the following, Section 3 describes the two-phase adaptive region-growing algorithm for carotid artery segmentation and then Section 4 provides a navigation path generation method for virtual angioscopy. Section 5 describes the mathematical phantom for validation of segmentation and navigation path generation. In Section 6, experimental results of segmentation and navigation path generation on both phantom and patient data sets are presented. Finally, Section 7 provides the discussions and concluding remarks.
3. Two-phase adaptive region growing The 3D time-of-flight (TOF) carotid artery (CA) MRA images are used to automatically generate the navigation
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
path for virtual angioscopy based on two-phase adaptive region-growing segmentation algorithm. The carotid artery is segmented as the following sequence: (1) image analysis, (2) initial thresholding for seed selection, (3) 2D region growing, (4) bifurcation detection, and (5) 3D region growing. Fig. 1 shows the overall segmentation process of carotid artery.
341
3.1. Image analysis The pairs of carotid and vertebral artery in the neck mainly deliver the blood to the brain about 80%, 20% of blood flow rate, respectively. Anatomically, the carotid artery is classified into common carotid artery (CCA), internal carotid artery (ICA), and external carotid artery (ECA). The CCA is branched off ICA and ECA at the position of thyroid cartilage and then connected to Circle of Willis through anterior and posterior communication artery. In this paper, the image analysis is defined as preliminary tasks to automatically segment the carotid artery. These consist of two image examination sub-tasks, which are: (1) the gray level distribution between objects and background and (2) the anatomic structure of carotid artery. Fig. 2(a) shows the pairs of CCA and vertebral artery (VA) before branch; Fig. 2(b) shows the ICA and ECA just after bifurcation from CCA. As shown in Fig. 2, CA and VA are brighter than surrounding tissues and background. This is the nature of MRA scanning protocol that focus on the blood flow area and suppresses the signal intensity of surrounding tissues. Fig. 3 shows the histogram of Fig. 2 that is plotted the frequency of intensity occurrence in the image [15–17]. Compared to background area and surrounding tissues in Fig. 2, blood flow area is small portion of the image. As indicated in the histogram, the frequency of gray level occurrence is biased to the left side. Consequently, the dominant background pixels provide uni-modal shape of histogram. While global and local scheme of thresholding is difficult to segment the object or region, region growing is suitable for the delineation of small structures such as tumors [18], lesions, and blood vessels in the medical imagery. Regarding to the anatomic structure of carotid artery, a priori knowledge is used as the following: (1) spatial information between VA and CA for seed selection of region growing, (2) spatial information and branch direction of ICA and ECA for navigation path generation. 3.2. Initial thresholding for seed selection
Fig. 1. Block diagram of carotid artery segmentation process.
Image segmentation is the process that partitions a digital image into disjoint connected sets of pixels, each of which corresponds to an object or region. The various image processing algorithms could be combined to segment the objects on the image space. In the case of pipeline operation, thresholding technique could be used as an initial or intermediate processing step to completely segment the specific objects or regions. Thresholding is a pixel classification process to identify the pixels of a given image into two classes: those pertaining to objects and those pertaining to background. While one class includes pixels with gray values that are below or equal to a certain threshold value, other class includes those with gray values above the threshold. Thresholding is the transformation of an input image f to an output binary image g as follows:
342
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
Fig. 2. Carotid artery MRA axial image: (a) before and (b) after bifurcation.
g(x, y) =
1
for f (x, y) ≥ T
0
for f (x, y) ≺ T
(1)
where T is the threshold, g(x,y) = 1 for image elements of objects, and g(x,y) = 0 for image elements for the background. In order to determine the region-growing seed in the sequence images, the initial segmentation is performed with threshold value T = (2h − 1)/2, where h is the number of bits for representation of each pixel. Fig. 4 shows the initial segmentation result that is the binary version of Fig. 2(a). The CA is isometric and almost symmetrically resides in left and right side of the neck, even though anatomic structures are varied among the subjects. By applying a priori knowledge of CA anatomic structure, the region R is divided into two separate sub-regions to locate two initial seeds on R1 and R2 for pairs of CA segmentation, respectively. The sub-region R1 and R2 are as follows: R1 = {(x, y), 1 ≤ x ≤ xm , 1 ≤ y ≤ yn /2}
(2)
R2 = {(x, y), 1 ≤ x ≤ xm , yn /2 + 1 ≤ y ≤ yn }
(3)
although the (row, column) orientation used in matrices is also often used in digital image processing. These sub-regions R1 and R2 are used to recognize the left and right carotid artery for region-growing seeds, respectively. In consideration of the spatial relationship between VA and CA as shown in Fig. 4, the every pixel, which have value of g1 (x,y) = 1, is examined from the top left corner toward the horizontal direction on each sub-region. By using eight connectivity features, the connected neighboring pixels are labeled as carotid artery starting at the first located pixel that have value of g1 (x,y) = 1. The labeled component of each sub-region is the binary version of carotid artery segmentation result in the first sequence image that can be used as seed region for further region growing. 3.3. 2D region growing The image can be modeled by a continuous function of two or three variables; in the 2D case arguments are coordinates (x,y) in a plane, while if images are expended to 3D a third variable might be added. The image sequence function F is as follows:
where xm , yn represent maximal image coordinates. The customary orientation of coordinates in an image is in the Cartesian fashion (horizontal x-axis, vertical y-axis),
Fig. 3. Histogram of carotid artery MRA axial image (Fig. 2).
Fig. 4. Binary version image (Fig. 2(a)).
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
343
Fig. 5. Candidate requirements of start pixel selection for region growing: (a) segmented image (current slice), (b) unsegmented image that have smaller region (next slice) than segmented current slice image (a), and (c) unsegmented image that have larger region (next slice) than segmented current slice image (a).
F=
n
fi (x, y)
(4)
i=1
where n is image sequence number, x and y are the spatial coordinates. CT and MRI scans produce series of axial slice images. These volumetric data could be used as input image sequence for segmentation of specific object or region. Based on Eqs. (1) and (4), the segmented output image g1 is a binary form of the input image f1 . In order to apply the region growing on the first sequence image, the intensity values are averaged only for CA region in image f1 , which is corresponding region of labeled component in image g1 . This averaged value is used as homogeneity criterion for region growing on the first sequence image f1 starting at arbitrary selected pixel with eight-connectivity feature. After the region growing is finished on the first sequence image f1 , homogeneity criterion is adjusted for segmentation of next sequence image f2 . Based on the connectivity-preserving features of sequence image, region growing is extended to the next slice image with adaptive homogeneity criterion that is re-calculated in the current slice image. In order to select the start pixel of region growing for next slice segmentation, the pixel must meet the following criteria. These conditions are: (1) pixel must be labeled in the current slice, (2) pixel must not be labeled in the next slice, (3) coordinates must be same on the current and next slice, (4) pixel of next slice must be passed the homogeneity test, and (5) next slice must not be last slice. Fig. 5 shows the candidate requirements of start pixel selection for next slice image segmentation. These example images (10 × 10) provide detail explanation of top three conditions that are described above to select the start pixel for region growing. Fig. 5(a) shows the segmentation result on current slice image. Fig. 5(b) shows the un-segmented next slice image that have smaller region; Fig. 5(c) shows the un-segmented image that have larger region than segmented current slice image. After start pixel is selected, pointer for current slice number is incremented and region growing is performed on next slice image. The 2D region-growing process is performed until next slice number is not last slice image. 3.4. Bifurcation detection After segmentation is finished on sequence image f2 , the face connection test is performed to detect bifurcation
between current slice image f2 and previous slice image f1 . In order to find the branch in the current slice, every pixel in the current slice must meet the following criteria. These conditions are: (1) pixel must be labeled in the previous slice, (2) pixel must not be labeled in the current slice, (3) coordinates must be same on the previous and current slice, and (4) pixel of current slice must be passed the homogeneity test. If the bifurcation is detected in the current slice image, the current slice number and homogeneity criterion are saved into the stack S and retrieved from the stack S to segment the other branch (ICA or ECA). In addition, stack S3D is duplicated from the stack S for second-phase 3D region growing. The stack S and S3D are identical that have slice number and homogeneity criterion in the case of bifurcation is detected. The stack S3D is only used to save in first-phase region growing and to retrieve the slice number and homogeneity criterion in second-phase region growing. The branch detection process is performed until next slice number is not last slice image. 3.5. 3D region growing The region growing is a procedure that group pixels or sub-regions into larger regions based on predefined criteria. The basic approach is to start with a set of seed points and from these grow regions by appending to each seed those neighboring pixels that have properties similar to the seed, such as specific ranges of gray level or color [19–21]. In similar, region growing could be expanded to 3D image space with voxel connectivity. The 3D region-growing method perform the homogeneous test from the start voxel (or 3D block) to the neighbor voxel (or 3D block) using graylevel, texture, color as an acceptance criterion. According to the homogeneous test result, the neighbor voxel (or 3D block) is included or excluded until termination condition is satisfied. In addition, the voxel connectivity is defined as the following: (1) six-connectivity with joint faces, (2) 18connectivity with joint edges, and (3) 26-connectivity with joint corners. If the next slice is last slice image and if stack is not empty, slice number and homogeneity criterion are retrieved from the stack S for the other branch segmentation. The 2D region-growing process is performed until next slice is last slice image and if stack is empty. After first-phase 2D adaptive region growing is terminated including branch
344
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
Fig. 7. (a)–(e) Five models of circumscribed quadrangle.
The circumscribed quadrangle is connected to the four tangent lines that run parallel with x- and y-axis in the image, respectively. The four tangent points are located outermost in the segmented carotid artery from left, right, top, and bottom directions, respectively. Fig. 7 shows the five models of circumscribed quadrangle including one circular and four ellipse shapes. The dashed lines represent the bisector of each circumscribe quadrangle. The generated centerline on each axial slice of carotid artery is used as a fly-through path for virtual angioscopy. 4.2. Rendering of endoscopic view
Fig. 6. Two-phase adaptive region-growing algorithm.
segmentation, the second-phase 3D region growing is performed as the following. These process are: (1) the intensity values of labeled pixels, which are from 2D region growing on whole sequence image, are averaged and used as an homogeneity criterion, (2) 26-connectivity feature is used, (3) only apply to the unlabeled voxel for preserving the labeled pixel in the first-phase region growing, and (4) use the duplicated stack S3D for branch handling. To summarize the carotid artery segmentation process, the algorithm for the two-phase adaptive region growing is specified in Fig. 6.
4. Navigation path generation 4.1. Navigation path By using a priori knowledge such as spatial information and branch direction of carotid artery, common carotid artery and internal carotid artery are segmented to determine the navigation path for virtual camera before and after the bifurcation, respectively. In order to visualize the inside of internal carotid artery, a fly-through path is pre-computed to automatically move the virtual camera from common carotid artery to internal carotid artery. Our method for positioning the virtual camera is based upon the intersecting coordinates of two bisectors on the circumscribed quadrangle of segmented carotid artery.
The endoscopic view of the internal carotid artery can be achieved by either surface or volume rendering computer graphic techniques. Surface rendering is an indirect method of obtaining an image from a volume data set. The surfaces are produced by mapping data values onto a set of geometric primitives in a process known as isosurfacing. These isosurfaces can then be rendered into a displayable image using standard computer graphics techniques. On the other hand, volume rendering is a more direct method for reconstruction of 3D structures. Volume rendering represents 3D objects as a collection of cube-like building blocks called voxels. The main advantage of this type of rendering is its ability to preserve the integrity of the original data throughout the visualization process. This technique, however, is generally expensive than conventional surface rendering technique. Although surface rendering lacks rendering quality and information beyond surface in terms of the usage of partial volume data, the surface rendering technique provides interactive rendering speed and usability of standard graphic hardware. Marching cubes is an algorithm for generating isosurfaces from volumetric data known as surface rendering technique. This algorithm uses the information at the corners of a voxel to construct a surface that approximates the original surface. The virtual angioscopy is implemented to provide the endoscopic view of internal carotid artery by using both Marching cubes algorithm [22] and volume rendering algorithm with perspective projection technique.
5. Phantom generation For the purpose of the performance quantification for segmentation method and specific medical application such as proposed automatic navigation path generation algorithm, validation experiments are necessary. Validation is typically
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
Fig. 8. Segmentation result of: (a) common carotid artery before bifurcation and (b) internal and external carotid artery after bifurcation.
Fig. 9. 3D reconstruction view of: (a) phantom and (b) clinical data set using segmentation results.
performed using one of two different types of truth models. The most straightforward approach to validation is by comparing the automated segmentations with manually obtained segmentations [23,24]. This approach, besides suffering from the drawbacks outlined of laborious and time consuming, does not guarantee a perfect truth model since an operator’s performance can also be flawed. The other common approach to validating segmentation methods is by use of physical phantom [25] or computational phantoms which consisting of digital (voxel-based) phantom [26,27] and mathematical (analytical) phantom [28]. Physical phantoms provide an accurate depiction of the image acquisition process but typically do not present a realistic representation of anatomy. Digital phantoms are mainly derived from segmented tomographic images of the human anatomy obtained by either CT or MRI. Mathematical phantoms consist of regularly shaped continuous objects defined by combinations of mathematical geometries such as spheres, ellipsoids, cylinders, and cones. Although anatomically less realistic than phantoms
345
Fig. 11. Projected centerline for: (a) phantom and (b) clinical data set.
derived from CT or MR images of patients, the mathematical phantom has the advantage that it can be easily modified to simulate a wide variety of patient anatomies. By considering the anatomical structure of blood vessel, the cylinder shape of mathematical phantom is devised to simulate the carotid artery bifurcation. A carotid artery phantom is implemented that consist of two cylinders model without branches (simulation of vertebral artery) and two cylinders model with bifurcation for each cylinder (simulation of carotid artery). The ellipse objects are drawn on 2D image space to simulate the blood vessels with a graphic package. To simulate the tomographic image sequence, the ellipse objects are located at the corresponding coordinates that are determined on clinical data set. For the purpose of intensity value assignment to blood vessels and background, the inside and outside of each ellipse objects are recognized on 2D image space. Based on the intensity assessment, the initial threshold value T = (2h − 1)/2, where h is the number of bits for representation of each pixel, is used to assign the intensity value for blood vessel and background. The intensity above the threshold value T is assigned to blood vessels and below the threshold value T is assigned to background based upon the examination result of intensity distribution between objects and background in the clinical data set. A random number generator also used to determine the specific intensity value both blood vessels and background after range are checked with threshold value T. Regarding to the determination of blood vessel size, the diameter of common carotid artery and internal carotid artery are designed to same radius; external carotid artery and vertebral artery are also implemented with same radius. According to the geometrical structure of carotid artery, isometrically simulated carotid and vertebral artery are located in the sub-regions R1 and R2 , respectively. The branch positions determined on clinical data
Fig. 10. (a)–(f) Center coordinates of segmented carotid artery based upon intersecting coordinates of two bisectors on circumscribed quadrangle.
346
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
Fig. 12. Endoscopic view of carotid artery where virtual camera was positioned at: (a)–(c) common carotid artery and (d)–(f) internal carotid artery by using ray casting algorithm.
set are also used to implement the bifurcations of the carotid artery phantom.
6. Experimental results The Intel P4 processor with MS Windows, Visual C++ and Matlab were used to segment the carotid artery, to generate the navigation path, and to implement the virtual angioscopy system. Our algorithm was applied on six clinical data sets of 3D time-of-flight (TOF) MRA images from GE Signa Horizon Echospeed scanner (TR = 25 ms, TE = 6.9 ms, FA = 25◦ , FOV = 220 mm). Each case consisted of 136 axial images and the dimension of image was 512 × 512 pixels, pixel spacing was 0.43 mm and slice thickness was 1.4 mm. The phantom data set was also used to validate our algorithm. Fig. 8(a) and (b) shows the segmentation result of CCA before bifurcation and ECA and ICA after bifurcation by using two-phase adaptive regiongrowing algorithm, respectively. Fig. 9(a) and (b) presents the fully rotatable external view of 3D reconstructed phantom and clinical data set using segmentation results. As a result of navigation path generation, Fig. 10 shows the center coordinates of segmented carotid artery based upon intersecting coordinates of two bisectors on the circumscribed quadrangle. Fig. 11 also shows the projected centerlines of phantom and clinical data set by using maximum intensity projection algorithm. The times to compute the navigation paths for 12 carotid arteries in six clinical data sets ranged from 19 to 21 s
depending on branch location of each carotid artery including segmentation times. For the purpose of visualization of carotid artery, Marching cubes algorithm was used to generate the endoscopic view. This algorithm provided interactive rendering speed at 15 frames per second, so effective navigation of carotid artery was possible. On the other hand, ray casting is one of the volume rendering techniques using backward projection. Rays are cast from each pixel of the image plane into the volume data. At locations along each ray, a sample value and a surface normal approximation are calculated using values of surrounding voxels. Using the sample value and normal, a sample opacity and color are dynamically assigned by a lookup table or in a preprocessing phase. Then a local shading model is applied and the samples along the ray are composed into a pixel value of the final image [29,30]. Fig. 12(a)–(c) and (d)–(f) shows the endoscopic view of carotid artery where the virtual camera was positioned at common carotid artery and internal carotid artery by using ray casting algorithm, respectively. This algorithm required lots of computing time, so only provided as an off-line navigation of internal carotid artery.
7. Discussion and conclusions Based on two-phase adaptive region-growing algorithm, we have presented a fast and automated navigation path generation method to visualize the carotid artery. This method
D.-Y. Kim et al. / Medical Engineering & Physics 28 (2006) 339–347
was validated by phantom data set and applied on clinical data sets. This algorithm is more effective than key framing and topological thinning method in terms of automated features and computing time. This algorithm is also applicable to generate the centerline of cylinder shape structures and tree-like blood vessels such as renal artery and coronary artery. The limitation of this algorithm occurs on occlusion or near occlusion cases resulted in segmentation failure of carotid artery. Caused by further problem on navigation path generation, it is impossible to visualize the inside of carotid artery. In addition, this algorithm is applicable to convex shape of blood vessel. The central axis of carotid artery, which was determined by using two-phase region growing and bounding box method, could be used as navigation path for virtual camera. In the case of that blood vessels run horizontally across the scan plane, however, the images provide non-circular shape of blood vessel. With this limitation, the navigation path cannot be easily determined by proposed algorithm. Hence, the further works are required to generate the oblique slices which scanning planes are perpendicular to the central axis of blood vessel. These slices provide the circular shape of blood vessel even though artery runs horizontally across the scan plane. Furthermore, it is beneficial that circular objects are restored to original shape, especially for distorted objects caused by stenosis, scanning artifacts and possible noise. In conclusion, newly developed virtual angioscopy system aims at supporting physicians in determining the precise spatial location and shape of the carotid artery stenosis that are the critical indicators of the transient ischemic stroke. Combining with a computerized quantification method for carotid artery stenosis [31], virtual angioscopy system will provide better diagnostic results of carotid artery disease.
References [1] Geiger B, Kikinis R. Simulation of endoscopy. In: Proceedings of the AAAI spring symposium series: applications of computer vision in medical image processing. 1994. p. 138–40. [2] Rubin GD, Beaulieu CF, Argiro V, Ringl H, Norbash AM, Feller JF, et al. Perspective volume rendering of CT and MR images: applications for endoscopic imaging. Radiology 1996;199:321–30. [3] Wickham JEA. Minimally invasive surgery: feature developments. Br Med J 1994;308:193–6. [4] Rusinek H, Mourino MR, Firooznia H, Weinreb JC, Chase NE. Volumetric rendering of MR images. Radiology 1989;171:269–72. [5] Fishman EK, Jeffrey RB. Spiral CT: principles, techniques, and applications. Philadelphia: Lippincott, Williams & Wilkins; 1998. [6] Robb RA. Virtual (computed) endoscopy: development and evaluation using the visible human datasets. In: Proceedings of the visible human project conference. 1996. p. 19–20. [7] Jolesz FA, Lorensen WE, Shinmoto H, Atsumi H, Nakajima S, kavanaugh P, et al. Interactive virtual endoscopy. Am J Roentgenol 1997;169(5):1229–35.
347
[8] He T, Hong L, Chen D, Liang Z. Reliable path for virtual endoscopy: ensuring complete examination of human organs. IEEE Trans Visual Comput Graph 2001;7(4):333–42. [9] Bitter I, Kaufman AE, Sato M. Penalized distance volumetric skeleton algorithm. IEEE Trans Visual and Comput Graph 2001;7(3):195–206. [10] Jain R, Kasturi R, Schunck BG. Machine vision. New York: McGraw-Hill; 1995. [11] Kim DY, Park JW. Visualization of the internal carotid artery using MRA images. Magn Reson Imag 2005;23(1):27–33. [12] Zhou Y, Toga AW. Efficient skeletonization of volumetric objects. IEEE Trans Visual Comput Graph 1999;5(3):196–209. [13] Hong L, Muraki S, Kaufman AE, Bartz D, He T. Virtual voyage: interactive navigation in the human colon. In: Proceeding of ACM SIGGRAPH. 1997. p. 27–34. [14] Paik DS, Beaulieu CF, Jeffery RB, Rubin GD, Napel S. Automated flight path planning for virtual endoscopy. Med Phys 1998;25(5):629–37. [15] Gonzalez RC, Woods RE. Digital image processing. New York: Addision-Wesley; 2000. [16] Parker JR. Algorithm for image processing and computer vision. New York: John Wiley & Sons; 1997. [17] Calotto MJ. Histogram analysis using a scale-space approach. IEEE Trans Pattern Anal Mach Intell 1987;9:121–9. [18] Kim DY, Park JW. Computer aided detection of kidney tumor on abdominal computed tomography scans. Acta Radiol 2004;45(7):791–5. [19] Sonka M, Hlavac V, Boyle R. Image processing, analysis, and machine vision. Pacific Grove: PWS Publishing; 1999. [20] Adams R, Bischof L. Seeded region growing. IEEE Trans Pattern Anal Mach Intell 1994;16:641–7. [21] Mehnert A, Jackway P. An improved seeded region growing algorithm. Pattern Recogn Lett 1997;18:1065–71. [22] Lorensen W, Cline H. Marching cubes: a high resolution 3D surface construction algorithm. In: Proceeding of ACM SIGGRAPH. 1997. p. 163–9. [23] Vaidyanathan M, Clarke LP, Hall LO, Heidtman C, Velthuizen R, Gosche K, et al. Monitoring brain tumor response to therapy using MRI segmentation. Magn Reson Imag 1997;15:323– 34. [24] Wust P, Gellermann J, Beier J, Wegner S, Tilly W, Troger J, et al. Evaluation of segmentation algorithms for generation of patient models in radiofrequency hyperthermia. Phys Med Biol 1998;43:3295–307. [25] Lei T, Sewchand W. Statistical approach to X-ray CT imaging and its applications in image analysis. IEEE Trans Med Imag 1992;11: 62–9. [26] Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, et al. Design and construction of a realistic digital brain phantom. IEEE Trans Med Imag 1998;17:463–8. [27] Zubal IG, Harrell CR, Smith EO, Rattner Z, Gindi G, Hoffer PB. Computerized 3-dimensional segmented human anatomy. Med Phys 1994;21:299–302. [28] Zaidi H. Relevance of accurate Monte Carlo modeling in nuclear medical imaging. Med Phys 1999;26:574–608. [29] Watt A. 3D computer graphics. Boston: Addison-Wesley; 1993. [30] Noordmans HJ, Voort H, Smeulders A. Spectral volume rendering. IEEE Trans Visual Comput Graph 2000;6(3):196–207. [31] Kim DY, Park JW. Computerized quantification of carotid artery stenosis using MRA axial images. Magn Reson Imag 2004;22(3):353–9.