CIELAB color difference

CIELAB色差计算比较详细的介绍：
https://en.wikipedia.org/wiki/Color_difference
CIE2000的计算步骤：
https://www.colortell.com/3424.html
MATLAB代码实现：

function E00 = delta_E00(img,src,bor) %返回两张图片的色差计算结果
imgt = rgb2lab(img);
srct = rgb2lab(src);
[height,width,deepth] = size(img);
fph=bor+1;
lph=bor;
fpv=bor+1;
lpv=bor;
E = 0;
kL = 1;
kC = 1;
kH = 1;
for a=fpv:height-lpv
    for b=fph:width-lph
        mean_cab = (CaiDu(imgt(a,b,2),imgt(a,b,3))...,
            + CaiDu(srct(a,b,2),srct(a,b,3))) / 2;
        mean_cab_pow7 = mean_cab^7;
        G = 0.5*(1-sqrt(mean_cab_pow7/(mean_cab_pow7 + 25^7)));
        LL1 = imgt(a,b,1);
        aa1 = imgt(a,b,2) * (1+G);
        bb1 = imgt(a,b,3);
        LL2 = srct(a,b,1);
        aa2 = srct(a,b,2) * (1+G);
        bb2 = srct(a,b,3);
        CC1 = CaiDu(aa1,bb1);
        CC2 = CaiDu(aa2,bb2);
        hh1 = SeDiaoJiao(aa1,bb1);
        hh2 = SeDiaoJiao(aa2,bb2);
        delta_LL = LL1 - LL2;
        delta_CC = CC1 - CC2;
        delta_hh = SeDiaoJiao(aa1,bb1) - SeDiaoJiao(aa2,bb2);
        delta_HH = 2 * sin(pi * delta_hh / 360) * sqrt(CC1 * CC2);
        mean_LL = (LL1 + LL2) / 2;
        mean_CC = (CC1 + CC2) / 2;
        mean_hh = (hh1 + hh2) / 2;
        SL = 1 + 0.015 * (mean_LL-50)^2 / sqrt(20 + (mean_LL-50)^2);
        SC = 1 + 0.045 * mean_CC;
        T = 1 - 0.17*cos((mean_hh-30)*pi/180)...,
            + 0.24*cos(2*mean_hh*pi/180)...,
            + 0.32*cos((3*mean_hh+6)*pi/180)...,
            - 0.2*cos((4*mean_hh-63)*pi/180);
        SH = 1 + 0.015*mean_CC*T;
        mean_CC_pow7 = mean_CC^7;
        RC = 2*sqrt(mean_CC_pow7/(mean_CC_pow7+25^7));
        delta_xita = 30*exp(-((mean_hh-275)/25)^2);
        RT = -sin((2*delta_xita)*pi/180)*RC;
        L_item = delta_LL/(kL*SL);
        C_item = delta_CC/(kC*SC);
        H_item = delta_HH/(kH*SH);
        E(a-bor,b-bor) = sqrt(L_item^2 + C_item^2....,
            + H_item^2 +RT*C_item*H_item); %CIE2000计算公式
    end
end
E00 = nanmean(mean(E)); %统计平均值，mean会产生NaN计算结果，所以用nanmean，欢迎指正
end

function cab = CaiDu(a,b)
cab = sqrt(a^2 + b^2);
end

function hab = SeDiaoJiao(a,b)
h = (180/pi)*atan(b/a);
if(a>0 && b>0)
    hab = h;
elseif(a<0 && b>0)
	hab = 180 + h;
elseif(a<0 && b<0)
	hab = 180 + h;
else
    hab = 360 + h;
end
end