% function to compute the sCIELAB error using the specific parameters of
% the example 8.1. Thus, all parameters are fixed which is not the
% case for a general sCIELAB error calculator.
% Eric Dubois 2018-05-07
function scielab_error = mdsp_scielab_error(IMG1,IMG2);
%compute the scielab error between two color images in gamma-corrected sRGB. They must be
%the same size.
%For each image, convert to OZW, apply the filters, then convert to LAB.
%fixed parameters
%1931 CIE XYZ to sRGB -- matrix from Poynton
A_XYZtosRGB = [3.240479 -1.537150 -0.498535; ...
              -0.969256  1.875992  0.041556; ...
               0.055648 -0.204043  1.057311];
%OZW of Zhang and Wandell
A_XYZtoOZW = [  0.279  0.720 -0.107;  ...
              -0.449 0.290 -0.077;  ...
              0.086 -0.590 0.501];
A_sRGBtoOZW = A_XYZtoOZW*inv(A_XYZtosRGB);
A_OZWtoXYZ = inv(A_XYZtoOZW);
%get the ZW filters
load('mdsp_ZWfilters.mat');
Hcat=cat(3,H11,H22,H33);
bound = 20;
%image 1
IMG1_linsRGB = sRGB_gamma(IMG1); %convert to  linear sRGB
% convert to OZW basis
for C = 1:3
    IMG1_OZW(:,:,C) = A_sRGBtoOZW(C,1)*IMG1_linsRGB(:,:,1) + ...
        A_sRGBtoOZW(C,2)*IMG1_linsRGB(:,:,2) + ...
        A_sRGBtoOZW(C,3)*IMG1_linsRGB(:,:,3);
end

%Filter the three components with the three filters
for C = 1:3
   IMG1_OZW_filt(:,:,C) = imfilter(IMG1_OZW(:,:,C),Hcat(:,:,C),'symmetric','same');
end
%convert to XYZ
for C=1:3
IMG1_XYZ_filt(:,:,C) = A_OZWtoXYZ(C,1)*IMG1_OZW_filt(:,:,1) + ...
    A_OZWtoXYZ(C,2)*IMG1_OZW_filt(:,:,2) + ...
    A_OZWtoXYZ(C,3)*IMG1_OZW_filt(:,:,3);
end
%convert to LAB
IMG1_LAB_filt = reshape(XYZ2LAB(reshape(IMG1_XYZ_filt,[],3)),size(IMG1));
%image 2
IMG2_linsRGB = sRGB_gamma(IMG2); %convert to  linear sRGB
% convert to OZW basis
for C = 1:3
    IMG2_OZW(:,:,C) = A_sRGBtoOZW(C,1)*IMG2_linsRGB(:,:,1) + ...
        A_sRGBtoOZW(C,2)*IMG2_linsRGB(:,:,2) + ...
        A_sRGBtoOZW(C,3)*IMG2_linsRGB(:,:,3);
end
%Filter the three components with the three filters
for C = 1:3
   IMG2_OZW_filt(:,:,C) = imfilter(IMG2_OZW(:,:,C),Hcat(:,:,C),'symmetric','same');
end
%convert to XYZ
for C=1:3
IMG2_XYZ_filt(:,:,C) = A_OZWtoXYZ(C,1)*IMG2_OZW_filt(:,:,1) + ...
    A_OZWtoXYZ(C,2)*IMG2_OZW_filt(:,:,2) + ...
    A_OZWtoXYZ(C,3)*IMG2_OZW_filt(:,:,3);
end
%convert to LAB
IMG2_LAB_filt = reshape(XYZ2LAB(reshape(IMG2_XYZ_filt,[],3)),size(IMG2));
%compute the Euclidean error excluding the boundary
IMGerr = (IMG1_LAB_filt -IMG2_LAB_filt ).^2;
IMGerr_win = IMGerr(bound+1:size(IMG1,1)-bound,bound+1:size(IMG1,2)-bound,:);
scielab_error = sqrt(sum(sum(sum(IMGerr_win)))/(size(IMGerr_win,1)*size(IMGerr_win,2)));