- Email: [email protected]
Skeleton extraction based on anisotropic partial differential equation
Optik 126 (2015) 3692–3697
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.de/ijleo
Skeleton extraction based on anisotropic partial differential equation Fang Zhang a , Zhitao Xiao a,∗ , Lei Geng a , Jun Wu a , Zhenbei Xu a , Danyu Wang a , Jiangtao Xi b a b
School of Electronics and Information Engineering, Tianjin Polytechnic University, Tianjin 300387, China School of Electrical, Computer and Telecommunications Engineering, University of Wollongong, Wollongong, NSW 2522, Australia
a r t i c l e
i n f o
Article history: Received 2 September 2014 Accepted 25 August 2015 Keywords: Skeleton extraction Partial differential equation Gradient vector field Divergence
a b s t r a c t Skeleton extraction is of great significance for image recognition processes. In this paper, a novel skeleton extraction method is proposed based on the anisotropic diffusion of the gradient vector field. Firstly, the gradient vector is calculated for the original image. Then, the gradient vector field is updated by an anisotropic partial differential equation to strengthen the source point or sink point of the image energy. At last, according to the divergence property of the final gradient vector field, the skeleton of the initial image is obtained. The experimental results demonstrate that the proposed method is robust for noise and can give the exact skeleton for images. © 2015 Elsevier GmbH. All rights reserved.
1. Introduction
conditions of the vector field (u, v), (u, v) is diffused with time by the following:
With the development of the image processing technology, image processing has been widely applied in many fields, such as biomedical [1,2]. As an enabling technology, skeleton extraction of images has been an active topic of research in recent years, both in 2D and 3D [3,4]. Many different methods have been proposed on this topic, such as the topological thinning method [5], the hierarchical approach based on Voronoi diagrams [6,7], the techniques based on the distance maps [8,9] and so on. However, all these methods require pre-processing of the initial images, e.g. enhancement, filtering, segmentation and binarization, and hence the effectiveness of these approaches is significantly impacted by the quality of the initial images and the performance of these preprocessing operations. Recently, partial differential equation (PDE) has been found as an effective tool for image processing and computer vision [10,11]. The basic idea of PDE-based methods is to set the image to be processed as the initial condition of a PDE with respect to time t, to deform the image by solving the PDE, and to extract the steady state solution of the PDE as the results of the processed image [12]. For the task of skeleton extraction, Yu and Bajaj [13] proposed a method based on anisotropic vector diffusion. In [13] f(x, y) = || ∇ G (x, y) ∗ I0 (x, y)||2 denotes the edge strength map of the original image I0 (x, y), where G (x, y) stands for a Gaussian kernel. With ∇ f(x, y) being the initial
⎧ ∂u ⎪ ⎨ ∂t = · div(g(˛) · ∇ u) − (u − fx ) fx2 + fy2
∗ Corresponding author. Tel.: +86 02283955639. E-mail address: [email protected] (Z. Xiao). http://dx.doi.org/10.1016/j.ijleo.2015.08.189 0030-4026/© 2015 Elsevier GmbH. All rights reserved.
⎪ ⎩ ∂v = · div(g(˛) · ∇ v) − (v − f ) f 2 + f 2 x x y
(1)
∂t
where ∇ = (∂/∂x)i + (∂/∂y)j, fx and fy are the first order spatial derivatives, g(·) is a decreasing function and ˛ is the angle between the central vector and the surrounding vectors. Using non-maximal suppression and double-threshold, the final skeletons are traced from the skeleton strength map [7], which is calculated from the diffused vector field (u, v) and provides a measure on the possibility for each pixel falling on the skeletons. Direkoglu [14] proposed to diffuse image I in the dominance of direction normal to the feature boundaries while allowing tangential diffusion to contribute slightly, i.e. ∂I = cI + I ∂t
(2)
where c is a constant, denotes the direction normal to the feature boundary (the gradient direction), and is the perpendiculargradient direction, which is shown in Fig. 1 and given by the following: =
(−Iy , Ix ) ∇⊥I = |∇ I| Ix2 + Iy2
=
∇I = |∇ I|
(I , Iy )
x
Ix2 + Iy2
(3)
F. Zhang et al. / Optik 126 (2015) 3692–3697
3693
2.2. The diffusion of the gradient vector field The GVF should be updated with the aim to smooth the noise. This function is similar to various techniques of image smoothing [15]. In order to remove the noise and simultaneously keep the image edges, an anisotropic PDE is introduced to update u0 (x, y) and v0 (x, y) respectively. Firstly, an inner orthogonal coordinate system is introduced, which is characterized by the local gradient direction () and the perpendicular-gradient direction (), respectively (shown in Fig. 1). Then, based on such coordinate system, an anisotropic PDE can be obtained as follows: Fig. 1. The inner orthogonal coordinate system.
∂I = I + g(|∇ G ∗ I|)I ∂t
where I and I are the second order spatial derivatives given as: I = |∇ I|div
∇I
|∇ I|
I = ∇ 2 I − |∇ I|div
∇I
∇ I
(4)
Then the skeleton strength map (SSM), which provides the likelihood to be a skeleton point, is obtained by computing the mean curvature of level-sets as follows: SSM = ∇ ·
∇I
|∇ I|
(5)
Although both methods in [13,14] work directly on the grayscale image, the quality of the skeleton extracted is unsatisfactory. In [13], Eq. (1) is used to update the two gradient components of the edge strength map of the original image, and the skeleton is extracted based on the updated gradient components. The final result will depend on the initial edge and hence usually exhibits discontinuities. In [14], the initial image is diffused by an anisotropic PDE (i.e., Eq. (2)) first, and noise can be removed by the grayvalue diffusion. As the skeleton is obtained based on the smoothed image, it is usually continuous, but many redundant branches can be resulted, especially at the end of the skeletons. In this paper, a novel method is proposed for the extraction of the skeletons of biomedicine images. The proposed method can be applied directly on to the images, and thus pre-processing is not requires. In contrast to existing works, for example Refs. [13,14], the proposed method is advantageous by skeleton continuity and no branches. The paper is organized as follows. In Section 2, the novel approach is described in detail. Then, in Section 3, skeleton obtained by the proposed method is presented in comparison with other existing approaches. Finally, conclusion is given in Section 4.
(7)
where I(x, y) is to be updated and I0 is its initial value. I and I are the second order spatial derivatives along and respectively and can be expressed as I =
Ixx Iy 2 − 2Ix Iy Ixy + Iyy Ix 2 2
Ix + Iy
2
I =
Ixx Ix 2 + 2Ix Iy Ixy + Iyy Iy 2 Ix 2 + Iy 2
(8)
Also in (7) G is the Gaussian filter, g(·) is a smooth nonincreasing function given by g(s) = k2 /(k2 + s2 ) [16] which will influence the diffusion speed, and k is a constant. The diffusion process is carried out based on the following ideas: when (x, y) is located at the edge of the image, ∇ G * I(x, y) will exhibit a large value, and g(| ∇ G ∗ I(x, y)|) will be a small value, and thus the image edge will be protected; if (x, y) is corrupted by isolated noise, the value of ∇ G * I(x, y) will be small, resulting in a large value of g(| ∇ G ∗ I(x, y)|), and hence the noise can be smoothed. In other words, fast diffusion will take place around noisy pixels, thus suppressing the noise, and slow diffusion occurs along the local gradient direction, and hence edges can be kept. Then, according to Eq. (7), u0 and v0 can be updated by
⎧ ∂u ⎪ = u + g(|∇ G ∗ u|)u u(x, y, 0) = u0 (x, y) ⎨ ∂t
⎪ ⎩ ∂v = v + g(|∇ G ∗ u|)v ∂t
(9)
v(x, y, 0) = v0 (x, y)
and the final GVF is y) = u(x, y)i + v(x, y)j F(x,
(10)
2.3. Skeleton extraction based on the divergence analysis Divergence measures the magnitude of the source or sink of a vector field at a given point (x, y) [17]. The divergence of a contin uously differentiable vector field A(x, y) = P(x, y)i + Q (x, y)j is:
2. Skeleton extraction based on the anisotropic partial differential equation The main process of the proposed method includes three steps, i.e. initialization of the gradient vector field (GVF) calculation from original image, GVF update by an anisotropic PDE, and skeleton extraction through divergence analysis on the updated GVF. 2.1. The initial gradient vector field Let I0 : R2 → R represents a gray-level image to be processed, and I0 (x, y) is the gray-level of pixel (x, y). The initial GVF F0 can be obtained by F0 (x, y) = u0 (x, y)i + v0 (x, y)j
I(x, y, 0) = I0 (x, y)
(6)
where u0 (x, y) = ∂I0 (x, y)/∂x and v0 (x, y) = ∂I0 (x, y)/∂y are the two components of F0 (x, y), i and j are the unit coordinate vectors.
∇ · A(x, y) =
∂Q (x, y) ∂P(x, y) + ∂x ∂y
(11)
When the divergence of a point is greater than zero, i.e. ∇ · y) > 0, this point has a positive source which emits flux in the A(x, vector field. When the divergence is less than zero, i.e. ∇ · A(x, y) < 0, this point has a negative source and the flux is sunk in the vector field. When the divergence is equal to zero, i.e. ∇ · A(x, y) = 0, there is no source at this position. Let us consider the following gradient vector field: ∂u(x, y) ∂v(x, y) ϕ(x, y) = ∇ · F(x, y) = + ∂x ∂y
(12)
The above can be used to determine the source points and sink points, corresponding to the lower-gray region and the higher-gray
3694
F. Zhang et al. / Optik 126 (2015) 3692–3697
Fig. 2. Experimentation with respect to noise in the binary image. Four groups of images from up to down give the skeleton extraction results on the noiseless image, the image with Gaussian noise (mean = 0, standard deviation = 0.2), the image with salt and pepper noise (density = 0.1), the image with salt and pepper noise (density = 0.1) and Gaussian noise (mean = 0, standard deviation = 0.2), respectively. The images from the left are the original images, the superimposition of the skeletons for brighter regions onto the initial images, the initial GVFs and the updated GVFs by PDE respectively.
region, respectively. By means of thresholding to ϕ(x, y), two binary images can be obtained:
3. Experimental results and analysis 3.1. Discretization of the evolution equation
ϕD (x, y) =
1
ϕ(x, y) > 0
0
ϕ(x, y) < 0
ϕL (x, y) =
0
ϕ(x, y) > 0
1 ϕ(x, y) < 0
,
(13)
Then, the single pixel-wide skeleton for the dark object in the original image can be obtained through thinning for ϕD (x, y) (morphologic thinning for example). Similarly, the skeleton for the light object can be extracted by thinning to the binary image ϕL (x, y). Based on the above, the proposed method can be summarized as follows:
Step 1: Calculate the initial GVF of the given image F0 (x, y) using Eq. (6); Step 2: Implement the anisotropic diffusion PDE in Eq. (9) on the two components u0 and v0 , and obtain the updated GVF F(x, y); y) using Eq. (12); Step 3: Calculate the divergence of F(x, Step 4: Get the skeleton according to the divergence property of y) by means of thresholding using Eq. (13) and thinning. F(x,
In order to solve Eq. (9) numerically, the equation has to be discretized. The initial image I0 (x, y) is represented by M × N matrix of intensity values. Let ui,j , vi,j denote u(x, y), v(x, y). The evolution equations give images at tn = nt (n is iterations). In gradient vector field, denote u(i, j, tn ) by uni,j and v(i, j, tn ) by vni,j . The time derivative ut (ut = ∂u/∂t) at (i, j, tn ) is approximated by the forward difference
(ut )ni,j =
un+1 − uni,j i,j
(14)
t
Here t is step-size in time domain. The spatial derivatives for u(i, j, tn ) are
⎧ uni+1,j − uni−1,j n ⎪ (ux )i,j = ⎪ ⎪ 2 ⎨
n
(uy )i,j =
n
(uxx )i,j = uni+1,j − 2uni,j + uni−1,j
⎪ ⎪ ⎪ ⎩ (u
n xy )i,j
=
2
(uyy )i,j = uni,j+1 − 2uni,j + uni,j−1
uni+1,j+1 − uni−1,j+1 − uni+1,j−1 + uni−1,j−1 4
n
uni,j+1 − uni,j−1
(15)
F. Zhang et al. / Optik 126 (2015) 3692–3697
3695
Then, using Eq. (8), (u )ni,j and (u )ni,j can be expressed as
⎧ n n 2 n n n n n 2 (uxx )i,j [(uy )i,j ] − 2(ux )i,j (uy )i,j (uxy )i,j + (uyy )i,j [(ux )i,j ] ⎪ n ⎪ ⎨ (u )i,j = 2 2 n n [(ux )i,j ] + [(uy )i,j ]
⎪ ⎪ ⎩ (u )ni,j =
n
n
2
n
n
n
n
n
(uxx )i,j [(ux )i,j ] + 2(ux )i,j (uy )i,j (uxy )i,j + (uyy )i,j [(uy )i,j ] n
2
n
[(ux )i,j ] + [(uy )i,j ]
2
(16)
2
(vt )ni,j , (v )ni,j and (v )ni,j can be expressed in a similar way to Eqs. (14)–(16). n denote g(| ∇ G * u(i, j, t )|) = k2 /(k2 + | ∇ G * u(i, Let ˛ni,j and ˇi,j n
j, tn )|2 ) and g(|∇ G ∗ v(i, j, tn )|) = k2 /(k2 + |∇ G ∗ v(i, j, tn )|2 ), respectively. Then the discrete form of Eq. (9) is gained as un+1 = uni,j + t[(u )ni,j + ˛ni,j (u )ni,j )] i,j n vn+1 = vni,j + t[(v )ni,j + ˇi,j (v )ni,j )] i,j
(17)
For solving Eq. (9) discretely, three parameters are needed, including discrete time step t, iteration time n, and diffusion speed parameter k. In the following, the values of these parameters are determined empirically. 3.2. Skeleton results and analysis The experimentation with respect to noise in the binary image is first presented. Fig. 2(a-1) is a binary image of size 50 × 50 pixels without any noise from Ref. [14], and Fig. 2(b-1), (c-1) and (d-1) show the noisy images with Gaussian noise, with salt and pepper noise, and with both salt and pepper noise and Gaussian noise respectively. For each of the cases, four images are presented from the left to right, namely the original image, the superimposition of the skeleton for brighter regions onto the initial image, the initial GVFs and the updated GVFs by the anisotropic PDE (Eq. (9)) with t = 0.2, n = 50, and k = 10, respectively. From Fig. 2, it can be seen that although the initial GVFs are corrupted by noise in Fig. 2(b-3), (c-3) and (d-3), the updated GVFs shown in Fig. 2(b-4), (c-4) and (d-4) are regular and very similar to the updated GVF in Fig. 2(a-4) which is obtained from the noiseless image. To visual inspection, the extracted skeletons of the subjects in Fig. 2(b-2), (c-2) and (d-2) are detected without obvious dislocation and without missing any part of the subject. Even for Fig. 2(d-1), in which the input image has very poor quality, the skeleton (Fig. 2(d-2)) is extracted with only small perturbations. The proposed method is also applied to extract the skeletons for gray images with the results shown in Figs. 3(a) and 4(a) are two MRI brain images. Figs. 3(b) and 4(b) show the skeletons obtained for the brighter regions by our method, where Fig. 3(b) is obtained with t = 0.2, n = 120, and k = 10 and Fig. 4(b) is extracted with t = 0.2, n = 100, and k = 10. For comparison, the results in [13] are also shown in Fig. 3(c) with many disconnected skeletons. The results in [14] are depicted in Fig. 4(c) containing many redundant branches, especially at the end of the skeletons. In contrast, the proposed method performs much better to extract the skeletons for these complex images of brains. The proposed method differs from the method in [13] in two main aspects. One is that the initial value of the updated components (u, v), i.e. the initial conditions of the PDE are different. Let I0 denote the original image. (u, v) equals to ∇ || ∇ I0 || in Yu’s method at the beginning, while (u, v) is initialized with ∇ I0 in our method. By computing ∇ || ∇ I0 || and ∇ I0 for Fig. 2(a-1) and (c-1) respectively, the results are obtained shown in Fig. 5, where the images from the left are the original images I0 , the GVF of I0 (i.e. ∇ I0 ) and the GVF of || ∇ I0 || (i.e. ∇ || ∇ I0 ||), respectively. It can be seen that both the two initialization form can represent the main feature of the original image, and initializing (u, v) with ∇ || ∇ I0 || results in
Fig. 3. Skeletonization of a MRI brain image. (a) The initial image; (b) the superimposition of the obtained skeletons by the proposed method for brighter regions onto the initial image; (c) the superimposition of the skeletons extracted by Yu et al., which comes from Ref. [13].
higher noise (see Fig. 5(b-3)). Therefore, initializing (u, v) with ∇ I0 is sufficient and more appropriate for GVF diffusion. The other difference is that the divergence of updated GVF is employed by the proposed method to assess the skeletons, while non-maximal suppression and double-threshold were utilized in [13] to the skeleton strength map for the same purpose. Such difference has resulted in elimination of discontinuities in the skeletons extracted. The proposed method can be considered as the improvement of the method in Ref. [14]. However, the proposed method also differs from Ref. [14] in two major aspects. Firstly, PDE is employed to update the GVF, rather than the image itself as in [14]. In other words, the proposed method computes the first order derivative of the updated object, while the method in [14] calculates the second order derivative of the updated object. The later magnifies the
3696
F. Zhang et al. / Optik 126 (2015) 3692–3697
Fig. 5. The initial values of the updated components (u, v) in the proposed method and in Yu’s method. Two groups of images from up to down correspond the noiseless image and the image with salt and pepper noise (density = 0.1), respectively. The images from the left are the original images, the initial values of (u, v) in the proposed method (i.e. (u, v) is initialized with ∇ I0 ) and Yu’s method (i.e. (u, v) is initialized with ∇ || ∇ I0 |), respectively.
method. Therefore, the anisotropic PDE proposed can effectively remove noise and protect feature of the GVF. The experiments show that the method proposed can be applied in the gray-scale image directly without any pre-processing, which can avoid information loss of the original image. 4. Conclusion In this paper, a novel skeleton extraction method based on anisotropic partial differential equation is proposed. The new method is implemented in two main steps, namely, update of the gradient vector field of the original image by an anisotropic partial differential equation and calculation of the divergence of the updated gradient vector field. The proposed method is tested by MRI images, and good results are obtained. The proposed technique does not need any pre-processing of the images, such as enhancement filtering, segmentation and binarization for the original image, thus overcoming the limitation of information loss in the traditional image skeleton extraction method. In addition, the proposed method is robust in suppressing noises. Experiments show that our method is superior to others by skeleton continuity and no branches. Acknowledgments
Fig. 4. Skeletonization of a MRI brain image. (a) The initial image; (b) the superimposition of the obtained skeletons by the proposed method for brighter regions onto the initial image; (c) the superimposition of the skeletons extracted by Direkoglu, which comes from Ref. [14].
This work was supported by National Nature Science Foundation of China (NSFC) under grant No. 61102150, Tianjin Science and Technology Supporting Projection under grant No. 13ZCZDGX02100 and No. 14ZCZDGX00033. References
image details, including noises, and thus leads to many branches at the end of the skeletons. Secondly, the anisotropic PDEs used for image diffusion are different. The PDE used in [14] allows tangential diffusion to contribute slightly and the diffusion speed is constant. In contrast, in the proposed method, the tangential diffusion plays a main role and the diffusion along the normal direction is adjustable, which is controlled by an adaptive gene in the proposed
[1] E. Zhang, F. Wang, Y. Li, X. Bai, Automatic detection of microcalcifications using mathematical morphology and a support vector machine, Bio-Med. Mater. Eng. 24 (2014) 53–59. [2] S. Hu, C. Xu, W. Guan, Y. Tang, Y. Liu, Texture feature extraction based on wavelet transform and gray-level co-occurrence matrices applied to osteosarcoma diagnosis, Bio-Med. Mater. Eng. 24 (2014) 129–143. [3] J. Xu, A generalized morphological skeleton transform using both internal and external skeleton points, Pattern Recognit. 47 (2014) 2607–2620.
F. Zhang et al. / Optik 126 (2015) 3692–3697 [4] L. Rossi, A. Torsello, Coarse-to-fine skeleton extraction for high resolution 3D meshes, Comput. Vis. Image Underst. 118 (2014) 140–152. [5] L. Lam, Thinning methodologies – a comprehensive survey, IEEE Trans. Pattern Anal. Mach. Intell. 14 (1992) 869–885. [6] J.W. Brandt, V.R. Algazi, Continuous skeleton computation by Voronoi diagram, CVGIP: Image Underst. 55 (1992) 329–338. [7] R.L. Ogniewicz, O. Kübler, Hierarchic Voronoi skeletons, Pattern Recognit. 28 (1995) 343–359. [8] G. Malandain, S. Fernández-Vidal, Euclidean skeletons, Image Vis. Comput. 16 (1998) 317–327. [9] M. Couprie, Topological maps and robust hierarchical Euclidean skeletons in cubical complexes, Comput. Vis. Image Underst. 117 (2013) 355–369. [10] E. Nadernejad, M. Nikpour, Image denoising using new pixon representation based on fuzzy filtering and partial differential equations, Digit. Signal Process. 22 (2012) 913–922. [11] X. Liu, L. Huang, Z. Guo, Adaptive fourth-order partial differential equation filter for image denoising, Appl. Math. Lett. 24 (2011) 1282–1288.
3697
[12] C. Tang, L. Han, H. Ren, T. Gao, Z. Wang, K. Tang, The oriented-couple partial differential equations for filtering in wrapped phase patterns, Opt. Express 17 (2009) 5606–5617. [13] Z. Yu, C. Bajaj, A segmentation-free approach for skeletonization of gray-scale images via anisotropic vector diffusion, in: IEEE International Conference on Computer Vision and Pattern Recognition, vol. 1, 2004, pp. 415–420. [14] C. Direkoglu, Feature extraction via heat flow analogy, Dissertation, University of Southampton, 2009. [15] F. Zhang, Z. Xiao, J. Wu, L. Geng, H. Li, J. Xi, J. Wang, Anisotropic coupled diffusion filter and binarization for the electronic speckle pattern interferometry fringes, Opt. Express 20 (2012) 21905–21916. [16] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (1990) 629–639. [17] G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review, Courier Dover Publications, 2000.
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.de/ijleo
Skeleton extraction based on anisotropic partial differential equation Fang Zhang a , Zhitao Xiao a,∗ , Lei Geng a , Jun Wu a , Zhenbei Xu a , Danyu Wang a , Jiangtao Xi b a b
School of Electronics and Information Engineering, Tianjin Polytechnic University, Tianjin 300387, China School of Electrical, Computer and Telecommunications Engineering, University of Wollongong, Wollongong, NSW 2522, Australia
a r t i c l e
i n f o
Article history: Received 2 September 2014 Accepted 25 August 2015 Keywords: Skeleton extraction Partial differential equation Gradient vector field Divergence
a b s t r a c t Skeleton extraction is of great significance for image recognition processes. In this paper, a novel skeleton extraction method is proposed based on the anisotropic diffusion of the gradient vector field. Firstly, the gradient vector is calculated for the original image. Then, the gradient vector field is updated by an anisotropic partial differential equation to strengthen the source point or sink point of the image energy. At last, according to the divergence property of the final gradient vector field, the skeleton of the initial image is obtained. The experimental results demonstrate that the proposed method is robust for noise and can give the exact skeleton for images. © 2015 Elsevier GmbH. All rights reserved.
1. Introduction
conditions of the vector field (u, v), (u, v) is diffused with time by the following:
With the development of the image processing technology, image processing has been widely applied in many fields, such as biomedical [1,2]. As an enabling technology, skeleton extraction of images has been an active topic of research in recent years, both in 2D and 3D [3,4]. Many different methods have been proposed on this topic, such as the topological thinning method [5], the hierarchical approach based on Voronoi diagrams [6,7], the techniques based on the distance maps [8,9] and so on. However, all these methods require pre-processing of the initial images, e.g. enhancement, filtering, segmentation and binarization, and hence the effectiveness of these approaches is significantly impacted by the quality of the initial images and the performance of these preprocessing operations. Recently, partial differential equation (PDE) has been found as an effective tool for image processing and computer vision [10,11]. The basic idea of PDE-based methods is to set the image to be processed as the initial condition of a PDE with respect to time t, to deform the image by solving the PDE, and to extract the steady state solution of the PDE as the results of the processed image [12]. For the task of skeleton extraction, Yu and Bajaj [13] proposed a method based on anisotropic vector diffusion. In [13] f(x, y) = || ∇ G (x, y) ∗ I0 (x, y)||2 denotes the edge strength map of the original image I0 (x, y), where G (x, y) stands for a Gaussian kernel. With ∇ f(x, y) being the initial
⎧ ∂u ⎪ ⎨ ∂t = · div(g(˛) · ∇ u) − (u − fx ) fx2 + fy2
∗ Corresponding author. Tel.: +86 02283955639. E-mail address: [email protected] (Z. Xiao). http://dx.doi.org/10.1016/j.ijleo.2015.08.189 0030-4026/© 2015 Elsevier GmbH. All rights reserved.
⎪ ⎩ ∂v = · div(g(˛) · ∇ v) − (v − f ) f 2 + f 2 x x y
(1)
∂t
where ∇ = (∂/∂x)i + (∂/∂y)j, fx and fy are the first order spatial derivatives, g(·) is a decreasing function and ˛ is the angle between the central vector and the surrounding vectors. Using non-maximal suppression and double-threshold, the final skeletons are traced from the skeleton strength map [7], which is calculated from the diffused vector field (u, v) and provides a measure on the possibility for each pixel falling on the skeletons. Direkoglu [14] proposed to diffuse image I in the dominance of direction normal to the feature boundaries while allowing tangential diffusion to contribute slightly, i.e. ∂I = cI + I ∂t
(2)
where c is a constant, denotes the direction normal to the feature boundary (the gradient direction), and is the perpendiculargradient direction, which is shown in Fig. 1 and given by the following: =
(−Iy , Ix ) ∇⊥I = |∇ I| Ix2 + Iy2
=
∇I = |∇ I|
(I , Iy )
x
Ix2 + Iy2
(3)
F. Zhang et al. / Optik 126 (2015) 3692–3697
3693
2.2. The diffusion of the gradient vector field The GVF should be updated with the aim to smooth the noise. This function is similar to various techniques of image smoothing [15]. In order to remove the noise and simultaneously keep the image edges, an anisotropic PDE is introduced to update u0 (x, y) and v0 (x, y) respectively. Firstly, an inner orthogonal coordinate system is introduced, which is characterized by the local gradient direction () and the perpendicular-gradient direction (), respectively (shown in Fig. 1). Then, based on such coordinate system, an anisotropic PDE can be obtained as follows: Fig. 1. The inner orthogonal coordinate system.
∂I = I + g(|∇ G ∗ I|)I ∂t
where I and I are the second order spatial derivatives given as: I = |∇ I|div
∇I
|∇ I|
I = ∇ 2 I − |∇ I|div
∇I
∇ I
(4)
Then the skeleton strength map (SSM), which provides the likelihood to be a skeleton point, is obtained by computing the mean curvature of level-sets as follows: SSM = ∇ ·
∇I
|∇ I|
(5)
Although both methods in [13,14] work directly on the grayscale image, the quality of the skeleton extracted is unsatisfactory. In [13], Eq. (1) is used to update the two gradient components of the edge strength map of the original image, and the skeleton is extracted based on the updated gradient components. The final result will depend on the initial edge and hence usually exhibits discontinuities. In [14], the initial image is diffused by an anisotropic PDE (i.e., Eq. (2)) first, and noise can be removed by the grayvalue diffusion. As the skeleton is obtained based on the smoothed image, it is usually continuous, but many redundant branches can be resulted, especially at the end of the skeletons. In this paper, a novel method is proposed for the extraction of the skeletons of biomedicine images. The proposed method can be applied directly on to the images, and thus pre-processing is not requires. In contrast to existing works, for example Refs. [13,14], the proposed method is advantageous by skeleton continuity and no branches. The paper is organized as follows. In Section 2, the novel approach is described in detail. Then, in Section 3, skeleton obtained by the proposed method is presented in comparison with other existing approaches. Finally, conclusion is given in Section 4.
(7)
where I(x, y) is to be updated and I0 is its initial value. I and I are the second order spatial derivatives along and respectively and can be expressed as I =
Ixx Iy 2 − 2Ix Iy Ixy + Iyy Ix 2 2
Ix + Iy
2
I =
Ixx Ix 2 + 2Ix Iy Ixy + Iyy Iy 2 Ix 2 + Iy 2
(8)
Also in (7) G is the Gaussian filter, g(·) is a smooth nonincreasing function given by g(s) = k2 /(k2 + s2 ) [16] which will influence the diffusion speed, and k is a constant. The diffusion process is carried out based on the following ideas: when (x, y) is located at the edge of the image, ∇ G * I(x, y) will exhibit a large value, and g(| ∇ G ∗ I(x, y)|) will be a small value, and thus the image edge will be protected; if (x, y) is corrupted by isolated noise, the value of ∇ G * I(x, y) will be small, resulting in a large value of g(| ∇ G ∗ I(x, y)|), and hence the noise can be smoothed. In other words, fast diffusion will take place around noisy pixels, thus suppressing the noise, and slow diffusion occurs along the local gradient direction, and hence edges can be kept. Then, according to Eq. (7), u0 and v0 can be updated by
⎧ ∂u ⎪ = u + g(|∇ G ∗ u|)u u(x, y, 0) = u0 (x, y) ⎨ ∂t
⎪ ⎩ ∂v = v + g(|∇ G ∗ u|)v ∂t
(9)
v(x, y, 0) = v0 (x, y)
and the final GVF is y) = u(x, y)i + v(x, y)j F(x,
(10)
2.3. Skeleton extraction based on the divergence analysis Divergence measures the magnitude of the source or sink of a vector field at a given point (x, y) [17]. The divergence of a contin uously differentiable vector field A(x, y) = P(x, y)i + Q (x, y)j is:
2. Skeleton extraction based on the anisotropic partial differential equation The main process of the proposed method includes three steps, i.e. initialization of the gradient vector field (GVF) calculation from original image, GVF update by an anisotropic PDE, and skeleton extraction through divergence analysis on the updated GVF. 2.1. The initial gradient vector field Let I0 : R2 → R represents a gray-level image to be processed, and I0 (x, y) is the gray-level of pixel (x, y). The initial GVF F0 can be obtained by F0 (x, y) = u0 (x, y)i + v0 (x, y)j
I(x, y, 0) = I0 (x, y)
(6)
where u0 (x, y) = ∂I0 (x, y)/∂x and v0 (x, y) = ∂I0 (x, y)/∂y are the two components of F0 (x, y), i and j are the unit coordinate vectors.
∇ · A(x, y) =
∂Q (x, y) ∂P(x, y) + ∂x ∂y
(11)
When the divergence of a point is greater than zero, i.e. ∇ · y) > 0, this point has a positive source which emits flux in the A(x, vector field. When the divergence is less than zero, i.e. ∇ · A(x, y) < 0, this point has a negative source and the flux is sunk in the vector field. When the divergence is equal to zero, i.e. ∇ · A(x, y) = 0, there is no source at this position. Let us consider the following gradient vector field: ∂u(x, y) ∂v(x, y) ϕ(x, y) = ∇ · F(x, y) = + ∂x ∂y
(12)
The above can be used to determine the source points and sink points, corresponding to the lower-gray region and the higher-gray
3694
F. Zhang et al. / Optik 126 (2015) 3692–3697
Fig. 2. Experimentation with respect to noise in the binary image. Four groups of images from up to down give the skeleton extraction results on the noiseless image, the image with Gaussian noise (mean = 0, standard deviation = 0.2), the image with salt and pepper noise (density = 0.1), the image with salt and pepper noise (density = 0.1) and Gaussian noise (mean = 0, standard deviation = 0.2), respectively. The images from the left are the original images, the superimposition of the skeletons for brighter regions onto the initial images, the initial GVFs and the updated GVFs by PDE respectively.
region, respectively. By means of thresholding to ϕ(x, y), two binary images can be obtained:
3. Experimental results and analysis 3.1. Discretization of the evolution equation
ϕD (x, y) =
1
ϕ(x, y) > 0
0
ϕ(x, y) < 0
ϕL (x, y) =
0
ϕ(x, y) > 0
1 ϕ(x, y) < 0
,
(13)
Then, the single pixel-wide skeleton for the dark object in the original image can be obtained through thinning for ϕD (x, y) (morphologic thinning for example). Similarly, the skeleton for the light object can be extracted by thinning to the binary image ϕL (x, y). Based on the above, the proposed method can be summarized as follows:
Step 1: Calculate the initial GVF of the given image F0 (x, y) using Eq. (6); Step 2: Implement the anisotropic diffusion PDE in Eq. (9) on the two components u0 and v0 , and obtain the updated GVF F(x, y); y) using Eq. (12); Step 3: Calculate the divergence of F(x, Step 4: Get the skeleton according to the divergence property of y) by means of thresholding using Eq. (13) and thinning. F(x,
In order to solve Eq. (9) numerically, the equation has to be discretized. The initial image I0 (x, y) is represented by M × N matrix of intensity values. Let ui,j , vi,j denote u(x, y), v(x, y). The evolution equations give images at tn = nt (n is iterations). In gradient vector field, denote u(i, j, tn ) by uni,j and v(i, j, tn ) by vni,j . The time derivative ut (ut = ∂u/∂t) at (i, j, tn ) is approximated by the forward difference
(ut )ni,j =
un+1 − uni,j i,j
(14)
t
Here t is step-size in time domain. The spatial derivatives for u(i, j, tn ) are
⎧ uni+1,j − uni−1,j n ⎪ (ux )i,j = ⎪ ⎪ 2 ⎨
n
(uy )i,j =
n
(uxx )i,j = uni+1,j − 2uni,j + uni−1,j
⎪ ⎪ ⎪ ⎩ (u
n xy )i,j
=
2
(uyy )i,j = uni,j+1 − 2uni,j + uni,j−1
uni+1,j+1 − uni−1,j+1 − uni+1,j−1 + uni−1,j−1 4
n
uni,j+1 − uni,j−1
(15)
F. Zhang et al. / Optik 126 (2015) 3692–3697
3695
Then, using Eq. (8), (u )ni,j and (u )ni,j can be expressed as
⎧ n n 2 n n n n n 2 (uxx )i,j [(uy )i,j ] − 2(ux )i,j (uy )i,j (uxy )i,j + (uyy )i,j [(ux )i,j ] ⎪ n ⎪ ⎨ (u )i,j = 2 2 n n [(ux )i,j ] + [(uy )i,j ]
⎪ ⎪ ⎩ (u )ni,j =
n
n
2
n
n
n
n
n
(uxx )i,j [(ux )i,j ] + 2(ux )i,j (uy )i,j (uxy )i,j + (uyy )i,j [(uy )i,j ] n
2
n
[(ux )i,j ] + [(uy )i,j ]
2
(16)
2
(vt )ni,j , (v )ni,j and (v )ni,j can be expressed in a similar way to Eqs. (14)–(16). n denote g(| ∇ G * u(i, j, t )|) = k2 /(k2 + | ∇ G * u(i, Let ˛ni,j and ˇi,j n
j, tn )|2 ) and g(|∇ G ∗ v(i, j, tn )|) = k2 /(k2 + |∇ G ∗ v(i, j, tn )|2 ), respectively. Then the discrete form of Eq. (9) is gained as un+1 = uni,j + t[(u )ni,j + ˛ni,j (u )ni,j )] i,j n vn+1 = vni,j + t[(v )ni,j + ˇi,j (v )ni,j )] i,j
(17)
For solving Eq. (9) discretely, three parameters are needed, including discrete time step t, iteration time n, and diffusion speed parameter k. In the following, the values of these parameters are determined empirically. 3.2. Skeleton results and analysis The experimentation with respect to noise in the binary image is first presented. Fig. 2(a-1) is a binary image of size 50 × 50 pixels without any noise from Ref. [14], and Fig. 2(b-1), (c-1) and (d-1) show the noisy images with Gaussian noise, with salt and pepper noise, and with both salt and pepper noise and Gaussian noise respectively. For each of the cases, four images are presented from the left to right, namely the original image, the superimposition of the skeleton for brighter regions onto the initial image, the initial GVFs and the updated GVFs by the anisotropic PDE (Eq. (9)) with t = 0.2, n = 50, and k = 10, respectively. From Fig. 2, it can be seen that although the initial GVFs are corrupted by noise in Fig. 2(b-3), (c-3) and (d-3), the updated GVFs shown in Fig. 2(b-4), (c-4) and (d-4) are regular and very similar to the updated GVF in Fig. 2(a-4) which is obtained from the noiseless image. To visual inspection, the extracted skeletons of the subjects in Fig. 2(b-2), (c-2) and (d-2) are detected without obvious dislocation and without missing any part of the subject. Even for Fig. 2(d-1), in which the input image has very poor quality, the skeleton (Fig. 2(d-2)) is extracted with only small perturbations. The proposed method is also applied to extract the skeletons for gray images with the results shown in Figs. 3(a) and 4(a) are two MRI brain images. Figs. 3(b) and 4(b) show the skeletons obtained for the brighter regions by our method, where Fig. 3(b) is obtained with t = 0.2, n = 120, and k = 10 and Fig. 4(b) is extracted with t = 0.2, n = 100, and k = 10. For comparison, the results in [13] are also shown in Fig. 3(c) with many disconnected skeletons. The results in [14] are depicted in Fig. 4(c) containing many redundant branches, especially at the end of the skeletons. In contrast, the proposed method performs much better to extract the skeletons for these complex images of brains. The proposed method differs from the method in [13] in two main aspects. One is that the initial value of the updated components (u, v), i.e. the initial conditions of the PDE are different. Let I0 denote the original image. (u, v) equals to ∇ || ∇ I0 || in Yu’s method at the beginning, while (u, v) is initialized with ∇ I0 in our method. By computing ∇ || ∇ I0 || and ∇ I0 for Fig. 2(a-1) and (c-1) respectively, the results are obtained shown in Fig. 5, where the images from the left are the original images I0 , the GVF of I0 (i.e. ∇ I0 ) and the GVF of || ∇ I0 || (i.e. ∇ || ∇ I0 ||), respectively. It can be seen that both the two initialization form can represent the main feature of the original image, and initializing (u, v) with ∇ || ∇ I0 || results in
Fig. 3. Skeletonization of a MRI brain image. (a) The initial image; (b) the superimposition of the obtained skeletons by the proposed method for brighter regions onto the initial image; (c) the superimposition of the skeletons extracted by Yu et al., which comes from Ref. [13].
higher noise (see Fig. 5(b-3)). Therefore, initializing (u, v) with ∇ I0 is sufficient and more appropriate for GVF diffusion. The other difference is that the divergence of updated GVF is employed by the proposed method to assess the skeletons, while non-maximal suppression and double-threshold were utilized in [13] to the skeleton strength map for the same purpose. Such difference has resulted in elimination of discontinuities in the skeletons extracted. The proposed method can be considered as the improvement of the method in Ref. [14]. However, the proposed method also differs from Ref. [14] in two major aspects. Firstly, PDE is employed to update the GVF, rather than the image itself as in [14]. In other words, the proposed method computes the first order derivative of the updated object, while the method in [14] calculates the second order derivative of the updated object. The later magnifies the
3696
F. Zhang et al. / Optik 126 (2015) 3692–3697
Fig. 5. The initial values of the updated components (u, v) in the proposed method and in Yu’s method. Two groups of images from up to down correspond the noiseless image and the image with salt and pepper noise (density = 0.1), respectively. The images from the left are the original images, the initial values of (u, v) in the proposed method (i.e. (u, v) is initialized with ∇ I0 ) and Yu’s method (i.e. (u, v) is initialized with ∇ || ∇ I0 |), respectively.
method. Therefore, the anisotropic PDE proposed can effectively remove noise and protect feature of the GVF. The experiments show that the method proposed can be applied in the gray-scale image directly without any pre-processing, which can avoid information loss of the original image. 4. Conclusion In this paper, a novel skeleton extraction method based on anisotropic partial differential equation is proposed. The new method is implemented in two main steps, namely, update of the gradient vector field of the original image by an anisotropic partial differential equation and calculation of the divergence of the updated gradient vector field. The proposed method is tested by MRI images, and good results are obtained. The proposed technique does not need any pre-processing of the images, such as enhancement filtering, segmentation and binarization for the original image, thus overcoming the limitation of information loss in the traditional image skeleton extraction method. In addition, the proposed method is robust in suppressing noises. Experiments show that our method is superior to others by skeleton continuity and no branches. Acknowledgments
Fig. 4. Skeletonization of a MRI brain image. (a) The initial image; (b) the superimposition of the obtained skeletons by the proposed method for brighter regions onto the initial image; (c) the superimposition of the skeletons extracted by Direkoglu, which comes from Ref. [14].
This work was supported by National Nature Science Foundation of China (NSFC) under grant No. 61102150, Tianjin Science and Technology Supporting Projection under grant No. 13ZCZDGX02100 and No. 14ZCZDGX00033. References
image details, including noises, and thus leads to many branches at the end of the skeletons. Secondly, the anisotropic PDEs used for image diffusion are different. The PDE used in [14] allows tangential diffusion to contribute slightly and the diffusion speed is constant. In contrast, in the proposed method, the tangential diffusion plays a main role and the diffusion along the normal direction is adjustable, which is controlled by an adaptive gene in the proposed
[1] E. Zhang, F. Wang, Y. Li, X. Bai, Automatic detection of microcalcifications using mathematical morphology and a support vector machine, Bio-Med. Mater. Eng. 24 (2014) 53–59. [2] S. Hu, C. Xu, W. Guan, Y. Tang, Y. Liu, Texture feature extraction based on wavelet transform and gray-level co-occurrence matrices applied to osteosarcoma diagnosis, Bio-Med. Mater. Eng. 24 (2014) 129–143. [3] J. Xu, A generalized morphological skeleton transform using both internal and external skeleton points, Pattern Recognit. 47 (2014) 2607–2620.
F. Zhang et al. / Optik 126 (2015) 3692–3697 [4] L. Rossi, A. Torsello, Coarse-to-fine skeleton extraction for high resolution 3D meshes, Comput. Vis. Image Underst. 118 (2014) 140–152. [5] L. Lam, Thinning methodologies – a comprehensive survey, IEEE Trans. Pattern Anal. Mach. Intell. 14 (1992) 869–885. [6] J.W. Brandt, V.R. Algazi, Continuous skeleton computation by Voronoi diagram, CVGIP: Image Underst. 55 (1992) 329–338. [7] R.L. Ogniewicz, O. Kübler, Hierarchic Voronoi skeletons, Pattern Recognit. 28 (1995) 343–359. [8] G. Malandain, S. Fernández-Vidal, Euclidean skeletons, Image Vis. Comput. 16 (1998) 317–327. [9] M. Couprie, Topological maps and robust hierarchical Euclidean skeletons in cubical complexes, Comput. Vis. Image Underst. 117 (2013) 355–369. [10] E. Nadernejad, M. Nikpour, Image denoising using new pixon representation based on fuzzy filtering and partial differential equations, Digit. Signal Process. 22 (2012) 913–922. [11] X. Liu, L. Huang, Z. Guo, Adaptive fourth-order partial differential equation filter for image denoising, Appl. Math. Lett. 24 (2011) 1282–1288.
3697
[12] C. Tang, L. Han, H. Ren, T. Gao, Z. Wang, K. Tang, The oriented-couple partial differential equations for filtering in wrapped phase patterns, Opt. Express 17 (2009) 5606–5617. [13] Z. Yu, C. Bajaj, A segmentation-free approach for skeletonization of gray-scale images via anisotropic vector diffusion, in: IEEE International Conference on Computer Vision and Pattern Recognition, vol. 1, 2004, pp. 415–420. [14] C. Direkoglu, Feature extraction via heat flow analogy, Dissertation, University of Southampton, 2009. [15] F. Zhang, Z. Xiao, J. Wu, L. Geng, H. Li, J. Xi, J. Wang, Anisotropic coupled diffusion filter and binarization for the electronic speckle pattern interferometry fringes, Opt. Express 20 (2012) 21905–21916. [16] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (1990) 629–639. [17] G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review, Courier Dover Publications, 2000.