Multi-Label Segmentation via Residual-Driven Adaptive Regularization
Byung-Woo Hong Ja-Keoung Koo Stefano Soatto








Multi-Label Segmentation via Residual-Driven Adaptive Regularization
    Byung-Woo Hong    Ja-Keoung Koo

Chung-Ang University, Korea

{hong,jakeoung}@cau.ac.kr









Stefano Soatto

University of California Los Angeles, U.S.A.

soatto@cs.ucla.edu

    
======================================================================================================================================================================================





We present a variational multi-label segmentation algorithm based on a robust Huber loss for both the data and the regularizer, minimized within a convex optimization framework. We introduce a novel constraint on the common areas, to bias the solution towards mutually exclusive regions.
We also propose a regularization scheme that is adapted to the spatial statistics of the residual at each iteration, resulting in a varying degree of regularization being applied as the algorithm proceeds: the effect of the regularizer is strongest at initialization, and wanes as the solution increasingly fits the data. This minimizes the bias induced by the regularizer at convergence. We design an efficient convex optimization algorithm based on the alternating direction method of multipliers using the equivalent relation between the Huber function and the proximal operator of the one-norm. We empirically validate our proposed algorithm on synthetic and real images and offer an information-theoretic derivation of the cost-function that highlights the modeling choices made. 





§ INTRODUCTION


To paraphrase the statistician Box, there is no such thing as a wrong segmentation. Yet, partitioning the image domain into multiple regions that exhibit some kind of homogeneity is useful in a number of subsequent stages of visual processing. So much so that segmentation remains an active area of research, with its own benchmark datasets that measure how right a segmentation is, often in terms of congruence with human annotators (who themselves are often incongruent).

The method of choice is to select a model, or cost function, that tautologically defines what a right segmentation is, and then find it via optimization. Thus, most segmentation methods are optimal, just with respect to different criteria.

Classically, one selects a model by picking a function(al) that measures data fidelity, which can be interpreted probabilistically as a log-likelihood, and one that measures regularity, which can be interpreted as a prior, with a parameter that trades off the two. In addition, since the number of regions is not only unknown, but also undefined (there could be any number of regions between one and the number of pixels in a single image), typically there is a complexity cost that drives the selection of a model among many.

Specifically for the case of multi-label, multiply-connected region-based segmentation, there is a long and illustrious history of contributions too long to list here, but traceable to <cit.> and references therein.



77pt





In each and every one of these works, to the best of our knowledge, the trade-off between data fidelity and regularization is determined from the distribution of the optimization residual and assumed constant over all the partitioning regions, leading to a trade-off or weighting parameter that is constant both in space (, on the entire image domain <cit.>) and in time, , during the entire course of the (typically iterative) optimization.

Neither is desirable: consider Fig. <ref>. Panels (c) and (d) show the residual and the variance, respectively, for each  region shown in (b), into which the image (a) is partitioned. Clearly, neither the residual, nor the variance (shown as a gray-level: bright is large, dark is small), are constant in space. Thus, we need a spatially adapted regularization, beyond static image features as studied in  <cit.>, or local intensity variations <cit.>. While regularization in these works is space-varying, the variation is tied to the image statistics, and therefore constant throughout the iteration. Instead, we propose a  spatially-adaptive regularization scheme that is a function of the residual, which changes during the iteration, yielding an automatically annealed schedule whereby the changes in the residual during the iterative optimization gradually guide the strength of the prior, adjusting it both in space, and in time/iteration.

To the best of our knowledge, we are the first to present an efficient scheme that uses the Huber loss for both data fidelity and regularization, in a manner that includes standard models as a special case, within a convex optimization framework. While the Huber loss <cit.> has been used before for regularizations <cit.>, to the best of our knowledge we are the first to use it both in the data and regularization terms.

Furthermore, to address the phenomenon of proliferation of multiple overlapping regions that plagues most multi-label segmentation schemes, we introduce a constraint that penalizes the common area of pairwise combinations of partitions. The  constraints often used to this end are ineffective in a convex relaxation <cit.>, which often leads to the need for user interaction <cit.>.

To boot, we present an efficient convex optimization algorithm within in the alternating direction method of multipliers (ADMM) framework <cit.> with a variable splitting technique that enables us to simplify the constraint <cit.>.




 §.§ Related Work


One of the most popular segmentation models relies on the bi-partition/piecewise-constant assumption <cit.>, which has been re-cast as the optimization of a convex functional based on the thresholding theorem <cit.>. 
In a discrete graph representation, global optimization techniques are developed based on min-cut/max-flow <cit.>, and there is an approach that has been applied to general multi-label problems <cit.>. Continuous convex  relaxation techniques have been applied to multi-label problems in a variational framework <cit.>, where the minimization of Total Variation (TV) is performed using a primal-dual algorithm. In minimizing TV, a functional lifting technique has been applied to the multi-label problem <cit.>.
Most convex relaxation approaches for the multi-label problems have been based on TV regularization while different data fidelity terms have been used by L_1 norm <cit.> or L_2 norm <cit.>. 
For the regularization term, the Huber norm has been used for TV in order to avoid undesirable staircase effects <cit.>. There have been adaptive algorithms proposed to improve the accuracy of the boundary using an edge indicator function <cit.>, generalized TV <cit.>, or an anisotropic structure tensor <cit.>. 
A local variation of the image intensity within a fixed size window has been also applied to consider local statistics into the regularization <cit.>, and the regularization parameter has been chosen based on the noise variance <cit.>.
Most adaptive regularization algorithms have considered spatial statistics that are constant during the iteration, irrespective of the residual. 

It has been known that most multi-label models suffer from inaccurate or duplicate partitioned regions when used with a large number of labels <cit.>, which forces user interactions including bounding boxes <cit.>, contours <cit.>, scribbles <cit.>, or points <cit.>. 
Alternatively, side information such as depth has been applied to overcome the difficulties that stem from uncertainty in characterization of regions in the multi-label problem <cit.>.

The words “deep learning” appear nowhere in this paper but this sentence: we believe there are deep connections between the dynamic data-driven regularization we have proposed and a process to design models that best exploit the informative content of the data, learning which can inform more useful models moving forward. Sec. <ref> is a first step in this direction.




 §.§ Contributions


Our first contribution is a multi-label segmentation model that adapts locally, in space and time/iteration, to the (data-driven) statistics of the residual (Sec. <ref>). 

The second contribution is to introduce a Huber functional as a robust penalty for both the data fidelity and the regularization terms, which are turned into proximal operators of the L_1 norm, allowing us to conduct the optimization efficiently via Moreau-Yosida regularization (Sec. <ref>). 

Third, unlike most previous algorithms that were ineffective at dealing with a large number of labels, we introduce a constraint on mutual exclusivity of the labels, which penalizes the common area of pairwise combination of labels so that their duplicates are avoided (Sec. <ref>).

Finally, we give an information-theoretic interpretation of our cost function, which highlights the underlying assumption and justifies the adaptive regularization scheme in a way that the classical Bayesian interpretation is not suited for the derivation of our model (Sec. <ref>).

We validate our model empirically in Sec. <ref>, albeit with the proviso that existing benchmarks are just one representative choice for measuring how correct a segmentation is, as hinted at in the introduction, whereas our hope is that our method will be useful in a variety of settings beyond the benchmarks themselves, for which purpose we intend to make our code available open-source upon completion of the anonymous review process.




§ MULTI-LABEL VIA ADAPTIVE REGULARIZATION





 §.§ General Multi-Label Segmentation
 

Let f : Ω→ be a real valued[Vector-valued images can also be handled, but we consider scalar for ease of exposition.] image with Ω⊂ℝ^2 its spatial domain.
Segmentation aims to divide the domain Ω into a set of n pairwise disjoint regions Ω_i where Ω = ∪_i=1^n Ω_i and Ω_i ∩Ω_j = ∅ if i ≠ j.  
The partitioning is represented by a labeling function l : Ω→Λ where Λ denotes a set of labels with | Λ | = n. The labeling function l(x) assigns a label to each point x ∈Ω such that Ω_i = { x   |   l(x) = i }.  
Each region Ω_i is indicated by the characteristic function χ_i : Ω→{ 0, 1 } defined by:


    χ_i(x) = 
       
       1     : l(x) = i,
    
       0     : l(x) ≠ i.


Segmentation of the image f(x) can be cast as an energy minimization problem in a variational framework, by seeking for regions {Ω_i } that minimize an energy functional with respect to a set of characteristic functions {χ_i }:


    ∑_i ∈Λ{λ 𝒟 (χ_i) + (1 - λ)  ℛ (D χ_i) },  ∑_i ∈Λχ_i(x) = 1,


where 𝒟 measures the data fidelity and ℛ measures the regularity of D χ_i, and D denotes an operator for the distributional derivative.
The trade-off between the data fidelity 𝒟 and the regularization ℛ is controlled by the relative weighting parameter λ∈ [0, 1]. 
A simple data fidelity can be designed to measure the homogeneity of the image intensity based on a piecewise constant model with an additional noise process: f(x) = c_i + η_i(x),   x ∈Ω_i where c_i ∈ and η_i(x) is assumed to follow a certain distribution independently in i ∈Λ with a specified parameter σ∈.
The regularization is designed to penalize the boundary length of region Ω_i that is preferred to have a smooth boundary. 
The classical form of the data fidelity and the regularization is:


    𝒟 = ∫_Ωχ_i(x)   | f(x) - c_i |^p  x,  ℛ = ∫_Ω | ∇χ_i(x) |  x


where p is given depending on the noise distribution assumption (e.g. p=1 for Laplacian and p=2 for Gaussian).
The control parameter λ in (<ref>) is related to the parameter σ of the noise distribution leading to be constant for all i since σ is assumed to be the same for all i.
The energy formulation in (<ref>) with the terms defined in (<ref>) is non-convex due to the integer constraint of the characteristic function χ_i(x) ∈{0, 1}, and the control parameter λ is given to be constant in Ω for all i due to the assumption that the associated parameter σ with the noise distribution is constant.
We present a convex energy formulation with a novel constraint on the set of partitioning functions in the following section.




 §.§ Energy with Mutually Exclusive Constraint
 

We derive a convex representation of the energy functional in (<ref>) following classical convex relaxation methods presented in <cit.> as follows:


    ∑_i ∈Λ{∫_Ωλ ρ(u_i(x))  x + ∫_Ω (1 - λ)  γ( ∇ u_i(x))  x }, 
       subject to  u_i(x) ∈ [0, 1],  ∑_i ∈Λ u_i(x) = 1,  ∀ x ∈Ω,


where ρ(u_i) represents the data fidelity, γ(∇ u_i) represents the regularization, and their relative weight is determined by λ about which will be discussed in the following section.  
The characteristic function χ_i in (<ref>) is replaced by a continuous function u_i ∈ BV(Ω) of bounded variation and its integer constraint is relaxed into the convex set u_i ∈ [0, 1]. 
In the determination of the energy functional in (<ref>), we employ a robust penalty estimator using the Huber loss function ϕ_η with a threshold parameter η > 0 as defined in <cit.>:


    ϕ_η(x) =
       1/2 η x^2     : |x| ≤η,
    
       |x| - η/2    : |x| > η.


We define the data fidelity ρ(u_i) and the regularization γ(∇ u_i) using the Huber loss function as follows:


    ρ(u_i(x); c_i)    ϕ_η (f(x) - c_i)   u_i(x), 
    γ(∇ u_i(x))    ϕ_μ (∇ u_i(x)),


where c_i ∈ is an approximate of f to estimate within the region Ω_i, and the threshold parameters η, μ > 0 are related to the selection of the distribution model depending on the residual. 
The data fidelity ρ(u_i; c_i) is defined by following the piecewise image constant model, however it can be generalized to ρ(u_i; c_i)  - log p_i(f(x)) where p_i(f(x)) is the probability that the observation f(x) fits a certain model distribution p_i. 
The advantage of using the Huber loss in comparison to the L_2 norm is that geometric features such as edges are better preserved while it has continuous derivatives whereas the L_1 norm is not differentiable leading to staircase artifacts.
In addition, the Huber loss enables efficient convex optimization algorithm due to its equivalence to the proximal operator of L_1 norm, which will be discussed in Sec. <ref>. 

The  regions are desired to be pairwise disjoint, namely Ω_i ∩Ω_j = ∅ if i ≠ j, however the summation constraint ∑_i ∈Λ u_i(x) = 1,  ∀ x ∈Ω in (<ref>) is ineffective to this end, especially  with a large number of labels <cit.>.
Thus, we introduce a novel constraint to penalize the common area of each pair of combinations in regions Ω_i in such a way that ∑_i ≠ j u_i u_j is minimized for all i, j ∈Λ. Then, we arrive at the following energy functional with the proposed mutually exclusive constraint:


    ∑_i ∈Λ{∫_Ωλ ρ(u_i(x); c_i) + τ( ∑_i ≠ j u_j(x) ) u_i(x)  x 
       + ∫_Ω (1 - λ)  γ( ∇ u_i(x))  x }, 
       subject to  u_i(x) ∈ [0, 1],  ∑_i ∈Λ u_i(x) = 1,  ∀ x ∈Ω,


where τ > 0 is a weighting parameter for the constraint of the mutual exclusivity in segmenting regions, and λ determines the trade-off between the data fidelity and the regularization. 
The optimal partitioning functions u_i and the approximates c_i are computed by minimizing the energy functional in (<ref>) in an alternating way. The desired segmentation results are obtained by the optimal set of partitioning functions u_i as given by:


    l(x) = _i u_i(x),    i ∈Λ, x ∈Ω.


In the optimization of the energy functional in (<ref>), we propose a novel regularization scheme that is adaptively applied based on the local fit of data to the model for each label as discussed in the following section.




 §.§ Residual-Driven Regularization



The trade-off between the data fidelity ρ(u_i; c_i) and the regularization γ(∇ u_i) in (<ref>) is determined by λ based on the noise distribution in the image model.
We assume that the diversity parameter in the probability density function of the residual that is defined by the difference between the observed and predicted values varies in label.
We propose a novel regularization scheme based on the adaptive weighting function λ_i depending on the data fidelity ρ(u_i; c_i) as defined by:


    ν_i(x)    = exp( - ρ(u_i(x); c_i)/β), 
    λ_i(x)    = min_λ1/2ν_i(x) - λ_2^2 + α λ_1,


where β > 0 is a constant parameter that is related to the variation of the residual, and α > 0 is a constant parameter for the sparsity regularization.
The relative weight λ_i(x) between the data fidelity and the regularization is adaptively applied for each label i ∈Λ and space x ∈Ω depending on ν_i(x) determining the local fit of data to the model.
The adaptive regularity scheme based on the weighting function λ_i(x) is designed so that regularization is stronger when the residual is large, equivalently ν_i(x) is small, and weaker when the residual is small, equivalently ν_i(x) is large, during the energy optimization process.
We impose sparsity on the exponential measure of the negative residual ν_i(x) to obtain the weighting function λ_i(x) as a solution of the Lasso problem <cit.> defined in (<ref>).
Such a solution with the identity operator as a predictor matrix can be efficiently obtained by the soft shrinkage operator 𝒯(ν | α) <cit.>:

    𝒯 (ν  |  α)    = 
    ν - α    : ν > α
    
    0     : ν_1 ≤α
    ν + α    : ν < -α

leading to the solution λ_i(x) = 𝒯(ν_i(x) | α).
The Lagrange multiplier α > 0 in the Lasso problem in (<ref>) where 0 < ν_i ≤ 1 restricts the range of λ_i to be [0, 1-α] so that the regularization is imposed everywhere, which leads to well-posedness even if ρ(u_i; c_i) = 0.
The constant α is related to the overall regularity on the entire domain.
The final energy functional for our multi-label segmentation problem with the adaptive regularity scheme reads:


    ∑_i ∈Λ{∫_Ωλ_i(x)  ρ(u_i(x); c_i) + τ( ∑_i ≠ j u_j(x) ) u_i(x)  x 
       + ∫_Ω (1 - λ_i(x))  γ( ∇ u_i(x))  x }, 
       subject to  u_i(x) ∈ [0, 1],  ∑_i ∈Λ u_i(x) = 1,  ∀ x ∈Ω,


where λ_i is obtained as the solution of (<ref>).
The optimization algorithm to minimize the energy functional in (<ref>) is presented in Sec. <ref> and the supplementary material.





 §.§ Information-Theoretic Interpretation



The energy functional for our basic model incorporating the adaptive regularization, simplified after removing auxiliary variables used in the optimization and the additional constraint on the mutual exclusivity of regions, consists of a point-wise sum, which could be interpreted probabilistically by assuming that the image f is a sample from an IID random field. Under that interpretation, we have that ρ(u_i(x)) = - log p(u(x) | f(x), c_i),  and γ(∇ u_i(x)) = - log p(u(x) | c_i). For simplicity, we indicate those as -log p(u | f) and -log p(u) respectively, and even further p ≐ p(u|f) and q ≐ p(u). Then  λ_i(x) ∝ p, and the overall cost function (<ref>) to be minimized, ∫λρ + (1 - λ) γ dx, can be written as




    ℰ^ min_p,q = - 𝔼(p  logp/q) + 𝔼(log q),

where the expectation 𝔼 is the sum with respect to the values taken by u and f on x. The first term is the Kullback-Leibler divergence between the chosen prior and the data-dependent posterior. The second term is constant once the prior is chosen. Therefore, the model chosen for the posterior is the one that maximizes the divergence between prior and posterior,



where the influence of the prior q wanes as the solution becomes an increasingly better fit of the data, without the need for manual tuning of the annealing, and without the undesirable bias of the prior on the final solution. The model chosen is therefore the one that, for a fixed prior, selects the posterior to be as divergent as possible, so as to make the data as informative as possible (in the sense of uncertainty reduction).

Compare this with the usual Bayesian interpretation of variational segmentation, whereby the function to be maximized is 

    ℱ^ min_p,q = ∫log p + β log q  x,

for some fixed β and a prior q whose influence does not change with the data. If we wanted it to change, we would have to introduce an annealing parameter λ, so 

    ℱ^ max_p,q = p^λ  q^(1-λ),

with no guidance from Bayesian theory on how to choose it or schedule it. Clearly choosing λ = p yields a form that is not easily interpreted within Bayesian inference. Thus the information-theoretic setting provides us with guidance on how to choose the parameter λ, whereby the cost function is given as:

    ℰ^ max_p,q = 𝕂𝕃(p   ||   q),


where 𝕂𝕃 denotes the Kullback-Leibler divergence. We obtain λ = p and consequently an automatic, data-driven annealing schedule.




§ ENERGY OPTIMIZATION
 

In this section, we present an efficient convex optimization algorithm in the framework of alternating direction method of multipliers (ADMM) <cit.>. The detailed derivations of our optimization algorithm are further provided in the supplementary material.

The energy (<ref>), that is convex in u_i with fixed c_i, is minimized with respect to the partitioning functions u_i and the intensity estimates c_i in an alternating way. 
We modify the energy functional in (<ref>) using variable splitting <cit.> introducing a new variable v_i with the constraint u_i = v_i as follows:


    ∑_i ∈Λ{∫_Ωλ_i  ρ(u_i; c_i) + τ( ∑_i ≠ j u_j ) u_i  x 
       + ∫_Ω (1 - λ_i)  γ( ∇ v_i )  x + θ/2 u_i - v_i + y_i _2^2 }, 
       subject to  u_i(x) ≥ 0,  ∑_i ∈Λ v_i(x) = 1,  ∀ x ∈Ω,


where y_i is a dual variable for each equality constraint u_i = v_i, and θ > 0 is a scalar augmentation parameter.
The original constraints u_i ∈ [0, 1] and ∑_i u_i = 1 in (<ref>) are decomposed into the simpler constraints u_i ≥ 0 and ∑_i v_i = 1 in (<ref>) by variable splitting u_i = v_i.
In the computation of the data fidelity and the regularization, we employ a robust estimator using the Huber loss function defined in (<ref>). 
An efficient procedure can be performed to minimize the Huber loss function ϕ_η following the equivalence property of Moreau-Yosida regularization of a non-smooth function | · | as defined by <cit.>:


    ϕ_η(x) = inf_r { | r | + 1/2 η (x - r)^2 },


which replaces the data fidelity ρ(u_i; c_i) and the regularization γ(∇ v_i) in (<ref>) with the regularized forms ρ(u_i; c_i, r_i) and γ(∇ v_i; z_i), respectively as follows:


    ρ(u_i; c_i, r_i)    = ( | r_i | + 1/2 η (f - c_i - r_i)^2 ) u_i, 
    γ(∇ v_i; z_i)    =   z_i _1 + 1/2 μ∇ v_i - z_i ^2_2,


where r_i and z_i are the auxiliary variables to be minimized in alternation.
The constraints on u_i and v_i in (<ref>) can be represented by the indicator function δ_A(x) of a set A defined by:


    δ_A(x)    = 
       
       0     : x ∈ A
    ∞    : x ∉ A


The constraint u_i ≥ 0 is given by δ_A(u_i) where A = { x | x ≥ 0 }, and the constraint ∑_i v_i = 1 is given by δ_B( { v_i } ) where B = {{x_i} | ∑_i x_i = 1 }.
The Moreau-Yosida regularization for the Huber loss function and the constraints represented by the indicator functions lead to the following unconstrained energy functional ℒ_i for label i:


    ℒ_i = ∫_Ωλ_i  ρ(u_i; c_i, r_i) + τ( ∑_i ≠ j u_j ) u_i  x + δ_A(u_i) 
       + ∫_Ω (1 - λ_i)  γ( ∇ v_i; z_i )  x + θ/2 u_i - v_i + y_i _2^2,


and the final unconstrained energy functional ℰ reads:


    ℰ({u_i, v_i, y_i, c_i, r_i, z_i}) = ∑_i ∈Λℒ_i + δ_B({v_i}).


The optimal set of partitioning functions { u_i } is obtained by minimizing the energy functional ℰ using ADMM, minimizing the augmented Lagrangian ℒ_i in (<ref>) with respect to the variables u_i, v_i, c_i, r_i, z_i, and applying a gradient ascent scheme with respect to the dual variable y_i followed by the update of the weighting function λ_i and the projection of { v_i } onto the set B in (<ref>). 
The alternating optimization steps using ADMM are presented in Algorithm <ref>, where k is the iteration counter. 








The technical details regarding the optimality conditions and the optimal solutions for the update of variable at each step in Algorithm <ref> are provided in the supplementary material.

































































































§ EXPERIMENTAL RESULTS
 

We demonstrate the effectiveness and robustness of our proposed algorithm for the multi-label segmentation problem using the images in the Berkeley segmentation dataset <cit.> and synthetic images.
The numerical experiments are designed to demonstrate the robustness of the energy functional that uses a Huber function as a penalty estimator, the effectiveness of the mutually exclusive constraint on the partitioning functions, and the effectiveness of the adaptive regularity scheme based on the local fit of data to the model. Note that we use a random initialization for the initial labeling function for all the algorithms throughout the experiments.

Robustness of the Huber-Huber (H^2) Model:





87pt




We empirically demonstrate the robustness of the proposed energy functional that uses a robust Huber loss function for both the data fidelity and the regularization.
We consider a bi-partitioning segmentation on the images that are suited for a bi-partitioning image model in the Berkeley dataset in order to demonstrate the robustness of the proposed H^2 model in comparison to TV-L_1 and TV-L_2 ignoring the effect of the constraint on the common area of the pairwise region combinations. 
We apply the ADMM optimization algorithm for our model with a variable splitting technique and the primal-dual algorithm is applied for TV-L_1 and TV-L_2 models.
It is shown that our model yields better and faster results as shown in Fig. <ref> where (a) F-measure and (b) error are presented for each iteration. The parameters for each algorithm are chosen fairly so accuracy and the convergence rate are optimized. The experimental results indicate that the proposed H^2 model with the presented optimization algorithm has a potential to be applied for a variety of imaging tasks.   

Effectiveness of Mutually Exclusive Constraint:







52pt
70pt
5pt




We demonstrate the effectiveness of our proposed model with a constraint of the mutual exclusivity in comparison with other multi-label segmentation algorithms.
We qualitatively compare the segmentation results with different number of labels on the classical junction test (cf. Fig. 12 of <cit.> or Fig. 5-8 of <cit.>), whereby the number of labels is fixed to one-less than the number of regions in the input image. The algorithm is then forced to “inpaint” the central disc with labels of surrounding regions. 
The segmentation results on the junction prototype images with different number of regions are shown in Fig. <ref> where the input junction images have 5 (top), 7 (middle), 9 (bottom) regions as shown in (a). 
We compare (f) our Huber-Huber (H^2) model without the mutual exclusivity constraint in (<ref>) and (g) our full H^2 model with the constraint in (<ref>) to the algorithms including: (b) fast-label (FL) <cit.>, (c) convex relaxation based on Total Variation using the primal-dual (TV) <cit.>, (d) vectorial Total Variation using the Dogulas-Rachford (VTV) <cit.>, (e) paired calibration (PC) <cit.>.

This experiment is particularly designed to demonstrate a need for the constraint of the mutual exclusivity, thus the input images are made to be suited for a precise piecewise constant model so that the underlying image model of the algorithm under comparison is relevant. The illustrative results are presented in Fig. <ref> where the performance of most algorithms degrades as the number of regions increases (top-to-bottom), while our algorithm yields consistently better results.

Effectiveness of Adaptive Regularity:




We empirically demonstrate the effectiveness of our proposed adaptive regularization using a representative synthetic image with four regions, each exhibiting spatial statistics of different dispersion, in Fig. <ref> (a). The artificial noises are added to the white background, the red rectangle on the left, the green rectangle on the middle, and the blue rectangle on the right with increasing degree of noises in order.

To preserve sharp boundaries, one has to manually choose a large-enough λ so that regularization is small; however, large intensity variance in some of the data yields undesirably irregular boundaries between regions, with red and blue scattered throughout the middle and right rectangles (d), all of which however have sharp corners. On the other hand, to ensure homogeneity of the regions, one has to crank up regularization (small λ), resulting in a biased final solution where corners are rounded (e), even for regions that would allow fine-scale boundary determination (red). Our approach with the adaptive regularization (b), however, naturally finds a solution with a sharp boundary where the data term supports it (red), and let the regularizer weigh-in on the solution when the data is more uncertain (blue). The zoom in images for the marked regions in (b) and (e) are shown in (c) and (f), respectively in order to highlight the geometric property of the solution around the corners.



75pt
10pt




Multi-Label Segmentation on Real Images:





68pt
1pt

exp5
62096
4

124084
3

15062
4

238011
3

12003
3

29030
5





We perform a comparative analysis for the multi-label segmentation based on the Berkeley dataset <cit.>.
We compare our algorithm to the existing state-of-the-art techniques of which the underlying model assumes the piecewise constant image for fair comparison, and consider the algorithms: FL <cit.>, TV <cit.>, VTV <cit.>, PC <cit.>. 






 








We provide the qualitative evaluation in Fig. <ref> where the input images are shown in (a) and the segmentation results are shown in (b)-(f) where the same number of labels is applied to all the algorithms for each image.
The algorithm parameters for the algorithms under comparison are optimized with respect to the accuracy while we set the parameters: η=0.5, μ=0.5, α=0.01, β=10, τ=0.5, θ=1.
While our method yields better labels than the others, the obtained results may seem to be imperfect in general, which is due to the limitation of the underlying image model in particular in the presence of texture or illumination changes.

The quantitative comparisons are reported in terms of precision and recall with varying number of labels in Tables <ref>, <ref>. 



The computational cost as a baseline excluding the use of special hardware (e.g. multi-core GPU/CPU) and image processing techniques (e.g. image pyramid) is provided for 481 × 321 × 3 RGB images in Table <ref>.













§ CONCLUSION


We have introduced a multi-label segmentation model that is motivated by Information Theory. The proposed model has been designed to make use of a prior, or regularization functional, that adapts to the residual during convergence of the algorithm in both space, and time/iteration. This results in a natural annealing schedule, with no adjustable parameters, whereby the influence of the prior is strongest at initialization, and wanes as the solution approaches a good fit with the data term. All this is done within an efficient convex optimization framework using the ADMM.

Our functional has been based on a Huber-Huber penalty, both for data fidelity and regularization. This is made efficient by a variable splitting, which has yielded faster and more accurate region boundaries in comparison to the TV-L_1 and TV-L_2 models. 

To make all this work for multi-label segmentation with a large number of labels, we had to introduce a novel constraint that penalized the common area of the pairwise region combinations. 

The proposed energy formulation and optimization algorithm can be naturally applied to a number of imaging tasks such as image reconstruction, motion estimation and motion segmentation. Indeed, segmentation can be improved by a more sophisticated data term integrated into our algorithm.












§ FROM BAYES-TIKHONOV TO INFORMATION-DRIVEN REGULARIZATION
 

This section offers an alternative motivation for our choice of adaptive regularization.

The cost function (<ref>), simplified after removing auxiliary variables used in the optimization, consists of a point-wise sum, which could be interpreted probabilistically by assuming that the image f is a sample from an IID random field. 

In a Bayesian setting, one would choose a prior q based on assumptions about the state of nature. For instance, assuming natural images to be piecewise smooth, one can choose q ≐ p(u(x) | c_i) ∝exp(- γ(∇ u_i(x))). Once chosen a prior, one would then choose or learn a likelihood model ℓ = p(f(x) | u(x), c_i), or equivalently a posterior p = ℓ q where, for instance, p = p(u(x) | f(x), c_i) ∝exp(-ρ(u_i(x))). The natural inference criterion in a Bayesian setting would be to maximize the posterior:

    ψ^ max_ℓ,q = ℓ q = p.

Here the prior exerts its influence with equal strength throughout the inference process, and is generally related to standard Tikhonov regularization. Instead, we would like a natural annealing process that starts with the prior as our initial model, and gradually increases the influence of the data term until, at convergence, the influence of the prior is minimized and the data drives convergence to the final solution. This could be done by an annealing parameter λ that changes from 0 to 1 according to some schedule, yielding (an homotopy class of cost functions mapping the prior to the likelihood as λ goes from 0 to 1)

    ψ^ max_ℓ,q(λ) = ℓ^λ q^(1-λ)≠ p

Ideally, we would like to not pick this parameter by hand, but instead have the data guide the influence of the prior during the optimization. Bayesian inference here not only does not provide guidance on how to choose the annealing schedule, but it does not even allow a changing schedule, for that would not be compatible with a maximum a-posteriori inference criterion. 

Therefore, we adopt an information-theoretic approach: We first choose the prior q, just like in the Bayesian case, but then we seek for a likelihood ℓ, or equivalently a posterior p, that makes the data as informative as possible. This is accomplished by choosing the posterior p to be as divergent as possible from the prior, in the sense of Kullback-Leibler. That way, the data is as informative as possible (if the posterior was chosen to be the same as the prior, the system would not even look at the data).  This yields a natural cost criterion

    ϕ^ max_p,q = 𝕂𝕃(p || q) +  const.

choosing the constant to be 𝔼(log q), we obtain the equivalent cost function to be minimized:

    ϕ^ min_p,q = - 𝔼(p logp/q) + 𝔼(log q)

where the expectation is with respect to the variability of the data sampled on the spatial domain. If such samples are assumed IID, we have that 

    ϕ^ min_p,q = - logψ_ℓ, q^ max(p)

from (<ref>). Thus the information-theoretic approach suggests a scheduling of the annealing parameter λ that is data-driven, and proportional to the posterior p.

This yields a model with adaptive regularization that automatically adjusts during the optimization: The influence of the prior q wanes as the solution becomes an increasingly better fit of the data, without the undesirable bias of the prior on the final solution (see Fig. <ref>). At the same time, it benefits from heavy regularization at the outset, turning standard regularization into a form of relaxation, annealing, or homotopy continuation method. 

In practice, the most tangible advantages of this choice are an algorithm that is easy to tune, and never again having to answer the question: “how did you choose λ?”






§ DETAILED DERIVATION OF ENERGY OPTIMIZATION
 

The energy (<ref>) is minimized with respect to the partitioning functions u_i and the intensity estimates c_i in an alternating way. Since it is convex in u_i with fixed c_i, we use an efficient convex optimization algorithm in the framework of alternating-direction method of multipliers (ADMM) <cit.>.
We modify the energy functional in (<ref>) using variable splitting <cit.> introducing a new variable v_i with the constraint u_i = v_i as follows:


    ∑_i ∈Λ{∫_Ωλ_i  ρ(u_i; c_i) + τ( ∑_i ≠ j u_j ) u_i  x 
       + ∫_Ω (1 - λ_i)  γ( ∇ v_i )  x + θ/2 u_i - v_i + y_i _2^2 }, 
       subject to  u_i(x) ≥ 0,  ∑_i ∈Λ v_i(x) = 1,  ∀ x ∈Ω,


where y_i is a dual variable for each equality constraint u_i = v_i, and θ > 0 is a scalar augmentation parameter.
The original constraints u_i ∈ [0, 1] and ∑_i u_i = 1 in (<ref>) are decomposed into the simpler constraints u_i ≥ 0 and ∑_i v_i = 1 in (<ref>) by variable splitting u_i = v_i.
In the computation of the data fidelity and the regularization, we employ a robust estimator using the Huber loss function defined in (<ref>). 
An efficient procedure can be performed to minimize the Huber loss function ϕ_η following the equivalence property of Moreau-Yosida regularization of a non-smooth function | · | as defined by <cit.>:


    ϕ_η(x) = inf_r { | r | + 1/2 η (x - r)^2 },


which replaces the data fidelity ρ(u_i; c_i) and the regularization γ(∇ v_i) in (<ref>) with the regularized forms ρ(u_i; c_i, r_i) and γ(∇ v_i; z_i), respectively as follows:


    ρ(u_i; c_i, r_i)    = ( | r_i | + 1/2 η (f - c_i - r_i)^2 ) u_i, 
    γ(∇ v_i; z_i)    =   z_i _1 + 1/2 μ∇ v_i - z_i ^2_2,


where r_i and z_i are the auxiliary variables to be minimized in alternation.
The constraints on u_i and v_i in (<ref>) can be represented by the indicator function δ_A(x) of a set A defined by:


    δ_A(x)    = 
       
       0     : x ∈ A
    ∞    : x ∉ A


The constraint u_i ≥ 0 is given by δ_A(u_i) where A = { x | x ≥ 0 }, and the constraint ∑_i v_i = 1 is given by δ_B( { v_i } ) where B = {{x_i} | ∑_i x_i = 1 }.
The Moreau-Yosida regularization for the Huber loss function and the constraints represented by the indicator functions lead to the following unconstrained energy functional ℒ_i for label i:


    ℒ_i = ∫_Ωλ_i  ρ(u_i; c_i, r_i) + τ( ∑_i ≠ j u_j ) u_i  x  + δ_A(u_i)
       + ∫_Ω (1 - λ_i)  γ( ∇ v_i; z_i )  x + θ/2 u_i - v_i + y_i _2^2,


and the final unconstrained energy functional ℰ reads:


    ℰ({u_i, v_i, y_i, c_i, r_i, z_i}) = ∑_i ∈Λℒ_i + δ_B({v_i}).


The optimal set of partitioning functions { u_i } is obtained by minimizing the energy functional ℰ using ADMM, minimizing the augmented Lagrangian ℒ_i in (<ref>) with respect to the variables u_i, v_i, c_i, r_i, z_i, and applying a gradient ascent scheme with respect to the dual variable y_i followed by the update of the weighting function λ_i and the projection of { v_i } onto the set B in (<ref>). 
The alternating optimization steps using ADMM are presented in Algorithm <ref>, where k is the iteration counter.







The update of the estimate c_i^k+1 in (<ref>) is obtained by the following step:


    c_i^k+1∫_Ωλ_i^k (f - r_i^k) u_i^k  x/∫_Ωλ_i^k u_i^k  x.

The update for the auxiliary variable r_i^k+1 in (<ref>) is obtained by the following optimality condition:

    0    ∈∂ | r_i^k+1 | - 1/η ( f - c_i^k+1 - r_i^k+1 ),

where ∂ denotes the sub-differential operator. 
The solution for the optimality condition in (<ref>) is obtained by the proximal operator of the L_1 norm <cit.> as follows:


    r_i^k+1( f - c_i^k+1 |  η  g ),


where g(x) =  x _1. 
The proximal operator (v   |  η  g) of the weighed L_1 norm η  g with a parameter η > 0 at v is defined by:


    (v   |  η  g) _x ( 1/2 x - v _2^2 + η  g(x) ).


The solution of the proximal operator of the L_1 norm is obtained by the soft shrinkage operator 𝒯 (v   |  η) defined in (<ref>).
Thus, the solutions of the proximal problem in (<ref>) is given by:


    r_i^k+1𝒯( f - c_i^k+1 |  η).


In the same way, the update of the auxiliary variable z_i^k+1 in (<ref>) is obtained by:


    z_i^k+1𝒯( ∇ v_i^k  |  μ).


For the primal variable u_i^k+1 in (<ref>), we employ the intermediate solution ũ_i^k+1 and its optimality condition is given by:


    0 ∈λ_i^k+1 d_i^k+1 + τ∑_i ≠ j u_j^k + θ (ũ_i^k+1 - v_i^k + y_i^k), 
    
    d_i^k+1 | r_i^k+1 | + 1/2 η (f - c_i^k+1 - r_i^k+1)^2,

leading to the following update:

    ũ_i^k+1 v_i^k - y_i^k - λ_i^k+1/θ d_i^k+1 - τ/θ∑_i ≠ j u_j^k.


Given the intermediate solution ũ_i^k+1, the positivity constraint is imposed to obtain the solution for the update of u_i^k+1 as follows:


    u_i^k+1Π_A( ũ_i^k+1 ) = max{ 0, ũ_i^k+1},


where the orthogonal projection operator Π_A on a set A = { x   |   x ≥ 0 } is defined by:


    Π_A(x) = min_y ∈ A y - x _2.


We also employ the intermediate solution ṽ_i^k+1 in (<ref>) for the update of the primal variable v_i^k+1 in (<ref>). 
The optimality condition for the update of ṽ_i^k+1 reads:


    0 ∈1 - λ_i^k+1/μ∇^* ( ∇ṽ_i^k+1 - z_i^k+1 ) - θ( u_i^k+1 - ṽ_i^k+1 + y_i^k ),


where ∇^* denotes the adjoint operator of ∇, leading to the following linear system of equation with ξ_i = 1 - λ_i^k+1/μθ:


    ṽ_i^k+1 -  ξ_i Δṽ_i^k+1 = u_i^k+1 + y_i^k - ξ_i z_i^k+1,


where - ∇^* ∇ = Δ is the Laplacian operator, and - ∇^* = is the divergence operator. 
We use the Gauss-Seidel iterations to solve the linear system in (<ref>). 
Given the set of intermediate solution {ṽ_i^k+1}, the solution for the update of the variable v_i^k+1 is obtained by the orthogonal projection of the intermediate solution to the set B in Eq. (<ref>):


    v_i^k+1 = ṽ_i^k+1 - 1/n( V - 1 ),   
    V = ∑_i ∈Λṽ_i^k+1.


The gradient ascent scheme is applied to update the dual variable y_i^k+1 in (<ref>).




ieee




