An efficient level set method for simultaneous intensity inhomogeneity correction and segmentation of MR images

Accepted Manuscript Title: An Efficient Level Set Method for Simultaneous Intensity Inhomogeneity Correction and Segmentation of MR Images Author: Tatyana Ivanovska Ren´e Laqua Lei Wang Andrea Schenk Jeong Hee Yoon Katrin Hegenscheid Henry V¨olzke Volkmar Liebscher PII: DOI: Reference:

S0895-6111(15)00177-9 http://dx.doi.org/doi:10.1016/j.compmedimag.2015.11.005 CMIG 1419

To appear in:

Computerized Medical Imaging and Graphics

Received date: Revised date: Accepted date:

8-7-2015 21-10-2015 30-11-2015

Please cite this article as: Tatyana Ivanovska, Ren´e Laqua, Lei Wang, Andrea Schenk, Jeong Hee Yoon, Katrin Hegenscheid, Henry Vddotolzke, Volkmar Liebscher, An Efficient Level Set Method for Simultaneous Intensity Inhomogeneity Correction and Segmentation of MR Images, Computerized Medical Imaging and Graphics (2015), http://dx.doi.org/10.1016/j.compmedimag.2015.11.005 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Highlights (for review)

Ac

ce pt

ed

M

an

us

cr

ip t

Highlights • We propose a novel method for simultaneous intensity inhomogeneity correction and segmentation. • Energy functional includes total variation regularization for bias field and level set function. • The minimization procedure searches the global minimum and does not depend on initialization. • The method was tested on synthetic and real images and evaluated qualitatively and quantitatively. • The approach produces results, similar or superior in quality compared to other techniques.

Page 1 of 34

*Manuscript Click here to view linked References

ip t

An Efficient Level Set Method for Simultaneous Intensity Inhomogeneity Correction and Segmentation of MR Images

cr

Tatyana Ivanovskaa,∗, Ren´e Laquaa , Lei Wangb , Andrea Schenkb , Jeong Hee Yoonc , Katrin Hegenscheida , Henry V¨olzkea , Volkmar Liebschera a

us

Ernst-Moritz-Arndt University, Greifswald, Germany Fraunhofer Institute for Medical Image Computing MEVIS, Bremen, Germany c Department of Radiology,Seoul National University Hospital, Seoul, Korea

an

b

Abstract

te

d

M

Intensity inhomogeneity (bias field) is a common artefact in magnetic resonance (MR) images, which hinders successful automatic segmentation. In this work, a novel algorithm for simultaneous segmentation and bias field correction is presented. The proposed energy functional allows for explicit regularization of the bias field term, making the model more flexible, which is crucial in presence of strong inhomogeneities. An efficient minimization procedure, attempting to find the global minimum, is applied to the energy functional. The algorithm is evaluated qualitatively and quantitatively using a synthetic example and real MR images of different organs. Comparisons with several state-of-the-art methods demonstrate the superior performance of the proposed technique. Desirable results are obtained even for images with strong and complicated inhomogeneity fields and sparse tissue structures.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Keywords: MRI, segmentation, bias field, intensity inhomogeneity, level set.

∗

Corresponding author Email address: [email protected] (Tatyana Ivanovska)

Preprint submitted to Computerized Medical Imaging and Graphics

October 21, 2015

Page 2 of 34

1. Introduction

te

d

M

an

us

cr

ip t

Intensity inhomogeneity is a common artefact that often occurs in medical images and is caused by imperfect image acquisition. For example, radio frequency field inhomogeneities appear due to either a non-uniform field itself or a non-uniform sensitivity of the receiver and transmitter coils. Such artefacts lead to undesired intensity variations in the same tissue types across the image. They hinder successful automated image segmentation, especially if the segmentation algorithms rely only on image intensities [15]. Therefore, efficient and reliable techniques for intensity inhomogeneity correction are required. In general, there are two main directions in this area [47, 18]. First, the correction can be done prospectively based on adjustments of hardware and acquisition methods, such as phantom-based calibration and multichannel transmit scanning. Such techniques can remove scanner-dependent inhomogeneities. Nonetheless, they are not applicable for already acquired data. Second, the retrospective correction is done on a post-processing level. Such methods are rather general and mainly rely on the information from acquired images and prior knowledge about imaged anatomies, including methods based on filtering, histogram analysis, surface fitting, statistical modeling, and segmentation [18]. Approaches that allow for simultaneous image segmentation and bias field correction are especially attractive, since the processes of segmentation and correction are interleaved to benefit from each other [47]. Level set methods, which are based on calculus of variations and partial differential equations, are naturally able to represent contours with complex topology and to change their topology during evolution. In the last decade, level set methods have become increasingly popular in image processing for segmentation [42, 22, 43], denoising, and enhancement [45, 23]. Region-based level set models [22, 37] aim to evolve the contour such that a certain region of interest is identified. A classical example is the well-known Chan and Vese model [7], which uses a region homogeneity assumption. Other methods based on piecewise smooth models might overcome the problem of intensity inhomogeneity. However, they are usually more computationally expensive and sensitive to initialization [27]. In this paper, we propose a novel region-based level set method, which allows for simultaneous intensity inhomogeneity correction and two-phase image segmentation. Using a generally accepted model of images with inten-

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2

Page 3 of 34

M

an

us

cr

ip t

sity inhomogeneity, we apply a local intensity clustering property proposed by Li et al. [26] and extend the energy functional to enable an explicit regularization of the bias field component. The functional consists of two parts related to segmentation and bias field correction. We reformulate the level set functional using a convexification procedure [6] and minimize each part with the Split Bregman algorithm[14, 3], which attempts to find the global minimum and does not depend on initialization. The algorithm is implemented with CUDA [36], a general-purpose computing platform on graphics processing units (GPGPU), to guarantee fast computations. The paper is organized as follows. First, a review of related works is given in Section 2. In Section 3, we shortly introduce the model and the basic level set formulation. In Section 4, we present our extended functional. The level set formulation and the minimization procedure are described in Section 4.3. The implementation is elaborated in Section 4.4. The results and findings are discussed in Section 5. This paper is summarized in Section 6. 2. Related Work

te

d

Recently, multiple retrospective algorithms for bias field correction applied to MR data have been proposed [47, 18]. These algorithms have become an essential part of segmentation pipelines. Such approaches as the nonparametric non-uniformity normalization (N3) method of Sled et al. [44] and its extension, the N4 approach, by Tustison et al. [46] have become popular due to their performance. They are applied in numerous studies, where the intensity inhomogeneity correction is introduced in a preprocessing step [49, 16, 40, 24] and the segmentation of the target organ (for instance, brain structures, breast, or liver tissue) is done thereafter. However, in certain applications prior bias field correction can not fully erase the inhomogeneity, which degrades the segmentation quality or requires semi-automatic post-processing, as mentioned, for example, in [32, 40]. Therefore, techniques, which jointly perform bias field correction and simultaneous segmentation, are required. Wells et al. [50], Guillemaud and Brady [17] presented a statistical approach based on the expectation maximization (EM) algorithm [15]. Their methods required prior knowledge on intensity distributions of different tissue classes. If there was no good initial guess for either the bias field or the classification estimate, the EM algorithm could lead to unsatisfactory results. Makarau et al. [34] developed a method based on the Mean shift approach [10] for breast images. It was tested on

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3

Page 4 of 34

te

d

M

an

us

cr

ip t

images with extracted regions of interest, i. e., without background. The method could lead to misclassifications, since it relied also on successful extraction of breast tissue. Pham and Prince added a multiplicative term to the objective function in the fuzzy C-Means (FCM) algorithm [2] to account for multiplicative intensity inhomogeneities [38]. However, the algorithm was rather sensitive to noise and computationally expensive. There are multiple extensions of this approach. Ahmed et al. proposed the labeling of a pixel (voxel) to be influenced by the labels in its immediate neighborhood [1], which made the algorithm insensitive to salt and pepper noise. Nevertheless, they did not consider that the algorithm might lead to unsatisfactory results without regularization terms. Li et al. also included the neighborhood statistics in the functional [30] and kept both regularization terms of the bias field, proposed by Pham and Prince [38]. They obtained promising results on images with a high noise level. However, these methods were tested on either phantom or brain images, and they stayed rather computationally demanding. Chen and Giger proposed a fast FCM based method [8]. The bias field was estimated in each iteration on the current tissue class centroids and the membership values of voxels. Then the estimated field was smoothed with an iterative low pass filter. The authors tested their algorithm on clinical breast MR images. That algorithm enhanced the signal intensities of fatty breast tissues, but the intensity of some fibroglandular tissue was also enhanced. This might change the overall contrast of images leading to wrong segmentation results, as reported by Lin et al. [32]. Lin et al. proposed a combined procedure for bias field correction and segmentation, consisting of the nonparametric non-uniformity normalization (N3) method [44] and a modified FCM [32]. In the first step of their algorithm, the breast region was extracted from the body in a semiautomatic manner. The algorithm was compared to the coherent local intensity clustering (CLIC), another FCM based method proposed by Li et al. [28], which showed promising results in removing strong inhomogeneities and, generally, did not require stripping of the region of interest. The study showed that both algorithms produced results comparable in quality. Cui et al. presented a localized FCM based method, motivated by the CLIC method of Li et al. [28], with the objective function including local and global information, which allowed for noise suppression [11]. The method was tested on synthetic images, X-ray vessel images, and MR brain images. Li et al. extended their CLIC method and converted the proposed functional to the level set formulation to improve performance. The results were promising even without additional (usually

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4

Page 5 of 34

te

d

M

an

us

cr

ip t

manual or semiautomatic) extraction of the region of interest [26]. Zhang et al. presented an extension of this approach, integrating a maximum likelihood function over the entire image domain to the level set energy [55, 54]. Zhan et al. [53] introduced a Gaussian distribution with a bias field as a local region descriptor. The variational model included three-phase level set model, and the method was tested on brain images. Nonetheless, the original level set method and its extensions are often sensitive to initialization in real-world images and rather computationally expensive. Moreover, that approach of Li et al. [26] sometimes produced misclassifications for images that consisted of two main classes, where the intensities of one tissue class and the background were similar. Ivanovska et al. [20] introduced a mask-based method to further improve results. However, the approach found only a local minimum, depended on initialization, and its computational costs made it practically inapplicable to massive data processing. Wang and Pan presented a level set method with an image-guided regularization term, to restrict the level set function [48]. The method was evaluated only on brain images. The bias field was also formulated in level set terms. Nevertheless, the level set procedure was dependent on the initial countour and might get stuck in a local minimum. The idea to convexify an energy functional and search the global minimum is becoming popular also for intensity inhomogeneity correction and segmentation algorithms. Li et al. [25] extended the basic model and introduced a method, which jointly performed bias field estimation and the tissue membership functions in an FCM-like energy minimization process to optimize two multiplicative intrinsic components of an MR image. The bias field was presented as a linear combination of a given set of smooth basis functions. The introduced energy functional is convex and, therefore, the minimization seeks the global minimum. The method was evaluated on brain images. Moreno et al. [35] proposed a multi-phase convex level set procedure for brain image segmentation using Chambolle’s dual minimization approach [5]. However, no bias field component was taken into account here. Yang and Wu formulated the local and global intensity fitting energy model and applied the Split Bregman method to it [51, 52]. The bias field regularization remained implicit though. Moreover, the authors reported that their method was sensitive to weighting parameter selection [51]. Ivanovska et al. proposed a more flexible model [21], which allowed for an explicit regularization of the bias field component. An efficient minimization procedure, namely, the Split Bregman algorithm [14] was applied both for contour

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

5

Page 6 of 34

an

us

cr

ip t

evolution and bias field correction processes. This paper presents an extended version of our preliminary work [21], namely, a fast novel level set based method for simultaneous bias field correction and segmentation, which allows for explicit bias field manipulation and uses a minimization procedure that searches the global minimum. Here, we provide a more thorough and accurate description of the model, level set functional derivation and minimization, introduce a full GPU-based implementation, test the approach both on synthetic and multiple MR data of different organs, and present quantitative result evaluation. The method requires no prior stripping of the region of interest. The utilized level set procedure requires no (re-)initialization and is less computationally demanding than the classical level set approaches. 3. Theoretical Background

d

M

3.1. Image Model The 2D image Im : Ω → R is defined on a continuous domain Ω, such that Im = bJ + n, (1)

te

where J is the true image, b is the intensity inhomogeneity component, and n is the additive zero-mean Gaussian noise. This is a standard multiplicative model that originates from MR physics [4] and is used in multiple works on intensity inhomogeneity correction [47]. Similarly as in the model presented by Li et al. [26], we assume that the bias field component b : Ω → R is slowly varying and essentially constant within a small circular neighborhood Oy (b(x) ≈ b(y), for x ∈ Oy ).

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3.2. CLIC Energy Functional and Level Set Minimization Here, we shortly describe the original CLIC energy formulation and the level set evolution procedure presented by Li et al. [26]. 3.2.1. Energy Functional The true image J is assumed to be piecewise constant, i. e., it consists of N regions Ω1 , . . . ΩN with intensities c = (c1 , . . . , cN ) and the bias field b. The noise Rn in the neighborhood Oy is minimized by the following fidelity term PN 2 i=1 Ωi ∩Oy (b(y)ci − Im(x)) dx, where y denotes the center of the circular neighbourhood Oy with a radius ρ: Oy = {x : |x − y| ≤ ρ}. 6

Page 7 of 34

i=1

K(y − x)(b(y)ci − Im(x))2 dx. Ωi ∩Oy

(2)

cr

Ey =

N Z X

ip t

Thereafter, a truncated Gaussian function K is introduced to account for the piecewise smoothness of the model. The clustering criterion for classifying the intensities in Oy is

Ey =

N Z X i=1

us

Since K(y − x) = 0 for x 6∈ Oy , the clustering criterion can be written as K(y − x)(b(y)ci − Im(x))2 dx. Ωi

Ωi

Ωi

te

with

d

Then, the order of integrations is exchanged: ! Z X Z X N N Z 2 ei (x)dx, K(y − x)(b(y)ci − Im(x)) dy dx = E= i=1

ei (x) =

Z

(3)

(4)

M

i=1

an

The clustering energy in the entire domain Ω is then formulated as ! Z Z X N Z K(y − x)(b(y)ci − Im(x))2 dx dy. E = Ey dy =

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

K(y − x)(b(y)ci − Im(x))2 dy.

(5)

i=1

(6)

The level set membership function Mi (φ(x)) is added as a multiplicative term to the clustering energy, and it is defined as follows Z X N ei (x)Mi (φ(x))dx, (7) E(φ, c, b) = i=1

where φ is a level set function. Two phase segmentation considered, where N = 2, M1 (φ) = H(φ) and M2 (φ) = 1 − H(φ). H(φ) is the Heaviside function and its derivative is the Dirac delta function δ: d H(x) = δ(x) dx  1 x≥0 H(x) = 0 x<0 7

Page 8 of 34

F(φ, c, b) = E(φ, c, b) + νL(φ) + µRp (φ),

ip t

3.2.2. Level Set Minimization The energy from Equation 7 is taken as the data term in the variational level set formulation

(8)

te

d

M

an

us

cr

Rwhere µ and ν are weighting factors, H is the Heaviside function, L(φ) = R |∇H(φ)|dx is the regularization term for contour smoothing. Rp (φ) = p(|∇φ|)dx is the distance regularization term, which allows for maintaining the signed distance property of φ without additional reinitialization [29], and p(s) = 0.5(s − 1)2 . The minimization of the energy F(φ, c, b) results in the image segmentation given by the level set function φ and the estimation of the bias field b. It is achieved by an iterative process, where on each iteration the energy is minimized with respect to each of its variables φ, c, b. The minimization is obtained by solving following equations (in each equation the other two variables are fixed):   ∂φ ∂F ∇φ (9) =− = −δ(φ)(e1 − e2 ) + νδ(φ)div ∂t ∂φ |∇φ| + µdiv (dp (|∇φ|)∇φ) , R (b ∗ K)ImMi (φ)dx cˆi = R 2 , (10) (b ∗ K)Mi (φ)dx PN ci Mi (φ)) ∗ K ˆb = (Im , (11) PN i=12 ( i=1 ci Mi (φ)) ∗ K

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where ∗ denotes the convolution operator. The finite difference scheme proposed by Li et al. [29] is applied for the level set evolution, and fixed space steps △x = △y = 1 are taken. For application, the Heaviside function H and the Dirac delta δ (the Dirac delta function) are replaced by their smooth approximations of H and δ, denoted as Hǫ and δǫ , as ǫ → 0:   1 2 x Hǫ (x) = 1 + arctan , (12) 2 π ǫ ǫ δǫ (x) = , (13) 2 π(ǫ + x2 ) with ǫ = 1, i. e. ǫ is equal to the fixed space steps. 8

Page 9 of 34

4. Method

E 1 = µb

Z X N

an

us

cr

ip t

4.1. Proposed Model with Explicit Bias Field Smoothing Parameter In the original CLIC model, the truncated Gaussian function K serves for an implicit smoothing of the bias field term b. The level of smoothness is controlled by the standard deviation σ of the Gaussian kernel K, i. e., one can manipulate b only by changing σ, which can be insufficient in presence of complicated intensity inhomogeneities and organ structures. Moreover, such a smoothing term requires computation of four convolutions in each iteration, which is rather computationally expensive. We propose to formulate the clustering energy with explicit bias field regularization, which requires no convolutions, as follows 2

(b(x)ci − Im(x)) dx +

Z

|∇b|T V dx,

(14)

M

i=1

R

te

d

where |∇b|T V dx is the total variation [41] of b, and µb is a constant, which denotes the weight between the data term and the regularization term. The total variation term serves for explicit smoothing of the bias field. The level set membership function Mi (φ(x)) is added as a multiplicative term to the fidelity term of the clustering energy: Z X N

(b(x)ci − Im(x))2 Mi (φ(x))dx,

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(15)

i=1

where φ is a level set function. We consider only the two phase cases, where N = 2, M1 (φ) = H(φ), and M2 (φ) = 1−H(φ) with the Heaviside function H. Let e∗i (x) = (b(x)ci − Im(x))2 . This term serves both as the data component for region-based level set function φ and the fidelity component for b. Finally, our clustering energy is formulated as follows: Enew (φ, c, b) = µb

Z X N

e∗i (x)Mi (φ(x))dx

i=1

+

Z

|∇b(x)|T V dx,

(16)

where µb is a weighting constant.

9

Page 10 of 34

Fnew (φ, c, b) = µφ Enew (φ, c, b) + L(φ),

ip t

4.2. Level Set Formulation and Convexification The variational level set formulation is formulated as (17)

us

cr

where µφ is a factor, which controls the relative influence of the data term and of the smoothing term. We omit the regularization term Rp (φ) (cf. 8), since no standard level set minimization procedure will be used. The variations of energy in Equation 17 with respect to the level set function φ for fixed c, b lead to the following gradient descent scheme:

an

∂φ ∇φ = −δ(φ)γ1 (e∗1 − e∗2 ) + δ(φ)div( ), ∂t |∇φ|

(18)

te

d

M

where the Dirac delta function is the derivative of the Heaviside function H, γ1 = µφ µb . In standard numerical implementation, the smooth approximations Hǫ and δǫ from Equation 12 are used. Equation 18 is similar to the gradient scheme for φ for the Chan-Vese model [6], therefore, the convexification procedure described by Chan et al. [6] can be utilized here. The level set functional is reformulated as follows. First, the approximate Heaviside function is omitted altogether, which gives   ∂φ ∇φ ∗ ∗ = −γ1 (e1 − e2 ) + div . (19) ∂t |∇φ|

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

This equation is the gradient descent for the following energy Z Z |∇φ|T V dx + γ1 (e∗1 − e∗2 )φdx.

(20)

According to Theorem 2 from [6], a global minimizer can be found by carrying out the following convex minimization: Z Z min |∇φ|T V dx + γ1 (e∗1 − e∗2 )φdx, (21) 0≤φ≤1

and then setting {x : φ(x) ≥ ν} for ν ∈ [0, 1]. We take ν = 1 in our implementation. For more details of the convexification procedure, we refer to the original work by Chan et al. [6].

10

Page 11 of 34

us

cr

ip t

Second, we transform our initial energy functional (cf. Equation 17) as follows. Z FConvex (φ, c, b) = µφ µb [φ(x)(b(x)c1 − Im(x))2 + (22) Z (1 − φ(x))(b(x)c2 − Im(x))2 ]dx + µφ |∇b(x)|T V dx Z + |∇φ(x)|T V dx,

M

an

4.3. Efficient Level Set Minimization FConvex (φ, c, b) is minimized with respect to three parameters, namely, φ, c, and b, given the other two updated in previous iteration. We have noticed that the efficient Split Bregman algorithm [14, 13, 3] can be utilized for minimization with respect to φ and b. Below we describe the procedures in more detail for each variable.

te

d

4.3.1. Minimization with respect to φ The main steps of the algorithm are as follows. First, the auxiliary variable d = ∇φ is introduced. λφ is a weighting factor. The constraint d = ∇φ is enforced using the Bregman iteration approach [3]: Z  k+1 k+1   (φ , d ) = arg min (|d| + γ1 (e∗1 − e∗2 )φ   φ∈[0,1],d  λφ (23) + |d − ∇φ − bv k |2 )dx, k≥0    2   k+1 bv = bv k + ∇φk+1 − dk+1

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The L1 and L2 components of this functional are split [14] so that the minimization is performed efficiently with respect to φ and d separately. These steps are: φk+1 = min γ1 (e∗1 − e∗2 )φ + φ

dk+1 = min |d| + d

λφ k |d − ∇φ − bv k |2 , 2

λφ k |d − ∇φk+1 − bv k |2 . 2

(24) (25)

The minimizing solution φk+1 is characterized by the optimality condition: ω △ φ = γ1 (e∗1 − e∗2 ) + λφ div(bv k − dk ),

φ ∈ [0, 1].

(26)

11

Page 12 of 34

dk+1 =

ip t

A fast approximate solution for φ is provided with a Gauss-Seidel iterative scheme. Finally, the minimizing solution for dk+1 is given by softthresholding: (∇φk+1 + bv k ) 1 max(|∇φk+1 + bv k | − , 0) k+1 k |∇φ + bv | λφ

(27)

us

cr

4.3.2. Minimization with respect to b Now we minimize the functional FConvex (φ, c, b) with respect to b. We omit all additive parts of FConvex , which do not depend on b and disappear with the differentiation, as well as the integrals for the sake of shortness. Now, the minimization problem looks as follows: min µb ((bc1 − Im)2 φ + (bc2 − Im)2 (1 − φ)) + |∇b|T V .

an

b

(28)

te

d

M

This equation is similar to the total variation minimization problem [41], therefore, the Split Bregman algorithm [14] can be applied here. We redefine M1 = φ and M2 = (1 − φ), µb = µ2 and take c1 and c2 out of the brackets: µ min ((b − Im/c1 )2 M1 c21 + (b − Im/c2 )2 M2 c22 ) + |∇bx | + |∇by |. (29) b 2 Then, the auxiliary variable d = ∇b is introduced. This yields µ min ((b − Im/c1 )2 M1 c21 + (b − Im/c2 )2 M2 c22 ) b,dx ,dy 2 + |dx | + |dy | λb λb (30) + |dx − ∇x b|2 + |dy − ∇y b|2 . 2 2 Finally, the Bregman iteration is applied to strictly enforce the constraints: µ min ((b − Im/c1 )2 M1 c21 + (b − Im/c2 )2 M2 c22 ) b,dx ,dy 2 λb + |dx | + |dy | + |dx − ∇x b − bvxk |2 2 λb + |dy − ∇y b − bvyk |2 , (31) 2 where the values bvxk and bvyk are chosen through the Bregman iteration. The subproblem with respect to b is µ bk+1 = min ((b − Im/c1 )2 M1 c21 + (b − Im/c2 )2 M2 c22 ) b 2 λb λb (32) + |dx − ∇x b − bvxk |2 + |dy − ∇y b − bvyk |2 , 2 2

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

12

Page 13 of 34

which has the optimality condition + λb ∇Tx (dkx − bvxk ) + λb ∇Ty (dky − bvyk ),

(33)

µ(M1 c21

(µ1 I − λ △ I)bk+1 = µ2 Im

+

M2 c22 ),

cr

where I is the identity matrix. We redefine µ1 = µ(M1 c1 + M2 c2 ). Then

ip t

(µ(M1 c21 + M2 c22 )I − λb △)bk+1 = µIm(M1 c1 + M2 c2 )

us

+ λb ∇Tx (dkx − bvxk ) + λb ∇Ty (dky − bvyk ).

µ2 =

(34)

an

A fast approximate solution for b can be found with Gauss-Seidel iterations. Similarly as for φ, the optimal solution for d is found using the shrinkage operator [14]: 1 (∇x bk+1 + bvxk ) max(|∇x bk+1 + bvxk | − , 0) k+1 k |∇x b + bvx | λb k+1 k (∇y b + bvy ) 1 = max(|∇y bk+1 + bvyk | − , 0) k+1 k |∇y b + bvy | λb

dk+1 y

(35)

M

dk+1 = x

The Bregman vector is updated as

d

bvxk+1 = bvxk + (∇x bk+1 − dk+1 x ), =

bvyk

te

bvyk+1

+ (∇y b

k+1

−

(36)

dk+1 y ).

4.3.3. Minimization with respect to c For fixed φ and b, the optimal c that minimizes the energy FConvex is obtained explicitly by the following equations: R b(x)Im(x)φ(x)dx R c1 = (37) b2 (x)φ(x)dx R b(x)Im(x)(1 − φ(x))dx R , c2 = b2 (x)(1 − φ(x))dx

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where φ(x) ∈ [0, 1].

4.4. Implementation Although the three partial minimization processes are of equal rank, we keep their historical order in the current implementation. Namely, the minimization for φ is the outer one, and the minimizations with respect to b and c are nested in it. The overall scheme is shown in Algorithm 1. We leave the investigation of influence of the minimization orders on the segmentation and bias field correction results for future work. 13

Page 14 of 34

ALGORITHM 1: Minimization scheme.

M

an

us

cr

ip t

Initialize φ with Image Im ∈ [0, 1] Compute initial b and c while |φOld − φN ew | > εφ AND outIter < maxOutIter do while innerGSIter < maxGSIter do Gauss-Seidel Iterations for φ (cf. Equation 26) end while Update d and Bregman vector bv for φ (cf. Equations 23, 27) Update c (Equation 37) while |bOld − bN ew | > εb AND outIterb < maxOutIterb do Gauss-Seidel Iterations for b (cf. Equation 34) Update d and Bregman vector bv for b (cf. Equations 36, 35) end while end while

te

d

CUDA Version The algorithm fits into a single instruction multiple-data (SIMD) programming model [33]. We replace the Gauss-Seidel iterations with the redblack successive overrelaxation (SOR) [39] and set the overrelaxation parameter ω = 1.5. Every part of the algorithm is implemented in CUDA as a separate kernel. The calculation process is started from the CPU side, whereas each calculation step is executed on the GPU side. The process of CPU-GPU interactions is shown in Algorithm 2. To calculate the cluster centers c, we use the reduction operation [36] to compute the sums efficiently.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

5. Experimental Results and Discussion

The parameter selection in our model is rather straightforward. We have two data fidelity terms µφ and µb and two Bregman factors λφ and λb . Both data fidelity terms affect the segmentation results for the φ process, and only µb affects the bias field correction b process. The values of the fidelity terms define the amount of smoothing, i. e., the smaller the fidelity terms, the more impact has the other term (|∇φ(x)| in the level set φ process) and the smoother the contours will be. Generally, a good strategy is to fix one of the fidelity terms and manipulate the other one. The Bregman factors serve for the Split Bregman minimization procedures. Exact values of these . parameters usually can be described by the following formulae: λφ = µφ 10−k 14

Page 15 of 34

cr

d

M

an

us

CPU-GPU: Initialize φ, c, b while CPU: F lagOU T Iterationsφ do while CPU: F lagIN Iterationsφ do GPU: Red-Black SOR for φ end while GPU: Update d and Bregman vector bv for φ GPU Update c while CPU: F lagOU T Iterationsb do while CPU: F lagIN Iterationsb do GPU: Red-Black SOR for b end while GPU: Update d and Bregman vector bv for b GPU: Check F lagOU T Iterationsb end while GPU: Check F lagOU T Iterationsφ end while CPU-GPU: Return φ, c, b

ip t

ALGORITHM 2: The scheme of CPU-GPU interactions.

te

Table 1: State-of-the-art methods for simultaneous segmentation and bias field correction for comparison.

Author Li et al. [26] Li et al. [28] Zhang et al. [55, 54]

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Method Level set FCM Level set

Implemented C. Li T. Ivanovska K. Zhang

Language Matlab Matlab Matlab

. and λb = µb 10+l , where k and l must be defined empirically. In our examples, we fixed µb = 1, selected µφ ≫ µb and k, l ∈ [3, 6]. We also observe that the Bregman factor λb influences th smoothness level of the extracted bias field. The parameters of other methods are selected according to the recommendations of the authors in original papers. We demonstrate the performance of our method as well as compare it to other state-of-the-art techniques in terms of speed and accuracy. The methods that do not require any prior extraction of the tissue of interest have been chosen. Performance is tested on a Intel(R)Xeon(R) CPU [email protected] with 24Gb RAM and NVIDIA Tesla C2070. The state-of-the-art methods are listed in Table 1. 15

Page 16 of 34

ip t

The estimated bias field b is used for intensity inhomogeneity correction, and the bias corrected image is computed as Im/(b + ǫ), where ǫ = 0.00001.

M

an

us

cr

5.1. Synthetic Example We tested the performance of the methods on a synthetic image shown in Figure 1a. The image imitates a real MR slice: it has homogeneous background with intensity level 0 and foreground with a local gradient. The image is affected by Gaussian noise with σ = 10.

(a)

(c)

(d)

te

d

(b)

(f)

(e)

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(g)

Figure 1: a) A synthetic 128 × 128 image with additive Gaussian noise (σ = 10). b)-g) Results: b) The proposed method. c)The FCM method of Li et al. [28]. d)The level set method of Li et al. [26]. e) The level set method of Zhang et al. [55]. f) The level set method of Li et al. [26] with another initialization. g) The level set method of Zhang et al. [55] with another initialization. The level set methods produce more accurate results with smoother contours, but they are rather sensitive to initialization. The proposed method overcomes this limitation.

For our method, we set µb = 1, µφ = 106 , λb = 300, λφ = 103 . We evaluated the segmentation results (cf. Figure 1) using a region-based measure, 2|S1 ∩ S2 | , namely, the Dice coefficient [12], which is computed as DSC = |S1 | + |S2 | where S1 and S2 are two segmentation results (in our case, the one provided by the automatic method and the one provided by the expert readings). The 16

Page 17 of 34

DSC 0.989 0.969 0.82 0.95

cr

Method Presented method Level set method of Li et al. [26] FCM method of Li et al. [28] Level set method of Zhang et al. [55, 54]

ip t

Table 2: DSC coefficients for the synthetic image.

te

d

M

an

us

results are presented in Table 2. As it can be observed in Figure 1 and Table 2, level set based methods produce more accurate results with smoother contours, whereas the FCM-based approach undersegmented the circular object. Nevertheless, the level-set approaches have several issues. First, classical gradient descent based procedure used by level set approaches of Li et al. [26] and Zhang et al. [55, 54] tend to converge to a local minimum and, therefore, the result depends on the location of initial contours. For example, we change the location of initial contours and run the algorithms with exactly the same parameters as in the “successful cases” (cf. Figure 1), on the example image, which lead to unsatisfactory results as documented in Figures 1f and 1g. To our experience, the method of Zhang et al. [55] is more sensitive to the initial contour location. The proposed method overcomes this limitation, since it does not depend on initial contours, and the minimization procedure uses the Split Bregman method, which seeks the global minimum. The result is shown in Figure 1b. Second, the level set based procedures are rather computationally expensive. Of course, running time of a method is heavily dependent on implementation, and one can not compare the above-mentioned methods implemented in Matlab to an approach, which runs on a GPU. In our previous work [21], we reimplemented the method of Li et al. [26] in C++ and compared sequential CPU implementations of two methods. It appeared that our method in that implementation was approximately 15 times faster. Here, to make the approaches comparable, we reimplement the method of Li et al. [26] in CUDA [20] and observe that our method is slightly more computationally expensive, which is explained by two evolution processes (for b and φ) there. For example, a call of 1 iteration for the GPU implementation of Li et al. [26] with σ = 4 takes 0.15 seconds on a 512 × 512 image, whereas our method requires 0.23 seconds. Our algorithm requires approximately 10 seconds for a 512 × 512 image on NVIDIA Tesla C2070, whereas the algorithm of Li et

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

17

Page 18 of 34

Table 3: Detailed information about example MR images.

ip t

Image Whole body TIRM (2.1 × 1.6 × 5 mm) T1 VIBE (1.48 × 1.48 × 3 mm) T1 TWIST native (0.9 × 0.7 × 1.5 mm)

cr

Organ A: Ankles B: Liver C1,2: Breast

us

al. [26] takes 5−10 seconds for the same image dependent on the convergence speed.

te

d

M

an

5.2. MR Examples Real world images are clearly more complicated than any synthetic examples due to a combination of different artefacts, such as intensity inhomogeneity, noise, and partial volume effect. Moreover, different organ structures often represent a significant challenge. Since our method is designed for two class segmentation, we select several examples, where the image can be roughly separated in two classes (dark and bright), apply the methods and discuss the results. In Figure 2, example slices from different MR sequences are shown. The ankle and breast data were acquired with 1.5-T MR Siemens Avanto imager, and the liver data were acquired with 3-T MR Siemens Trio Tim imager. The detailed image information is given in Table 3.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

A

C1

B

C2

Figure 2: Example slices from different MR sequences depicting ankles (A), liver (B), and breast(C).

The first example is an image with depicted ankles. The image exhibits intensity inhomogeneities on outer ankle edges, where intensities are significantly higher than in the rest of the tissue. This example is rather straightforward to segment, since the objects are big and have quite homogeneous structure. All methods produce here satisfactory segmentation results (cf. Figure 3), the FCM-based method slightly undersegments the lower region, 18

Page 19 of 34

us

cr

ip t

and the intensity inhomogeneity is not fully corrected. The level set based correction results are visually more homogeneous with no bias field present.

(b) Level set of Li et al. [26].

d

M

an

(a) The presented global level set method.

te

(c) FCM of Li et al. [28].

(d) Level set of Zhang et al. [55, 54].

Figure 3: Segmentation results (red contours) overlaid with the initial images and bias field corrected images. The initial contours are marked in blue.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The liver example is a cut-out of a slice from a 3 Tesla T1-weighted dynamic sequence. Here, there is a clear intensity gradient within the liver, where the upper inner organ’s part is significantly darker than the rest tissue. The highest intensities are reached in lower and outer liver parts. As one can observe in Figure 4, all methods produce decent results on this example. Strong intensity inhomogeneity is removed in the whole organ apart from a small region near the upper liver boundary, which did not hinder a successful segmentation. The native breast examples have a rather sophisticated bias field. Here, intensity is gradually decreasing from breast tops, where nipples are located, to the region, where heart and lungs are located. Within breast tissue itself, there is also another intensity gradient. The regions, which are close to the breast-air boundary, have higher intensity values than other regions. More19

Page 20 of 34

ip t cr

(b) Level set of Li et al. [26].

(c) FCM of Li et al. [28].

(d) Level set of Zhang et al. [55, 54].

M

an

us

(a) The presented global level set method.

d

Figure 4: Segmentation results (red contours) overlaid with the initial images and bias field corrected images. The initial contours are marked in blue.

te

over, breast parenchyma (dark regions) is sparse structured and interconnected with fatty tissue (light grey regions). Hence, it is not straightforward for the algorithms to properly estimate intensity inhomogeneity and not to misclassify parenchymal tissue, assigning it to the fatty tissue class. The first breast example (C1) contains relatively low amount of parenchyma, which is sparse, and the breast muscle is not well separated from the breast fatty tissue there. The second breast example (C2) contains more parenchymal regions and the breast muscle tissue is clearly visible. The results are presented in Figures 5 and 6. The FCM-based method undersegments parenchymal regions in both examples and fails to accurately detect the sophisticated breast muscle region in the example C1. Triggering the method’s parameters has not lead to result improvement. The method of Zhang et al. [54] has not been able to converge to a stable solution in a reasonable number of iterations. We show intermediate results for this method after 700 iterations. The approach of Li et al. [26] produces plausible results with some undersegmentation of the breast parenchyma and muscle regions. However, this method is based on a gradient descent minimization scheme, and it becomes

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

20

Page 21 of 34

ip t cr

(b) Level set of Li et al. [26].

(c) FCM of Li et al. [28].

(d) Level set of Zhang et al. [55, 54]. Intermediate result after 700 iterations.

M

an

us

(a) The presented global level set method.

te

d

Figure 5: Segmentation results (yellow contours) overlaid with the initial images and bias field corrected images. The initial contours are marked in red.

really sensitive to initialization in complex real-world images, such as breast data. In Figure 7, the segmentation results of the level set approach by Li et al. [26] with different initializations are shown. The results are obtained with the same parameter values, as in Figure 6b, namely, with ν = 255, µ = 1, σ = 15, dt = 0.1. As in the synthetic case (cf. Figure 1), the location of initial contours influences segmentation results. Specifically, the sparse parenchymal regions that have been detected in Figure 6b are missing in Figure 7. Unlike other methods, our approach overcomes the above mentioned problems due to its flexible model with a term for direct smoothing bias field component and the global optimization procedure. Hence, the sparse parenchymal regions as well as the breast muscle boundary were segmented appropriately, as documented in Figures 5 and 6. We have shown that breast images pose a significant challenge for intensity inhomogeneity correction methods. Now we confirm our observations quantitatively using inhomogeneity correction and segmentation met-

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

21

Page 22 of 34

ip t cr

(b) Level set of Li et al. [26].

(c) FCM of Li et al. [28].

(d) Level set of Zhang et al. [55, 54]. Intermediate result after 700 iterations.

M

an

us

(a) The presented global level set method.

te

d

Figure 6: Segmentation results (yellow contours) overlaid with the initial images and bias field corrected images. The initial contours are marked in red.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 7: The segmentation results for the level set procedure [26] depend on the position of the initial contour (marked red).

rics. The results are compared to the method of Li et al. [26]. Since we do not have the “true” inhomogeneity field to compare with the extracted one, we only rely on the intensities of corrected images and segmentation results. Therefore, we apply intensity variation-based and 22

Page 23 of 34

ip t

segmentation-based (region-based) metrics. For comparisons, we use manually annotated ground truth provided by our radiology experts.

te

d

M

an

us

cr

Intensity Variation Based Metrics. The results of inhomogeneity correction are usually evaluated with two measures: coefficient of variation (CV) and coefficient of joint variation (CJV) [47]. The coefficient of variation is the quotient between standard deviation (σ) and the mean value (µ) of a tissue class σ(C) (C): CV (C) = . CV is invariant to multiplicative uniform intensity µ(C) transformation, but sensitive to the additive part (brightness). Moreover, CV can be computed for one tissue class, but does not provide any information about the overlap between intensity distributions of different tissue classes. The coefficient of joint variation overcomes the above mentioned drawbacks of CV [31], since it estimates the overlap between two tissue classes (C1 σ(C1 ) + σ(C2 ) . Additionally, and C2 ). CJV is defined as CJV (C1 , C2 ) = |µ(C1 ) − µ(C2 )| we also compute the relative change of CJV before and after the correction |CJVaf ter − CJVbef ore | × 100%. CJVbef ore We measure CV for the whole breast tissue (fatty and parenchymal tissues together), and CJV as well as the relative CJV for fatty and parenchymal tissues in original and bias field corrected data (cf. Table 4). One can observe that the presented approach reduces the inhomogeneity in both images. The relative variation change is larger than 63%. The level set method of Li et al. [26] yields the CV and CJV values higher than for the proposed method.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Table 4: Coefficients of variation (CV) for breast tissue and coefficients of joint variation (CJV) for fatty and parenchymal tissues in the original and corrected, and relative CJV for original and corrected data.

Image

C1 C2

Original

CV 0.44 0.46

CJV 5.41 1.71

Corrected Our Method Level set of [26] CV CJV Rel CJV % CV CJV Rel CJV % 0.11 1.16 78.1 0.11 1.47 72.8 0.15 0.61 63.8 0.21 0.905 47.07

Segmentation Based Metrics. Since the intensity variation based metrics are indirect [9] and their usage can be unreliable we also evaluated the segmen23

Page 24 of 34

Table 5: Dice’s coefficients for automatic results compared to the manual segmentation.

Level set of [26] Breast Parenchyma 0.93 0.31 0.981 0.52

cr

C1 C2

Our Method Breast Parenchyma 0.97 0.65 0.989 0.89

ip t

Image

te

d

M

an

us

tation results using the Dice coefficient. Of course, the segmentation results require a range of basic post-processing steps, as extracting connected components of breast tissue, smoothing breastpectoral muscle and breast-air boundaries, and analysing 3D connected component. We applied the proposed method for the automated breast density estimation [19] and validated it using 37 3D MR mammograms. The average Dice’s Similarity Coefficient (DSC) values were 0.96±0.0172 and 0.83±0.0636 for breast and parenchymal volumes. Here, for a proof of concept we apply the postprocessing steps to the example 2D slices processed by our approach and the approach of Li et al. [26] and compare the breast and parenchymal regions to manually obtained masks. In Table 5, the results of the slice-wise comparison are presented. One can observe that the method of Li et al. [26] slightly misclassifies the breastmuscle boundary (cf. Figure 5), the DSC for the image C1 is lower than for our method, namely, 0.93, whereas our method achieves a DSC of 0.97. For the example C2, where the breast muscle boundary is clearly visible (cf. Figure 6), the methods produce high quality results with a DSC value of 0.98. In parenchymal tissue segmentation, the method of Li et al. [26] achieves quite low DSCs, when compared to the proposed approach, which confirms our qualitative observations.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

5.3. Further Evaluation on MR Breast Data Here, we extend our test set by randomly selecting 30 MR images (5 slices from 6 datasets) and apply the proposed method, the level set method of Li et al. [26], and a combined pipeline, consisting of the N4 technique [46] and a subsequent intensity clustering with the classical fuzzy C-Means approach [2] into three classes. The segmentation performance is evaluated using the Dice coefficients. The box plots for breast volume and parenchymal tissue segmentation results are presented in Figure 8. The proposed method is superior than the other techniques both for the breast volume and parenchymal tissue segmentation. The N4 method could 24

Page 25 of 34

DSC Coefficients for Breast Parenchyma

0.9

0.8

0.8

0.6

0.7

0.4

0.6

0.2

0.5

cr

DSC

1

ip t

DSC Coefficients for Breast Volume 1

DSC

0 Proposed Method

Li et al. Levelset

N4+FCM

(a) Breast tissue segmentation results.

Proposed Method

Li et al. Levelset

us

N4+FCM

(b) Breast parenchyma segmentation results.

an

Figure 8: Segmentation results for three approaches. The presented approach is superior both in breast and parenchymal tissue segmentation.

te

d

M

not erase the intensity inhomogeneities completely, which is reflected in segmentation results, where the parts of breast tissue are assigned to different intensity classes. The approach of Li et al. [26] usually detects the breast tissue rather accurately, but assigns the parenchymal tissue parts to the breast tissue class. The manual segmentations are cut approximately in the sternum position (cf. Figure 9a). For quantitative evaluation the automatically obtained results are cut in the same line after processing. Bias field correction and segmentation results are documented in Figures 9b, 9c, 9d.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

5.4. Limitations of the Proposed Method The proposed approach in its current implementation is a 2D two-phase level set method. This fact poses certain limitations on its application areas. First, there is the intensity range problem. Although the segmentation results as a stack of 2D slices can be combined and post-processed as a 3D volume, as it was done in [19], the intensity ranges of bias field correction results are not correlated with each other. This makes it difficult to process the bias field correction results. Second, there is a two-phase level set limitation. The proposed algorithm extracts two classes from an image, which is desirable for certain application areas, as breast or ankle segmentation, but will lead to inappropriate results for the data, where more intensity classes are present, such as brain structures. The multi-phase implementation is left for future research. 25

Page 26 of 34

ip t cr

M

an

us

(a) Left: original image, Right: man- (b) Results for N4 (left) and Fuzzy Cual segmentation result, which is cut ap- Means results (right). DSCBreast = 0.9 proximately in the sternum location. and DSCP arenchyma = 0.41.

te

d

(c) Intensity inhomogeneity correction (left) and segmentation results (right) for the proposed method. DSCBreast = 0.963 and DSCP arenchyma = 0.75.

(d) Intensity inhomogeneity correction (left) and segmentation results (right) for the method of Li et al. [26]. DSCBreast = 0.95 and DSCP arenchyma = 0.31.

Figure 9: Initial slice and a groundtruth produced by radiologists and intensity inhomogeneity correction and segmentation results from three automated approaches. N4 does not fully erase the intensity inhomogeneity, the clustering misclassifies the breast tissue (cf. Figure 9b). The method of Li et al. [26] accurately detects the breast volume, but undersegments the parenchymal tissue (cf. Figure 9d). The proposed method is superior than the other approaches.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Third, there is a computational time limitation. Although the presented approach works rather fast on GPU, processing of thousands of images (hundreds of 3D datasets) stays computationally expensive. In current implementation, the processing of a 3D dataset with the resolution of 512 × 512 × 128 takes approximately 30 minutes including data transfer to GPU and result saving. We will investigate further optimization of the presented approach to reduce the computational costs.

26

Page 27 of 34

6. Conclusions and Future Work

te

d

M

an

us

cr

ip t

We present a variational level set method for simultaneous segmentation and intensity inhomogeneity correction. We formulate an energy by applying a generally accepted model for images with intensity inhomogeneities and a local coherent clustering criterion. We extend this criterion to include an explicit manipulation of the bias field component. The novel functional is minimized using the efficient Split Bregman method. It attempts to find the global minimum and is not sensitive to initialization for both contour segmentation and bias field extraction. Our method has been tested on synthetic and real examples and compared to state-of-the-art methods. The proposed approach with the flexible bias field component model produces promising results, which are similar to the results from the other methods or superior in quality, even for images with strong and complicated inhomogeneity gradients and sparse tissue structures. We observe that level set based approaches usually produce smooth and rather accurate contours. However, the results and convergence speed of the level set based procedures depend on the location of initial contours. Moreover, due to the classical gradient descent based minimization, level set methods usually find a local minimum. Our method overcomes these limitations. As future work, we will extend the method to three dimensions and consider multi-phase level sets, since the main limitations of the presented approach are that it operates only in 2D and is able to segment just two classes. Moreover, we also plan to investigate the influence of different versions of the minimization procedure on the results.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Acknowledgment

The authors would like to thank Dr. C. Li, Dr. X. Bresson, Dr. D.J. Kroon, and Dr. T. Goldstein for providing a research source code of the methods. This work was supported by the Ministry for Education, Science and Culture and the European Social Fund (Grant UG 11 035A). MR imaging examinations in SHIP are supported by the Federal Ministry of Education and Research (grant no. 03ZIK012) and a joint grant from Siemens Healthcare, Erlangen, Germany and the Federal State of Mecklenburg West Pomerania.

27

Page 28 of 34

ip t

[1] Ahmed, M. N., Yamany, S. M., Mohamed, N., Farag, A. A., Moriarty, T., 2002. A modified fuzzy c-means algorithm for bias field estimation and segmentation of mri data. IEEE Transactions on Medical Imaging 21 (3), 193–199.

cr

[2] Bezdek, J., Ehrlich, R., Full, W., 1984. Fcm: The fuzzy c-means clustering algorithm. Computers and Geosciences 10 (2), 191–203.

us

[3] Bresson, X., 2009. A short guide on a fast global minimization algorithm for active contour models.

an

[4] Brown, R. W., Cheng, Y.-C. N., Haacke, E. M., Thompson, M. R., Venkatesan, R., 2014. Magnetic resonance imaging: physical principles and sequence design. John Wiley & Sons.

M

[5] Chambolle, A., 2004. An algorithm for total variation minimization and applications. Journal of Mathematical imaging and vision 20 (1-2), 89– 97.

d

[6] Chan, T. F., Esedoglu, S., Nikolova, M., 2006. Algorithms for finding global minimizers of image segmentation and denoising models. SIAM Journal of Applied Mathematics 66 (5), 1632–1648.

te

[7] Chan, T. F., Vese, L. A., 2001. Active contours without edges. IEEE Transactions on Image Processing 10 (2), 266–277. [8] Chen, W., Giger, M. L., 2004. A fuzzy c-means (fcm) based algorithm for intensity inhomogeneity correction and segmentation of mr images. In: Biomedical Imaging: Nano to Macro, 2004. IEEE International Symposium on. IEEE, pp. 1307–1310.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[9] Chua, Z. Y., Zheng, W., Chee, M. W., Zagorodnov, V., 2009. Evaluation of performance metrics for bias field correction in mr brain images. Journal of Magnetic Resonance Imaging 29 (6), 1271–1279.

[10] Comaniciu, D., Meer, P., 2002. Mean shift: A robust approach toward feature space analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (5), 603–619. [11] Cui, W., Wang, Y., Fan, Y., Feng, Y., Lei, T., 2013. Localized fcm clustering with spatial information for medical image segmentation and bias field estimation. Journal of Biomedical Imaging 2013, 13. 28

Page 29 of 34

ip t

[12] Dice, L. R., 1945. Measures of the amount of ecologic association between species. Ecology 26 (3), 297–302.

cr

[13] Goldstein, T., Bresson, X., Osher, S., 2010. Geometric applications of the split bregman method: Segmentation and surface reconstruction. Journal of Scientific Computing 45, 272–293.

us

[14] Goldstein, T., Osher, S., 2009. The split bregman method for l1regularized problems. SIAM Journal of Imaging Sciences 2 (2), 323–343. [15] Gonzalez, R., Woods, R. E., 2002. Digital Image Processing. Prentice Hall International.

M

an

[16] Gubern-M´erida, A., Kallenberg, M., Mann, R. M., Marti, R., Karssemeijer, N., 2015. Breast segmentation and density estimation in breast mri: a fully automatic framework. Biomedical and Health Informatics, IEEE Journal of 19 (1), 349–357.

d

[17] Guillemaud, R., Brady, M., 1997. Estimating the bias field of mr images. IEEE Transactions on Medical Imaging 16 (3), 238–251.

te

[18] Hou, Z., 2006. A review on mr image intensity inhomogeneity correction. International Journal of Biomedical Imaging 1, 1–11. [19] Ivanovska, T., Laqua, R., Wang, L., Liebscher, V., V¨olzke, H., Hegenscheid, K., 2014. A level set based framework for quantitative evaluation of breast tissue density from mri data. PloS one 9 (11), e112709.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[20] Ivanovska, T., Laqua, R., Wang, L., V¨olzke, H., Hegenscheid, K., 2013. Fast implementations of the levelset segmentation method with bias field correction in mr images: Full domain and mask-based versions. In: Pattern Recognition and Image Analysis. Springer, pp. 674–681. [21] Ivanovska, T., Wang, L., Laqua, R., Hegenscheid, K., V¨olzke, H., Liebscher, V., 2013. A fast global variational bias field correction method for mr images. In: Image and Signal Processing and Analysis (ISPA), 2013 8th International Symposium on. IEEE, pp. 667–671. [22] Kimmel, R., 2004. Numerical Geometry of Images. Theory, Algorithms, and Applications. Springer-Verlag New York Inc.

29

Page 30 of 34

ip t

[23] Li, B. N., Chui, C. K., Chang, S., Ong, S. H., 2011. Integrating spatial fuzzy clustering with level set methods for automated medical image segmentation. Computers in Biology and Medicine 41 (1), 1–10.

cr

[24] Li, B. N., Chui, C. K., Chang, S., Ong, S. H., 2012. A new unified level set method for semi-automatic liver tumor segmentation on contrastenhanced ct images. Expert Systems with Applications 39 (10), 9661– 9668.

an

us

[25] Li, C., Gore, J. C., Davatzikos, C., 2014. Multiplicative intrinsic component optimization (mico) for mri bias field estimation and tissue segmentation. Magnetic resonance imaging 32 (7), 913–923.

M

[26] Li, C., Huang, R., Ding, Z., et al., 2011. A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI. IEEE Transactions on Image Processing 20, 2007–2016.

d

[27] Li, C., Kao, C., Gore, J. C., Ding, Z., 2008. Minimization of regionscalable fitting energy for image segmentation. IEEE Transactions on Image Processing 17 (10), 1940–1949.

te

[28] Li, C., Xu, C., Anderson, A. W., Gore, J. C., 2009. Mri tissue classification and bias field estimation based on coherent local intensity clustering: A unified energy minimization framework. In: Information Processing in Medical Imaging (IPMI). Springer, pp. 288–299.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[29] Li, C., Xu, C., Gui, C., Fox, M. D., 2010. Distance regularized level set evolution and its application to image segmentation. IEEE Transactions on Image Processing 19 (12), 3243–3254. [30] Li, X., Li, L., Lu, H., Chen, D., Liang, Z., 2003. Inhomogeneity correction for magnetic resonance images with fuzzy c-mean algorithm. In: Medical Imaging 2003. International Society for Optics and Photonics, pp. 995–1005. [31] Likar, B., Viergever, M. A., Pernus, F., 2001. Retrospective correction of mr intensity inhomogeneity by information minimization. IEEE Transactions on Medical Imaging 20 (12), 1398–1410.

30

Page 31 of 34

ip t

[32] Lin, M., Chan, S., Chen, J., Chang, D., et al., 2011. A new bias field correction method combining n3 and fcm for improved segmentation of breast density on mri. Medical physics 38, 5.

cr

[33] Luebke, D., Harris, M., Govindaraju, N., Lefohn, A., Houston, M., et al., 2006. GPGPU: general-purpose computation on graphics hardware. In: Proceedings of the 2006 ACM/IEEE conference on Supercomputing. ACM, pp. 208–208.

an

us

[34] Makarau, A., Huisman, H., Mus, R., Zijp, M., Karssemeijer, N., 2010. Breast mri intensity non-uniformity correction using mean shift. In: Proc. SPIE 7624, Medical Imaging 2010: Computer-Aided Diagnosis.

M

[35] Moreno, J. C., Prasath, V. S., Proen¸ca, H., Palaniappan, K., 2014. Fast and globally convex multiphase active contours for brain mri segmentation. Computer Vision and Image Understanding 125, 237–250. [36] Nvidia, 2008. Nvidia cuda: Programming guide.

d

[37] Osher, S., Paragios, N. (Eds.), 2003. Geometric Level Set Methods in Imaging, Vision, and Graphics. Springer.

te

[38] Pham, D. L., Prince, J. L., 1999. An adaptive fuzzy c-means algorithm for image segmentation in the presence of intensity inhomogeneities. Pattern Recognition Letters 20 (1), 57–68.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[39] Press, W. H., 2007. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press. [40] Razavi, M., Wang, L., Gubern-M´erida, A., Ivanovska, T., Laue, H., Karssemeijer, N., Hahn, H. K., 2015. Towards accurate segmentation of fibroglandular tissue in breast mri using fuzzy c-means and skin-folds removal. In: Image Analysis and ProcessingICIAP 2015. Springer, pp. 528–536. [41] Rudin, L., Osher, S., Fatemi, E., 1992. Nonlinear rotal variation based noise removal algorithms. Physica D 60, 259–268. [42] Sapiro, G., 2009. Geometric Partial Differential Equations and Image Analysis. Cambridge University Press.

31

Page 32 of 34

ip t

[43] Sethian, J. A., 1999. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. Cambridge University Press.

cr

[44] Sled, J. G., Zijdenbos, A. P., Evans, A. C., 1998. A nonparametric method for automatic correction of intensity nonuniformity in mri data. IEEE Transactions on Medical Imaging 17 (1), 87–97.

an

us

[45] Tsai, A., Yezzi Jr, A., Willsky, A. S., 2001. Curve evolution implementation of the mumford-shah functional for image segmentation, denoising, interpolation, and magnification. Image Processing, IEEE Transactions on 10 (8), 1169–1186.

M

[46] Tustison, N. J., Avants, B. B., Cook, P. A., Zheng, Y., Egan, A., Yushkevich, P. A., Gee, J. C., 2010. N4itk: improved n3 bias correction. Medical Imaging, IEEE Transactions on 29 (6), 1310–1320.

d

[47] Vovk, U., Pernus, F., Likar, B., 2007. A review of methods for correction of intensity inhomogeneity in mri. IEEE Transactions on Medical Imaging 26 (3), 405–421.

te

[48] Wang, L., Pan, C., 2014. Image-guided regularization level set evolution for mr image segmentation and bias field correction. Magnetic Resonance Imaging 32 (1), 71–83. [49] Wang, L., Shi, F., Li, G., Gao, Y., Lin, W., Gilmore, J. H., Shen, D., January 2014. Segmentation of neonatal brain mr images using patchdriven level sets. NeuroImage 84, 141158.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[50] Wells III, W. M., Grimson, W. E. L., Kikinis, R., Jolesz, F. A., 1996. Adaptive segmentation of mri data. IEEE Transactions on Medical Imaging 15 (4), 429–442. [51] Yang, Y., Wu, B., 2012. Convex image segmentation model based on local and global intensity fitting energy and split bregman method. Journal of Applied Mathematics 2012.

[52] Yang, Y., Wu, B., 2012. Fast multiphase image segmentation model for images with inhomogeneity. Journal of Electronic Imaging 21 (1).

32

Page 33 of 34

ip t

[53] Zhan, T., Zhang, J., Xiao, L., Chen, Y., Wei, Z., 2013. An improved variational level set method for mr image segmentation and bias field correction. Magnetic Resonance Imaging 31 (3), 439 – 447.

cr

[54] Zhang, K., Zhang, L., Lam, K., Zhang, D., 2013. A local active contour model for image segmentation with intensity inhomogeneity. arXiv preprint arXiv:1305.7053.

te

d

M

an

us

[55] Zhang, K., Zhang, L., Zhang, S., 2010. A variational multiphase level set approach to simultaneous segmentation and bias correction. In: 17th IEEE International Conference on Image Processing (ICIP). IEEE, pp. 4105–4108.

Ac ce p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

33

Page 34 of 34