- Email: [email protected]
Digital subtraction CT angiography based on efficient 3D registration and refinement
Computerized Medical Imaging and Graphics 28 (2004) 391–400 www.elsevier.com/locate/compmedimag
Digital subtraction CT angiography based on efficient 3D registration and refinement Sung Min Kwona, Yong Sun Kima, Tae-Sung Kimb, Jong Beom Raa,* a
Department of Electrical Engineering and Computer Science, Korea Advanced Institute of Science and Technology, 373-1, Guseongdong, Yuseonggu, Daejeon 305-701, South Korea b Department of Diagnostic Radiology, College of Medicine, Yonsei University, Seoul, South Korea Received 30 January 2004; accepted 1 June 2004
Abstract A novel method for fast, automatic 3D digital subtraction CT angiography (DS-CTA) is presented to generate artifact-free angiograms. The proposed method consists of two steps: 3D registration to align a CT image to the CT angiography (CTA) image and subtraction-andrefinement to extract blood vessels only. For efficient and accurate 3D registration in the first step, an normalized mutual information (NMI) based algorithm is adopted, and its fast version is developed by introducing a new measure. To further improve the subtracted image quality in the second step, a novel 3D refinement algorithm is suggested to effectively remove unwanted residuals. Experimental results of seven clinical CT/CTA head datasets demonstrate that cerebral vessels are well extracted from CTA images with almost no loss. The typical processing time is 3–9 min depending on the image size in a PC with a 2.4 GHz CPU. q 2004 Elsevier Ltd. All rights reserved. Keywords: Digital subtraction; Angiography; CTA; DS-CTA; Registration; Mutual information; Refinement; Mask; Cerebral vessels
1. Introduction Recently, angiographic investigation using highresolution multi-detector (MD) CT images has become a useful tool for the sensitive assessment of vascular structures. In the head, for example, most of blood supplied to brain passes the special region called the Circle of Willis. Hence, its branching portions are continuously stressed and may undergo aneurismal change. Consequently, they are apt to be ruptured and a devastating subarachnoid hemorrhage may occur [1]; therefore, the accurate diagnosis of cerebral arterial aneurysm is critical. Since a CTA image describes blood vessels with high intensity, it is known to be very helpful for diagnosing this kind of vessel disease; however, since blood vessel parts of interest are often located near the high intensity skull base, general CTA with bones may not be quite helpful for accurate diagnosis, so a sophisticated
* Corresponding author. Tel.: C82-42-869-3434; fax: C82-42-8698360. E-mail address: [email protected] (J.B. Ra). 0895-6111/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2004.06.003
procedure to separate only blood vessels from bone is needed. Several methods to extract blood vessels from CTA images have been tried, based on segmentation or tracking [2,3]; however, those methods tend to leave some artifacts of bone near vessels and lose thin vessels. Meanwhile, a couple of vessel extraction techniques have been attempted on the basis of the CT and CTA image registration, namely, the subtraction of the CT image from a CTA image [4,5] and bone-elimination by replacing the intensities of bone with a soft tissue intensity [6]. In Ref. [4], during CT/CTA acquisition, a patient’s head is immobilized in a headholder by using forehead and chin straps. In this case, the patient feels uncomfortable during the data acquisition and slight motion due to imperfect fixation is a strong limitation to this subtraction technique. A feature-based registration for correcting patient motion between CT/CTA images releases the patient from the discomfort of the head-holder [5]; however, a subtle mis-registration of the two images may cause some bone fractions near the bone boundary in the subtracted image. Therefore, for a satisfactory vessel extraction, an accurate, robust 3D registration is
392
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
indispensable. Maximization of mutual information (MI) has been recommended as a powerful intrinsic criterion for accurate multi-modal or mono-modal image registration [7–10]. However, MI requires long processing time, especially for a recent large size 3D volume dataset obtained from a high resolution CT. Practically, there are many urgent cases in the diagnosis using DS-CTA, such as rupturing cerebral aneurysms; therefore, the importance of fast 3D registration cannot be overestimated, and its processing time should be improved. Even with the accurate registration between CT and CTA images by using an advanced scheme, their subtraction and bone-elimination results still contain unwanted residuals at boundaries between bone and tissue, which is due to and the limited resolution of CT and CTA images, and these residuals cause annoying artifacts in the angiogram. In this paper, we develop a novel DS-CTA scheme by adopting an accurate NMI-based registration technique. By solving the two problems in the scheme, namely, the calculation complexity of the NMI-based registration and the residual artifact due to subtle mis-registration, we propose a fast automatic 3D DS-CTA method based on CT and CTA images. The scheme is easily applicable to a PC environment. In Section 2, the proposed fast registration method is described, and in Section 3, the 3D refinement process to remove unwanted artifacts is described in detail. Experimental results are presented in Section 4, and the conclusions are drawn in Section 5.
2. The proposed effective 3D registration based on NMI The CT and CTA head images acquired from a patient include the same information of bone and soft-tissue with a certain mis-alignment, while the blood vessel information is included only in the CTA image. Therefore, the blood vessel information can be extracted using CT and CTA images simultaneously. Fig. 1 shows the overall block diagram of
Fig. 1. Overall block diagram of the adopted 3D DS-CTA algorithm.
Fig. 2. Block diagram of an NMI-based registration algorithm.
the 3D DS-CTA algorithm we adopt. As shown in the figure, the algorithm first registers a CT image to the CTA image, subtracts it from the CTA image, and finally refines the subtracted image by removing unwanted residuals or artifacts. As mentioned in Section 1, the precise registration is very important for a satisfactory subtraction result. Hence, we adopt the NMI-based scheme for the registration of CT and CTA images due to its superior performance, such as robust and overlap-invariant characteristics [10]. In the following sub-sections, we first briefly describe the NMIbased scheme and a multi-resolution approach, which we adopt to our algorithm for accurate and fast registration. Then, we propose a volume mask based registration scheme to further accelerate the processing speed. 2.1. Conventional NMI-based multi-resolution registration Fig. 2 shows the block diagram for a conventional NMIbased registration. Registration is a procedure to find the optimal transformation matrix T projecting a floating volume into the reference frame as shown in Fig. 3, i.e. a procedure to match the two volumes best using a similarity
Fig. 3. Registration procedure of a floating volume to the reference frame by transformation matrix T.
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
measure. NMI is an overlap-invariant entropy measure that correctly describes a dispersion of the joint histogram for two random variables and measures their statistical dependence. Hence, NMI has been extensively investigated as a similarity measure [7–16] and has been proven to perform well in 3D medical image alignment [10]. NMI is defined as YðA; BÞ Z
HðAÞ C HðBÞ HðA; BÞ
complexity of multi-resolution registration is represented as CC Z
LK1 X ð the number of trialsÞl !ðcomputational lZ0
complexity of a similarity measureÞl f
LK1 X
ð of trialsÞl
lZ0
!ðthe number of sub-sampled floating (1)
where H(A) and H(B) are the marginal entropy of random variable A and B, respectively, and H(A, B) is the joint entropy of A and B. These three entropies are calculated from the joint probability distribution of A and B using the Shannon measure of entropy, and the joint probability distribution can be estimated from the joint histogram as shown in Fig. 4; therefore, the NMI is obtained if the joint histogram is available. During the registration procedure in Fig. 2, the elements of a transformation matrix T are updated to maximize the NMI in each trial. Accordingly, the corresponding projected position of the floating volume is also refined. In order to generate a joint histogram based on intensity pairs in each trial (see Fig. 4), we need intensity values at the reference volume positions corresponding to the projected positions of all the voxels in the floating volume. In obtaining those intensity values, tri-linear interpolation is to be performed in the reference volume. This extensive tri-linear interpolation processing for all floating voxels in each trial requires massive computation, which is the major reason for the long processing time in NMI-based registration. To overcome this problem, a multi-resolution approach is adopted by using a hierarchical structure of input data [11–13]. For multi-resolution image registration, a floating volume is sub-sampled with an interval of 2LKlK1 voxels at the lth level (lZ0,1,.,LK1). Then, based on the similarity measure that is obtained by using the subsampled floating volume at each level, the transformation matrix T is searched and refined starting from the initial level. As the level increases, the sub-sampling interval decreases and the registration result becomes more accurate with higher computational complexity. The computational
Fig. 4. Joint histogram of intensity pairs.
393
volume voxelsÞl (2) Even though this multi-resolution approach can considerably reduce the computational complexity, we still need to reduce it. In Eq. (2), the majority of the computation results from a calculation of a similarity measure at the final level because the number of sub-sampled floating volume voxels increases exponentially as the level becomes higher. Hence we will focus the reduction of the number of voxels to be used in the calculation of the similarity measure at the higher levels. 2.2. The proposed fast registration A fast 3D registration algorithm based on the multiresolution scheme is given in Fig. 5. Here, we introduce a floating volume mask to the conventional NMI-based multiresolution registration scheme described above so that a reduced number of voxels are used to evaluate NMI in the repeated trials to find the best transformation matrix at each level. Fig. 6 demonstrates a mask setting procedure at level l.
Fig. 5. Block diagram of the proposed 3D multi-resolution registration algorithm.
394
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
the registration procedure. Note here that Thfv is properly decided according to the bin size of the joint histogram.
3. 3D refinement of the subtracted image
Fig. 6. The procedure to obtain a floating volume mask at resolution level l.
As shown in Fig. 1, the aligned CT image is subtracted from the CTA image to extract blood vessels. Even though the alignment result by registration is fairly good, the subtracted image still contains unwanted residuals from a subtle mis-alignment at bone and air boundaries due to abrupt and considerable intensity changes; moreover, the residuals hinder proper examining blood vessels because their intensity values are still comparable with blood vessel intensity values. To get rid of these unwanted residuals, we introduce and apply a masking procedure after
This procedure begins with checking the convergence rate (CR) of the previous level, CRlK1, which is defined as follows vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 8 u1 X CRlK1 Z t ! jðtransformed position of vi with T1;lK1 Þ K ðtransformed position of vi with TK;lK1 Þj2 8 iZ1 where vi represents eight corner voxels of a floating volume, and T1,lK1 and TK,lK1 denote the initially and finally obtained transformation matrices at level lK1, respectively. If CRlK1, is smaller than the minimum voxel-interval dmin among the ones along the x-, y-, and z-axes, we set a floating volume mask Mfv($) that consists of the voxels having an MIDV (maximum intensity difference of a voxel) larger than a threshold. Otherwise, the whole floating volume is set to the mask. Therefore, at higher levels where the number of subsampled voxels increases and the convergence rate becomes slow, we can partially select the floating volume voxels, reducing the calculation time for a similarity measure, NMI. Now, we consider a measure MIDV for the proper selection of a mask. As described in Section 2.1, during the registration procedure, the elements of T are refined to the direction increasing the NMI value. If the intensity change of an intensity pair for a voxel between two consecutive trials is too little to affect the joint histogram, the corresponding voxel can be considered as not contributing to any change to the value of NMI; therefore, those voxels may be removed from the registration procedure without degrading the registration accuracy. To discriminate those voxels, we introduce a measure MIDV for each floating volume voxel, which is defined as dmin MIDV Z max jIðvi Þ K Iðvc Þj ! ; i Z 1; .; 26 : (4) d vi Here, vi represents 26 neighboring voxels of a center voxel vc, dvi denotes the distance between vi and vc, and I($) denotes the intensity of a voxel. Note that MIDV represents the degree of intensity change from a voxel vc. If the MIDV of a floating volume voxel is smaller than a threshold value Thfv, the voxel is excluded from a mask, thereby from
(3)
the subtraction, i.e. Vðx; y; zÞ Z Mba ðx; y; zÞ !Sðx; y; zÞ; where
(
Sðx; y; zÞ Z
(5)
f ðx; y; zÞ K gðx; y; zÞ
if f ðx; y; zÞO gðx; y; zÞ;
0
otherwise:
(6) Here, V(x,y,z), Mba(x,y,z), S(x,y,z), f(x,y,z), and g(x,y,z) represent the final vessel image, mask, subtraction image, CTA image, and aligned CT image, respectively. Since unwanted residuals as well as blood vessels are contained in S(x,y,z), S(x,y,z) is the target image for masking operation. The purpose of masking is to effectively remove unwanted residuals by using a mask function that represents definite bone and air areas excluding vessels. The procedure to obtain the mask, Mba(x,y,z) is shown in Fig. 7. We first consider the bone mask Mb that represents definite bone areas, and then an initial bone mask is generated in the aligned CT image by properly selecting a threshold value of bone, Thb, based on the intensity histogram of CT image, i.e. ( 1 if gðx; y; zÞO Thb ; Mb ðx; y; zÞ Z (7) 0 otherwise:
Fig. 7. The procedure to obtain the final mask Mba.
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
395
Table 1 Dimension and resolution of test image sets CT/CTA datasets
Dimension
Resolution (mm3)
Phantom Patient_1 Patient_2 Patient_3 Patient_4 Patient_5 Patient_6 Patient_7
512!512!180 512!512!185 512!512!185 512!512!163 512!512!215 512!512!147 512!512!201 512!512!180
0.300!0.300!0.300 0.268!0.268!0.300 0.268!0.268!0.300 0.326!0.326!0.500 0.301!0.301!0.300 0.289!0.289!0.500 0.275!0.275!0.300 0.254!0.254!0.300
All the datasets include CT and CTA images having the same dimension and resolution.
Note that the initial bone mask is made not in the CTA image but in the aligned CT image where the vessel intensity is much lower than that of bone; therefore, reliable bone areas can be selected as the initial bone mask. Since we assume the inherent subtle mis-alignment between registered CT and CTA images, however, the initial bone mask defined on the registered CT image is not suitable for completely removing unwanted bone residuals from the subtraction image. Hence, we refine the initial bone mask using the CTA image. Since the initial bone mask obtained through simple thresholding does not cover all the bone areas, the mask should be expanded until it covers all the bone areas so that the final bone mask may remove all the unwanted bone fractions from the subtraction image. It is obvious in the CTA image that the gradient is steep at the boundary between bone and tissue; hence, in the proposed algorithm, the bone mask is expanded voxel by voxel toward the outside of its boundary if the intensity value of the corresponding CTA voxel differs enough from any one of its 26 neighboring voxel intensity values. In other words, a background voxel adjacent to the mask boundary is included as part of the bone mask, if it is assured to belong to a transition region from bone to other tissue. The bone mask updating process can be represented as follows Mb ðx C i; y C j; z C kÞ Z 1; if Mb ðx; y; zÞ Z 1; & Mb ðx C i; y C j; z C kÞ Z 0;
(8)
& f ðx C i; y C j; z C kÞ! f ðx; y; zÞ K ThG ; where ThG denotes a threshold value for the intensity difference in the CTA. Here, i, j, kZK1, 0, 1 except for
Fig. 8. Volume rendered images of a phantom dataset. (a) CT and (b) CTA.
Fig. 9. Volume rendered images of a dataset, Patient_1. (a) CT and (b) CTA.
the case of iZjZkZ0. This procedure is repeated until the expansion stops naturally at the end of the bone boundary where the gradient becomes slow. It is interesting to note that in the case where a vessel is located near bone, the mask expansion based on Eq. (8) can prohibit the mask from intruding into the vessel because the gradient becomes slow at the boundary between bone and vessel. An initial air mask is also generated to get rid of unwanted residuals between soft tissue (or bone) and air. Since there is no vessel in direct contact with air and the air provides a fairly constant intensity value, the initial air mask is generated on the basis of the intensity value rather than the gradient value as in the bone mask. The initial air mask is defined as Ma ðx; y; zÞ Z 1 if gðx; y; zÞ! Tha ;
(9)
and the final air mask Ma is simply refined on the basis of CTA image intensities as follows Ma ðx C i; y C j; z C kÞ Z 1; if Ma ðx; y; zÞ Z 1 & Ma ðx C i; y C j; z C kÞ Z 0
(10)
& f ðx C i; y C j; z C kÞ! Tha ; where Tha is the threshold value based on the intensity histogram of CT image, and i, j, kZK1, 0, 1, except iZjZ kZ0. Then, the final mask is obtained by simply combining the bone and air masks, i.e. Mba ðx; y; zÞ Z ðMb ðx; y; zÞ C Ma ðx; y; zÞÞ:
(11)
This mask is used for excluding definite bone and air areas.
Fig. 10. Volume rendered images of a dataset, Patient_3. (a) CT and (b) CTA.
396
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
Fig. 11. Volume rendered images of a dataset, Patient_4. (a) CT and (b) CTA.
Fig. 12. Floating volume mask. (a) A slice image of CT data for floating volume of Patient_1 and (b) its selected mask with Thfv of 8 HU.
4. Experimental results and discussions The proposed DS-CTA method is applied to one CT/CTA phantom dataset and seven intra-patient head CT/CTA images. Here, the CT/CTA phantom dataset simulates bone and vessels in the human head and is generated by using a 3D Gaussian blurring model for the spiral CT system [17]. Dimensions and resolutions of all the test images are as shown in Table 1, and some of the images are shown in Figs. 8–11, respectively. To provide the 3D structure intuitively, the figures show volume rendered images for the corresponding datasets. Notice that these volume rendered images are acquired from the original source data without any preprocessing or
segmentation. Usually the sampling interval in the z-direction is sparser than in x- and y-directions in conventional CT images. In this paper, the anisotropic datasets are directly used without preprocessing. The bin size of joint histogram is chosen as 16 HU!16 HU, and the threshold value of MIDV, Thfv, is set to a half of the bin size, or 8 HU. An example of a selected floating volume mask is shown in Fig. 12. Here, the mask including boundary features of the bone and air corresponds to about 31% of the volume of original CT data, and this tendency is consistent for all the datasets we use. The threshold values of bone and air, or Thb and Tha, are chosen as 176 HU and K24 HU, for the intensity value range of CT images from K1024 HU to 2112 HU. The threshold value used in expanding the bone mask, ThG, is set to 120 HU by examining the 1D profile of the CTA image. First, the registration accuracy of the proposed fast algorithm is examined by comparing it with the conventional one. For this examination, we use the CTA image of the phantom and the corresponding CT image that is intentionally moved from the CTA one by (Tx, Ty, Tz, Rx, Ry, Rz) of (3.1, 2.5, K7.8, K0.88, K0.38, 0.28). In Table 2, we can note that the determined registration parameters are the same in both the conventional and proposed fast NMI-based multi-resolution registration methods and coincide with the intentionally applied ones with tiny absolute error. Here, the absolute error is defined as a RMS (root mean squared) value of position errors at the eight corner points of the 3D floating volume (or the CT image), which is normalized by the minimum voxel-interval among those in x-, y-, and z-directions. In the clinical datasets of CT and CTA, however, the registration accuracy cannot be examined quantitatively because the position error between CT and CTA is not known. To effectively demonstrate the registration performance, the two volume datasets of CTA and aligned CT are mixed into one volume dataset by alternating voxel cubes as follows
Table 2 Simulation results using the phantom data with the known mis-alignment parameters
Conventional method
Proposed method
Multi-resolution step
Tx (voxels)
Ty (voxels)
Tz (voxels)
Rx
Ry
Rz
Absolute error (voxels)
(16:16:16)
K0.14
2.90
K7.74
K0.658
1.018
K0.828
9.54
(8:8:8) (4:4:4) (2:2:2) (1:1:1) (16:16:16)
3.08 3.10 3.10 3.10 K0.14
2.49 2.49 2.50 2.50 2.90
K7.80 K7.80 K7.80 K7.80 K7.74
K0.808 K0.808 K0.808 K0.808 K0.658
K0.298 K0.308 K0.308 K0.308 1.018
0.218 0.208 0.208 0.208 K0.828
0.08 0.02 0.01 0.01 9.54
3.08 3.10 3.10 3.10
2.49 2.49 2.50 2.50
K7.80 K7.80 K7.80 K7.80
K0.808 K0.808 K0.808 K0.808
K0.298 K0.308 K0.308 K0.308
0.218 0.208 0.208 0.208
0.08 0.02 0.01 0.01
(8:8:8) (4:4:4) (2:2:2) (1:1:1)
Here, both the conventional and proposed methods are based on NMI and a multi-resolution approach, while the latter uses a floating volume mask Mfv($).
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
397
Fig. 13. Registration results using (a) the conventional registration method and (b) proposed fast registration method for Patient_1 and Patient_4, respectively. To examine the registration accuracy, cubes of CTA and aligned CT images are alternatively displayed into an image.
VMixed ðx;y;zÞ 8
x y z < C C is even f ðx;y;zÞ if Z Wcube Wcube Wcube : gðx;y;zÞ otherwise (12) where VMixed(x,y,z) and Wcube represent the mixed volume data and the width of cube, respectively. In Fig. 13, registration results are given for a couple of clinical datasets. Even in the alternating cube display, the volume rendered image of mixed volume provides a smooth skull base, which proves the accurate registration. Meanwhile, scattering in the vessel regions is caused by the alternation of two volumes with and without high intensity vessels, and the rough volume boundary is due to the lack of data in one of the CT and CTA images; hence, those are not relevant to the accuracy of registration. Based on the quantitative analysis for the phantom data in Table 2 and the qualitative analysis for clinical data in Fig. 13, we may conclude that the proposed fast registration algorithm does not lose any registration accuracy compared to the conventional method. Since the registration time is the major part of the total computing time for DS-CTA, we examine it for all the datasets by using a PC with an Intel Pentium IV CPU of 2.4 GHz (see Table 3). As shown in the table, the speed-up ratio of the proposed fast algorithm is about 1.74w3.01 compared with the conventional registration method.
The registration computing time of 2–8 min, which is reduced from 4–18 min, seems to be acceptable in many clinical applications, since the additional computing time for subtraction and masking operation requires about a minute. The results of the proposed DS-CTA method are shown in Figs. 14–17. Figs. 14 and 15 show MIP (maximal intensity projection) images for clinical datasets Patient_1 and Patient_3. It is clear in the figures that simple subtraction is not enough to produce a satisfactory extraction result even after the fine 3D registration. Meanwhile, the proposed novel 3D refinement scheme based on adaptive masking provides results good enough for Table 3 Registration time for 8 datasets CT/CTA datasets
Conventional method (s)
Proposed method (s)
Speed-up ratio
Volume of mask (%)
Phantom Patient_1 Patient_2 Patient_3 Patient_4 Patient_5 Patient_6 Patient_7
254 850 1049 825 218 652 956 909
118 389 349 475 112 365 454 364
2.15 2.19 3.01 1.74 1.95 1.79 2.10 2.50
17 31 24 52 47 40 38 32
Intel P4 of 2.4 GHz is used for experiment.
398
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
Fig. 14. Vessel extraction results from dataset Patient_1. All the images are represented by MIP. (a) CTA image, (b) brute-force subtraction image after registration and (c) refined image by the proposed adaptive masking scheme.
Fig. 15. Vessel extraction results from dataset Patient_3. All the images are represented by MIP. (a) CTA image, (b) brute-force subtraction image after registration and (c) refined image by the proposed adaptive masking scheme.
diagnosis. This can be easily observed through the 1D profiles in Fig. 16. As shown in the graphs, simple subtraction provides considerable residuals near the bone (voxel positions of around 245, 285, and 310), while the proposed mask successfully excludes unwanted residuals and identifies vessels only.
Fig. 17 shows the volume rendered images of CTA and aligned CT data of Patient_4 and the MIP and volume rendered images of the DS-CTA result. The MIP and volume rendered images by the proposed vessel extraction scheme in Fig. 17(c) and (d) clearly show all the vessels including the hidden ones, without artifacts due to bone fragments.
Fig. 16. 1D profiles along the white line in Fig. 14.
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
Fig. 17. Vessel extraction results from dataset Patient_4. (a) CT image and (b) CTA image represented by volume rendering; and extracted vessel images represented by (c) MIP and (d) volume rendering, respectively.
399
registration algorithm selects a sub-volume out of floating volume data (or the CT data) using a newly defined measure, MIDV, and calculates NMI only inside the subvolume. Thereby, the computational complexity is reduced to 33–58% while maintaining the registration accuracy. After the registration, the blood vessel image, obtained by subtracting the aligned CT from CTA, still contained unwanted residuals. Hence, to remove those residuals, we developed a novel 3D refinement algorithm by using a mask, which represents definite bone and air areas. The proposed 3D DS-CTA method was applied to a phantom CT/CTA head dataset and seven clinical CT/CTA head datasets and was found to be fast, accurate, and robust. The proposed method well extracts intracranial vessels from the CTA image with almost no loss. The method also exposes the blood vessels hidden by bone in CTA image. The total execution time is short enough for many clinical applications.
5. Conclusions
References
In this paper, a fast and accurate method is proposed for the 3D DS-CTA using CT and CTA images. The proposed DS-CTA method is composed of the fast registration algorithm and the refinement algorithm after subtraction. The fast 3D registration algorithm reduces the computational complexity to 33–58% while maintaining the registration accuracy, and the novel 3D refinement algorithm after subtraction is found to be robust and efficiently removes unwanted bone residuals. Therefore, the proposed DS-CTA method well extracts intracranial blood vessels from the CTA image with almost no loss, and exposes the blood vessels that are hidden by bone in the CTA image. The total execution time in a PC is short enough for many clinical applications.
[1] Yoshimoto T, Mizoi K. Importance of management of unruptured cerebral aneurysms. Surg Neurol 1997;47:522–6. [2] Abrahams JM, Saha PK, Hurst RW, LeRoux PD, Udupa JK. Threedimensional bone-free rendering of the cerebral circulation by use of computed tomographic angiography and fuzzy connectedness. Neurosurgery 2002;51(1):264–9. [3] Yan C, Hirano S, Hata Y. Extraction of blood vessel in CT angiography image aided by fuzzy logic. Proc IEEE ICSP 2000;2: 926–9. [4] Hunter GJ, Hamberg LM, Ponzo JA, Huanghellinger FR, Morris PP, Rabinov J, Farkas J, Lev MH, Schaefer PW, Ogilvy CS, Schwamm L, Buonanno FS, Koroshetz WJ, Wolf GL, Gonzalez RG. Assessment of cerebral perfusion and arterial anatomy in hyperacute stroke with three-dimensional functional CT: early clinical results. AJNR 1998; 19:29–37. [5] Yeung MM, Yeo BL, Liou SP, Banihashemi A. Three-dimensional image registration for spiral CT angiography. Proc IEEE CVPR 1994;423–9. [6] Venema HW, Hulsmans FJH, Heeten GJD. CT angiography of the circle of willis and intracranial internal carotid arteries: maximum intensity projection with matched mask bone elimination-feasibility study. Radiology 2001;218(3):893–8. [7] Maes F, Collignon A, Vandermeulen D, Marchal G, Suetens P. Multimodality image registration by maximization of mutual information. IEEE T Med Imaging 1997;16(2):187–98. [8] Well III WM, Viola P, Atsumi H, Nakajima S, Kikins R. Multimodal volume registration by maximization of mutual information. Med Image Anal 1996;1(1):35–51. [9] Studholme C, Hill DLG, Hawkes DJ. Automated 3D registration of MR and CT images in the head. Med Image Anal 1996;1(2):163–75. [10] Studholme C, Hill DLG, Hawkes DJ. An overlap invariant entropy measure of 3D medical image alignment. Pattern Recogn 1999;32(1): 71–86. [11] Maes F, Vandermeulen D, Suetens P. Comparative evaluation of multiresolution optimization strategies for multimodality image registration by maximization of mutual information. Med Image Anal 1999;3(4):373–86. [12] Thevenaz P, Unser M. Optimization of mutual information for multiresolution image registration. IEEE T Image Process 2000;9(12): 2083–99.
6. Summary Recently, angiographic investigation using high-resolution multi-detector CT images has become a useful tool for the sensitive assessment of vascular pathological changes. Despite the high resolution of CTA, this representation provides the blood vessel information merged with bone information and does not allow a straightforward understanding of vascular structure. In the human head, for example, the accurate diagnosis of cerebral arterial aneurysm is important, but since some blood vessel parts of interest can be located near the high intensity skull base region, a complicated procedure is needed to separate blood vessels from bone. We have developed a new, fast, and accurate method for the 3D DS-CTA using CT and CTA images acquired from the same patient. The effective multi-resolution 3D
400
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
[13] Capek M, Mroz L, Wegenkittl R. Robust and fast medical registration of 3D-multi-modality data sets Technical Report. Vol. TR-VRVIS-2001-008: Research center for virtual reality and visualization 2001. [14] Pluim JPW, Maintz JBA, Viergever MA. Interpolation artefacts in mutual information-based image registration. Comput Vis Image Und 2000;77(2):211–32. [15] Shuqian L, Xiang L. Implementation of mutual information based multi-modality medical image registration. Proc IEEE Eng Med Biol 2000;2:1447–50. [16] West J, Fitzpatrick JM, Wang MY, Dawant BM, Maurer CR, Kessler RM, Maciunas RJ, Barillot C, Lemoine D, Collignon A, Maes F, Suetens P, Vandermeulen D, van den Elsen PA, Napel S, Sumanaweera TS, Harkness B, Hemler PF, Hill DLG, Hawkes DJ, Studholme C, Maintz JBA, Viergever MA, Malandain G, Pennec X, Noz ME, Maguire GQ, Pollack M, Pelizzari CA, Robb RA, Hanson D, Woods RP. Comparison and evaluation of retrospective intermodality brain image registration techniques. J Comp Assist Tomogr 1997; 21(4):554–66. [17] Wang G, Vaanier MW, Skinner MW, Cavalcanti MGP, Harding GW. Spiral CT image deblurring for cochlear implantation. IEEE T Med Imaging 1998;17(2):251–61.
Sung Min Kwon received the BS and MS degrees in Electrical Engineering from KAIST, Korea in 1997 and 1999, respectively. He is working toward PhD degree in the Department of Electrical Engineering and Computer Science (EECS), KAIST, Korea. His research interests include medical image processing, 3D visualization, and volume registration.
Yong Sun Kim received the BS and MS degrees in Electrical Engineering from KAIST, Korea in 2002 and 2004, respectively. He is a PhD candidate of the Department of EECS, KAIST. His research interests include medical image processing and registration.
Tae-Sung Kim received his medical degree from Yonsei University, Seoul, Korea, in 2000. He completed his internship and residency in Diagnostic Radiology at Yonsei University Medical Center, in 2001 and 2005, respectively. His research interests are medical image processing and medical informatics.
Jong Beom Ra received the BS degree in Electronic Engineering from Seoul National University, Korea, in 1975, and the MS and PhD degrees in Electrical Engineering from KAIST, Korea, in 1977 and 1983, respectively. From 1983 to 1987, he served as an associate research scientist in Radiology, Columbia University, New York, USA. He joined the Department of EECS in KAIST, in July 1987, where he is now a professor and the director of medical imaging research center, KAIST. His research interests are digital image and video processing, medical imaging, and 3-D visualization. He is a senior member of IEEE.
Digital subtraction CT angiography based on efficient 3D registration and refinement Sung Min Kwona, Yong Sun Kima, Tae-Sung Kimb, Jong Beom Raa,* a
Department of Electrical Engineering and Computer Science, Korea Advanced Institute of Science and Technology, 373-1, Guseongdong, Yuseonggu, Daejeon 305-701, South Korea b Department of Diagnostic Radiology, College of Medicine, Yonsei University, Seoul, South Korea Received 30 January 2004; accepted 1 June 2004
Abstract A novel method for fast, automatic 3D digital subtraction CT angiography (DS-CTA) is presented to generate artifact-free angiograms. The proposed method consists of two steps: 3D registration to align a CT image to the CT angiography (CTA) image and subtraction-andrefinement to extract blood vessels only. For efficient and accurate 3D registration in the first step, an normalized mutual information (NMI) based algorithm is adopted, and its fast version is developed by introducing a new measure. To further improve the subtracted image quality in the second step, a novel 3D refinement algorithm is suggested to effectively remove unwanted residuals. Experimental results of seven clinical CT/CTA head datasets demonstrate that cerebral vessels are well extracted from CTA images with almost no loss. The typical processing time is 3–9 min depending on the image size in a PC with a 2.4 GHz CPU. q 2004 Elsevier Ltd. All rights reserved. Keywords: Digital subtraction; Angiography; CTA; DS-CTA; Registration; Mutual information; Refinement; Mask; Cerebral vessels
1. Introduction Recently, angiographic investigation using highresolution multi-detector (MD) CT images has become a useful tool for the sensitive assessment of vascular structures. In the head, for example, most of blood supplied to brain passes the special region called the Circle of Willis. Hence, its branching portions are continuously stressed and may undergo aneurismal change. Consequently, they are apt to be ruptured and a devastating subarachnoid hemorrhage may occur [1]; therefore, the accurate diagnosis of cerebral arterial aneurysm is critical. Since a CTA image describes blood vessels with high intensity, it is known to be very helpful for diagnosing this kind of vessel disease; however, since blood vessel parts of interest are often located near the high intensity skull base, general CTA with bones may not be quite helpful for accurate diagnosis, so a sophisticated
* Corresponding author. Tel.: C82-42-869-3434; fax: C82-42-8698360. E-mail address: [email protected] (J.B. Ra). 0895-6111/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2004.06.003
procedure to separate only blood vessels from bone is needed. Several methods to extract blood vessels from CTA images have been tried, based on segmentation or tracking [2,3]; however, those methods tend to leave some artifacts of bone near vessels and lose thin vessels. Meanwhile, a couple of vessel extraction techniques have been attempted on the basis of the CT and CTA image registration, namely, the subtraction of the CT image from a CTA image [4,5] and bone-elimination by replacing the intensities of bone with a soft tissue intensity [6]. In Ref. [4], during CT/CTA acquisition, a patient’s head is immobilized in a headholder by using forehead and chin straps. In this case, the patient feels uncomfortable during the data acquisition and slight motion due to imperfect fixation is a strong limitation to this subtraction technique. A feature-based registration for correcting patient motion between CT/CTA images releases the patient from the discomfort of the head-holder [5]; however, a subtle mis-registration of the two images may cause some bone fractions near the bone boundary in the subtracted image. Therefore, for a satisfactory vessel extraction, an accurate, robust 3D registration is
392
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
indispensable. Maximization of mutual information (MI) has been recommended as a powerful intrinsic criterion for accurate multi-modal or mono-modal image registration [7–10]. However, MI requires long processing time, especially for a recent large size 3D volume dataset obtained from a high resolution CT. Practically, there are many urgent cases in the diagnosis using DS-CTA, such as rupturing cerebral aneurysms; therefore, the importance of fast 3D registration cannot be overestimated, and its processing time should be improved. Even with the accurate registration between CT and CTA images by using an advanced scheme, their subtraction and bone-elimination results still contain unwanted residuals at boundaries between bone and tissue, which is due to and the limited resolution of CT and CTA images, and these residuals cause annoying artifacts in the angiogram. In this paper, we develop a novel DS-CTA scheme by adopting an accurate NMI-based registration technique. By solving the two problems in the scheme, namely, the calculation complexity of the NMI-based registration and the residual artifact due to subtle mis-registration, we propose a fast automatic 3D DS-CTA method based on CT and CTA images. The scheme is easily applicable to a PC environment. In Section 2, the proposed fast registration method is described, and in Section 3, the 3D refinement process to remove unwanted artifacts is described in detail. Experimental results are presented in Section 4, and the conclusions are drawn in Section 5.
2. The proposed effective 3D registration based on NMI The CT and CTA head images acquired from a patient include the same information of bone and soft-tissue with a certain mis-alignment, while the blood vessel information is included only in the CTA image. Therefore, the blood vessel information can be extracted using CT and CTA images simultaneously. Fig. 1 shows the overall block diagram of
Fig. 1. Overall block diagram of the adopted 3D DS-CTA algorithm.
Fig. 2. Block diagram of an NMI-based registration algorithm.
the 3D DS-CTA algorithm we adopt. As shown in the figure, the algorithm first registers a CT image to the CTA image, subtracts it from the CTA image, and finally refines the subtracted image by removing unwanted residuals or artifacts. As mentioned in Section 1, the precise registration is very important for a satisfactory subtraction result. Hence, we adopt the NMI-based scheme for the registration of CT and CTA images due to its superior performance, such as robust and overlap-invariant characteristics [10]. In the following sub-sections, we first briefly describe the NMIbased scheme and a multi-resolution approach, which we adopt to our algorithm for accurate and fast registration. Then, we propose a volume mask based registration scheme to further accelerate the processing speed. 2.1. Conventional NMI-based multi-resolution registration Fig. 2 shows the block diagram for a conventional NMIbased registration. Registration is a procedure to find the optimal transformation matrix T projecting a floating volume into the reference frame as shown in Fig. 3, i.e. a procedure to match the two volumes best using a similarity
Fig. 3. Registration procedure of a floating volume to the reference frame by transformation matrix T.
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
measure. NMI is an overlap-invariant entropy measure that correctly describes a dispersion of the joint histogram for two random variables and measures their statistical dependence. Hence, NMI has been extensively investigated as a similarity measure [7–16] and has been proven to perform well in 3D medical image alignment [10]. NMI is defined as YðA; BÞ Z
HðAÞ C HðBÞ HðA; BÞ
complexity of multi-resolution registration is represented as CC Z
LK1 X ð the number of trialsÞl !ðcomputational lZ0
complexity of a similarity measureÞl f
LK1 X
ð of trialsÞl
lZ0
!ðthe number of sub-sampled floating (1)
where H(A) and H(B) are the marginal entropy of random variable A and B, respectively, and H(A, B) is the joint entropy of A and B. These three entropies are calculated from the joint probability distribution of A and B using the Shannon measure of entropy, and the joint probability distribution can be estimated from the joint histogram as shown in Fig. 4; therefore, the NMI is obtained if the joint histogram is available. During the registration procedure in Fig. 2, the elements of a transformation matrix T are updated to maximize the NMI in each trial. Accordingly, the corresponding projected position of the floating volume is also refined. In order to generate a joint histogram based on intensity pairs in each trial (see Fig. 4), we need intensity values at the reference volume positions corresponding to the projected positions of all the voxels in the floating volume. In obtaining those intensity values, tri-linear interpolation is to be performed in the reference volume. This extensive tri-linear interpolation processing for all floating voxels in each trial requires massive computation, which is the major reason for the long processing time in NMI-based registration. To overcome this problem, a multi-resolution approach is adopted by using a hierarchical structure of input data [11–13]. For multi-resolution image registration, a floating volume is sub-sampled with an interval of 2LKlK1 voxels at the lth level (lZ0,1,.,LK1). Then, based on the similarity measure that is obtained by using the subsampled floating volume at each level, the transformation matrix T is searched and refined starting from the initial level. As the level increases, the sub-sampling interval decreases and the registration result becomes more accurate with higher computational complexity. The computational
Fig. 4. Joint histogram of intensity pairs.
393
volume voxelsÞl (2) Even though this multi-resolution approach can considerably reduce the computational complexity, we still need to reduce it. In Eq. (2), the majority of the computation results from a calculation of a similarity measure at the final level because the number of sub-sampled floating volume voxels increases exponentially as the level becomes higher. Hence we will focus the reduction of the number of voxels to be used in the calculation of the similarity measure at the higher levels. 2.2. The proposed fast registration A fast 3D registration algorithm based on the multiresolution scheme is given in Fig. 5. Here, we introduce a floating volume mask to the conventional NMI-based multiresolution registration scheme described above so that a reduced number of voxels are used to evaluate NMI in the repeated trials to find the best transformation matrix at each level. Fig. 6 demonstrates a mask setting procedure at level l.
Fig. 5. Block diagram of the proposed 3D multi-resolution registration algorithm.
394
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
the registration procedure. Note here that Thfv is properly decided according to the bin size of the joint histogram.
3. 3D refinement of the subtracted image
Fig. 6. The procedure to obtain a floating volume mask at resolution level l.
As shown in Fig. 1, the aligned CT image is subtracted from the CTA image to extract blood vessels. Even though the alignment result by registration is fairly good, the subtracted image still contains unwanted residuals from a subtle mis-alignment at bone and air boundaries due to abrupt and considerable intensity changes; moreover, the residuals hinder proper examining blood vessels because their intensity values are still comparable with blood vessel intensity values. To get rid of these unwanted residuals, we introduce and apply a masking procedure after
This procedure begins with checking the convergence rate (CR) of the previous level, CRlK1, which is defined as follows vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 8 u1 X CRlK1 Z t ! jðtransformed position of vi with T1;lK1 Þ K ðtransformed position of vi with TK;lK1 Þj2 8 iZ1 where vi represents eight corner voxels of a floating volume, and T1,lK1 and TK,lK1 denote the initially and finally obtained transformation matrices at level lK1, respectively. If CRlK1, is smaller than the minimum voxel-interval dmin among the ones along the x-, y-, and z-axes, we set a floating volume mask Mfv($) that consists of the voxels having an MIDV (maximum intensity difference of a voxel) larger than a threshold. Otherwise, the whole floating volume is set to the mask. Therefore, at higher levels where the number of subsampled voxels increases and the convergence rate becomes slow, we can partially select the floating volume voxels, reducing the calculation time for a similarity measure, NMI. Now, we consider a measure MIDV for the proper selection of a mask. As described in Section 2.1, during the registration procedure, the elements of T are refined to the direction increasing the NMI value. If the intensity change of an intensity pair for a voxel between two consecutive trials is too little to affect the joint histogram, the corresponding voxel can be considered as not contributing to any change to the value of NMI; therefore, those voxels may be removed from the registration procedure without degrading the registration accuracy. To discriminate those voxels, we introduce a measure MIDV for each floating volume voxel, which is defined as dmin MIDV Z max jIðvi Þ K Iðvc Þj ! ; i Z 1; .; 26 : (4) d vi Here, vi represents 26 neighboring voxels of a center voxel vc, dvi denotes the distance between vi and vc, and I($) denotes the intensity of a voxel. Note that MIDV represents the degree of intensity change from a voxel vc. If the MIDV of a floating volume voxel is smaller than a threshold value Thfv, the voxel is excluded from a mask, thereby from
(3)
the subtraction, i.e. Vðx; y; zÞ Z Mba ðx; y; zÞ !Sðx; y; zÞ; where
(
Sðx; y; zÞ Z
(5)
f ðx; y; zÞ K gðx; y; zÞ
if f ðx; y; zÞO gðx; y; zÞ;
0
otherwise:
(6) Here, V(x,y,z), Mba(x,y,z), S(x,y,z), f(x,y,z), and g(x,y,z) represent the final vessel image, mask, subtraction image, CTA image, and aligned CT image, respectively. Since unwanted residuals as well as blood vessels are contained in S(x,y,z), S(x,y,z) is the target image for masking operation. The purpose of masking is to effectively remove unwanted residuals by using a mask function that represents definite bone and air areas excluding vessels. The procedure to obtain the mask, Mba(x,y,z) is shown in Fig. 7. We first consider the bone mask Mb that represents definite bone areas, and then an initial bone mask is generated in the aligned CT image by properly selecting a threshold value of bone, Thb, based on the intensity histogram of CT image, i.e. ( 1 if gðx; y; zÞO Thb ; Mb ðx; y; zÞ Z (7) 0 otherwise:
Fig. 7. The procedure to obtain the final mask Mba.
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
395
Table 1 Dimension and resolution of test image sets CT/CTA datasets
Dimension
Resolution (mm3)
Phantom Patient_1 Patient_2 Patient_3 Patient_4 Patient_5 Patient_6 Patient_7
512!512!180 512!512!185 512!512!185 512!512!163 512!512!215 512!512!147 512!512!201 512!512!180
0.300!0.300!0.300 0.268!0.268!0.300 0.268!0.268!0.300 0.326!0.326!0.500 0.301!0.301!0.300 0.289!0.289!0.500 0.275!0.275!0.300 0.254!0.254!0.300
All the datasets include CT and CTA images having the same dimension and resolution.
Note that the initial bone mask is made not in the CTA image but in the aligned CT image where the vessel intensity is much lower than that of bone; therefore, reliable bone areas can be selected as the initial bone mask. Since we assume the inherent subtle mis-alignment between registered CT and CTA images, however, the initial bone mask defined on the registered CT image is not suitable for completely removing unwanted bone residuals from the subtraction image. Hence, we refine the initial bone mask using the CTA image. Since the initial bone mask obtained through simple thresholding does not cover all the bone areas, the mask should be expanded until it covers all the bone areas so that the final bone mask may remove all the unwanted bone fractions from the subtraction image. It is obvious in the CTA image that the gradient is steep at the boundary between bone and tissue; hence, in the proposed algorithm, the bone mask is expanded voxel by voxel toward the outside of its boundary if the intensity value of the corresponding CTA voxel differs enough from any one of its 26 neighboring voxel intensity values. In other words, a background voxel adjacent to the mask boundary is included as part of the bone mask, if it is assured to belong to a transition region from bone to other tissue. The bone mask updating process can be represented as follows Mb ðx C i; y C j; z C kÞ Z 1; if Mb ðx; y; zÞ Z 1; & Mb ðx C i; y C j; z C kÞ Z 0;
(8)
& f ðx C i; y C j; z C kÞ! f ðx; y; zÞ K ThG ; where ThG denotes a threshold value for the intensity difference in the CTA. Here, i, j, kZK1, 0, 1 except for
Fig. 8. Volume rendered images of a phantom dataset. (a) CT and (b) CTA.
Fig. 9. Volume rendered images of a dataset, Patient_1. (a) CT and (b) CTA.
the case of iZjZkZ0. This procedure is repeated until the expansion stops naturally at the end of the bone boundary where the gradient becomes slow. It is interesting to note that in the case where a vessel is located near bone, the mask expansion based on Eq. (8) can prohibit the mask from intruding into the vessel because the gradient becomes slow at the boundary between bone and vessel. An initial air mask is also generated to get rid of unwanted residuals between soft tissue (or bone) and air. Since there is no vessel in direct contact with air and the air provides a fairly constant intensity value, the initial air mask is generated on the basis of the intensity value rather than the gradient value as in the bone mask. The initial air mask is defined as Ma ðx; y; zÞ Z 1 if gðx; y; zÞ! Tha ;
(9)
and the final air mask Ma is simply refined on the basis of CTA image intensities as follows Ma ðx C i; y C j; z C kÞ Z 1; if Ma ðx; y; zÞ Z 1 & Ma ðx C i; y C j; z C kÞ Z 0
(10)
& f ðx C i; y C j; z C kÞ! Tha ; where Tha is the threshold value based on the intensity histogram of CT image, and i, j, kZK1, 0, 1, except iZjZ kZ0. Then, the final mask is obtained by simply combining the bone and air masks, i.e. Mba ðx; y; zÞ Z ðMb ðx; y; zÞ C Ma ðx; y; zÞÞ:
(11)
This mask is used for excluding definite bone and air areas.
Fig. 10. Volume rendered images of a dataset, Patient_3. (a) CT and (b) CTA.
396
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
Fig. 11. Volume rendered images of a dataset, Patient_4. (a) CT and (b) CTA.
Fig. 12. Floating volume mask. (a) A slice image of CT data for floating volume of Patient_1 and (b) its selected mask with Thfv of 8 HU.
4. Experimental results and discussions The proposed DS-CTA method is applied to one CT/CTA phantom dataset and seven intra-patient head CT/CTA images. Here, the CT/CTA phantom dataset simulates bone and vessels in the human head and is generated by using a 3D Gaussian blurring model for the spiral CT system [17]. Dimensions and resolutions of all the test images are as shown in Table 1, and some of the images are shown in Figs. 8–11, respectively. To provide the 3D structure intuitively, the figures show volume rendered images for the corresponding datasets. Notice that these volume rendered images are acquired from the original source data without any preprocessing or
segmentation. Usually the sampling interval in the z-direction is sparser than in x- and y-directions in conventional CT images. In this paper, the anisotropic datasets are directly used without preprocessing. The bin size of joint histogram is chosen as 16 HU!16 HU, and the threshold value of MIDV, Thfv, is set to a half of the bin size, or 8 HU. An example of a selected floating volume mask is shown in Fig. 12. Here, the mask including boundary features of the bone and air corresponds to about 31% of the volume of original CT data, and this tendency is consistent for all the datasets we use. The threshold values of bone and air, or Thb and Tha, are chosen as 176 HU and K24 HU, for the intensity value range of CT images from K1024 HU to 2112 HU. The threshold value used in expanding the bone mask, ThG, is set to 120 HU by examining the 1D profile of the CTA image. First, the registration accuracy of the proposed fast algorithm is examined by comparing it with the conventional one. For this examination, we use the CTA image of the phantom and the corresponding CT image that is intentionally moved from the CTA one by (Tx, Ty, Tz, Rx, Ry, Rz) of (3.1, 2.5, K7.8, K0.88, K0.38, 0.28). In Table 2, we can note that the determined registration parameters are the same in both the conventional and proposed fast NMI-based multi-resolution registration methods and coincide with the intentionally applied ones with tiny absolute error. Here, the absolute error is defined as a RMS (root mean squared) value of position errors at the eight corner points of the 3D floating volume (or the CT image), which is normalized by the minimum voxel-interval among those in x-, y-, and z-directions. In the clinical datasets of CT and CTA, however, the registration accuracy cannot be examined quantitatively because the position error between CT and CTA is not known. To effectively demonstrate the registration performance, the two volume datasets of CTA and aligned CT are mixed into one volume dataset by alternating voxel cubes as follows
Table 2 Simulation results using the phantom data with the known mis-alignment parameters
Conventional method
Proposed method
Multi-resolution step
Tx (voxels)
Ty (voxels)
Tz (voxels)
Rx
Ry
Rz
Absolute error (voxels)
(16:16:16)
K0.14
2.90
K7.74
K0.658
1.018
K0.828
9.54
(8:8:8) (4:4:4) (2:2:2) (1:1:1) (16:16:16)
3.08 3.10 3.10 3.10 K0.14
2.49 2.49 2.50 2.50 2.90
K7.80 K7.80 K7.80 K7.80 K7.74
K0.808 K0.808 K0.808 K0.808 K0.658
K0.298 K0.308 K0.308 K0.308 1.018
0.218 0.208 0.208 0.208 K0.828
0.08 0.02 0.01 0.01 9.54
3.08 3.10 3.10 3.10
2.49 2.49 2.50 2.50
K7.80 K7.80 K7.80 K7.80
K0.808 K0.808 K0.808 K0.808
K0.298 K0.308 K0.308 K0.308
0.218 0.208 0.208 0.208
0.08 0.02 0.01 0.01
(8:8:8) (4:4:4) (2:2:2) (1:1:1)
Here, both the conventional and proposed methods are based on NMI and a multi-resolution approach, while the latter uses a floating volume mask Mfv($).
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
397
Fig. 13. Registration results using (a) the conventional registration method and (b) proposed fast registration method for Patient_1 and Patient_4, respectively. To examine the registration accuracy, cubes of CTA and aligned CT images are alternatively displayed into an image.
VMixed ðx;y;zÞ 8
x y z < C C is even f ðx;y;zÞ if Z Wcube Wcube Wcube : gðx;y;zÞ otherwise (12) where VMixed(x,y,z) and Wcube represent the mixed volume data and the width of cube, respectively. In Fig. 13, registration results are given for a couple of clinical datasets. Even in the alternating cube display, the volume rendered image of mixed volume provides a smooth skull base, which proves the accurate registration. Meanwhile, scattering in the vessel regions is caused by the alternation of two volumes with and without high intensity vessels, and the rough volume boundary is due to the lack of data in one of the CT and CTA images; hence, those are not relevant to the accuracy of registration. Based on the quantitative analysis for the phantom data in Table 2 and the qualitative analysis for clinical data in Fig. 13, we may conclude that the proposed fast registration algorithm does not lose any registration accuracy compared to the conventional method. Since the registration time is the major part of the total computing time for DS-CTA, we examine it for all the datasets by using a PC with an Intel Pentium IV CPU of 2.4 GHz (see Table 3). As shown in the table, the speed-up ratio of the proposed fast algorithm is about 1.74w3.01 compared with the conventional registration method.
The registration computing time of 2–8 min, which is reduced from 4–18 min, seems to be acceptable in many clinical applications, since the additional computing time for subtraction and masking operation requires about a minute. The results of the proposed DS-CTA method are shown in Figs. 14–17. Figs. 14 and 15 show MIP (maximal intensity projection) images for clinical datasets Patient_1 and Patient_3. It is clear in the figures that simple subtraction is not enough to produce a satisfactory extraction result even after the fine 3D registration. Meanwhile, the proposed novel 3D refinement scheme based on adaptive masking provides results good enough for Table 3 Registration time for 8 datasets CT/CTA datasets
Conventional method (s)
Proposed method (s)
Speed-up ratio
Volume of mask (%)
Phantom Patient_1 Patient_2 Patient_3 Patient_4 Patient_5 Patient_6 Patient_7
254 850 1049 825 218 652 956 909
118 389 349 475 112 365 454 364
2.15 2.19 3.01 1.74 1.95 1.79 2.10 2.50
17 31 24 52 47 40 38 32
Intel P4 of 2.4 GHz is used for experiment.
398
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
Fig. 14. Vessel extraction results from dataset Patient_1. All the images are represented by MIP. (a) CTA image, (b) brute-force subtraction image after registration and (c) refined image by the proposed adaptive masking scheme.
Fig. 15. Vessel extraction results from dataset Patient_3. All the images are represented by MIP. (a) CTA image, (b) brute-force subtraction image after registration and (c) refined image by the proposed adaptive masking scheme.
diagnosis. This can be easily observed through the 1D profiles in Fig. 16. As shown in the graphs, simple subtraction provides considerable residuals near the bone (voxel positions of around 245, 285, and 310), while the proposed mask successfully excludes unwanted residuals and identifies vessels only.
Fig. 17 shows the volume rendered images of CTA and aligned CT data of Patient_4 and the MIP and volume rendered images of the DS-CTA result. The MIP and volume rendered images by the proposed vessel extraction scheme in Fig. 17(c) and (d) clearly show all the vessels including the hidden ones, without artifacts due to bone fragments.
Fig. 16. 1D profiles along the white line in Fig. 14.
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
Fig. 17. Vessel extraction results from dataset Patient_4. (a) CT image and (b) CTA image represented by volume rendering; and extracted vessel images represented by (c) MIP and (d) volume rendering, respectively.
399
registration algorithm selects a sub-volume out of floating volume data (or the CT data) using a newly defined measure, MIDV, and calculates NMI only inside the subvolume. Thereby, the computational complexity is reduced to 33–58% while maintaining the registration accuracy. After the registration, the blood vessel image, obtained by subtracting the aligned CT from CTA, still contained unwanted residuals. Hence, to remove those residuals, we developed a novel 3D refinement algorithm by using a mask, which represents definite bone and air areas. The proposed 3D DS-CTA method was applied to a phantom CT/CTA head dataset and seven clinical CT/CTA head datasets and was found to be fast, accurate, and robust. The proposed method well extracts intracranial vessels from the CTA image with almost no loss. The method also exposes the blood vessels hidden by bone in CTA image. The total execution time is short enough for many clinical applications.
5. Conclusions
References
In this paper, a fast and accurate method is proposed for the 3D DS-CTA using CT and CTA images. The proposed DS-CTA method is composed of the fast registration algorithm and the refinement algorithm after subtraction. The fast 3D registration algorithm reduces the computational complexity to 33–58% while maintaining the registration accuracy, and the novel 3D refinement algorithm after subtraction is found to be robust and efficiently removes unwanted bone residuals. Therefore, the proposed DS-CTA method well extracts intracranial blood vessels from the CTA image with almost no loss, and exposes the blood vessels that are hidden by bone in the CTA image. The total execution time in a PC is short enough for many clinical applications.
[1] Yoshimoto T, Mizoi K. Importance of management of unruptured cerebral aneurysms. Surg Neurol 1997;47:522–6. [2] Abrahams JM, Saha PK, Hurst RW, LeRoux PD, Udupa JK. Threedimensional bone-free rendering of the cerebral circulation by use of computed tomographic angiography and fuzzy connectedness. Neurosurgery 2002;51(1):264–9. [3] Yan C, Hirano S, Hata Y. Extraction of blood vessel in CT angiography image aided by fuzzy logic. Proc IEEE ICSP 2000;2: 926–9. [4] Hunter GJ, Hamberg LM, Ponzo JA, Huanghellinger FR, Morris PP, Rabinov J, Farkas J, Lev MH, Schaefer PW, Ogilvy CS, Schwamm L, Buonanno FS, Koroshetz WJ, Wolf GL, Gonzalez RG. Assessment of cerebral perfusion and arterial anatomy in hyperacute stroke with three-dimensional functional CT: early clinical results. AJNR 1998; 19:29–37. [5] Yeung MM, Yeo BL, Liou SP, Banihashemi A. Three-dimensional image registration for spiral CT angiography. Proc IEEE CVPR 1994;423–9. [6] Venema HW, Hulsmans FJH, Heeten GJD. CT angiography of the circle of willis and intracranial internal carotid arteries: maximum intensity projection with matched mask bone elimination-feasibility study. Radiology 2001;218(3):893–8. [7] Maes F, Collignon A, Vandermeulen D, Marchal G, Suetens P. Multimodality image registration by maximization of mutual information. IEEE T Med Imaging 1997;16(2):187–98. [8] Well III WM, Viola P, Atsumi H, Nakajima S, Kikins R. Multimodal volume registration by maximization of mutual information. Med Image Anal 1996;1(1):35–51. [9] Studholme C, Hill DLG, Hawkes DJ. Automated 3D registration of MR and CT images in the head. Med Image Anal 1996;1(2):163–75. [10] Studholme C, Hill DLG, Hawkes DJ. An overlap invariant entropy measure of 3D medical image alignment. Pattern Recogn 1999;32(1): 71–86. [11] Maes F, Vandermeulen D, Suetens P. Comparative evaluation of multiresolution optimization strategies for multimodality image registration by maximization of mutual information. Med Image Anal 1999;3(4):373–86. [12] Thevenaz P, Unser M. Optimization of mutual information for multiresolution image registration. IEEE T Image Process 2000;9(12): 2083–99.
6. Summary Recently, angiographic investigation using high-resolution multi-detector CT images has become a useful tool for the sensitive assessment of vascular pathological changes. Despite the high resolution of CTA, this representation provides the blood vessel information merged with bone information and does not allow a straightforward understanding of vascular structure. In the human head, for example, the accurate diagnosis of cerebral arterial aneurysm is important, but since some blood vessel parts of interest can be located near the high intensity skull base region, a complicated procedure is needed to separate blood vessels from bone. We have developed a new, fast, and accurate method for the 3D DS-CTA using CT and CTA images acquired from the same patient. The effective multi-resolution 3D
400
S.M. Kwon et al. / Computerized Medical Imaging and Graphics 28 (2004) 391–400
[13] Capek M, Mroz L, Wegenkittl R. Robust and fast medical registration of 3D-multi-modality data sets Technical Report. Vol. TR-VRVIS-2001-008: Research center for virtual reality and visualization 2001. [14] Pluim JPW, Maintz JBA, Viergever MA. Interpolation artefacts in mutual information-based image registration. Comput Vis Image Und 2000;77(2):211–32. [15] Shuqian L, Xiang L. Implementation of mutual information based multi-modality medical image registration. Proc IEEE Eng Med Biol 2000;2:1447–50. [16] West J, Fitzpatrick JM, Wang MY, Dawant BM, Maurer CR, Kessler RM, Maciunas RJ, Barillot C, Lemoine D, Collignon A, Maes F, Suetens P, Vandermeulen D, van den Elsen PA, Napel S, Sumanaweera TS, Harkness B, Hemler PF, Hill DLG, Hawkes DJ, Studholme C, Maintz JBA, Viergever MA, Malandain G, Pennec X, Noz ME, Maguire GQ, Pollack M, Pelizzari CA, Robb RA, Hanson D, Woods RP. Comparison and evaluation of retrospective intermodality brain image registration techniques. J Comp Assist Tomogr 1997; 21(4):554–66. [17] Wang G, Vaanier MW, Skinner MW, Cavalcanti MGP, Harding GW. Spiral CT image deblurring for cochlear implantation. IEEE T Med Imaging 1998;17(2):251–61.
Sung Min Kwon received the BS and MS degrees in Electrical Engineering from KAIST, Korea in 1997 and 1999, respectively. He is working toward PhD degree in the Department of Electrical Engineering and Computer Science (EECS), KAIST, Korea. His research interests include medical image processing, 3D visualization, and volume registration.
Yong Sun Kim received the BS and MS degrees in Electrical Engineering from KAIST, Korea in 2002 and 2004, respectively. He is a PhD candidate of the Department of EECS, KAIST. His research interests include medical image processing and registration.
Tae-Sung Kim received his medical degree from Yonsei University, Seoul, Korea, in 2000. He completed his internship and residency in Diagnostic Radiology at Yonsei University Medical Center, in 2001 and 2005, respectively. His research interests are medical image processing and medical informatics.
Jong Beom Ra received the BS degree in Electronic Engineering from Seoul National University, Korea, in 1975, and the MS and PhD degrees in Electrical Engineering from KAIST, Korea, in 1977 and 1983, respectively. From 1983 to 1987, he served as an associate research scientist in Radiology, Columbia University, New York, USA. He joined the Department of EECS in KAIST, in July 1987, where he is now a professor and the director of medical imaging research center, KAIST. His research interests are digital image and video processing, medical imaging, and 3-D visualization. He is a senior member of IEEE.