PDNet: Semantic Segmentation integrated with a Primal-Dual Network for Document binarization
Kalyan Ram Ayyalasomayajula Filip Malmberg Anders Brun
stat.ML cs.LG



empty
                                                             


empty






empty



















1]Kalyan Ram Ayyalasomayajulacor1 
[cor1]Corresponding author: 
  
kalyan.ram@it.uu.se
1]Filip Malmberg
1]Anders Brun

[1]Division of Visual Information and Interaction, Dept. of Information Technology, Uppsala University, Uppsala, 751 05, Sweden.


1 May 2013
10 May 2013
13 May 2013
15 May 2013
S. Sarkar



Binarization of digital documents is the task of classifying each pixel in an image of the document as belonging to the background (parchment/paper) or foreground (text/ink). Historical documents are often subjected to degradations, that make the task challenging. In the current work a deep neural network architecture is proposed that combines a fully convolutional network with an unrolled primal-dual network that can be trained end-to-end to achieve state of the art binarization on four out of seven datasets. Document binarization is formulated as an energy minimization problem. A fully convolutional neural network is trained for semantic segmentation of pixels that provides labeling cost associated with each pixel. This cost estimate is refined along the edges to compensate for any over or under estimation of the foreground class using a primal-dual approach. We provide necessary overview on proximal operator that facilitates theoretical underpinning required to train a primal-dual network using a gradient descent algorithm. Numerical instabilities encountered due to the recurrent nature of primal-dual approach are handled. We provide experimental results on document binarization competition dataset along with network changes and hyperparameter tuning required for stability and performance of the network. The network when pre-trained on synthetic dataset performs better as per the competition metrics. 





Binarization, Semantic Segmentation, Convolutional Neural Networks, Energy Minimization, Primal-Dual scheme.
 











§ INTRODUCTION



The process of binarizing digital documents deals with classifying each pixel as belonging to the background (parchment/paper) or foreground (text/ink) while preserving most of the relevant visual information in the image. Binarization is a common pre-processing step in most tasks performed on document images, such as word spotting and transcription where a high-quality and accurate binarization significantly simplifies the task at hand. In addition to the challenges due to uneven illumination and artifacts introduced by capturing devices, historical documents may have other degradations such as; bleed through; fading or paling of the ink in some areas; smudges, stains and blots covering the text; textured background and handwritten documents with heavy-feeble pen strokes for cursive or calligraphic effects to name a few. In general, this makes the task of document binarization very challenging as shown in Fig.<ref>. The task is often subjective with multiple acceptable outcomes encountered in corner cases such as considering ink blot as being part of foreground or background. The problem of document binarization garners interest in the field, which has led to the document image binarization content (DIBCO) <cit.>, for automatic methods with minimum parameters to tune.







The task of document binarization borrows techniques from denoising, background removal, image segmentation and image in-painting, hence there exist several successful methods with individual strengths. The classical approaches have tried to separate the pixels into two classes using a single global threshold or series of finer local thresholds. The approach by <cit.>, which tries to maximize the gray level separation between foreground (FG) and background (BG) classes by maximizing the inter class variance to separate the classes. However, local intensity variations and other artifacts introduced when creating a digital image have led to  the success of locally adaptive techniques, such as the methods from <cit.>, <cit.>. The techniques discussed in all these classical methods are generic and applicable to any image in general. However, developing an approach specific to document images has been the trend in winning entries of DIBCO in the past. These methods seek to improve binarization through modeling properties of FG/BG in documents images specifically. <cit.>, have for instance modeled background using polynomial smoothing followed by local thresholding on detected text strokes, <cit.>, iteratively grow FG and BG within a 7 × 7 window.


The recent success of deep neural methods in vision related tasks has been broadly due to their ability to effectively encode the spatial dependencies and redundancy in an image. A fully convolutional neural network (FCNN) <cit.>, is best suited for semantic labeling of pixels, which is the primary objective in segmentation. The crucial idea is to use skip connections to combine the coarse features from deep layers with fine features from shallow layers to improve the final segmentation. Training such models on text images, however, often result in loss of finer details along edges. Hence a post processing step such as a graph-cut <cit.>, often improves the results as shown in our previous work <cit.>. Here, we improve upon our previous work by incorporating the energy minimization step directly in the network to facilitate joint end-to-end training of both the semantic segmentation and energy minimization steps. To this end, we adopt a primal-dual update (𝒫𝒟Update scheme <cit.>). See Fig.<ref>(a). This framework helps in training the unary cost associated with pixel labeling, pairwise cost <cit.> associated with smoothness of neighboring pixels and the cost of overlooking an edge when merging regions into a single framework, resulting in an optimally segmented image. Our contributions in the present method can be summarized as:




 
  * A stable framework that allows end-to-end training of an energy minimization function along with a semantic labeling network termed as Primal-Dual Net (PDNet).
 
  * Improved segmentation output from a semantic labeling network that is lightweight in terms of trainable weights.
 
  * A numerically stable, unrolled 𝒫𝒟Update scheme when formulating binarization as a total-variation problem that can be extended to generic image based segmentation with multiple classes.
 
  * Improved gradient propagation from 𝒫𝒟Update, using modified class weighting in loss function.


This paper is divided into five main sections describing the approach and an appendix that provides a summary of the necessary theory. An overview of document binarization with relevant background is covered in the introduction section. This is followed by a section reviewing some of the related work, highlighting some of the challenges that were overcome in our approach. An overview of the network architecture and functionality of its basic blocks are discussed in the methodology section, covering all the details on architectural changes made in building the network. This is followed by an experimental section that covers the results on DIBCO dataset for the architectural choices discussed previously. The article is concluded with contributions in the current work and possible directions for future research. All the necessary details on the theoretical framework is covered in an appendix towards the end.



§ RELATED WORK

 
The basic algorithm proposed in this paper draws motivations from other ideas that have used a high level loss function as an energy associated with binarization. This loss is then optimized by minimizing the said energy over the image; typical examples include the use of Markov random fields (MRF) for binarization <cit.> and intensity variation reduction using Laplacian kernel <cit.>. These methods take both the global and the local aspects of the image into consideration to label the pixels. The former uses the Laplacian of the image to obtain invariance in BG intensity, followed by a graph-cut with suitable source-sink priors (seed points for FG-BG, respectively) and edge estimates required to build an image graph. Although the fundamental idea of using a defined loss as employed in Howe's method has been explored previously as separate methods, combining them into an energy function proved particularly effective. Further improvement of Howe's approach was proposed in our previous method <cit.> by defining a 3D binarization space comprising of the intensity, horizontal and vertical derivatives at a pixel. A hierarchical clustering exploiting the inherent topology in this space, led to effective detection of seeds for the source and sink estimates and refinement of edges to improve the binarization result.


Before proceeding with the details in the method, we would like to motivate the reader towards an end-to-end trainable model as proposed in the current approach. The methods winning DIBCO competitions, 2009-2016 tune their parameters by training on the labeled data available from the previous competitions. Though this approach is common to all the supervised training methods in learning based approaches, these binarization methods fail to support their claim to generalizability through an exhaustive k-fold cross validation where the model is trained on all the labeled data from other competitions years and testing on data for the DIBCO competition year under consideration. We do compare our results with these state-of-the-art methods in Table.<ref> however, it must be noted that these DIBCO competition winners are different for each year (<cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>).


As the current approach is based on a deep neural network we focus on FCNNs based approaches for further discussion. To our knowledge there are two methods based on FCNNs that provide an exhaustive cross-validation results on each of the DIBCO datasets for document binarization namely, <cit.> and <cit.>. We used a post-processing based on graphcut in <cit.>, to improve the segmentation output of a FCNN, however this approach is not end-to-end trainable. The proposed approach is an improvement over the former, as it is an end-to-end trained energy minimization approach formulated as a total-variation scheme on the FCNN output. As shown in Fig.<ref>(b) output from the later network compensates for the overestimated FG along the thick stroke boundaries and preserves the under estimated FG that is missed along the thin strokes by the former network. 

We conclude this sections by discussing some of the short comings in the other deep network proposed for document image binarization. As introduced in the previous section binarization can be solved through various formulations involving pixel classification. The approach that is relevant in all the FCNN based binarization methods is that of semantic segmentation. In the semantic segmentation literature <cit.> integrating the context information surrounding a pixel showed an improved segmentation output. The proposed algorithm and <cit.> broadly fall into this category of refining the output from a FCNN. The former uses a total-variation framework and the later uses feature augmentation along with a loss function tailored for DIBCO metrics to meet this end. The performance of <cit.> depends two aspects. First, the binarization output is obtained from an ensemble of 5 networks trained on the DIBCO data instead of single trained network. Second, using feature augmentation such as image Laplacian <cit.> or Relative Brightness feature, which effect the binarization quality. <cit.>, have shown that using an ensemble is always better than using a single trained network and <cit.> shown the advantage of feature augmentation on classification output of deep networks. Both these tricks though often used in training deep networks do not alter the fundamental segmentation output of the FCNN, instead improve the final output of segmentation after a suitable loss function. However, improving the fundamental output of FCNN is known to improve segmentation output <cit.>. PDNet is generic network for semantic segmentation based on this idea of CRF unrolling <cit.> that is aimed at improving the fundamental segmentation output from the FCNN, that achieves either state of the art or close to state of the art results on a single trained network instance instead of an ensemble.



§ METHODOLOGY







The basic architecture of the proposed end-to-end binarization network PDNet is shown in Fig.<ref>(a). The network is built of three basic blocks:



  – Unary network: This is a semantic segmentation network that is capable of classifying each pixel in a given image into respective classes. Ideally such a network is quite capable of segmenting a given image by itself. Underlying such a classification is a cost associated with labeling each pixel as a particular class. We use the efficient neural network (ENet) proposed by <cit.>. The motivation for such an architecture is presented in the following section. However, as shown in our previous work <cit.> and <cit.>, text segmentation is sensitive to edge artifacts. Instead of using this network output directly we use the output prior to a typical softmax-like classification layer as the cost term for each pixel that can be further refined to improve the segmentation result.


  – Primal-Dual Update: The design of this network is inspired from previous methods that made use of conditional random field (CRF) such as <cit.>, as a post processing layer in segmentation. This allows for a way to incorporate structural information among neighboring pixels in segmentation. However, we wanted to extend this idea further in text images, where the label propagation between neighbors should be encouraged but also restricted along the edges as shown in Fig.<ref>(b). Use of 𝒫𝒟Update scheme in tasks that involve total-variation formulation of energy function is already explored in depth super resolution <cit.>, and multi-class labeling problem <cit.>. We extend these ideas into a more stable architecture that permits end-to-end training within the intended theoretical framework of the underlying proximal operator, eliminating the exploding gradient problem due to its recurrent structure.


  – Loss function:
A typical document image has a lot of background pixels as opposed to written text, which naturally leads to class imbalance between the two classes in the training data. The loss function used in the network is a weighted Spatial Cross Entropy loss <cit.>, often used to counter any class imbalance in the training samples. However, due to the redistribution of pixel labels resulting from the total-variation regularization as a part of the 𝒫𝒟Update these weights need to be re-adjusted. We propose an empirical approach to achieve this intended outcome.




  §.§.§ ENet architecture

The ENet architecture is inspired from scene parsing CNNs based on probabilistic auto-encoders <cit.>, where two separate neural networks are combined as an encoder-decoder pair. The encoder is trained to classify an input through downsampling and the decoder is used to up sample the encoder output. The ENet architecture also tries to make use of the bottleneck module that was introduced in the ResNets <cit.>. A bottleneck structure consists of a main branch that is separated from an extension consisting of convolutions filters. These two branches are later merged using elementwise addition as shown in Fig.<ref>(b). The convolution layer , is either a regular, dilated or full convolution and if the bottleneck is downsampling then a max pooling layer is added to the main branch.




We conclude the section by discussing a few key aspects of ENet  architecture as shown in Table.<ref>. The encoder constitutes of the bottleneck sections 1-3 and sections 4,5 are part of the decoder. (indicated by BNs.l; s for section and l for layer within the section). ENet has some important architectural details that improve the speed of training keeping the parameters quite low. The projection layers do not include bias terms, to reduce the number of kernel calls and overall memory operations without effecting the accuracy. The network architecture heavily reduces the input size in the first two blocks allowing for small feature maps to be retained for further processing. This is because visual information can be highly compressed due to its inherent spatial redundancy. ENet opts for a large encoder with smaller decoder as opposed to a more symmetric design. This is motivated by the fact that encoder must be able to operate on smaller resolution data, reducing the role of the decoder to that of simple upsampling.

ENet further exploits the redundancy in convolution weights by strategic replacement of n × n convolutions with two n × 1 and 1 × n convolutions filters as discussed in <cit.> and <cit.>. Strong downsampling of feature space needs to be compensated with equally adept up-sampling.  Convolutions layer were alternated with dilated convolutions to have a wide receptive field, while at the same time avoiding overly downsampling the input. Parametric rectifier linear Units (PReLU) <cit.>, were used to learn the negative slope of non-linearities through an additional learnable parameter. Most of these aspects on receptive fields, non-linear activation functions and concerned limitation in text segmentation were raised in our previous work <cit.>, making ENet architecture worth exploring for text binarization. 





  §.§.§  Primal-Dual Update scheme

The primal-dual update is built on three basic concepts 


  – Total variation formulation for segmentation.

  – Proximal operator approach to decouple the underlying total variation and gradient operations.

  – Bregman functions based proximal operator smoothing to make proximal calculations differentiable.

as discussed in the Appendix section. The primal-dual formulation of the segmentation problem is given by

    min_ u=(u_l)_l=1^kmax_ p=(p_l)_l=1^k( ∑_l=1^k( ∇ u_lp_l + u_lf_l) ) + δ_U(u) - δ_P(p)

with the primal and dual updates given by 

    û = Π_U (u̅ - τ∇^T p - τ f ) 
       p̂ = Π_P (p̅ + σ∇ u )

where are u=(u_l)_l=1^k, p=(p_l)_l=1^k, f=(f_l)_l=1^k are the primal, dual and cost vectors for 1,⋯,k classes respectively. δ_U(u), δ_P(p) are the indicator functions for the primal and dual variables u,p corresponding to the constraints sets U,P, respectively defined in Eqs.<ref>,<ref>. The orthogonal projections on to U,P are given by Π_U, Π_P, respectively . One approach to obtain a closed form representation for the projections in Eq.<ref> satisfying the constraints implicitly, is to use the Bregman proximity functions. The updates for the p,u are given in Eqs.<ref>,<ref>,  respectively. For the segmentation result to converge, the primal and dual updates need to be iterated over. PDNet has primal-dual updates unrolled over five times. A 𝒫𝒟Update with thrice such unrolled iteration is shown in Fig.<ref>, The p̅_i^k,u̅_i^k,  and  u_sum are initialized to 0,1/k and 1/k, respectively. The overall output of the network can be interpreted as a perturbation of the unary cost using primal-dual updates to give a more controlled segmentation. The final segmentation is obtained by using a weighted cross entropy loss on the final u_sum.



  §.§.§  Loss function

A cross entropy criterion (CEC) combines the logistic-softmax over the class with classwise negative log likelihood criterion to obtain the final classification. A common problem in classification is imbalance in the samples over classes. This problem can be countered by adjusted the weights associating with each class, often estimated from the class histograms. However, due to the smoothing introduced by the 𝒫𝒟Update the weights for CEC loss estimated from the histogram alone overcompensate for the imbalance. We propose a power law over the histogram based weight calculation. The weighing used is inverse of square-root of class histograms. This power law is determined by the computing the training loss with 𝒫𝒟Update being part of the network iterated over various exponents as shown in Fig.<ref>(a).


The source code for the implementation of the PDNet used for binarization is made publicly available at <https://github.com/krayyalasomayajula/pdNet.git>



 §.§ Architectural modifications

In the subsections below, we summarize the architectural changes from the basic network blocks.


  §.§.§  E-Net architecture initial block changes

When carrying out experiments on the DIBCO dataset we experimented with both color and gray scale images as well as using both color and gray channels. Although state of the art results were obtained using gray-scale images, we would like to highlight some extreme cases as shown in Fig.<ref>(a). When using color and gray channels the initial block was modified to include RGB and two gray channels with the max pooling applied to RGB channels alone; results for this case were similar to using RBG channels alone. When training the network for RGB and gray-scale alone the original ENet architecture was used without any changes.




  §.§.§ Clamped primal-dual updates

The typical primal-dual updates when properly initialized are usually stable. However, when training a 𝒫𝒟Update scheme on image data in deep networks exploding gradient problem is commonly encountered. Implementations by <cit.> and  <cit.>, have dealt with this issue by gradient clipping to specific bounds during back-propagation. We resorted to another approach of clamping the values in primal and dual updates as shown in Fig.<ref>(b). This approach has two advantages:


  – The clamping as shown in Fig.<ref>(b) resets the pixel where instability was encountered to their initial values 0, 1/k for p̅_i^k, u̅_i^k, respectively. The cost estimates for these pixels can be refined in further iterations, thus the resulting scheme is more faithful to the theoretical primal-dual approach.

  – The gradients in this approach are not clipped during back-propagation thus leading to a faster training of the network making the loss converge within 10 epochs as opposed to 30 epochs in a network without clamping.




  §.§.§ Choice of labeling

The convention often followed in semantic segmentation is to allow an unknown class, to neglect objects such as background scenery or to include classes that need to be ignored. Generally, presence of such classes does not cause any change in the result, but do increase the training time for the PDNet. The  unknown class is usually labeled 1 and then other classes are labeled incrementally. We begin by relabeling the FG, BG and unknown classes to 1,2 and 3, respectively. Once the unary network is trained, the cost vector can be truncated to include just the FG and BG costs. Training PDNet on these costs leads to faster convergence as the size of p̅_i^k, u̅_i^k and further vectors used in computations are reduced by 30%.



 §.§ Tuning hyperparameters

The binarization network has two hyperparameters to tune: The index to be used in the power law for the class weights used in the CEC loss function and the weights to be associated with the edges along the class boundaries. As mentioned previously the exponent in the power law for the weights was determined through an exhaustive search carried out on the training error for various exponents. The weights associated with the edges is a shared parameter over all the unrolled loops and multiplied with τ weight of each loop unroll. Training for both the edge weights and τ can cause instability. We instead initialized the edge weight to 1.0 and trained on a small training-set to determine the convergence. This value was later used in the network. Although the edge weight determined by this approach is suboptimal, it is compensated for by learning τ values through back-propagation in the network over the complete training sets. This method of training for the edge weight produced a better validation loss, though this mode of pre-training is not critical for the network performance it can save time required for the over all training of the network.





§ EXPERIMENTAL RESULTS



  §.§.§ DIBCO dataset

The experiments were conducted on the DIBCO datasets by <cit.>, for binarization competition from years 2009-2016 consisting of 76 images in total. The images and ground truth images were augmented by applying an identical deformation field transformation to both. These augmented images were then converted into 128 × 256 pixels of cropped images with 25% overlap horizontally and vertically to create more data for training. The cropped size 128 × 256 was selected based on the encoder requirement for the width and height to be equal to 2^n for some n and to fit the training batches into memory in Torch7 framework <cit.>. The training set was picked by including the all the augmented images, except those from competition year under evaluation. The validation set was made by randomly picking 3000 original cropped images from the DIBCO datasets excluding the year under evaluation. The trained network was then used to produce binary output on the dataset of DIBCO competition year. The images were then stitched to the original size to compute the DIBCO evaluation metrics. Three main evaluation metrics for comparison are F-Measure, Peak signal to noise ratio (PSNR) and Distance Reciprocal Distortion Metric (DRD), definitions of which can be found in <cit.>, <cit.>.







  §.§.§ Synthetic Dataset

To have better weight initialization and prevent over fitting the network to data, the network was pre-trained on synthetic data. Documents resembling historical handwritten and printed material were generated synthetically. Various filters were applied to resemble background textures and degradations in the parchment. The text was generated using handwriting and machine printed fonts from the Google™<cit.>. Fig.<ref>(b) shows few cropped images from the synthetic dataset. The results from binarization on DIBCO dataset using the network pre-trained on the synthetic dataset are presented in Table-<ref>.



  §.§.§ Training

The network was trained in three stages: 


  – Pre-training unary network on synthetic dataset.

  – Training unary network on DIBCO dataset.

  – Combined training of unary and 𝒫𝒟Update scheme on DIBCO dataset.


The encoder of the unary network was pre-trained on patches of sizes 128×256 with batch-size 30 on the Synthetic dataset. The pre-trained encoder was then used to initialize the decoder weights and then the decoder was per-trained on the Synthetic dataset using the same patch size and batch size. This gives the pre-trained unary model. This model was then used to initialize the weights of the unary model that is to be trained on the augmented DIBCO dataset. The gradients were trained using ADAM <cit.>, with the  learning rate, weight decay, and momentum set to 5 × 10^-4, 2 × 10^-4 and  0.9, respectively. The unary models were trained for 20 epochs and the best model with least loss was picked in further steps. When training the unary and 𝒫𝒟Update the learning rate was set to 5 × 10^-4 for the first 10 epochs and decayed to 2 × 10^-4 and  1 × 10^-4 between 11 - 15 and 16 - 20 epochs, respectively. The results on DIBCO dataset, with the proposed binarization network is summarized in Table.<ref>.


PD_g, PD_c are the outputs from primal-dual networks trained on grayscale and color images, respectively. GC is the result obtained from a graph-cut based approach as discussed in our previous work <cit.>. It takes the output from the unary model as seed points along with the output from the classification layer acting as costs for pixels and Canny edges as boundaries estimates for segmentation. This approach requires three parameters a) cost associated with pixels labeling as obtained from the unary network, 2) weight associated to edges and Canny threshold to estimates the boundaries. It then employees an external graph-cut method <cit.>, to obtain the final segmentation. In contrast to this approach the current methods does all these parameter estimations and the energy minimization in end-to-end manner in a single framework. TM are the results from another CNN based approach developed by <cit.>, which augments the segmentation result with relative darkness feature <cit.>, to aid in binarization. DBC are results from winning entries in DIBCO competitions using various classical approaches in  as given in <cit.>, <cit.>.



§ CONCLUSION

In the current work we have extended our previous approach of using graph-cuts as post-processing to a FCNN, to obtain binarized documents images. In this paper, we propose an architecture that combines the use of an energy minimization function and FCNN based feature learning. The energy minimization framework is flexible in imposing constraints on the desired segmentation. This primal-dual network imposes a total variation based energy minimization, with an unrolled primal-dual update scheme in a gradient descent trainable CNN architecture. The unary CNN, on the other hand, learns pixel labeling costs and the combined network learns all the associated parameters within a single end-to-end trainable framework.



The numerical instability caused by the recurrent nature of the primal-dual steps was solved by using a clamping function on the primal and dual updates within the 𝒫𝒟Update block. Propagation of gradients from the 𝒫𝒟Update was further facilitated by employing a weighted cross entropy loss adjusted by a power law. The final primal-dual architecture improves the binarization results, compared to the FCNN output alone, even when combined post-processing such as graph-cut or feature augmentation. The binarization achieves state of the art results on four out of seven DIBCO datasets, with comparable to best result on the rest. Further investigation in dropout layers in the unary model and architectural variations in 𝒫𝒟Update are planned for future work. The results from using a trained network to other historical document images and cases of transfer learning are of practical interest. Using more sophisticated unary and primal-dual schemes could also yield an even better result on document binarization.





§ ACKNOWLEDGMENTS

This project is a part of q2b, From quill to bytes, an initiative sponsored by the Swedish Research Council "Vetenskapsrådet D.Nr 2012-5743) and Riksbankens Jubileumsfond (R.Nr NHS14-2068:1) and Uppsala university. The authors would like to thank Tomas Wilkinson of Dept. of Information Tech., Uppsala University for discussions on debugging and performance optimization in Torch framework.

model2-names




§ APPENDIX


Total variation formulation: Formulating image segmentation as a saddle point problem and applying proximal operator to the primal-dual variables is a well studied problem in fixed point analysis of convex functions <cit.>. An overview of theory is provided here to make the discussion comprehensive and to provide intuition for 𝒫𝒟Update discussed in the method. We consider segmentation of image into k-pairwise  disjoint regions as a total variation on the segmentation image as: 

    min_(R_l)_l=1^k, (c_l)_l=1^k1/2∑_l=1^k Per(R_l,Ω) + λ/2∑_l=1^k∫_R_l |g(x)-c_l|^2 dx,

where Per(R_l,Ω) is the perimeter of the region R_l in a domain Ω, g: Ω→ℝ is the input image, c_l ∈ℝ are the optimal mean values and the regions (R_l)_l=1^k form a partition of Ω that is, R_l ∩ R_m = ∅, l  m and ∪_l=1^k R_l = Ω, λ is a regularization weight. This is an optimization problem between the data fitting term |g(x)-c_l| and the length term
 Per(R_l,Ω) where the ideal mean values of the region c_l = ∫_R_l g(x)/|R_l| are unknown a-priori as they depend on the partition we seek.

 
By introducing a labeling function u = (u_l)_l=0^k: Ω→ℝ and viewing the data fitting term as a weighing function f_l = λ/2|g(x)-c_l| Eq(<ref>) can be generalized as:

    min_ u=(u_l)_l=1^k J(u) + ∑_l=1^k∫_Ω u_l f_l dx, 
     s.t    ∑_l=1^k u_l(x)=1, u_l ≥ 0, ∀ x ∈Ω,

where J(u) is the relaxation term. 


Proximal operator:
This basic structure in Eq.<ref> is of the form 

    min_x ∈ X F(Kx) + G(x),

involving a linear map K: X → Y with the induced norm ||K|| = max{ ||Kx||: x ∈ X  with  ||x|| ≤ 1 } where X,Y are the primal and dual spaces, respectively. The corresponding dual formulation of this equation is a generic saddle-point problem 

    min_x ∈ Xmax_y ∈ YKxy +G(x) - F^*(y),

where Kxy is the inner product induced by the vector space Y and F^* is the convex conjugate of F.


The advantage of such an approach is discussed further down with the segmentation example, for now we can observe that this structure readily presents a computationally tractable algorithm. The dual variable y ∈ Y acts like bounded slack variables introduced to ease the solution in the resulting dual space Y. Introducing the dual variable y relieves the initial composition of F(Kx) making computations involving Kxy independent of F^*(y). As per the structure of segmentation problem, F is typically an indicator function for capturing the constraints on x, which translates to F^* being an indicator function for capturing the constraints on its dual variable y. Since the primal space has ||x|| ≤ 1 if the dual variable is bounded, which is most often the case then iterating repeatedly between the two variables should converge to a solution.


The solution takes a form involving the proximal operator or gradient of the functions F,G depending on them being convex or convex as well as differentiable, respectively. The basic idea behind a proximal operator is a generalization of projection on to a vector space. This makes it an ideal operator that can be used in a gradient descent algorithm where the iteration involves taking a suitable step towards the solution along the gradient direction. But since the function need not be differentiable the gradient need not necessarily exist and hence the question of uniqueness along the gradient direction does not arise. This results in a set of permissible vectors that though not strictly a gradient can act as one at a given point x, such a set of permissible vectors is called subgradient. The set ∂ F is the subgradient it is also the set of underestimators of F at x. A closely related set is the resolvent operator with the property

    x = (I + τ∂ F )^-1 (y) = min_x{||x-y||^2/2τ +F(x) }.

It can be interpreted as the closest point x in the set under consideration to the given point y under an associated error F(x). The primal-dual formulation allows for an algorithm that iterates between the primal and dual variables x,y, respectively in this case leading to convergence according to the forward-backward algorithm

    y^n+1=prox_x(y^n) = (I + σ∂ F^* )^-1 (y^n + σ K x̅^n); 
       x^n+1=prox_y(x^n) = (I + τ∂ G )^-1 (x^n - τ K^* y^n+1);
       x̅^n+1 = x^n+1 + θ (x^n+1 - x^n),

where τ, σ are step lengths along dual and primal subgradients and θ is the relaxation parameter in iterating the primal variable.


Considering the image discretised over the Cartesian grid of size M × N as { (ih,jh) : 1 ≤ i ≤ M, 1 ≤ j ≤ N }, where h is the size spacing and (i,j) the indices in discrete notation. X is a vector space in ℝ^MN equipped with standard inner product uv for u,v ∈ X. The gradient is defined as ∇: X → Y, ∇ u = (u_i+1,j-u_i,j/h,u_i,j+1-u_i,j/h) with Y=X × X equipped with the inner product defined as,

    pq_Y = ∑_i,j p_i,j^1 q_i,j^1 + p_i,j^2 q_i,j^2, p=(p^1,p^2), q=(q^1,q^2) ∈ Y.


Applying the above framework to Eq.<ref> with J(u) = 1/2∑_l=1^k∫_Ω|∇ u_l|, we have

    min_ u=(u_l)_l=1^k1/2∑_l=1^k( ∫_Ω|∇ u_l| + u_lf_l) + δ_U(u),

where G(u) = δ_U(u) is the indicator function for the unit simplex,

    U = { u ∈ X^k : ∑_l=1^k u_l(x)=1, u_l ≥ 0 }.

f=(f_l)^k_l=1∈ X^k is the discretized weighting function or the cost per pixel, u=(u_l)^k_l=1 is the primal variable and X^k is the extension of the vector space for k classes.
Considering 

    1/2∑_l=1^k(∫_Ω|∇ u_l| + u_lf_l) = 
       max_ p=(p_l)_l=1^k( ∑_l=1^k( ∇ u_lp_l + u_lf_l) ) - δ_P(p),

where p ∈ Y^k is the dual variable, with Y^k is the extension of the gradient vector space for k classes, δ_P(p) is the indicator function for p ∈ P defined as 

    P = { p ∈ Y^k : ||p_l||_∞≤1/2}.

We have the primal-dual formulation as

    min_ u=(u_l)_l=1^kmax_ p=(p_l)_l=1^k( ∑_l=1^k( ∇ u_lp_l + u_lf_l) ) + δ_U(u) - δ_P(p),

with u,p related as ∇ up_Y = -udiv p_X. This result is a consequence of applying Guass Divergence theorem on the scalar function u and vector field p. A corollary of the fore mentioned result is the relation -div=∇^* where div, ∇^* are divergence in Y;  and the conjugate of gradient ∇, respectively. Further since ∇^* = -∇^T it turns out that div=∇^T.


Bregman Proximity Functions:
The sets considered so far are unit simplex and unit ball (or more strictly a ball of radius 1/2) as defined by U,P in Eqs.<ref>,<ref>, respectively. The resolvent of these sets are orthogonal projections on unit simplex and point projection onto unit ball, respectively. However, in the case of segmentation when using more sophisticated relaxation that yield better delineation along edges, like paired calibrations is used in Eq.<ref> given by,

    J(u) = ∫_ΩΨ(Du); 
    
    s.t    Ψ(a) = sup_b{∑_l=1^ka_lb_m : |a_l - b_m| ≤ 1, 1≤ l ≤ m ≤ k },

where a=(a_1,⋯,a_k), b=(b_1,⋯,b_k). The corresponding set for the dual variables is no longer a unit ball, but intersection of unit balls given by,

    P ={ p ∈ Y^k : |p_l - p_m|_∞≤ 1, 1≤ l ≤ m ≤ k }.

The resolvent of which, is an orthogonal projection on to such an intersection of unit balls. As the relaxations get more sophisticated the corresponding resolvent set becomes more complex and orthogonal projections on to them get computationally more involved. One approach towards getting a solution that can be used in a computational algorithm is to use Bregman functions. Suppose we have a convex function ψ(x) that is continuously differentiable on the interior of its domain; int(X) and continuous on its closure; cl(X), we use x̅ to denote a point from int(X). A Bregman proximity function D_ψ : X ×int(X) →ℝ generated by ψ is defined as 

    D_ψ(x,x̅) = ψ(x)-ψ(x̅) - ∇ψ(x̅)x-x̅.


In iterative algorithms, the Bregman proximity function can be used with the proximity operator for a convex function g : X →ℝ as

    prox_α g^ψ(x̅) = min_x ∈ Xα g(x) + D_ψ(x,x̅).


In image segmentation problem the basic class of functions of interest are of the form g(x) = xc + δ_X(x) as seen in Eq.<ref>. The associated proximal operator is 

    prox_α g^ψ(x̅) = min_x ∈ Xαxc + D_ψ(x,x̅).

The necessary and sufficient condition for optimality, which has a unique solution 
for Eq.<ref> is

    ∇ψ(x̅) - c = ∇ψ(x).

This constraint is implicitly taken care by the Bregman proximity function. For further details on Bregman functions one may refer <cit.>. In image segmentation the dual variables belong to the intersection of unit balls as shown in Eq.<ref> so each coordinate of the dual variable p should satisfy -1 ≤ p_j ≤ 1 and solve the dual problem

    max_ p=(p_l)_l=1^k∑_l=1^k∇ u_lp_l  - δ_P(p).

A suitable Bregman proximity functions that encode the dual variable constraints and the corresponding proximity solution along each coordinate i are given by

    ψ(x)= 1/2[ (1+x)log (1+x) + (1-x)log (1-x) ]; 
       ( prox_α g^ψ(x̅) )_i = exp( -2 α c_i ) - 1-x̅_i/1+x̅_i/exp( -2 α c_i ) + 1-x̅_i/1+x̅_i,

where c_i = (∇ u)_i is the i-component of the ∇ u. Similarly, the primal problem deals with 

    min_ u=(u_l)_l=1^ku_lf_l  + δ_U(u),

with the primal variables u restricted by u_i ≥ 0 along each coordinate. The Bregman function encoding these constraints and the corresponding proximal solution are given by

    ψ(x)= x log x, 
       ( prox_α g^ψ(x̅) )_i = x_i exp( -2 α c_i ).

To satisfy ∑_l=1^k u_l(x)=1 in Eq.<ref> it is sufficient to normalize it as 

    ( prox_α g^ψ(x̅) )_i = x_i exp( -2 α c_i )/∑_j=1^Kx_j exp( -2 α c_j ),

where c_i = (∇^T p)_i - f_i, f_i denoting the cost associated with the i^th class. Eqs.<ref> and <ref>, can now be used in Eq.<ref> to converge to a solution for image segmentation.

