Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) 3905–3926
Physics in Medicine & Biology doi:10.1088/0031-9155/60/10/3905
Quantitative evaluation of noise reduction and vesselness filters for liver vessel segmentation on abdominal CTA images Ha Manh Luu, Camiel Klink, Adriaan Moelker, Wiro Niessen and Theo van Walsum Biomedical Imaging Group Rotterdam, Departments of Radiology and Medical Informatics, Erasmus MC, Dr. Molewaterplein 50/60, Rotterdam, The Netherlands E-mail: [email protected] Received 11 January 2015, revised 18 March 2015 Accepted for publication 24 March 2015 Published 24 April 2015 Abstract
Liver vessel segmentation in CTA images is a challenging task, especially in the case of noisy images. This paper investigates whether pre-filtering improves liver vessel segmentation in 3D CTA images. We introduce a quantitative evaluation of several well-known filters based on a proposed liver vessel segmentation method on CTA images. We compare the effect of different diffusion techniques i.e. Regularized Perona–Malik, Hybrid Diffusion with Continuous Switch and Vessel Enhancing Diffusion as well as the vesselness approaches proposed by Sato, Frangi and Erdt. Liver vessel segmentation of the pre-processed images is performed using a histogram-based region grown with local maxima as seed points. Quantitative measurements (sensitivity, specificity and accuracy) are determined based on manual landmarks inside and outside the vessels, followed by T-tests for statistic comparisons on 51 clinical CTA images. The evaluation demonstrates that all the filters make liver vessel segmentation have a significantly higher accuracy than without using a filter (p λ2 > λ3. Second, λc = min(−λ2,−λ3) is the smaller value of the two sign-inverse smallest eigenvalues λ2 and λ3. The response is formulated as: ⎧ ⎛ −λ12 ⎞ ⎪ exp ⎜ ⎟.λc if λ1 ⩽ 0, λc ≠ 0 ⎪ ⎝ 2(α1λc )2 ⎠ , V (2) S(σ ) = ⎨ ⎛ −λ12 ⎞ ⎪ ⎪ exp ⎜⎝ 2(α λ )2 ⎟⎠.λc if λ1 > 0, λc ≠ 0 2 c ⎩
where α1 λ2 > λ3. Then the vesselness function is given by: ⎛2 ⎞ V (3) E (σ ) = K .⎜ λ1 − λ2 − λ3⎟, ⎝3 ⎠
where K = 1 −
λ2 − λ3 λ2 + λ3
. In the case of vessel structures, where λ1 ≈ 0 and λ2 ≈ λ3 ≪ 0,
VE (σ) gives a high response. Otherwise, it gives a low response. The K factor acts as a modulator which approaches 1 when λ2 ≈ λ3, in the case of a vessel-like structure and goes to 0 if the difference between λ2 and λ3 becomes large. This helps in reducing very bright, non-tube-like structures. This filter only has the scale of the Gaussian as the input parameter. 2.2. Diffusion methods
Diffusion filters solve the partial differential equation ut = div (D. ∇u) where ∇u is the gradient of the image and D is the diffusion tensor, which steers the diffusion. If the diffusion tensor D is replaced by a scalar-valued diffusivity g, the diffusion will be isotropic. In case D is an anisotropic tensor, it can model anisotropic diffusion. RPM is thus an isotropic diffusion method, as it varies only the amount of smoothing based on the local gradient magnitude, whereas VED and HDCS are anisotropic diffusion methods, which not only locally change the magnitude but also the direction of smoothing by adapting the diffusion tensor D. Each of the diffusion methods is described in more detail below. (a) RPM: Perona and Malik (1990) introduced an isotropic nonlinear diffusion described by ut = div (g (∣∇u∣). ∇u). The scalar-valued diffusivity g (∣∇u∣) is a function of the gradient magnitude ∣∇u ∣ , causing filtering in homogenous areas while retaining edges with a high gradient. Catte et al (1992) proposed the following scalar-valued diffusivity function for the non-linear diffusion using a Gaussian derivative at scale σ: −C
g(4) ( ∇ uσ ) = 1 − e( ∇ uσ 2 /λ2)4 ,
where C = 3.1488 and λ is the contrast parameter. The contrast parameter λ acts as a threshold scale for the gradient magnitude ∣∇uσ ∣. If the gradient is large compared to the contrast parameter, i.e ∣∇uσ ∣2 ≫ λ2, this results in g (∣∇uσ ∣) ≈ 0, reducing the amount of diffusion. Therefore strong edges, where ∣∇uσ ∣ is large, are preserved. Parameter σ denotes the scale of the Gaussian used to calculate the gradient. The value of σ should be chosen based on the variance of the noise and the size of the small structures we want to retain. (b) HDCS: Mendrik et al (2009) introduced HDCS as a combination of two other diffusion filters: CED and EED Weickert (1998). Both CED and EED use the structure tensor to derive a diffusion tensor. While CED is suitable for filtering tube-like structures, EED works well with flat areas and edges. The main idea of HDCS is to use a voting criterion to decide whether a local structure is tubular or non-tubular. The structure classifier is defined as: μ1 μ2 ξ(5) = − , α + μ 2 α + μ3 3909
H Manh Luu et al
Phys. Med. Biol. 60 (2015) 3905
where α = 0.001 and μ1 > μ2 > μ3 are eigenvalues of tensor product Jρ(∇ uσ ) = K ρ∗(∇ uσ ∇ uσT ) and Kρ is a Gaussian convolution kernel which acts as a smoothing factor in the tensor product Jρ. The classifier ξ ≪ 0 when the structure is tubular, ξ ≈ 0 when the structure is sphere-like (noise) and ξ ≫ 0 when the structure is plate-like (background). The HDCS diffusion tensor DH is then given by DH = Q.Λ.QT, where Q is the matrix of eigenvectors and Λ, a diagonal matrix with λ hi on the diagonal (i = 1, 2, 3 in case of 3D data), is a combination of the eigenvalues of EED (λ ei ) and CED (λ ci ): λ hi = (1 − ε )λ ci + ελ ei, (6) ⎛ μ ((λ 2ξ− ξ ) − 2μ ) ⎞ 2 3 h ⎟⎟, ε(7) = exp ⎜⎜ 4 2 λ ⎝ ⎠ h
where λh is a contrast parameter. When the local structure is tubular, ε → 0 and the diffusion is CED-like, for other structures ε → 1 and the diffusion is EED-like. (c) VED: Manniesing et al (2006) used the multi-scale Hessian filter response to drive the diffusion. The main idea is that by using eigenvalues of the Hessian matrix, we can define a diffusion tensor DV that depends on the local curvature. Let V ∈ [0, 1] be the output of a multiscale scale vesselness filter. V should be around 1 inside tubular structures and 0 elsewhere. If ∣λ1 ∣
Physics in Medicine & Biology doi:10.1088/0031-9155/60/10/3905
Quantitative evaluation of noise reduction and vesselness filters for liver vessel segmentation on abdominal CTA images Ha Manh Luu, Camiel Klink, Adriaan Moelker, Wiro Niessen and Theo van Walsum Biomedical Imaging Group Rotterdam, Departments of Radiology and Medical Informatics, Erasmus MC, Dr. Molewaterplein 50/60, Rotterdam, The Netherlands E-mail: [email protected] Received 11 January 2015, revised 18 March 2015 Accepted for publication 24 March 2015 Published 24 April 2015 Abstract
Liver vessel segmentation in CTA images is a challenging task, especially in the case of noisy images. This paper investigates whether pre-filtering improves liver vessel segmentation in 3D CTA images. We introduce a quantitative evaluation of several well-known filters based on a proposed liver vessel segmentation method on CTA images. We compare the effect of different diffusion techniques i.e. Regularized Perona–Malik, Hybrid Diffusion with Continuous Switch and Vessel Enhancing Diffusion as well as the vesselness approaches proposed by Sato, Frangi and Erdt. Liver vessel segmentation of the pre-processed images is performed using a histogram-based region grown with local maxima as seed points. Quantitative measurements (sensitivity, specificity and accuracy) are determined based on manual landmarks inside and outside the vessels, followed by T-tests for statistic comparisons on 51 clinical CTA images. The evaluation demonstrates that all the filters make liver vessel segmentation have a significantly higher accuracy than without using a filter (p λ2 > λ3. Second, λc = min(−λ2,−λ3) is the smaller value of the two sign-inverse smallest eigenvalues λ2 and λ3. The response is formulated as: ⎧ ⎛ −λ12 ⎞ ⎪ exp ⎜ ⎟.λc if λ1 ⩽ 0, λc ≠ 0 ⎪ ⎝ 2(α1λc )2 ⎠ , V (2) S(σ ) = ⎨ ⎛ −λ12 ⎞ ⎪ ⎪ exp ⎜⎝ 2(α λ )2 ⎟⎠.λc if λ1 > 0, λc ≠ 0 2 c ⎩
where α1 λ2 > λ3. Then the vesselness function is given by: ⎛2 ⎞ V (3) E (σ ) = K .⎜ λ1 − λ2 − λ3⎟, ⎝3 ⎠
where K = 1 −
λ2 − λ3 λ2 + λ3
. In the case of vessel structures, where λ1 ≈ 0 and λ2 ≈ λ3 ≪ 0,
VE (σ) gives a high response. Otherwise, it gives a low response. The K factor acts as a modulator which approaches 1 when λ2 ≈ λ3, in the case of a vessel-like structure and goes to 0 if the difference between λ2 and λ3 becomes large. This helps in reducing very bright, non-tube-like structures. This filter only has the scale of the Gaussian as the input parameter. 2.2. Diffusion methods
Diffusion filters solve the partial differential equation ut = div (D. ∇u) where ∇u is the gradient of the image and D is the diffusion tensor, which steers the diffusion. If the diffusion tensor D is replaced by a scalar-valued diffusivity g, the diffusion will be isotropic. In case D is an anisotropic tensor, it can model anisotropic diffusion. RPM is thus an isotropic diffusion method, as it varies only the amount of smoothing based on the local gradient magnitude, whereas VED and HDCS are anisotropic diffusion methods, which not only locally change the magnitude but also the direction of smoothing by adapting the diffusion tensor D. Each of the diffusion methods is described in more detail below. (a) RPM: Perona and Malik (1990) introduced an isotropic nonlinear diffusion described by ut = div (g (∣∇u∣). ∇u). The scalar-valued diffusivity g (∣∇u∣) is a function of the gradient magnitude ∣∇u ∣ , causing filtering in homogenous areas while retaining edges with a high gradient. Catte et al (1992) proposed the following scalar-valued diffusivity function for the non-linear diffusion using a Gaussian derivative at scale σ: −C
g(4) ( ∇ uσ ) = 1 − e( ∇ uσ 2 /λ2)4 ,
where C = 3.1488 and λ is the contrast parameter. The contrast parameter λ acts as a threshold scale for the gradient magnitude ∣∇uσ ∣. If the gradient is large compared to the contrast parameter, i.e ∣∇uσ ∣2 ≫ λ2, this results in g (∣∇uσ ∣) ≈ 0, reducing the amount of diffusion. Therefore strong edges, where ∣∇uσ ∣ is large, are preserved. Parameter σ denotes the scale of the Gaussian used to calculate the gradient. The value of σ should be chosen based on the variance of the noise and the size of the small structures we want to retain. (b) HDCS: Mendrik et al (2009) introduced HDCS as a combination of two other diffusion filters: CED and EED Weickert (1998). Both CED and EED use the structure tensor to derive a diffusion tensor. While CED is suitable for filtering tube-like structures, EED works well with flat areas and edges. The main idea of HDCS is to use a voting criterion to decide whether a local structure is tubular or non-tubular. The structure classifier is defined as: μ1 μ2 ξ(5) = − , α + μ 2 α + μ3 3909
H Manh Luu et al
Phys. Med. Biol. 60 (2015) 3905
where α = 0.001 and μ1 > μ2 > μ3 are eigenvalues of tensor product Jρ(∇ uσ ) = K ρ∗(∇ uσ ∇ uσT ) and Kρ is a Gaussian convolution kernel which acts as a smoothing factor in the tensor product Jρ. The classifier ξ ≪ 0 when the structure is tubular, ξ ≈ 0 when the structure is sphere-like (noise) and ξ ≫ 0 when the structure is plate-like (background). The HDCS diffusion tensor DH is then given by DH = Q.Λ.QT, where Q is the matrix of eigenvectors and Λ, a diagonal matrix with λ hi on the diagonal (i = 1, 2, 3 in case of 3D data), is a combination of the eigenvalues of EED (λ ei ) and CED (λ ci ): λ hi = (1 − ε )λ ci + ελ ei, (6) ⎛ μ ((λ 2ξ− ξ ) − 2μ ) ⎞ 2 3 h ⎟⎟, ε(7) = exp ⎜⎜ 4 2 λ ⎝ ⎠ h
where λh is a contrast parameter. When the local structure is tubular, ε → 0 and the diffusion is CED-like, for other structures ε → 1 and the diffusion is EED-like. (c) VED: Manniesing et al (2006) used the multi-scale Hessian filter response to drive the diffusion. The main idea is that by using eigenvalues of the Hessian matrix, we can define a diffusion tensor DV that depends on the local curvature. Let V ∈ [0, 1] be the output of a multiscale scale vesselness filter. V should be around 1 inside tubular structures and 0 elsewhere. If ∣λ1 ∣
