d = double(imread('jupiter1.tif'));
figure(1); imagesc(d); colormap(gray); axis equal
h = d(70:101,220:251);
h = h./sum(sum(h));
small = find(h<0.05*max(max(h)));
h(small) = 0;
D = fft2(d,256,256);
H = fft2(h,256,256);
%
% Draw an L-curve
%
llist = logspace(-4,0,51);
l = length(llist);
C = zeros(l,1);
S = zeros(l,1);
h = waitbar(0,'Scanning lambda');
for k = 1:l
    lambda = llist(k);
    F = D.*conj(H)./(abs(H).^2+lambda.^2);
    S(k) = norm(F(:));
    C(k) = norm(D(:)-H(:).*F(:));
    waitbar(k/l);
end
close(h);
figure(2); loglog(C,S,'x');
for k = 1:2:l
    text(C(k),S(k),num2str(llist(k)))
end    
lambda = 0.1;
F = D.*conj(H)./(abs(H).^2+lambda.^2);
f = real(ifft2(F));
figure(3); imagesc(f); colormap(gray); axis equal