J Med Syst (2014) 38:10 DOI 10.1007/s10916-014-0010-8
RESEARCH ARTICLE
Multiresolution Image Registration in Digital X-Ray Angiography with Intensity Variation Modeling Mansour Nejati & Hossein Pourghassem
Received: 13 September 2013 / Accepted: 10 January 2014 / Published online: 28 January 2014 # Springer Science+Business Media New York 2014
Abstract Digital subtraction angiography (DSA) is a widely used technique for visualization of vessel anatomy in diagnosis and treatment. However, due to unavoidable patient motions, both externally and internally, the subtracted angiography images often suffer from motion artifacts that adversely affect the quality of the medical diagnosis. To cope with this problem and improve the quality of DSA images, registration algorithms are often employed before subtraction. In this paper, a novel elastic registration algorithm for registration of digital X-ray angiography images, particularly for the coronary location, is proposed. This algorithm includes a multiresolution search strategy in which a global transformation is calculated iteratively based on local search in coarse and fine sub-image blocks. The local searches are accomplished in a differential multiscale framework which allows us to capture both large and small scale transformations. The local registration transformation also explicitly accounts for local variations in the image intensities which incorporated into our model as a change of local contrast and brightness. These local transformations are then smoothly interpolated using thin-plate spline interpolation function to obtain the global model. Experimental results with several clinical datasets demonstrate the effectiveness of our algorithm in motion artifact reduction.
Keywords Digital subtraction angiography . elastic image registration . multiresolution . thin-plate spline
M. Nejati : H. Pourghassem (*) Department of Electrical Engineering, Najafabad Branch, Islamic Azad University, 517, Najafabad, Isfahan, Iran e-mail: [email protected]
Introduction Digital subtraction angiography (DSA) is a widely used X-ray technique for vascular imaging and treatment [1, 2]. In this technique, a sequence of X-ray projection images is acquired during the passage of the injected radiopaque contrast material through the vessels of interest. The first few images of the sequence are pre-contrast or mask images taken prior to the injection of the contrast material and therefore, blood vessels are not visible in them. Other images that are acquired during the passage of the contrast material are often referred to as live images. By subtracting the mask image from each live image, interfering background structures in the live images are largely removed and the visualization of vessels is improved. However, the mask and the live images are acquired at different times, and the position of the background tissues around the vessels often changes due to the organ movement of the patient and unavoidable motions such as cardiac and respiratory motions. Thus, direct subtracting of these images would produce motion artifacts in resulting DSA images. An example of motion artifacts in the DSA images is shown in Fig. 1. These artifacts may prevent proper interpretation of the images and even lead to misdiagnosis [2]. To handle this problem and improve the diagnostic value of the DSA, registration of the mask and live images is required prior to the image subtraction. Image registration is the process of finding a transformation to align two or more images of the same scene taken at different times, from different viewpoints, and/or by different sensors [3]. Image registration has been an active area of research for more than three decades. Survey and classification of the image registration methods may be found in the literatures, such as Brown [4], Elsen et al. [5], Maintz et al. [6],
10, Page 2 of 10
J Med Syst (2014) 38:10
Fig. 1 Example of motion artifacts in angiography images. a A mask image and b a live image from a coronary angiographic image sequence. c The difference image (DSA image) that shows considerable cardiac motion artifacts
and Zitova et al. [3]. A collection of algorithms particularly suitable for registration of the medical images has been edited into a book by Hajnal et al. [7]. During the past 20 years, many efforts have been made to decrease the motion artifacts in DSA images using image registration techniques [2, 8–24]. However, the nonrigid motions of the tissues inside human body are intricate and there are local dissimilarities between the live and mask images. Thus, the registration of these images is very difficult. The most presented registration algorithms in the DSA application use a limited set of corresponding control points or landmarks to define the global transformation for aligning images [2, 8–22]. In the easiest way, the control points can be taken to lie on a regular grid [8, 19]. In this way, however, the reliability of the displacement estimations for those control points located in homogeneous regions is low. To solve this problem, an exclusion technique can be used to take into account the amount of contrast variations within a window around the control points and discard those points that are located in regions with insufficient contrast variation [10–12]. More sophisticated algorithms use image features such as vessel center lines [13], strong edges in the live or mask images [2, 14–18, 21], high-curvature points (corners) [22], or highly structured regions [20] for control point selection. After obtaining the corresponding control points of live and mask images, the parameters of registration model are estimated. This model can be a simple transformation such as affine [10–12], piecewise linear [2, 8, 21], or even a higher order model such as B-spline [9, 19], or thin-plate spline transformation [13, 15–17, 22]. Due to the large and complex nonrigid motions in the coronary or abdominal angiographic images, the aforementioned single-step algorithms cannot effectively eliminate these motion artifacts. In order to improve the image registration performance, a multiresolution registration algorithm has been developed by Yang et al. [23]. In this algorithm, the mask image is decomposed into coarse-to-fine sub-image blocks which rigidly registered to the live image. This algorithm, however, has some drawbacks. In this algorithm, Euclidean transformation is used to register the sub-image blocks that the search space is restricted to 2-D rotation and translation while a higher order model such as affine transformation can be used
to improve the performance of nonrigid image registration. Also, the optimization method used for registering the subimage blocks dose not have a closed-form, many iterations is required for convergence and easily traps in the local minimums. Moreover, in the presented algorithm in [23], if the parameters of local transformations were converged incorrectly to false values, there was not any method or strategy for detecting and compensating this shortcoming on the global transformation. To cope with these problems, we propose a new multiresolution registration algorithm. In our algorithm, the both mask and live images are decomposed to coarse and fine sub-image blocks iteratively, and search is accomplished in a multiscale framework for the affine registration of sub-image block pairs in each decomposition level [24]. The multiscale registration leads us to fast convergence and avoids the pitfalls of iterative minimization scheme. Intensity variations between mask and live images are also explicitly modeled with local changes in the brightness and the contrast of sub-image blocks. The local transformations are then smoothly interpolated using the thin-plate spline (TPS) interpolation function to obtain a global model. In this approach, the registration is refined iteratively so that the coarse-to-fine search strategy allows us to capture both large and small scale transformations. The rest of this paper is organized as follows: The proposed multiresolution registration algorithm is described in detail in “Proposed Registration Algorithm” section. It is followed by the experimental results in “Experimental Results and Discussion” section. In this section, the proposed algorithm is evaluated on simulated and real clinical datasets and then compared with two competing registration algorithms. Finally, conclusions are presented in “Conclusions” section.
Proposed Registration Algorithm In our algorithm, a multiresolution search strategy is used to register the angiography images. In this algorithm, the decomposition of the mask and live images is performed from coarse to fine level, iteratively. The registration calculations start from the coarsest level of decomposition or the first search
J Med Syst (2014) 38:10
step, where the size of sub-image block equals to the whole mask image and the geometrical difference between mask and live images is modeled coarsely with a global affine transformation. Intensity variations are also taken into account explicitly by introducing two additional parameters into our model for regarding the change of contrast and brightness. To estimate the parameters of this model, we apply mean square error (MSE) measure on the intensity values as an error function and least-squares method for minimization [25]. This procedure is accomplished upon a multiscale framework on image pyramid, which yields a fast convergence without the pitfalls of iterative minimization scheme. In multiscale estimation, a Gaussian pyramid with four levels is built for both mask and live images. The model parameters including affine, contrast, and brightness parameters simultaneously are estimated at the coarsest level. These parameters are used to warp the mask image in the next level of pyramid as an initial guess. A new estimate is computed at this level, and the process repeated through each level of the pyramid. The transformations at each level of the pyramid are accumulated to yield a final single transformation (see Fig. 2). After estimating the final affine parameters, this transformation is applied to the mask image to produce a new warped mask. This new mask image is compared with the live image using a suitable similarity measure and if the similarity is improved, the current mask image is replaced by the new warped mask image; otherwise the current mask image is retained. In the next decomposition level or the second search step, both mask and live images are decomposed to four nonoverlapping sub-image blocks and an affine registration based on aforementioned multiscale framework is performed for each pair of the corresponding blocks. Then, the calculated affine transformation is applied to each block of the mask
Fig. 2 Multiscale representation of image (Gaussian pyramid)
Page 3 of 10, 10
image. If the similarity between transformed block and its corresponding block in the live image is decreased, this calculated affine transformation will be replaced by its identity. After calculating the transformation parameters for all mask image blocks, center point of block and four points with the equal distance (a quarter of the block size) to block corners for all original blocks together with their transformed points are selected as two sets of corresponding control points for the mask and live images, respectively. The position of selected control points on the mask image in the second decomposition level is shown in Fig. 3. Given the two sets of corresponding control points, the parameters of the nonrigid global transformation are calculated and elastic warping of the mask image is performed by means of thin-plate spline (TPS) interpolation function [26]. Similar to the previous decomposition level, the similarity between the newly warped mask image and the live image is evaluated and if this similarity is improved, the current mask image will be replaced by its warped one; otherwise the current mask image is retained. The same procedure is performed in the sequential decomposition levels where at level i of decomposition (in other words, i-th search step) the mask and live images are decomposed as i×i sub-image blocks of the equal size. This strategy continues until the registration result is acceptable or the size of the sub-image block reaches to a user-specified threshold. Block diagram of our multiresolution registration algorithm with three decomposition levels is illustrated in Fig. 4. The source and target images in this block diagram refer to the mask and live images, respectively. Also, outline of updating procedure of the source image is depicted in Fig. 5. In the following, the main parts of our registration algorithm are described in detail.
Fig. 3 The position of selected control points (black dots) in the second decomposition level. Corresponding control points are obtained by applying the calculated affine transform of each block to the control points of the same block
10, Page 4 of 10
J Med Syst (2014) 38:10
Fig. 4 Overview of the proposed multiresolution registration algorithm with three decomposition levels. In this multiresolution refinement strategy, the registration result is iteratively refined toward the finer decomposition levels
Affine Registration As previously mentioned, the corresponding control point pairs at each decomposition level of registration algorithm are obtained using the affine registration of corresponding sub-image blocks in the mask and live images. To estimate the affine parameters, we employ MSE measure on the intensity values as an error function which is minimized by using least-squares method [25]. To obtain fast convergence and avoiding the pitfalls of iterative minimization method, this procedure is accomplished in a multiscale framework that allows us to capture both large and small scale transformations.
Let f(x,y,t) and f(x,y,t−1) denote the source and target images, respectively. Assume that the motion between the source and target images can be modeled by an affine transform: f ðx; y; t Þ ¼ f ðm1 x þ m2 y þ m5 ; m3 x þ m4 y þ m6 ; t−1Þ; ð1Þ
where m1, m2, m3 and m4 are the linear affine parameters, and m5, m6 are the translation parameters in x- and y-direction, respectively. To estimate these parameters, we use the following quadratic error function on the intensity values to be minimized: X E ! m ¼ ½ f ðx; y; t Þ− f ðm1 x þ m2 y þ m5 ; m3 x þ m4 y þ m6 ; t−1Þ2 x;y
ð2Þ
where ! m ¼ ðm1 ; m2 ; …; m6 ÞT is unknown parameter vector. In order to simplify the minimization, this error function is approximated using a first-order truncated Taylor series expansion as follows: X E ! m ≈ ð f ðx; y; t Þ−½ f ðx; y; t Þ þ ðm1 x þ m2 y þ m5 −xÞf x ðx; y; t Þþ x;y
ðm3 x þ m4 y þ m6 −yÞf y ðx; y; t Þ− f t ðx; y; t Þ
i2
;
ð3Þ
Fig. 5 Block diagram of the updating procedure of the source image for i-th decomposition level of our proposed algorithm
where fx (⋅) and fy (⋅) are the spatial derivatives of f(⋅) in x-direction and y-direction, and ft (⋅) is the temporal derivative of f(⋅). This approximated error can be represented more compactly as [25]:
J Med Syst (2014) 38:10
Page 5 of 10, 10
2 X T E ! m ¼ k−! c ! m ;
ð4Þ
m7 : f ðx; y; t Þ þ m8 ¼ f ðm1 x þ m2 y þ m5 ; m3 x þ m4 y þ m6 ; t−1Þ; ð7Þ
x;y
where the scalar k and vector ! c are given as: T k ¼ f t þ xf x þ y f y ; ! c ¼ xf x ; y f x ; xf y ; y f y ; f x ; f y ; ð5Þ
where m7 and m8 are two parameters to model a change in the contrast and brightness between the source and target images, respectively. Like it was previously mentioned, the MSE function is approximated with a first-order truncated Taylor series which is expressed in vector form same to Eq. (4), but, the scalar k and vector ! c are now given as: k ¼ f t −f þ xf x þ y f y ;
The error function minimization problem has a closed-form solution which is obtained by differentiating function in Eq. (4) with respect to its variable ! m and setting this result ! equal to zero. Then we have m as follows [25]: # " #−1 " X X T ! ! ! c! c ck ð6Þ m ¼ x;y
T and ! c ¼ x f x ; yf x ; x f y ; y f y ; f x ; f y ; − f ; −1 :
ð8Þ As before, the model parameters ! m are obtained by minimization of the error function E ! m . The solution takes the same form as Eq. (6), with defined k and ! c in above.
x;y
Affine Parameter Validation In practice, to estimate the affine parameters ! m , an iterative scheme is employed: Let T and ΔT denote two affine transformation matrices that the matrix T is initialized by the estimated parameters ! m from the closed-form solution in Eq. (6). At each iteration, the transformation T is applied to the original source image, and then the new affine parameters are calculated for the transformed source image and target image using Eq. (6). These new parameters are put into the matrix ΔT and the result of T×ΔT is replaced by the current affine matrix T. This process is then iterated until the affine parameters converge. Note that spatial/temporal derivatives of image in Eq. (5) affect the convergence to true affine parameters. Simple differences between neighboring pixel values are poor approximation to image derivatives. In [27], a set of filters are presented which are specifically designed for differentiation of multidimensional signals. These filters are exploited in our algorithm to compute derivatives.
In the X-ray angiographic images, we deal with nonrigid motions of tissues and local dissimilarity caused by inflow of contrast material [2]. Thus, the iterative affine registration may converge to incorrect parameters. To reduce the undesirable effect of such incorrect affine registrations on the global elastic registration, we replace the affine transformations which decrease the similarity of corresponding sub-image blocks in the mask and live images by identity transformation. We use entropy of difference image (EDI) as a similarity measure for validating the affine parameters which is the most suitable similarity measure for the DSA registration according to the presented quantitative evaluation in [20]. The similarity between two images A and B in EDI measure is defined as the entropy of normalized gray-value histogram of the difference image: S ðA; BÞ ¼ −
g max X pg log pg
ð9Þ
g¼gmin
Intensity Variations In calculating the parameters of model presented in the previous section, we assumed that the image intensities between the source and target images are unchanged, i.e., the assumption of brightness constancy. This assumption is likely to fail in the case of X-ray angiography images. This fact is mainly due to the changes of local densities as a result of the contractions and expansions of tissues, and also fluctuations in the intensity of X-rays [18]. To account for intensity variations, an explicit change of contrast and brightness is incorporated into our model:
where gmin and gmax are the minimum and maximum grayvalues of difference image, respectively, and pg is the fraction of pixels in this image with gray-value g. When the similarity between two images is increased, the gray-value histogram of difference image becomes more clustered and, therefore, its entropy is decreased. Thin-Plate Spline Warping After calculating the affine parameters for all sub-image blocks in each decomposition level, two sets of the
10, Page 6 of 10
J Med Syst (2014) 38:10
corresponding control points, P={pi} with pi=(xi,yi), i=1,2, …,n, and Q={qi} with qi=(x′i,y′i), i=1,2,…,n, are obtained. The elastic warping of the mask image is then performed by applying the TPS interpolation function on these two sets of corresponding control points. In other words, in order to carry out the warping of the mask image, it is required to have a dense correspondence between the mask and live images which can be estimated through interpolation of the corresponding control points P={pi} and Q={qi}. The TPS is a well-known method to smooth interpolation of bivariate scattered data which was first proposed as a pointbased elastic registration algorithm for medical images by Bookstein [26]. The TPS interpolation function is defined from two sets of n corresponding control pairs as
x0 0 y
¼
a11 a21
t x a12 þ x þ : ty y a22
Xn Xni¼1 i¼1
wx;i ϕðkðxi ; yi Þ−ðx; yÞkÞ wy;i ϕðkðxi ; yi Þ−ðx; yÞkÞ
! ;
ð10Þ
where a11, a12, a21, a22, tx and ty represent the global affine transformation parameters, ‖.‖ denotes the Euclidean norm, φ(r)=r2 ln(r) is the TPS radial basis function, and parameters wx,i and wy,i representing the weights of TPS function φ. Given the two sets P and Q of n corresponding control point pairs, these parameters can be obtained by solving the following linear system: K P W Q ¼ ð11Þ PT O a o
where K is a n×n matrix with entries Kij=φ(||(xi,yi)−(xj,yj)||), O and o are 3×3 and 3×2 matrices of zeros, respectively. P is n× 3 matrix which its i-th row is (1,xi,yi), Q is a n×2 matrix formed from coordinate of the control points qi, and the coefficient matrices W and a are defined as 2 3 2 3 wx;1 wy;1 ty tx ⋮ 5; a ¼ 4 a11 a21 5 W ¼4 ⋮ ð12Þ wx;n wy;n a12 a22
The complete set of parameters, defining the TPS registration transformation, is then used to map each point (x,y) in the mask image to the corresponding point (x′,y′) in the live image and consequently the elastic warping of the mask image.
Experimental Results and Discussion Our registration algorithm was implemented in Matlab and tested on AMD processor 2.3GHz with 1GB RAM. To
evaluate the performance of the proposed registration algorithm in enhancing the DSA image sequences, we carried out experiments on various datasets. All datasets that were used in the experiments are clinical coronary angiography image sequences, acquired from a Siemens AXIOM Artis FC DSA imaging system consisted of 512×512 images with grayvalue resolution of 8 bits. Registration Results on Simulated Images For quantitative evaluation of the proposed registration algorithm, firstly, it is employed to register simulated image pairs. Each test dataset consists of a live image taken from a cardiac angiography image sequence and its elastically warped version as a simulated live image. To generate a simulated live image, we selected a set of the control points randomly distributed on the live image. For each control point, an uncorrelated random shift to both x- and y- direction ranging from −10 to 10 pixels was added and a new set of random control points were thus formed. Using the generated control point pairs and TPS interpolation method described in “Thin-Plate Spline Warping” section, a new simulated live image was produced by elastic warping of the live image. In this way, the true transformation parameters and therefore true correspondences between these two images are known. Thus, by calculating the root-mean-squares (RMS) error of the control point pairs after registration of the simulated image, the algorithm can be evaluated objectively. The RMS error, δrms, is computed as: !!0:5 n 2 1 X
0 0
δrms ¼ ð13Þ
xi ; yi − bxi ; byi n i¼1
where ðbxi ; byi Þ , i=1,2,…,n are corresponding points of a set of n randomly selected control points on the simulated live image which are calculated using the known TPS transformation, and (x′i,y′i), i=1,2,…,n, are corresponding points estimated after registration of the live and simulated live images. Moreover, the registration performance may be evaluated according to the RMS error reduction ratio that is defined as, R¼
δ0 −δ % δ0
ð14Þ
where δ0 denotes the initial RMS error between the corresponding control points of the live and simulated images, and δ is the value of RMS error after registration. We applied our proposed registration algorithm to several simulated image pairs and compared the registration results with the results of two competing algorithms introduced
J Med Syst (2014) 38:10
Page 7 of 10, 10
Table 1 The comparison between the proposed registration algorithm and other algorithms on ten simulated image pairs Image Pairs
No. 1 No. 2 No. 3 No. 4 No. 5 No. 6 No. 7 No. 8 No. 9 No. 10 Average
Initial error before registration
After registration Single step algorithm [2]
Multiresolution algorithm [23]
Proposed algorithm
RMS (δ0)
RMS
R%
Elapsed time (s)
RMS
R%
Elapsed time (s)
RMS
R%
6.098 6.892 6.442 6.201 6.363 6.270 6.813 6.920 7.399 7.494 –
3.903 2.725 2.659 3.621 3.965 2.233 3.717 2.995 3.898 3.266 –
36.0 60.5 58.7 41.6 37.7 64.4 45.4 56.7 47.3 56.4 50.5
34.4 33.1 34.8 35.9 33.4 35.8 33.5 36.1 33.9 33.5 34.4
1.511 1.015 0.870 1.603 1.254 1.542 1.532 1.272 1.784 2.331 –
75.2 85.3 86.5 74.1 80.3 75.4 77.5 81.6 75.9 68.9 78.1
8954 9958 8745 9180 9553 9210 8768 9003 9548 9346 9226
1.124 0.781 0.875 0.882 1.154 1.316 1.191 1.020 1.337 1.037 –
81.6 88.7 86.4 85.8 81.9 79.0 82.5 85.2 81.9 86.2 83.9
% % % % % % % % % % %
respectively by Meijering et al. [2], and Yang et al. [23]. The first one [2] is a single step algorithm in which a set of control points are extracted from strong edges of the mask image. The displacement of each control point is then computed by means of the template matching using the energy of difference histogram as similarity measure. Also, a hill-climbing approach is used as an optimization method. At last, the image warping is performed using a Delaunay triangulation and piecewise linear interpolation. The second algorithm [23] is a multiresolution registration algorithm in which the mask image is decomposed to coarse and fine sub-image blocks Fig. 6 Example of the registration results on simulated images. a Original live image from a coronary angiography sequence. b Simulated live image. It is warped from the original live image with 30 random control point pairs using thin-plate spline. c Subtraction image corresponding to the difference between a and b. d Subtraction image after registration by single step algorithm [2]. e Subtraction image after registration by multiresolution algorithm [23]. f Subtraction image after registration by our proposed registration algorithm
% % % % % % % % % % %
Elapsed time (s) % % % % % % % % % % %
285 336 319 323 304 308 325 322 314 305 314
iteratively and each block is rigidly registered to the live image using Powell optimization method [28]. The calculated registration transform in each decomposition level is used as an initial guess for the next finer level. The comparison between our proposed algorithm and the single step [2] and multiresolution [23] algorithms is accomplished according to the registration performance and computational time. The performance is evaluated based on RMS pixel error between the true corresponding control points and the ones estimated by registration of the live and simulated live image pairs. Table 1 lists the achieved registration results
10, Page 8 of 10 Table 2 Evaluation of affine parameter validation procedure on the registration performance
J Med Syst (2014) 38:10
Test image pairs
No. 1 No. 2 No. 3 No. 4 No. 5 No. 6 No. 7 No. 8 No. 9 No. 10 Average
Registration results of proposed algorithm With affine parameter validation
Without affine parameter validation
RMS
R%
RMS
R%
1.124 0.781 0.875 0.882 1.154 1.316 1.191 1.020 1.337 1.037 –
81.6 88.7 86.4 85.8 81.9 79.0 82.5 85.2 81.9 86.2 83.9
1.208 2.191 1.364 1.956 1.669 1.729 1.743 2.114 1.682 1.450 –
80.19 68.21 78.83 68.46 73.77 72.42 74.41 69.45 77.26 80.65 74.37
by our proposed algorithm and two presented algorithms in [2] and [23] for the ten simulated image pairs. The results of our proposed algorithm in Table 1 are obtained at level 8 of the image decomposition. For the ten test image pairs, the reduction ratio of RMS error (R) of 83.9 % has been obtained by our proposed registration algorithm that is better than the corresponding value for multiresolution algorithm in [23], while our algorithm has considerably lower computational time. The computational time of the single step algorithm in [2] is lower than computational time of our proposed algorithm, however, the single step algorithm results in poor registration performance and our algorithm considerably outperforms for all image pairs. Figure 6 shows the
Fig. 7 Multiresolution registration results for three coronary angiography images. Columns from left to right correspond to mask images, live images, the original subtraction of the respective mask and the live images, and the subtraction after correction of motion artifacts using the proposed registration algorithm with eight decomposition levels, respectively
% % % % % % % % % % %
% % % % % % % % % % %
registration results of these algorithms for one of the simulated image pairs. As, it can be seen in Fig. 6d, e and f, the extensive artifacts still remain in the subtraction images when single step algorithm [2] or multiresolution algorithm [23] are applied. As previously mentioned in “Affine Parameter Validation” section, in our application, we employ an affine parameter validation procedure to reduce the undesirable effect of incorrect affine registration of blocks on the global elastic transformation. To evaluate this procedure on the registration performance, we also tested our algorithm with and without affine parameter validation on the simulated image pairs. Table 2 demonstrates the significant improvement in the registration results by using affine parameter validation.
J Med Syst (2014) 38:10
Registration Results on Clinical Images We carried out several experiments on clinical datasets to evaluate the performance of proposed multiresolution registration algorithm in the motion artifact reduction of DSA images. By several experiments on the coronary angiography images with different levels of decomposition, it was concluded that the best results are often achieved in the level 8 where the size of sub-image blocks is equal to 64×64 pixels. In the multiscale affine registration of sub-image blocks, a Gaussian pyramid is constructed for both the source and target images which the number of pyramid levels is selected such that the size of the coarsest level is greater than or equal to 64×64 pixels (because the coarser levels significantly lose the gradient patterns of original image which can degrade initial affine estimation). Figure 7 shows the multiresolution registration results for the three cardiac angiography images of datasets. In this figure, columns from left to right correspond to the mask images, the live images, original subtraction of the mask and live images, and the subtraction after correction of motion artifacts, using the proposed multiresolution registration algorithm, respectively. These results demonstrate that the motion artifacts in the original DSA images are largely removed after registration.
Conclusions In this paper, a new elastic registration algorithm for motion artifact reduction in DSA was presented. This method is based on a multiresolution search strategy in which both the mask and live images are decomposed to coarse and fine sub-image blocks iteratively. All the corresponding block pairs at each decomposition level are affinely registered using least-squares optimization which for fast convergence and avoiding the pitfalls of iterative optimization method, the entire procedure is done in a multiscale framework. To account for intensity variations, an explicit change of the contrast and brightness is also incorporated into our block registration model. Furthermore, in order to detect the probable incorrect registrations, a validation procedure is performed after affine registration of each corresponding block pairs based on the similarity assessment. The estimated affine parameters of registered sub-image blocks in each decomposition level are then smoothly interpolated by mean of the thin-plate spline interpolation function to generate a global transformation at each level. To evaluate our algorithm, we tested it with various clinical and simulated datasets. The obtained registration results demonstrate that our proposed algorithm outperforms other presented algorithms in the literature in terms of visual and quantitative evaluations. The experiments on the coronary angiography images also show the effectiveness of our method in the motion artifact reduction.
Page 9 of 10, 10 Acknowledgments The authors would like to thank Dr. Hashemi of Sina Heart Hospital in Isfahan, Iran, for providing us the datasets used in the experiments.
References 1. Brody, W. R., Digital subtraction angiography. IEEE Trans. Nucl. Sci. 29(3):1176–1180, 1982. 2. Meijering, E. H. W., Zuiderveld, K. J., and Viergever, M. A., Image registration for digital subtraction angiography. Int. J. Comput. Vis. 31(2/3):227–246, 1999. 3. Zitova, B., and Flusser, J., Image registration methods: A survey. Image Vis. Comput. 21(11):977–1000, 2003. 4. Brown, L. G., A survey of image registration techniques. ACM Comput. Surv. 24(4):325–276, 1992. 5. Elsen, V. D., Pol, P. A., and Viergever, E.-J. D., Medical image matching: A review with classification. IEEE Eng. Med. Biol. Magazine 12(1):26–39, 1993. 6. Maintz, J. B. A., and Viergever, M. A., A survey of medical image registration. Med. Image Anal. 2(1):1–36, 1998. 7. Hajnal, J. V., Hill, D. L. G., and Hawkes, D. J., Medical image registration. CRS Press, Baton Rouge, Florida, 2001. 8. Cox, G. S., De Jager, G. “Automatic registration of temporal image pairs for digital subtraction angiography,” In: Loew M. H. (ed.), Proceedings of the SPIE, Medical imaging: image processing, vol. 2167, pp. 188–199, Bellingham, WA, 1994. 9. Nejati, M., Sadri, S., and Amirfattahi, R., “Nonrigid Image Registration in Digital Subtraction Angiography Using Multilevel B-Spline”, BioMed Research International, pp. 1–12, vol. 2013. 10. Buzug, T. M., Weese, J., Fassnacht, C., and Lorenz, C., Using an entropy similarity measure to enhance the quality of DSA images with an algorithm based on template matching. In: Höhne, K.-H., and Kikinis, R. (Eds.), VBC’96, Lecture Notes in Computer Science, vol. 1131. Springer, Berlin, pp. 235–240, 1996. 11. Buzug, T. M., Weese, J., Lorenz, C., and Beil, W., Histogram-based image registration for digital subtraction angiography. In: Del Bimbo, A. (Ed.), ICIAP’97, Lecture Notes in Computer Science, vol. 1311. Springer, Berlin, pp. 380–387, 1997. 12. Buzug, T. M., and Weese, J., Image registration for DSA quality enhancement. Comput. Med. Imaging Graph. 22(2):103–113, 1998. 13. Cao, Z., Liu, X., Peng, B., Moon, Y.-S. “DSA image registration based on multiscale Gabor filters and mutual information”. In: Proceedings of the IEEE ICIA’05, Hong Kong and Macau, China, pp. 105–110, 2005. 14. Taleb, N., and Jetto, L., Image registration for applications in digital subtraction angiography. Control. Eng. Pract. 6(2):227–238, 1998. 15. Bentoutou, Y., Taleb, N., Chikr El Mezouar, M., Taleb, M., and Jetto, L., An invariant approach for image registration in digital subtraction angiography. Pattern Recogn. 35(12):2853–2865, 2002. 16. Bentoutou, Y., and Taleb, N., Automatic extraction of control points for digital subtraction angiography image enhancement. IEEE Trans. Nucl. Sci. 52(1):238–246, 2005. 17. Bentoutou, Y., and Taleb, N., A 3-D space-time motion detection for an invariant image registration approach in digital subtraction angiography. Comput. Vis. Image Underst. 97(1):30–50, 2005. 18. Meijering, E. H. W., Niessen, W. J., and Viergever, M. A., Retrospective motion correction in digital subtraction angiography: a review. IEEE Trans. Med. Imaging 18(1):2–21, 1999. 19. Ping, L., Hong, N., Ye, S. “An efficient method for image registration in DSA,” In: Proceedings of the ISISE’08, Shanghai, China, vol. 2, pp. 551–554, 2008. 20. Nejati, M., Amirfattahi, R., Sadri, S.“A fast image registration algorithm for digital subtraction angiography”, In: Proceedings of 17th
10, Page 10 of 10
21.
22.
23.
24.
Iranian Conference of Biomedical Engineering (ICBME2010), Isfahan, Iran, pp. 1–4, November 2010. Chu, Y., Bai, N., Ji, Z., Chen, S., Mou, X. “Registration for DSA image using triangle grid and spatial transformation based on stretching”, In: Proceedings of the 8th International Conference on Signal Processing, Beijing, China, 16–20 Nov. 2006. Wang, J., Zhang, J. Q. “An iterative refinement DSA image registration algorithm using structural image quality measure”, In: Proceedings of the 5th IIH-MSP, Kyoto, Japan, pp. 973–976, 2009. Yang, J., Wang, Y., Tang, S., Zhou, S., Liu, Y., and Chen, W., Multiresolution elastic registration of X-ray angiography images using thin-plate spline. IEEE Trans. Nucl. Sci. 54(1):152–166, 2007. Nejati, M., Pourghassem, H. “Multiresolution search strategy for elastic registration of x-ray angiography images”, In Proceedings of International Conference on Intelligent Computation and Bio-
J Med Syst (2014) 38:10
25.
26.
27.
28.
Medical Instrumentation (ICBMI 2011), pp. 216–219, Wuhan, China, 14–17 Dec. 2011. Periaswamy, S., and Farid, H., Elastic registration in the Presence of Intensity Variations. IEEE Trans. Med. Imaging 22(7):865–874, 2003. Bookstein, F. L., Principal warps: thin-plate splines and the decomposition of deformations. IEEE Transact. Pattern Anal. Mach. Intell. 11(6):567–585, 1989. Farid, H., Simoncelli, E. P.“Optimally rotation-equivariant directional derivative kernels”, In Proceedings of International Conference on Computer Analysis of Images and Patterns, pp. 207–214, Berlin, Germany, Sept. 1997. Powell, M. J. D., An efficient method for finding the minimum of a function of several variables without calculating derivatives. Comput. J. 7(2):155–162, 1964.
RESEARCH ARTICLE
Multiresolution Image Registration in Digital X-Ray Angiography with Intensity Variation Modeling Mansour Nejati & Hossein Pourghassem
Received: 13 September 2013 / Accepted: 10 January 2014 / Published online: 28 January 2014 # Springer Science+Business Media New York 2014
Abstract Digital subtraction angiography (DSA) is a widely used technique for visualization of vessel anatomy in diagnosis and treatment. However, due to unavoidable patient motions, both externally and internally, the subtracted angiography images often suffer from motion artifacts that adversely affect the quality of the medical diagnosis. To cope with this problem and improve the quality of DSA images, registration algorithms are often employed before subtraction. In this paper, a novel elastic registration algorithm for registration of digital X-ray angiography images, particularly for the coronary location, is proposed. This algorithm includes a multiresolution search strategy in which a global transformation is calculated iteratively based on local search in coarse and fine sub-image blocks. The local searches are accomplished in a differential multiscale framework which allows us to capture both large and small scale transformations. The local registration transformation also explicitly accounts for local variations in the image intensities which incorporated into our model as a change of local contrast and brightness. These local transformations are then smoothly interpolated using thin-plate spline interpolation function to obtain the global model. Experimental results with several clinical datasets demonstrate the effectiveness of our algorithm in motion artifact reduction.
Keywords Digital subtraction angiography . elastic image registration . multiresolution . thin-plate spline
M. Nejati : H. Pourghassem (*) Department of Electrical Engineering, Najafabad Branch, Islamic Azad University, 517, Najafabad, Isfahan, Iran e-mail: [email protected]
Introduction Digital subtraction angiography (DSA) is a widely used X-ray technique for vascular imaging and treatment [1, 2]. In this technique, a sequence of X-ray projection images is acquired during the passage of the injected radiopaque contrast material through the vessels of interest. The first few images of the sequence are pre-contrast or mask images taken prior to the injection of the contrast material and therefore, blood vessels are not visible in them. Other images that are acquired during the passage of the contrast material are often referred to as live images. By subtracting the mask image from each live image, interfering background structures in the live images are largely removed and the visualization of vessels is improved. However, the mask and the live images are acquired at different times, and the position of the background tissues around the vessels often changes due to the organ movement of the patient and unavoidable motions such as cardiac and respiratory motions. Thus, direct subtracting of these images would produce motion artifacts in resulting DSA images. An example of motion artifacts in the DSA images is shown in Fig. 1. These artifacts may prevent proper interpretation of the images and even lead to misdiagnosis [2]. To handle this problem and improve the diagnostic value of the DSA, registration of the mask and live images is required prior to the image subtraction. Image registration is the process of finding a transformation to align two or more images of the same scene taken at different times, from different viewpoints, and/or by different sensors [3]. Image registration has been an active area of research for more than three decades. Survey and classification of the image registration methods may be found in the literatures, such as Brown [4], Elsen et al. [5], Maintz et al. [6],
10, Page 2 of 10
J Med Syst (2014) 38:10
Fig. 1 Example of motion artifacts in angiography images. a A mask image and b a live image from a coronary angiographic image sequence. c The difference image (DSA image) that shows considerable cardiac motion artifacts
and Zitova et al. [3]. A collection of algorithms particularly suitable for registration of the medical images has been edited into a book by Hajnal et al. [7]. During the past 20 years, many efforts have been made to decrease the motion artifacts in DSA images using image registration techniques [2, 8–24]. However, the nonrigid motions of the tissues inside human body are intricate and there are local dissimilarities between the live and mask images. Thus, the registration of these images is very difficult. The most presented registration algorithms in the DSA application use a limited set of corresponding control points or landmarks to define the global transformation for aligning images [2, 8–22]. In the easiest way, the control points can be taken to lie on a regular grid [8, 19]. In this way, however, the reliability of the displacement estimations for those control points located in homogeneous regions is low. To solve this problem, an exclusion technique can be used to take into account the amount of contrast variations within a window around the control points and discard those points that are located in regions with insufficient contrast variation [10–12]. More sophisticated algorithms use image features such as vessel center lines [13], strong edges in the live or mask images [2, 14–18, 21], high-curvature points (corners) [22], or highly structured regions [20] for control point selection. After obtaining the corresponding control points of live and mask images, the parameters of registration model are estimated. This model can be a simple transformation such as affine [10–12], piecewise linear [2, 8, 21], or even a higher order model such as B-spline [9, 19], or thin-plate spline transformation [13, 15–17, 22]. Due to the large and complex nonrigid motions in the coronary or abdominal angiographic images, the aforementioned single-step algorithms cannot effectively eliminate these motion artifacts. In order to improve the image registration performance, a multiresolution registration algorithm has been developed by Yang et al. [23]. In this algorithm, the mask image is decomposed into coarse-to-fine sub-image blocks which rigidly registered to the live image. This algorithm, however, has some drawbacks. In this algorithm, Euclidean transformation is used to register the sub-image blocks that the search space is restricted to 2-D rotation and translation while a higher order model such as affine transformation can be used
to improve the performance of nonrigid image registration. Also, the optimization method used for registering the subimage blocks dose not have a closed-form, many iterations is required for convergence and easily traps in the local minimums. Moreover, in the presented algorithm in [23], if the parameters of local transformations were converged incorrectly to false values, there was not any method or strategy for detecting and compensating this shortcoming on the global transformation. To cope with these problems, we propose a new multiresolution registration algorithm. In our algorithm, the both mask and live images are decomposed to coarse and fine sub-image blocks iteratively, and search is accomplished in a multiscale framework for the affine registration of sub-image block pairs in each decomposition level [24]. The multiscale registration leads us to fast convergence and avoids the pitfalls of iterative minimization scheme. Intensity variations between mask and live images are also explicitly modeled with local changes in the brightness and the contrast of sub-image blocks. The local transformations are then smoothly interpolated using the thin-plate spline (TPS) interpolation function to obtain a global model. In this approach, the registration is refined iteratively so that the coarse-to-fine search strategy allows us to capture both large and small scale transformations. The rest of this paper is organized as follows: The proposed multiresolution registration algorithm is described in detail in “Proposed Registration Algorithm” section. It is followed by the experimental results in “Experimental Results and Discussion” section. In this section, the proposed algorithm is evaluated on simulated and real clinical datasets and then compared with two competing registration algorithms. Finally, conclusions are presented in “Conclusions” section.
Proposed Registration Algorithm In our algorithm, a multiresolution search strategy is used to register the angiography images. In this algorithm, the decomposition of the mask and live images is performed from coarse to fine level, iteratively. The registration calculations start from the coarsest level of decomposition or the first search
J Med Syst (2014) 38:10
step, where the size of sub-image block equals to the whole mask image and the geometrical difference between mask and live images is modeled coarsely with a global affine transformation. Intensity variations are also taken into account explicitly by introducing two additional parameters into our model for regarding the change of contrast and brightness. To estimate the parameters of this model, we apply mean square error (MSE) measure on the intensity values as an error function and least-squares method for minimization [25]. This procedure is accomplished upon a multiscale framework on image pyramid, which yields a fast convergence without the pitfalls of iterative minimization scheme. In multiscale estimation, a Gaussian pyramid with four levels is built for both mask and live images. The model parameters including affine, contrast, and brightness parameters simultaneously are estimated at the coarsest level. These parameters are used to warp the mask image in the next level of pyramid as an initial guess. A new estimate is computed at this level, and the process repeated through each level of the pyramid. The transformations at each level of the pyramid are accumulated to yield a final single transformation (see Fig. 2). After estimating the final affine parameters, this transformation is applied to the mask image to produce a new warped mask. This new mask image is compared with the live image using a suitable similarity measure and if the similarity is improved, the current mask image is replaced by the new warped mask image; otherwise the current mask image is retained. In the next decomposition level or the second search step, both mask and live images are decomposed to four nonoverlapping sub-image blocks and an affine registration based on aforementioned multiscale framework is performed for each pair of the corresponding blocks. Then, the calculated affine transformation is applied to each block of the mask
Fig. 2 Multiscale representation of image (Gaussian pyramid)
Page 3 of 10, 10
image. If the similarity between transformed block and its corresponding block in the live image is decreased, this calculated affine transformation will be replaced by its identity. After calculating the transformation parameters for all mask image blocks, center point of block and four points with the equal distance (a quarter of the block size) to block corners for all original blocks together with their transformed points are selected as two sets of corresponding control points for the mask and live images, respectively. The position of selected control points on the mask image in the second decomposition level is shown in Fig. 3. Given the two sets of corresponding control points, the parameters of the nonrigid global transformation are calculated and elastic warping of the mask image is performed by means of thin-plate spline (TPS) interpolation function [26]. Similar to the previous decomposition level, the similarity between the newly warped mask image and the live image is evaluated and if this similarity is improved, the current mask image will be replaced by its warped one; otherwise the current mask image is retained. The same procedure is performed in the sequential decomposition levels where at level i of decomposition (in other words, i-th search step) the mask and live images are decomposed as i×i sub-image blocks of the equal size. This strategy continues until the registration result is acceptable or the size of the sub-image block reaches to a user-specified threshold. Block diagram of our multiresolution registration algorithm with three decomposition levels is illustrated in Fig. 4. The source and target images in this block diagram refer to the mask and live images, respectively. Also, outline of updating procedure of the source image is depicted in Fig. 5. In the following, the main parts of our registration algorithm are described in detail.
Fig. 3 The position of selected control points (black dots) in the second decomposition level. Corresponding control points are obtained by applying the calculated affine transform of each block to the control points of the same block
10, Page 4 of 10
J Med Syst (2014) 38:10
Fig. 4 Overview of the proposed multiresolution registration algorithm with three decomposition levels. In this multiresolution refinement strategy, the registration result is iteratively refined toward the finer decomposition levels
Affine Registration As previously mentioned, the corresponding control point pairs at each decomposition level of registration algorithm are obtained using the affine registration of corresponding sub-image blocks in the mask and live images. To estimate the affine parameters, we employ MSE measure on the intensity values as an error function which is minimized by using least-squares method [25]. To obtain fast convergence and avoiding the pitfalls of iterative minimization method, this procedure is accomplished in a multiscale framework that allows us to capture both large and small scale transformations.
Let f(x,y,t) and f(x,y,t−1) denote the source and target images, respectively. Assume that the motion between the source and target images can be modeled by an affine transform: f ðx; y; t Þ ¼ f ðm1 x þ m2 y þ m5 ; m3 x þ m4 y þ m6 ; t−1Þ; ð1Þ
where m1, m2, m3 and m4 are the linear affine parameters, and m5, m6 are the translation parameters in x- and y-direction, respectively. To estimate these parameters, we use the following quadratic error function on the intensity values to be minimized: X E ! m ¼ ½ f ðx; y; t Þ− f ðm1 x þ m2 y þ m5 ; m3 x þ m4 y þ m6 ; t−1Þ2 x;y
ð2Þ
where ! m ¼ ðm1 ; m2 ; …; m6 ÞT is unknown parameter vector. In order to simplify the minimization, this error function is approximated using a first-order truncated Taylor series expansion as follows: X E ! m ≈ ð f ðx; y; t Þ−½ f ðx; y; t Þ þ ðm1 x þ m2 y þ m5 −xÞf x ðx; y; t Þþ x;y
ðm3 x þ m4 y þ m6 −yÞf y ðx; y; t Þ− f t ðx; y; t Þ
i2
;
ð3Þ
Fig. 5 Block diagram of the updating procedure of the source image for i-th decomposition level of our proposed algorithm
where fx (⋅) and fy (⋅) are the spatial derivatives of f(⋅) in x-direction and y-direction, and ft (⋅) is the temporal derivative of f(⋅). This approximated error can be represented more compactly as [25]:
J Med Syst (2014) 38:10
Page 5 of 10, 10
2 X T E ! m ¼ k−! c ! m ;
ð4Þ
m7 : f ðx; y; t Þ þ m8 ¼ f ðm1 x þ m2 y þ m5 ; m3 x þ m4 y þ m6 ; t−1Þ; ð7Þ
x;y
where the scalar k and vector ! c are given as: T k ¼ f t þ xf x þ y f y ; ! c ¼ xf x ; y f x ; xf y ; y f y ; f x ; f y ; ð5Þ
where m7 and m8 are two parameters to model a change in the contrast and brightness between the source and target images, respectively. Like it was previously mentioned, the MSE function is approximated with a first-order truncated Taylor series which is expressed in vector form same to Eq. (4), but, the scalar k and vector ! c are now given as: k ¼ f t −f þ xf x þ y f y ;
The error function minimization problem has a closed-form solution which is obtained by differentiating function in Eq. (4) with respect to its variable ! m and setting this result ! equal to zero. Then we have m as follows [25]: # " #−1 " X X T ! ! ! c! c ck ð6Þ m ¼ x;y
T and ! c ¼ x f x ; yf x ; x f y ; y f y ; f x ; f y ; − f ; −1 :
ð8Þ As before, the model parameters ! m are obtained by minimization of the error function E ! m . The solution takes the same form as Eq. (6), with defined k and ! c in above.
x;y
Affine Parameter Validation In practice, to estimate the affine parameters ! m , an iterative scheme is employed: Let T and ΔT denote two affine transformation matrices that the matrix T is initialized by the estimated parameters ! m from the closed-form solution in Eq. (6). At each iteration, the transformation T is applied to the original source image, and then the new affine parameters are calculated for the transformed source image and target image using Eq. (6). These new parameters are put into the matrix ΔT and the result of T×ΔT is replaced by the current affine matrix T. This process is then iterated until the affine parameters converge. Note that spatial/temporal derivatives of image in Eq. (5) affect the convergence to true affine parameters. Simple differences between neighboring pixel values are poor approximation to image derivatives. In [27], a set of filters are presented which are specifically designed for differentiation of multidimensional signals. These filters are exploited in our algorithm to compute derivatives.
In the X-ray angiographic images, we deal with nonrigid motions of tissues and local dissimilarity caused by inflow of contrast material [2]. Thus, the iterative affine registration may converge to incorrect parameters. To reduce the undesirable effect of such incorrect affine registrations on the global elastic registration, we replace the affine transformations which decrease the similarity of corresponding sub-image blocks in the mask and live images by identity transformation. We use entropy of difference image (EDI) as a similarity measure for validating the affine parameters which is the most suitable similarity measure for the DSA registration according to the presented quantitative evaluation in [20]. The similarity between two images A and B in EDI measure is defined as the entropy of normalized gray-value histogram of the difference image: S ðA; BÞ ¼ −
g max X pg log pg
ð9Þ
g¼gmin
Intensity Variations In calculating the parameters of model presented in the previous section, we assumed that the image intensities between the source and target images are unchanged, i.e., the assumption of brightness constancy. This assumption is likely to fail in the case of X-ray angiography images. This fact is mainly due to the changes of local densities as a result of the contractions and expansions of tissues, and also fluctuations in the intensity of X-rays [18]. To account for intensity variations, an explicit change of contrast and brightness is incorporated into our model:
where gmin and gmax are the minimum and maximum grayvalues of difference image, respectively, and pg is the fraction of pixels in this image with gray-value g. When the similarity between two images is increased, the gray-value histogram of difference image becomes more clustered and, therefore, its entropy is decreased. Thin-Plate Spline Warping After calculating the affine parameters for all sub-image blocks in each decomposition level, two sets of the
10, Page 6 of 10
J Med Syst (2014) 38:10
corresponding control points, P={pi} with pi=(xi,yi), i=1,2, …,n, and Q={qi} with qi=(x′i,y′i), i=1,2,…,n, are obtained. The elastic warping of the mask image is then performed by applying the TPS interpolation function on these two sets of corresponding control points. In other words, in order to carry out the warping of the mask image, it is required to have a dense correspondence between the mask and live images which can be estimated through interpolation of the corresponding control points P={pi} and Q={qi}. The TPS is a well-known method to smooth interpolation of bivariate scattered data which was first proposed as a pointbased elastic registration algorithm for medical images by Bookstein [26]. The TPS interpolation function is defined from two sets of n corresponding control pairs as
x0 0 y
¼
a11 a21
t x a12 þ x þ : ty y a22
Xn Xni¼1 i¼1
wx;i ϕðkðxi ; yi Þ−ðx; yÞkÞ wy;i ϕðkðxi ; yi Þ−ðx; yÞkÞ
! ;
ð10Þ
where a11, a12, a21, a22, tx and ty represent the global affine transformation parameters, ‖.‖ denotes the Euclidean norm, φ(r)=r2 ln(r) is the TPS radial basis function, and parameters wx,i and wy,i representing the weights of TPS function φ. Given the two sets P and Q of n corresponding control point pairs, these parameters can be obtained by solving the following linear system: K P W Q ¼ ð11Þ PT O a o
where K is a n×n matrix with entries Kij=φ(||(xi,yi)−(xj,yj)||), O and o are 3×3 and 3×2 matrices of zeros, respectively. P is n× 3 matrix which its i-th row is (1,xi,yi), Q is a n×2 matrix formed from coordinate of the control points qi, and the coefficient matrices W and a are defined as 2 3 2 3 wx;1 wy;1 ty tx ⋮ 5; a ¼ 4 a11 a21 5 W ¼4 ⋮ ð12Þ wx;n wy;n a12 a22
The complete set of parameters, defining the TPS registration transformation, is then used to map each point (x,y) in the mask image to the corresponding point (x′,y′) in the live image and consequently the elastic warping of the mask image.
Experimental Results and Discussion Our registration algorithm was implemented in Matlab and tested on AMD processor 2.3GHz with 1GB RAM. To
evaluate the performance of the proposed registration algorithm in enhancing the DSA image sequences, we carried out experiments on various datasets. All datasets that were used in the experiments are clinical coronary angiography image sequences, acquired from a Siemens AXIOM Artis FC DSA imaging system consisted of 512×512 images with grayvalue resolution of 8 bits. Registration Results on Simulated Images For quantitative evaluation of the proposed registration algorithm, firstly, it is employed to register simulated image pairs. Each test dataset consists of a live image taken from a cardiac angiography image sequence and its elastically warped version as a simulated live image. To generate a simulated live image, we selected a set of the control points randomly distributed on the live image. For each control point, an uncorrelated random shift to both x- and y- direction ranging from −10 to 10 pixels was added and a new set of random control points were thus formed. Using the generated control point pairs and TPS interpolation method described in “Thin-Plate Spline Warping” section, a new simulated live image was produced by elastic warping of the live image. In this way, the true transformation parameters and therefore true correspondences between these two images are known. Thus, by calculating the root-mean-squares (RMS) error of the control point pairs after registration of the simulated image, the algorithm can be evaluated objectively. The RMS error, δrms, is computed as: !!0:5 n 2 1 X
0 0
δrms ¼ ð13Þ
xi ; yi − bxi ; byi n i¼1
where ðbxi ; byi Þ , i=1,2,…,n are corresponding points of a set of n randomly selected control points on the simulated live image which are calculated using the known TPS transformation, and (x′i,y′i), i=1,2,…,n, are corresponding points estimated after registration of the live and simulated live images. Moreover, the registration performance may be evaluated according to the RMS error reduction ratio that is defined as, R¼
δ0 −δ % δ0
ð14Þ
where δ0 denotes the initial RMS error between the corresponding control points of the live and simulated images, and δ is the value of RMS error after registration. We applied our proposed registration algorithm to several simulated image pairs and compared the registration results with the results of two competing algorithms introduced
J Med Syst (2014) 38:10
Page 7 of 10, 10
Table 1 The comparison between the proposed registration algorithm and other algorithms on ten simulated image pairs Image Pairs
No. 1 No. 2 No. 3 No. 4 No. 5 No. 6 No. 7 No. 8 No. 9 No. 10 Average
Initial error before registration
After registration Single step algorithm [2]
Multiresolution algorithm [23]
Proposed algorithm
RMS (δ0)
RMS
R%
Elapsed time (s)
RMS
R%
Elapsed time (s)
RMS
R%
6.098 6.892 6.442 6.201 6.363 6.270 6.813 6.920 7.399 7.494 –
3.903 2.725 2.659 3.621 3.965 2.233 3.717 2.995 3.898 3.266 –
36.0 60.5 58.7 41.6 37.7 64.4 45.4 56.7 47.3 56.4 50.5
34.4 33.1 34.8 35.9 33.4 35.8 33.5 36.1 33.9 33.5 34.4
1.511 1.015 0.870 1.603 1.254 1.542 1.532 1.272 1.784 2.331 –
75.2 85.3 86.5 74.1 80.3 75.4 77.5 81.6 75.9 68.9 78.1
8954 9958 8745 9180 9553 9210 8768 9003 9548 9346 9226
1.124 0.781 0.875 0.882 1.154 1.316 1.191 1.020 1.337 1.037 –
81.6 88.7 86.4 85.8 81.9 79.0 82.5 85.2 81.9 86.2 83.9
% % % % % % % % % % %
respectively by Meijering et al. [2], and Yang et al. [23]. The first one [2] is a single step algorithm in which a set of control points are extracted from strong edges of the mask image. The displacement of each control point is then computed by means of the template matching using the energy of difference histogram as similarity measure. Also, a hill-climbing approach is used as an optimization method. At last, the image warping is performed using a Delaunay triangulation and piecewise linear interpolation. The second algorithm [23] is a multiresolution registration algorithm in which the mask image is decomposed to coarse and fine sub-image blocks Fig. 6 Example of the registration results on simulated images. a Original live image from a coronary angiography sequence. b Simulated live image. It is warped from the original live image with 30 random control point pairs using thin-plate spline. c Subtraction image corresponding to the difference between a and b. d Subtraction image after registration by single step algorithm [2]. e Subtraction image after registration by multiresolution algorithm [23]. f Subtraction image after registration by our proposed registration algorithm
% % % % % % % % % % %
Elapsed time (s) % % % % % % % % % % %
285 336 319 323 304 308 325 322 314 305 314
iteratively and each block is rigidly registered to the live image using Powell optimization method [28]. The calculated registration transform in each decomposition level is used as an initial guess for the next finer level. The comparison between our proposed algorithm and the single step [2] and multiresolution [23] algorithms is accomplished according to the registration performance and computational time. The performance is evaluated based on RMS pixel error between the true corresponding control points and the ones estimated by registration of the live and simulated live image pairs. Table 1 lists the achieved registration results
10, Page 8 of 10 Table 2 Evaluation of affine parameter validation procedure on the registration performance
J Med Syst (2014) 38:10
Test image pairs
No. 1 No. 2 No. 3 No. 4 No. 5 No. 6 No. 7 No. 8 No. 9 No. 10 Average
Registration results of proposed algorithm With affine parameter validation
Without affine parameter validation
RMS
R%
RMS
R%
1.124 0.781 0.875 0.882 1.154 1.316 1.191 1.020 1.337 1.037 –
81.6 88.7 86.4 85.8 81.9 79.0 82.5 85.2 81.9 86.2 83.9
1.208 2.191 1.364 1.956 1.669 1.729 1.743 2.114 1.682 1.450 –
80.19 68.21 78.83 68.46 73.77 72.42 74.41 69.45 77.26 80.65 74.37
by our proposed algorithm and two presented algorithms in [2] and [23] for the ten simulated image pairs. The results of our proposed algorithm in Table 1 are obtained at level 8 of the image decomposition. For the ten test image pairs, the reduction ratio of RMS error (R) of 83.9 % has been obtained by our proposed registration algorithm that is better than the corresponding value for multiresolution algorithm in [23], while our algorithm has considerably lower computational time. The computational time of the single step algorithm in [2] is lower than computational time of our proposed algorithm, however, the single step algorithm results in poor registration performance and our algorithm considerably outperforms for all image pairs. Figure 6 shows the
Fig. 7 Multiresolution registration results for three coronary angiography images. Columns from left to right correspond to mask images, live images, the original subtraction of the respective mask and the live images, and the subtraction after correction of motion artifacts using the proposed registration algorithm with eight decomposition levels, respectively
% % % % % % % % % % %
% % % % % % % % % % %
registration results of these algorithms for one of the simulated image pairs. As, it can be seen in Fig. 6d, e and f, the extensive artifacts still remain in the subtraction images when single step algorithm [2] or multiresolution algorithm [23] are applied. As previously mentioned in “Affine Parameter Validation” section, in our application, we employ an affine parameter validation procedure to reduce the undesirable effect of incorrect affine registration of blocks on the global elastic transformation. To evaluate this procedure on the registration performance, we also tested our algorithm with and without affine parameter validation on the simulated image pairs. Table 2 demonstrates the significant improvement in the registration results by using affine parameter validation.
J Med Syst (2014) 38:10
Registration Results on Clinical Images We carried out several experiments on clinical datasets to evaluate the performance of proposed multiresolution registration algorithm in the motion artifact reduction of DSA images. By several experiments on the coronary angiography images with different levels of decomposition, it was concluded that the best results are often achieved in the level 8 where the size of sub-image blocks is equal to 64×64 pixels. In the multiscale affine registration of sub-image blocks, a Gaussian pyramid is constructed for both the source and target images which the number of pyramid levels is selected such that the size of the coarsest level is greater than or equal to 64×64 pixels (because the coarser levels significantly lose the gradient patterns of original image which can degrade initial affine estimation). Figure 7 shows the multiresolution registration results for the three cardiac angiography images of datasets. In this figure, columns from left to right correspond to the mask images, the live images, original subtraction of the mask and live images, and the subtraction after correction of motion artifacts, using the proposed multiresolution registration algorithm, respectively. These results demonstrate that the motion artifacts in the original DSA images are largely removed after registration.
Conclusions In this paper, a new elastic registration algorithm for motion artifact reduction in DSA was presented. This method is based on a multiresolution search strategy in which both the mask and live images are decomposed to coarse and fine sub-image blocks iteratively. All the corresponding block pairs at each decomposition level are affinely registered using least-squares optimization which for fast convergence and avoiding the pitfalls of iterative optimization method, the entire procedure is done in a multiscale framework. To account for intensity variations, an explicit change of the contrast and brightness is also incorporated into our block registration model. Furthermore, in order to detect the probable incorrect registrations, a validation procedure is performed after affine registration of each corresponding block pairs based on the similarity assessment. The estimated affine parameters of registered sub-image blocks in each decomposition level are then smoothly interpolated by mean of the thin-plate spline interpolation function to generate a global transformation at each level. To evaluate our algorithm, we tested it with various clinical and simulated datasets. The obtained registration results demonstrate that our proposed algorithm outperforms other presented algorithms in the literature in terms of visual and quantitative evaluations. The experiments on the coronary angiography images also show the effectiveness of our method in the motion artifact reduction.
Page 9 of 10, 10 Acknowledgments The authors would like to thank Dr. Hashemi of Sina Heart Hospital in Isfahan, Iran, for providing us the datasets used in the experiments.
References 1. Brody, W. R., Digital subtraction angiography. IEEE Trans. Nucl. Sci. 29(3):1176–1180, 1982. 2. Meijering, E. H. W., Zuiderveld, K. J., and Viergever, M. A., Image registration for digital subtraction angiography. Int. J. Comput. Vis. 31(2/3):227–246, 1999. 3. Zitova, B., and Flusser, J., Image registration methods: A survey. Image Vis. Comput. 21(11):977–1000, 2003. 4. Brown, L. G., A survey of image registration techniques. ACM Comput. Surv. 24(4):325–276, 1992. 5. Elsen, V. D., Pol, P. A., and Viergever, E.-J. D., Medical image matching: A review with classification. IEEE Eng. Med. Biol. Magazine 12(1):26–39, 1993. 6. Maintz, J. B. A., and Viergever, M. A., A survey of medical image registration. Med. Image Anal. 2(1):1–36, 1998. 7. Hajnal, J. V., Hill, D. L. G., and Hawkes, D. J., Medical image registration. CRS Press, Baton Rouge, Florida, 2001. 8. Cox, G. S., De Jager, G. “Automatic registration of temporal image pairs for digital subtraction angiography,” In: Loew M. H. (ed.), Proceedings of the SPIE, Medical imaging: image processing, vol. 2167, pp. 188–199, Bellingham, WA, 1994. 9. Nejati, M., Sadri, S., and Amirfattahi, R., “Nonrigid Image Registration in Digital Subtraction Angiography Using Multilevel B-Spline”, BioMed Research International, pp. 1–12, vol. 2013. 10. Buzug, T. M., Weese, J., Fassnacht, C., and Lorenz, C., Using an entropy similarity measure to enhance the quality of DSA images with an algorithm based on template matching. In: Höhne, K.-H., and Kikinis, R. (Eds.), VBC’96, Lecture Notes in Computer Science, vol. 1131. Springer, Berlin, pp. 235–240, 1996. 11. Buzug, T. M., Weese, J., Lorenz, C., and Beil, W., Histogram-based image registration for digital subtraction angiography. In: Del Bimbo, A. (Ed.), ICIAP’97, Lecture Notes in Computer Science, vol. 1311. Springer, Berlin, pp. 380–387, 1997. 12. Buzug, T. M., and Weese, J., Image registration for DSA quality enhancement. Comput. Med. Imaging Graph. 22(2):103–113, 1998. 13. Cao, Z., Liu, X., Peng, B., Moon, Y.-S. “DSA image registration based on multiscale Gabor filters and mutual information”. In: Proceedings of the IEEE ICIA’05, Hong Kong and Macau, China, pp. 105–110, 2005. 14. Taleb, N., and Jetto, L., Image registration for applications in digital subtraction angiography. Control. Eng. Pract. 6(2):227–238, 1998. 15. Bentoutou, Y., Taleb, N., Chikr El Mezouar, M., Taleb, M., and Jetto, L., An invariant approach for image registration in digital subtraction angiography. Pattern Recogn. 35(12):2853–2865, 2002. 16. Bentoutou, Y., and Taleb, N., Automatic extraction of control points for digital subtraction angiography image enhancement. IEEE Trans. Nucl. Sci. 52(1):238–246, 2005. 17. Bentoutou, Y., and Taleb, N., A 3-D space-time motion detection for an invariant image registration approach in digital subtraction angiography. Comput. Vis. Image Underst. 97(1):30–50, 2005. 18. Meijering, E. H. W., Niessen, W. J., and Viergever, M. A., Retrospective motion correction in digital subtraction angiography: a review. IEEE Trans. Med. Imaging 18(1):2–21, 1999. 19. Ping, L., Hong, N., Ye, S. “An efficient method for image registration in DSA,” In: Proceedings of the ISISE’08, Shanghai, China, vol. 2, pp. 551–554, 2008. 20. Nejati, M., Amirfattahi, R., Sadri, S.“A fast image registration algorithm for digital subtraction angiography”, In: Proceedings of 17th
10, Page 10 of 10
21.
22.
23.
24.
Iranian Conference of Biomedical Engineering (ICBME2010), Isfahan, Iran, pp. 1–4, November 2010. Chu, Y., Bai, N., Ji, Z., Chen, S., Mou, X. “Registration for DSA image using triangle grid and spatial transformation based on stretching”, In: Proceedings of the 8th International Conference on Signal Processing, Beijing, China, 16–20 Nov. 2006. Wang, J., Zhang, J. Q. “An iterative refinement DSA image registration algorithm using structural image quality measure”, In: Proceedings of the 5th IIH-MSP, Kyoto, Japan, pp. 973–976, 2009. Yang, J., Wang, Y., Tang, S., Zhou, S., Liu, Y., and Chen, W., Multiresolution elastic registration of X-ray angiography images using thin-plate spline. IEEE Trans. Nucl. Sci. 54(1):152–166, 2007. Nejati, M., Pourghassem, H. “Multiresolution search strategy for elastic registration of x-ray angiography images”, In Proceedings of International Conference on Intelligent Computation and Bio-
J Med Syst (2014) 38:10
25.
26.
27.
28.
Medical Instrumentation (ICBMI 2011), pp. 216–219, Wuhan, China, 14–17 Dec. 2011. Periaswamy, S., and Farid, H., Elastic registration in the Presence of Intensity Variations. IEEE Trans. Med. Imaging 22(7):865–874, 2003. Bookstein, F. L., Principal warps: thin-plate splines and the decomposition of deformations. IEEE Transact. Pattern Anal. Mach. Intell. 11(6):567–585, 1989. Farid, H., Simoncelli, E. P.“Optimally rotation-equivariant directional derivative kernels”, In Proceedings of International Conference on Computer Analysis of Images and Patterns, pp. 207–214, Berlin, Germany, Sept. 1997. Powell, M. J. D., An efficient method for finding the minimum of a function of several variables without calculating derivatives. Comput. J. 7(2):155–162, 1964.
