gfp.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A MATLAB code to demonstrate the GFP-SOR method for 
% Total Variation denoising model (ROF model)
%
% Min { alpha*Tv(u) + (1/2)*int(u-z)^2 } (*)
%  u
%
% This is equaivalent to solve the following PDE:
%
% -alpha*div[Grad(u)/|Grad(u)|]+(u-z) = 0
%
% Here u is the image to be recovered 
%      z is the observed image corrupted by the Gaussian noise 
% and  
%      alpha > 0 is the regularization parameter.
% 
% Reference: 
%  (1) Rudin, Osher and Fatemi, 'Nonlinear total variation based 
%      noise removal algorithms', Phys. D 60, pp. 259-268, 1992.
%
%  (2) C.R. Vogel and M.E. Oman,'Iterative methods for total variation 
%      denoising, SIAM J. Sci. Comput. 17, pp. 227-238, 1996.
%
% LAST MODIFIED: 2018-January-04
%
% Programed by (for SEAMS School2018)
%
% Asistant Profesor Dr. Noppadol Chumchob
% Department of Mathematics,
% Silpakorn University, 
% Nakhon-Pathom, 73000, THAILAND.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [u] = gfp(u,z,alpha,beta,omega,GSiter) 

[n,m] = size(u);
h = 1;
D = CoefficientD(u,beta);
G = z;

for iter=1:GSiter

u(1,1) = (1-omega)*u(1,1)+...
    omega*(G(1,1)+(alpha/h^2)*(D(1,1)*(u(2,1)+u(1,2))))...
    /(1+(alpha/h^2)*(2*D(1,1)));

if n>2
for i=2:n-1
    
u(i,1) = (1-omega)*u(i,1)+...
    omega*(G(i,1)+(alpha/h^2)*(D(i,1)*(u(i+1,1)+u(i,2))+D(i-1,1)*u(i-1,1)))...
    /(1+(alpha/h^2)*(2*D(i,1)+D(i-1,1)));

end
end

u(n,1) = (1-omega)*u(n,1)+...
    omega*(G(n,1)+(alpha/h^2)*(D(n,1)*(u(n,2))+D(n-1,1)*u(n-1,1)))...
    /(1+(alpha/h^2)*(D(n,1)+D(n-1,1)));

if m>2
for j=2:m-1

u(1,j) = (1-omega)*u(1,j)+...
    omega*(G(1,j)+(alpha/h^2)*(D(1,j)*(u(2,j)+u(1,j+1))+D(1,j-1)*u(1,j-1)))...
    /(1+(alpha/h^2)*(2*D(1,j)+D(1,j-1)));

for i=2:n-1
    
u(i,j) = (1-omega)*u(i,j)+omega*(G(i,j)+(alpha/h^2)*...
    (D(i,j)*(u(i+1,j)+u(i,j+1))+D(i-1,j)*u(i-1,j)+D(i,j-1)*u(i,j-1)))/...
    (1+(alpha/h^2)*(2*D(i,j)+D(i-1,j)+D(i,j-1)));

end
u(n,j) = (1-omega)*u(n,j)+...
    omega*(G(n,j)+(alpha/h^2)*(D(n,j)*(u(n,j+1))+D(n-1,j)*u(n-1,j)+D(n,j-1)...
    *u(n,j-1)))/(1+(alpha/h^2)*(D(n,j)+D(n-1,j)+D(n,j-1)));

end
end

u(1,m) = (1-omega)*u(1,m)+...
    omega*(G(1,m)+(alpha/h^2)*(D(1,m)*(u(2,m))+D(1,m-1)...
    *u(1,m-1)))/(1+(alpha/h^2)*(D(1,m)+D(1,m-1)));

if n>2
for i=2:n-1

u(i,m) = (1-omega)*u(i,m)+...
    omega*(G(i,m)+(alpha/h^2)*(D(i,m)*(u(i+1,m))+D(i-1,m)*u(i-1,m)+D(i,m-1)...
    *u(i,m-1)))/(1+(alpha/h^2)*(D(i,m)+D(i-1,m)+D(i,m-1)));

end
end

u(n,m) = (1-omega)*u(n,m)+...
    omega*(G(n,m)+(alpha/h^2)*(D(n-1,m)*u(n-1,m)+D(n,m-1)...
    *u(n,m-1)))/(1+(alpha/h^2)*(D(n-1,m)+D(n,m-1)));

end