Inverse problem with quadratic data-fidelity and total variation regularity

A solution of problem \eqref{eq:primal-dual-inverse} can be numerically computed using the algorithm of Chambolle-Pock [3], which is designed to solve the generic saddle-point problem

The algorithm is briefly presented on this page.

Taking \(X=\RR^\Omega\), \(Y=\left(\RRRR{\Omega} \right) \times \RR^\omega\), \(x=u\), \(y=(p,q)\), \(F^\star(p,q) = \delta(p) + \tfrac{1}{2}\|q+u_0\|_2^2\), \(G(u) = 0\) and \(Ku = \left(\lambda \nabla u , Au \right)\) into problem \eqref{eq:cp.primal-dual} allows to recover problem \eqref{eq:primal-dual-inverse}.

We note \(\divergence = -\nabla^*\) the opposite of the adjoint of the operator \(\nabla\) (in analogy with the continuous setting) and \(A^*\) the adjoint of \(A\), we have then \(K^* = -\lambda \divergence + A^*\). The update of the dual variable (here the tuple (p, q)) in the Chambolle-Pock algorithm (see Algorithm 1 in on this page) can been split into two independant updates thanks to the additive separability with respect to \(p\) and \(q\) of \((p,q) \to \langle K u , (p,q) \rangle - F^\star(p,q)\), hence applying the Chambolle and Pock algorithm to problem \eqref{eq:cp.primal-dual} boils down to Algorithm 1.

The dual and primal updates (Equations \eqref{cp.dual-p}, \eqref{cp.dual-q} and \eqref{cp.primal}) can be computed in closed-form.

Notice that the convergence of Algorithm 1 toward the (here unique) solution of problem \eqref{eq:primal-dual-inverse} is ensured when \(\theta = 1\) and \(\tau \sigma {|||K|||}^2 < 1\). Here we have \({|||K|||}^2 \leq \lambda^2 {|||\nabla|||}^2 + {|||A|||}^2\), this upper bound depends on the choice of the operator \(A\) and the discretization scheme used for \(\nabla\).

Let \(n_x\) and \(n_y\) denote the width and height of the initial (noisy) image \(u_0\), we note \(\Omega = \left\{ 1, \dots , n_x \right\} \times \left\{ 1, \dots , n_y \right\}\) the discrete domain of \(u_0\). We propose here to use the following finite difference scheme for the gradient operator

$$ \forall u \in \RR^\Omega, \quad \forall (i,j) \in \Omega, \quad \nabla u (i,j) = \left( \nabla_x u (i,j), \nabla_y u (i,j) \right)\;, $$

where

This operator is easily implementable in Scilab langage,

function [Dx,Dy] = grad(u)
  [ny,nx] = size(u);
  Dx = u(:,[2:nx,nx]) - u;
  Dy = u([2:ny,ny],:) - u;
endfunction

and we can prove that the corresponding divergence operator (the opposite of the adjoint of \(\nabla\)) is given by

$$ \forall p = (p_x,p_y) \in \RRRR{\Omega}, \quad \forall (i,j) \in \Omega, \quad \divergence{(p)}(i,j) = \text{div}_x (p_x)(i,j) + \text{div}_y (p_y)(i,j) \;, $$

where

This operator is also easily implementable in Scilab langage.

function d = div(px,py) // compute d = div(p) where p = (px,py) and size(px) == size(py)
  [ny,nx] = size(px);
  div_x      =  px - px(:,[1,1:(nx-1)]);
  div_x(:,1) =  px(:,1);
  div_x(:,nx) = -px(:,nx-1);
  div_y      =  py - py([1,1:(ny-1)],:);
  div_y(1,:) =  py(1,:);
  div_y(ny,:) = -py(ny-1,:);
  d = div_x + div_y;
endfunction

With the choice of discretization scheme described in the previous section, we can show that the induced \(\ell^2\) norm of the operator \(\nabla\) is less than \(\sqrt{8}\) (that is \({|||\nabla|||} \leq \sqrt{8}\)). Let us note \(L_A\) an upper bound of \({|||A|||}\) (will be precised for each application proposed), we have \({|||K|||} \leq L := \sqrt{8\lambda^2 + L_A^2}\), therefore taking \(\tau = \sigma = \tfrac{0.99}{L}\) into Algorithm 1 ensures that \(\tau \sigma \, {|||K|||}^2 < 1\) (which, in addition to the setting \(\theta=1\) is sufficient to ensure the convergence of the algorithm toward the solution of problem \eqref{eq:primal-inverse}).

We propose below a Scilab implementation of Algorithm 1 that can be used to solve the inverse problem \eqref{eq:primal-inverse} once the operator \(A\), its adjoint \(A^*\) and an upper bound \(L_A\) of \(|||A|||\) are set into Scilab memory.

More precisely, the user must fill according to the intended application, the following code and load it into Scilab memory

//* define operator A *//
function v = A(u)
  //* enter your code here *//
  ...
endfunction

//* define adjoint of A *//
function u = adj_A(v)
  //* enter your code here *//
  ...
endfunction

//* define an upperbound for the induced l2 norm of A *//
L_A = ... //* enter the value here *//

optionally, the user may define an initial value for the primal variable \(u\) (initializing the algorithm with a more or less good guess of the solution may drastically accelerate the convergence). For that simply load into Scilab memory the initialization of your choice with the name init. If this init is not defined the primal variable \(u\) will be initialized to the null image (image full of zeros).

init = ... //* enter the value here *//

Once the setting is done, you can use the following implementation of Algorithm 1 to solve the inverse problem with the Chambolle and Pock algorithm.

 //          Solver for inverse problems with quadratic data-fidelity
 //                     and total variation regularity
 //
 // u0 = input image
 // lambda = regularity parameter (TV weight)
 // niter = number of iterations for the Chambolle-Pock algorithm
 // u = restored image
 // E = energy (E(i) = energy computed at iteration i)
 //
 // Warning: the bound L_A and the operators A, adj_A must have been
//           defined into Scilab memory.
 function [u,E] = inverse_solver_tvreg(u0,lambda,niter)

   E = zeros(1,niter);

   // retrieve dimensions of domains omega & Omega (recall
   // that operator A is from R^Omega to R^omega).
   [ny0,nx0] = size(u0);      // domain omega
   [ny,nx] = size(adj_A(u0)); // domain Omega

   // initialize primal and dual variables
   if(exists("init")) then
     if((size(u,1) == ny)&(size(u,2) == nx)) then u = init;
     else u = zeros(ny,nx); end
   else u = zeros(ny,nx); end
   ubar = u;
   px = zeros(u);
   py = zeros(u);
   q = zeros(ny0,nx0);

   // set Chambolle & Pock algorithm time-step parameters
   L = sqrt(8*lambda^2+L_A^2);
   tau = 0.99/L;
   sigma = 0.99/L;
   theta = 1;

   //**************************** main loop ****************************//
   for i = 1:niter
     //* compute (Dx,Dy) = grad(ubar) *//
     [Dx,Dy] = grad(ubar);
     //* update p = (px,py) *//
     px  = px + sigma*lambda*Dx;
     py  = py + sigma*lambda*Dy;
     nrm = px.*px + py.*py;
     id = nrm > 1; nrm = sqrt(nrm(id));
     px(id) = px(id) ./ nrm;
     py(id) = py(id) ./ nrm;
     //* update q *//
     q = (q + sigma*(A(ubar)-u0))/(1+sigma);
     //* compute d = div(p) *//
     d = div(px,py);
     //* update u and ubar *//
     ubar = -theta * u;
     u    = u + tau*(lambda*d-adj_A(q));
     ubar = ubar + (1+theta)*u;
     //* compute energy of u *//
     [Dx,Dy] = grad(u);
     tv = sum(sqrt(Dx.^2 + Dy.^2));
     E(i) = 0.5*sum((A(u)-u0).^2) + lambda*tv;
     printf("iteration %d : E = %.10g\n",i,E(i));
   end
   //************************* end of main loop *************************//
 endfunction

We first propose some very basic tools dedicated to image manipulation in Scilab (opening, visualization).

Image viewer

function imview(u)
  [height,width] = size(u);
  fg = figure();
  set(fg,"axes_size",[width,height],"infobar_visible","on","menubar_visible","on","toolbar_visible","off","auto_resize","off","color_map",graycolormap(256));
  set(fg.children,"margins",zeros(1,4));
  Matplot(255/(%eps + max(u)-min(u)) * (u-min(u)),"010",[1,1,width,height]);
endfunction

Image opening (format PMG)

// open a PGM RAW image (8 bits)
// usage: u = readpgm('lena.pgm');
function image=readpgm(filename)
  //* local function *//
  function l=getline(u)
    h=[]
    while %t
      c=mget(1,'uc',u)
      if c==10 then break, end
      h=[h c]
    end
    l=ascii(h)
  endfunction
  //* main *//
  [u,err]=mopen(filename,'rb')
  if err<>0 then error('Impossible to open file '+filename), end
  if getline(u)~='P5' error('Unrecognized format'), end
  z=getline(u), while part(z,1)=='#', z=getline(u), end
  execstr('n=['+z+']')
  getline(u)
  image=matrix(mget(n(1)*n(2),'uc',u),n)'
  mclose(u)
endfunction

In the case of image deconvolution, the linear operator \(A\) in \eqref{eq:primal-inverse} is the convolution with a point spread function (modeling for instance some blurring phenomena such as diffraction, defocus or motion, etc.). Notice that a convolution can only model uniform phenomena while in practice optical devices generally suffer from more complicated distortions such as chromatic aberrations, stigmatism and coma, vignetting, etc. The correction of such non-uniform distortions is therefore out of the scope of the restoration application detailed here.

Let us consider a discrete convolution kernel \(k_A \in \RR^{\omega_A}\) with domain \(\omega_A \subset \ZZ^2\) that we can assume to be symmetric (\((a,b) \in \omega_A \Leftrightarrow (-a,-b) \in \omega_A\)) without loss of generality. We define \(A : \RR^\Omega \mapsto \RR^\omega\) by

and \(\omega\) denotes the subset of \(\Omega\) made of all the pixels \((x,y) \in \Omega\) such as \((x,y)-\omega_A \subset \Omega\). Remark that it is also possible to define the convolution with kernel \(k_A\) as an operator \(A : \RR^\Omega \mapsto \RR^\Omega\) at the cost of an extension of \(u\) outside of \(\Omega\), usually a periodic or nearest-neighbour recopy, or a zero-extension of the image \(u\) is considered, which is of course nonrealistic in most situations.

In order to use the Scilab implementation inverse_solver_tvreg, we need to compute the adjoint \(A^*\) and an upper bound of \(|||A|||\).

We can now implement those operators in Scilab, actually thanks to the built-in conv2 function, the code is straightforward. The (motion) blur kernel used for this experiment can be downloaded here in sod format (Scilab variable storing format).

load('motion_kernel.sod') ; // load the variable "ker" (a motion blur kernel) into Scilab memory...
// ...or define your own kernel, for instance: "ker = ones(5,5)/25;"
imview(ker.*.ones(10,10)); // visualize kernel (zoom x10)

//* set operator A = convolution with kernel "ker" *//
function v = A(u)
  v = conv2(u,ker,"valid");
endfunction

//* define the adjoint of A *//
ker_sym = ker($:-1:1,$:-1:1);
function u = adj_A(v)
  u = conv2(v,ker_sym,"full");
endfunction

//* set an upper bound for |||A||| *//
L_A = sum(abs(ker));

motion kernel (size 27 \(\times\) 27)

Let us now compute a blurry and noisy observation of a reference image. The reference image fuji.pgm used in this section is downloadable in PGM RAW format here (licence: CC0 public domain, source).

ref = readpgm('fuji.pgm'); // change the path according to the location of the image on your disk.
imview(ref);
sigma = 2;
u0 = A(ref);
u0 = u0 + grand(u0,'nor',0,sigma);
imview(u0);

reference image (domain \(\Omega\))	\(u_0\) blurry and noisy observation (domain \(\omega\))
\(~\)

Before running the algorithm, we suggest to use an initialization (this is optional but drastically reduce the number of iteration necessary to reach convergence). The input (\(u_0\)) of the algorithm has domain \(\omega\), the output (\(u\)) has domain \(\Omega\). We propose to extend \(u_0\) to the domain \(\Omega\) with a Neumann boundary condition (that is pixel recopy at the frontier of \(\omega\)) and use this extended image as initializer.

// compute init = u0 extended to domain Omega with pixel recopy
// boundary condition (drastically improve convergence speed)
[ny0,nx0] = size(u0);
dx = floor((size(ker,2)-1)/2);
dy = floor((size(ker,1)-1)/2);
init = u0([ones(1,dy),1:ny0,ny0*ones(1,dy)],[ones(1,dx),1:nx0,nx0*ones(1,dx)]);

We are now ready to use inverse_solver_tvreg to solve problem \eqref{eq:primal-inverse} with the Chambolle and Pock algorithm.

lambda = 0.2;
niter = 200;
[u,E] = inverse_solver_tvreg(u0,lambda,niter);
imview(u); // display output (restored)
figure(); plot(1:niter,E); // plot energy evolution
xlabel("iteration"); ylabel("energy");
title("Energy decrease");

reference	blurry and noisy (\(\sigma = 2\))
\(~\)
restored (\(200\) iterations, \(\lambda=0.2\))	energy decrease
\(~\)

To illustrate how we benefit from a good initialization, let us run the algorithm taking the null image as initializer for the primal variable \(u\),

init = zeros(ref);
lambda = 0.2;
v200  = inverse_solver_tvreg(u0,lambda,200);  // niter =  200
v2000 = inverse_solver_tvreg(u0,lambda,2000); // niter = 2000
v5000 = inverse_solver_tvreg(u0,lambda,5000); // niter = 5000

null initialization, 200 iterations	null initialization, 2000 iterations
\(~\)
null initialization, 5000 iterations	null initialization, 10000 iterations
\(~\)

The inverse problem \eqref{eq:primal-inverse} can also be used to perform image zooming (see [6,3]). For such application, the operator \(A\) is oftenly assumed to be a blurring kernel followed by a subsampling procedure.

We consider two discrete domains, a (small) domain \(\omega = \left[ 0, M \right) \times \left[ 0, N \right) \cap \ZZ^2\), and a (bigger) domain \(\Omega = \left[ 0, z M \right) \times \left[ 0, z N \right) \cap \ZZ^2\), where \(z\) denotes a given positive integer (zoom factor). We set \(\omega_A = \left[0,z\right) \times \left[0, z\right) \cap \ZZ^2\) and define \(A: \RR^\Omega \mapsto \RR^\omega\) by

Remark that \(Au\) is nothing more than a zero order unzoom with factor $z of the image \(u\), it consists in partitioning \(\Omega\) into square cells \(\omega_A\) and to average \(u\) over each one of the cells. The operator \(A\) can also be seen as a discrete approximation of a captor integration procedure, wich models the photon count averaging process that is done over the (in reality continuous) square domain \(\left[ 0, z \right) \times \left[ 0, z \right)\) of each captor covering the focal plane of the image acquisition system (for instance for CCD or CMOS cameras).

Let us now define those operators into Scilab memory.

//* set A operator = unzoom with factor z *//
z = 4; // set unzoom factor
function v = A(u)
  [ny,nx] = size(u);
  nx0 = nx/z; // warning: must be an integer
  ny0 = ny/z; // warning: must be an integer
  idx = 1+z*(0:(nx0-1));
  idy = 1+z*(0:(ny0-1));
  v = zeros(ny0,nx0);
  for dx = 0:(z-1)
    for dy = 0:(z-1)
      v = v + u(dy+idy,dx+idx);
    end
  end
  v = v/z^2;
endfunction

//* define the adjoint of A *//
k = ones(z,z)/z^2;
function u = adj_A(v)
  u = v .*. k;
endfunction

//* set an upper bound for |||A||| *//
L_A = 1/z;

Let us now unzoom with factor \(4\) a reference image, we will then rezoom it by solving problem \eqref{eq:primal-inverse}. The reference image eye.pgm used in this section is downloadable in PGM RAW format here (licence: CC0 public domain, source).

ref = readpgm('eye.pgm'); // change the path according to the location of the image on your disk.
imview(ref);
sigma = 2;
u0 = A(ref);
u0 = u0 + grand(u0,'nor',0,sigma);
imview(u0);

reference	unzoomed (\(\times 4\)) and noisy (\(\sigma = 2\))

Now let us use inverse_solver_tvreg to solve problem \eqref{eq:primal-inverse}, we will use as initializer a nearest neighbor zoom of \(u_0\).

niter = 2000;
lambda = 0.2;
init = u0.*.ones(z,z); // set init = nearest neighbor zoom of factor z of u0.
[u,E] = inverse_solver_tvreg(u0,lambda,niter);
imview(u); // display output (restored)
figure(); plot(1:niter,E); // plot energy evolution
xlabel("iteration"); ylabel("energy");
title("Energy decrease");

reference / \(u_0\)	details of \(u_0\) (nearest-neighbor zoom of \(u_0\))
\(~\)
energy decrease	restored (\(\lambda = 0.2\))
\(~\)

Thanks to Scilab built-in functions, we can easily rezoom \(u_0\) using the bicubic spline interpolation.

//* zoom u0 using bicubic spline interpolation *//
[ny0,nx0] = size(u0);
x = 0:(nx0-1);
y = 0:(ny0-1);
C = splin2d(y,x,u0,"natural");
xx = linspace(0,nx0-1,z*nx0);
yy = linspace(0,ny0-1,z*ny0);
[XX,YY] = ndgrid(xx,yy);
z_bicubic = (interp2d(YY,XX, y, x, C, "natural"))';
imview(z_bicubic);

bicubic spline zoom (\(\times 4\)) of \(u_0\)	\(\TV\) restored (\(\lambda = 0.2\))
\(~\)

The total variation based zooming yields much sharper image but clearly suffers to much from staircasing effect (that is the creation of flat region with artificial boundaries in the restored image) and from a strong block effect. We will see that the use of the Huber total variation allows to get rid of this undesirable artifact, but will not reduce the appearance of blocks. Those blocks are favored by the fidelity term in \eqref{eq:primal-inverse} (due to the particular structure of the operator \(A\) that we try to inverse, indeed \(Au\) is obtained by averaging \(u\) over non-overlapping square cells of size \(z \times z\)) and are not penalized by the \(\TV\) regularizer which is not very sensitive to sharp horizontal and vertical edges.

The spectrum extrapolation method has been first introduced by Guichard and Malgouyres in [5]. Let us assume that the spectrum (understand here its discrete Fourier transform) of the image \(u\) is known on the (small) domain \(\omega\), we want to extrapolate it to the bigger domain \(\Omega\).

We define \(A = \tfrac{|\omega|}{|\Omega|} \cdot \mathcal{F}^{-1}_\omega \circ R_{\omega} \circ \mathcal{F}_\Omega\) in order to have

$$A u = u_0 \quad \Leftrightarrow \quad \tfrac{|\omega|}{|\Omega|} \cdot R_\omega \circ \mathcal{F}_\Omega (u) = \mathcal{F}_\omega(u_0) \;,$$

where \(R_{\omega}\) denotes the restriction operator to the set \(\omega\) (that is \(R_{\omega} v = v_{|\omega}\), for any \(v \in \RR^\Omega\)), \(\mathcal{F}_\Omega\) and \(\mathcal{F}_\omega\) (resp. \(\mathcal{F}^{-1}_\omega\)) denote the direct (resp. inverse) discrete Fourier transforms of two-dimensional signals in \(\RR^\Omega\) and \(\RR^\omega\). The multiplicative factor \(\tfrac{|\omega|}{|\Omega|}\) is used to ensure that \(u\) and \(u_0\) share the the same mean.

Notice that in [5], Guichard and Malgouyres are interested in computing the image \(u\) with smallest total variation among those satisfying the constraint \(Au = u_0\), that is

The resolution of such problem is adressed in this page dedicated to the resolution of constrained minimization of the total variation). Instead of \eqref{eq:tvextrap} we consider here the relaxed version given by problem \eqref{eq:primal-inverse}, which amounts to compute the image that satisfy a good compromise between two requirements: "having a small total variation" and "verifying Au ≈ u₀".

Before going further, we need to discuss an important point that is usually ignored in the litterature. In general, when at least one of the dimensions of \(\omega\) is even, the use of \(\mathcal{F}_\omega\) and \(\mathcal{F}_\omega^{-1}\) in our image processing algorithms generates nonreal values because of the dissymmetry of \(\omega\) which induces in general non-Hermitian symmetry for the spectrum of manipulated images. In particular here, if \(\omega\) has not two odd dimensions, we obtain in general images with nonreal values when computing \(Au\) or \(A^* v\) (for \(u\in\RR^\Omega\) and \(v \in \RR^\omega\)). See this page to better understand where this phenomenon comes from and how it can be correctly handled.

For shake of simplicity, we will here simply restrict the use of the application to images \(u_0\) having odd \(\times\) odd dimensions.

Let us now configure our inverse_solver_tvreg algorithm to perform spectrum extrapolation. The reference image hummingbird.pgm used in this section is downloadable in PGM RAW format here (licence: CC0 public domain, source).

//* load the reference image *//
ref = readpgm('hummingbird.pgm'); // change the path according to the location of the image on your disk.

//* precompute Omega and omega dimensions *//
z = 3; // set the square root of the ratio |Omega|/|omega|
[ny,nx] = size(ref);
nx0 = round(nx/z);
ny0 = round(ny/z);

//* ensures omega has odd x odd dimensions *//
if(modulo(nx0,2)==0) then nx0 = nx0-1; end
if(modulo(ny0,2)==0) then ny0 = ny0-1; end
z = sqrt((nx*ny)/(nx0*ny0)); // update z = sqrt(|Omega|/|omega|)


//* define omega_x, omega_y = indexes of frequencies omega into Omega. *//
//* this allows to restrict the Fourier domain Omega to omega.         *//
omega_x = 1+modulo(nx + ifftshift((0:nx0-1)-floor(nx0/2)), nx);
omega_y = 1+modulo(ny + ifftshift((0:ny0-1)-floor(ny0/2)), ny);

//* define operator A *//
function v = A(u)
  dft_u = fft(u);
  dft_v = dft_u(omega_y,omega_x);
  v = 1/z^2 * real(ifft(dft_v));
endfunction

//* define the adjoint of A *//
function u = adj_A(v)
  dft_u = zeros(ny,nx);
  dft_u(omega_y,omega_x) = fft(v);
  u = real(ifft(dft_u));
endfunction

//* set an upper bound for |||A||| *//
L_A = 1;

We compute now a noisy observation \(u_0\) of the reference image.

sigma = 2; // set the level of noise
u0 = A(ref) + grand(ny0,nx0,'nor',0,sigma);
imview(ref); // display reference image
imview(u0); // display observed image

reference image	observed image \(u_0\) (\(\sigma=2\))

Let us run our inverse_solver_tvreg algorithm.

lambda = 0.1;
niter = 300;
clear init; // we will use default initializer (null image)
// but you can try: "init = z^2*adj_A(u0);" to set init = zero-padding zoom of u0.
[u,E] = inverse_solver_tvreg(u0,lambda,niter);
imview(u); // display output (restored)
imview(fftshift(log(1+abs(fft(u))))); // display log(1+Fourier-modulus) of the restored image
imview(fftshift(log(1+abs(fft(ref))))); // display log(1+Fourier-modulus) of the reference image
figure(); plot(1:niter,E); // plot energy evolution
xlabel("iteration"); ylabel("energy");
title("Energy decrease");

reference image	restored image (\(\lambda=0.1\))
Log(1+Fourier-modulus)	Log(1+Fourier-modulus)

energy decrease

Let us resample \(u_0\) over \(\Omega\) using the bicubic spline interpolation.

//* resample u0 over domain Omega using a bicubic spline interpolation *//
x = 0:(nx0-1);
y = 0:(ny0-1);
C = splin2d(y,x,u0,"natural");
xx = linspace(0,nx0-1,nx);
yy = linspace(0,ny0-1,ny);
[XX,YY] = ndgrid(xx,yy);
z_bicubic = (interp2d(YY,XX, y, x, C, "natural"))';
imview(z_bicubic);
imview(fftshift(log(1+abs(fft(z_bicubic))))); // display log(1+Fourier-modulus)

bicubic spline resampling of \(u_0\) over \(\Omega\)	TV restored image (\(\lambda=0.1\))
Log(1+Fourier modulus)	Log(1+Fourier-modulus)

The image obtained by mean of bicubic interpolation of \(u_0\) is too smooth and presents some ringing artifacts (some edges are replicated into vanishing parallel edges spaced by 2 pixels). Besides, the looking at the modulus of its Fourier transform, we see that the spectrum was completed in a very poor and unrealistic way. The total variation based spectrum extrapolation method yields a better quality image, with a much richer and natural spectrum.