CIEDE2000 色差公式C++及MATLAB实现

主要功能包括：

输入处理：

接收两个LAB颜色值：[L1, a1, b1] 和 [L2, a2, b2]

LAB颜色空间模拟人眼感知，包含：

L：亮度（0-100）

a：红绿轴（负值绿，正值红）

b：黄蓝轴（负值蓝，正值黄）

核心计算：

色度调整：通过G因子补偿人眼对中性色的感知差异

色调角处理：处理角度跨越360°的特殊情况

差异计算：

ΔL’：亮度差异

ΔC’：调整后的色度差异

ΔH’：色调差异（转换为欧几里得距离）

感知权重：

Sl：亮度差异的感知权重

Sc：色度差异的感知权重

Sh：色调差异的感知权重

旋转项：处理蓝色区域的特殊感知特性

输出结果：

del_C00：色相差分量（仅包含色度和色调差异）

del_E00：总色差值（包含亮度、色度和色调差异）

#include <vector>
#include <cmath>// 定义圆周率常量
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endifstd::vector<double> CIEDE2000Computation(const std::vector<double>& arr) {// =============================================// 步骤1: 提取输入参数// =============================================// 输入数组包含两个LAB颜色值: [L1, a1, b1, L2, a2, b2]double L1 = arr[0];  // 第一个颜色的L分量（亮度）double a1 = arr[1];  // 第一个颜色的a分量（红绿轴）double b1 = arr[2];  // 第一个颜色的b分量（黄蓝轴）double L2 = arr[3];  // 第二个颜色的L分量double a2 = arr[4];  // 第二个颜色的a分量double b2 = arr[5];  // 第二个颜色的b分量// =============================================// 步骤2: 初始化常量// =============================================double kL = 1.0;  // 亮度权重因子double kC = 1.0;  // 色度权重因子double kH = 1.0;  // 色调权重因子// =============================================// 步骤3: 计算初始色度值// =============================================// 计算两个颜色的原始色度（a和b分量的欧几里得距离）double c1 = std::sqrt(a1*a1 + b1*b1);double c2 = std::sqrt(a2*a2 + b2*b2);double c_bar = (c1 + c2) / 2.0;  // 平均色度// =============================================// 步骤4: 计算色度调整因子G// =============================================// G因子用于补偿人眼对中性色的感知差异double G = 0.5 * (1 - std::sqrt(std::pow(c_bar, 7) / (std::pow(c_bar, 7) + std::pow(25.0, 7)));// =============================================// 步骤5: 计算调整后的a分量// =============================================// 应用G因子调整a分量，以补偿感知非线性double a1p = (1 + G) * a1;  // 调整后的第一个颜色a分量double a2p = (1 + G) * a2;  // 调整后的第二个颜色a分量// =============================================// 步骤6: 计算调整后的色度值和色调角// =============================================// 基于调整后的a分量重新计算色度double c1p = std::sqrt(a1p*a1p + b1*b1);  // 调整后的第一个颜色色度double c2p = std::sqrt(a2p*a2p + b2*b2);  // 调整后的第二个颜色色度// 使用atan2函数计算色调角（弧度制）double h1p = std::atan2(b1, a1p);  // 第一个颜色的色调角double h2p = std::atan2(b2, a2p);  // 第二个颜色的色调角// 确保色调角在[0, 2π)范围内if (h1p < 0) h1p += 2 * M_PI;if (h2p < 0) h2p += 2 * M_PI;// =============================================// 步骤7: 计算色调角差值// =============================================double dhp = h2p - h1p;  // 原始色调角差// 处理特殊情况：当任一颜色为中性色时（色度接近0）if (c1p * c2p == 0) {dhp = 0;  // 中性色没有色调差异}// 处理色调角跨越360°边界的情况else if (std::abs(dhp) > M_PI) {if (dhp > 0) {dhp -= 2 * M_PI;  // 正向跨越边界} else {dhp += 2 * M_PI;  // 负向跨越边界}}// =============================================// 步骤8: 计算ΔL', ΔC', ΔH'// =============================================double del_L_p = L2 - L1;  // 亮度差异double del_C_p = c2p - c1p;  // 调整色度差异// 色调差异（转换为色度空间中的距离）double del_H_p = 2 * std::sqrt(c1p * c2p) * std::sin(dhp / 2);// =============================================// 步骤9: 计算平均色调角// =============================================double h_bar = 0;  // 平均色调角初始化// 仅当两个颜色都有色调时计算平均色调角if (c1p * c2p != 0) {if (std::abs(h2p - h1p) <= M_PI) {// 当色调角差小于180°时直接平均h_bar = (h1p + h2p) / 2;} else {// 处理跨越360°边界的平均计算if (h1p + h2p < 2 * M_PI) {h_bar = (h1p + h2p + 2 * M_PI) / 2;} else {h_bar = (h1p + h2p - 2 * M_PI) / 2;}}}// =============================================// 步骤10: 计算中间平均值// =============================================double L_bar_p = (L1 + L2) / 2;    // 平均亮度double C_bar_p = (c1p + c2p) / 2;  // 平均调整色度// =============================================// 步骤11: 计算色调相关参数T// =============================================// T因子表示色调角位置的权重调整double T = 1.0 - 0.17 * std::cos(h_bar - M_PI/6.0)+ 0.24 * std::cos(2*h_bar)+ 0.32 * std::cos(3*h_bar + M_PI/30.0)- 0.20 * std::cos(4*h_bar - 63.0*M_PI/180.0);// =============================================// 步骤12: 计算Δθ和Rc// =============================================double h_bar_deg = h_bar * 180.0 / M_PI;  // 弧度转角度// 计算基于色调位置的旋转因子double del_theta_deg = 30.0 * std::exp(-std::pow((h_bar_deg - 275.0)/25.0, 2));// 计算色度压缩因子double Rc = 2.0 * std::sqrt(std::pow(C_bar_p, 7) / (std::pow(C_bar_p, 7) + std::pow(25.0, 7)));// =============================================// 步骤13: 计算权重因子// =============================================// 亮度差异的权重因子（对中间亮度更敏感）double Sl = 1.0 + (0.015 * std::pow(L_bar_p - 50, 2)) / std::sqrt(20 + std::pow(L_bar_p - 50, 2));// 色度差异的权重因子double Sc = 1.0 + 0.045 * C_bar_p;// 色调差异的权重因子double Sh = 1.0 + 0.015 * C_bar_p * T;// =============================================// 步骤14: 计算旋转项Rt// =============================================// 该项处理蓝色区域的感知非线性double Rt = -std::sin(2 * del_theta_deg * M_PI / 180.0) * Rc;// =============================================// 步骤15: 计算最终结果// =============================================// 计算总色差ΔE00 (CIEDE2000色差)double term1 = std::pow(del_L_p / (kL * Sl), 2);  // 亮度差异项double term2 = std::pow(del_C_p / (kC * Sc), 2);  // 色度差异项double term3 = std::pow(del_H_p / (kH * Sh), 2);  // 色调差异项double cross_term = Rt * (del_C_p / (kC * Sc)) * (del_H_p / (kH * Sh));  // 交叉项double del_E00 = std::sqrt(term1 + term2 + term3 + cross_term);  // 总色差// 计算色相差分量ΔC00 (不包括亮度差异)double del_C00 = std::sqrt(term2 + term3 + cross_term);  // 色相+色度差异// 返回结果向量 [ΔC00, ΔE00]return {del_C00, del_E00};
}

MATLAB

function [rlst] = CIEDE2000Computation(arr)
%CIEDE2000COMPUTATION 此处显示有关此函数的摘要,
%   此处显示详细说明 L1,a1,b1,L2,a2,b2
%   Constants for CIEDE2000 calculation
L1 = arr(1);
a1 = arr(2);
b1 = arr(3);
L2 = arr(4);
a2 = arr(5);
b2 = arr(6);
kL = 1; kC = 1; kH = 1;%compute C1,C2,h'1,h'2
c1 = sqrt(a1^2 + b1^2);%compute C1_ab
c2 = sqrt(a2^2 + b2^2);%compute C2_ab
c_bar = (c1 + c2) / 2;%compute Cave_ab
G = 0.5*(1 - sqrt((c_bar^7)/(c_bar^7 + 25^7)));%compute G
a1p = (1 + G)*a1;%computer a1'
a2p = (1 + G)*a2;%computer a2'
c1p = sqrt(a1p^2 + b1^2);%computer C1'
c2p = sqrt(a2p^2 + b2^2);%computer C2'
%c_bar_p = (c1p + c2p) / 2;%%%compute C'_bar
h1p = atan2(b1,a1p);%compute h1'
h2p = atan2(b2,a2p);%compute h2'
%%h_bar_p = (h1p + h2p) / 2;%%%not clear%compute ΔL',ΔC',ΔH'
del_L_p = L2 - L1;%compute ΔL'
del_C_p = c2p - c1p;%compute ΔC’
dhp = h2p - h1p;%compute h2' - h1'
if c1p*c2p == 0; dhp = 0;%if C1' × C2' = 0,Δh' = 0;
elseif abs(dhp) <= pi; dhp = dhp;%if |Δh' ≤ 180°|,Δh' = dhp;
elseif dhp > pi; dhp = dhp - 2*pi;%if Δh'＞180°,Δh' = Δh' - 360°;
elseif dhp < -pi; dhp = dhp + 2*pi;%if Δh'＜-180°,Δh' = Δh' +360°;
end
del_H_p = 2 * (c1p*c2p)^(1/2)*sin(dhp/2);%%compute h'_bar
h_bar = h1p + h2p;%compute h'_sum
if c1*c2 == 0; h_bar = 0;
elseif abs(h1p - h2p) <= pi; h_bar = h_bar / 2;
elseif ((abs(h2p - h1p) > pi)&&((h1p + h2p) < 2*pi)); h_bar = ((h2p + h1p) + 2*pi)/2;
elseif ((abs(h2p - h1p) > pi)&&((h1p + h2p) >= 2*pi)); h_bar = ((h2p + h1p) - 2*pi)/2;
end
%%endL_bar_p = (L1 + L2) / 2;%compute Lave'
C_bar_p = (c1p + c2p) / 2;%compute Cave'
Tmp = 1 - 0.17*cos(h_bar - pi/6) + 0.24*cos(2*h_bar) + 0.32*cos(3*h_bar + pi/30) - 0.20*cos(4*h_bar - 63*pi/180);%compute T% Calculate delta C00
del_theta = 30 * exp(-((h_bar-275*180/pi)/25)^2);
Rc = 2 * (C_bar_p^7/(C_bar_p^7+25^7))^(1/2);
Sl = 1 + 0.015*(L_bar_p - 50)^2/(20+(L_bar_p-50)^2)^(1/2);
Sc = 1 + 0.045 * C_bar_p;
Sh = 1 + 0.015 * C_bar_p * Tmp;
Rt = -sin(2 * del_theta) * Rc;
del_E00 = sqrt((del_L_p/(kL*Sl))^2 + (del_C_p/(kC*Sc))^2 + (del_H_p/(kH*Sh))^2 + Rt*((del_C_p/(kC*Sc))*(del_H_p/(kH*Sh))));
del_C00 = sqrt((del_C_p/(kC*Sc))^2 + (del_H_p/(kH*Sh))^2 + Rt*((del_C_p/(kC*Sc))*(del_H_p/(kH*Sh))));rlst = [del_C00,del_E00];
end