% inpaint1_template.m
% template for 2D in-painting

%% 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+[0:4],ny/4+[0:1]) = 1; xtrue(nx/2:3*nx/4,ny/2:3*ny/4) = 1;

%% ideal noiseless data (just the true image in a denoising problem)
rng(0)
samp = rand(nx,ny) > 0.5; % sample only this fraction of the pixels
A = Gdiag(samp);
ytrue = A * xtrue;

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

im plc 2 2
clim = [-0.1 1.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 = Cfd2(nx,ny); % 2D finite differences

%% run steepest descent

reg = 1.1; % todo: experiment with this to find improved SNR
step = 1 / (1 + reg * 8); % todo: is this 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)
	% todo: modify to use edge-preserving regularization
	grad = A' * (A * x - y) + reg * (C' * (C * x));
%	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(:))))
