Stabilised bias field: segmentation with intensity inhomogeneity

Abstract

Automatic segmentation in the variational framework is a challenging task within the field of imaging sciences. Achieving robustness is a major problem, particularly for images with high levels of intensity inhomogeneity. The two-phase piecewise-constant case of the Mumford-Shah formulation is most suitable for images with simple and homogeneous features where the intensity variation is limited. However, it has been applied to many different types of synthetic and real images after some adjustments to the formulation. Recent work has incorporated bias field estimation to allow for intensity inhomogeneity, with great success in terms of segmentation quality. However, the framework and assumptions involved lead to inconsistencies in the method that can adversely affect results. In this paper we address the task of generalising the piecewise-constant formulation, to approximate minimisers of the original Mumford-Shah formulation. We first review existing methods for treating inhomogeneity, and demonstrate the inconsistencies with the bias field estimation framework. We propose a modified variational model to account for these problems by introducing an additional constraint, and detail how the exact minimiser can be approximated in the context of this new formulation. We extend this concept to selective segmentation with the introduction of a distance selection term. These models are minimised with convex relaxation methods, where the global minimiser can be found for a fixed fitting term. Finally, we present numerical results that demonstrate an improvement to existing methods in terms of reliability and parameter dependence, and results for selective segmentation in the case of intensity inhomogeneity.

Introduction

The task of partitioning an image into multiple regions each sharing certain characteristics – such as texture, intensity, shape or colour – is called segmentation, and is an important aspect of image processing. Given an image z(x) in a bounded domain $\Omega \subset ℝ2$ , we look for an edge Γ that partitions Ω into regions $\{\Omega i,i=1,2,\dots ,N\}$ in $\Omega \setminus \Gamma$ . The particular case of selecting objects in an image based on intensity has been widely studied over the last 20 years,^1–3 and is particularly challenging in cases of intensity inhomogeneity.

We consider the variational approach to this problem, and in particular region-based active contour models (ACMs). These are based on the introduction of the minimisation of the Mumford-Shah functional, ³ given by

$E\mathrm{MS}(\Gamma ,w)=\left|\Gamma \right|+\lambda \int \Omega (z-w)2\text{d}x+\mu \int \Omega \setminus \Gamma |\nabla w|\text{d}x,$

(1)

where $\mu ,\lambda >\text{F}0$ are weighting parameters and $\left|\Gamma \right|$ denotes the length of the edge curve Γ, the boundary between regions Ω _i . Mumford and Shah ³ demonstrated that, theoretically, the existence and regularity of minimisers of this functional can be achieved, and Tsai et al. ⁴ and Vese and Chan ⁵ used the variational level set method of Osher et al. ⁶ and Zhao et al. ⁷ to minimise (equation (1)). Chan and Vese ² introduced a functional that was a particular case of equation (1), the two-phase piecewise-constant example, i.e. $|\nabla w|=0$ and N = 2 in equation (1).

$E\mathrm{CV}(\Gamma ,c1,c2)=\left|\Gamma \right|+\lambda (\int \Omega 1(z-c1)2\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\int \Omega 2(z-c2)2\text{d}x)$

(2)

The Chan-Vese model (CV) has been widely used in segmentation applications since its introduction in 2001. Its framework has been generalised by the introduction of new fitting terms to incorporate extensive intensity inhomogeneity, such as Li et al.^8,9 who introduced a region scalable fitting energy and local cluster method. Jung et al.^10,11 introduced a nonlocal ACM utilising distance functions. Brox and Cremers ¹² and Lanktona and Tannenbaum ¹³ introduced new local models, incorporating Gaussian kernel functions. Recent work related to this area includes the work of Ali et al., ¹⁴ who form fitting terms using multiplicative and difference image data, and L₀ regularisation for simultaneous bias correction and segmentation by Duan et al. ¹⁵

A drawback of the CV approach is a lack of convexity in the level-set-based minimisation. The method is based on using the Heaviside function to represent the two regions Ω₁ and Ω₂ with respect to a level set function, φ, and computing the Euler-Lagrange equation. The solution of the corresponding partial differential equation is often a local minimum, which reduces the reliability of the results. The introduction of the idea of convex relaxation by Chan et al. ¹ demonstrated that a global minimum of equation (2) with respect to Γ, and consequently the regions Ω₁ and Ω₂, can be found. The idea is to represent the two regions by an indicator function, and relax the constraint such that both the functional and the constraint sets are convex. Further work has been done by Bresson et al. ¹⁶ and Chambolle et al. ¹⁷ and extended to the multiphase framework, i.e. N > 2, by Lellmann et al., ¹⁸ Bae et al. ¹⁹ and Gu et al. ²⁰

The recent work of Chen et al.,²³ combines the idea of convex relaxation and segmentation with intensity inhomogneity, and our work focuses on aiming to improve their formulation. The ‘true’ image data are formulated^22,23 as

$T=\sum _{i}ci\chi i,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}i=1,2,\dots ,N,$

(3)

where c_i are intensity constants, and χ_i are characteristic functions of the regions Ω _i . It is based on the idea that the image can be modelled as

$z(x)=B(x)T+\eta ,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}x\in \Omega i,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}i=1,2,\dots ,N,$

(4)

where η is additive noise. Here, as with CV, we consider the two-phase case, i.e. N = 2. Chen et al. ²¹ aim to estimate the bias field B and recover the ‘true’ image T. However, a lack of convergence of the bias field and the intensity constants means that recovering an accurate T is not possible. We propose introducing a constraint on these variables to correct this inconsistency.

The paper is organised as follows. In “VMS model” we detail the Variant Mumford-Shah (VMS) model, ²¹ briefly introduced above, and discuss the problems with recovering the ‘true’ image, in particular the lack of convergence of c₁ and c₂ due to the formulation. In “Stabilised bias field” we detail the introduction of a constraint to the work of Chen et al. ²¹ in order to automatically establish feasible intensity constants and ensure the convergence of all variables being minimised. We discuss how this alters the minimisation of the bias field, how the functional is iteratively minimised, and details of the numerical implementation. We also highlight the link the proposed method provides between Mumford-Shah ³ and Chan-Vese. ² In “Results” we include experimental results that measure the accuracy of the proposed method compared to VMS by using the Tanimoto Coefficient, ²⁴ and demonstrate the convergence of the intensity constants for examples used in Chen et al. ²¹ We extend this idea to selective segmentation in “Selective segmentation with SBF”. We consider incorporating a distance selection term from a recent selection model, ²⁵ and include experimental results for one challenging case. We discuss the benefits of the proposed method in “Conclusions”.

VMS model

The VMS model by Chen et al. ²¹ is formulated as follows

$E\mathrm{VMS}(\Gamma ,c1,c2,B)=\left|\Gamma \right|+\lambda \int \Omega ((z-\mathrm{Bc}1)2\chi 1\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+(z-\mathrm{Bc}2)2\chi 2)\text{d}x+\mu \int \Omega |\nabla B|2\text{d}x,$

(5)

where λ and μ are weighting parameters. The idea is that the intensity constants represent the ‘true’ image, and the bias field B varies such that their combination gives a piecewise-smooth approximation of z (with respect to the Mumford-Shah ³ formulation (equation (1.1))), given by

$w\mathrm{VMS}=\mathrm{Bc}1\chi 1+\mathrm{Bc}2\chi 2.$

(6)

The functional equation (5) is minimised iteratively by the following steps. Step (1): For fixed characteristic functions χ₁ and χ₂, and intensity constants c₁ and c₂, minimise equation (5) with respect to bias field estimator B. Based on the work of Nielsen et al. ²⁶ and Brox and Cremers, ¹² the exact minimiser can be well approximated. Step (2): For fixed characteristic functions χ₁ and χ₂, and bias field estimator B, minimise equation (5) with respect to intensity constants c₁ and c₂. These can be computed precisely. Step (3): For fixed intensity constants c₁ and c₂, and bias field estimator B, minimise equation (5) with respect to $\chi 1,\chi 2$ . Based on the work of Chan et al., ¹ minimising two characteristic functions can be achieved by using an indicator function, u(x), which is allowed to take intermediate values.

Convergence behaviour of VMS

In Figure 1, we demonstrate a result for VMS that is also used in Chen et al., ²¹ and is of comparable quality. However, the question remains: based on the image model described above (equation (3)), what is the ‘true’ image? Whilst the joint minimisation of equation (5) with respect to c₁, c₂, B and Γ is nonconvex, and therefore we cannot determine the correct c₁ and c₂ precisely, there is a problem with the current framework, which we will now discuss. In Figure 1 (after 1000 iterations), we show that the values of the intensity constants continually rise, such that $c1=9.1\times 104$ and $c2=6.4\times 104$ . The convergence of w^VMS (equation (6)) comes from the reduction in scale of B. To demonstrate this, after the same number of iterations $\left|\right|B\left|\right|F=4.7\times 10-4$ (where $\left|\right|·\left|\right|F$ is the Frobenius norm). This motivates our proposal for modifying VMS, in the form of an additional constraint, which can automatically control the scale of c₁, c₂ and B.

Figure 1.

Convergence behaviour. The first row, from left to right, shows the lack of convergence for the intensity constants, giving $c1=9.1\times 104$ and $c2=6.4\times 104$ , and the scale of the bias field B(x). The bottom row, from left to right, shows the image z(x), and the computed contour $\Gamma *$ , and the progression of the energies E^B and E^F (equation (7)). (a) c₁ Progression, (b) c₂ Progression, (c) Bias Field, B(x), (d) z(x) ∈ [0; 1]; Γ*, (e) E^B and (f) E^f.

To explain this phenomenon let us examine the VMS functional equation (5). The smoothness penalty included, which we denote E^B, is similar to the penalty enforced in the Mumford-Shah functional (equation (1)), except that it applies throughout the domain. We denote the fitting energy E^F, and it is again similar to the Mumford-Shah fitting energy

$EB=\int \Omega |\nabla B|2\text{d}x,EF=\int \Omega ((z-\mathrm{Bc}1)2\chi 1+(z-\mathrm{Bc}2)2\chi 2)\text{d}x.$

(7)

However, crucially, the Mumford-Shah fitting term only involves one variable, w. The VMS fitting term involves the products Bc₁ and Bc₂. This means that a change in one variable does not necessarily alter the energy, as long as the other variable changes accordingly. In practice that means that the minimum of the VMS functional is attained when $EB\to 0$ , despite the convergence of E^F. This is due to the lack of convergence of the intensity constants c₁ and c₂, but this contradicts the assumptions of recovering the ‘true’ image (equation (3)) discussed in “Introduction”. This is demonstrated in Figure 1.

Stabilised bias field

VMS produces a piecewise-smooth approximation of the image (in the Mumford-Shah sense ³ ), given by w^VMS (equation (6)). However, it does not give values for c₁ and c₂ that are consistent with the observed image. It is possible to manually rescale these without changing w^VMS, but this is not a sensible approach as these values are unknown by definition. The immediate question is: is it possible to incorporate constraints into the formulation in a reliable way, i.e. can we use information in the image to automatically restrict the scale of B,c₁ and c₂? There are two obvious approaches. The first is to constrain the values of c₁ or c₂. The situation when the optimal intensity constants are not known a priori has been studied by Brown et al. ²⁷ in the piecewise-constant case, but not in cases of intensity inhomogeneity. It is worth considering how this method could be incorporated in the presence of a bias field function; however, we do not discuss this here. The second is to control the scale of the bias field, B. We therefore consider how to introduce a constraint in such a way that it provides a link between the piecewise-constant and piecewsie-smooth approximations of z that are consistent with the image, which we will return to later.

With VMS, B is encouraged to be close to 0, which leads to the lack of convergence for c₁ and c₂. To prevent this we propose a new model we call stabilised bias field (SBF), with the introduction of an additional constraint that encourages B to be close to a positive constant. However, this alters the minimisation step for the bias field from Chen et al. ²¹ We now consider how to obtain this with the addition of this constraint. To distinguish between the two methods we refer to the bias field in SBF as $\tilde{B}$ . The new formulation is given as follows

$E\mathrm{SBF}(\Gamma ,c1,c2,\tilde{B})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=\left|\Gamma \right|+\lambda \int \Omega ((z-\tilde{B}c1)2\chi 1+(z-\tilde{B}c2)2\chi 2)\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\mu \int \Omega |\nabla \tilde{B}|2\text{d}x+\gamma \int \Omega (\tilde{B}-\mathbf{s})2\text{d}x,$

(8)

where $\mathbf{s},\gamma$ are positive parameters. We intend to use the framework of VMS to approximate the exact minimiser of equation (8) for $\tilde{B}$ . With this in mind, the previous formulation is equivalent to

$E\mathrm{SBF}(\Gamma ,c1,c2,\tilde{B})=\lambda \int \Omega ([(z-\tilde{B}c1)2+\frac{\gamma }{\lambda }(B-\mathbf{s})2]\chi 1\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+[(z-\tilde{B}c2)2+\frac{\gamma }{\lambda }(\tilde{B}-\mathbf{s})2]\chi 2)\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\left|\Gamma \right|+\mu \int \Omega |\nabla \tilde{B}|2\text{d}x,$

where the new constraint has been incorporated into the fitting term. We can reformulate this as follows, first looking at the χ₁ term

$[(z-\tilde{B}c1)2+\frac{\gamma 2}{\lambda }(\tilde{B}-\mathbf{s})2]\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=\tilde{B}2(c_1^2+\tilde{\gamma })-2\tilde{B}(c1z+\tilde{\gamma }\mathbf{s})+(z2+\tilde{\gamma }\mathbf{s}2)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=(c_1^2+\tilde{\gamma })[\tilde{B}-\frac{c1z+\tilde{\gamma }\mathbf{s}}{c_1^2+\tilde{\gamma }}]2+f1(z,c1,\mathbf{s},\tilde{\gamma })\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=[\frac{c1z+\tilde{\gamma }\mathbf{s}}{\sqrt{c_1^2+\tilde{\gamma }}}-\tilde{B}\sqrt{c_1^2+\tilde{\gamma }}]2+f1(z,c1,\mathbf{s},\tilde{\gamma }),$

where $\tilde{\gamma }=\frac{\gamma }{\lambda }$ and $f1(z,c1,\mathbf{s},\tilde{\gamma })=\frac{z2+\tilde{\gamma }\mathbf{s}2}{c_1^2+\tilde{\gamma }}$ . In a similar way, for the χ₂ term

$[(z-\tilde{B}c2)2+\frac{\gamma 2}{\lambda }(\tilde{B}-\mathbf{s})2]=[\frac{c2z+\tilde{\gamma }\mathbf{s}}{\sqrt{c_2^2+\tilde{\gamma }}}-\tilde{B}\sqrt{c_2^2+\tilde{\gamma }}]2\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+f2(z,c2,\mathbf{s},\tilde{\gamma }),$

where $f2(z,c2,\mathbf{s},\tilde{\gamma })=\frac{z2+\tilde{\gamma }\mathbf{s}2}{c_2^2+\tilde{\gamma }}$ . Therefore equation (8) is equivalent to

$E\mathrm{SBF}(\Gamma ,c1,c2,\tilde{B})=\lambda \int \Omega ([\frac{c1z+\tilde{\gamma }\mathbf{s}}{\sqrt{c_1^2+\tilde{\gamma }}}-\tilde{B}\sqrt{c_1^2+\tilde{\gamma }}]2\chi 1\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+[\frac{c2z+\tilde{\gamma }\mathbf{s}}{\sqrt{c_2^2+\tilde{\gamma }}}-\tilde{B}\sqrt{c_2^2+\tilde{\gamma }}]2\chi 2)\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\int \Omega f1(z,c1,\mathbf{s},\tilde{\gamma })\chi 1\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\int \Omega f2(z,c2,\mathbf{s},\tilde{\gamma })\chi 2\text{d}x+\left|\Gamma \right|\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\mu \int \Omega |\nabla \tilde{B}|2\text{d}x.$

Minimising $E\mathrm{SBF}(\Gamma ,c1,c2,\tilde{B})$ with respect to $\tilde{B}$ is given by

$min\tilde{B}\{\lambda \int \Omega ([z1-\tilde{B}\tilde{c}1]2\chi 1+[z2-\tilde{B}\tilde{c}2]2\chi 2)\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+\mu \int \Omega |\nabla \tilde{B}|2\text{d}x\}$

(9)

since $f1(z,c1,\mathbf{s},\tilde{\gamma }),f2(z,c2,\mathbf{s},\tilde{\gamma })$ , and $\left|\Gamma \right|$ are not dependent on $\tilde{B}$ . Here

$z1=\frac{c1z+\tilde{\gamma }\mathbf{s}}{\sqrt{c_1^2+\tilde{\gamma }}},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\tilde{c}1=\sqrt{c_1^2+\tilde{\gamma }},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}z2=\frac{c2z+\tilde{\gamma }\mathbf{s}}{\sqrt{c_2^2+\tilde{\gamma }}},\text{and}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\tilde{c}2=\sqrt{c_2^2+\tilde{\gamma }}.$

In the same way as VMS, ²¹ we can approximate the exact minimiser of equation (9) with a Gaussian G

$\tilde{B}=\frac{\tilde{c}1z1\chi 1+\tilde{c}2z2\chi 2}{{\tilde{c}}_1^2\chi 1+{\tilde{c}}_2^2\chi 2}*G.$

(10)

Relationship to Chan-Vese and Mumford-Shah

We now discuss how the proposed model relates to the two important works discussed in the Introduction. The SBF functional is given as

$E\mathrm{SBF}(\Gamma ,c1,c2,\tilde{B})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=\left|\Gamma \right|+\lambda \int \Omega ((z-\tilde{B}c1)2\chi 1+(z-\tilde{B}c2)2\chi 2)\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\mu \int \Omega |\nabla \tilde{B}|2\mathrm{dx}+\gamma \int \Omega (\tilde{B}-\mathbf{s})2\text{d}x.$

(11)

It relates to Mumford-Shah in the same sense that VMS does. That is, we can compute a piecewise-smooth approximation of the image

$w\mathrm{SBF}=\tilde{B}c1\chi 1+\tilde{B}c2\chi 2,$

(12)

except that the values computed correspond to the observed image $z\in [0,1]$ , and the variables converge reliably. However, it also relates to the Chan-Vese functional. If s = 1 and $\gamma \to \infty$ we have the CV formulation (equation (2))

$E\mathrm{CV}(\Gamma ,c1,c2)=\left|\Gamma \right|+\lambda \int \Omega ((z-c1)2\chi 1+(z-c2)2\chi 2)\text{d}x.$

Iterative minimisation of SBF formulation

We now detail how to minimise the functional equation (8), in line with the method of Chen et al., ²¹ in order to effectively compare our proposed method against VMS. The SBF model is given as follows

$min\Gamma ,c1,c2,\tilde{B}\{E\mathrm{SBF}(\Gamma ,c1,c2,\tilde{B})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=\left|\Gamma \right|+\lambda \int \Omega ((z-\tilde{B}c1)2\chi 1+(z-\tilde{B}c2)2\chi 2)\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\mu \int \Omega |\nabla \tilde{B}|2\mathrm{dx}+\gamma \int \Omega (\tilde{B}-\mathbf{s})2\text{d}x\}.$

This is minimised iteratively (e.g. the iterative process method, Li et al. ⁹ ) by the following steps:

For fixed characteristic functions χ₁ and χ₂, and intensity constants c₁ and c₂, minimise equation (11) with respect to bias field estimator $\tilde{B}$ .

For fixed characteristic functions χ₁ and χ₂, and bias field estimator $\tilde{B}$ , minimise equation (11) with respect to intensity constants c₁ and c₂.

For fixed intensity constants c₁ and c₂, and bias field estimator $\tilde{B}$ , minimise equation (11) with respect to characteristic functions χ₁ and χ₂.

We provide a summary of how each step is minimised in the following:

Step (1): detailed in the previous section. It can be approximated, according to the work of Nielsen et al. ²⁶ and Brox and Cremers ¹² and as discussed by Chen et al., ²¹ by

$\tilde{B}=\frac{\tilde{c}1z1\chi 1+\tilde{c}2z2\chi 2}{{\tilde{c}}_1^2\chi 1+{\tilde{c}}_2^2\chi 2}*G.$

Step (2): minimising with respect to c₁ and c₂ gives

$c1=\frac{\int \Omega z(x)B(x)\chi 1\text{d}x}{\int \Omega B2(x)\chi 1\text{d}x},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}c2=\frac{\int \Omega z(x)B(x)\chi 2\text{d}x}{\int \Omega B2(x)\chi 2\text{d}x}.$

(13)

Step (3): achieved by the following minimisation

$min\chi 1,\chi 2\{\left|\Gamma \right|+\lambda \int \Omega ((z-\tilde{B}c1)2\chi 1+(z-\tilde{B}c2)2\chi 2)\text{d}x\}$

(14)

Minimising two characteristic functions can be achieved by using an indicator function, u(x), which is allowed to take intermediate values. This is based on the work of Chan et al. ¹

$min0\le u\le 1\{\int \Omega |\nabla u(x)|\text{d}x+\lambda \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\times \int \Omega ((z-\tilde{B}c1)2-(z-\tilde{B}c2)2)u(x)\text{d}x\},$

(15)

where the length term $\left|\Gamma \right|$ in equation (14) has been replaced by the total variation (TV) of u. The dual formulation method of Chambolle ²⁸ consists of introducing a new variable v and alternating between minimising u and v. Alternative minimisation methods for this type of problem include Goldstein et al. ²⁹ and an improved additive operator splitting method ²⁵ that is applicable to problems of this type. By splitting in this way, the minimisation of u concentrates on the TV term and the minimisation of v satisfies the fitting and constraint requirements

$minu,v\{\int \Omega |\nabla u(x)|\text{d}x+\frac{1}{2\theta }\int \Omega (u(x)-v(x))2\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\int \Omega \lambda r(x)v(x)+\alpha \psi (v)\text{d}x\},$

(16)

where $r(x)=(z-\tilde{B}c1)2-(z-\tilde{B}c2)2$ , and $\psi (v)=max\{0,2|v-\frac{1}{2}|-1\}$ . Two parameters are introduced here: θ > 0 is a small parameter and $\alpha >\frac{\lambda }{2}\left|\right|r(x)\left|\right|L\infty (\Omega )$ ensures the constraints on the indicator function u(x) in equation (15) are met. The minimisation of u and v, with fixed $c1,c2,\tilde{B}$ , can be achieved iteratively by the following steps:

• Step (3a): With fixed v

  $minu\{\int \Omega |\nabla u(x)|\text{d}x+\frac{1}{2\theta }\int \Omega (u(x)-v(x))2\text{d}x\}$
  which can be solved by ²⁸
  $u(x)=v(x)-\theta \nabla ·\rho (x),$

  (17)
  where $\rho =(\rho 1,\rho 2)$ is the solution of
  $\nabla (\theta \nabla ·\rho -v)-|\nabla (\theta \nabla ·\rho -v)|\rho =0,$

  (18)
  which can be solved by a fixed point method ρ⁰= 0 and
  $\rho n+1=\frac{\rho n+\tau \nabla (\nabla ·\rho n-v/\theta )}{1+\tau |\nabla \rho n-v/\theta |}.$

  (19)

• Step (3b): With fixed u

  $minv\{\frac{1}{2\theta }\int \Omega (u(x)-v(x))2\text{d}x+\int \Omega \lambda r(x)v(x)+\alpha \psi (v)\text{d}x\},$
  given, based on the work of Bresson et al., ¹⁶ by
  $v(x)=min\{max\{u(x)-\theta \lambda r(x),0\},1\}.$

  (20)

Numerical implementation

We now provide details of implementing the three steps above. We follow the work of Chen at al., ²¹ who use slight variations on the formulation, in order to be consistent and ensure a fair comparison between VMS and SBF. The intensity constants are computed using smooth region descriptors $H_{\varsigma }^{(1)}$ and $H_{\varsigma }^{(2)}=1-H_{\varsigma }^{(1)}$ instead of characteristic functions χ₁ and χ₂, respectively. This descriptor is defined as follows

$H_{\varsigma }^{(1)}(\varphi (x))=\frac{1}{2}(1+\frac{2}{\pi }\arctan (\varphi *G\varsigma )),\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}x\in \Omega .$

The variable $\varphi (x)$ is given by

$\varphi (x)=\{c,\text{for}x\in \Omega :u(x)>\epsilon ,-c,\text{for}x\in \Omega :u(x)\le \epsilon ,$

where $\epsilon \in [0,1],c=2$ and $\varsigma =1$ . This adjusts the computation of the intensity constants (equation (13)) to

$c1=\frac{\int \Omega z(x)\tilde{B}(x)H_{\varsigma }^{(1)}\text{d}x}{\int \Omega \tilde{B}2(x)H_{\varsigma }^{(1)}\text{d}x},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}c2=\frac{\int \Omega z(x)\tilde{B}(x)H_{\varsigma }^{(2)}\text{d}x}{\int \Omega \tilde{B}2(x)H_{\varsigma }^{(2)}\text{d}x}.$

(21)

The Gaussian kernel, $G\sigma$ , is truncated as a $\varrho \times \varrho$ mask, where ϱ is the smallest odd number greater than $4\sigma +1$ (σ is the standard deviation of the Gaussian kernel ⁹ ). Other parameters mentioned in “Iterative minimisation of SBF formulation” are set as follows $\tau =1/8,\theta =1/3,\epsilon =1/2$ .

The primary motivation of this model is to have convergence of the intensity constants c₁ and c₂. The value of s determines the size of these values; as $\mathbf{s}\to 0,c1,c2\to \infty$ , as discussed in “Convergence behaviour of VMS”. For consistency, it is desirable that $c1,c2\in [0,1]$ given $z\in [0,1]$ . With that in mind, a natural selection is $\mathbf{s}=1$ given that the intensity constants are then related to the average value of z(x) inside and outside the contour. Additionally, SBF is then clearly related to Chan-Vese ² as detailed in “Relationship to Chan-Vese and Mumford-Shah”. The choice of σ is related to each example as with VMS. The choice of λ is also related to each example in the same way as other segmentation problems of a similar type.^1,2

Algorithm 1 Stabilised Bias Field: ${\Omega }_1^{*}\to \mathrm{SBF}(z,\mathrm{maxit},\delta ,\lambda ,\sigma ,\mathbf{s},\gamma )$

1: Initialise $u(0)$ , estimate $c_1^{(1)},c_2^{(1)},\tilde{B}(1)$ .

2: for $\ell \to 1:\mathrm{maxit}$ do

3: if $\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}mod(\ell ,2)=0$ then

4:  Calculate $\tilde{B}(\ell )$ using equation (10).

5:  Calculate $c_1^{(\ell )},c_2^{(\ell )}$ using equation (21).

6:   $\tilde{B}(\ell +1)$  =  $\tilde{B}(\ell ),c_1^{(\ell +1)}=c_1^{(\ell )},c_2^{(\ell +1)}=c_2^{(\ell )}$ .

7: end if

8: while $max\{||u(\ell )-u(\ell -1)||,||v(\ell )-v(\ell -1)||\}>\delta$ do

9:  Find $\rho (\ell )$ , solution of equation (18), using equation (19).

10:  Get $u(\ell )$ using equation (17), and $v(\ell )$ using equation (20).

11: end while

12: end for

13: $\tilde{B}*=\tilde{B}(\ell ),c_1^{*}=c_1^{(\ell )},c_2^{*}=c_2^{(\ell )},{\Omega }_1^{*}=u(\ell )>\epsilon ,\Gamma *=\partial {\Omega }_1^{*}$ .

Results

This section is in two parts. First we will test SBF using images from the VMS tests in Chen et al., ²¹ intending to show that the proposed method retains the segmentation quality of VMS, whilst demonstrating the convergence of the intensity constants. Another aspect of the success of SBF is what $c_1^{*}$ and $c_2^{*}$ are; we can check whether the computed values are feasible, i.e. $c1,c2\in [0,1]$ , whilst maintaining the quality of the segmentation. Secondly, we investigate other advantages of SBF over VMS. In particular, we look at the segmentation accuracy depending on the fitting parameter λ, and how the piecewise-smooth approximations w^SBF and w^VMS compare for a model example.

Set 1: Convergence behaviour

We test four examples (Images 1–4) in Figure 2, all used in Chen et al. ²¹ In Figure 3, we present the results for each case. We set $\mathbf{s}=1$ , and vary the constraint parameter for each case. Given this value, for an 128 × 128 image, we expect $\left|\right|\tilde{B}\left|\right|F\approx 128$ with the addition of the constraint. It is worth noting that in some sense the value of s is arbitrary as the constraint is also dependent on γ. Experimentally, similar results as presented below can be attained for different values of s, if γ is adjusted accordingly. However, setting $\mathbf{s}=1$ and varying γ is the most intuitive approach to take.

Figure 2.

Images tested with SBF and compared to results of Chen et al. ²¹

Figure 3.

Set 1 Results. SBF convergence behaviour Rows 1–4 are for Images 1–4, respectively. From left to right: z(x) and $\Gamma *$ computed with SBF, bias field $\tilde{B}(x)$ , and the progression of c₁ values (vertical axis) against iterations (horizontal axis.) Similar behaviour for c₂ values is also observed.

For Image 1 $\gamma =0.1$ , and the intensity constants converge to $c_1^{*}=0.352,c_2^{*}=0.240$ . The bias field, $\tilde{B}$ , also converges and we compute $\left|\right|\tilde{B}\left|\right|F=129.8$ . This result is of a similar quality to VMS, shown in Figure 1 of Chen et al. ²¹ For Image 2 $\gamma =0.2$ , and the intensity constants converge to $c_1^{*}=0.807,c_2^{*}=0.660$ . The bias field, $\tilde{B}$ , also converges and we compute $\left|\right|\tilde{B}\left|\right|F=67.0$ . However, row 2 of Figure 3 demonstrates that the convergence of c₁ and c₂ is quite slow, taking over 500 iterations which is much more than for the convergence of w_VMS in VMS. Image 3 is an ultrasound image, containing intensity variation in the background. For this example $\gamma =0.1$ , and the intensity constants converge to $c_1^{*}=0.446,c_2^{*}=0.383$ . The bias field, $\tilde{B}$ , also converges and we compute $\left|\right|\tilde{B}\left|\right|F=120.9$ . Image 4 is another example of vessel segmentation, where intensity varies smoothly throughout the vessel. For this example $\gamma =0.1$ , and the intensity constants converge to $c_1^{*}=0.352,c_2^{*}=0.239$ . The bias field, B, also converges and we compute $\left|\right|\tilde{B}\left|\right|F=122.2$ . We see results of comparable quality to rows 3 and 4 in Figure 3, in Figures 9 and 4 of Chen et al., ²¹ respectively, except that there is no convergence for c₁ and c₂. This demonstrates a clear improvement for these examples. Despite the slow convergence of SBF in the case of Image 2, we have fast convergence in the other examples, meaning the additional constraint generally does not slow down the computation of a solution. Also, the results are not sensitive to the constraint parameter, γ. For Image 2, it was adjusted to 0.2, but for all other cases it was set at 0.1, and for all examples, $c_1^{*},c_2^{*}\in [0,1]$ , showing that the method produces results consistent with the image.

Set 2: Comparison to VMS

With Result Set 1, we have successfully demonstrated that SBF achieves the intended goal: the convergence of the intensity constants within a feasible range, and the computation of a stabilised bias field. However, we now intend to examine the success of the proposed method in another way: how does the method affect the accuracy of the final segmentation. With this in mind, we can quantifiably measure the solution of each model (VMS and SBF) against this using the Tanimoto Coefficient ²⁴

$\mathrm{TC}=\frac{N(\mathrm{GT}\cap {\Omega }_1^{*})}{N(\mathrm{GT}\cup {\Omega }_1^{*})},$

(22)

where $N(·)$ is the number of pixels in the enclosed region, GT is the ground truth, and ${\Omega }_1^{*}$ is the result computed with VMS or SBF. However, without the ground truth data for Images 1–4 (Figure 2) we cannot measure this for the examples used by Chen et al. ²¹ Instead, we test one model example shown in Figure 4 where the ground truth is known precisely. We observe two things with this example. Firstly, how w^VMS and w^SBF compare visually with each other. In Figure 4 we can see that around $\Gamma *$ there are significant differences between the two approximations. SBF appears to produce a sharper approximation of the image, dealing with the discontinuity in the intensity more effectively. Secondly, we have tested each model with a large range of the fitting parameter λ. Our results are promising in the sense that an optimal result can be computed for a wider range of λ with SBF over VMS. It is worth noting that we have not observed such pronounced results with Images 1–4, although an advantage is still present.

Figure 4.

Set 2 Results. SBF compared to VMS (a) Successful segmentation of the image, z(x), given by contour $\Gamma *$ . (b) Computed stabilised bias field, $\tilde{B}(x)$ . (c) Convergence of c₁ values (50 iterations). Similar behaviour for c₂ values is also observed. (d) Piecewise-smooth approximation of z(x) with SBF. (e) Piecewise-smooth approximation of z(x) with VMS. (f) Difference between SBF and VMS approximations, demonstrating significant differences around $\Gamma *$ . (g) The TC measure for VMS (dotted red) and SBF (blue), demonstrating the segmentation quality falls away for VMS with large values of λ. (h) The TC measure for $\lambda \in [5,100]$ shows that an optimal $\Gamma *$ can be computed for a larger range of λ with SBF than VMS.

Selective segmentation with SBF

Selective segmentation is the task of extracting one particular object of interest, from a foreground with similar characteristics. We now consider the problem of selecting objects in images that contain significant intensity inhomogeneity, which is beyond recent work on selective segmentation.^30,31 By incorporating the proposed SBF idea into a current selective segmentation model we aim to demonstrate the flexibility of SBF as a fitting term. We recently proposed a selection method in a piecewise-constant framework using a polygon formed by user input, called convex distance selective segmentation (CDSS), ²⁵ which we now introduce briefly. The formulation is given as follows

$E\mathrm{CDSS}(\Gamma ,c1,c2)=\left|\Gamma \right|+\lambda \int \Omega ((z-c1)2\chi 1+(z-c2)2\chi 2)\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\theta \int \Omega Pd(x)\chi 1\text{d}x,$

(23)

where $Pd(x)$ is the normalised Euclidean distance of each point $x\in \Omega$ from its nearest point in the user-defined polygon. Further details are given in Spencer and Chen. ²⁵ Whilst results demonstrate this approach is robust, even in quite difficult cases, it is limited by the piecewise-constant assumption it relies on. We therefore extend this idea to incorporate bias field estimation, which we call selective SBF

$E\mathrm{SSBF}(\Gamma ,c1,c2,\tilde{B})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=\left|\Gamma \right|+\lambda \int \Omega ((z-\tilde{B}c1)2\chi 1+(z-\tilde{B}c2)2\chi 2)\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\mu \int \Omega |\nabla \tilde{B}|2\mathrm{dx}+\gamma \int \Omega (\tilde{B}-\mathbf{s})2\text{d}x\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\theta \int \Omega Pd(x)\chi 1\text{d}x.$

(24)

We minimise this functional equation (24) as outlined in “Iterative minimisation of SBF formulation” and “Numerical implementation” above, except that for Step (3) we use an improved additive operator splitting method from CDSS. ²⁵

Results

For selective SBF we test one image that involves significant intensity inhomogeneity in the foreground and background, shown in Figure 5. The foreground consists of a series of distinct objects that could conceivably be of interest, and was chosen as it is clearly beyond the scope of the piecewise-constant framework used in CDSS. ²⁵ By using just four markers to loosely define the shape of the target object, as well as its location and size, we define a distance selection term P_d that is capable of excluding unwanted artefacts. We demonstrate that we get a successful result for this example, both in terms of the computed contour $\Gamma *$ and the convergence of the intensity constants c₁ and c₂. A particularly challenging aspect of this type of problem can be highlighted by noting that the intensity constants computed are very close: $c_1^{*}=0.35$ and $c_2^{*}=0.33$ . The role of the stabilised bias field, $\tilde{B}$ , is particularly important here.

Figure 5.

Selective SBF Results (a) Successful selective segmentation of the image, z(x), given by contour $\Gamma *$ . (b) Distance selection function, $Pd(x)$ , with user markers. (c) Convergence of c₁ values (200 iterations). (d) Convergence of c₂ values (200 iterations). (e) Computed stabilised bias field, $\tilde{B}(x)$ . (f) Computed indicator function, u(x).

It is worth considering two alternatives for this example, which we now briefly discuss. Firstly, what would the performance of CDSS be like for the image in Figure 5? For brevity, we do not include those results here. As might be expected for a model that relies on a piecewise-constant framework, the results for this image are inadequate as the segmentation favours exterior artefacts to the target object that are of a similar intensity value. Secondly, what does SBF contribute here, i.e. what would Selective VMS (γ = 0 in equation (24)) results be like? Again, we do not include results here, but Selective VMS is capable of achieving a successful segmentation, although as expected c₁ and c₂ do not converge. However, we observe a similar effect as observed in Results Set 2 given by Figure 4. That is, with all other parameters fixed and varying the selection parameter θ, there is a successful result for a wider range of values. We do not know the ground truth for this case, which makes quantifying differences between methods difficult, but we aim to further investigate this phenomenon with different examples.

Conclusions

We have proposed the introduction of a constraint to the VMS model, ²¹ although it applies to any model using bias field correction in this way. It is a framework that provides a link between the Mumford-Shah functional ³ and the piecewise-constant functional of Chan and Vese. ² This constraint does not affect the computation time as we have shown how the exact minimiser can be well approximated in a similar way to Chen et al. ²¹ It is an improvement over current methods in the sense that the intensity constants reliably converge and are feasible in relation to the image. This allows for a meaningful representation of the data by the definition of the image model (equation (3)). We also observe possible advantages with this framework in terms of the quality of the piecewsie-smooth approximation of the image, and a model less reliant on the fitting parameter. We have successfully extended the proposed method to a selective segmentation model, ²⁵ to allow for selection in the presence of intensity inhomogeneity, and have again observed an improvement in terms of parameter dependence. This is a potentially important finding, as this ‘stabilisation’ of the bias field appears to allow for more parameter variation thus improving the reliability of the models. We will investigate this idea further in the future, and attempt to accurately quantify an improvement.