This example shows how to recover a blurred image by solving a large-scale bound-constrained linear least-squares optimization 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; imagesc(P); colormap(gray); axis off image; title([int2str(m) ' x ' int2str(n) ' (' int2str(mn) ') pixels'])
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)')
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(:); imagesc(reshape(G,m,n)); axis off 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); imagesc(reshape(G(:),m,n)); axis off image; title('G, the blurred image'); subplot(1,2,2); imagesc(reshape(P(:),m,n)); axis off image; title('P, the original image');
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')
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<=P<=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''');
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); imagesc(reshape(x,m,n)) axis off image
Let's compare the blurred and deblurred pictures.
subplot(1,2,1); imagesc(reshape(G,m,n)); axis image; axis off; title('Blurred'); subplot(1,2,2); imagesc(reshape(x,m,n)); axis image; axis off; title('De-Blurred');