%This is an "edited" implementation of the algorithm for calculating the
%Visual Information Fidelity (VIF) measure (may also be known as the Sheikh
%-Bovik Index) between two images. Please refer
%to the following paper:
%
%H. R. Sheikh and A. C. Bovik "Image Information and Visual Quality"
%IEEE Transactios on Image Processing, in publication, May 2005.
%Download manuscript draft from http://live.ece.utexas.edu in the
%Publications link
%
% For questions and/or comments please email
% Hisham Rahman (hishamrahman@gmail.com).

% Compute the pixel wise self-information I(D,D'|S=s) for a given subband
% coefficients.
% imgRefSb is HxW in [0, 255] double-precision.

function [imgSelfInfo] = selfinfoQMapFunc_PW(imgRefSb)
    % Tunable parameters 
    sigma_nsq = 0.4;
    M = 3; %Block Size
    
    % Calculate parameters of reference image
    sizey=size(imgRefSb);
    
    % Collect MxM blocks. Rearrange each block into an
    % M^2 dimensional vector and collect all such vectors.
    % Collece ALL possible MXM blocks (even those overlapping) from the subband
    temp=[];
    for j = 1:M
        for k = 1:M
            temp=cat(1, temp, reshape(imgRefSb(k: end-(M-k), j: end-(M-j)), 1, []));
        end
    end
    
    % estimate mean and covariance
    mcu = mean(temp')';
    cu = ((temp-repmat(mcu, 1, size(temp, 2)))*(temp-repmat(mcu, 1, size(temp, 2)))')./size(temp, 2); % covariance matrix for U
    
    % calculate the S field
    % Pad (M-1)/2 elements to each side using mirror boundary conditions.
    pl = (M-1)/2;
    imgRefSb_pad = padarray(imgRefSb, [pl, pl], 'symmetric');
    temp = zeros(M*M, size(imgRefSb, 1)*size(imgRefSb, 2));
    for i = 1:size(imgRefSb, 2)
        for j = 1:size(imgRefSb, 1)
            temp(:, (i-1)*size(imgRefSb, 1)+j) = reshape(imgRefSb_pad(j:j+M-1, i:i+M-1), 1, [])';
        end
    end
    
    % Calculate the S field
    ss = (inv(cu)*temp);
    ss = sum(ss.*temp,1)./(M*M);
    ss = reshape(ss,sizey);
    
    % Eigen-decomposition
    [v, d] = eig(cu);
    larr(1, :) = diag(d)';
    
    % Compute self information.
    lambda = larr(1, :);

    j = 1;
    imgSelfInfo = log2(1+ss.*lambda(j)./(sigma_nsq));
    for j = 2:length(lambda)
        imgSelfInfo = imgSelfInfo + (log2(1+ss.*lambda(j)./(sigma_nsq)));
    end
end

