www.gusucode.com > optim 案例源码 matlab代码程序 > optim/optdeblur.m

    %% Large-Scale Constrained Linear Least-Squares
% This example shows how to recover a blurred image by solving a large-scale 
% bound-constrained linear least-squares optimization problem.

% Copyright 1990-2016 The MathWorks, Inc.

%% The Problem
% We will add motion blur to a photo of Mary Ann and Matthew sitting in Joe's
% car, then try to restore the original. Our starting image is this black and
% white image, contained the m x n matrix P. Each element in the matrix
% represents a pixel's gray intensity between black and white (0 and 1).

load optdeblur P
[m,n] = size(P);
mn = m*n;
imshow(P);
colormap(gray);
axis off image;
title([int2str(m) ' x ' int2str(n) ' (' int2str(mn) ') pixels']) 

%% Add Motion
% We can simulate the effect of vertical motion blurring by averaging each pixel
% with the 5 pixels above and below. We construct a sparse matrix D, that will 
% do this with a single matrix multiply.
%
% Create D.

blur = 5;  mindex = 1:mn;  nindex = 1:mn;
for i = 1:blur
  mindex=[mindex i+1:mn 1:mn-i];
  nindex=[nindex 1:mn-i i+1:mn];
end
D = sparse(mindex,nindex,1/(2*blur+1));

%%
% Draw a picture of D.

cla
axis off ij
xs = 25;
ys = 15;
xlim([0,xs+1]);
ylim([0,ys+1]);
[ix,iy] = meshgrid(1:(xs-1),1:(ys-1));
l = abs(ix-iy)<=5;
text(ix(l),iy(l),'x')
text(ix(~l),iy(~l),'0')
text(xs*ones(ys,1),1:ys,'...');
text(1:xs,ys*ones(xs,1),'...');
title('Blurring Operator D ( x = 1/11)')

%% Apply Blurring Operator
% We multiply the image by this operator to create the blurred image. P is the
% original image, D is the operator, and G is the blurred image.

G = D*P(:);
imshow(reshape(G,m,n));

%% The Blurred Image
% Now, let's pretend Joe took this blurred picture G from a moving elevator.
% Assume we know how fast the elevator is moving, so we know the blurring
% operator D. How well can we remove the blur and recover the original image P?
% 
% The simplest approach is to solve the least squares problem:
% 
% $$ \min( \| D P(:) - G\|^2 )$$

subplot(1,2,1);
imshow(reshape(G(:),m,n)); 
title('G, the blurred image');
subplot(1,2,2);
imshow(reshape(P(:),m,n));
title('P, the original image');

%% Regularization Term
% In practice the results obtained with this simple approach tend to be noisy.
% To compensate for this, a regularization term is added:
%
% $$ 0.0004*\| L P \|^2 $$
%
% L is the discrete Laplacian, which relates each pixel to those surrounding it. 
% Since we know we are looking for a gray intensity, we also impose the constraint 
% that the elements of P must fall between 0 and 1.
%
% Create L.

L = sparse( [1:mn,2:mn,1:mn-1],  [1:mn,1:mn-1,2:mn], ...
   [4*ones(1,mn) -1*ones(1,2*(mn-1))]  );

%%
% Draw a picture of L.

subplot(1,1,1) ;
delete(gca);
axis ij
axis off;
xs=11;
ys=11;
xlim([0,xs+1]);
ylim([0,ys+1]);
[ix,iy]=meshgrid(1:(xs-1),1:(ys-1));
four=(ix==iy);
one=(abs(ix-iy)==1);
text(ix(one),iy(one),'-1')
text(ix(four),iy(four),'+4')
text(ix(~four & ~one),iy(~four & ~one),' 0')
text(xs*ones(ys,1),1:ys,'...');
text(1:xs,ys*ones(xs,1),'...');
title('L, The Discrete Laplacian')

%% Set Up Problem to Deblur Image
% To obtain the deblurred picture we want to solve for P:
%
% $$ \min( \| D P(:) - G(:) \|^2 + 0.0004 \| L P(:) \|^2 ) $$
%
% We can simplify this expression by defining A and b:
%
%    A = [D ; 0.02*L]; b = [ G(:) ; zeros(mn,1) ];
%
% which changes the last equation to:
%
% $$ \min( \| A P(:) - b \|^2 ) $$
%
% subject to $0\le P\le 1$. Both matrices D and L relate each pixel to a few of its
% neighbors. This makes A structured and sparse.

A = [D ; 0.02*L]; b = [ G(:) ; zeros(mn,1) ];
spy(A');
axis equal tight
title('Structure of Matrix A''');

%% Solve Optimization Problem
% Because A is sparse, we can use a large-scale algorithm to solve this linear
% least squares optimization problem. We call LSQLIN with A, b, lower bounds,
% upper bounds, and options.

options = optimoptions('lsqlin','Algorithm','trust-region-reflective',...
    'Display','off');
x = lsqlin(A, b, [], [], [], [], zeros(mn,1), ones(mn,1), [], options);
imshow(reshape(x,m,n))

%% Compare Images
% Let's compare the blurred and deblurred pictures.

subplot(1,2,1);
imshow(reshape(G,m,n)); 
title('Blurred');

subplot(1,2,2);
imshow(reshape(x,m,n)); 
title('De-Blurred');