% Deblurs the image J = E*i + N, where i = I(:),
% using the Optimization Toolbox function lsqlin.
% J, E, N must be defined in the workspace.

% h = ones(3, 3)/9;
% B = imfilter(J, h);  % Apply an averaging filter to approximate B = J - N.
B = J-N;
b = B(:);

[m,n] = size(J);
mn = m*n;

options = optimset('LargeScale', 'on', 'Display', 'off');
i = lsqlin(E, b, [], [], [], [], zeros(mn, 1), 255*ones(mn, 1), [], options);
I = reshape(i, m, n);

subplot(2,1,1)
imagesc(J);
colormap(gray)
axis image
title('Blurred Image J')

subplot(2,1,2)
imagesc(I);
colormap(gray)
axis image
title('Deblurred Image I')

clear h B b m n mn options i 
