function [G1,G2] = grad(U)
[n1,n2] = size(U);
G1 = [U(2:n1,:)-U(1:n1-1,:); zeros(1,n2)];
G2 = [U(:,2:n2)-U(:,1:n2-1), zeros(n1,1)];
endfunction

function Div = div(F1,F2)
[n1,n2] = size(F1);
Div = -[zeros(1,n2); F1(1:n1-1,:)] + [F1(1:n1-1,:); zeros(1,n2)] - [zeros(n1,1), F2(:,1:n2-1)] + [F2(:,1:n2-1), zeros(n1,1)];
endfunction

function I = proxROF(D,lambda,tau,n)
[n1,n2] = size(D);
F1 = zeros(n1,n2);
F2 = F1; G1 = F1; G2 = F1;
One = ones(n1,n2);
for k=1:n
   K = div(F1,F2);
   [G1,G2] = grad(tau*(K-D/lambda));
   Den = One + sqrt(G1.*G1+G2.*G2);
   F1 = (F1 + G1)./Den;
   F2 = (F2 + G2)./Den;
end
I = D-lambda*div(F1,F2);
endfunction

I = im2double(imread('barbara.pgm'));
D = imnoise(I,'gaussian');
J = proxROF(D,0.1,0.25,100);
imshow([D J])


h = im2double(imread('PSF.pgm'));
D = im2double(imread('blurred.pgm'));

function I = deblur(D,h,lambda,gamma,n,m)
[n1,n2] = size(D); [m1,m2] = size(h);
mm1 = floor(m1/2); mm2 = floor(m2/2);
hh = zeros(n1,n2);
hh(1:mm1+1,1:mm2+1) = h(mm1+1:m1,mm2+1:m2);
hh(1:mm1+1,n2-mm2+1:n2) = h(mm1+1:m1,1:mm2);
hh(n1-mm1+1:n1,1:mm2+1) = h(1:mm1,mm2+1:m2);
hh(n1-mm1+1:n1,n2-mm2+1:n2) = h(1:mm1,1:mm2);
Fh = fft(hh); Fh = Fh/Fh(1,1); FD = fft(D); I = D;
for i=1:n
  I = I - gamma*real(ifft(conj(Fh).*(Fh.*fft(I)-FD)));
  I = proxROF(I,lambda*gamma,0.25,m)
end
endfunction

J = deblur(D,h,.01,1,10,40);
imshow([D J])