Reconstruction of Three-Dimensional Blood Vessel Model Using Fractal Interpolation

Fractal method is used in the image processing and studying the irregular and the complex shapes in the image. It is also used in the reconstruction and smoothing of one-, two-, and three-dimensional data. In this chapter, we present an interpolating fractal algorithm to reconstruct 3D blood vessels. Firstly, the proposed method determines the blood vessel centerline from the 2D retina image, and then it uses the Douglas-Peucker algorithm to detect the control points. Secondly, we use the 3D fractal interpolation and iterated function systems for the visualization and reconstruction of these blood vessels. The results showed that the obtained reduction rate is between 71 and 94% depending on the tolerance value. The 3D blood vessels model can be reconstructed efficiently by using the 3D fractal interpolation method.


Introduction
Since a long time, the interpolation is used to model and visualize data. Indeed, a way to analyze experimental data is to represent them on a 2D or 3D graph; it fills this gap by adding intermediate points that could simulate a more detailed experiment. Our proposed approach is to achieve a 3D fractal interpolation, more adapted to many natural phenomena. We use the fractal interpolation to perform the image restoration. The principle of this technique is explained hereafter [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18].
In this research, we are interested in the 3D geometrical models of human organs [2][3][4][5][6][7][8][9][10][11]; the three-dimensional reconstruction of the vascular networks presents a major medical interest for the diagnostic and prognostic monitoring of several diseases such as atheromatous disease. However, the disadvantage of threedimensional reconstruction approaches is too costly in terms of storage capacity and transmission time. The objective of this research is to create a technique that allows generating with very little information at the beginning, a large amount of data with few errors. To accomplish this goal, we propose an original method based on mathematical morphology in order to provide a 3D reconstruction taking into account all the data of the problem [11][12][13][14][15][16][17][18].
A lot of available works focuses on 3D reconstruction by the 3D fractal interpolation. These approaches were used to generate a 3D model with a high degree of realism.
Guedri et al. [11] have proposed a method for the 3D model blood vessel reconstruction by using the 3D iterated function systems (IFS) and 3D fractal interpolation, and the approach is within the framework of 2D/3D techniques. Sun [12,13] has presented the principles of the multifractal interpolation surface, and he described the methods of dividing local interpolation neighborhoods and determining multiple vertical compression ratios. Then, he gave a practical MATLAB program for the multifractal interpolation surface, and he explained the main parameters in his proposed program. Chen and Bi [14] have created fascinating fractal scenes by using the 3D IFS, and they proposed an efficient coloring, lighting, and mist effect scheme. Guérin et al. [15] have used a fractal model called projected iterated function system (IFS) model that allows the extension of the iteration space to a barycentric space R n2 by enriching the classical IFS model with a set of control points for approximating smooth or rough surfaces defined in R 3 . He et al. [16] have proposed a probabilitybased method to speed up the fractal interpolation execution to three-dimensional terrain reconstruction by a few sparse points in the digital elevation model (DEM). Huang and Chen [17] have adopted two sets of real-world 3D terrain profile data to precede data reducing and then reconstruct them through 3D fractal reconstruction. Xiong [18] has analyzed the reasons that fractal can be used in 3D terrain surface reconstruction and introduced a fractional Brownian motion. Then, they presented an interpolating algorithm to reconstruct 3D terrain surface.
The chapter is organized as follows. In the first step, we give the image source and introduce its features. Then, we present preprocessing of the 2D image datasets and the 3D reconstruction model. Then, a linear simplification method is provided, the Douglas-Peucker method, in order to detect the control points and reduce the data. Thereafter, we develop a reconstruction algorithm which adds new data points by 3D fractal interpolation; this algorithm generates more data points to restore the 3D original model. To evaluate the results of the fractal interpolation method, some blood vessel samples extracted from a real retinal image are used to make the fractal interpolation.

Algorithm
The aim of this research is to reconstruct a 3D blood vessel model from a retinal image by using the 3D fractal interpolation. Figure 1 illustrates the different steps of our approach.
1. First, after giving the source of the used images, we will present the different vessel process detections (preprocessing of the 2D image): a. Binarization image b.Centerline detection c. Point classification 2. Then, we will detect the characteristic points by using the Douglas-Peucker algorithm.
3. Subsequently, we will present the 3D reconstruction model. 4.Finally, the fractal interpolation will be presented.

Image source
There are several databases of available retinal images. Most popular databases used for blood vessel segmentation approach and containing low-resolution images are the Digital Retinal Images for Optic Nerve Segmentation (DRIONS) Database [19], Structured Analysis of the Retina (STARE) [20], and Digital Retinal Images for Vessel Extraction (DRIVE) [21]. On the other hand, there are other databases such as MESSIDOR Digital Retinal Images (MESSIDOR) [22], Retinal Identification Database (RIDB) [23], Glaucoma Database (DB) [24], and high-resolution fundus (HRF) with higher resolutions [25]. In this chapter, we used images illustrated from the databases STARE and HRF.

High-resolution fundus (HRF)
The database contains three image sets: 15 images of healthy patients, 15 images of patients with diabetic retinopathy, and 15 images of glaucomatous patients. All images are captured using a fundus camera CANON CF-60UVi equipped with CANON EOS 20D digital cameras; the capture angle is 60°(FOV). The resolution is 3504 Â 2336 pixels. In addition, this database provided binary gold standard vessel segmentation for each image. Also the masks determining the field of view (FOV) are provided for particular datasets [25]. Figure 2 shows an example of an image from the database.

Vessel process detection
In the first step, our approach starts with the extraction of vessel centerlines, i.e., the preprocessing phase.

Binarization phase
The binarization phase is used to convert the gray image to a binary image using the thresholding technique [28][29][30][31][32]. The intensity Io (x, y) of the binary image is transformed using the following equation: where Ii (x, y) is the intensity of the initial image.

Retinal vessel centerline detection
One of the very useful applications in image processing and in shape recognition is skeletonization. The skeletons represent interconnected lines at the center of an object. They do not only describe the shape but also some mathematical properties of objects, such as length or surface [33][34][35]. The skeletons of the blood vessels are produced by a thinning technique [35][36][37][38]. In this study, the skeleton makes it possible to represent in a compact manner the properties of the blood vessels of the human retina in its particular shape. Its purpose is to describe each vessel by a set of infinitely fine lines. To achieve this objective, we have used the thinning method; the basic idea of this method is to reduce an object iteratively. The algorithm of this method consists of two main steps: the first step is detecting the contour of the object, and in the second step, the contour points are deleted only if they do not change the object's topology.

Feature point extraction
In order to extract the shape characteristic points of the curve, we used a matrix 3 Â 3, and we calculate the point neighbor numbers (PNN) of a pixel Pi containing the information (Pi and the neighboring points must be equal to 1) [39][40][41][42]: • If PNN = 0, P is an isolated point.
• If PNN = 2, P is a connection point.
• If PNN = 3, P is a bifurcation point.
• If PNN ≥ 4; P is a crossing point.

Determination of the control points
Instead of using all the points of blood vessel centerline, it is worthwhile to be contented only with the control points.
Several algorithms are available to archive this reduction [43][44][45][46]. In this work, we will use the Douglas-Peucker algorithm since it gives the best results according to the works of White [43]. This algorithm is still used in many applications.
The Douglas-Peucker method is an automated algorithm, and its principle is given by the following steps: We connect the start pixel to the end pixel of the blood vessel curve (Figure 4a) by a straight line (denoted by S; see Figure 4b).
Then, we compute the perpendicular distance (PD) of each point of the curve to the line S ( Figure 4b).
We compare PD with tolerance ε, and we determine the maximal perpendicular distance MPD (Figure 4b).
If MPD > ε, then we take this point as a starting point of the new interval and repeat step 1 (Figure 4b and c).
This process is repeated until all the distances become smaller than a tolerance value ε (Figure 4b-d, and f). Figure 4 illustrates the evolution of the progressive determination of the control points.

3D reconstruction
The proposed method relied on on the natural shape of the blood vessel, it has a deformable cylindrical shape, which can be represented as a combination of an axis and a surface [8,47,48]; more precisely, we used the approach the right generalized cylinder state model (RGC-sm); then to create the cylindrical shape, we need a set of successive 3D circles with its axes being the blood vessel center (as shown in Figure 5). The RGC parametric equations include a coupling of an ordinary axis H  and a regular surface S (radius parameter rand w is the surface azimuthal parameter) to represent the cylinder as a volumetric object: Figure 6a shows the control points obtained for the blood vessel represented in Figure 5.
It is clear that the number of control points is greatly reduced, and therefore the image transmission time will be reduced with a limited information loss.
Actually, the 3D reconstruction is not performed to all the centerline, but it is limited only to the control points ( Figure 6b).

Fractal interpolation
The principle of fractal interpolation studied by Barnsley [49,50] is to construct a continuous fractal function, passing by a number of giving points of the shape {(x n ; F n ): n = 0, 1, 2,…, N} with x0 < x1 < x2 < … < xN. An interpolation function corresponding to this data set is a continuous function f: [x0; xN] passing with the interpolation points (x n ; F n ) and checking f(x n ) = F n with n = 0, 1, 2,…, N. In this section, we will study constructions of 3D IFS, whose attractors, f, are the graphs of continuous functions, with affine transformation W n [50][51][52][53]. We use the same principle to reconstruct our 3D model by using the 3D fractal interpolation method.

Affine transformation
The affine transformations affect the rotations, translations, scaling, and shear of a data set. He generates each point in a space R n to another point in the R n . The affine transformation contains two major parameters, A and t, where he represented in the AP + T form. The affine transformations on a data set in 2D space are is defined by the following equation [50][51][52][53]: , and P i x i y i the data set of control points (2) The equation system provides five equations for five parameters, so d i , the vertical scaling factor, is computed by using the fractal dimension (DF) calculated by the box-counting algorithm [54,55]. We can solve the above equations for a i , c i , d i , e i , F i which are defined as where N is the affine transformation map number.

Moving from 2D to 3D
In the 2D affine transformation, it only requires a single equation in AP + t form. On the other hand, the 3D affine transformation composes two 2D affine transformations, where this is represented by W1(P) = A1P+ t1 and the other is represented by W2(P) = A2P + t2; all the elements of a 2D transformation matrix could be solved and used to generalize the 3D transformations. Indeed, 3D space (x-y-z) is divided into two 2D spaces (x-y space and x-z space). To determine the 3D affine transformation parameters, we calculate, at first, the parameters of W1 in the space x-y. For that, the parameter bi (bi = 0) is eliminated, so that the y-axis has no term of rotation. The function W1 is defined as follows: Similarly for the second affine transformation W2 in the space x-z, we use the same approach:

Iterated function system (IFS)
Considering W N is the 3D affine transformation (W 1 ,…, W n ), we select an initial set circle which can be selected at random; then we compute iteratively the new data sets. In practice, we start with an the initial circle set and a whole group of affine functions to generate the first set circle. Then, using this circle set and a whole group of affine functions to generate the second circle set, we continue to generate circle set until the generated circle set conforms the requirements of the desired shape, as can be seen in Figure 7 [53].

Results
In order to test the method presented in this work, a demo has been designed and coded in MATLAB. This program has been run on a workstation equipped with an Intel Pentium B960 CPU at 2.20 GHz and 4GB of RAM processor. The reconstruction time takes into account the time required for image segmentation and fractal interpolation. Figure 8 provides two examples of blood vessel skeletons; they are created from the retinal images shown in Figures 2 and 3. These images have been transformed to binary image using thresholding algorithm (Figure 8a and c). Then, the thinning technique is used to produce a blood vessel skeleton that preserves all geometric and topological features of these vessels (Figure 8b and d).

Image segmentation
After performing the skeleton of the retinal vasculature, we use the proposed method for the detection and extraction of the feature points. Figure 9 illustrates an example for the detection of endpoints and bifurcation points. In these examples, the red pixels correspond to the endpoints, and the blue pixels correspond to the bifurcation points.
In addition, Table 1 presents some typical results of endpoints and bifurcation point detection for two test images (image from STARE database and image from

Database
Images Endpoints Bifurcation points Table 1.
Detection of endpoints and bifurcation points for image from high-resolution image database (HRF) and image from STARE database.
HRF database). According to the results, we can see that the images from HRF database give a number of endpoints and bifurcation points higher than the images from the other used database; it is due to the image quality and the presence of a higher number of blood vessels in these images.

Douglas-Peucker algorithm
In the first step, the aim is to compress data by using the Douglas-Peucker method to calculate the data-reduction rate (DR). We use the following equation:

3D interpolation fractal
The control point's density is large in the regions having a large curvature. The original shape of the blood vessel is preserved. These control points are used to reconstruct the 3D blood vessel via fractal interpolation. The affine transformations W n are performed out. The execution time, the number of interpolated points, and the error are calculated versus the iteration number.
In addition, to analyze the performance of our proposed interpolation algorithm, we calculate the error rate between the original image and the interpolated image in the equation below: where NEP is the number of erroneous pixels and TP is the total pixel. The number of erroneous pixels is the pixel numbers that do not appear in the interpolated image and the total pixels are the numbers of pixels in the original image. Table 3 shows the evolution of the iteration number, as well as the error between original and interpolated image and the number of iterations for two test images used; for a small iteration number (N = 50), the error is about 7%. Moreover, For N= 400 iterations, the error decreases to 2.3% for the first image 12 (im0139.vk. jpg) and 2% for the second image used (02_h.tif) For N = 1000, the minimum error rate (0.3%) is obtained from the two test images.
The algorithm described in this document takes a lower time of 2 s to find the curve for a blood vessel and the 3D reconstruction. Moreover, the remnants of execution times are devoted to 3D fractal interpolation. The execution times were between 4 and 14 s. For a small iteration number (N = 50), the execution times is about 4 s. By looking at Table 3, we can see that the value of the execution time increases in the last three tests; this is justified by the increased number of interpolated point, for example, for N = 1000 iterations, the execution time is about 14 s.
Increasing the number of iterations reduces the error but dramatically increases the number of interpolated points and execution times. A trade-off should be looked from these two parameters.

Discussion
The aim of this work is not only to make a 3D reconstruction for visualization but also to describe the geometric shape-based fractal interpolation. The benefits of the proposed method to perform 3D reconstruction by 3D fractal interpolation method lies mainly in reducing storage memory for 3D blood vessel images. Indeed, the data-reduction rate ranges from 72 to 94% for the first test image and between 73 and 95% for the second image used, when the tolerance value varies from ε = 0.5 to 2 pixels. This allows us to reduce costs and transmission time. Another advantage in terms of execution time is obtained. This execution time range is varying from 4 to 14 s, according to the iteration number. Furthermore, an accurate form of description is reached since the algorithm does not only interpolate the original points but also determines the interconnection between these points. However, the weakness of our algorithm is the space complexity of algorithm calculation. It requires a specialized computer to minimize the execution time.
In previous works, the result obtained by Guedri et al. [11] offers an innovative approach for 3D reconstruction with fractal calculation time between 4 and 40 s. The execution time of the present method is about 4-14 s for the same number of iterations. Thus, we can note that the present method is faster than other methods.
Regarding the error evolution of the curvature by using the 3D fractal interpolation, Guedri et al. [11] obtained an error between 0.8 and 8%. The new method is indeed better than the other method in terms of error. The error value is between 0.27 and 6.7% and between 0.25 and 6.5% for the first test image and for the second test image, respectively. This comparison shows that our proposed method can correctly follow all the 3D branches of the blood vessels and reduces the execution time while having a minimal error.

Conclusion
In this work, we propose a new method for reconstructing a 3D vascular tree model from the human retina image by fractal interpolation. Firstly, we describe the method to extract the human retina vascular tree, and then we present a 3D reconstruction algorithm of the blood vessel. In the second step, we applied the Douglas-Peucker algorithm of simplification to detect control points and to compress the data in this 3D model. In the last step, the 3D fractal interpolation and control points are used to generate new data points and restore the original 3D model. From the obtained result, it is found that the Douglas-Peucker method has a high reduction ratio between 72 and 94% for the first test image and between 73 and 95% for the second image. According to the obtained results after the application of the proposed 3D fractal interpolation, we find that the error is small and ranging between 0.2 and 2% for the iteration number greater than 400 iterations. On the other hand, the error increases while decreasing the iteration number. The advantage of this method is that there are more interpolated data points than the original data points; therefore, the reconstruction model is powerful and provides valuable information.