deltaE1(Ls,as,bs,Lr,ar,br)

function DE = deltaE1(Ls,as,bs,Lr,ar,br)
% this function calculate the Delta E according to CIE 2000 formula
%%L*a*b* sample -L*a*b* reference;
%% this is the correct algrithm for CIE 2000
%%In RUO and Rigg's paper, they were using L'a'b'. In this progamm, L*a*b*
%%were used. Be aware of this difference when comparing the results.

Lbp = (Ls+Lr)/2;
Cr = sqrt(ar^2 + br^2);
Cs = sqrt(as^2 + bs^2);
Cb = (Cr +Cs)/2;
G = (1- sqrt(Cb^7/(Cb^7 + 25^7)))/2;
arp = ar*(1+G);
asp = as*(1+G);
brp = br;
bsp = bs;
Crp = sqrt(arp^2 + brp^2);
Csp = sqrt(asp^2 + bsp^2);
Cbp = (Crp + Csp)/2;

if  arp>=0 & brp>=0           %% h ranges from 0 to 360 degrees
    hrp = 180*atan(brp/arp)/pi; %% atan returns radians, so multiply by 180/pi 
elseif arp>=0 & brp<0
    hrp = 180*atan(brp/arp)/pi + 360;
else
    hrp = 180*atan(brp/arp)/pi + 180;       
end
    

if  asp>=0 & bsp>=0           %% h ranges from 0 to 360 degrees
    hsp = 180*atan(bsp/asp)/pi; %% atan returns radians, so multiply by 180/pi 
elseif asp>=0 & bsp<0
    hsp = 180*atan(bsp/asp)/pi + 360;
else
    hsp = 180*atan(bsp/asp)/pi + 180;       
end
    

if abs(hsp-hrp)>180
    Hbp = (hrp+hsp-360)/2;
else
    Hbp = (hrp+hsp)/2;
end

    
T = 1- 0.17*cos((Hbp-30)*pi/180) + 0.24* cos(2*Hbp*pi/180) +  ...
     0.32*cos((3*Hbp+6)*pi/180) - 0.20* cos((4*Hbp - 63)*pi/180);
 
if abs(hsp- hrp)<= 180
    dhp = hsp - hrp; %%%abs why??
elseif abs(hsp - hrp)<= 180 & hsp <= hrp
    dhp = hsp - hrp + 360;
else
    dhp = hsp - hrp -360;
end
 
     
dLp = Ls - Lr;
dCp = Csp - Crp;

dHp = 2*sqrt(Crp*Csp)*sin((dhp/2)*pi/180);

SL = 1+ 0.015*(Lbp-50)^2/sqrt(20+(Lbp-50)^2);
SC = 1+ 0.045*Cbp;
SH = 1+ 0.015*Cbp*T;

dtheta = 30* exp(-((Hbp-275)/25)^2);

RC = sqrt(Cbp^7/(Cbp^7+25^7));
RT = -2*RC*sin(2*dtheta*pi/180);
KL = 1;  %% default values for KL, KC, KH
KC = 1;
KH = 1;

%% calculate the delta E using CIE 2000
DE = sqrt((dLp/(KL*SL))^2 + (dCp/(KC*SC))^2 + (dHp/(KH*SH))^2 ...
       + RT*(dCp/(KC*SC))*(dHp/(KH*SH)));