% =================================================================
%  Sparse Feature Fidelity (SFF) Version 2.0
%  Copyright(c) 2013  Hua-wen Chang
%  All Rights Reserved.
%  changhuawen@gmail.com
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is here
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
% ----------------------------------------------------------------------
% Please refer to the following paper
%
% Hua-wen Chang, Hua Yang, Yong Gan, and Ming-hui Wang, "Sparse Feature Fidelity
% for Perceptual Image Quality Assessment", IEEE Transactions on Image Processing,
% vol. 22, no. 10, pp. 4007-4018, October 2013
% ----------------------------------------------------------------------
% INPUT variables:
% Ir:       reference image
% Id:       distorted image
% W:        feature detector (in the file 'W.mat'
%                    or can be generated by running 'TrainW' in \Training.)
% 
% Example:
% load('W.mat'); % load feature detector W
% score = SFF(refImage, disImage, W); 
% =================================================================

function SFFscore = SFF(Ir, Id, W)

patchSize = size(W,1); %8
patchDim = (patchSize^2)*3;
Cm = 0.001;
C = 0.08;
Lamda = 0.8;
Tv = 0.4;
Tx = 1;
Tm = 1;

%%%%%% DIVIDING EACH IMAGE INTO BLOCKS %%%%%%
startPosition = 1;
sizeY = size(Ir,1); sizeX = size(Ir,2);
gridY = startPosition : patchSize : sizeY-patchSize; %
gridX = startPosition : patchSize : sizeX-patchSize; %
Y = length(gridY);  X = length(gridX);
Xr = zeros(patchDim, Y*X);
Xd = zeros(patchDim, Y*X);
ij = 0;
for i = gridY;
    for j = gridX
        ij = ij+1;
        Xr(:,ij) = reshape( Ir(i:i+patchSize-1, j:j+patchSize-1, 1:3), [patchDim 1] );
        Xd(:,ij) = reshape( Id(i:i+patchSize-1, j:j+patchSize-1, 1:3), [patchDim 1] );
    end
end


Xr = double(Xr);        Xd = double(Xd);
mXr = mean(Xr);         mXd = mean(Xd);     % mean value
% remove mean value
Xr = Xr-ones(size(Xr,1),1)*mXr;
Xd = Xd-ones(size(Xd,1),1)*mXd;  


% ==========  [FEATURE SIMILARITY]  ============
% ---------- SELECTION OF PATCH PAIRS ----------
Xe = mean(abs(Xr-Xd));
medXe = Tx*median(Xe);
Y_dis = Xd(:,Xe>=medXe);
Y_ref = Xr(:,Xe>=medXe);

% ------------- FEATURE EXTRACTION ------------
A = W*Y_ref; B = W*Y_dis;

% ------------- VISUAL THRESHOLD --------------
responseIntensity = sum(A.^2);
VT = Tv*mean(responseIntensity);
A_th = A(:,responseIntensity>VT);
B_th = B(:,responseIntensity>VT);

% ------------- FEATURE SIMILARITY ------------
C1 = A_th.*B_th;
C2 = A_th.^2+B_th.^2;
SFF_f = mean2((2*C1+C)./(C2+C));



% =========  [LUMINANCE CORRELATION]  ==========
% ------ SELECTION OF MEAN VALUE PAIRS ------
mXe = abs(mXr-mXd);
mmXe = Tm*median(mXe);
mXd = mXd(:,mXe>=mmXe);
mXr = mXr(:,mXe>=mmXe);

% ---------- LUMINANCE CORRELATION ----------
meanXr = mXr-mean(mXr);
meanXd = mXd-mean(mXd);
SFF_m = (sum(meanXr.*meanXd)+Cm) / (sqrt(sum(meanXr.^2)*sum(meanXd.^2))+Cm);



%%%%%%%%%%%% QUAITY SCORE %%%%%%%%%%%
SFFscore = Lamda*SFF_m+(1-Lamda)*SFF_f;
