IMAGE SEGMENTATION WITH DYNAMIC ARTIFACTS DETECTION AND BIAS CORRECTION

,

1. Introduction.Image segmentation is the task of partitioning the image domain Ω into homogeneous regions corresponding to individual objects, Ω = i Ω i , or by duality, to find the contours Γ that define the boundaries ∂Ω i of these objects.Image segmentation is most commonly formulated variationally as a minimization problem, where one optimizes a parameterization of the regions or their contours towards specific segmentation criteria encoded in the objective functional.Fundamental image segmentation models underlying the two respective segmentation goals are snakes [19] or geodesic active contours (GAC) [4] for edge-based segmentation, and active contours without edges (ACWE, a.k.a.Chan-Vese or CV model) [5] for region-based segmentation.For both GAC and CV, the segmentation can conveniently be parameterized implicitly using the level set method [37], and various efficient optimization schemes have been presented, e.g.[2,7,15,27].
Here, we are interested in region-based image segmentation, and improvements to the CV model in particular, to segment images that are affected by both image artifacts (single or grouped outliers such as scars, occlusions, scratches) and illumination bias (gain inhomogeneity, shadows).Examples of such images are found in atomic force microscopy (AFM), where these image imperfections are due both to the samples and the image acquisition, in infrared imaging, which suffers from noise and gain inhomogeneity, as well as in medical imaging modalities such as magnetic resonance imaging (MRI; bias field inhomogeneity, coil sensitivity, lesions).In some applications, artifacts (and outliers) are not only a nuisance but the actual object of interest, such as in lesion and tumor detection in brain MRI [38], exploration in geochemistry [16], or domain boundaries of molecular monolayers [9,44].
The CV model is initially derived from the Mumford-Shah functional (MS) [35], which from an input image I 0 : Ω ⊂ R n → R recovers an image I : Ω → R to be smooth almost everywhere: (1) where the three terms aim towards data fidelity, smoothness, and short boundary length; Ω ⊂ R n is the entire image domain and Γ ⊂ Ω the boundary set of Hausdorff dimension n − 1 where the approximation I exhibits discontinuities.Indeed, the CV model is defined as the cartoon-limit (α → ∞) of the MS functional, where the recovered image is required to be piecewise constant, and regions are thus characterized by a single representative color each, µ i : (2) , with λ i and β being parameters.The representative colors µ i are to be determined along with the optimal image domain partitioning Ω = i Ω i .The Chan-Vese active contour model has seen a number of specific modifications.For example, while the initial model could only segment greyscale images into two regions, the authors have extended it to adapt to vector-valued images [8], and they have also illustrated an extension of the model to multi-phase segmentations [42].Unfortunately, all of these models often fail to segment images which feature damage, such as artifacts and bias.Occlusion of objects, scars, and other artifacts on the image, and images featuring regional color inhomogeneity challenge CV and CV-like segmentation models.Several models propose partial solutions.
Regarding images affected by artifacts, Jung, Kang, and Kang [18] created an L 1 variant of the CV scheme, which is more robust to outliers, along with a numerical scheme to minimize it.Indeed, the original quadratic fidelity term in the CV model corresponds to a Gaussian additive noise prior on top of the piecewise constant image model [3].Impulsive noise or other outliers, in particular due to image artifacts, therefore overly impact the functional and skew the image statistics µ i , against which the L 1 -term is more robust.The same idea was used in [28] more recently.
To deal with illumination bias, one could remove regional color inhomogeneity as a preprocessing step, or build a model which explicitly incorporates non-constant piecewise regions.An example of the former is an L 1 -Retinex model [29,47,48], which gathers the piecewise-constant (or relatively constant) structure of an image by removing the smooth bias field.Retinex models remove shadows, but do not guide their corrections to the segmentation.On the other hand, Li et al. [24] propose a region-scalable fitting model which replaces the piecewise constant requirement of the CV model with a local averaging fidelity term which decreases penalty from distant points.This is effectively a partial roll-back from the CV cartoon-limit of the MS-model insofar as regions are not required to be piecewise flat, but piecewise smooth, only.Wang et al. [43] apply this model to medical imaging with success, while Yang et al. [46] use the split Bregman method to expedite segmentation.As illustrated in their papers, the model was able to identify inhomogeneous regions, but, in our experience, the algorithm has trouble when objects are sparse and disconnected.Moreover, the region-scalable fitting model is unable to effectively deal with artifacts.
Another (unspecific) approach is to use some prior knowledge of the object region's shape to guide the segmentation and overcome difficulties with pure intensity cues, e.g.[6,10,11,39].In all cases, these models are able to use shape information to guide the segmentation, but require an accurate shape prior to work with.Moreover, aligning this shape prior can be computationally expensive, or require supervision.
In this paper, we propose a new model that dynamically identifies artifacts and corrects shadows during the segmentation.Our goal is to formulate a single variational model to deal with both artifact detection and illumination bias correction.We introduce a binary artifact label X, that marks individual pixels as artifacts if they violate the two-phase piece-wise constancy assumption of the CV segmentation model.This approach effectively corresponds to statistical hypothesis testing; a pixel is classified as artifact if it fails the Gaussian null-hypothesis by virtue of a z−score more extreme than a given threshold.This threshold parameter thus directly controls the statistical significance level associated with the artifact classification, or the expected false-positive rate, respectively.We call this the CV+X model.On the other hand, to deal with illumination bias, we include a Retinexlike image decomposition.Starting from a simple additive image formation model, I = B + S, we decompose the input image I into a structure part S and a smooth bias B. The structure part is expected to be two-phase piecewise constant, as modeled by the CV-terms.We call this second variant the CV+B model.Combining both parts performs image segmentation coupled with dynamic artifacts detection and bias estimation as intended.The resulting full CV+XB model yields as output: A binary segmentation u, the region characteristics µ 1 and µ 2 , an artifact map X, a nearly two-phase structure estimate S, and a smooth bias field B. We also propose a numerical scheme to minimize this model efficiently, based on MBO-like threshold dynamics [13,31,32], which is much faster than the original levelset approach [5].
The remainder of this paper is structured as follows.First, in §2 we provide some more complementary background on the CV model as relevant to our work.We then address the artifact detection ( §3) and bias correction ( §4) problem separately, before we put them together in a single combined variational model ( §5).We present results and a discussion of our model in §6 and conclude in §7.

2.
More background on the CV model.In the following paragraphs, we review some aspects of and modifications to the CV model [5], both in terms of model formulation and algorithms for its optimization, that are going to guide the development of our proposed artifact-and bias-resistant segmentation model.

2.1.
Levelset formulation and gradient descent.The reduction of the segmentation model (2) to two phases and using the level-set representation of these two phases leads to the following classical formulation [5]: (3) where the levelset function φ is positive in object regions, negative in background regions, zero on the object boundaries, and H is the Heaviside function.The last term represents the total variation of the characteristic function of the object, H(φ), and is the co-area-formula equivalent of the boundary length term.The original and immediate strategy to solve (3) is gradient descent: , which is to say that µ 1 and µ 2 are the average colors of their respective regions, and (5) In practice, the distributional Heaviside H and its derivative δ need to be smoothly approximated by analytical H ε and δ ε = H ε (φ), such as the logistic function and its derivative.Irrespective of the chosen Heaviside approximation H ε , the CV-PDE ( 5) is hard to solve, due to the presence of the divergence term associated with the interface perimeter.Further, many implementations of the Chan-Vese model use a modified descent without the Dirac-masking (corresponding to H ε (φ) = φ and δ ε (φ) = 1).Such a gradient descent corresponds to an energy functional homogeneous of degree 1 in φ; for minimizers to exist the level-set function needs to be restricted, e.g., to a signed distance function [15,25,26] or to φ : Ω → [0, 1] (phase-field approximation) [7].
2.2.Threshold dynamics for mean curvature motion.An important contribution towards more efficient CV optimization stems from the phase-field approach outlined in [13], and goes back to the diffusion-threshold scheme for approximating motion by mean curvature proposed by Merriman, Bence, and Osher [32].The fundamental idea is to reproduce the motion by mean curvature due to the divergence term in the CV-PDE (5) by more efficient means.For this, let us replace the total variation of the phase field u ∈ W 1,2 (Ω), by the real Ginzburg-Landau functional [33] (we use the variable u instead of φ to highlight Inverse Problems and Imaging Volume 11, No. 3 (2017), 577-600 its phase-field nature rather than just being a generic level set function): where W (s) is a double well potential with two equal minima at s = 0 and s = 1, for example W (s) := s 2 (1 − s) 2 .Minimizing this functional yields a phase field that is smooth and tends to be binary.The GL-functional Γ-converges to the total variation functional of binary classifiers (u : Ω → {0, 1} and ε → 0) [33]: where σ(W ) is a surface tension term depending on the double well potential.The minimizing flow of this functional for ε → 0 + produces motion by mean curvature of the interface, which is exactly what one needs in the CV model minimization.
The gradient descent of ( 6) is the Allen-Cahn equation ( 8) The MBO scheme can be motivated by a two-step splitting scheme for the Allen-Cahn equation: where the first step is simply diffusion by the heat equation, and the second step is an ordinary differential equation minimizing the double well potential.The implicit heat step is carried out using the fast Fourier transform, which diagonalizes the Laplace operator in a spectral method [40].The two-step method above is then modified again to produce the MBO-scheme [32], involving alternating between a heat equation step and thresholding (projection onto the binary set {0, 1}): 2.3.CV model with threshold dynamics.Starting from the phase-field approximation H(φ) ≈ u, where u : Ω → [0, 1], Esedoglu and Tsai devise the following modified GL-based diffuse interface approximation to the CV-model: They propose several different time splitting variants to incorporate the image segmentation terms into a threshold dynamics scheme for mean curvature motion [13].Their main scheme is to combine heat equation and data-fidelity forcing term in one step, and perform thresholding in the second step.Here, we are more inspired by the variant that iterates over three steps, namely propagation of the image fidelity term ODE, followed by heat diffusion, and eventually thresholding.For a more complete discussion of threshold dynamics and energy functions that they minimize, see [12].
3. Image segmentation with dynamic artifact detection.The main issue to be addressed in this section is the presence of outliers (individual or small groups of pixels) in images, that do not comply with the piecewise constant image model underlying the CV image segmentation.These outliers greatly affect the regions' intensity statistics µ 1 and µ 2 and may effectively lead to image segmentation failure.In some cases, artifacts may even be falsely detected as objects of interest in the image, because their presence outweighs the actual image contrast between background and object.This issue was partially addressed by [18,28] through the substitution of the quadratic data fidelity by the absolute fidelity term.However, this also represents a major difference in the underlying noise prior (Gaussian noise versus impulsive noise), and thus has consequences beyond dealing with artifacts.Also, in these models, artifacts are not actively dealt with and they remain part of the segmentation.
If artifacts were readily identifiable through preprocessing steps such as thresholding or other heuristics, then one simple approach is to eliminate the fidelity term in these regions.Here, however, we integrate the artifact detection into the very same variational model.
While in some specific applications artifacts may follow known statistical distributions in terms of shape, location, size or appearance, in the general case this is not true.Instead, here we describe artifacts by exclusion, only knowing what they are not.Indeed, we may characterize artifacts as isolated or small groups of pixels that fail to adhere to the two-phase piecewise constant image model, without any other prior regarding their intensity distribution or artifact shape regularity.Definition 3.1.We denote X : Ω → {0, 1} as an artifact indicator function.We accommodate this additional information in the CV model by using X as a mask on the data-fidelity term.To control the artifact detection rate, we add a penalty γ on the size of the artifact class, and optimize for X along with the actual segmentation.The CV image segmentation with dynamic artifact labeling model is defined as follows: The preliminary interpretation is that now the artifact indicator function X optimizes a trade-off between size penalty γ and outlier elimination.This will become much clearer once we look at the actual minimization scheme, further below.

Remark 1.
In support of this model, it is important to note that similar approaches are employed in other imaging tasks: in image sequences, Ayvaci, Raptis, and Soatto [1] and Estellers and Soatto [14] identify pixels overly violating the optical flow constraint as occlusions, while Yan [45] localizes individual damaged pixels for impulse-noise image in-painting.4. Image segmentation with bias correction.The consequence of inhomogeneous illumination is that the same material (object or background) does not have a uniform appearance over the entire image domain.More specifically, the image pixels of a region of interest are not all drawn from the same idealized Gaussian distribution, since the mean intensity of such a region varies across Ω due to inhomogeneous illumination, and a region cannot be characterized by its mean intensity and noise variance alone, anymore.More abstractly, the very two-phase piecewise constancy assumption of the Chan-Vese [5] cartoon limit to the Mumford-Shah [35] model is violated.One way to deal with this complication is to abandon the CV assumptions partially, and to model images to be segmented as piecewise smooth, instead.This path was chosen by Li, Kao, Gore and Ding, with their popular region-scalable fitting energy for image segmentation [24].In contrast, in [30] the authors proposed to multiply a spatially varying illumination field based on spline interpolation into the CV model.Here, we want to maintain the piecewise constancy assumption of the CV model, and deal with the illumination inhomogeneity by decomposing the image into its (piecewise constant) structure and an (illumination) bias part, inspired by Retinex theory.
Retinex is a theory on the human visual perception [21,22,23].It was an attempt at explaining how a combination of processes supposedly taking place both in the retina and the cortex is capable of adaptively coping with illumination that varies spatially [47,48].The fundamental observation is the insensitivity of human visual perception with respect to a slowly varying illumination on a Mondrian-like piecewise constant scene.
Starting from an additive image formation model2 , ( where the illumination b is supposed to vary smoothly, the spatial derivatives of the observed intensity i are mostly due to edges in the structure s.This is the core assumption to most Retinex implementations, such as [17,20,34].Here, we propose to integrate this Retinex model into the CV image segmentation model, by segmenting only the piecewise constant image part while extracting the smoothly varying illumination bias.The most immediate solution might consist in performing Retinex image decomposition as a preprocessing step to CV image segmentation.This could be an unwise strategy, however, since the CV image model is not just expecting a piecewise constant input image, but more specifically a two-phase piecewise constant image, which is a much stronger image prior, simply ignored in such Retinex preprocessing.
Instead, we propose to integrate the Retinex image decomposition and the CV image segmentation problem into a single variational model: We are interested in an image decomposition consisting of a smooth bias field and a two-phase piecewise constant structure part with regular phase interface.In this, rather than segment an image directly, the model seeks to segment the image's underlying structure.We optimize segmentation and bias correction simultaneously.denote the estimated bias and structure part, respectively.Then, we define the energy functional for joint image segmentation and bias correction as follows: where the first three terms are the classical CV model acting on the structure S instead of the input image I 0 , while the new term models bias smoothness.
Remark 2. Note that the proposed model ( 16) differs from a simple two-phase approximation of the MS-functional (1), or the region-scalable fitting model [24], in that the bias field, B, representing the smooth variations, is restricted to be smooth everywhere, including at the phase interface.

5.
Combined artifact detection and bias corrective image segmentation model.Having tuned the original CV model to deal with artifact or illumination bias corrupted images, respectively, we now have all ingredients ready to build a single variational, region-based image segmentation model that inherently detects local artifacts and corrects smoothly varying illumination bias.Having both components is crucial, because either type of corruption has the critical potential to prevent correctly handling the other.On the one hand, artifacts and outliers can exert a significant negative impact on most Retinex decomposition models, since their resulting high contrast cannot be accommodated in the bias field but would need to be attributed to the structure, distorting the intensity information even outside the actual artifact regions.Uneven illumination, on the other hand, has the potential to identify artifacts falsely if it pushes the intensities of some pixels outside their respective region's "comfort zone" in terms of distance from the mean intensity.Having artifact and bias terms in a single variational functional allows handling of these complications concurrently.
Remark 3. It is important to note that in this combined model, the artifact variable X only masks the CV terms on the structure part S, while interface perimeter T V (u) and bias field smoothness |∇B| 2 are not affected.and then solving the resulting unconstrained minimization problem (the Penalty method ).It is well-known that strictly enforcing the image formation model using this penalty, only, would require ρ → ∞, which renders the functional illconditioned.Therefore, if the constraint is to be imposed strictly, in addition to the quadratic penalty term, a Lagrangian multiplier can be employed while keeping ρ finite [36] (the Augmented Lagrangian method).

Optimization strategy.
Due to the newly introduced variables modulating the data-fidelity term, the resulting energy functional in (µ 1 , µ 2 , u, X, B, S) is nonconvex and the global minimizer results from [7] do not fully apply to the associated minimization problem.Due to the non-convex and non-linear nature, we propose solving the minimization problem in a coordinate descent and dual ascent approach.

5.3.1.
Minimization w.r.t.µ.The first problem in ( 21) is to update the region statistics.
Lemma 5.3.Let both inside and outside, non-artifact regions be non-empty: Then the updates in (21a) are found as , where we have omitted most occurrences of • n and • n+1 for simpler notation.
Proof.The relevant optimization problems in (21a) are Both objective functions are proper, convex, and differentiable in µ 1,2 .Their respective unique minimizer is found by letting the first variation vanish.We obtain 0 from which the lemma immediately results.Remark 4. The computation of the mean intensities takes the usual form as known from the CV-model [5], but is computed based on the structure part S and restricted to regions not being artifacts.Note that situations where these regions happen to be empty are degenerate and need to be handled separately.
As outlined in §2.2 and §2.3, one way of overcoming the stiffness of the gradient flow associated with the boundary term |∇u| is to use time splitting.The idea behind our proposed minimization scheme is to combine the threshold dynamics of Esedoglu and Tsai's [13] CV-model with the artifact detection and bias-structure decomposition.To this end, we perform time-splitting on the gradient descent

Inverse Problems and Imaging
Volume 11, No. 3 (2017), 577-600 scheme w.r.t. the phase-field u, to separate the data-fidelity ODE from the heatdiffusion PDE, followed by thresholding for rectification.
We look for a minimizer to (21b)/( 26) by evolving the phase-field u for some time τ u each along the following directions: 1.According to the CV-data term on the structure part, with v(x, 0) = u n (x) as initial condition, 2. According to the heat equation with the previously propagated phase-field w(x, 0) = v(x, τ u ) as initial condition and using appropriate boundary conditions.3. Rectify the diffused w by thresholding: Note that β replaces β as the parameter controlling the weight of the interface regularity prior; This β and the original β are non-trivially related through the propagation time τ u [13].Here, we propose to tune the new β directly.
Remark 5.In the absence of artifacts, the phase-field evolution is similar to the approach in [13].Where artifacts have been detected, however, the phase-field ignores the observed and allegedly unreliable structure information and is purely driven by the interface regularizing MBO terms.
5.3.3.Minimization w.r.t.X.We now consider the update of the artifact indicator function X, according to (21c).
Lemma 5.4.The update X n+1 is locally found as Proof.The update must solve the minimization problem where the local cost associated with X is Locally, the optimal X n+1 thus satisfies (30) For each pixel x ∈ Ω, only two cases need to be considered: ( 31) Inverse Problems and Imaging Volume 11, No. 3 (2017), 577-600 Thus X(x) = 1 is optimal if c(x) ≤ 0 from which the lemma immediately results.
Remark 6.We include this proof despite being trivial, because it highlights the impact of the artifact indicator X.Indeed, let us assume that the parameters λ 1 and λ 2 represent the inverse of the Gaussian noise variance from a Bayesian interpretation of the CV model [3].Without loss of generality, consider now a pixel x in the object region, u = 1.The relevant term in ( 27) is thus which is the squared z-score (standard score) of the local intensity under a Gaussian distribution of corresponding region statistics.This squared z-score is then compared against the threshold γ.The immediate interpretation is thus the following: the artifact classification is effectively a concealed statistical hypothesis z-test of the pixel intensity with a Gaussian distribution N (µ i , λ i ) as null-hypothesis, and a pixel is classified as an artifact if the z-score of its intensity is more extreme than √ γ.The model parameter γ is thus intimately related to the level of statistical significance attached to the artifact classification and its expected false positives rate.Proof.These EL-equations directly result from simple calculus of variations on the three relevant terms of (21d) involving B: Optimality requires the first variation to vanish: from which ( 32) is immediately obtained by rearrangement.
The update rule for B n+1 in ( 32) is simply the heat equation, which can be efficiently solved, e.g., spectrally using the Fourier transform F{•}(ω) (thus implicitly assuming periodic boundary conditions): Remark 7. The immediate interpretation of the B-update is as follows: B strives to optimize the trade-off between image formation constraint I 0 −S and the smoothness prior imposed by |∇B| 2 , and B is updated as a lowpass filtered version of what I 0 −S suggests.As will be seen farther below, the Lagrangian multiplier λ accumulates the error on the reconstruction constraint, and will eventually enforce that constraint strictly.
5.3.5.Minimization w.r.t. S. Similarly, let us now determine the EL-equations governing the update of the structure S.
Lemma 5.6.The update rule for S reads Proof.We consider the relevant minimization problem (21e), (37) S n+1 = arg min The EL equation of the quadratic, differentiable problem ( 36) is immediately obtained by having the first variation in S vanish (38) from which the lemma immediately follows by rearrangement.
Remark 8.This update for the structure variable S again admits an immediate and intuitive interpretation: S balances the trade-off between the image decomposition rule I −B and compliance with the piecewise constancy terms due to the CV model, achieved by averaging the two contributions.In the presence of artifacts, however, the structure is exclusively derived from the image formation model.Again, the error-accumulating Lagrangian multiplier adds the error back into the updates.
5.3.6.Maximization w.r.t.λ.Finally, we note that the Lagrangian multiplier is simply updated through gradient ascent (dual ascent), through propagation of the following ODE (39) from which the error-accumulating nature of the Lagrangian multiplier becomes evident.

5.4.
Algorithm.Putting all previous results together as described, we propose the following algorithm for the CV+XB model: Algorithm 1.The joint image segmentation, artifact detection, and bias correction problem can be solved by minimizing the constrained CV+XB energy (17), viz.by solving the associated Augmented Lagrangian saddle point problem (20).Starting from some initialization u 0 , X 0 = 0, S 0 + B 0 = I 0 , and λ 0 = 0 we iterate these steps until convergence: 1. Update µ 1 and µ 2 according to lemma 5.3 2. Evolve phase field through time-split threshold dynamics as described in §5.3.2 3. Update the artifact classification X pointwise, according to lemma 5.4 4. Update the bias field estimate B spectrally, according to lemma 5.5 5. Update the structure estimate S pointwise, according to lemma 5.6 6. Update the Lagrangian multiplier λ by dual ascent (21f) where we define convergence as the situation where subsequent phase fields produce the same segmentation.
Remark 9. Note that due to their nature as exact coordinate minimizations, the update steps 1 and 3-5 are non-increasing in the Augmented Lagrangian (19), while the dual ascent step 6 is non-decreasing in the AL (strictly increasing outside the feasible set).The behavior of the phase-field evolution in step 2 is more convoluted, and thus overall convergence of the scheme to a stationary point is not immediately obvious.The convergence behavior of the isolated threshold dynamics scheme is discussed in [12].In our numerical experiments, three different types of behavior have been observed in numerical experiments: 1) most often, convergence to a stationary point (local optimum); 2) rarely, convergence to a degenerate case (trivial solutions); and 3) very rarely, convergence to a limit cycle.
Practical results obtained with this algorithm are presented in §6, below.5.5.Parameter influence and algorithm variants.Here, we briefly summarize the various model parameters and their qualitative impact on the expected optimization outcome, before we consider a few limit cases and sub-models.5.5.1.Model parameters.The model in (17) contains five main parameters with the following roles: λ i : As noted in remark 6, these can be interpreted as the inverse of the noise variance of the region i.Larger values of λ i lead to increased data-fidelity.β: This parameter controls the regularity of the segmentation interface.In the algorithm we replace it by β , which directly tunes the amount of diffusion in the threshold-dynamics evolution of the phase field.Larger values of β (β ) lead to shorter domain interfaces.γ: The parameter γ controls the false positives rate of the artifact detection (see also remark 6).Larger values of γ result in more conservative artifact classifications.α: This parameter controls the smoothness of the bias field estimate.Larger values of α lead to smoother bias field estimations.In addition to these model parameters, the Augmented Lagrangian (19) and the algorithm introduce three more computational parameters: ρ: This defines the relative importance of the quadratic image formation penalty; the larger ρ, the less slack between S and B during their evolution.τ λ : Time step of the Lagrange multiplier update.Typical values 0 ≤ τ λ ≤ ρ.For τ λ = 0, the Lagrange multiplier becomes inactive and the image formation constraint is not strictly enforced.τ u : This is the time step parameter of the phase-field evolution through threshold dynamics.The parameter balances precision (τ u small) with speed (τ u big); Too small parameter choice leads to early freezing or pinning of the phase-field [31,32,41].

Limit cases.
Here we briefly describe how particular limit-cases for parameter choice lead to reduced models: α → ∞: Letting the bias smoothness parameter α go to infinity results in a biasstructure decomposition where the bias field estimate is flat (constant).The structure part S therefore captures the entire contrast of the input image.As a result, the model then corresponds to (14) and performs Chan-Vese image segmentation with artifact detection.In the experiments of §6, we call this model CV+X.The corresponding algorithm is most simply obtained from algorithm 1 by replacing the B and S updates by B = 0 and S = I 0 , respectively.Indeed, B = 0 results from Inverse Problems and Imaging Volume 11, No. 3 (2017), 577-600 (35) with α → ∞ without loss of generality, and consequently the image formation constraint results in S = I 0 .γ → ∞: Choosing the artifact penalty parameter γ infinitely large essentially leads to X = 0 everywhere, which removes the artifact detection from the model.This reduced model then corresponds to (16) and performs Chan-Vese segmentation on the structure part of the bias-structure decomposed input image.We call this model CV+B.We obtain the corresponding algorithm from algorithm 1 immediately by simply replacing the artifact detection update by X = 0, which is the natural consequence of ( 27) with threshold γ → ∞. α, γ → ∞: Letting both parameters go to infinity removes both the artifact detection and the bias-field estimation from the segmentation model, and we find ourselves with the phase-field parameterization of the original Chan-Vese (CV) image segmentation model (3).
6. Results and discussion.We have implemented the complete CV+XB model according to the algorithm outlined in algorithm 1 of §5.2, in MATLAB3 .This full model includes the sub-models CV, CV+X, and CV+B, easily selected by simply suppressing some or all of the sub-updates to X, S, and B listed in algorithm 1.First, we use test images provided4 by the authors of the region-scalable fitting model [24].These five images mainly illustrate the use of bias field estimation incorporated into the segmentation problem.In addition, for a specific case for the artifact detection model, we test our model on a synthetic test image of ellipses 5 .
All six test images have been processed with all four models, and the results are illustrated in figures 1-3.Images and results are described in detail in the figure captions.Slight tuning of model parameters was necessary between images (but not between models for a single image).This is natural since images greatly vary in noise level, bias smoothness, outlier presence, etc.The success of the segmentation models depends on the chosen model and image considered; naturally, not all of the images are equally affected by bias and artifacts, and only appropriate priors and model choices fit.Most images worked impressively well, and the algorithm converged rapidly, typically within 10 and 50 iterations, amounting to fractions of a second on standard personal computing equipment (2011 laptop with Intel Core i7 2.80 GHz (dual core), 4GB RAM, running 64-bit Fedora 20 Linux OS, and MATLAB 2011b), as reported in table 1.
For comparison with a state-of-the-art bias-resistant image segmentation algorithm, we provide results obtained using the code provided in [24] for the first 5 test images, as well, shown in figure 4. As can be seen from this figure, the method presented in [24] achieves reasonable bias-estimation, but requires careful initialization, as the authors noted.Also, the gradient descent code provided for this method takes on the order of 1-2 seconds per image, which is almost an order of magnitude slower than our proposed method.Note that the method in [24] can be somewhat accelerated by making use of splitting techniques [46]; however, the algorithm run times reported therein suggest that our proposed model still converges faster.
We further apply the proposed models to real images from atomic force microscopy (AFM) and fluorescence microscopy.In figures 5 and 6 we present both the classical CV result as well as the most appropriate proposed model for each case Table 1.Algorithm convergence: number of iterations and required computer time.The algorithm is deterministic.Convergence is defined as the phase-field u not changing during its update.Computation time is statistical depending on CPU scheduling; here, we report average numbers over 50 repetitions.The extra cost of artifact is negligible compared to basic CV, in particular since it may speed up convergence in appropriate images.Bias correction roughly doubles the computational load per iteration, which is an acceptable price for its benefits when appropriate.).In all these images, our models allow significant improvement of the segmentation, while also providing additional information about artifact location and/or structure-bias decomposition, where applicable.

7.
Conclusions.In this paper, we have presented an energy functional that describes joint region-based image segmentation, dynamic artifact classification, and bias field estimation.Combining these tasks within a single functional allows simultaneous solution of the different problem aspects in a common variational approach, as opposed to more error-prone sequential processing.Indeed, joint optimization enables making full use of the few priors employed in the model, while sequential processing would each time focus on subsets of priors, only.We show relations between the proposed model and state-of-the-art techniques in image processing: the proposed artifact handling is related to recent results in occlusion detection in optical flow modeling [1,14] and outlier handling for impulsenoise denoising [45].In our case, we can make an elegant connection between the artifact classification performed by our model and statistical hypothesis testing: the model classifies pixels as artifacts if their z−score under a Gaussian noise model is more extreme than the parameter γ, which means that we have direct control over the expected error rate on artifact labeling.The bias field estimation, on the other hand, is intimately tied to the Retinex model for image decomposition.However, using a joint model, our decomposition makes direct use of the two-phase piecewise constancy assumption of the CV-model, which is a much stronger prior on the structure part than typical Retinex models employ.
In addition to formulating a joint model for these tasks, we also devise an efficient algorithm for the optimization of the proposed variational problem.Our algorithm makes use of the phase-field parameterization and threshold-dynamics [32] that have successfully been introduced to the CV-model in [7].The algorithm solves the saddle-point problem associated with the constrained functional optimization by iterating several steps of coordinate descent and dual ascent.The resulting The goal is to separate the two black ellipses from the gray background, considering the white ring to be an occlusion artifact.CV, however, groups the white ring with the light background.CV+X successfully identifies the ring as artifact, and closes the black ellipses thanks to the interface regularization.The bias correction is misleading, since much of the white ring will be considered overly illuminated background (CV+B, CV+XB), the corners being captured as transient artifacts (CV+XB).
subproblems have closed form solutions (region statistics µ 1 , µ 2 , artifacts map X, and structure and bias field decomposition S and B, respectively), or admit efficient gradient descent/ascent steps (segmentation phase-field u, through modified MBO, Input/Init RSF [24] Input/Init RSF [24] i  Results by region-scalable-fitting (RSF) [24].Left: Starting from generic initialization (as used in our method), RSF fails to capture the correct bias/segmentation result.Right: Starting from a tuned initialization (as provided in [24]), RSF produces the desired bias resistant region-based image segmentation result.Note that for optimally chosen number of iterations, i (as provided in [24]), the computation time is about an order of magnitude slower than our proposed method, as reported in table 1.
and dual update λ).The proposed algorithm typically converges in about 10 to 50 iterations and takes fractions of a second on standard computing equipment.Practical results from disparate imaging modalities illustrate the range of images that are successfully modeled by our joint energy functional.Finally, the MATLAB code implementing our algorithm and that produces the results presented in this paper is made available at http://www.math.montana.edu/dzosso/codeand at MATLAB Central.

5. 1 .
Combined model.Formulating this single variational model involving both dynamic artifact detection and bias correction is now merely a matter of combining the respective definitions of CV+X(14) and CV+B (16): Definition 5.1.Let X : Ω → {0, 1} denote the artifact indicator function, and let B, S : Ω → R be the bias and structure decomposition of the input image I 0 : Ω → R. We define the CV+XB functional for joint image segmentation, artifact classification, and bias correction as the following energy in terms of the phase-field u : Ω → [0, 1]:(17)

Figure 1 .Figure 2 .Figure 3 .
Figure 1.Example cases 1 & 2. Top: Coronal MRI slice.Input image and initial segmentation contour are shown in the top-right corner.The image is heavily affected by intensity bias, such that the classical CV model fails.The superior parts of the white-matter are undersegmented, while the inferior regions are markedly oversegmented.CV+X is not very helpful, here.In contrast, CV+B fixes the problem: the extracted structure is nearly flat (the brain is essentially two-phase piecewise constant), while CV+XB marks some non-brain pixels as outliers.Bottom: This synthethic image is a combination of piecewise constant regions affected by strong noise and oscillating bias.Again, CV fails, and artifacts detection not appropriate.Bias correction greatly improves the segmentation, but errors persist since from the simple initialization the algorithm converges to a wrong local minimum (upper part of left hand structure).Inverse Problems and ImagingVolume 11, No. 3 (2017), 577-600

Figure 4 .
Figure 4.Results by region-scalable-fitting (RSF)[24].Left: Starting from generic initialization (as used in our method), RSF fails to capture the correct bias/segmentation result.Right: Starting from a tuned initialization (as provided in[24]), RSF produces the desired bias resistant region-based image segmentation result.Note that for optimally chosen number of iterations, i (as provided in[24]), the computation time is about an order of magnitude slower than our proposed method, as reported in table 1.

Figure 5 .
Figure 5. Microscopy example cases.Top: This AFM image is severely affected by inhomogeneity compared to pattern contrast.As a result, classical CV segmentation fails.The proposed CV+XB model is able to capture some of the bias as such, and correctly segments some of the actual pattern.The central-square contour initialization, however, provokes an incorrect bias field estimate at early stages of the optimization, and leads to an incorrect local minimum, misclassifying the central portions.Middle: Starting from a near-optimal contour initialization of a single pattern element, the CV model still fails entirely.The proposed CV+XB model is not misled into incorrect minima, anymore, and successfully separates bias from actual pattern contrast.Bottom: The CV+XB model captures most of the inhomogeneity present in this AFM sample image and leads to reasonable segmentation of the diamond pattern.(Note: The inversion of foreground/background between CV and CV+XB models is arbitrary and triggered by domain size.)