Contents

%  image_deblurring.m
%
%  Set up a discretization of a convolution integral operator K with a
%  Gaussian kernel. Generate a true solution and convolve it with the
%  kernel K. Then add random noise to the resulting data.
%
%  Matthias Heinkenschloss. September 14, 2011
%  Modified by Yin Zhang.   September 12, 2013

set(0, 'defaultaxesfontsize',14,'defaultaxeslinewidth', 1.25,...
    'defaultlinelinewidth',1.25,'defaultpatchlinewidth',1.25,...
    'defaulttextfontsize',14);

clear all

Set up parameters.

n       =  80; % input(' No. of grid points = ');
sig     = .05; % input(' Kernel width sigma = ');
err_lev = 2;   % input(' Percent error in data = ');

Set up grid.

h = 1/n;            % interval length
x = (h/2:h:1-h/2)'; % interval midpoints

Compute matrix K corresponding to convolution with Gaussian kernel

kernel = (1/sqrt(pi)/sig) * exp(-(x-h/2).^2/sig^2);
K = toeplitz(kernel)*h;

Set up true solution f_true and data d = K*f_true + error.

f_true = 0.75*(0.1<x&x<0.25) + 0.25*(0.25<x&x<0.32) + 0.9*(0.5<x&x<0.8);
Kf     = K*f_true;
rng('default')
eta    = err_lev/100 * norm(Kf) * randn(n,1)/sqrt(n); % error
d      = Kf + eta;                                    % measurement

Display the data.

figure(1)
mesh(K)
title('Mesh Plot Representation of K')

figure(2)
plot(x,f_true,'-', x,d,'o',x,Kf,'--')
xlabel('x axis')
title('True Solution and Discrete Noisy Data')
legend('true image', 'blurred image', 'blurred image + noise', 'Location','NorthWest')
%print -depsc deblurr_data.eps

Least Squares approximation

fprintf( '\nL2 approximation \n' )
f_lls = pinv(K)*d;

figure(3)
plot(x,f_true,'b-',x,f_lls,'r--')
axis([0 1 -5 5])
xlabel('x axis')
title('True Solution and Least Squares Estimate')
legend('true image', 'LLS estimate', 'Location','NorthWest')
%print -depsc2 deblurr_lls.eps

L2 approximation 

L1 approximation

fprintf( '\nL1 approximation \n' )
[m,n] = size(K);
K = sparse(K);
Im = speye(m);

A   = [  K  -Im
        -K  -Im];
c   = [zeros(n,1); ones(m,1)];
b   = [d; -d];

[f_l1,fval,exitflag] = linprog(c,A,b);

fprintf( 'linprog terminated with exitflag = %d \n', exitflag )
fprintf( '        and objective function value   %12.6e \n', fval )

f_l1 = f_l1(1:n);

figure(4)
plot(x,f_true,'b-',x,f_l1,'r--')
axis([0 1 -5 5])
xlabel('x axis')
title('True Solution and L1 Estimate')
legend('true image', 'L1 estimate', 'Location','NorthWest')
%print -depsc2 deblurr_l1.eps

L1 approximation 
Exiting: Maximum number of iterations exceeded; increase options.MaxIter.
linprog terminated with exitflag = 0 
        and objective function value   6.795403e-01 

L1 approximation with nonnegativity constraint

fprintf( '\nL1 approximation with nonnegativity constraint\n' )
[f_l1,fval,exitflag] = linprog(c,A,b,[],[],zeros(n+m,1),[]);

fprintf( 'linprog terminated with exitflag = %d \n', exitflag )
fprintf( '        and objective function value   %12.6e \n', fval )

f_l1 = f_l1(1:n);

figure(5)
plot(x,f_true,'b-',x,f_l1,'r--')
xlabel('x axis')
title('True Solution and L1 Estimate + nonnegativity')
legend('true image', 'L1 estimate', 'Location','NorthWest')
%print -depsc2 deblurr_l10.eps

L1 approximation with nonnegativity constraint
Optimization terminated.
linprog terminated with exitflag = 1 
        and objective function value   5.675778e-01 

L1 approximation with TV regularization

fprintf( '\nL1 approximation with TV regularization\n' )
alpha = 2;  % regularization parameter
e = ones(n-1,1);
Diff = spdiags([e -e],0:1,n-1,n);
Zmn1 = sparse(m,n-1);
In1  = speye(n-1);

A   = [  K      -Im     Zmn1
        -K      -Im     Zmn1
         Diff  -Zmn1'   -In1
        -Diff  -Zmn1'   -In1 ];
c   = [ zeros(n,1); ones(m,1); alpha*ones(n-1,1)];
b   = [d; -d; zeros(n-1,1); zeros(n-1,1)];

[f_l1,fval,exitflag] = linprog(c,A,b);

fprintf( 'linprog terminated with exitflag = %d \n', exitflag )
fprintf( '        and objective function value   %12.6e \n', fval )

f_l1 = f_l1(1:n);

figure(6)
plot(x,f_true,'b-',x,f_l1,'r--')
xlabel('x axis')
title('True Solution and L1 Estimate with TV regularization')
legend('true image', 'L1 estimate', 'Location','NorthWest')
axis([0,1,-0.1,1.1])
%print -depsc2 deblurr_l1_TV2.eps

L1 approximation with TV regularization
Exiting: One or more of the residuals, duality gap, or total relative error
 has grown 100000 times greater than its minimum value so far:
         the dual appears to be infeasible (and the primal unbounded).      
         (The primal residual < TolFun=1.00e-08.)
linprog terminated with exitflag = -3 
        and objective function value   2.991629e+00 

L1 approximation with nonnegativity constraint and TV regularization

fprintf( '\nL1 approximation with TV regularization and nonnegativity\n' )
lB = zeros(2*n+m-1,1);
[f_l1,fval,exitflag] = linprog(c,A,b,[],[],lB);

fprintf( 'linprog terminated with exitflag = %d \n', exitflag )
fprintf( '        and objective function value   %12.6e \n', fval )

f_l1 = f_l1(1:n);

figure(7)
plot(x,f_true,'b-',x,f_l1,'r--')
xlabel('x axis')
title('True Solution and L1 Estimate with TV regularization + nonnegativity')
legend('true image', 'L1 estimate', 'Location','NorthWest')
axis([0,1,-0.1,1.1])
%print -depsc2 deblurr_l1_TV20.eps

L1 approximation with TV regularization and nonnegativity
Optimization terminated.
linprog terminated with exitflag = 1 
        and objective function value   7.302192e+00 

Published with MATLAB® R2013b