*This work was supported by the Project of the National Fundamental Research 973 Program of China under Grant

No. 2003CB716105.
DSA Image Registration based on Multiscale Gabor
Filters and Mutual Information
Zhiguo Cao , Xiaoxiao Liu and Bo Peng Yiu-Sang Moon
Key Laboratory of Ministry of Education for Image Processing
and Intelligent Control
Department of Computer Science and Engineering
Institute for Pattern Recognition and Artificial Intelligence Chinese University of Hong Kong
Huazhong University of Science and Technology Shatin N.T.
Wuhan 430074, China Hong Kong, China
{zgcao & sharonLiu}@hust.edu.cn [email protected]
Abstract - In clinical practice, digital subtraction
angiography (DSA) is a powerful technique for the visualization
of blood vessels in X-ray image sequences. Different with
traditional DSA image registration processes, in our proposed
image registration method, the control points are selected from
the vessel centerlines using multiscale Gabor filters, and mutual
information (MI) is then taken as the similarity criterion to find
the correspondences. Experimental results demonstrate our
algorithm efficiently yields satisfying registration result for DSA
images.
Index Terms - Image registration; Digital Subtraction
Angiography; Multiscale Gabor Filters; Mutual information;
I. INTRODUCTION
In modern medical clinical practice, Digital Subtraction
Angiography (DSA) is a widely used imaging tool for
visualizing the blood vessels. Especially, the coronary
angiography is the gold-standard image technique for the
visualization of the motion and the morphology of the arteries.
During the DSA imaging process, a successive series of
images is taken to show the passage of a bolus of injected
contrast medium through the interested vessels. So-called
mask images are taken prior to the arrival of the contrast
medium, and the live images are taken when the contrast
medium is full of the vessels. The background is largely
removed by subtracting a mask image from the live image.
Most ideally, only the high-contrast vessels are left in the
picture.
However, due to the unavoidable and complex motion of
the tissues inside the human body, the image sequences are
often full of artefacts, including peristaltic motion of the
patient as well as the superimposed spine, catheter and bowels.
It is not likely to directly conduct subtraction without
alignment. Instead of developing new imaging devices, the
easier and cheaper way is to take image registration techniques
by calculating the correspondence between pixels in the
successive images.
Lots of image registration strategies have been carried
out during the past decade [1]. A survey specific in medical
field can be found in [2]. In [3], the authors provide a detailed
and extensive framework of DSA image registration.
Generally speaking, the registration involves two steps: (a)
find out the correspondence between the pixels in the live and
mask images, and (b) warp one of the images according to the
correspondence with respect to the other. To reduce the
computation price, only the control points are taken into
account for the first step. Later, Bentoutou etc. [4] pointed out
that the control points should be selected from the edges of the
live images, instead of mask images [3], to ensure the vascular
structures are included in order to better register the dissimilar
regions. They take the combined invariant moments as the
similarity measure to compute the displacements of each
control point, and the registration results are very well.
However, there is still some room for further improving the
registration process. First, the common edge detection filters
may not behave very well when encountering the low contrast
vessel images with noisy background. Second, the robustness
of the similarity criterion (invariant moments) towards the
dissimilarity of vessel information in the live image is not
fully investigated. Also, the computation cost of calculating
the 5 moments at each position in the searching window to get
the displacement of each control point still may not be cheap
enough for clinical practice.
The purpose of this paper is to develop an effective and
fast method of DSA image registration to improve the
visualization of blood vessels. The rest of the paper is
organized as follows. In section , control points selection is
developed using multiscale Gabor filtering and the
displacement of the each control point is calculated by MI
criterion, illustrated in section . The detailed approach of a
smoothed thin-plate spline warping strategy can be found in
section . And some experiments results are shown in section
. Finally, a conclusion is given in section.
II. CONTROL POINTS SELECTION
To explicitly compute the correspondence for each pixel
might be desirable but unpractical, due to the expensive
computation cost. The possible alternative is to merely
calculate the correspondence for a selected set of so-called
control points. And the overall alignment can then be obtained
by interpolation. Usually, the control points can be selected
manually, or by taking a regular grid, or by auto extracting
image features such as strong edges [3]. Considering the
particular case of DSA images, where artefacts are observable
only in the areas with dynamic objects, we extract the control
points from centerlines of these moving objects, like vessels
and catheters, using multiscale Gabor filtering [5].
0-7803-9303-1/05/$20.00 2005 IEEE
105
Proceedings of the 2005 IEEE
International Conference on Information Acquisition
June 27 - July 3, 2005, Hong Kong and Macau, China
Compared with the edge points, the skeleton points, as
the important feature of the tubular objects, can better control
the deformation of the whole image, and the reduced number
of the control points will help to curtail the computation cost.
In addition, common edge detectors often fail to detect the
vessel area with low contrast, while Gabor filters with
different directions and frequency bands and can well enhance
the centreline signals of tubular objects and conduct noise
suppression
A. Centerline extraction by Multiscale Gabor filters
Here, we exploit the 2-D even Gabor filter (1) for
extracting control points from the live image.
( ) ( ) ( )
' ' '
2 cos , , Fx y x g y x h =
(1)
where ( ) ( ) cos sin , sin cos ,
' '
y x y x y x + + =
,
denoting the orientation of the filter. Any desired
orientation can be achieved via a rigid rotation in the x-y
plane. F is the spatial central frequency, ) , ( y x g is the 2-D
Gaussian function as in (2)
( )

(
(

|
|
.
|

\
|
+
|
|
.
|

\
|
=
2
2
2
1
exp
2
1
,
y x y x
y x
y x g

. (2)
Where
x

and
y

correspond the horizontal and vertical

and
y

are
related to the filters frequency bandwidth and orientation
bandwidth in the way illustrated in (3) and (4):
1 2
1 2 1
2
2 ln

+
=
F
F
B
B
x
F

(3)
( ) 2 tan
1 1
2
2 ln

B F
y
=
(4)
F
B
is the spatial frequency bandwidth and

B
is the
orientation bandwidth.
Fig.1 Two dimensional even Gabor filter in spatial domain
The following steps are taken to extract the centerlines.
An example is given in Fig 2.to help illustration:
1. Conduct convolution of the original DSA image with
a group of pre-designed Gabor filter templates. According to
the width of the tubular objects in the live image, we modulate
the filters into specific 5 frequency band (the 5 scales
corresponds to the 5 degrees of width), and 12 directions in
each scale. Extracting the maximum response value among
those convolution results in different frequencies and
orientations, we get the enhanced picture of the vessel
centrelines (Fig 2. (b)).
2. Extract the ridge points or local maximums from the
maximum response map, using non-maximum suppression
method, guided by the maximum response direction.
3. Threshold at the 90% of the histogram is able to get
rid of most of the background ridges, which are not supposed
to belong to the strong moving objects.
4. Utilizing the continuousness of the scales and
directions along the centerline curve, the smooth centerlines of
interest are obtained by further decline short and broken lines.
(Fig. 2(d))
(a) Original DSA image (b) Maximum response Map
(c) Ridge detection result (d) Extracted centreline (white )
Fig.2 Centerline extraction by Mutiscale Gabor filters
B. Control points sampling
A slipping window of m m is taken to make a
sampling of the centerline points so as to ensure one point in
each window at most. The window size m is related to the
size of the image itself. As for a specific image, the larger the
size m, the smaller the number of the control points. However,
to choose the best number of the control points is to make a
best balance between the computation cost and the registration
effect, which has no rigorous standards. Here we set the ratio
of control points out of the sum of the pixels to be 0.1% ~
0.5%, and our experimental results demonstrate that the small
group of accurately selected centerline points is able to well
control the deformation during the TPS interpolation.
CONTROL POINTS DISPLACEMENT CALCULATION
After selecting the control points
{ ( , )}
i
V v x y =
,
1,... i N =
from the live image, we now calculate the
displacement of each point, viz., to search for the
correspondences. There are two major categories of
displacement calculation methods: gradient-based optic-flow
techniques and the template-matching based techniques. We
choose the latter, which can be made more robust against the
inflow of contrast than optic-flow techniques, by a proper
choice for the similarity criterion [3]. Then the similarity
criterion comes to be the crucial factor of the whole
registration process.
A number of similarity measurements have been
proposed to determine the correspondence between regions
106
[3]; however, most of them are not robust towards the regions
with intensity changes, especially where obvious
dissimilarities appear, thus not suited to DSA applications.
The major dissimilarity between the live image and the mask
image is the blood vessel information: the mask image has no
vessel information at all, while the live images contain vessels
in different sizes. Therefore, we need a similarity criterion,
which is not sensitive to the vessel signal in order to find the
correspondence. Besides, the ideal similarity criterion should
also be invariant to the local transformations in certain degree,
which can be simplified here to be rotation, translation and
blur. The estimated displacements of the correspondence in
DSA images are in small scales [4]. Reasonably, we assume
that the largest deformation parameters would not exceed 5
degrees angle for rotation, 10 pixels shift for translation and
Gaussian blur with the
1.5 sigma
. The above expectations
form the standard for choosing the ideal similarity
measurement in our DSA application.
In this paper, the mutual information (MI) or relative
entropy is used to describe the similarity of the DSA image
pairs. MI is a very powerful tool to measure the statistical
dependence between two random variables or the amount of
information that one variable contains about the other. It has
been applied on a large variety of applications for its
robustness to measure image intensities, especially in
multimodality image registration [6]. We, therefore, not only
utilize MI to measure the final effect of our registration
algorithms, but also take it as the similarity criterion of the
correspondence search. And we further discover the
robustness characteristic of MI by testing it with technically
designed synthetic experiments for DSA images.
A. Mutual Information theory
Generally, the entropy of random variable A is defined
as
( ) ( ) log ( )
A A
a
H A P a P a =
_
. (5)
Where a A and
( )
A
P a
is the possibility distribution of A .
Similarly, the joint entropy of two random variables A and B
is defined as
,
( , ) ( , ) log ( , )
AB AB
a b
H A B P a b P a b =
_
. (6)
Where a A b B and
( . )
ab
P a b
is the joint possibility
distribution of two variables. If
( | ) H A B
stands for the
condition entropy of A when B is known, then the difference
of
( ) H A
and
( | ) H A B
expresses the entropy of A contained
in B, i.e. MI. So the MI between two systems can be presented
by (7):
( , ) ( ) ( ) ( , )
( ) ( | )
( ) ( | )
I A B H A H B H A B
H A H A B
H B H B A
= +
=
=
. (7)
Obviously, the MI, which measures how much one image
represents the information in another image, will reach
maximum if the spatial position of two images is the same.
Therefore, MI may be used as similarity measure between two
images and is usually estimated by the following expression:
( ) ( )
( )
( ) ( ) ,
,
, , log
AB
AB
a b A B
P a b
I A B P a b
P a P b
=
_

(8)
B. Synthetic experiments for MI measure
In order to test the robustness of MI measure in the
dissimilar area, where vessel information is visible, we
synthesize a precisely registered sub-image pairs from original
live image for tests (Fig.3): First clip a small 5 51
rectangular region, filled with contrast information, from the
original 165 165 live image, and then paste it onto the center
of the 51 51 mask template, which is arbitrarily selected from
the background. So we get the 51 51 template pair, which can
be regarded as strictly registered with each other.

165 165 51 51 5 51 51 51
Fig.3 The formation of the synthesized image pairs
The following 4 small experiments are carried out to test
the robustness of the MI similarity criterion towards
dissimilarities, using the above synthetic image pairs.
Test 1: Taking a 21 21 sub-template from the live
template, calculate the MI value at each position within
1111 neighborhood (translation from -5 to +5 pixels at both
x and y directions) at the corresponding mask template (Fig.
4). We can see the peak appears in the template center, which
tells us that the grey information in the 2121 window is rich
enough to determine the similarity that the vessel information
in the live image has little influence. Certainly, if we take a
smaller window, where the vessel information is too much to
be ignored, the peak will be weaker and even destroyed. Since
we hope the window size to be as small as possible in order to
reduce the computing cost, it is necessary to find out how the
MI measure is related to the window size, which leads us to
the Test 2.
Fig.4 Peak appears at the matching position in MI measure space
Test 2: Change the sub-template size n to see how the
MI measure works towards the vessel information. We find
107
that the peaks are becoming less strong as the decrease of
template size n . And when the area of the vessel part (5 pixels
wide) takes up almost half of the template (
9,11 n =
), the
difference of MI value are too small to accurately detect the
maximum. Considering the deformations in real DSA images,
we take the templates with the size 21 21 (big enough to
counteract the deformations, and computationally acceptable)
in our later algorithms to deal with the vessels with width less
than 5 pixels. Experimentally speaking, the template size n is
supposed to be as large as about 4 times of the estimated
biggest vessel width.
Another two tests are used to discover the invariant
ability of MI method towards transformation (rotation and blur
separately) in small scale in DSA data.
The tests demonstrate that with a proper template size,
the maximum of MI is always at the right correspondent
position, invariant to rotation and blur in a certain range,
which is robust enough to deal with the local transformation in
DSA image sequences.
By searching the maximal value in the 10 10
neighborhood of control point set
{ ( , )}
i
V v x y =
of the MI
measure space, we can obtain the best matching corresponding
point set { ( ', ')}
i
U u x y = ,
1,... i N =
.
C. Sub-pixel accuracy
Even sub-pixel misalignments may produce significant
artefacts in the subtraction images. Instead of translating the
mask image by sub-pixel increments around the found best
match position, which is computationally expensive, here, we
utilize the mutual information value to do an analytical match-
interpolation [7]. The need for greater precision requires that
the coordinates for the corresponding peak in the match space
of mutual information.
At each matched position
( ', ') x y
, given the MI values of
its discrete neighbourhood, the estimated continuous
coordinates ( , ) x y
of the match peak are calculated as:
( , ) ( , )
( , ) ( , ) ( , )
( , ) ( , ) ( , )
( ', ') ( ' 1, ')
1
'
2 2 ( ', ') ( ' 1, ) ( ' 1, ')
1 ( ', ') ( ', ' 1)
' .
2 2 ( ', ') ( ', 1) ( ', ' 1)
x y x y
x y x y x y
x y x y x y
I x y I x y
x x
I x y I x y I x y
I x y I x y
y y
I x y I x y I x y

= +
+

= +
+
(9)
Where
( , )
( ', ')
x y
I x y
is the MI value between the live
template centred at
( , ) x y
and the mask template centred
at
( ', ') x y
.
. IMAGE WARPING BASED ON SMOOTHING TPS
The Thin-plate Spline (TPS) functions provide us an
effective way to generate a smoothly interpolated spatial
mapping with adherence to two sets of landmark points [8].
As a natural non-rigid extension of the affine map, the TPS
can be cleanly decomposed into affine and non-affine
subspaces while minimizing a bending energy based on the
second derivative of the spatial mapping [9]. In this paper, the
smoothing TPS is utilized to decrease the influence of wrong
correspondences.
Given the control point pairs { ( , )}
i
U u x y =
and
{ ( , )}
i
V v x y =
,
1,... i N =
, the TPS energy function is
expressed as
2 2 2
2
2 2 2
( ) 2 2
1
( ) ( ) 2( ) ( )
N
TPS f i i
i
f f f
E u f v dxdy
x x y y

=
(
= + + +
(


_
))
(10)
Where f is the mapping function between U and V , and is
a regularization parameter, which adjusts the degree of
smoothness of the TPS function. Large values of greatly
limit the range of non-rigidity of the transformation. On the
other hand, the transformation can turn out to be very flexible
at small values of . For each fixed , there exists a unique
minimizer f which is denoted by:

( ) ( ) f v v A v W = +
(11)
Where A is a 3 3 matrix, representing the affine
transformation and W is a 3 N warping coefficient matrix,
representing the non-affine deformation. The vector ( ) v is a
1 K vector for each point
i
v , where each
entry
2
( ) log
j i j i j i
v c v v v v =
.
When we substitute f in () using (), the TPS
energy function becomes:
2
( , ) ( )
T
TPS
E A W U VA w trace W W = +
(12)
Where is a N N matrix made up by
( )
i
v
.
Then,
QR
decomposition is used to separate the affine
and non-affine warping space:
1 2
[ ]
0
R
V QQ
| |
=
|

\ .
(13)
Where
1
Q
and
2
Q
are orthonormal matrices, R is upper
triangular. The final solution for W and A are:
1
2 2 2 ( 3) 2

( )
T T
N
W Q Q Q I Q U

= +
(14)

1
1

( );
T
A R Q V W

=
(15)
Then, using function with the calculated parameters, we
apply this warping to each pixel of the mask image
( , ) M x y
to
obtain a new registered mask image
*
( , ) M x y
. During the
warping, however, the calculated new coordinates might not
be integers. If we simply cut out the decimal parts to get
integers, put the inaccuracy resulted by round operation aside,
many positions can not get gray value or get values time after
time, especially in areas where larger deformation happens. To
improve the accuracy of TPS warping result, we define the
following interpolation function (16) to calculate
*
( , ) M x y
:
108
2
1 *
2
1
2 2
( , )
( , ) ,
1
( ) ( ) ,
1
1,
K
i i i
i i
K
i i
i i i i
i i
i i
M x y
D
M x y
D
x x y y x x and y y
D D
x x and y y
=
=
=

+
= ,

= =

_
_
(16)
For each position ( , ) x y , , x y N there are some calculated
coordinates ( , )
i i
x y ,
, , 1...
i i
x y R i K =
, which are all close
to ( , ) x y that the distance
i
D
between them are less than 1
pixel. We distribute the K gray values by giving them each a
weight according to their distance to the integer coordinates.
Thus we get smoothing and accurately registered new mask
image.
. EXPERIMENTAL RESULTS
We carried out experiments on various sets of DSA
images, most of the artefacts can be removed and the weak
vessels can be better visualized after registration.
Two examples are given in this section to illustrate our
registration techniques. Fig.5 illustrates the registration
process of the cerebral image pairs in detail, with the middle
products shown orderly. Fig.5(c) is the maximum response
map after multiscale Gabor filtering, vessel structure gets well
enhanced. The extracted centerlines are superposed on to the
original live image in Fig.5(c), the extracted strong tubular
structures are supposed to be the key moving objects that need
to be registered. Calculated the displacement vectors of each
control points are depicted on Fig.5 (f) by the arrows. Then
TPS function warping deforms the original mask image both
globally and locally and produces the new mask image in
Fig.5(g).The increased MI value between mask image
(
*
( , ) ( , ) M x y and M x y
separately) and the original live image
( , ) L x y , from 2.90 to 4.70, demonstrates the effectiveness of
our algorithms. The final better vision of the vessel structures
(Fig.5 (i)) is obtained by simply subtracting the new mask
image from the original live image. Compared with the direct
subtraction result before registration, the improvement is
obvious. Fig.6 shows the registration result of peripheral
image pairs and the MI value increases from 1.60 to 2.00.
(a) Original live image ( , ) L x y
(b) Original mask image
( , ) M x y
(c) Maximum response map (d) Centerline points
(f) Displacements calculation
(g) Registered mask image
*
( , ) M x y
(h) Direct subtraction result (i) The result of our approach
Fig.5 Registration of cerebral images
109
(a) Original live image (b) Original mask image
(c) Control points on the centerlines
(d) Displacement calculation
(e) Direct subtraction result (f) The result of our approach
Fig.6 Registration of peripheral images
CONCLUSIONS
The techniques we propose for the automatic and fast
DSA image registration successfully remove most of the
motion artefacts thus providing a better difference image for
blood vessel visualization. The centerline points selected by
multiscale Gabor filtering controls the image transformation
very well, and the similarity criterion of Mutual Information is
robust enough to detect the correspondences regardless of the
dissimilarity of vessel information in the live image. Tied with
the advantage of thin-plate spine in non-rigid deformation, the
proposed algorithm yields good global and local registration.
ACKNOWLEDGEMENT
Many thanks to Dr. N.Table of Dpartement d
Electronique at universit de Sidi Bel Abbes for providing us
the DSA images used in our paper. Those image sequences are
originally from Boston City Hospital in Boston, USA and
Saint Quentin Hospital, in Saint Quentin, France.
REFERENCES
[1] L.G. Brown, A survey of image registration techniques, ACM
Computing Surveys, 1992, vol.24, no.4, pp.325-376
[2] J.B.A. Maintz and M.A. Viergever, A survey of medical image
registration, Medical Image Analysis, 1998, vol.2, no.1, pp.1-36
[3] Erik H.W. Meijering, Karel J Zuiderveld and Max A.Viergever, Image
registration for digital subtraction angiography, International Journal of
Computer Vision, 1999, vol.31, no.2/3, pp.227-246
[4] Y. Bentoutou, N. Taleb, M. Chikr El Mezouar, M.taleb and L. Jetto, An
invariant approach for image registration in digital subtraction
angiography, Pattern Recognition, 2002, vol.35, pp. 2853-2865
[5] Nong Sang, Qiling Tang, Xiaoxiao Liu, Wenjie Weng, Multiscale
centerlie extraction of angiogram vessels using gabor filters, CIS2004,
LNCS 3314, pp.570-575
[6] Frederik Maes, Andr Collignon, Dirk Vandermeulen, Guy Marchal and
Paul Suetens, Multimodality Image Registration by Maximization of
Mutual Information, IEEE Transactions on Medical Imaging, 1997,
vol.16, pp.187-198
[7] G.S.Cox and G.de Jager, Automatic registration of temporal image
pairs fro digital subtraction angiography, in: M.H. Loew
(Ed.),Proc.SPIE, Medical Imaging:Image Processing, 1994, vol.2167,
Bellignham, Washington, USA, pp.188-199
[8] Fred L. Bookstein, Principal Warps: Thin plate splines and the
decomposition of deformations, IEEE Transactions on Pattern Analysis
and Machine Intelligence, 1989, vol.11, pp.567-585
[9] Haili Chui, Anand Rangarajan, A new point matching algorithm for
non-rigid registration, Computer Vision and Image Understanding,
2003, vol.89, pp.114-141
110