CAM16UCS_deltaE.m

function [dE,Jab1,Jab2,name] = CAM16UCS_deltaE(rgb1,rgb2,varargin)
% Calculate perceptual color difference (deltaE) between sRGB colors using CAM16-UCS.
%
%%% Syntax %%%
%
%   dE = CAM16UCS_deltaE(rgb1,rgb2)
%   dE = CAM16UCS_deltaE(rgb1,rgb2,...,coords)
%   dE = CAM16UCS_deltaE(rgb1,rgb2,...,output)
%   dE = CAM16UCS_deltaE(rgb1,rgb2,...,<opts>)
%   [dE,Jab1,Jab2] = CAM16UCS_deltaE(...)
%
%% Examples %%
%
%   >> rgb1 = [64,128,255]/255;
%   >> rgb2 = [128,64,200]/255;
%   >> dE = CAM16UCS_deltaE(rgb1,rgb2)
%   dE =
%        10.7747
%
%   >> rgb1 = uint8([64,128,255; 100,150,200]);
%   >> rgb2 = uint8([128,64,200]);
%   >> dE = CAM16UCS_deltaE(rgb1,rgb2)
%   dE =
%        10.7747
%        12.6697
%
%   >> dE = CAM16UCS_deltaE(rgb1,rgb2,'JCh')
%   dE =
%        10.7747
%        12.6697
%
%   >> dE = CAM16UCS_deltaE(rgb1,rgb2,'Jab','LCD','Euclid')
%   dE =
%        33.1934
%        42.6910
%
%% Input Arguments (**==default) %%
%
%   rgb1 = Numeric array of sRGB values. Floating point values must be
%          0<=rgb<=1, integer must be 0<=rgb<=intmax(class(rgb)).
%          Size Nx3 or RxCx3, the last dimension encodes the R,G,B values.
%   rgb2 = Numeric array of sRGB values, same format as <rgb1>. Note that
%          either <rgb1> or <rgb2> may consist of one color, which will be
%          implicitly expanded to match the colors of the other rgb array.
%   coords = StringScalar or CharRowVector, either of the following:
%          'Jab'** / 'JCh', which selects the CAM16 deltaE coordinates:
%          Jab   cartesian coordinates: lightness - red/green - yellow/blue
%          JCh cylindrical coordinates: lightness - chroma - hue angle
%   output = StringScalar or CharRowVector, either of the following:
%          'Euclid' / 'Huang'**, which select the output transform:
%          Euclid: no power transform, dE is euclidean in CAM16 colorspace.
%          Huang: applies the CAM16-UCS power transform from Li et al.
%          "Comprehensive color solutions: CAM16, CAT16, and CAM16-UCS",
%          based on the power-function approach of Huang et al.
%   <opts> = all remaining inputs are passed to CAM16UCS_parameters.
%
%% Output Arguments %%
%
%   dE = Column vector of color differences in CAM16-UCS space.
%        Size Nx1, where N is the max number of colors in <rgb1> or <rgb2>.
%   Jab1 = Numeric Nx3 matrix of the <rgb1> colors in CAM16 colorspace.
%   Jab2 = Numeric Nx3 matrix of the <rgb2> colors in CAM16 colorspace.
%
%% Dependencies %%
%
% * MATLAB R2009b or later.
% * sRGB_to_CAM16UCS.m from <https://github.com/DrosteEffect/CIECAM16>
%
% See also SRGB_TO_CAM16UCS CAM16UCS_TO_SRGB CAM16UCS_PARAMETERS

%% Input Wrangling %%
%
ddE = 'HUANG'; % default deltaE
dcs = 'JAB'; % default coordinates
for k = numel(varargin):-1:1
	try %#ok<TRYNC>
		tmp = upper(varargin{k});
		switch tmp
			case {'JAB','JCH'}
				varargin(k) = [];
				dcs = tmp;
			case {'EUCLID','HUANG'}
				varargin(k) = [];
				ddE = tmp;
		end
	end
end
%
% Convert to CAM16-UCS with isd=true for deltaE calculations.
% The conversion function performs input checking on the sRGB.
[Jab1,  ~] = sRGB_to_CAM16UCS(rgb1,true,varargin{:});
[Jab2,csn] = sRGB_to_CAM16UCS(rgb2,true,varargin{:});
%
% Reshape to Nx3
Jab1 = reshape(Jab1,[],3);
Jab2 = reshape(Jab2,[],3);
%
n1 = size(Jab1,1);
n2 = size(Jab2,1);
%
% Check size compatibility
assert(n1==n2 || n1==1 || n2==1,...
	'SC:CAM16UCS_deltaE:rgb:IncompatibleSizes',...
	'Inputs <rgb1> and <rgb2> must have the same number of rows, or one must have a single row.')
%
%% Calculate DeltaE %%
%
% Calculate Euclidean distance based on selected colorspace
switch dcs
	case 'JAB'
		% Already in J'a'b' format - use bsxfun for compatibility
		diff = bsxfun(@minus, Jab1, Jab2);
		dE = hypot(hypot(diff(:,1),diff(:,2)),diff(:,3));
	case 'JCH'
		% Convert to J'C'h
		J1 = Jab1(:,1); a1 = Jab1(:,2); b1 = Jab1(:,3);
		J2 = Jab2(:,1); a2 = Jab2(:,2); b2 = Jab2(:,3);
		C1 = hypot(a1,b1);
		C2 = hypot(a2,b2);
		h1 = atan2(b1,a1);
		h2 = atan2(b2,a2);
		dh = bsxfun(@minus, h1, h2);
		dh = mod(dh+pi,2*pi)-pi;
		dH = 2*sqrt(bsxfun(@times, C1, C2)).*sin(dh/2);
		dE = hypot(hypot(...
			bsxfun(@minus, J1, J2), ...
			bsxfun(@minus, C1, C2)), dH);
end
%
if strcmpi(ddE,'HUANG')
	switch csn
		case 'CAM16UCS'
			dE = 1.41 .* dE.^0.63;
		otherwise
			error('SC:CAM16UCS_deltaE:NoPowerTransform',...
				'The Huang power transform is only defined for UCS. Specify option "Euclid" instead.')
	end
end
%
name = strcat(csn,'_d',dcs,'_',ddE);
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CAM16UCS_deltaE
% Copyright (c) 2026 Stephen Cobeldick
%
% Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%license