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');