EFFICIENT NUMERICAL SCHEMES FOR GRADIENT VECTOR FLOW ...
EFFICIENT NUMERICAL SCHEMES FOR GRADIENT VECTOR FLOW Djamal Boukerroui Universit´e de Technologie de Compi`egne, CNRS UMR 6599 Heudiasyc BP 20529 - 60205 Compigne Cedex, France. ABSTRACT
2. BACKGROUND MATERIALS
Since its publication 10 years ago, the GVF has been used and adapted to various models and problems. Its popularity is certainly due to its effectiveness. The main drawback of GVF and its generalization, is their expensive computation load and its consequence on the capture range. We propose and compare in this work different efficient numerical schemes to solve the GVF and its generalizations.
The GVF approach extends the gradient map using a diffusion process. GVF is defined as a vector field v(x ) : Ω → IRm that minimizes the following energy function [6]:
1. INTRODUCTION Active deformable contours are based on the minimization of a generally two terms energy function (internal and external) [1]. The internal encodes our prior about the shape of the contour and the external is data dependent and is designed to attract the solution to the desired boundary. Traditionally, the external energy is based only on edge information and is define as −f (·), where f (·) is an edge map of the image I(·), usually defined as |∇Gσ (·) ∗ I(·)|2 . Vectors of the edge map gradient ∇f are normal to the edge, and pointing towards it. ∇f is known as the external potential force field that pulls the contour towards the desired edges. These gradient vectors however have a large magnitude only in the immediate vicinity of the boundary. Several researchers pointed out this limitations and proposed alternatives solutions to increase the capture range of the external force (see eg. [2, 3]). Recent works have also the ability to capture concavities. For instance, the vector field convolution (VFC) [4] where the external force field is calculated using a simple convolution of the edge map with a user-defined vector field kernel. This simplicity combined with its robustness to noise makes of the VFC a very promising technique. The most popular and first solution is, however, the Gradient Vector Flow (GVF) [5] and its generalization, the GGVF [6]. Since its publication 10 years ago, the GVF has been used and adapted to various models and problems (eg. segmentation [5, 7, 8], tracking [9], skeletonisation [10]). The popularity of Xu’s GVF is certainly due to its effectiveness. A major drawback of GVF and its generalization GGVF, is their expensive computation load and its consequence on the capture range. Recently the multigrid method was applied with success to speed up the computation [11]. In this paper we propose and discuss other efficient alternatives.
Z
g(x)|∇v|2 + h(x)|v − ∇f |2 dx ,
E=
(1)
Ω
g(x) and h(x) are nonnegative functions on Ω. The basic idea
behind this variational formulation is that of making the solution smooth where there is no data, and to be as faithful as possible to the data where reliable measurements are possible. Using calculus of variations the solution of (1) can be obtained using the following diffusion equation [6] ∂t v = ∇ · g(x)∇v − h(x)(v − ∇f ) ,
(2)
with v(x, 0) = ∇f as an initial state and reflective boundary conditions. Depending on the choices of the two weighting functions g(x) and h(x), Xu and Prince proposed two vector flow fields. The GVF is defined for the special case where g(x) = µ and h(x) = |∇f |2 . µ is a positive constant parameter governing the amount of the smoothness term in (1). Equation (2) is then simplified to: ∂t v = µ∇2 v − h(|∇f |)(v − ∇f ) .
(3)
In order to avoid smoothing in the proximity of detected edges and increase the amount of regularization in areas where information is constant, the following choices are adopted by Xu and Prince [6] g(x) = e−
|∇f | K
, h(x) = 1 − g(|∇f |) .
(4)
They decided, however, to approximate the evolution equation (2) by ignoring the variations on g. This defines the GGVF as the solution of ∂t v = g(|∇f |)∇2 v − h(|∇f |)(v − ∇f ) .
(5)
We will denote by True Generalized Gradient Vector Flow (TGGVF), the solution of (2) using the above weighting functions. In the remaining of the paper, we limit ourselves to two dimensional domains and give details only for the
component u of the vector field v in the x direction. From (5), we get: ∂t u = g∇2 u − hu + h1 ,
with
h1 = h∂x f .
(6)
3. NUMERICAL IMPLEMENTATIONS We will denote by uki,j the discrete version of u(i∆x, j∆y, kτ ) where τ is the time step and ∆x and ∆y are the spatial grid sizes. It is natural and hence common in the image processing community, to use a uniform spatial grid (i.e ∆x = ∆y = h). Finite difference methods are generally used in this context in order to approximate the different derivatives which lead to a discrete set of equations. For instances, the forward difference scheme is used for the time derivative and central differences are used for spatial derivatives. a) Explicit scheme: The explicit scheme of the GGVF is given by [5]: k k k k k uk+1 i,j = ri,j (ui+1,j + ui,j+1 + ui−1,j + ui,j−1 − 4ui,j ) +
(1 − τ hi,j )uki,j + τ h1i,j ,
(11)
where I is the unit matrix, u ∈ IR2 is represented as a vector. The matrices Al corresponds to the derivatives along the lth coordinate axis, i.e the matrix-vector multiplication Al u is the 1 h(x)u which is discrete approximation of ∂xl g(x)∂xl u − m simply given by [12, 13] X gj + gi 1 (uj − ui ) − hi ui m 2h2
(12)
j∈N (i)
2h2 . maxi,j (h2 hi,j + 8gi,j )
(8)
Explicit schemes are very popular for solving partial differential equations, mainly because of their simplicity. This simplicity in counter part imposes very small and restrictive time steps to unsure stability, leading to slow convergence. b) Alternating Direction Explicit Scheme: The ADE is a 2 steps method and is known to be unconditionally stable. The two steps for the evolution equation (6) are given by: k+1 ri,j (uk+1 i−1,j +ui,j−1 ) = 1+2ri,j +τ hi,j k 1 ri,j (uk +u ) τ h i+1,j i,j+1 + 1+2ri,ji,j 1+2ri,j +τ hi,j +τ hi,j ,
(1−2ri,j )uk i,j 1+2ri,j +τ hi,j
n+2 ri,j (un+2 i+1,j +ui,j+1 ) un+2 − = i,j 1+2ri,j +τ hi,j k+1 1 ri,j (uk+1 +u ) τ h i−1,j i,j−1 + 1+2ri,ji,j 1+2ri,j +τ hi,j +τ hi,j .
(1−2ri,j )uk+1 i,j 1+2ri,j +τ hi,j
uk+1 − i,j
(I − τ Ax − τ Ay ) uk+1 = uk + τ h1 ,
(7)
where r = τ h−2 g. The above numerical scheme is called explicit as the solution at iteration k + 1 is explicitly given knowing the solution at the previous iteration k. The study of the stability of the above discrete schemes results to the following sufficient condition (see App. A)1 : τ ≤
numerical method suffers from anisotropy when very large time steps are used. In order to reduces such effect, left-right flip of the domain is recommended between iterations, We will denote by ADES this variant. c) Operating Splitting Schemes: Inspecting the evolution equation (2), we notice that the diffusion term takes exactly the form of the one in nonlinear diffusion filtering [12, eq. 1]. Splitting methods were used successfully to solve the later problem. In our case, there is an additional extra term “ h(x)(v − ∇f ) ”. Such techniques are know as Approximation Factorization methods [13]. For notational convenience, we rewrite the implicit scheme of (2) in a matrix vector notation:
where N (i) is the set of two neighbors of pixel i along the lth axis (boundary pixel have only one neighbor). Thus, the elements of Al are given by gj +gi (j ∈ N (i)), 2h2 P gi +gn h i (al )ij = − m − n∈N (i) 2h2 (13) (j = i), 0 (else). Barash et al. [13] analyzed different splitting operator to solve the diffusion equation that could be used to solve equation (11). Namely the Additive Operator Splitting (AOS), the Locally one dimensional (LOD) and the Additive Multiplicative Operator Splitting (AMOS) given respectively by:
uk+1 = +
2 1X (I − 2τ Al )−1 (uk + τ h1 ) , 2
(14)
l=1
(9) uk+1 =
2 Y
(I − τ Al )−1 (uk + τ h1 ) ,
(15)
1h (I − τ A1 )−1 (I − τ A2 )−1 + 2 i
(16)
l=1
+
(10)
Although it is not obvious, the above scheme is fully explicit. Indeed, in equation (9), the calculation proceeds in a raster scan and equation (10) proceeds in an inverse raster scan. This 1 Our stability condition is tighter than Xu and Prince’s τ ≤ ∆x∆y/4gmax [6, 5]. We believe that there was an error in Xu and Prince derivation. Fortunately the 2 conditions are equivalent for the GGVF and for the special case of the weighting functions given in (4) when h = 1;
uk+1 =
(I − τ A2 )−1 (I − τ A1 )−1 (uk + τ h1 ) .
The LOD is known to be the most efficient, but not symmetric. The AOS has the advantage of being be symmetric and the AMOS method is the symmetrized version of the LOD. All the above approximations however do not account for the Factorization Error (i.e difference between the true operator in (11) and the approximation). Tested on our problem, all
produced resonantly good results for small time steps. Analyzing the factorization error of the LOD approximation we get: EAF = τ 2 A1 A2 uk+1 . It is a second time order error and involves mixed forth-order derivatives of the solution at step k + 1. This explains one of our observations about the error behavior of the LOD scheme. The LOD method has tendency to miss estimated the gradient vector field near edge locations. As an alternative we propose the below approximation factorization, (I − τ A1 )u∗ =(I + τ A2 )(uk + τ h1 ) (I − τ A2 )u
k+1
∗
=u − τ A2 u
(17)
k
which has an EAF = τ 2 A1 A2 (uk − uk+1 ). In contrast to the LOD scheme, the error depends on the residual not on the solution. Hence, we expect it to decrease as the method converges. We will denoted this scheme by the AFI scheme. Given the above form of Al , all the above approximation schemes can be solved efficiently by the Thomas algorithm [12]. Notice that the above iterative solutions (14)-(17) solves the TGGVF and without any extra computational cost as the matrices Al are calculated only once. 4. RESULTS AND DISCUSSION We have presented several alternatives to the basic explicit scheme, that alow larger time steps. All schemes have a second order accuracy in space and are first order in time and are unconditionally stable. All methods methods have an O(N ) complexity. The computational time for one iteration of the ADE schemes is the same as for the Explicit scheme, and is 2 to 3 times slower for the ADES, LOD, AOS and AFI shemes. The AMOS runs about 5 times slower. This results are based on an average of 600 iterations runs and for 8 different image sizes. This imposes a minimum time step between 0.25 to 1.5 for the proposed methods to be faster than the explicit one. Lager time steps however may introduce a loss of accuracy. In section, we propose a way to quantify their accuracy when larger time steps are used. We have to keep also in mind our final objective for the use of GVF. Indeed, generally, in most uses of the GVF tech-
nique, the important information is the orientations of the vectors. The amplitude is rarely used in segmentation applications for example. To this end, we synthesized an edge map image composed of a perfect circle for which we could compute the theoretical optimal orientation map of the ideal vector field (see Fig. 1). Two criterions are then used to evaluate accuracy: The mean absolute angular error of the estimated vector field in comparison to the ground truth and the value of the energy function E given by equation (1). The results for different values of the time step τ are shown in Fig. 2. The first observation is the expected behavior of all methods for being very close to the explicit one for small τ values and then divert as τ increases. An important observation is that the proposed measure is very informative as its behavior is not always correlated with the energy. We notice that the curves of the mean absolute angular error are within 1◦ for τ values less than 10 except for the ADE, LOD and AOS schemes. The worst behavior are obtained by the LOD and the AOS methods (the highest error approaches 5◦ at time t = 300 for the LOD method with τ = 25). This confirms the observations made by Papandreou and Marogos when used in active contour models[14]. The results of the proposed AFI and ADES are very satisfying. The angular error is of the same order as the explicit method for t > 50 with τ = 10 and t > 100 for τ = 25; meaning after 4 to 5 iterations. Notice, as expected, the AFI method improves considerably the performances of approximation methods. It showed however a lower rate of convergence. The ADES method is probably the most appropriate. We have also implemented the the Full MultiGrid solver, proposed in [11]. The application of the FMG-GS(1,2,2) scheme gives 1.9◦ angular error. The ADES errors is lower than 1.9 after 4 iterations with τ = 15. The processing time is about the same for the two methods with our implementation. Our experiences suggest however to use of about 10 iterations. Further experimentations should be carried on. This experiments suggest, however, that the proposed scheme may be a promising alternative. 5. SUMMARY AND CONCLUSION Since its introduction by Xu and Prince, Gradient Vector Flow has become very popular and is extensively used in contour
Fig. 2. (left) evolution of the GGVF energy function (1) versus Fig. 1. (left) edge map image; (right) theoretical optimal orientation map of the ideal vector field.
time for different schemes; (right) Evolution of the mean absolute angular error in degree versus time. The small images are a zoom around t = 300.
based segmentation as an external energy force. Probably the main drawback of such a force is its high computation burden. Efficient schemes for the minimization of the GVF and its generalization are proposed and analyzed. We have shown that the ADES scheme may be a good alternative to the Multigrid technique. The proposed AFI scheme also showed its superiority on the others approximation factorization methods. Although included, we did not discuss much the AMOS method as its behavior may be expected, and also because it is the most computationally expensive. Finally we may emphasis the fact that the Alternative Direction Explicit schemes are not well know by the image processing community. They may be also be of interest to solve other PDE problems for image processing. 6. REFERENCES [1] M Kass, A Witkin, and D Terzopoulos, “Snake : Active contour model,” Int J. Computer Vision, vol. 1, pp. 321– 331, 1987. [2] L.D Cohen and I Cohen, “Finite element method for active contour models and balloons for 2D and 3D images,” IEEE PAMI, vol. 15, no. 11, pp. 1131–1147, 1993. [3] B Leroy, I Herlin, and L.D Cohen, “Multi-resolution algorithm for active contour models,” In 12th Int. Conf. Anal. Opt. System, vol. 15, no. 11, pp. 58–65, 1996. [4] Bing Li and Scott T Acton, “Active contour external force using vector field convolution for image segmentation,” IEEE Image Process., vol. 16, no. 8, pp. 2096– 2106, 2007. [5] Chenyang Xu and J L Price, “Snakes, shapes, and gradient vector flow,” IEEE Image Process., vol. 7, no. 3, pp. 359–369, Mar. 1998. [6] Chenyang Xu and J L Price, “Generalized gradient vector flow external forces for active contours,” Signal Processing, vol. 71, pp. 131–139, 1998.
[10] M. Sabry Hassouna and Aly A. Farag, “On the extraction of curve skeletons using gradient vector flow,” in ICCV, Brazil, Oct. 14-20, 2007, pp. 1–8. [11] X. Han, C. Xu, and J.L. Prince, “Fast numerical scheme for gradient vector flow computation using a multigrid method,” IET Image Processing, vol. 1, no. 1, pp. 48 – 55, Mar. 2007. [12] Joachim Weickert, Bart M. ter Haar Romeny, and Max A. Viergever, “Efficient reliable schemes for nonlinear diffusion filtering,” IEEE Image Process., vol. 7, no. 3, pp. 398–410, 1998. [13] D. Barash, T. Schlick, M. Israeli, and R. Kimmel, “Multiplicative operator splitting in nonlinear diffusion: from spatial splitting to multiple timesteps,” Journal of Mathematical Imaging and Vision, vol. 19, no. 1, pp. 33 – 48, Mar. 2003. [14] Papandreou G and Maragos P., “Multigrid geometric active contour models.,” IEEE Trans Image Process., vol. 16, no. 1, pp. 229–40, Jan. 2007. [15] G Aubert and P. Kornprobst, Mathematical Problems in Image Processing, vol. 147 of Applied Mathematical Sciences, Springer Verlag, Nov. 2001. A. STABILITY OF A DIFFERENCE SCHEME The stability of a difference can be studied by means of the discrete von Neumann criterion for stability (see eg. [15]). This approach considers a discrete Fourier mode for the problem. In our case given by uki,j = ξ k exp{ıπh(im+jn)}, where, the subscript on the ξ term is a multiplicative exponent. Inserting this general Fourier mode in the difference scheme, a necessary condition for stability is obtained by restricting h and τ so that |ξ| ≤ 1. Doing so for the explicit scheme given in equation (7) we obtain: ξ = 1 − τ hi,j − 4ri,j + 2ri,j (cos(πmh) + cos(πnh)) .
Let us denote by γ = cos(πmh) + cos(πnh) and replace ri,j by τ h−2 gi,j , and after some algebraic manipulations we get the following condition:
[7] N. Paragios and O. Mellina-Gottardoand V. Ramesh, “Gradient vector flow fast geometric active contours,” IEEE PAMI, vol. 26, no. 3, pp. 402–407, Mar. 2004. [8] O. Colliot, O. Camara, and I. Bloch, “Integration of fuzzy spatial relations in deformable models– Application to brain MRI segmentation,” Pattern Recogn., vol. 39, no. 8, pp. 1401–1414, Aug. 2006. [9] N Ray and S T Acton, “Motion gradient vector flow: An external force for tracking rolling leukocytes with shape and size constrained active contours,” IEEE Med. Imag., vol. 23, no. 12, pp. 1466–1478, 2004.
0≤τ
2gi,j hi,j + 2 (2 − γ) ≤ 2 . h
Since γ is the sum of two cosines, then 0 ≤ (2 − γ) ≤ 4. This leads to 0≤τ
2gi,j 8gi,j hi,j + 2 (2 − γ) ≤ τ hi,j + 2 ≤2 . h h
The stability condition for the explicit scheme is then given by: τ ≤
2h2 maxi,j (h2 hi,j
+ 8gi,j )
(A.1)
2. BACKGROUND MATERIALS
Since its publication 10 years ago, the GVF has been used and adapted to various models and problems. Its popularity is certainly due to its effectiveness. The main drawback of GVF and its generalization, is their expensive computation load and its consequence on the capture range. We propose and compare in this work different efficient numerical schemes to solve the GVF and its generalizations.
The GVF approach extends the gradient map using a diffusion process. GVF is defined as a vector field v(x ) : Ω → IRm that minimizes the following energy function [6]:
1. INTRODUCTION Active deformable contours are based on the minimization of a generally two terms energy function (internal and external) [1]. The internal encodes our prior about the shape of the contour and the external is data dependent and is designed to attract the solution to the desired boundary. Traditionally, the external energy is based only on edge information and is define as −f (·), where f (·) is an edge map of the image I(·), usually defined as |∇Gσ (·) ∗ I(·)|2 . Vectors of the edge map gradient ∇f are normal to the edge, and pointing towards it. ∇f is known as the external potential force field that pulls the contour towards the desired edges. These gradient vectors however have a large magnitude only in the immediate vicinity of the boundary. Several researchers pointed out this limitations and proposed alternatives solutions to increase the capture range of the external force (see eg. [2, 3]). Recent works have also the ability to capture concavities. For instance, the vector field convolution (VFC) [4] where the external force field is calculated using a simple convolution of the edge map with a user-defined vector field kernel. This simplicity combined with its robustness to noise makes of the VFC a very promising technique. The most popular and first solution is, however, the Gradient Vector Flow (GVF) [5] and its generalization, the GGVF [6]. Since its publication 10 years ago, the GVF has been used and adapted to various models and problems (eg. segmentation [5, 7, 8], tracking [9], skeletonisation [10]). The popularity of Xu’s GVF is certainly due to its effectiveness. A major drawback of GVF and its generalization GGVF, is their expensive computation load and its consequence on the capture range. Recently the multigrid method was applied with success to speed up the computation [11]. In this paper we propose and discuss other efficient alternatives.
Z
g(x)|∇v|2 + h(x)|v − ∇f |2 dx ,
E=
(1)
Ω
g(x) and h(x) are nonnegative functions on Ω. The basic idea
behind this variational formulation is that of making the solution smooth where there is no data, and to be as faithful as possible to the data where reliable measurements are possible. Using calculus of variations the solution of (1) can be obtained using the following diffusion equation [6] ∂t v = ∇ · g(x)∇v − h(x)(v − ∇f ) ,
(2)
with v(x, 0) = ∇f as an initial state and reflective boundary conditions. Depending on the choices of the two weighting functions g(x) and h(x), Xu and Prince proposed two vector flow fields. The GVF is defined for the special case where g(x) = µ and h(x) = |∇f |2 . µ is a positive constant parameter governing the amount of the smoothness term in (1). Equation (2) is then simplified to: ∂t v = µ∇2 v − h(|∇f |)(v − ∇f ) .
(3)
In order to avoid smoothing in the proximity of detected edges and increase the amount of regularization in areas where information is constant, the following choices are adopted by Xu and Prince [6] g(x) = e−
|∇f | K
, h(x) = 1 − g(|∇f |) .
(4)
They decided, however, to approximate the evolution equation (2) by ignoring the variations on g. This defines the GGVF as the solution of ∂t v = g(|∇f |)∇2 v − h(|∇f |)(v − ∇f ) .
(5)
We will denote by True Generalized Gradient Vector Flow (TGGVF), the solution of (2) using the above weighting functions. In the remaining of the paper, we limit ourselves to two dimensional domains and give details only for the
component u of the vector field v in the x direction. From (5), we get: ∂t u = g∇2 u − hu + h1 ,
with
h1 = h∂x f .
(6)
3. NUMERICAL IMPLEMENTATIONS We will denote by uki,j the discrete version of u(i∆x, j∆y, kτ ) where τ is the time step and ∆x and ∆y are the spatial grid sizes. It is natural and hence common in the image processing community, to use a uniform spatial grid (i.e ∆x = ∆y = h). Finite difference methods are generally used in this context in order to approximate the different derivatives which lead to a discrete set of equations. For instances, the forward difference scheme is used for the time derivative and central differences are used for spatial derivatives. a) Explicit scheme: The explicit scheme of the GGVF is given by [5]: k k k k k uk+1 i,j = ri,j (ui+1,j + ui,j+1 + ui−1,j + ui,j−1 − 4ui,j ) +
(1 − τ hi,j )uki,j + τ h1i,j ,
(11)
where I is the unit matrix, u ∈ IR2 is represented as a vector. The matrices Al corresponds to the derivatives along the lth coordinate axis, i.e the matrix-vector multiplication Al u is the 1 h(x)u which is discrete approximation of ∂xl g(x)∂xl u − m simply given by [12, 13] X gj + gi 1 (uj − ui ) − hi ui m 2h2
(12)
j∈N (i)
2h2 . maxi,j (h2 hi,j + 8gi,j )
(8)
Explicit schemes are very popular for solving partial differential equations, mainly because of their simplicity. This simplicity in counter part imposes very small and restrictive time steps to unsure stability, leading to slow convergence. b) Alternating Direction Explicit Scheme: The ADE is a 2 steps method and is known to be unconditionally stable. The two steps for the evolution equation (6) are given by: k+1 ri,j (uk+1 i−1,j +ui,j−1 ) = 1+2ri,j +τ hi,j k 1 ri,j (uk +u ) τ h i+1,j i,j+1 + 1+2ri,ji,j 1+2ri,j +τ hi,j +τ hi,j ,
(1−2ri,j )uk i,j 1+2ri,j +τ hi,j
n+2 ri,j (un+2 i+1,j +ui,j+1 ) un+2 − = i,j 1+2ri,j +τ hi,j k+1 1 ri,j (uk+1 +u ) τ h i−1,j i,j−1 + 1+2ri,ji,j 1+2ri,j +τ hi,j +τ hi,j .
(1−2ri,j )uk+1 i,j 1+2ri,j +τ hi,j
uk+1 − i,j
(I − τ Ax − τ Ay ) uk+1 = uk + τ h1 ,
(7)
where r = τ h−2 g. The above numerical scheme is called explicit as the solution at iteration k + 1 is explicitly given knowing the solution at the previous iteration k. The study of the stability of the above discrete schemes results to the following sufficient condition (see App. A)1 : τ ≤
numerical method suffers from anisotropy when very large time steps are used. In order to reduces such effect, left-right flip of the domain is recommended between iterations, We will denote by ADES this variant. c) Operating Splitting Schemes: Inspecting the evolution equation (2), we notice that the diffusion term takes exactly the form of the one in nonlinear diffusion filtering [12, eq. 1]. Splitting methods were used successfully to solve the later problem. In our case, there is an additional extra term “ h(x)(v − ∇f ) ”. Such techniques are know as Approximation Factorization methods [13]. For notational convenience, we rewrite the implicit scheme of (2) in a matrix vector notation:
where N (i) is the set of two neighbors of pixel i along the lth axis (boundary pixel have only one neighbor). Thus, the elements of Al are given by gj +gi (j ∈ N (i)), 2h2 P gi +gn h i (al )ij = − m − n∈N (i) 2h2 (13) (j = i), 0 (else). Barash et al. [13] analyzed different splitting operator to solve the diffusion equation that could be used to solve equation (11). Namely the Additive Operator Splitting (AOS), the Locally one dimensional (LOD) and the Additive Multiplicative Operator Splitting (AMOS) given respectively by:
uk+1 = +
2 1X (I − 2τ Al )−1 (uk + τ h1 ) , 2
(14)
l=1
(9) uk+1 =
2 Y
(I − τ Al )−1 (uk + τ h1 ) ,
(15)
1h (I − τ A1 )−1 (I − τ A2 )−1 + 2 i
(16)
l=1
+
(10)
Although it is not obvious, the above scheme is fully explicit. Indeed, in equation (9), the calculation proceeds in a raster scan and equation (10) proceeds in an inverse raster scan. This 1 Our stability condition is tighter than Xu and Prince’s τ ≤ ∆x∆y/4gmax [6, 5]. We believe that there was an error in Xu and Prince derivation. Fortunately the 2 conditions are equivalent for the GGVF and for the special case of the weighting functions given in (4) when h = 1;
uk+1 =
(I − τ A2 )−1 (I − τ A1 )−1 (uk + τ h1 ) .
The LOD is known to be the most efficient, but not symmetric. The AOS has the advantage of being be symmetric and the AMOS method is the symmetrized version of the LOD. All the above approximations however do not account for the Factorization Error (i.e difference between the true operator in (11) and the approximation). Tested on our problem, all
produced resonantly good results for small time steps. Analyzing the factorization error of the LOD approximation we get: EAF = τ 2 A1 A2 uk+1 . It is a second time order error and involves mixed forth-order derivatives of the solution at step k + 1. This explains one of our observations about the error behavior of the LOD scheme. The LOD method has tendency to miss estimated the gradient vector field near edge locations. As an alternative we propose the below approximation factorization, (I − τ A1 )u∗ =(I + τ A2 )(uk + τ h1 ) (I − τ A2 )u
k+1
∗
=u − τ A2 u
(17)
k
which has an EAF = τ 2 A1 A2 (uk − uk+1 ). In contrast to the LOD scheme, the error depends on the residual not on the solution. Hence, we expect it to decrease as the method converges. We will denoted this scheme by the AFI scheme. Given the above form of Al , all the above approximation schemes can be solved efficiently by the Thomas algorithm [12]. Notice that the above iterative solutions (14)-(17) solves the TGGVF and without any extra computational cost as the matrices Al are calculated only once. 4. RESULTS AND DISCUSSION We have presented several alternatives to the basic explicit scheme, that alow larger time steps. All schemes have a second order accuracy in space and are first order in time and are unconditionally stable. All methods methods have an O(N ) complexity. The computational time for one iteration of the ADE schemes is the same as for the Explicit scheme, and is 2 to 3 times slower for the ADES, LOD, AOS and AFI shemes. The AMOS runs about 5 times slower. This results are based on an average of 600 iterations runs and for 8 different image sizes. This imposes a minimum time step between 0.25 to 1.5 for the proposed methods to be faster than the explicit one. Lager time steps however may introduce a loss of accuracy. In section, we propose a way to quantify their accuracy when larger time steps are used. We have to keep also in mind our final objective for the use of GVF. Indeed, generally, in most uses of the GVF tech-
nique, the important information is the orientations of the vectors. The amplitude is rarely used in segmentation applications for example. To this end, we synthesized an edge map image composed of a perfect circle for which we could compute the theoretical optimal orientation map of the ideal vector field (see Fig. 1). Two criterions are then used to evaluate accuracy: The mean absolute angular error of the estimated vector field in comparison to the ground truth and the value of the energy function E given by equation (1). The results for different values of the time step τ are shown in Fig. 2. The first observation is the expected behavior of all methods for being very close to the explicit one for small τ values and then divert as τ increases. An important observation is that the proposed measure is very informative as its behavior is not always correlated with the energy. We notice that the curves of the mean absolute angular error are within 1◦ for τ values less than 10 except for the ADE, LOD and AOS schemes. The worst behavior are obtained by the LOD and the AOS methods (the highest error approaches 5◦ at time t = 300 for the LOD method with τ = 25). This confirms the observations made by Papandreou and Marogos when used in active contour models[14]. The results of the proposed AFI and ADES are very satisfying. The angular error is of the same order as the explicit method for t > 50 with τ = 10 and t > 100 for τ = 25; meaning after 4 to 5 iterations. Notice, as expected, the AFI method improves considerably the performances of approximation methods. It showed however a lower rate of convergence. The ADES method is probably the most appropriate. We have also implemented the the Full MultiGrid solver, proposed in [11]. The application of the FMG-GS(1,2,2) scheme gives 1.9◦ angular error. The ADES errors is lower than 1.9 after 4 iterations with τ = 15. The processing time is about the same for the two methods with our implementation. Our experiences suggest however to use of about 10 iterations. Further experimentations should be carried on. This experiments suggest, however, that the proposed scheme may be a promising alternative. 5. SUMMARY AND CONCLUSION Since its introduction by Xu and Prince, Gradient Vector Flow has become very popular and is extensively used in contour
Fig. 2. (left) evolution of the GGVF energy function (1) versus Fig. 1. (left) edge map image; (right) theoretical optimal orientation map of the ideal vector field.
time for different schemes; (right) Evolution of the mean absolute angular error in degree versus time. The small images are a zoom around t = 300.
based segmentation as an external energy force. Probably the main drawback of such a force is its high computation burden. Efficient schemes for the minimization of the GVF and its generalization are proposed and analyzed. We have shown that the ADES scheme may be a good alternative to the Multigrid technique. The proposed AFI scheme also showed its superiority on the others approximation factorization methods. Although included, we did not discuss much the AMOS method as its behavior may be expected, and also because it is the most computationally expensive. Finally we may emphasis the fact that the Alternative Direction Explicit schemes are not well know by the image processing community. They may be also be of interest to solve other PDE problems for image processing. 6. REFERENCES [1] M Kass, A Witkin, and D Terzopoulos, “Snake : Active contour model,” Int J. Computer Vision, vol. 1, pp. 321– 331, 1987. [2] L.D Cohen and I Cohen, “Finite element method for active contour models and balloons for 2D and 3D images,” IEEE PAMI, vol. 15, no. 11, pp. 1131–1147, 1993. [3] B Leroy, I Herlin, and L.D Cohen, “Multi-resolution algorithm for active contour models,” In 12th Int. Conf. Anal. Opt. System, vol. 15, no. 11, pp. 58–65, 1996. [4] Bing Li and Scott T Acton, “Active contour external force using vector field convolution for image segmentation,” IEEE Image Process., vol. 16, no. 8, pp. 2096– 2106, 2007. [5] Chenyang Xu and J L Price, “Snakes, shapes, and gradient vector flow,” IEEE Image Process., vol. 7, no. 3, pp. 359–369, Mar. 1998. [6] Chenyang Xu and J L Price, “Generalized gradient vector flow external forces for active contours,” Signal Processing, vol. 71, pp. 131–139, 1998.
[10] M. Sabry Hassouna and Aly A. Farag, “On the extraction of curve skeletons using gradient vector flow,” in ICCV, Brazil, Oct. 14-20, 2007, pp. 1–8. [11] X. Han, C. Xu, and J.L. Prince, “Fast numerical scheme for gradient vector flow computation using a multigrid method,” IET Image Processing, vol. 1, no. 1, pp. 48 – 55, Mar. 2007. [12] Joachim Weickert, Bart M. ter Haar Romeny, and Max A. Viergever, “Efficient reliable schemes for nonlinear diffusion filtering,” IEEE Image Process., vol. 7, no. 3, pp. 398–410, 1998. [13] D. Barash, T. Schlick, M. Israeli, and R. Kimmel, “Multiplicative operator splitting in nonlinear diffusion: from spatial splitting to multiple timesteps,” Journal of Mathematical Imaging and Vision, vol. 19, no. 1, pp. 33 – 48, Mar. 2003. [14] Papandreou G and Maragos P., “Multigrid geometric active contour models.,” IEEE Trans Image Process., vol. 16, no. 1, pp. 229–40, Jan. 2007. [15] G Aubert and P. Kornprobst, Mathematical Problems in Image Processing, vol. 147 of Applied Mathematical Sciences, Springer Verlag, Nov. 2001. A. STABILITY OF A DIFFERENCE SCHEME The stability of a difference can be studied by means of the discrete von Neumann criterion for stability (see eg. [15]). This approach considers a discrete Fourier mode for the problem. In our case given by uki,j = ξ k exp{ıπh(im+jn)}, where, the subscript on the ξ term is a multiplicative exponent. Inserting this general Fourier mode in the difference scheme, a necessary condition for stability is obtained by restricting h and τ so that |ξ| ≤ 1. Doing so for the explicit scheme given in equation (7) we obtain: ξ = 1 − τ hi,j − 4ri,j + 2ri,j (cos(πmh) + cos(πnh)) .
Let us denote by γ = cos(πmh) + cos(πnh) and replace ri,j by τ h−2 gi,j , and after some algebraic manipulations we get the following condition:
[7] N. Paragios and O. Mellina-Gottardoand V. Ramesh, “Gradient vector flow fast geometric active contours,” IEEE PAMI, vol. 26, no. 3, pp. 402–407, Mar. 2004. [8] O. Colliot, O. Camara, and I. Bloch, “Integration of fuzzy spatial relations in deformable models– Application to brain MRI segmentation,” Pattern Recogn., vol. 39, no. 8, pp. 1401–1414, Aug. 2006. [9] N Ray and S T Acton, “Motion gradient vector flow: An external force for tracking rolling leukocytes with shape and size constrained active contours,” IEEE Med. Imag., vol. 23, no. 12, pp. 1466–1478, 2004.
0≤τ
2gi,j hi,j + 2 (2 − γ) ≤ 2 . h
Since γ is the sum of two cosines, then 0 ≤ (2 − γ) ≤ 4. This leads to 0≤τ
2gi,j 8gi,j hi,j + 2 (2 − γ) ≤ τ hi,j + 2 ≤2 . h h
The stability condition for the explicit scheme is then given by: τ ≤
2h2 maxi,j (h2 hi,j
+ 8gi,j )
(A.1)
