function [deltaE] = getdeltaE(LABcomp,LABspec)

%% This function computes the deltaEab and deltaEcmc values
%% between the inputed LAB coordinates.
%% The input values should be of shape [3 x n].
%%
%% created: HPN 3/12/00
%% modified: JSB 3/16/00

%%%% deltaEab:
deltaEab = ((LABcomp(1,:)-LABspec(1,:)).^2 + ...
	    (LABcomp(2,:)-LABspec(2,:)).^2 + ...
	    (LABcomp(3,:)-LABspec(3,:)).^2 ).^.5;




%%%% deltaEcmc:

% Cstar (chroma):
cstar_comp = (LABcomp(2,:).^2 + LABcomp(3,:).^2 ).^.5;
cstar_spec = (LABspec(2,:).^2 + LABspec(3,:).^2 ).^.5;
deltacstar = (cstar_comp)-(cstar_spec);

% deltaL
deltaL = LABcomp(1,:) - LABspec(1,:);

% deltaHab (hue difference):
deltaHab = (deltaEab.^2 - deltaL.^2 - deltacstar.^2).^.5;

%% S values for comp
% sL_comp

[n m] = size(LABcomp(1,:));

sL_comp = zeros(1,m);
for i = 1:m
	if LABcomp(1,i) < 16 
		sL_comp(i)=.511; 
	else
		sL_comp(i)=(LABcomp(1,i)*.040957)./(1+LABcomp(1,i)*.01765);
	end
end


% sC_comp
sC_comp = ( cstar_comp*.0638./(1+cstar_comp*.0131) )+.638;

%sH_comp
f_comp = (cstar_comp.^4./( cstar_comp.^4+1900)).^.5;
h_comp = atan(LABcomp(2,:)./LABcomp(3,:));

[n m] = size(h_comp);
t_comp = zeros(1,m);

for i = 1:m
	if (h_comp(i) > (deg2rad(164))) & (h_comp(i) < (deg2rad(345)))  
		t_comp(i) = .56 + abs(cos(h_comp(i)+168).*.2);
	else
		t_comp(i) = .36 + abs(cos(h_comp(i)+35).*.4);
	end
end

sH_comp = (f_comp.*t_comp+1-f_comp).*sC_comp;

% deltaEcmc value

deltaEcmc = ((deltaL./sL_comp).^2 + (deltacstar./sC_comp).^2 ...
		+ (deltaHab./sH_comp).^2).^.5;

deltaEcmc = abs(deltaEcmc);

%%% Finally, we can create deltaE, composed of the deltaEab and the
%%% deltaEcmc values, to be returned by the function:
 
deltaE = [deltaEab; deltaEcmc];



