% This script is used to calculate the pearson linear correlation
%coefficient and root mean sqaured error after regression

%get the objective scores computed by the IQA metric and the subjective
%scores provided by the dataset
%% INPUT
% x = (differential) mean opinion scores (MOS)
% score = objective scores
%
%  OUTPUTS
% res = [PLCC  SROCC KROCC RMSError MAE ]
% and plots of all performance parameters
%%

mos = x;
y = score;

for i=1:size(y,2)
objectiveValues = y(:,i);

%initialize the parameters used by the nonlinear fitting function
% The following beta were used in VIS, TIP 2014

beta = [ max(mos), min(mos), mean(objectiveValues), 0.1, 0.1 ];


%plot objective-subjective score pairs
figure;
p = plot(objectiveValues,mos,'o');
set(p,'Color','blue','LineWidth',3);

%% fitting a curve using the data
[bayta ehat,J] = nlinfit(objectiveValues,mos,@logistic,beta); 
%given a ssim value, predict the correspoing mos (ypre) using the fitted curve
[ypre junk] = nlpredci(@logistic,objectiveValues,bayta,ehat,J);

RMSError(i) = sqrt(sum((ypre - mos).^2) / length(mos)); %root meas squared error
    PLCC(i) = corr(mos, ypre, 'type','Pearson');        %pearson linear coefficient
   SROCC(i) = corr(mos, ypre, 'type','Spearman');      %pearson linear coefficient
   KROCC(i) = corr(mos, ypre, 'type','Kendall');       %pearson linear coefficient
     MAE(i) = meanAbsError(mos,ypre); 

%draw the fitted curve
t = min(objectiveValues):(max(objectiveValues)-min(objectiveValues))/1000:max(objectiveValues);
[ypre junk] = nlpredci(@logistic,t,bayta,ehat,J);
hold on;
p = plot(t,ypre);
set(p,'Color','black','LineWidth',3);
legend('Test dataset ','Curve fitted with logistic function', 'Location','SouthWest');
xlabel('Obj');
ylabel('Sub');
% hold off;

end

res = [PLCC  SROCC KROCC RMSError MAE ]
