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 = proxTV(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

