% denoise_2d_template.m
% template for 2D denoising

%% noiseless original image
nx = 32; ny = 20; % image size (must be a multiple of 4 for next line)
xtrue = zeros(nx,ny); xtrue(nx/4,ny/4) = 2; xtrue(nx/2:3*nx/4,ny/2:3*ny/4) = 1;

%% ideal noiseless data (just the true image in a denoising problem)
A = 1; % trick: for denoising can use "1" instead of identity matrix in code
	% an alg
ytrue = A * xtrue;

%% add noise
rng(0)
snr = 12; % specify desired SNR in dB
sig = 10^(-snr/20) * norm(ytrue(:)) / sqrt(numel(ytrue));
y = ytrue + sig * randn(size(ytrue));

im plc 2 2
clim = [-0.1 2.1]; % display all images on same gray scale
im(1, xtrue, clim)
im(3, ytrue, clim)
im(4, y, clim)
xlabelf('SNR = %4.1f dB', 20*log10(norm(xtrue(:)) / norm(y(:) - xtrue(:))))
drawnow

C = 1; % todo: replace this with an appropriate 2D fatrix

%% run steepest descent

reg = 1.1; % todo: experiment with this to find improved SNR
step = 0.1; % todo: replace with appropriate step size!
niter = 200; % max # of iter
tol = 1e-5; % stop tolerance for change in x

x = A' * y; % initial guess is just the noisy image itself
for iter=1:niter
	ticker(mfilename, iter, niter)
	grad = A' * (A * x - y) + reg * x; % todo: modify to consider C
%	pr norm(grad(:))
	x = x - step * grad;
	if (norm(step*grad(:)) < norm(x(:)) * tol)
		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(:))))
