tvmm_a.m

Total variation image deconvolution_A majorization-minimization approach

function [img_estimated,energyi,ISNRi,energyo,ISNRo]=tvmm_a(img_noisy,h,varargin)
% function   [img_estimated]=tvmm_a(img_noisy,h,...
%            'optional_parameter_name1',value1,'optional_parameter_name2', value2,...);
%
% Total Variation-based image deconvolution with
% a majorization-minimization approach.
%  
% Written by: Joao Oliveira, Jose Bioucas-Dias, Mario Figueiredo
% email: [e-mail]
% SITE: www.lx.it.pt/~jpaos/tvmm
% Date: 21/11/2005
%  
%        
% ========================== INPUT PARAMETERS (required) =================
% Parameter     Values
% name          and description
% ========================================================================
% img_noisy		(double) Noisy blured image of size ny. 
% h             (double) Blur kernel.
%
% ======================== OPTIONAL INPUT PARAMETERS ====================
% Parameter     Values
% name          and description
% =======================================================================
% boa_iter      (double) The number of outer loop iterations. 
%               Default: 20 
% cg_iter       (double) The max number of inner CG iterations.
%               Default: 200
% cg_thrld      (double) Conjugate Gradient threshold.
%               Default: 1e-5;
% image         (double) Original image.
% displayIm     ({'yes','no'}) if 'yes' is passed  the restored imaged
%               is displayed along the CG iterations
%               Default: 'no'
% info_energyi  ({'yes','no'}) if 'yes' is passed the isnr and energy of
%               the objective function is displayed inside the inner loop.
%               Default: 'no'
% info_energyo  ({'yes','no'}) if 'yes' is passed the energy of
%               the objective function is displayed at each outter loop.
%               Default: 'no'
% info_ISNRi    ({'yes','no'}) if 'yes' is passed the isnr is displayed.
%               Default: 'no'
%               Requires: image
% info_ISNRo    ({'yes','no'}) if 'yes' is passed the isnr is displayed.
%               Default: 'no'
%               Requires: image
% x_0           (double) initial image iteration.
%               Default: Wiener filter
% method        ({'cg','soim'}) Selects the method used to solve the linear
%               system: cg=conjugate gradient; soim=second order iterative
%               method.
%               Default: 'cg'
% l_min         (double) Minimum eigenvalue of the system for the second
%               order iterative method.
%               Default: 1e-5
% l_max         (double) Maximum eigenvalue of the system for the second
%               order iterative method.
%               Default: 1
% 
%       ===================================================================
%       The following functions can be provided in order to overwrite the 
%       internal ones. Recall that, by default, all calculations are
%       performed with circular convolutions. (see Technical Report)
%       ===================================================================
%
% mult_H        (function_handle) Function handle to the function that 
%               performes H*x. This function must accept two parameters:
%                       x : image to apply the convolution kernel h
%                       h : Blur kernel
% mult_Ht       (function_handle) Function handle to the function that 
%               performes Ht*x (Ht = H transpose). This function must 
%               accept two parameters:
%                       x : image to apply the convolution kernel h
%                       h : Blur kernel
% diffh         (function_handle) Function handle to the function that 
%               computes the horizontal differences of an image x.
% diffv         (function_handle) Function handle to the function that 
%               computes the vertical differences of an image x.
% diffht         (function_handle) Function handle to the function that 
%               computes the horizontal differences (transposed) of an 
%               image x.
% diffh         (function_handle) Function handle to the function that 
%               computes the vertical differences (transposed) of an 
%               image x.
%
%
% ====================== Output parameters ===============================
% img_estimated     Estimated image
% energyi           Energy of inner loop iterations.
%                   Required input arguments: 'info_energyi' and 'image'.
% ISNRi             Improvement Signal-to-noise ratio of inner loop iterations.
%                   Required input arguments: 'info_ISNRi' and 'image'.
% energyo           Energy of outer loop iterations.
%                   Required input arguments: 'info_energyo' and 'image'.
% ISNRo             Improvement Signal-to-noise ratio of outer loop iterations.
%                   Required input arguments: 'info_ISNRo' and 'image'.
%
% ============= EXAMPLES =========================================
% Normal use:
% [img_estimated]=tvmm_a(img_noisy,h)
% 
% Display restored images along CG iterations:
% [img_estimated]=tvmm_a(img_noisy,h,...
%                       'displayIm','yes')
%
% Perform ISNR calculations of the outer loop iterations:
% [img_estimated,energyi,ISNRi,energyo,ISNRo]=tvmm_a(img_noisy,...
%                       h,'info_SNRo','yes','image',image);
%

% ======================= DEFAULT PARAMETERS ===================
global mH mHt diffh diffv diffht diffvt
boa_iter = 20;          % Number of outer loop iterations
cg_iter = 200;          % Number of inner CG iterations
cg_thrld = 1e-5;        % CG threshold
displayIm = 0;
image=[];
info_energyi='no';
info_energyo='no';
info_ISNRi='no';
info_ISNRo='no';
info_int = 0;
info_ext = 0;
energyi=0;
ISNRi=0;
energyo=0;
ISNRo=0;
mH=@conv2c;
mHt=@conv2c;
diffh=@f_diffh;
diffv=@f_diffv;
diffht=@f_diffht;
diffvt=@f_diffvt;
x_0=[];


% ====================== INPUT PARAMETERS ========================
% Test for number of required parametres
if (nargin-length(varargin)) ~= 2
     error('Wrong number of required parameters');
end

% Read the optional parameters
if (rem(length(varargin),2)==1)
    error('Optional parameters should always go by pairs');
else
    for i=1:2:(length(varargin)-1)
        % change the value of parameter
        switch varargin{i}
            case 'boa_iter'                     % Outer loop iterations
                boa_iter = varargin{i+1};
            case 'cg_iter'                      % Inner loop iterations
                cg_iter = varargin{i+1};
            case 'displayIm'                    % display sucessive xe estimates
                if (isequal(varargin{i+1},'yes'))
                    displayIm=1;
                end
            case 'cg_thrld'                     % CG threshold
                cg_thrld = varargin{i+1};
            case 'image'                        % Original image
                image = varargin{i+1};
            case 'info_ISNRi'                   % display ISNR inside inner loop
                info_ISNRi = varargin{i+1};
            case 'info_ISNRo'                   % display ISNR on outer loop
                info_ISNRo = varargin{i+1};
            case 'info_energyi'                 % display energy inside inner loop
                info_energyi = varargin{i+1};
            case 'info_energyo'                 % display energy on outer loop
                info_energyo = varargin{i+1};
            case 'x_0'                          % Initial image iteration
                x_0 = varargin{i+1};
            case 'mult_H'                       % Function handle to perform H*x
                mH = varargin{i+1};
            case 'mult_Ht'                      % Function handle to perform Ht*x
                mHt = varargin{i+1};
            case 'diffh'                   % External function to perform 
                diffh = varargin{i+1};     % horizontal differences    
            case 'diffv'                   % External function to perform 
                diffv = varargin{i+1};     % vertical differences    
            case 'diffht'                  % External function to perform 
                diffht = varargin{i+1};    % horizontal differences (transposed)   
            case 'diffvt'                  % External function to perform 
                diffvt = varargin{i+1};    % vertical differences (transposed)  
            case 'method'                  % External function to perform 
                method = varargin{i+1};    % vertical differences (transposed)  
                
            otherwise
                % Hmmm, something wrong with the parameter string
                error(['Unrecognized parameter: ''' varargin{i} '''']);
        end;
    end;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parameters pos-validations (after load optional parameters)

if (isequal(info_ISNRi,'yes')||isequal(info_ISNRo,'yes'))&&isempty(image)
    error('ISNR calculations: original image not found');
end

% Which internal info to display...
if isequal(info_energyi ,'yes')
    info_int=info_int+1;
end
if isequal(info_ISNRi,'yes')
    info_int=info_int+2;
end

% Which external info to display...
if isequal(info_energyo ,'yes')
    info_ext=info_ext+1;
end
if isequal(info_ISNRo,'yes')
    info_ext=info_ext+2;
end

% Initializations
[hy,hx]=size(h);
[imy,imx]=size(img_noisy);

% Initialization with Wiener filter

    h_aux=h;
    h_aux(imy,imx)=0;
    h_aux=imshift(h_aux,-floor(hy/2),-floor(hx/2));
    db = fft2(h_aux);               % eigenvalues of B
    varx = var(img_noisy(:));       % approximate var of img_noisy
    aux = diffh(img_noisy);
    sigma =  median(abs(aux(:)-median(aux(:))))/.6745/2;
    x_old = real(fft2(db./(db.^2+10*sigma^2/varx).*ifft2(img_noisy)));
    fprintf('Sigma estimation = %g\n',sigma);

% Initialization is given...    
if ~isempty(x_0)
    x_old=x_0;
end

    
% Initial Improve Signal-to-noise Ratio 
if bitand(info_ext,2)==2 || bitand(info_int,2)==2
    ISNRi=[20*log10(norm(img_noisy-image,'fro')/norm(x_old-image,'fro'))];
    ISNRo=ISNRi;
end

% Initial energy
if bitand(info_int,1)==1 || bitand(info_ext,1)==1
    energyi=[norm(img_noisy-mH(x_old,h),'fro')^2+lambda*sum(sum(sqrt(diffh(x_old).^2+diffv(x_old).^2)))];
    energyo=energyi;
    ti=clock;
end

%====================== Conjugate Gradient Method =========================

i=1;
ti=clock;
wij=diffh(x_old).^2+diffv(x_old).^2;
%     sigma=sqrt(1/imy/imx*norm(img_noisy-mH(x_old,h),'fro')^2)
while(i<=boa_iter)

    t1=clock;
    if (info_ext>0)
        fprintf('BOA %g of %g \n',i,boa_iter);
    end
    Wij=1./sqrt(wij);
    %         sigma=sqrt(1/imy/imx*norm(img_noisy-mH(x_old,h),'fro')^2)
    lambda=1/2*imy*imx*sigma.^2/sum(sum(sqrt(wij)))