Three-dimensional reconstruction of blood vessels extracted from retinal fundus images

M. Elena Martinez-Perez; Arturo Espinosa-Romero

doi:10.1364/OE.20.011451

1. Introduction

The quantitative assessment of retinal vasculature provides useful clinical information to assist in the diagnosis of various diseases. The detection and measurement of retinal vasculature can be used to quantify the severity of disease and the progression of therapy. Retinal blood vessel tree geometry, topology and vessel tracking have been studied by means of digital image processing mainly using retinal images which are also known as fundus images [1, 2]. However, the majority of these works have been carried out in 2D fundus images (Fig. 1(a) and 1(c)). A first effort to obtain a 3D view of fundus images was developed by [3,4]. Their work focuses on the reconstruction and display of 3D fundus patterns using branching vessel point correspondences between images of two and multiple views. Another approach to reconstruct the retinal fundus is presented by [5] where an epipolar geometry with a plane-and-parallax method is used. An approach to actually reconstruct the skeleton of a blood vessel tree is presented by [6] where a self-calibration based on essential matrix is used. We have presented before a previous work with a first approximation to the 3D reconstruction of blood vessel tree surfaces [7], some improvements will be presented and discussed herein.

Fig. 1 Fundus images, (a) and (c) are two different views from the same eyeball and (b) and (d) are their respective segmented binary images.

Download Full Size | PDF

Photogrammetric analysis of features in images of human ocular fundus is affected by various sources of error, like aberrations of the fundus camera and the eye optics. The magnification in a fundus image is equal to the focal length of the fundus camera divided by the focal length of the eye. This formula can only be applied under different ametropic conditions and changes in the camera position with respect to the subject eye. This difference in magnification between one fundus image and another from different subjects has been pointed out in [8]. As a result, it is not possible to make a direct comparison between measurements from different subjects. Consequently, the 3D metric reconstruction of retinal blood vessel trees is a great challenge.

We have applied well known computer vision techniques to get a nearly automatic approximation to a 3D reconstruction of fundus vasculature. The method considers the following procedures: 1) the estimation of the intrinsic camera parameters for each “fundus camera-eyeball” system, 2) the extraction of blood vessel trees from 2D fundus images, 3) the match of the correspondences of vessel skeleton tree points between images, 4) the computation of the fundamental matrix from the correspondences, 5) the computation of the camera projection matrices, 6) the computation of a point in the space via linear triangulation for each correspondence point, and finally, 7) construction of a 3D surface model of the blood vessel tree.

In our previous results [7] the intrinsic camera parameters were compute as in [9], radial distortions were not considered and correspondence points were selected by hand, furthermore fundus images where adquired with 30° filed of view. We have improved those results by using a pre-processing stage where: 1) 2 new sets of digital fundus images were acquire by a CCD camera attached to the fundus camera using a filed of view of 50°, in order to get larger vessel tree views; 2) the optical radial distortion of the fundus camera lens was automatically corrected, 3) the intrinsic parameters of the fundus camera were obtained by a classical calibration method with a regular pattern and were used later as the initial values for the auto-calibration technique to find the intrinsic parameters of the whole “fundus camera-eyeball” system; and finally, 4) both the correspondences and the elimination of outliers were performed automatically. Figure 2 shows a flowchart depicting the whole 3D reconstruction process.

Fig. 2 Flowchart of the whole 3D reconstruction process.

Download Full Size | PDF

The rest of the paper is organised as follows: Section 2 describes the techniques used for distortion lens correction, calibration of fundus camera and auto-calibration of the “fundus camera-eyeball” system. Section 3 deals with the blood vessel extraction from 2D fundus images and the analysis of the binary images resulted. Section 4 explains the computation of the projection matrices and the linear triangulation procedure to get 3D views of vessel tree skeletons. Section 5 depicts an approach to the surface reconstruction, and finally Sections 6 and 7 present some experimental results and conclusions.

2. System calibration

The calibration task we are dealing with is particularly difficult since we ignore the intrinsic and extrinsic parameters of the optical system formed by the cornea, lens and vitreous humour of the eyeball being examined. Thus, we need to use an auto-calibration technique for the whole “fundus camera-eyeball” system, which uses a number of specific correspondences of branching blood vessel points. These points are taken automatically from a set of retinal images captured from different views of the same eyeball. A pre-processing stage has been introduced first to eliminate the radial distortions of the fundus camera lens and second to get the fundus camera intrinsic parameters through a classical technique, in order to use them as the initial values for the auto-calibration technique.

2.1. Fundus camera distortion correction and calibration

The images obtained with the fundus camera are affected by radial and tangential distortion, an artifact produced by the camera optical system. Since the 3D reconstruction methods we use in our work assume a lineal pin-hole camera model, it is necessary to eliminate the distortion. For this effect we use a classical polynomial model adjusted by using a least mean square method [10].

For the calibration of the fundus camera, we need to find the camera calibration matrix, K, which is a 3 × 3 matrix having the well-known form:

K = [\begin{matrix} α_{u} & - α_{u} \cot θ & u_{0} \\ 0 & α_{v} / \sin θ & v_{0} \\ 0 & 0 & 1 \end{matrix}]

where α_u and α_v correspond to the focal distances in pixels along the axes of the image, θ is the angle between the two image axes and (u₀, v₀) are the coordinates of the principal point. The method we use to find the coordinates of the principal point and the focal distances is that of Kanatani [11], which will be briefly described in following Sections.

2.1.1. Principal point

The principal point of the fundus camera, or optical centre is calculated by finding the focus of expansion (FOE) of a set of spatial reference points while moving the camera along the optical axis. The FOE is defined as the vanishing points of the trajectories of points moving on the scene. If we move the camera along the optical axis, the trajectories of the reference scene points will be parallel to it and therefore will intersect the image plane at the principal point.

We used the corners of a frame as the spatial reference points, similar to the one shown in Fig. 3(a), which was held perpendicular to the camera main axis. Different pictures of the reference frame were captured while carefully moving the fundus camera forward along its main axis. By tracking the frame corners, finding the lines connecting corresponding corners in different images, and finding their intersections, we found an estimation of the camera principal point.

Fig. 3 Examples of the vanishing point estimation process. (a) shows an image of the calibration frame captured by the fundus camera, (b) shows the edges of the calibration frame obtained from the image in (a), and (c) shows the rectangle fitted to the calibration frame in red lines, and the two pairs of lines whose vanishing points are used to estimate the focal length which are shown in blue.

Download Full Size | PDF

2.1.2. Focal distance

In order to estimate the camera focal distance, also called focal length, we rely on the fact that two space lines are orthogonal to each other if and only if their vanishing points are conjugate on the image plane (see Fig. 3(c)). This method considers that the pixels of the camera sensor used to capture the images are square, hence α_u = α_v.

Two points in the image plane whose homogeneous coordinates are the vectors m and m′ are conjugate if and only if: m · m′ = 0. The homogeneous coordinates of the point m in the image plane can be defined as the vector [x,y, f]^⊤. In terms of Cartesian coordinates two points in the image plane with coordinates (a,b) and (a′,b′) are conjugate if and only if

a a^{'} + b b^{'} + f^{2} = 0

Since we do not know the value of the focal length, f, the homogeneous coordinates do not provide an explicit expression for the equivalent Cartesian coordinates on the image plane. For this reason we express conjugacy in terms of the change of a vector m as the focal length of the system changes. Since the homogeneous coordinates of the point m defined before, represents an image point with respect to the focal length f, the point m′ represents the same image point with respect to the focal length f′, and it is given by

m^{'} = {[m_{1}, m_{2}, (f^{'} / f) m_{3}]}^{⊤}

By substituting Eq. (3) into Eq. (2) we obtain the following expression: $m_{1} m_{1}^{'} + m_{2} m_{2}^{'} + \frac{{\hat{f}}^{2}}{f^{2}} m_{3} m_{3}^{'} = 0$ , where m and m′ are the vectors of the vanishing points of two mutually orthogonal space lines, f is the tentative focal length, and f̂ is the true focal length. From the last expression it follows that:

\hat{f} = f \sqrt{- \frac{m_{1} m_{1}^{'} + m_{2} m_{2}^{'}}{m_{3} m_{3}^{'}}}

2.2. Auto-calibration of the “fundus camera-eyeball” system

In order to introduce the eye optics into the calibration process we assume that an auto-calibration technique will include both the eye and the fundus camera optics together. The latter is because the feature points for calibration are taken directly from the fundus images.

Although the eyeball and the camera are separate units, each of them capable of independent movement, the variation on their relative position for the images captured in this article is small, due to the way the images are captured. The computation performed in the method described in this article depends on the number of correspondences present between the image pairs; this obliges to have a large overlap area between images. The procedure to capture the fundus images is as follows: the subject maintains the eye fixed on a beacon set in the space in front of the camera (a fixed flashing red light) by the camera operator (see Fig. 4). The operator moves the camera around the subject eye in order to frame the subject’s retina area to be captured. To be able to have large overlaping areas, the camera principal axis should not move more than certain degrees. The latter is estimated based on the fact that the FOV (field of view) used is 50°, we calculate while taking the images, that camera principal axis does not move more than a quarter of the FOV, this is approximately no more than 12.5°. Since movements are minimal and in search of simplicity, we consider safe to asume a black box “fundus camera-eyeball” system with fixed intrinsic parameters.

Fig. 4 Procedure to capture the fundus images.

Download Full Size | PDF

The method we applied for camera auto-calibration is based on the work reported by [9] which employs a simplification of the Kruppa equations. The method relies solely on the singular value decomposition (SVD) of the fundamental matrix that reduces the number of equations to be solved. It also avoids to recover noise-sensitive quantities such as epipoles, since its accurate estimation is difficult in the presence of noise and/or degenerate motions.

Kruppa equations can be considered as an epipolar matching constraint for the projections of quadratics or conics [12]. The image conics are identical when the images are captured with a camera having fixed intrinsic parameters, which is the case of any pair of views taken with the fundus camera to the same eyeball. Let ω be the projection of the absolute conic. The matching constraint can be expressed in the following form:

F ω^{*} F^{T} = {[e^{'}]}_{\times} ω^{*} {[e^{'}]}_{\times}

where F is the fundamental matrix of the two views, e′ is the second epipole. The equality is up to scale (ω^* = ω⁻¹, dual image of the absolute conic) and [e′]_× is the skew-symmetric matrix associated with the cross-product of e′. Equation (5) is called the Kruppa equations.

Using singular value decomposition (SVD) of matrix F = UDV^T, Eq. (5) is equivalent to:

\frac{u_{2}^{T} ω^{*} u_{2}}{r^{2} v_{1}^{T} ω^{*} v_{1}} = - \frac{u_{1}^{T} ω^{*} u_{2}}{r s v_{1}^{T} ω^{*} v_{2}} = \frac{u_{1}^{T} ω^{*} u_{1}}{s v_{2}^{T} ω^{*} v_{2}}

where r,s are the eigenvalues of the matrix FF^T ; u₁, u₂, u₃ are the column vectors of U ; v₁, v₂, v₃ are the column vectors of V, and U, V are two orthogonal matrices. The aim is to find ω^* from Eq. (6). ω^* can be expressed as well in terms of the calibration matrix, K, as:

ω^{*} = K K^{T}

therefore K is extracted from Eq. (7) by computing K^T employing the Cholesky decomposition of ω⁻¹, then it is transposed and finally inverted to yield K.

1. Initial solution and point correspondences. Unlike [9] algorithm, the initial camera parameters we used are taken from the previous stage of fundus camera classical calibration described in Section 2.1. Assuming we have M images, that have been acquired with constant camera intrinsic parameters, a total of N ≤ (M(M − 1))/2 fundamental matrices can be defined to be used for the auto-calibration method.

Feature points are selected from each 2D image automatically using the skeletons of the segmented images, looking for those points which connectivity is either a bifurcation or a vessel crossing point (see Fig. 5(a), 5(d) and [1] for more details). The correspondence between image feature points are obtained by a classical normalised correlation technique to get the putative matches, followed by the solution of a nonlinear minimisation problem base on the least-median-of-squares (LMedS) of the distances between putative points to eliminate the outliers.

Fig. 5 Process to extract a tree from binary image. (a) and (b) show skeleton and significant points marked over the original image, (b) and (e) show the automated tracking of the tree selected by the operator and (c) and (f) show the extraction of the tree segments from the binary image. Upper raw: subject set 1, lower raw: subject set 2.

Download Full Size | PDF

For every feature point in image 1 we extract a window of normalised data of size m×m and correlate it with a window of the same size corresponding to every feature point in image 2. The correlation matrix is as follow:

Cor = \frac{W_{1} * W_{2}}{\sqrt{\sum \sum W_{2}^{2}}}

where W₁ and W₂ are the windows of data from image 1 and 2 respectively which centres are the feature points. W₁ should be also a window with normalised data

W_{1} = W / (\sqrt{\sum \sum W^{2}})

. This correlation matrix holds the correlation strength of every point relative to every other point. The latter is needed since we want to find pairs of points that correlate maximally in both directions (rows and columns).

After finding the putative matches, we eliminate the outliers using the least median of squares (LMedS) of the distances between matched putative points. The LMedS method estimates the parameters by solving the nonlinear minimisation problem which in one dimension reduces to:

min_{\hat{θ}} {med}_{i} r_{i}^{2}

where r_i = (y_i − θ̂) is the residual of the difference between what is actually observed, y_i, and what is estimated, θ̂. The robust standard deviation is based on the minimal median and multiplied by a finite correction factor for the case of normal errors:

\hat{σ} = 1.4826 (1 + 5 / (n - p)) \sqrt{{med}_{i} r_{i}^{2} (\hat{θ})}

where n is the number of points and p is the dimension (p=1). This scale factor is then used to determine a weight k_i for the ith observation namely:

k_{i} = {\begin{array}{l} 1 & if | r_{i} / \hat{σ} | \\ 0 & otherwise \end{array}

The correspondences having k_i = 0 are outliers and should not be further taken into account. The reader is refer to [13] for further explanation about robust estimators. The LMedS method to eliminate outliers in this way works properly for our application since rotation between fundus images is minimal. In cases with a very large rotation this method does not work.

2. Non-linear optimisation. Matrix ω^*, from Eq. (5), is computed as the solution of a nonlinear least squares problem.

Let π_ij(S_F, ω^*) denote the differences of ratios ij and $σ_{i j}^{2} (S_{F}, ω^{*})$ the variance from Eq. (6). Matrix ω^* is computed as the solution of the non-linear least squares problem:

ω^{*} = {argmin}_{{\tilde{ω}}^{*}} \sum_{i = 1}^{N} \frac{π_{12}^{2} (S_{F_{i}}, {\tilde{ω}}^{*})}{σ_{π 12}^{2} (S_{F_{i}}, {\tilde{ω}}^{*})} + \frac{π_{13}^{2} (S_{F_{i}}, {\tilde{ω}}^{*})}{σ_{π 13}^{2} (S_{F_{i}}, {\tilde{ω}}^{*})} + \frac{π_{23}^{2} (S_{F_{i}}, {\tilde{ω}}^{*})}{σ_{π 23}^{2} (S_{F_{i}}, {\tilde{ω}}^{*})}

where S_F is a vector formed by parameters from the SVD of F. The method used to find F is based on the one proposed in [14]. For a detailed description of this method please refer to [9]).

3. Blood vessel extraction

Blood vessels are automatically segmented and extracted from 2D fundus images using a previously described algorithm based on multi-scale analysis [2]. Two geometrical features based upon the first and the second derivative of the intensity image, maximum gradient and principal curvature, are obtained at different scales by means of Gaussian derivative operators. A multiple pass region growing procedure is used which progressively segments the blood vessels using the feature information together with spatial information about the 8-neighbouring pixels. Figures 1(a) and 1(c) show two different views a fundus image and Fig. 1(b) and 1(d) their segmented binary images, respectively. The optic disc region in this pair of images is on the bottom-left, blood vessels are tracked from this area outwards.

3.1. Analysis of binary images

A semi-automatic method to measure and quantify geometrical and topological properties of continuous vascular trees in clinical fundus images was developed on a previous work [1]. These measurements are made from the segmented binary images. The skeletons of the segmented trees are produced by a thinning technique, branching and crossing points are identified and segments of the trees are labelled and stored as a chain code. Figures 5(a) and 5(b) show the skeleton lines and the significant points over imposed in the original image. The operator selects the tree to be measured and decides if it is an arterial or venous tree. An automatic process then tracks the tree and measures the lengths, areas and angles of individual vessel segments. Figures 5(b) and 5(e) show the tracking of the selected tree. Geometrical data such as diameter, length, branching angle, and the connectivity information between vessel segments from continuous retinal vessel trees are tabulated. A number of geometrical properties and topological indices are derived. Note that in Fig. 5(c) and 5(f) only the vessel segments are extracted, branching areas are ignore. From each tree an ASCII table of measurements is generated. The blood vessel width information will be used later for vessel tree surface modelling. First raw of Fig. 5 corresponds to subject set 1, whereas second raw corresponds to subject set 2.

4. 3D skeleton reconstruction

Suppose a set of points of correspondence x_i ↔ x′_i is given, such as the skeleton points of a tree. It is assumed that these correspondences come from a set of 3D points X_i, which are unknown. The reconstruction task is to find the camera projection matrices P and P′, as well as the 3D points X_i such that

\begin{array}{l} x_{i} = {P X}_{i} \\ x_{i}^{'} = P^{'} X_{i} & for all i \end{array}

The reconstruction process consist of: 1) compute the fundamental matrix F from the correspondences of tree skeleton points, 2) compute the camera projection matrices, P and P′, from fundamental matrix, and 3) for each pair of correspondence points x_i ↔ x′_i compute the point in the space that projects these two image points by linear triangulation [12, 15].

5. Surface modelling

Since the 3D reconstruction algorithm described in Section 4 is applied only to the skeleton points of each blood vessel tree, we did not recover 3D information from the actual shape of its surface. We considered a fair assumption to describe the blood vessel segment cross-section as a circle in order to generate a 3D surface. Each blood vessel segment was modelled by a generalised cylinder [16] which axis is the 3D skeleton estimated earlier. The radius of each circular cross-section is kept constant for each vessel segment, and it is obtained from the measures taken from the 2D fundus images [1]. The cylinders are scaled using the camera intrinsic parameters from the matrix K of the “fundus camera-eyeball” system.

6. Experimental results

Two sets of 10 fundus images each were taken belonging to two different subjects. Both sets of retinal images were taken using a fundus camera with a 50° field of view (Zeiss FF 450 IR). The first set of images were acquired with a CCD Sony camera (Sony Power HAD 3CCD Colour Video Camera) attached to the fundus camera. Images were 768 × 576 pixels in size. For the second set of images a CCD Canon camera (Canon EOS 5D Mark II) was used. Images were 783 × 522 pixels in size. For the purpose of these experiments, 8 different views were taken for the auto-calibration process and one different pair for the reconstruction example, for each set. The distortion of all the images captured was corrected and the fundus camera was manually calibrated using the procedure described in Section 2.1. Under the general assumption of having a sensor with square pixels and the angle between image axes θ = π/2, the following fundus camera calibration matrix, K_mi, was obtained as an initial condition for the auto-calibration technique. Since each set of data (for each subject) was taken with different CCD camera, we do compute a K_mi matrix for each one:

\begin{matrix} K_{m 1} = [\begin{array}{r} 909.68 & 0 & 391.28 \\ 0 & 909.68 & 295.61 \\ 0 & 0 & 1 \end{array}]; & K_{m 2} = [\begin{array}{r} 639.79 & 0 & 460.75 \\ 0 & 639.39 & 247.02 \\ 0 & 0 & 1 \end{array}] \end{matrix}

K_m1 for fundus camera with Sony camera and K_m2 for fundus camera with Canon camera.

6.1. Camera auto-calibration matrix

For each set of images, eight different views of the fundus images were chosen for auto-calibration purposes. The branching points were detected from fundus images automatically, from the skeletons as explained in Section 3.1, and were matched using a normalised correlation technique followed by a LMedS estimator as described in Section 2.2.

Since we have 8 images acquired with constant camera intrinsic parameters, a total of 28 fundamental matrices can be defined. However due to the high disparity of some of these combination of pairs we undertake a further pair selection by using only those pairs which ROI (region of interest) overlapped was larger than 70%. The latter is to ensure we have enough correspondence points. After this selection they turn out to be 12 pairs. The average overlap area for pairs in set 1 is 87.7%, and for pairs in set 2 is 90.9%.

The parameters for correlation consider for both experiments were: a window size of 25×25 and a distance for searching window dmax = 300 pixels. The choice of these parameters depend on the image and feature size. To compensate for brightness differences between images, before applying the correlation, we subtract a smoothed image with an averaging filter of same size as the correlation window from each of the images. Figures 6(a) and 6(b) show one of these image pairs with the feature points, Fig. 6(c) shows the putative points after correlation and Fig. 6(d) shows the inlier points after LMedS minimisation.

Fig. 6 Matching process: (a) feature points on image 1 (n=55), (b) feature points on image 2 (n=59), (c) putative matches after normalised correlation (n=36) and (d) inlier points after LMedS minimisation (n=29). Images taken from set 1.

Download Full Size | PDF

Finally, we assumed the initial parameters for auto-calibration equal to K_mi. Solving the system 5 we obtained the calibration matrix, K_i, for the whole “camera-eyeball” system as:

\begin{matrix} K_{1} = [\begin{array}{r} 998.81 & - 1 & 372.74 \\ 0 & 905.54 & 317.42 \\ 0 & 0 & 1 \end{array}]; & K_{2} = [\begin{array}{r} 462.89 & 0.0090 & 344.16 \\ 0 & 515.60 & 329.64 \\ 0 & 0 & 1 \end{array}] \end{matrix}

K₁ for the first set and K₂ for the second set.

Each auto-calibration matrix K_i shown above was the one with the smaller covariance matrix from those similar to the matrix K_mi. We presume that this difference between the fundus camera matrix and auto-calibrated matrix, might be due the optics of the eyeball, which is unknown and was not considered in the model of the system.

6.2. 3D skeleton reconstruction

One pair of images was selected and automatically segmented from each set of images [2]. Skeletons of the extracted blood vessel trees, as those shown in Fig. 5(c) a vein from set 1 and 5(f) a vein from set 2, are marked and measured. Measurements are kept in an ASCII table per tree for further morphological analysis [1].

Since continuous blood vessel tree skeletons are marked they can be tracked along the tree. Using these tables, the correspondences (x_i ↔ x′_i) of a continuous tree skeletons from a pair of images are extracted automatically. Figures 7(a) and 7(b) show the same tree extracted from images Fig. 1(b) and 1(d), taken from set 1, respectively. Figure 7(c) shows the corresponding matched tree skeletons. Note that a small branch at the top-centre of the tree, shown in Fig. 7(a), is missing in image 2 of the same eye shown in Fig. 7(b). This small branch is of course automatically discarded in the correspondences, as it is shown in Fig. 7(c). Figures 7(d) and 7(e) show the same tree taken from images of set 2 and Fig. 7(f) corresponds to the matched tree skeletons.

Fig. 7 (a) T2 extracted tree from image 1 from set 1, (b) T2 extracted tree from image 2 from set 1 and (c) matched skeleton points between the two images. (d), (e) and (f) correspond to the same example using a pair of images from set 2.

Download Full Size | PDF

A set of correspondences x_i ↔ x′_i for only one tree is needed in order to compute the two projection matrices P and P′, as described in Section 4, with the use of the camera auto-calibration matrix K_i calculated in Section 6.1. In Fig. 8 the epipolar lines corresponding to a subset of matching points for each of the two image pairs used to reconstruct the vascular trees are shown.

Fig. 8 (a) and (b) show the epipolar lines drawn on the image pair from subject 1, used to recover the 3D vascular tree, while (c) and (d) show the same from the image pair from subject 2.

Download Full Size | PDF

Once the projection matrices are obtained, the triangulation of the rest of correspondences for each skeleton tree extracted from the pair of images is performed. A total of 2 trees were reconstructed from each set. Correspondences of the skeleton points were smoothed by a Gaussian kernel and interpolated using splines in order to obtain a smoother and refined version of the reconstructed skeletons. Figure 9 shows different 3D views of the blood vessel skeletons extracted and plotted on the same coordinate system. First raw of Fig. 9 corresponds to set 1 and second raw to set 2.

Fig. 9 (a) and (b) different 3D views of skeleton trees reconstruction from set 1. (c) and (d) same example from set 2.

Download Full Size | PDF

Figures 9(b) and 9(d) show the curvature of the eyeball for each subject, respectively, these views correspond to about 50° field of view of the fundus camera. Figures 9(a) and 9(c) are just rotated versions of the formers.

6.3. Surface visualisation

Using the surface model described in Section 5, and based on the 3D skeleton points and vessel segment diameter, the surface of each vessel segment is generated. Each tree is visualised using Geomview (http://www.geomview.org/), a 3D viewing program. Each tree is built as a set of disconnected segments such as those shown in Fig. 7(a) and 7(d), each of them with a defined spatial position and orientation.

Each vessel segment surface is described by a mesh of square patches wrapped in one direction so that it forms a cylinder. The mesh vertices for each segment are calculated as follows: for each 3D skeleton point, i.e. the set that describes the generalised cylinder axis, a polygon or cross–section is considered. Each polygon is centred along the axis position in space, and its orientation is set so that the vector normal to the polygon centre is tangent to the skeleton axis at that point. The coordinates of the polygon vertices are used as the mesh vertices. Each vertex is connected to the two neighbouring vertices in the same cross–section and to the corresponding vertices in the neighbouring cross-sections. The vertices of the cross-section at the extreme points of each segment are only connected to one neighbouring cross-section. Three 3D views of different surface model of blood vessel trees are shown in Fig. 10.

Fig. 10 Blood vessel surfaces extracted from the pair of images. (a), (b) and (c) correspond to set 1; (d), (e) and (f) to set 2. Diameter measure of each vessel segment is obtained from the corresponded 2D binary images. (c) and (f) show both trees from each set respectively, where the curvature of each eyeball can be seen.

Download Full Size | PDF

Figure 10(a) corresponds to the surface reconstruction of the tree shown in Fig. 5(c) a vein. Figure 10(b) corresponds to the surface reconstruction of the tree shown in Fig. 7(a) an artery and Fig. 10(c) shows both trees together but rotated to show the eyeball curvature of about 50°. Note that in Fig. 10(a) and 10(b) the optic disc is at the top left of the coordinate system, and blood vessel segments coming from that area are wider than those from the tips of the tree, as vessels naturally tapered. Figure 10(d) is a vein shown inf Fig. 5(f), Fig. 10(e) is an artery from the same pair and Fig. 10(f) are both trees together. Figures 10(a), 10(b) and 10(c) correspond to the set 1 whereas Fig. 10(d), 10(e) and 10(f) correspond to the set 2.

7. Conclusions

We have presented a 3D model of a retinal blood vessel tree, reconstructed from two different views of fundus images through well known computer vision techniques: camera lens distortion correction, correspondence matching, auto-calibration procedure, 3D projection, triangulation and surface meshing.

By contrast, [3, 4] do not use an auto-calibration procedure but instead they employed two parallel planes in front of the fundus camera and used a transparent acrylic plate with a grid in order to calibrate the camera. Additionally, they projected the correspondences and the fundus pattern into a spherical surface which means that even a vessel crossing point belongs to the same surface. In the case of [5] they reconstructed the retinal fundus as a whole without distinguishing between the retinal fundus itself and blood vessels. Finally [6] used a self-calibration approach based on the essential matrix, where they only reconstruct the skeleton lines of blood vessels using a 30° FOV, therefore the eyeball curvature is not clear. A more realistic assumption, however can be seen from Fig. 10(c) and 10(f) in which all vessel trees surfaces follow the same eyeball curvature but they do not belong to the same spherical surface. Moreover, we reconstructed blood vessel trees rather than fundus patterns.

The core of the reconstruction presented in this paper is the auto-calibration technique which was based on the work of [9]. The authors have shown experimental results from extensive simulations and several image sequences which demonstrate the effectiveness of their method in accurately estimating the intrinsic calibration matrices. Their evaluation was done with syntetic and real data taken from buildings. Unlike [9], our initial camera parameters for the auto-calibration algorithm were calculated using a classical camera calibration algorithm [11], the latter makes our auto–calibration implementation faster to converge, since the initial value is nearer to the solution.

There are various questions that still have to be addressed in order to use this model for Euclidean measurable applications. The most significant is the evaluation of the current blood vessel surfaces, this is a difficult process because unlike buildings, retinal fundus has the difficulty in accessing grand truth measurements. A possible source of grand truth measurements can be obtained from ultrasound or optical coherence tomography (OCT) data, therefore future work will attempt to do so. Another issue to be addressed is that since trees are extracted from individual vessel segments separately, branching regions are not considered in the mesh of the current model. At the present Fig. 9 and 10 only show projective views of the reconstruction. Finally, larger eyeball views could be obtained by using multiple view geometry methods. Nonetheless our results are promising and, after clearing up these issues, the 3D model obtained could be used for realistic applications.

The representation of vascular trees with 3D models offers many advantages for ophthalmology: it can provide a more realistic view for clinical and for educational purposes, it can permit to extend all the geometrical and topological properties already measured in 2D images [1] to a 3D model and finally, it can be used for blood flow simulation.

Acknowledgments

This work was supported in part by CONACYT Grant CB-2007-83088 and in part by SEP PIFI 2011-31MSU0098J-15.

References and links

1. M. E. Martinez-Perez, A. D. Hughes, A. V. Stanton, S. A. Thom, N. Chapman, A. A. Bharath, and K. H. Parker, “Retinal vascular tree morphology: A semi-automatic quantification,” IEEE Trans. Biomed. Eng. 49, 912–917 (2002). [CrossRef]

2. M. E. Martinez-Perez, A. D. Hughes, S. A. Thom, A. A. Bharath, and K. H. Parker, “Segmentation of blood vessels from red-free and fluorescein retinal images,” Med. Image Anal. 11, 47–61 (2007). [CrossRef] [PubMed]

3. K. Deguchi, J. Noami, and H. Hontani, “3d fundus pattern reconstruction and display from multiple fundus images,” in Proceedings 15th International Conference on Pattern Recognition (IEEE, 2000), pp. 94–97. [CrossRef]

4. K. Deguchi, D. Kawamata, K. Mizutani, H. Hontani, and K. Wakabayashi, “3d fundus shape reconstruction and display from stereo fundus images,” IEICE Trans. Inf. Syst. E83-D, 1408–1414 (2000).

5. T. E. Choe, G. Medioni, I. Cohen, A. C. Walsh, and S. R. Sadda, “2-D registration and 3-D shape inference of the retinal fundus from fluorescein images,” Med. Image Anal. 12, 174–190 (2008). [CrossRef]

6. D. Liu, N. Wood, X. Xu, N. Witt, A. Hughes, and S. Thom, “3D reconstruction of the retinal arterial tree using subject-specific fundus images,” in Advances in Computational Vision and Medical Image Processing, J. M. R. S. Tavares and R. M. N. Jorge ed. (Springer, 2009), pp. 187–201. [CrossRef]

7. M. E. Martinez-Perez and A. Espinosa-Romero, “3D Reconstruction of Retinal Blood Vessels From Two Views,” in Proceedings of the 4th Indian Conference on Computer Vision, Graphics and Image Processing, B. Chanda, S. Chandran, and L. Davis, ed. (Indian Statistical Insitute, 2004), pp. 258–263.

8. J. Arnold, J. Gates, and K. Taylor, “Possible errors in the measurement of retinal lesions,” Invest. Ophthalmol. Vis. Sci. 34, 2576–2580 (1993). [PubMed]

9. M. I. A. Lourakis and R. Deriche, “Camera self-calibration using the singular value decomposition of the fundamental matrix: from point correspondences to 3d measurements,” Research Report 3748, INRIA Sophia-Antipolis, (1999).

10. M. Sonka, Image Processing, Analysis, and Machine Vision (Thomson, 2008).

11. K. Kanatani, Geometry Computation for Machine Vision (Oxford Science Publications, 1993).

12. R. Hartley and A. Zisserman, Multiple View Geometry in Computer Vision (Cambridge Uiversity Press, 2000).

13. P. J. Rousseeuw and A. M. Leroy, Robust Regression and Outilier Detection (John Wiley & Sons, 1987). [CrossRef]

14. Z. Zhang, R. Deriche, O. Faugeras, and Q.-T. Luong, “A robust technique for matching two uncalibrated images through the recovery of the unknown epipolar geometry,” Research Report 2273, INRIA Sophia-Antipolis, (1994).

15. R. I. Hartley, “Estimation of relative camera positions for uncalibrated cameras,” in Proceedings of the 2nd European Conference on Computer Vision, G. Sandini, ed. (Springer-Verlag, 1992), pp. 579–587.

16. J. Ponce, “Straight homogeneous generalized cylinders: Differential geometry and uniqueness results,” Int. J. Comput. Vis. 4, 79–100 (1990). [CrossRef]

Three-dimensional reconstruction of blood vessels extracted from retinal fundus images

Abstract

1. Introduction

2. System calibration

2.1. Fundus camera distortion correction and calibration

2.1.1. Principal point

2.1.2. Focal distance

2.2. Auto-calibration of the “fundus camera-eyeball” system

3. Blood vessel extraction

3.1. Analysis of binary images

4. 3D skeleton reconstruction

5. Surface modelling

6. Experimental results

6.1. Camera auto-calibration matrix

6.2. 3D skeleton reconstruction

6.3. Surface visualisation

7. Conclusions

Acknowledgments

References and links

Supplementary Material (6)

Cited By

Figures (10)

Equations (15)

Optics Express