function [J,T]=diffusion_timestep(I,sigma,iter,g,dt,h)
%% Total Variation denoising.
%% Example: J=chambolle1(I,iter,lam,I0)
%% Input: I    - image (double array gray level 1-256),
%%        sigma - standard deviation of image noise [sigma=0.1]
%%        iter - max. number of iterations [iter=500]
%%        g    - name of the m-file for minimization method [g='glag_TV']
%%        dt   - factor of real time step  delta_t=dt*h^2  [dt=0.5] 
%%        h    - spatial stepsize [h=1/max(size(I))]
%%       (default values are in [])
%% Output: J -  evolved image
%%         T -  evolved time
%%
%%   On the aforementioned example,
%%   about 10000 Iterations give precision 3-5 * 10^(-3)  with dt=0.5

[ny,nx]=size(I); 
I=double(I);

if ~exist('iter')
    iter=500
end
if ~exist('sigma')
    sigma=0.1
end
if ~exist('h')
    h=1/max([ny,nx])
end
if ~exist('dt')
   dt=0.5 
end
if ~exist('g')
    g='glag_TV'
end

k=0;
T=0;
I0=I;
c=sqrt(nx*ny)*sigma;

while (k<iter) & (norm(I-I0,'fro')<c)
    norm(I-I0,'fro')/c
    k=k+1;
    T=T+dt;
    % nabla and div are discretized as in Chambolle's dual method
    % estimate nabla(I)
    G1 = [I(:,[2:nx]),zeros(ny,1)]-[I(:,[1:nx-1]),zeros(ny,1)];
	G2 = [I([2:ny],:);zeros(1,nx)]-[I([1:ny-1],:);zeros(1,nx)];
    % include the stepsize h for computation of the diffusivity function
	arg_lag=1/h*sqrt(G1.^2+G2.^2);
    gval=feval(g,arg_lag); 
    G1=gval.*G1;
    G2=gval.*G2;
    % estimate div(g(|nabla I|)*(nabla I))
    G1 = [G1(:,[1:nx-1]),zeros(ny,1)]-[zeros(ny,1),G1(:,[1:nx-1])];
	G2 = [G2([1:ny-1],:);zeros(1,nx)]-[zeros(1,nx);G2([1:ny-1],:)];
	I = I+dt*(G1+G2);
%    imshow(I);
%    pause
end % for i
%% return image
%J=I*Imean/mean(mean(I)); % normalize to original mean
J=I;