function [J]=diffusion_lag_deblur(I,iter,lam,K,g,h)
%% Deblurring of image I by lagged diffusion with diffusivity function g and psf-kernel K
%%       [J]=diffusion_lag_deblur(I,iter,lam,K,g,h)
%%
%% Input: I    - image (double array gray level 0-255),
%%        iter - num of iterations,
%%        lam  - fidelity term lambda; choosing the correct value
%%               depends on the size of the image and the range of 
%%               pixel values; for the range issue, both commands
%%                 J1 = diffusion_lag(I,iter,dt,lam,G,I0)
%%                 J2 = diffusion_lag(I/a,iter,dt,a*lam,G,I0)
%%               give the same result; for the size issue,
%%               images whose pixel values are in [0..1] and 
%%               whose spatial stepsize is h,
%%               the value  lam=h*400/p  and TV-minimization
%%               (see lag_TV.m) will shrink the jump
%%               in an image with black left half (I=0) and 
%%               white right half (I=1) by precisely p percent.
%%        K    - sparse matrix of dimension prod(size(I)) 
%%               which defines the known psf of the blur.
%%               K should have row-sums equal to 1.
%%        g    - name of the m-file for minimization method [g='g_TV']
%%        h    - spatial stepsize [h=1/max(size(I))]
%%       (default values are in [])
%% Output: evolved image
%%
%%   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=5;
end
if ~exist('h')
    h=1/max([ny,nx]);
end
if ~exist('lam')
    lam=400; % shrink jumps in TV-denoising by at most 1 percent
end
if ~exist('K')
    K=speye(nx*ny,nx*ny); % no blurring
end
if ~exist('g')
    g='g_TV';
end

% right hand side
f=h^2*lam*K'*I(:);
% normalization of K
K=h*sqrt(lam)*K;  

% prepare for plots
% hFig=figure;
% res=get(0,'ScreenSize');
% set(hFig,'Position',[1 1 res(3) res(4)-64]); %Leave room for title bar

% prepare sparse matrix for difference method
i=zeros(1,8*nx*ny);
j=zeros(1,8*nx*ny);
s=zeros(1,8*nx*ny);
fak=zeros(1,8*nx*ny);
indk=zeros(1,8*nx*ny); % index into matrix g_lag
indl=zeros(1,8*nx*ny);
num=0;
% construct index vectors and signs for diffusion_lag entries in the
% matrix
for k=1:nx-1
    kk=(k-1)*ny;
    for l=1:ny-1
        kkl=kk+l; % row number for diagonal entry 2*g(l,k)
        num=num+1;
        i(num)=kkl;
        j(num)=kkl;
        fak(num)=2;
        indk(num)=k;
        indl(num)=l;
        num=num+1;
        i(num)=kkl+1;
        j(num)=kkl+1;
        fak(num)=1;
        indk(num)=k;
        indl(num)=l;
        num=num+1;
        i(num)=kkl+ny;
        j(num)=kkl+ny;
        fak(num)=1;
        indk(num)=k;
        indl(num)=l;
        num=num+1;
        i(num)=kkl;
        j(num)=kkl+1;
        fak(num)=-1;
        indk(num)=k;
        indl(num)=l;
        num=num+1;
        i(num)=kkl;
        j(num)=kkl+ny;
        fak(num)=-1;
        indk(num)=k;
        indl(num)=l;
        num=num+1;
        i(num)=kkl+1;
        j(num)=kkl;
        fak(num)=-1;
        indk(num)=k;
        indl(num)=l;
        num=num+1;
        i(num)=kkl+ny;
        j(num)=kkl;
        fak(num)=-1;
        indk(num)=k;
        indl(num)=l;
    end
end
numend=num;
i=i(1:numend);
j=j(1:numend);
s=s(1:numend);
fak=fak(1:numend);
indtotal=1:numend;
% convert indk,indl to single index 
indg=sub2ind([ny,nx],indl(indtotal),indk(indtotal));

for k=1:iter  % do iterations
    fprintf('  diffusion_lag_deblur: start iteration  %d\n',k)
    % nabla and div are discretized as in Chambolle's dual method
    % estimate nabla(I0)
    I_x = [I(:,[2:nx]),zeros(ny,1)]-[I(:,[1:nx-1]),zeros(ny,1)];
	I_y = [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(I_x.^2+I_y.^2);
    arg_lag=arg_lag(:).';
    g_lag=feval(g,arg_lag);    
    % construct sparse matrix for difference scheme
    s(indtotal)=fak(indtotal).*g_lag(indg);
    A=sparse(i,j,s,nx*ny,nx*ny,numend);
    [I,iter]=my_cg_AplusKstarK(A,K,f,I(:),1e-5);  % the new image 
    fprintf('                        %d Iterationen des cg-Verfahrens\n',iter)
    I=reshape(I,ny,nx);
   % imshow(I);
   % pause
end % for i
%% return image
%J=I*Imean/mean(mean(I)); % normalize to original mean
J=I;