Multigrid methods and automatic segmentation: an application to CT images of the liver

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Multigrid methods and automatic segmentation: an application to CT images of the liver J. F. Garamendi1, N. Malpica1, E. Schiavi2 1Dpto. Ingenier´ıa Telem´atica y Tecnolog´ıa Electr´onica, Univ. Rey Juan Carlos de Madrid {juanfrancisco.garamendi, norberto.malpica}@urjc.es. 2Dpto. de Matem´atica Aplicada, Univ. Rey Juan Carlos de Madrid [email protected]. Palabras clave: Quasilinear elliptic equations, Total variation minimization, Multigrid methods, binary segmentation, medical images, liver surgery. Resumen We consider a segmentation problem which arises in medical imaging and liver surgery. The model problem is based on an active contour without edges technique formulated in a level set dictionary. Previous work indicates that a feasible solution can be obtained solving the gradient descent equation associated to the original minimization problem but the convergence of the algorithm is too slow for practical clinical purposes. Here, we study the implementation of multigrid methods to the elliptic problem and the numerical results are compared with the parabolic approach. 1. Introduction Planning liver surgery requires accurate volumetric measures of the organ to be operated on. After acquiring images of the adequate modality for a certain organ, the main step of the process is the segmentation of the organ, which is generally carried out by manually using a software platform. Manual delineation of the liver is time-consuming and can lack repeatability among users. This leads to the need of developing automatic segmentation techniques. The automatic segmentation of the liver is quite challenging due to several structural reasons, among with, the fact that the Hounsfield units corresponding to the liver are the same as those of neighbouring organs, so the use of simple gray level segmentation methods give inaccurate results. Also, the shape and size of the liver can vary a lot between patients, making it difficult to impose a priori conditions on pure statistical based models. Recently, we proposed a PDE based method for the automatic 1 J. F. Garamendi, N. Malpica, E. Schiavi segmentation of the liver, which needs no a priori shape or size information and stops automatically [2]. It is based on Total Variation regularization, a very useful technique in image restoration [11]. The associated PDE is solved following a gradient descent method [2]. Results of the algorithm have been compared to manual segmentations by an expert with good correlation. The convergence of the gradient descent implementation, however, is too slow for pratical clinical use. In this work we study the possibility to accelerate the algorithm using multigrid techniques. Multigrid methods are generally accepted as being the fastest numerical methods for the solution of elliptic partial differential equations [4] but their perfomance and accuracy depend on the structural difficulties of the underlying PDE. This paper is organized as follows: in section 2 the model is described in detail and in section 3 the algorithm is shown. Section 4 shows some results on clinical datasets. Finally, we discuss our conclusions and some open problems which deserve a throughout study (section 5). 2. Model equations Efficient resolution of segmentation models is a basic problem in image processing and computer vision. We consider here the active contour method proposed by Chan and Vese, [7], which can be deduced from the piecewise Mumford-Shah model for image segmentation [5]. In its level set formulation this leads to consider a quasilinear elliptic equation for the Total Variation (TV) differential operator (which corresponds to the p−laplacian operator for p= 1). Basically the active contour is the set of discontinuities of the Mumford-Shah problem (the minimal partition problem) and the resulting functional defines gradient descent evolution equations to deform the active contour which provide the final segmentation. Considering that the original problem is an elliptic one, this amounts to solve parabolic equations until stabilization and it is typically slow, not feasible for medical image segmentation [2]. Moreover, these variational problems have multiple (local and global) minima and this makes the initial conditions for gradient descent critical. Some effort has been dedicated to speed up the algorithms involved in the resolution of TV operator in this context [8] and we explore here their suggestions. Reformulating the original, not convex minimization problem associated to the two-phase (binary) Chan-Vese model (CV) as a convex one we solve the Chan-Essedouglu-Nikolova model (CEN) in a multigrid framework. Chan-Vese model Let Ω ⊂IR2be an open, bounded domain (usually a rectangle) where (x, y)∈Ω denotes pixel location and I0(x, y) is a function representing the intensity image values. The Chan-Vese model for binary segmentation is based on the minimization of an energy functional expressed in terms of a level set dictionary which we briefly resume. Let ω⊂Ω be an open (eventually not connected) positive measured sub-region of the original domain. If the curve Crepresents the boundary of such a segmentation ωthen, in the level set formulation, the (free) boundary Cis the zero level set of a Lipschitz function φ: Ω →IR, that is: C={(x, y)∈Ω : φ(x, y) = 0},C=∂ω where ω={(x, y)∈Ω : φ(x, y)>0}, Ω\ω={(x, y)∈Ω : φ(x, y)<0}. Let H(φ) be the Heaviside function defined by H(φ) = 1 if φ≥0 and H(φ) = 0 otherwise. Notice that the term H(φ) parametrize 2 Multigrid methods for liver segmentation the level set in the class of binary functions. The level set function φ(the active contour) can be characterized as a minimum of the following energy functional: Jcv(φ) = µZΩ |∇H(φ)|dxdy +λin ZΩ H(φ)eindxdy +λout ZΩ (1 −H(φ))eoutdxdy (1) where µ,λin and λout are parameters which can be considered as weight factors which control the trade off between smoothness (µ) and the fidelity data terms (λin,λout). The functions ein,eout are defined as [2] ein(φ) =     ZD I0dxdy |D|−cin(φ)    2 , eout(φ) =     ZD I0dsdr |D|−cout(φ)    2 (2) where D=Dx,y ⊂Ω is a neighbour centered at pixel (x, y)⊂Ω. The local mean term in (2) is a slight modification of the functional proposed in [8] which we proposed in [2]. A two step algorithm which makes ein,eout,cin,cout explicit is implemented. The (constants) cin and cout are the mean value inside (cin) and outside (cout) the binary segmentation. Following variational calculus, the minimum of (1) corresponds to a solution of the Euler-Lagrange equation: 0 = δ(φ)µ∇ · ∇φ |∇φ|−λinein(φ) + λouteout(φ), a.e. (x, y)∈Ω (3) where δ(φ) is a non-compactly supported, smooth approximation of the Dirac delta function δ(φ) located at C[7]. The equation is complemented with (no flux) homogeneous Neumann boundary conditions. Chan-Esedo¯ glu-Nikolova model This model has been proposed in [8] and convexifies minimization of (1) in the class of binary functions (characteristic functions of sets which are a non-convex collection). The key idea is to extend the minimization to all functions and to transform the minimizers of the original (non-convex minimization problem) into minimizers of (1) by thresholding. As a result, equation (3) has the same solutions as: 0 = µ∇ · ∇φ |∇φ|−λinein +λouteout, a.e. (x, y)∈Ω which is the Euler-Lagrange of the (linear) energy functional: J(φ) = µZΩ |∇φ|dxdy +λin ZΩ φ ein(φ)dxdy −λout ZΩ φ eout(φ)dxdy which has to be constrained [8] to the level sets 0 ≤φ≤1 to fix the non-uniqueness of the representation of the level sets. After some observations it is claimed that for any given fixed cin,cout ∈IR (hence ein,eout by (2)) a global minimizer for the energy (1) can be found solving the convex, unconstrained minimization problem: Jcen(φ) = µZΩ |∇φ|dxdy +αZΩ ν(φ) + λin ZΩ φ eindxdy −λout ZΩ φ eoutdxdy (4) 3 J. F. Garamendi, N. Malpica, E. Schiavi where ν(ξ) :=max{0,2|ξ−1/2| − 1}is a penalty term and αhas to be properly chosen [8]. The (global) minima of the energy (4) correspond to a solution of the Euler-Lagrange equation: 0 = µ∇ · ∇φ |∇φ|−λinein +λouteout −αν0(φ), a.e. (x, y)∈Ω (5) which we discretize in the following section. 2.1. Numerical Implementation The fast and accurate numerical resolution of the TV operator is still a challenging problem and it has been typically solved by finite difference discrete formulations embedded in time-marching schemes. Interesting contributions are due, among many others, to Vogel and Oman, [12], with the so called lagged diffusivity method which corresponds to a fixed-point iteration, Chan et al. with their primal-dual method [6] and Oman for a multigrid approach [10] in the context of denoising and deblurring. In this section we describe and compare the algorithms we used to solve (5). A number of additional difficulties arise when considering (5) which are associated to the multivalued character of the subdifferential of the penalty term. We shall discuss this in the last section. Parabolic equation Instead of solving the elliptic equation (3) directly, we solve the associated parabolic equation: ∂φ ∂t =δ(φ)µ∇ · ∇φ |∇φ|−λinein +λouteout(6) This introduces an artificial time which diminishes the difficulties inherent to the TV elliptic operator. A simple L1-norm bound is selected as a stopping rule between consecutives iterations. The derivatives are implemented with a finite difference scheme, where D± ddenote the forward (backward) difference in the direction d. The parabolic equation (6) is solved by an explicit gradient descent method: φn+1 −φn ∆t=δ(φn) µD− x  D+ xφn q(D+ xφn)2+ (D+ yφn)2+1  + +µD− y  D+ yφn q(D+ xφn)2+ (D+ yφn)2+1  −λinein +λouteout  where ein and eout are computed using (2) and 0 < 11 is a small regularizing parameter avoiding the singularity due to the nondifferentiability of the TV functional when ∇φ= 0. This explicit scheme linearizes the highly nonlinear oscillating diffusion term in (6). 4 Multigrid methods for liver segmentation Elliptic equation Equation (5) is discretized following a finite difference scheme: 0 = µD− x  D+ xφ q(D+ xφ)2+ (D+ yφ)2+1  +µD− y  D+ yφ q(D+ xφ)2+ (D+ yφ)2+1   −λinein +λouteout −αν0 2(φ) (7) where ν0 2is a smooth approximation of ν0, which is a bounded maximal monotone graph. System of equations obtained by finite differences discretization is solved following a multigrid method. Multigrid methods for numerical analysis are algorithms for solving differential equations using a hierarchy of discretizations between coarser and finer grids. The key idea is that low frequencies in finer grids appear as high frequencies in coarser grids. They can be traced back to the pioneering work of A. Brandt. The typical application for multigrid methods is in the numerical solution of elliptic partial differential equations. Notably, they do not depend on the kind of discretization (finite difference or finite elements). Nevertheless it is known that highly oscillating coefficients in quasilinear equations can deteriorate the performance of this strategy [13]. A great deal of care is needed in order to assess the basic steps of the algorithm. Iterative solvers, such as Gauss-Seidel, have fast smoothing properties but, after few iterations the convergence slows down significantly. Multigrid algorithms take this idea and combine it with a hierarchycal model of equations systems that come from differents levels of detail in the problem discretization. There are several kinds of multigrid methods: multilevel methods, bidirectional multigrid and full multigrid. The most simple (multilevel methods) use an aproximated solutions, calculated from some level of discretization, as a seed for the iterative solver. However, the main idea of more advanced multigrid methods is to first approximate a solution in a level of discretization, then calculate the error of the solution from a coarser level and to correct the approximate solution with the computed error. In the next section we describe our algorithm and compare the results obtained with multilevel, bidirectional and full multigrid. 3. Algorithm Starting with the finite difference discretization scheme (5) we use the frozen coefficient method of Savage and Chen [14] to linearize the elliptic operator. Also a fixed point iteration is implemented to deal with the term ν0 2. The resulting system of equations can be written in form A(x) = b. We briefly describe the multigrid methods we used for comparison. Multilevel Method: The equation system A(x) = bobtained at the finest grid is solved by an iterative method using a seed calculated at a coarser level. This approach is implemented iteratively through several levels starting from the coarsest one. Bidirectional multigrid (FAS method): Let ˜xbe an approximated solution where x= ˜x+eis the exact solution of A(x) = band eis the error which can be estimated from the residual equation r=b−A(˜x) solving the residual system A(˜x+e)−A(˜x) = r. 5 J. F. Garamendi, N. Malpica, E. Schiavi This has the same complexity as the original system A(x) = b, but we can solve it at a coarser discretisation level, where the system is smaller and easier to solve. Once we have estimated eat the coarser level, we calculate ˜einterpolating efrom the coarser level to the finest level and later we do the correction ˜ ˜x= ˜x+ ˜e. The biridectional multigrid method is based on this scheme of prediction-correction using several levels of discretisation. This approach is also called V-cycle. Full multigrid: Full multigrids methods combine the multilevel approach and the bidirectional implementations. Basically, it is a multilevel method where a bidirectional multigrid is implemented at each level. The global procedure for liver segmentation is implemented using a multistep approach, consisting of three tightly interleaved steps: preprocessing step, segmentation step, and postprocessing step. The preprocessing and postprocessing step are carried out as in [2]. The segmentation step is as follows: 1. Initialize φ0with a segmentation by Otsu algorithm. 2. Compute cout(φn). 3. Solve the system (7) using a multigrid method, obtaining φn+1. 4. Reinitialize φto a binary function. 5. Go to 2 until |φn−φn−1|< tolerance 4. Results To test the three multigrid implementations, we used a set of CT images of the liver, acquired at Fundaci´on Hospital de Alcorc´on in Madrid. All images are 512x512 pixels. All images were manually segmented by an expert radiologist and automatically using the three versions of the multigrid algorithms with different settings. Results of the segmentation are shown in figure 1, where we considered the parabolic equation by gradient descent, parabolic equation by gradient descent in a multilevel scheme, the elliptic equation in a V-Cycle and the elliptic equation in the full multigrid scheme. The black contour in the top row is the manual segmentation performed by the radiologist and the white (active) contour is the automatic segmentation. Tables 1-4 show, for every method, the symmetric difference between the manual and automatic segmentations and the time taken by the method. Our preliminary results indicate that in terms of time consumption the best algorithm is the V-cycle with five levels which also minimizes the symmetric difference between the automatic and the manual segmentation. On the contrary, the full multigrid performs worst as we increase the number of levels. This, in a different problem but with the same elliptic operator has been observed also in [13]. Symmetric difference CPU Time 3.230599 321.156 Tabla 1: Results of gradient descent method 6 Multigrid methods for liver segmentation Figura 1: Top row shows manual segmentation (in black) and automatic segmentations (in white), from left to right, parabolic (gradient descent), elliptic with multilevel (2 levels), elliptic with V Cycle (5 levels) and elliptic with full multigrid (3 levels). At the bottom row, the symmetric differences between manual segmentation and automatic segmentation are shown in black. Coarse grid size # Levels Symmetric difference CPU Time 25x255 3.227273 410.688000 26x264 3.229491 410.297000 27x273 3.225056 401.421000 28x282 3.219513 390.516000 Tabla 2: Results of multilevel method 5. Conclusions and open problems Real applications (i.e. not phantoms o synthetic images) are always affected by some sort of unknown noise contamination or uncertainty in the grey levels which can deteriorate the final segmentation. Despite of this, we showed that the CV model (which is quite robust with respect to noise) and the low order statistics we introduced in [2] can attract the active contour towards the desired position. Our preliminary numerical results indicate that a iterative smoother can be used successfully but much care is need choosing the initial parameters λin,λout,µ,α. The best results in time have been obtained with the V-cycle (bidirectional) multigrid scheme (see tabla 3) where we can reach a quite coarse level of discretization without making the solution worse. This is not any more true when the Full multigrid is applied whereas the solution deteriorates in coarser levels of discretization. Interesting performances of the full multigrid and of the FAS method for the TV operator have been reported in Bruhn et al. [1] for the optic flow problem and this and this suggests to explore their indications in future work. Mathematically we believe that the difficulties associated to the minimization of (4) plead for a rigorous formulation of (5) in terms of a multivalued equation accounting for the maximal monotone graph ∂ν(u) which arise in the unconstrained formulation. This, together with a finite elements framework can lead to the definition of new duality methods [3] for the numerical resolution of (5) and shall be considered in future work. More work also has to be done with a view to hybrid methods, able to take into account natural variability and pathology. 7 J. F. Garamendi, N. Malpica, E. Schiavi Coarse grid size # Levels Symmetric difference CPU Time 25x255 3.395788 47 26x264 3.852551 49.312 27x273 3.968959 130.156 28x282 4.300444 230.766 Tabla 3: Results of V-cycle method Coarse grid size # Levels Symmetric difference CPU Time 25x255 34.718403 194.9680 26x264 23.050998 861.9240 27x273 4.863637 1444.6560 28x282 4.996675 464.8910 Tabla 4: Results of full multigrid method Referencias [1] A. Bruhn, J. Weickert, T. Kohlberger, C. Schn¨orr: Discontinuity-preserving computation of variational optic flow in real-time. In R. Kimmel, N. Sochen, J. Weickert (Eds.): Scale-Space and PDE Methods in Computer Vision. Lecture Notes in Computer Science, Vol. 3459, Springer, Berlin, 279- 290, 2005. [2] J.F. Garamendi, N. Malpica, J. Martel, E. Schiavi. Automatic Segmentation of the Liver Using Level Sets without Edges. IbPRIA 2007, Part I, LNCS 4477, pp. 161-168, 2007. [3] A. Chambolle: An Algorithm for Total Variation Minimization and Applications. Journal of Mathematical Imaging and Vision 20 (1-2): 89-97, January - March, 2004. [4] Ulrich Trottenberg, Cornelis Oosterlee, Anton sch¨uller: Multigrid Academic Press 2001. [5] D. Mumford and J. Shah, Optimal approximations by piecewise smooth functions and associated variational problems, Comm. Pure Appl. Math, Vol 42, pages 577-685, 1989. [6] Tony F. Chan and Gene H. Golub and Pep Mulet, A nonlinear primal-dual method for total variationbased image restoration, ICAOS ’96 (Paris, 1996), vol. 219, Springer-Verlag, 241–252, 1996, url: citeseer.ist.psu.edu/chan99nonlinear.html [7] Tony F. Chan, Luminita A. Vese: Active Contours Without Edges. IEEE Transactions on Image Processing, Vol. 10, No.2. IEEE, 2001. [8] Tony F. Chan, Selim Esedoglu and Mila Nikolova, Algorithms for Finding Global Minimizers of Image Segmentation and Denoising Models, UCLA-CAM-Report preprint, September 2004. [9] Freid, M. A Level Set Based Adaptive Finite Element Algorithm for Image Segmentation. Scicade 05. [10] M. E. Oman. Fast multigrid techniques in total variation-based image reconstruction, Seventh Copper Mountain Conference on Multigrid Methods, vol CP 3339, 649–659, 1996. [11] D. Strong, T. Chan. Exact Solutions to Total Variation Regularization Problems. UCLA Math Department CAM Report 96-41, October 1996. http://citeseer.ist.psu.edu/strong96exact.html [12] C. R. Vogel and M. E. Oman, Fast numerical methods for total variation minimization in image reconstruction, proceedings of SPIE 1995, San Diego, Advanced Signal Processing Algorithms, Vol. 2563, edited by F. T. Luk. http://citeseer.ist.psu.edu/vogel95fast.html [13] Vassilevski, P.S., Wade, J.G. A comparison of multilevel methods for Total Variation regularization. Electronic Transactions on Numerical Analysis. Vol 6, pp. 255-270, 1997. Kent State University. ISSN 1068-9613. [14] Joseph Savage and Ke Chen. An improved and accelerated nonlinear multigrid method for totalvariation denoising. IJCM, Vol. 82 (8), pp. 1001-1015. 8