function [ mImageApprox, mNCn, mCn, mICn, actualPSNR, processingTime, mError ] = ImageApproximation( mImage, mAtoms, algorithm, criteria, blockWidth )
% ImageApproximation Returns an approximation of an image generated by choosing atoms from
% dictionary using either thresholding or a greedy algorithm.
%
% Usage:            [ mImageApprox, mNCn, mCn, mICn, actualPSNR, processingTime, ...
%                       mError ] = ImageApproximation( mImage, mAtoms, algorithm, ...
%                       criteria, blockWidth );
%
% Inputs:
%   mImage          the matrix of pixel intensity values
%   mAtoms          matrix whos columns are the atoms from the dictionary we use to
%                   approximate the image with
%   algorithm       name of the algorithm to use
%   criteria        the target psnr between the original image and its approximation if using
%                   a greedy algorithm or the threshold if using Thresholding
%   blockWidth      width of the square blocks that the image will be processed in
%
% Outputs:
%   mImageApprox    matrix representing the approximated image
%   mNCn            matrix containing the number of coefficients used to represent each
%                   block of the image
%   mICn            array containing the index's of the retained coefficients
%   actualPSNR      psnr between the original image and its approximation
%   processingTime  time in seconds to process the image
%   mError          matrix containing the norm of the error between the original image and
%                   the approximated image for each block processed.

% See   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/

% Get the dimensions of the image
[y x] = size(mImage);
[szSpace nAtoms] = size(mAtoms);


%<<<<<<<<<<<<<<<<<<<<< Validate the input >>>>>>>>>>>>>>>>>>>>>>>>%

% If we are using a greedy algorithm
if ~strcmp(algorithm,'Thresholding')
    
    if ~Span(mAtoms, szSpace, nAtoms)
        
        error('Atoms do not span the space!!');
        
    end
    
    blockWidth = sqrt(szSpace);
    if ~IsInteger(blockWidth)
        
        error('Supplied dictionary is not of correct dimensions');
        
    end
    
    % If vectors in mAtoms are orthogonal then use Thresholding and exit as
    % the we have no way of knowing if the algorithmCriteria is correct.
    if ~IsOrthogonal(mAtoms, szSpace, nAtoms)
             
        % Calculate the per block abs error from the PSNR for the whole image
        targetAbsError = PSNRToAbsError( criteria, y, x );
        algorithmCriteria = ImageToBlockError( targetAbsError, y, x, blockWidth );        
                
    else
        
        fprintf('Supplied orthogonal matrix for use with %s!\n',algorithm); 
        error('Using thresholding instead!!');
        
    end
            
else
    
    % If we are Thresholding using the seperable product then we will have
    % nAtoms^2.
    if blockWidth == szSpace
        
        nAtoms = nAtoms^2;
        
    end
    
end

% Check that the image can be split into square blocks of width blockWidth
if ~IsInteger(y/blockWidth) || ~IsInteger(x/blockWidth)
    
    error('Image cannot be divided exactly into blocks');
    
end

% Check algorithm exists in current search path
if (exist(algorithm)) == 0
    
    error('Algorithm is not in search path!!')
    
end

fprintf('Path to algorithm: %s\n',which(algorithm));

%<<<<<<<<<<<<<<<<<<<<< End validation >>>>>>>>>>>>>>>>>>>>>>>>>>>>%


% Split the image into a cell array of square blocks
cImage = mat2cell(mImage, blockWidth*ones(1, (y/blockWidth)), blockWidth*ones(1, ...
    (x/blockWidth)));
[yC xC] = size(cImage);

% Initialize matrices
mICn = zeros(nAtoms,yC,xC);
mCn = zeros(nAtoms,yC,xC);
mNCn = zeros(yC,xC);
mError = zeros(yC,xC);
cImageApprox = cell(yC,xC);


%<<<<<<<<<<<<<<<<<<<< Approximate each block of the image >>>>>>>>>>>>>>>>%

tic;
for i = 1:yC
    
        for j = 1:xC
            
            [ mImageBlockApprox, vCn, mNCn(i,j), vICn,mError(i,j) ] = ProcessBlocks( ...
                cImage{i,j}, mAtoms, algorithm, algorithmCriteria);
            mICn(1:mNCn(i,j),i,j) = vICn(1:mNCn(i,j));
            mCn(1:mNCn(i,j),i,j) = vCn(1:mNCn(i,j));            
            cImageApprox{i,j} = mImageBlockApprox;
   
        end
        
end
processingTime = toc;

%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%

mImageApprox = cell2mat(cImageApprox);
actualPSNR = CalcPSNR( mImage, mImageApprox );

end      


%<<<<<<<<<<<<<<<<<<<<<<<<< Begin Sub Functions >>>>>>>>>>>>>>>>>>>>>>>>>>>%

function integer = IsInteger(x)

    integer = false;
    
    if (round(x) - x) == 0
        
        integer = true;
        
    end
    
end
        
function [ mImageApprox, vCn, nCn, vICn, error ] = ProcessBlocks( mImage, mAtoms, ...
    algorithm, approxCriteria )
                
    % Ensures that the same values get returned
    if strcmp(algorithm,'Thresholding')
        
        [ mImageApprox, vICn, vCn, nCn ] = ...
            Thresholding( mImage, mAtoms, criteria );
        
    else                
                
        % Initialize mImageApprox to return a square matrix
        mImageApprox = zeros(size(mImage));
        fHandle = str2func(algorithm);
        [ mChosenAtoms, vICn, beta, vCn, Q]=fHandle(mImage(:)',mAtoms,approxCriteria);
        nCn = numel(vCn);
        mImageApprox(:) = mChosenAtoms(:,1:nCn)*vCn(1:nCn)';
        
    end
    
    error = norm(mImage(:) - mImageApprox(:));

end


function spaceSpanned = Span(mAtoms, szSpace, nAtoms)

    spaceSpanned = true;
    
    if szSpace > nAtoms || rank(mAtoms) < szSpace
        
        spaceSpanned = false;         
        
    end
    
end


function orthogonal = IsOrthogonal(m, y, x)

    orthogonal = false;
    tol = 1e-10;
    
    if x <= y
        
        % check the inner product
        if norm(m*m' - eye(y,x)) < tol
            
            orthogonal = true;
        end
        
    end

end

    
function blockError = ImageToBlockError( imageError, y, x, blockWidth )

    c = (imageError^2)/(y*x);c = (imageError^2)/(y*x);
    blockError = sqrt(c*blockWidth^2);

end


function tol = PSNRToAbsError( PSNR, y, x, maxIntesity )

    % If not specified assume 8 bit intensiy:
    % min 2^0 - 1 = 0 - black
    % max 2^8 - 1 = 255 - white
    if nargin < 4
        
        maxIntensity = 255;
        
    end
    
    tol = sqrt((maxIntensity^2*y*x)/(10^(PSNR/10)));

end