Total Variation Regularization with Chambolle Algorihtm

The total variation prior is a non-smooth convex function of the image. To perform a minimization of an energy involving the total variation, one needs to use dedicated algorithm that can handle non-smooth priors.

Contents

Total Variation Prior

The total variation is a Banach norm. On the contrary to the Sobolev norm, it is able to take into account step edges.

First we load a clean image.

n = 256;
name = 'hibiscus';
f0 = load_image(name,n);
f0 = rescale( sum(f0,3) );

The total variation of a smooth image \(f\) is defined as \[J(f)=\int \|\nabla f(x)\| d x\]

It is extended to non-smooth images having step discontinuities.

The total variation of an image is also equal to the total length of its level sets. \[J(f)=\int_{-\infty}^{+\infty} L( S_t(f) ) dt\]

Where \(S_t(f)\) is the level set at \(t\) of the image \(f\) \[S_t(f)=\{ x \backslash f(x)=t \}\]

Exercice 1: (the solution is exo1.m) Compute the total variation of f0.

exo1;

The Gradient of the TV norm is \[ \text{Grad}J(f) = \text{div}\left( \frac{\nabla f}{\|\nabla f\|} \right) . \]

The gradient of the TV norm is not defined if at a pixel \(x\) one has \(\nabla f(x)=0\). This means that the TV norm is difficult to minimize, and its gradient flow is not well defined.

Total Variation Regularization

We consider here the denoising problem and its solution using variational methods.

Add some noise to the original image.

sigma = .1;
y = f0 + randn(n,n)*sigma;

We consider the solution of a variational TV regularization to denoise an image \(f_0\) from obervations \(y=f_0+w\) \[\min_{f} E(f)= \frac{1}{2} \|y-f\|^2 + \lambda J(f)\]

The parameter \(\lambda\) should be adjusted to match the noise level \(\epsilon = \|w\|\) (a larger \(\lambda\) implies a stronger denoising).

Set the regularization parameter.

lambda = .2;

The mapping f \mapsto E(f) is not smooth, and one cannot use a gradient-descent method to solve the problem.

Chambolle Dual Algorithm

The minimization is performed using the algorithm proposed by Antonin Chambolle in

The idea is to replace the optimization of the image \(f\) by the optimization of a vector field \(G\) that is related to \(f\) by \( f=y-\lambda\text{div}(G) \). The vector field is the one that minimizes \[ \min_G \|y-\lambda \text{div}(G)\| \]

Subject to the constraints that for all position \(x\), one has \( \|G_x\| \leq 1\).

Chambolle proposed to perform this minimization with a fixed point iteration. Another approach is simply to perform a projected gradient descent on \(G\).

An iteration of this projected gradient reads: \[ G^{(\ell+1)} = \text{Proj}\left( f^{(\ell)} - \tau \nabla (\text{div}(G^{(\ell)}) - y/\lambda) \right) \]

Where \(\text{Proj}\) is the orthogonal projector on the constraint \( \|G_x\| \leq 1\).

The descent step size should satisfy \[ \tau < \frac{2}{ \| \nabla \circ \text{div} \| } = \frac{1}{4}. \]

tau = 0.9/4;

We initialize the iteration with a zero vector field, that corresponds to \(G^{(0)}\).

G = zeros(n,n,2);

Compute the gradient of the functional \( \|y-\lambda \text{div}(G)\|^2 \).

dG = grad( div(G) - y/lambda );

Do a step of gradient descent.

G1 = G + tau*dG;

Do the projection to ensure that the entries are smaller than 1.

d = repmat( sqrt(sum(G1.^2,3)), [1 1 2] ); % norm of the vectors
G = G1 ./ max(d,ones(n,n,2));

Display.

clf;
imageplot( G1, 'Vector field before projection', 1, 2, 1 );
imageplot( G, 'Vector field after projection', 1, 2, 2 );

Exercice 2: (the solution is exo2.m) Perform Chambolle algorithm to find the solution the variational problem by iterating this gradient descent. Reconstruct the denoise image according to the fTV=y-lambda*div(G)

exo2;

Display the denoised image.

clf;
imageplot(clamp(y), 'Noisy image', 1,2,1);
imageplot(clamp(fTV), strcat(['TV regularized using \lambda=' num2str(lambda)]), 1,2,2);

Parameter Adaptation

The issue with the minimization is that it is difficult to set the value of \( \lambda \) for denoising. It is often simpler to solve the following constraint minimization. \[ \min_f J(f) \quad \text{subj. to} \quad \|f-y\| \leq \epsilon \]

Where \(\epsilon\) is the noise level, set to \(\epsilon = n \sigma\) where \(\sigma\) is the standard deviation of the image.

This constraint minimization can be solved by modifying the value of \(\lambda\) at each iteration of Chambolle's algorithm. If the solution, at a given step, is \( f^{(\ell)} = y-\lambda^{(\ell)} \text{div}(G^{(\ell)}) \), then the \(\lambda=\lambda^{(\ell)}\) is modified according to \[ \lambda^{(\ell+1)} = \lambda^{(\ell)} \frac{n\sigma}{ \|f-y\| } \]

Exercice 3: (the solution is exo3.m) Implement Chambolle's algorithm with this update rule for lambda. Display the evolution of lambda and norm( y-fTV,'fro')-n*sigma during the iterations.

exo3;

Display the result

clf;
imageplot(clamp(y), strcat(['Noisy, SNR=' num2str(snr(f0,y),3) 'dB']), 1,2,1 );
imageplot(clamp(fTV), strcat(['TV Regularization Denoised, SNR=' num2str(snr(f0,fTV),3) 'dB']), 1,2,2 );

Copyright (c) 2010 Gabriel Peyre