tvmm_a.m

Total variation image deconvolution_A majorization-minimization approach

    [x_new,n_it_cg,cost_i,isnri]=cg_tv(Wij,img_noisy,lambda,x_old,h,cg_iter,cg_thrld,displayIm,info_int,image);
    wij=diffh(x_new).^2+diffv(x_new).^2;
    t2=clock;
    % Keep track of internal data
    if bitand(info_int,1)==1
        energyi=[energyi cost_i];
    end
    if bitand(info_int,2)==2
        ISNRi=[ISNRi isnri];
    end

    % Display external simulation info
    if (info_ext>0)

        fprintf('BOA %g of %g: CG iterations = %g',i,boa_iter,n_it_cg);

        if bitand(info_ext,1)==1
            energyo=[energyo norm(img_noisy-mH(x_new,h),'fro')^2+lambda*sum(sum(sqrt(wij)))];
            fprintf(', energy = %g',energyo(end));
        end
        if bitand(info_ext,2)==2
            ISNRo=[ISNRo 10*log10(norm(img_noisy-image,'fro').^2/norm(x_new-image,'fro').^2)];
            fprintf(', ISNR = %g dB',ISNRo(end));
        end
        fprintf(', time = %g secs\n',etime(t2,t1));
    end

    x_old=x_new;
    i=i+1;
end



tf=clock;
if info_ext>0 || info_int>0
    fprintf('Simulation time = %g secs \n',etime(tf,ti));
end

img_estimated=x_new;


% =============    Auxiliary functions   =================

function sol=f_diffh(x)
h=[0 1 -1];
sol=conv2c(x,h);


function sol=f_diffht(x)
h=[-1 1];
sol=conv2c(x,h);


function sol=f_diffv(x)
h=[0 1 -1]';
sol=conv2c(x,h);


function sol=f_diffvt(x)
h=[-1 1]';
sol=conv2c(x,h);

function y = conv2c(x,h)
% Circular 2D convolution
x=wraparound(x,h);
y=conv2(x,h,'valid');


function y = wraparound(x, m)
% Extend x so as to wrap around on both axes, sufficient to allow a
% "valid" convolution with m to return the cyclical convolution.
% We assume mask origin near centre of mask for compatibility with
% "same" option.
[mx, nx] = size(x);
[mm, nm] = size(m);
if mm > mx | nm > nx
    error('Mask does not fit inside array')
end

mo = floor((1+mm)/2); no = floor((1+nm)/2);  % reflected mask origin
ml = mo-1;            nl = no-1;             % mask left/above origin
mr = mm-mo;           nr = nm-no;            % mask right/below origin
me = mx-ml+1;         ne = nx-nl+1;          % reflected margin in input
mt = mx+ml;           nt = nx+nl;            % top of image in output
my = mx+mm-1;         ny = nx+nm-1;          % output size

y = zeros(my, ny);
y(mo:mt, no:nt) = x;      % central region
if ml > 0
    y(1:ml, no:nt) = x(me:mx, :);                   % top side
    if nl > 0
        y(1:ml, 1:nl) = x(me:mx, ne:nx);            % top left corner
    end
    if nr > 0
        y(1:ml, nt+1:ny) = x(me:mx, 1:nr);          % top right corner
    end
end
if mr > 0
    y(mt+1:my, no:nt) = x(1:mr, :);                 % bottom side
    if nl > 0
        y(mt+1:my, 1:nl) = x(1:mr, ne:nx);          % bottom left corner
    end
    if nr > 0
        y(mt+1:my, nt+1:ny) = x(1:mr, 1:nr);        % bottom right corner
    end
end
if nl > 0
    y(mo:mt, 1:nl) = x(:, ne:nx);                   % left side
end
if nr > 0
    y(mo:mt, nt+1:ny) = x(:, 1:nr);                 % right side
end




function [x_new,i,iter_energy,isnri]=cg_tv(Wij,b,lambda,x_old,h,imax,cg_thrld,displayIm,info_int,image)
% Function cg_tv applies the conjugate gradient
% method in order to solve the linear system Ax=b, where
% A = (D'*V*D*lamda/2 +H'H)
%    D - toeplitz matrix
%    V - Diagonal matrix
%
% Inputs parameters are:
%    Wij : 
%    lambda  : lambda parameter
%    b  : particular solution b
%    x_old : starting value for x
%    h : blur kernel
%    imax : maximum CG iterations
%    cg_thrld : threshold for stoping CG
%    displayIm : display restored images along CG iterations
%    info_int : perform (or not) calculations for ISNR or energy
%    image : original image for ISNR calculations
% Output parameters are:
%    x_new : solution of conjudate gradient
%    i : number of iterations used
%    energy : energy function of all iterations 
%    isnri : Improve Signal-to-noise ratio of all iterations
%
% Written by: Joao Oliveira
% email: [e-mail]
% Date: 28/09/2005

global mH mHt diffh diffv diffht diffvt
% Inicialization
i=0;
iter_energy=[];
isnri=[];
decr=2*cg_thrld;

% Initial Conjugate gradient parameters
% System: Az=y, y=Ht*x

% y=conv2c(b,h);
y=mHt(b,h);
% Multiply x_old by [lambda/2*D'*V*D + Ht*H]
r=y-(lambda/2)*(diffht(Wij.*diffh(x_old))+diffvt(Wij.*diffv(x_old)))-mHt(mH(x_old,h),h);
d=r;
delta_new=r(:)'*r(:);
delta0=delta_new;

while(i<imax && decr>cg_thrld)
    
    i=i+1;
    % q=A*d; A=[lambda/2*D'*V*D + Ht*H]
    q=(lambda/2)*(diffht(Wij.*diffh(d))+diffvt(Wij.*diffv(d)))+mHt(mH(d,h),h);
    alfa=delta_new/(d(:)'*q(:));

    x_new=x_old+alfa*d;

    % Display internal simulation info
    if (info_int>0)
        fprintf('\t\tCG iteration %g',i);
        if bitand(info_int,1)==1
            iter_energy=[iter_energy norm(b-mH(x_new,h),'fro')^2+lambda*sum(sum(sqrt(diffh(x_new).^2+diffv(x_new).^2)))];
            fprintf(', energy = %g',iter_energy(end));
        end
        if bitand(info_int,2)==2
            isnri=[isnri 20*log10(norm(b-image,'fro')/norm(x_new-image,'fro'))];
            fprintf(', ISNR = %g dB',isnri(end));
        end
        fprintf('\n');
    end

    % Display restored image
    if (displayIm)
        figure(7);imagesc(x_new);colormap gray; drawnow;
    end
    
    r=r-alfa*q;
    delta_old=delta_new;
    delta_new=r(:)'*r(:);
    beta=delta_new/delta_old;
    d=r+beta*d;

    if i>10
        decr=abs(norm(x_new-x_old,'fro')/norm(x_old,'fro'));
    end
    x_old=x_new;
end


function  z = imshift(x,shx,shy)
%         z = imshift(x,shx,shy);
%
% circularly shift image by nx and ny 
%
% ----------- Input parameters -----------------
% x		    -> original 2D signal 
% shx		-> x-shift 
% shy		-> y-shift 
%
% ----------- Output parameters -----------------
% z 		-> Shifted image
   
[ny,nx] = size(x);
I=sqrt(-1);

% shift dft 
fs=exp(-I*2*pi/ny*shy*(0:(ny-1))')*exp(-I*2*pi/nx*shx*(0:(nx-1)));

% symmetric shift
z=real(ifft2(fft2(x).*fs));

 
function [z,i]=cgtv2(Wij,w,lambda,z,imax)
% Solve Az=Dw
% Inicialization
global diffh diffv diffht diffvt mH mHt
i=1;
[dimy,dimx]=size(z);


% Initial Conjugate gradient parameters
% System: Az=y
Dw=[diffh(w); diffv(w)];

Dt=diffht(z(1:(dimy/2),:))+diffvt(z(dimy/2+1:end,:));
r=Dw-[diffh(Dt); diffv(Dt)]-2/lambda*[Wij.*z(1:(dimy/2),:); Wij.*z(dimy/2+1:end,:)];
d=r;
delta_new=r(:)'*r(:);

% Iterations

while(i<=imax)
   
    % q=A*d;
    Dt=diffht(d(1:(dimy/2),:))+diffvt(d(dimy/2+1:end,:));
    q=[diffh(Dt); diffv(Dt)]+2/lambda*[Wij.*d(1:(dimy/2),:); Wij.*d(dimy/2+1:end,:)];
    alfa=delta_new/(d(:)'*q(:));
    z=z+alfa*d;
    r=r-alfa*q;
    delta_old=delta_new;
    delta_new=r(:)'*r(:);
    beta=delta_new/delta_old;
    d=r+beta*d;
    i=i+1;
end