% deblur2_race.m
% "racing" optimization algorithms for deblurring
% template for student solutions

%% read predefined data
if ~isvar('y')
	load deblur2_data
	[nx ny] = size(y);

	%% system matrix for blur, with non-periodic edge conditions
	A = Gblur(true([nx ny]), 'psf', psf);

	C = Cfd2(nx,ny); % 2D finite differences
end


im plc 2 2
clim = [-10 130]; % display all images on same gray scale
im(1, xtrue, clim)
im(3, pad_into_center(psf, 8), 'psf')
im(4, y, clim)
drawnow


%% image deblurring parameters

reg = 2^-1; % keep this fixed for fair comparisons
step = 1 / (1 + reg * 8); % because rho(A'A) = 1 and rho(C'C) = 2*4
niter = 300; % max # of iter

grad0 = A' * (A * (0*xtrue) - y); % gradient at "0" image to normalize tolerance
tol_grad = 8e-4 * norm(grad0(:)); % stop tolerance for gradient

delta = 2; % small compared to max(xtrue) - min(xtrue).  keep fixed
dpot = @(t) t ./ sqrt(1 + (t / delta).^2); % derivative of hyperbola potential

x = y; % initial guess is just the noisy image itself


%% do not change anything above this line! %%%%%%%%%%%%%%%%%%%%%%%


%% todo: replace GD with a faster optimization method of your own 
%% after you find a faster version, report your # of iterations to google doc
%% then try to find a still faster version (until time runs out).
%% report your # of iterations for each method to google doc

%% run gradient descent (requires 206 iterations)
for iter=1:niter
	ticker(mfilename, iter, niter)
	grad = A' * (A * x - y) + reg * (C' * dpot(C * x)); % gradient
	x = x - step * grad;
	if (norm(grad(:)) < tol_grad)
		printm('stopped at %d iterations', iter)
		break
	end
end
im(2, x, clim)
xlabelf('SNR = %4.1f dB', 20*log10(norm(xtrue(:)) / norm(x(:) - xtrue(:))))
