Remove the blurr using Tikhonov or TGSVD methods

This is my Matlab code, where its applied blurr to a image of a Mandril. The main ideia is using one of Tikhonov or TGSVD methods to remove the blurr.

clear all

close all

clc

load mandrill

clear caption map

N = 480;

X_exact = double(X(:,10+(1:N)));

x_exact = X_exact(:);

% Creation of point-spread-function

[I, J] = meshgrid(0.01 * (-15:15));

P = (1 + 1 / (0.25 * 0.35 - 0.1^2) * ((I / 0.25).^2 - 2 * 0.1 * I .* J + (J / 0.35).^2)).^(-1);

P = P / sum(P(:));

% Creation of blurred image

b_exact = Amult(P, x_exact);

e = randn(size(b_exact));

e = e / norm(e);

e = 0.01 * norm(b_exact) * e;

b = b_exact + e;

B = reshape(b, N, N);

% Display original and blurred images

figure(1)

subplot(1, 2, 1)

imshow(X, []);

title("Original Image")

subplot(1, 2, 2)

imshow(B, []);

title("Blurred Image")

I want to use one of those twos methods to remove the blurr in the image. Any ideia how?

I'm really lost, I'm sorry that I can't give more details.

This is the Tikhonov function

function [x_lambda,rho,eta] = tikhonov(U,s,V,b,lambda,x_0)

%TIKHONOV Tikhonov regularization.

%

% [x_lambda,rho,eta] = tikhonov(U,s,V,b,lambda,x_0)

% [x_lambda,rho,eta] = tikhonov(U,sm,X,b,lambda,x_0) , sm = [sigma,mu]

%

% Computes the Tikhonov regularized solution x_lambda. If the SVD

% is used, i.e. if U, s, and V are specified, then standard-form

% regularization is applied:

% min { || A x - b ||^2 + lambda^2 || x - x_0 ||^2 } .

% If, on the other hand, the GSVD is used, i.e. if U, sm, and X are

% specified, then general-form regularization is applied:

% min { || A x - b ||^2 + lambda^2 || L (x - x_0) ||^2 } .

%

% If x_0 is not specified, then x_0 = 0 is used

%

% Note that x_0 cannot be used if A is underdetermined and L ~= I.

%

% If lambda is a vector, then x_lambda is a matrix such that

% x_lambda = [ x_lambda(1), x_lambda(2), ... ] .

%

% The solution norm (standard-form case) or seminorm (general-form

% case) and the residual norm are returned in eta and rho.

% Per Christian Hansen, IMM, April 14, 2003.

% Reference: A. N. Tikhonov & V. Y. Arsenin, "Solutions of

% Ill-Posed Problems", Wiley, 1977.

% Initialization.

if (min(lambda)<0)

error('Illegal regularization parameter lambda')

end

m = size(U,1);

n = size(V,1);

[p,ps] = size(s);

beta = U(:,1:p)'*b;

zeta = s(:,1).*beta;

ll = length(lambda); x_lambda = zeros(n,ll);

rho = zeros(ll,1); eta = zeros(ll,1);

% Treat each lambda separately.

if (ps==1)

% The standard-form case.

if (nargin==6), omega = V'*x_0; end

for i=1:ll

if (nargin==5)

x_lambda(:,i) = V(:,1:p)*(zeta./(s.^2 + lambda(i)^2));

rho(i) = lambda(i)^2*norm(beta./(s.^2 + lambda(i)^2));

else

x_lambda(:,i) = V(:,1:p)*...

((zeta + lambda(i)^2*omega)./(s.^2 + lambda(i)^2));

rho(i) = lambda(i)^2*norm((beta - s.*omega)./(s.^2 + lambda(i)^2));

end

eta(i) = norm(x_lambda(:,i));

end

if (nargout > 1 & size(U,1) > p)

rho = sqrt(rho.^2 + norm(b - U(:,1:n)*[beta;U(:,p+1:n)'*b])^2);

end

elseif (m>=n)

% The overdetermined or square general-form case.

gamma2 = (s(:,1)./s(:,2)).^2;

if (nargin==6), omega = V\x_0; omega = omega(1:p); end

if (p==n)

x0 = zeros(n,1);

else

x0 = V(:,p+1:n)*U(:,p+1:n)'*b;

end

for i=1:ll

if (nargin==5)

xi = zeta./(s(:,1).^2 + lambda(i)^2*s(:,2).^2);

x_lambda(:,i) = V(:,1:p)*xi + x0;

rho(i) = lambda(i)^2*norm(beta./(gamma2 + lambda(i)^2));

else

xi = (zeta + lambda(i)^2*(s(:,2).^2).*omega)./...

(s(:,1).^2 + lambda(i)^2*s(:,2).^2);

x_lambda(:,i) = V(:,1:p)*xi + x0;

rho(i) = lambda(i)^2*norm((beta - s(:,1).*omega)./...

(gamma2 + lambda(i)^2));

end

eta(i) = norm(s(:,2).*xi);

end

if (nargout > 1 & size(U,1) > p)

rho = sqrt(rho.^2 + norm(b - U(:,1:n)*[beta;U(:,p+1:n)'*b])^2);

end

else

% The underdetermined general-form case.

gamma2 = (s(:,1)./s(:,2)).^2;

if (nargin==6), error('x_0 not allowed'), end

if (p==m)

x0 = zeros(n,1);

else

x0 = V(:,p+1:m)*U(:,p+1:m)'*b;

end

for i=1:ll

xi = zeta./(s(:,1).^2 + lambda(i)^2*s(:,2).^2);

x_lambda(:,i) = V(:,1:p)*xi + x0;

rho(i) = lambda(i)^2*norm(beta./(gamma2 + lambda(i)^2));

eta(i) = norm(s(:,2).*xi);

end

end

And this is the TSVD function

function [x_k,rho,eta] = tsvd(U,s,V,b,k)

%TSVD Truncated SVD regularization.

%

% [x_k,rho,eta] = tsvd(U,s,V,b,k)

%

% Computes the truncated SVD solution

% x_k = V(:,1:k)*inv(diag(s(1:k)))*U(:,1:k)'*b .

% If k is a vector, then x_k is a matrix such that

% x_k = [ x_k(1), x_k(2), ... ] .

%

% The solution and residual norms are returned in eta and rho.

% Per Christian Hansen, IMM, 12/21/97.

% Initialization.

[n,p] = size(V); lk = length(k);

if (min(k)<0 | max(k)>p)

error('Illegal truncation parameter k')

end

x_k = zeros(n,lk);

eta = zeros(lk,1); rho = zeros(lk,1);

beta = U(:,1:p)'*b;

xi = beta./s;

% Treat each k separately.

for j=1:lk

i = k(j);

if (i>0)

x_k(:,j) = V(:,1:i)*xi(1:i);

eta(j) = norm(xi(1:i));

rho(j) = norm(beta(i+1:p));

end

end

if (nargout > 1 & size(U,1) > p)

rho = sqrt(rho.^2 + norm(b - U(:,1:p)*beta)^2);

end