% deblur1_rms.m
% Image deblurring for several values of regularization parameter.
% In this simulation we know the true image so we can find the minimum RMSE
% value of the regularization parameter.

%% 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 (blurred true image)
A = Gblur(true(nx,ny), 'psf', ones(3)/9);
ytrue = A * xtrue;
nd = numel(ytrue);

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

% function for computing normalized RMS error in %:
nrmse_fun = @(x) 100 * norm(x(:) - xtrue(:)) / norm(xtrue(:));
snr_fun = @(x) 20*log10(norm(xtrue(:)) / norm(x(:) - xtrue(:)));

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

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


delta = 0.1; % small compared to max(xtrue) - min(xtrue)
dpot = @(t) t ./ sqrt(1 + (t / delta).^2); % derivative of hyperbola potential
niter = 200; % max # of iter
tol = 1e-5; % stop tolerance for change in x

%% try several values of regularization parameter and examine RMS errors

reglist = 2 .^ [-6:-1]; % working over a log scale is often helpful
nr = numel(reglist);
nrmse = zeros(nr,1); % store RMS error values in this array
rss = zeros(nr,1); % store RSS values in this array
nrmse_best = nrmse_fun(y);
x_best = y;
for ir=1:nr
	reg = reglist(ir);
	step = 1 / (1 + reg * 8); % todo: is this appropriate step size?

	xinit = y; % run GD with edge-preserving regularization
	x = gradient_descent1(A, y, xinit, reg, C, dpot, step, tol, niter);

	rss(ir) = norm(col(A*x) - y(:))^2 / sig^2; % RSS

	nrmse(ir) = nrmse_fun(x);
	if nrmse(ir) < nrmse_best % archive the lowest NRMSE case
		nrmse_best = nrmse(ir);
		x_best = x;
	end

	im(4, x, clim)
	xlabelf('NRMSE = %4.1f %%', nrmse(ir))

	im subplot 5
	plot(log2(reglist), rss, '-o')
	titlef 'RSS'

	im subplot 6
	plot(log2(reglist), nrmse, '-o')
	titlef 'NRMSE (%%)'
	xlabelf 'log_2(\beta)'
	drawnow
end

	im(4, x_best, clim)
	xlabelf('NRMSE = %4.1f %%', nrmse_best)
