IEEE TMI SPECIAL ISSUE MB0016 — [RUNNING ENHANCED CLASS V1.5]
1
Predictive (un)distortion model and 3D Reconstruction by Biplane Snakes *Cristina Cañero, Fernando Vilariño, Josepa Mauri and Petia Radeva
Abstract— This paper is concerned with the 3D reconstruction of coronary vessel centerlines and with how distortion of X-ray angiographic images affects it. Angiographies suffer from pincushion and other geometrical distortions, caused by the peripheral concavity of the Image Intensifier (II) and the non-linearity of electronic acquisition devices. In routine clinical practice, where a Field-Of-View (FOV) of cm is commonly used for the acquisition of coronary vessels, this distortion introduces a positional error of upto 7 pixels for an image matrix size of and a FOV of cm. This error increases with the size of the FOV. Geometrical distortions have a significant effect on the validity of the 3D reconstruction of vessels from these images. We show how this effect can be reduced by integrating a predictive model of (un)distortion into the biplane snakes formulation for 3D reconstruction. First, we prove that the distortion can be accurately modeled using a polynomial for each view. Also, we show that the estimated polynomial is independent of focal length, but not of changes in anatomical angles, as the II is influenced by the earth’s magnetic field. Thus, we decompose the polynomial into two components, namely the steady and the orientation-dependent component. We determine the optimal polynomial degree for each component, which is empirically determined to be 5 for the steady component, and 3 for the orientation-dependent one. This fact simplifies the prediction of the orientation-dependent polynomial, since the number of polynomial coefficients to be predicted is lower. The integration of this model into the biplane snakes formulation enables us to avoid image unwarping, which deteriorates image quality and therefore complicates vessel centerline feature extraction. Moreover, we improve the biplane snake behavior when dealing with wavy vessels, by means of using Generalized Gradient Vector Flow (GGVF). Our experiments show that the proposed methods in this paper decrease up to the reconstruction error obtained when geometrical distortion effects are ignored. Tests on imaged phantoms and real cardiac images are presented as well. Index Terms—3D reconstruction, angiography, geometrical distortion, snakes.
I. I NTRODUCTION IGITAL X-ray angiography provides high quality images of coronary vessels. It represents the most common image modality applied in clinical practice to assist coronary diagnosis and therapy. X-ray medical images are mainly used to
D
Manuscript received December 1, 2001. This work was supported in part by the catalonian Departament d’Universitats, Recerca i Societat de la Informació under the grant 1999FI 00753 APTIND CVC and by the spanish Ministerio de Ciencia y Tecnología under the projects TIC2000-1635-C04-04 and TIC20000399-C02-01. *C. Cañero, F. Vilariño and P. Radeva are with the Computer Vision Center, Edifici O, campus UAB, 08193 Bellaterra, Barcelona, Spain, and with the Department of Computer Science, Universitat Autònoma de Barcelona, Spain (e-mail:
[email protected]). J. Mauri is with the University Hospital "Germans Tries i Pujol", Badalona, Barcelona, Spain.
provide qualitative information (i.e. presence of vessel stenosis). However, medical images as a source of quantitative information about vessel lesions are not sufficiently explored. The percentage of stenosis assessment and lesion length measurement have major potential as Quantitative Coronary Angiography (QCA) applications. These measurements are obtained in clinical practice by comparing the size of the lesion to the width of the catheter appearance in the image, disregarding the imprecision due to the perspective projection. In order to obtain reliable vessel measurements, 3D reconstruction of vessels is necessary. Other applications that need 3D reconstruction of vessels from angiography are related to the determination of optimal views of the coronary tree to improve the visualization and the fusion of angiography with other medical image modalities such as intravascular ultrasound images to obtain a real volumetric 3D reconstruction of the vessel [1], [2]. In order to obtain an accurate three-dimensional reconstruction of vessels from angiography, X-ray image distortion must be taken into account [3]. The peripheral concavity of the Image Intensifier and non-linearity of the electronic devices cause pincushion and other geometrical distortions. Moreover, the influence of the earth’s magnetic field is quite perceptible: the distortion depends on the orientation of the detector [4]. Some investigators use a few small lead markers attached to the II screen to estimate the distortion [5], [6], [7]. Since the geometrical distortion is complex, in most cases a large number of lead markers [4] or a calibration grid attached to the II screen, is necessary [8], [9], [10], [11], [4], [12], [13], [14], [3]. Applying this procedure for each image sequence acquisition may be inconvenient in clinical practice; hence, some investigators have suggested to calibrate the system only once [13]. Thus, they correct only the radial or pincushion distortion. However, the influence of the II orientation on the geometrical distortion [11] and its effects on 3D reconstruction are too large to be neglected [4], [3]. Other authors propose to predict the distortion [11], [14] from the anatomical angles (rotation and angulation) of the II by interpolation, obtaining a mean error of about pixels and a maximum of about three pixels, for an image matrix size of pixels [14]. Therefore, we conclude that although much work has been done, the problem of correcting X-ray image distortion remains open. In this work, we use a predictive model of geometrical distortion based on anatomical angles, since the geometrical distortion varies with the extrinsic parameters. The main differences with the work presented in [14] are that we predict a global model, whereas Kerrien et al. predict positions of the grid nodes in the image, and also that we use an optimal degree polynomial for the prediction.
IEEE TMI SPECIAL ISSUE MB0016 — [RUNNING ENHANCED CLASS V1.5]
Once the X-ray system is calibrated, we proceed with the 3D reconstruction of the coronary vessels. The methods proposed in [15], [16], [17], [18], [19], [20], [13], [3] use corresponding points to obtain the 3D reconstruction. However, automatic correspondence determination is difficult and must deal with the image ambiguities due to the vessel projections. On the other hand, manual assignment of corresponding points becomes time-consuming and tedious if high accuracy is needed. We showed in [21] that both problems could be reduced by using biplane snakes, which are derived from the physics-based model of snakes [22]. Such a curve deforms in order to adapt its projections to the images, inherently solving the point correspondence problem, and obtaining higher precision compared to the approach based on user-defined corresponding points in both X-ray images [21]. The biplane snake deforms in space to capture the spatial position of the vessel part by minimizing the distance between its projections and the vessel appearance in the X-ray images. However, similar to most snake models, the biplane snake has difficulties converging onto concave shapes when reconstructing wavy vessels. The problem of active contour convergence to concave objects is addressed in [23], [24] by defining a new kind of external force, called the Gradient Vector Flow (GVF) and an improved version, Generalized Gradient Vector Flow (GGVF). In this article, we propose a further development of biplane snakes to deform on a GGVF to obtain a robust 3D reconstruction of the centerline of coronary vessels. The remainder of this paper is organized as follows: in section II, we propose a predictive model for geometrical distortion based on the parameters of the X-ray acquisition system. Section III describes the biplane snakes and the GGVF method for 3D reconstruction. Section IV presents experimental results of both calibration and 3D reconstruction on synthetic images and real cardiac images. The paper finishes with conclusions and future work. II. (U N )D ISTORTION METHOD There is a strong relationship between the camera model used in Computer Vision applications and the angiography acquisition system. Let us consider a pinhole camera model and note the differences with the angiography model. Figure 1(a) illustrates the angiography acquisition model and figure 1(b) shows the pinhole camera model (see [25], [26] for more details). Plane , called the focal plane, is placed in front of the image plane . Distance is the focal length. Let point be the optical center, then the optical axis is the line passing through perpendicular to . The optical axis intersects at point , called the principal point. It should be noted that, from a geometric viewpoint, in both models we deal with the perspective projection; the main difference between them is that the object in the X-ray image is magnified and not inverted compared to the pinhole camera image. The transformation from 3D world coordinates to camera pixel coordinates is modeled in [27] as a process of four steps. By changing some details, we obtain an equivalent procedure. First, we perform a linear transformation from 3D world coordinate system
to 3D camera coordinate system
:
2
Fig. 1. The angiography acquisition model (a) and the pinhole camera model (b)
#
% #
"!
$
%
(1)
where matrix is a & '& rotation matrix and is a translation vector. constitute the extrinsic parameters. Second, a perspective projection is applied, following the pinhole camera model, which is also applicable in the angiographic context: )(
-,.
10 2 , /
*+ 15 2
(2)
8 +77 6
when
43 ,
where matrix embodies the intrinsic parameters. A standard expression for the matrix is: 8 8 ;: 8
,9.
/< 8 < :=3 ;
where : is the inverse pixel size and />< ?, 3 < are the coordinates of the principal point on the image matrix. This step differs from Tsai’s approach [27], where matrix dependes only on . There, the pixel aspect ration is addressed at a fourth step by including the aspect factor @BA , and expressing in pixels/cm. We do not introduce any aspect factor, as we consider it8 as8 a part of image distortion. Another difference with the Tsai’s approach is that in [27], the image center is supposed to be at , which is defined at the center of the image matrix. Instead, we include the parameters />< ?3 < . Equations (1) and (2) can be combined as follows: )(
C-,) D
)!
*+
#
% #
$
%
/
02
43
(3)
52
8 +77 6
when
However, ignoring distortion is unacceptable when doing 3D measurement [27]. This is also true in the context of angiographic projection, since the pinhole camera model is not accurate. Therefore, a third step to model geometrical distortion
IEEE TMI SPECIAL ISSUE MB0016 — [RUNNING ENHANCED CLASS V1.5]
3
is needed. Thus, we perform a non-linear transformation using / functions and to obtain distorted pixel coordinates 3 : /
3
3
/
In order to obtain a unwarped image, this is the model to be applied, since for each coordinate of the target undistorted image, we obtain a unique corresponding distorted coordinate on the image to be unwarped. By using a undistortion model to unwarp images, we ought to deal with image holes or multiple pixel values. Nevertheless, image unwarping introduces noise in the image, and is computationally expensive. So, it is preferable to avoid it. Instead, we propose to estimate also a model to obtain the undistorted coordinates: / 3
3
/
Thus, to project a 3D point on the distorted image (as done with biplane snakes, as we will see later) we use the distortion model, and to reconstruct a 3D point from a pair of points on the images, we use the undistortion model. Thus, no image unwarping is needed. In Computer Vision, for real cameras, the geometrical distor tion is introduced by the imperfect lens shape/ [26], and and can be approximated as polynomials of and 3 (see [27], [25], [28], [26]). However, in the case of digital angiography geometrical distortion is generally described as the combination of two effects from different sources: First, the peripheral concavity of the Image Intensifier and the non-linearity of the electronic devices that causes pincushion and other geometrical distortions. Second, the influence of the earth’s magnetic field is quite perceptible and the distortion depends on the orientation of the detector; hence, and change depending of the orientation of the II [11], [4]. Other non-fixed distortion sources may include fields emitted by nearby monitors, or even non-linear distortions with a warming-up of the image intensifier. Since in our case the monitors remain at a fixed position, we assume that the distortion introduced by them can also be assumed as depending of the orientation of the detector. Regarding the warming-up of the II, authors in [4] claim that the distortion remains constant even during several months, always after a warming-up of a minimum of 2h from the switching on of the II. Our aim is to propose a model to predict , , and given the orientation of the detector, which is determined by a rotation angle and an angulation angle . Figure 2 illustrates the rotation and angulation angles in physical terms. These an8 gles are directly avalaible from the DICOM format with a precision of . A. Geometrical Distortion Estimate To estimate the distortion, a calibration grid with a spacing of 1.0 cm is placed against the II screen. We acquire the grid images from different rotation and angulation angles of the acquisition system. Thus, for each image : , we obtain subpixel
Fig. 2. Anatomical angles, which define the orientation of the detector. (a) illustrates de rotation, denoted by , and (b) the angulation, denoted by .
estimates of crossing points and their corresponding physical position (see appendix A for details). Our aim is to determine the ideal pixel coordinates of the grid for each image. The main drawback is that may have a different coordinate origin for each image. To solve this problem, the authors in [9] use the fact that near the center of the image the distortion is low; they propose to track among the images a crossing point placed near the image center and to assign the origin of coordinates to it. We can define the world coordinate system as placed on the calibration grid. Thus, the world coordinates of the crossing points on the grid will be:
!
$
#
"
(4)
Since we place the calibration grid against the II screen, and thus parallel to the image plane, we can assume that the external 8 parameters are :
&%(') *,+-) . /0* ) . /1 8 * %('8 ) *
.
8
32
.
6
25A4 +76
(5)
is the distance between the II screen and the where * 2 image 254 plane, describes a rotation around the axis, and ( A ) is the displacement parallel to the image plane between the coordinate origin of the grid and the one of the image. From (3), (4) and (5) we obtain the following:
%(') *
;:
+76
+
;:
+76
) . /1*
;:
+86
2
A
) . /0* %(') * 52 4 ;: ;: ;: +76 + 6 7 +76 :9; :9; ? 9; 9 9 9 2 2 /< 2 52 4 , A and = < = > < > =
/<
Let @ we can re-write equation (6) as follows:
@ BA
@ @
@
%(') *,+ ) . /0*?2 @ ) . /0* (% ')* 2 A @
< 3
3
<
(6) , then
C
Therefore, we can compute * the2 ideal undistorted coordinates 2 from parameters @ , , and . These parameters must minimize the following expression:
D
E
;
E FHG I + I K I + I 5J(L 5J
IEEE TMI SPECIAL ISSUE MB0016 — [RUNNING ENHANCED CLASS V1.5]
4
To minimize this expression, parameters and are introduced to linearise the problem as follows [29]:
%(') *
@
+
#
8
8
8
E
;
E F
E ;
#
%
E
# # #
G
J
E F ;
#
( #
G
E
*
I
I J I (J L
E ;
E F I
I I I I
is the total number of detected crossing points in all images. When there is high distortion, we can constrain the estimate of these parameters to the crossing points at the central part of the image. For a FOV of cm, all the crossing points can be used. Once the ideal pixel coordinates are obtained, we can model (and analogously ) as the polynomial 1 : #
G
!
L
#
#
#
#
#
% %
I + I +
For simplicity, we omit here the
%
% % % %
I + I +
%
I J
where corresponds to the static component of distortion and to the orientation-dependent one. Let be the estimated at an arbitrary rotation < and angulation < . Hence, the estimate of is obtained by minimizing the expression:
I + I + I +
E F
% % %
C
I 5 J
superindex, but this process is applied for each image. See appendixB for details on bi-variate polynomial fitting
%
I
I
5J
And for we minimize: %
#
.. .
%
We have chosen polynomials for , , and because most distortion / undistortion models present in Computer Vision literature can be reduced to a bi-variate polynomial [27], [25], [28], [26], [30]. The authors in [12] pay special attention to the suitability of a bi-variate polynomial. The degrees of these polynomials have been empirically estimated (see section IV-A.1). However, the estimated distortion is not always the same, since it is a combination of two components: one is fixed and the other depends on the orientation. We expect to model the distortion variation by a low degree polynomial. Thus, we re define the polynomials , , and as:
%
where is a vector containing the coefficients of the polynomial. These coefficients are determined by minimizing the following expression2 :
E F
E F
and
#
%
which should minimize the expression:
< < J J
L
#
%
C
+ I I K I I
#
!
#
#
%
J
E F ;
#
%
*
E F
%
; E
2 2
E F I
% #
E F
;
E
%
I
(
+
where
< < J #
The least squares approach leads to a set of four linear equations: 8 #
) . /1*
@
Analogously, the polynomial (as ) is given by:
E F
I + I + I +
J
As above, the degree of has been empirically determined (see section IV-A.2). Now our aim is to predict the coefficients of the polynomials , , and for each rotation angle and angulation angle !
. Hence, we acquire different X-ray images from different ; views.; Each image : ; ; ; ! is ; acquired with rotation angle and angulation angle
. For each view, we obtain the & " , represented by the coefpolynomials #" , ; $" , ; %" and ; ; ficient vectors '" , )( " , ' " and ) ( " . Each coefficient of the polynomials #" , $" , %" and $" has its own predictive polynomial. For instance, coefficient I (the * th coefficient of
IEEE TMI SPECIAL ISSUE MB0016 — [RUNNING ENHANCED CLASS V1.5]
5
However, the previous approach minimizes a criterion that does not have a good physical interpretation (see [26] for details). We use the obtained solution as starting point to minimize:
from projections
Fig. 3. Approximated reconstruction
# "
;
and
.
) is predicted by the following polynomial: # # #
I '" G
I '"
L
#
'"
#
I # #
<
<
% % % % %
#
!
!
@
I
I + I +
F I I
/ 3
I
I A
@
/ 3
I I
F I + I +
I
I ! ! $# I I where corresponds to the *% row of matrix . /
3
+
J
A problem which arises is that automatic correspondence determination may be computationally expensive and it must deal with ambiguities. On the other hand, manually marking corresponding points can become a tedious task if high accuracy is needed. We showed in [21] that both problems could be reduced by using biplane snakes.
I I A
,
I
First, the user manually determines two pairs of corresponding points near the beginning and the end of the vessel to be reconstructed. Then, an initial 3D curve is obtained from these corresponding points. Let @B be the 3D curve. Its deformation is guided by an energy minimization procedure. The general form of the energy to be minimized for a snake has the following expression:
&
*) F ' (& ' (& ',+ (& 6 F ' where (& is the internal energy and preserves smooth',+ (& is the external energy and attracts the snake ness, and to the features in the image. Applying the Euler-Lagrange equa
A
F
A
@
OF VESSELS
Now that we have corrected the distortion of the acquisition system we can perform the 3D reconstruction of the coronary vessels. The usual method to do the 3D reconstruction is based on corresponding points in both X-ray images [15], [16], [17], [18], [19], [20], [13], [3]. Theoretically, the 3D reconstruction of a point from two / / projections
?3 and
?3 corresponds J J J to the intersection of projective rays and J J , but in practice they may fail to intersect. Dumay et al. [17] propose as an approximated reconstruction (also called retroprojection) the point , which minimizes the distance to both projective lines (see figure 3). However,+ this method does+ not necessarily minimize the dis and and are the tances , where J J J projections of . Moreover, it is not clear how to apply it when more than two views are available. , D ! I We propose to compute as follows: let I I I be the projection matrix for view * . We can derive the followingF constraint for each view F using equation (3):
J
A. Biplane Snakes
C
where #" is a vector containing the coefficients of the predic and
5 . The tive polynomial, degree of these polynomials is also empirically determined (see section IV-A.3). Once geometrical distortion is calibrated, the biplane imaging geometry can be determined by any of the methods proposed in the literature [17], [18], [20], [13]. III. 3D R ECONSTRUCTION
I
I "! ! $#
+
%
F
%
J J
6
%
F
E F
where corresponds to the row of matrix . For any two views, we have 4 equations and 3 unknowns. Any linear least-squares technique can be used to solve this problem.
tion, we get [22]:
.- / &10 +
@
-8
-J 3 + 54 - J 2&100 @
A
In our case, the aim is 3D reconstruction of the curve so that its projections coincide with the vessels in the X-ray images. + To this end, we define the external force A of the 3D snake curve as a function of distances of the projected snake to the image features:
3 +
6 7 (&
:9 82+ (& < ; +6
