%%
clear
clc
format compact

%%
% We repeat the earlier process to blur the original image X to get Z. 
load('mandrill', 'X')
[m, n] = size(X)
v = [1/4 1/2 1/4];
B = spdiags(repmat(v,m,1), -1:1, m, m);
B(1, m) = v(3); % periodic BC
B(m, 1) = v(1);
C = spdiags(repmat(v,n,1), -1:1, n, n);
C(1, n) = v(3); % periodic BC
C(n, 1) = v(1);
blur = @(X, k) B^k * X * C^k;
nbl = 16;
Z = blur(X, nbl);

%%
% Now we imagine that X is unknown and that the blurred Z is given.
% We want to invert the blur transformation using the transformation itself.
% But we have to translate between vectors and images each time. 
vec = @(X) reshape(X, m*n, 1);
unvec = @(x) reshape(x, m, n);
T = @(x) vec(blur(unvec(x), nbl));

%%
% Now we apply gmres to the composite blurring transformation T.
y = gmres(T, vec(Z), 50, 1e-5);
Y = unvec(y);

subplot(131)
image(X)
colormap(gray(256))
axis equal
axis tight
title('Original')

subplot(132)
image(Z)
colormap(gray(256))
axis equal
axis tight
title('Blurred')

subplot(133)
image(Y)
colormap(gray(256))
axis equal
axis tight
title('Deblurred')