#include <stdio.h>
#include <math.h>
#define _USE_MATH_DEFINES
#define DEG (180 / M_PI)
#define RAD (M_PI / 180)
void XYZ2Lab(double X, double Y, double Z, double* L, double* a, double* b){
double fx, fy, fz;
double X_Xn, Y_Yn, Z_Zn;
double Xn = 95.046;
double Yn = 100.;
double Zn = 108.906;
X_Xn = X / Xn;
Y_Yn = Y / Yn;
Z_Zn = Z / Zn;
if (X_Xn > pow((6. / 29.), 3.)){
fx = pow(X_Xn, 1. / 3.);
}else{
fx = (pow((29. / 3.), 3.) * X_Xn + 16.) / 116.;
}
if (Y_Yn > pow((6. / 29.), 3.)){
fy = pow(Y_Yn, 1. / 3.);
}else{
fy = (pow((29. / 3.), 3.) * Y_Yn + 16.) / 116.;
}
if (Z_Zn > pow((6. / 29.), 3.)){
fz = pow(Z_Zn, 1. / 3.);
}else{
fz = (pow((29. / 3.), 3.) * Z_Zn + 16.) / 116.;
}
*L = 116. * fy - 16.;
*a = 500. * (fx - fy);
*b = 200. * (fy - fz);
}
void XYZ2sRGB(double X, double Y, double Z, double* R, double* G, double* B){
int i, j;
double R_tmp, G_tmp, B_tmp;
double M_inv[3][3] = {{3.240970, -1.537383, -0.498611},
{-0.969244, 1.875968, 0.041555},
{0.055630, -0.203977, 1.056972}};
R_tmp = M_inv[0][0] * X / 100. + M_inv[0][1] * Y / 100. + M_inv[0][2] * Z / 100.;
G_tmp = M_inv[1][0] * X / 100. + M_inv[1][1] * Y / 100. + M_inv[1][2] * Z / 100.;
B_tmp = M_inv[2][0] * X / 100. + M_inv[2][1] * Y / 100. + M_inv[2][2] * Z / 100.;
if(R_tmp <= 0.0031308){
*R = (12.92 * R_tmp) * 255.;
}else{
*R = (1.055 * pow(R_tmp, (1. / 2.4)) - 0.055) * 255.;
}
if(G_tmp <= 0.0031308){
*G = (12.92 * G_tmp) * 255.;
}else{
*G = (1.055 * pow(G_tmp, (1. / 2.4)) - 0.055) * 255.;
}
if(B_tmp <= 0.0031308){
*B = (12.92 * B_tmp) * 255.;
}else{
*B = (1.055 * pow(B_tmp, (1. / 2.4)) - 0.05) * 255;
}
}
void RGB2HTMLcolorcode(int R, int G, int B, char* code, int* dec){
sprintf(code, "0x%02x%02x%02x", R, G, B);
*dec = (R << 16) + (G << 8) + B;
}
double deltaEab(double L1, double a1, double b1, double L2, double a2, double b2){
return sqrt(pow(L2 - L1, 2) + pow(a2 - a1, 2) + pow(b2 - b1, 2));
}
double deltaE00(double L1, double a1, double b1, double L2, double a2, double b2){
double kL, kC, kH;
double dLp, Lbar;
double Cs_1, Cs_2, Cbar;
double ap_1, ap_2;
double Cp_1, Cp_2, Cpbar, dCp;
double hp_1, hp_2;
double dsmallhp, dHp, Hpbar;
double T;
double S_L, S_C, S_H;
double R_C, dtheta, R_T;
double term1, term2, term3, term4;
kL = 1.;
kC = 1.;
kH = 1.;
dLp = L2 - L1;
Lbar = (L1 + L2) / 2.;
Cs_1 = sqrt(pow(a1, 2.) + pow(b1, 2.));
Cs_2 = sqrt(pow(a2, 2.) + pow(b2, 2.));
Cbar = (Cs_1 + Cs_2) / 2.;
ap_1 = a1 + (a1 / 2.) * (1. - sqrt(pow(Cbar, 7.) / (pow(Cbar, 7.) + pow(25., 7.))));
ap_2 = a2 + (a2 / 2.) * (1. - sqrt(pow(Cbar, 7.) / (pow(Cbar, 7.) + pow(25., 7.))));
Cp_1 = sqrt(pow(ap_1, 2.) + pow(b1, 2.));
Cp_2 = sqrt(pow(ap_2, 2.) + pow(b2, 2.));
Cpbar = (Cp_1 + Cp_2) / 2.;
dCp = Cp_2 - Cp_1;
hp_1 = DEG * atan2(b1, ap_1);
hp_2 = DEG * atan2(b2, ap_2);
if (Cp_1 == 0. || Cp_2 == 0.){
dsmallhp = 0.;
}else if(fabs(hp_1 - hp_2) <= 180.){
dsmallhp = hp_2 - hp_1;
}else if((fabs(hp_1 - hp_2) > 180.) && (hp_2 <= hp_1)){
dsmallhp = hp_2 - hp_1 + 360.;
}else{
dsmallhp = hp_2 - hp_1 - 360.;
}
dHp = 2. * sqrt(Cp_1 * Cp_2) * sin(RAD * (dsmallhp/ 2.));
if (Cp_1 == 0. || Cp_2 == 0.){
Hpbar = hp_1 + hp_2;
}else if(fabs(hp_1 - hp_2) <= 180.){
Hpbar = (hp_1 + hp_2) / 2.;
}else if((fabs(hp_1 - hp_2) > 180.) && (hp_1 + hp_2 < 360)){
Hpbar = (hp_1 + hp_2 + 360.) / 2.;
}else{
Hpbar = (hp_1 + hp_2 - 360.) / 2.;
}
T = 1. - 0.17 * cos(RAD * (Hpbar - 30.)) + 0.24 * cos(RAD * (2. * Hpbar)) + 0.32 * cos(RAD * (3. * Hpbar + 6.)) - 0.20 * cos(RAD * (4. * Hpbar - 63.));
S_L = 1. + (0.015 * pow((Lbar - 50.), 2.)) / (sqrt(20. + pow((Lbar - 50.), 2.)));
S_C = 1. + 0.045 * Cpbar;
S_H = 1. + 0.015 * Cpbar * T;
R_C = 2. * sqrt(pow(Cpbar, 7.) / (pow(Cpbar, 7.) + pow(25., 7.)));
dtheta = 30. * exp(-1. * pow(((Hpbar - 275.) / 25.), 2.));
R_T = -1. * sin(RAD * (2 * dtheta)) * R_C;
term1 = pow(dLp / (kL * S_L), 2.);
term2 = pow(dCp / (kC * S_C), 2.);
term3 = pow(dHp / (kH * S_H), 2.);
term4 = R_T * (dCp / (kC * S_C)) * (dHp / (kH * S_H));
return sqrt(term1 + term2 + term3 + term4);
}
int main(void){
double _L, _a, _b;
double R, G, B;
char code[8];
int dec;
XYZ2Lab(30., 30., 30., &_L, &_a, &_b);
printf("%f\n%f\n%f\n", _L, _a, _b);
printf("\n");
XYZ2sRGB(30., 30., 30., &R, &G, &B);
printf("%f\n%f\n%f\n", R, G, B);
printf("\n");
RGB2HTMLcolorcode(15, 15, 15, code, &dec);
printf("%s\n", code);
printf("%d\n", dec);
printf("\n");
printf("%f\n", deltaEab(83, -50, -4, 55, -15, 7));
printf("%f\n", deltaE00(83, -50, -4, 55, -15, 7));
return 0;
}
EXCELでΔE00の計算シートも作成、Labの色も表示するマクロコード追加
DeltaE.xlsm - Google ドライブ
以下VBAコードも
Function D65wp()
Dim result(2) As Double
result(0) = 95.046
result(1) = 100#
result(2) = 108.906
D65wp = result
End Function
Function D65spectrum()
Dim D65(32) As Double
D65(0) = 49.9755
D65(1) = 54.6482
D65(2) = 82.7549
D65(3) = 91.486
D65(4) = 93.4318
D65(5) = 86.6823
D65(6) = 104.865
D65(7) = 117.008
D65(8) = 117.812
D65(9) = 114.861
D65(10) = 115.923
D65(11) = 108.811
D65(12) = 109.354
D65(13) = 107.802
D65(14) = 104.79
D65(15) = 107.689
D65(16) = 104.405
D65(17) = 104.046
D65(18) = 100#
D65(19) = 96.3342
D65(20) = 95.788
D65(21) = 88.6856
D65(22) = 90.0062
D65(23) = 89.5991
D65(24) = 87.6987
D65(25) = 83.2886
D65(26) = 83.6992
D65(27) = 80.0268
D65(28) = 80.2146
D65(29) = 82.2778
D65(30) = 78.2842
D65(31) = 69.7213
D65(32) = 71.6091
D65spectrum = D65
End Function
Function xyz_func_10deg()
Dim xyz(32, 3) As Double
xyz(0, 0) = 0.00016
xyz(1, 0) = 0.002362
xyz(2, 0) = 0.01911
xyz(3, 0) = 0.084736
xyz(4, 0) = 0.204492
xyz(5, 0) = 0.314679
xyz(6, 0) = 0.383734
xyz(7, 0) = 0.370702
xyz(8, 0) = 0.302273
xyz(9, 0) = 0.195618
xyz(10, 0) = 0.080507
xyz(11, 0) = 0.016172
xyz(12, 0) = 0.003816
xyz(13, 0) = 0.037465
xyz(14, 0) = 0.117749
xyz(15, 0) = 0.236491
xyz(16, 0) = 0.376772
xyz(17, 0) = 0.529826
xyz(18, 0) = 0.705224
xyz(19, 0) = 0.878655
xyz(20, 0) = 1.01416
xyz(21, 0) = 1.11852
xyz(22, 0) = 1.12399
xyz(23, 0) = 1.03048
xyz(24, 0) = 0.856297
xyz(25, 0) = 0.647467
xyz(26, 0) = 0.431567
xyz(27, 0) = 0.268329
xyz(28, 0) = 0.152568
xyz(29, 0) = 0.081261
xyz(30, 0) = 0.040851
xyz(31, 0) = 0.019941
xyz(32, 0) = 0.009577
xyz(0, 1) = 0.000017
xyz(1, 1) = 0.000253
xyz(2, 1) = 0.002004
xyz(3, 1) = 0.008756
xyz(4, 1) = 0.021391
xyz(5, 1) = 0.038676
xyz(6, 1) = 0.062077
xyz(7, 1) = 0.089456
xyz(8, 1) = 0.128201
xyz(9, 1) = 0.18519
xyz(10, 1) = 0.253589
xyz(11, 1) = 0.339133
xyz(12, 1) = 0.460777
xyz(13, 1) = 0.606741
xyz(14, 1) = 0.761757
xyz(15, 1) = 0.875211
xyz(16, 1) = 0.961988
xyz(17, 1) = 0.991761
xyz(18, 1) = 0.99734
xyz(19, 1) = 0.955552
xyz(20, 1) = 0.868934
xyz(21, 1) = 0.777405
xyz(22, 1) = 0.658341
xyz(23, 1) = 0.527963
xyz(24, 1) = 0.398057
xyz(25, 1) = 0.283493
xyz(26, 1) = 0.179828
xyz(27, 1) = 0.107633
xyz(28, 1) = 0.060281
xyz(29, 1) = 0.0318
xyz(30, 1) = 0.015905
xyz(31, 1) = 0.007749
xyz(32, 1) = 0.003718
xyz(0, 2) = 0.000705
xyz(1, 2) = 0.010482
xyz(2, 2) = 0.086011
xyz(3, 2) = 0.389366
xyz(4, 2) = 0.972542
xyz(5, 2) = 1.55348
xyz(6, 2) = 1.96728
xyz(7, 2) = 1.9948
xyz(8, 2) = 1.74537
xyz(9, 2) = 1.31756
xyz(10, 2) = 0.772125
xyz(11, 2) = 0.415254
xyz(12, 2) = 0.218502
xyz(13, 2) = 0.112044
xyz(14, 2) = 0.060709
xyz(15, 2) = 0.030451
xyz(16, 2) = 0.013676
xyz(17, 2) = 0.003988
xyz(18, 2) = 0
xyz(19, 2) = 0
xyz(20, 2) = 0
xyz(21, 2) = 0
xyz(22, 2) = 0
xyz(23, 2) = 0
xyz(24, 2) = 0
xyz(25, 2) = 0
xyz(26, 2) = 0
xyz(27, 2) = 0
xyz(28, 2) = 0
xyz(29, 2) = 0
xyz(30, 2) = 0
xyz(31, 2) = 0
xyz(32, 2) = 0
xyz_func_10deg = xyz
End Function
Function trans2XYZ(trans As range)
Dim result(2, 0) As Double
i = 0
teisuk = 0
For Each i In trans
teisuk = teisuk + D65spectrum(j) * xyz_func_10deg(j, 1)
j = j + 1
Next
teisuk = 100# / teisuk
i = 0
j = 0
X = 0
For Each i In trans
X = X + i * 0.01 * D65spectrum(j) * xyz_func_10deg(j, 0)
j = j + 1
Next
X = teisuk * X
i = 0
j = 0
Y = 0
For Each i In trans
Y = Y + i * 0.01 * D65spectrum(j) * xyz_func_10deg(j, 1)
j = j + 1
Next
Y = teisuk * Y
i = 0
j = 0
Z = 0
For Each i In trans
Z = Z + i * 0.01 * D65spectrum(j) * xyz_func_10deg(j, 2)
j = j + 1
Next
Z = teisuk * Z
result(0, 0) = X
result(1, 0) = Y
result(2, 0) = Z
trans2XYZ = result
End Function
Function XYZ2Lab(X, Y, Z)
Dim result(2, 0) As Double
Xn = D65wp(0)
Yn = D65wp(1)
Zn = D65wp(2)
X_Xn = X / Xn
Y_Yn = Y / Yn
Z_Zn = Z / Zn
If X_Xn > ((6# / 29#) ^ 3#) Then
fx = X_Xn ^ (1# / 3#)
Else
fx = ((29# / 3#) ^ 3# * X_Xn + 16#) / 116#
End If
If Y_Yn > ((6# / 29#) ^ 3#) Then
fy = Y_Yn ^ (1# / 3#)
Else
fy = ((29# / 3#) ^ 3# * Y_Yn + 16#) / 116#
End If
If Z_Zn > ((6# / 29#) ^ 3#) Then
fz = Z_Zn ^ (1# / 3#)
Else
fz = ((29# / 3#) ^ 3# * Z_Zn + 16#) / 116#
End If
L = 116# * fy - 16#
a = 500# * (fx - fy)
b = 200# * (fy - fz)
result(0, 0) = L
result(1, 0) = a
result(2, 0) = b
XYZ2Lab = result
End Function
Function trans2lab(trans As range)
Dim result(2, 0) As Double
X = trans2XYZ(trans)(0, 0)
Y = trans2XYZ(trans)(1, 0)
Z = trans2XYZ(trans)(2, 0)
result(0, 0) = XYZ2Lab(X, Y, Z)(0, 0)
result(1, 0) = XYZ2Lab(X, Y, Z)(1, 0)
result(2, 0) = XYZ2Lab(X, Y, Z)(2, 0)
trans2lab = result
End Function
Function deltaEab(L1, a1, b1, L2, a2, b2)
deltaEab = Sqr((L2 - L1) ^ 2 + (a2 - a1) ^ 2 + (b2 - b1) ^ 2)
End Function
Function deltaE00(L1, a1, b1, L2, a2, b2)
kL = 1#
kC = 1#
kH = 1#
dLp = L2 - L1
Lbar = (L1 + L2) / 2#
Cs_1 = Sqr(a1 ^ 2# + b1 ^ 2#)
Cs_2 = Sqr(a2 ^ 2# + b2 ^ 2#)
Cbar = (Cs_1 + Cs_2) / 2#
ap_1 = a1 + (a1 / 2#) * (1# - Sqr(Cbar ^ 7# / (Cbar ^ 7# + 25# ^ 7#)))
ap_2 = a2 + (a2 / 2#) * (1# - Sqr(Cbar ^ 7# / (Cbar ^ 7# + 25# ^ 7#)))
Cp_1 = Sqr(ap_1 ^ 2# + b1 ^ 2#)
Cp_2 = Sqr(ap_2 ^ 2# + b2 ^ 2#)
Cpbar = (Cp_1 + Cp_2) / 2#
dCp = Cp_2 - Cp_1
hp_1 = WorksheetFunction.Degrees(WorksheetFunction.Atan2(ap_1, b1))
hp_2 = WorksheetFunction.Degrees(WorksheetFunction.Atan2(ap_2, b2))
If ((Cp_1 = 0#) Or (Cp_2 = 0#)) Then
dsmallhp = 0#
ElseIf (Abs(hp_1 - hp_2) <= 180#) Then
dsmallhp = hp_2 - hp_1
ElseIf ((Abs(hp_1 - hp_2) > 180#) And (hp_2 <= hp_1)) Then
dsmallhp = hp_2 - hp_1 + 360#
Else
dsmallhp = hp_2 - hp_1 - 360#
End If
dHp = 2# * Sqr(Cp_1 * Cp_2) * Sin(WorksheetFunction.Radians(dsmallhp / 2#))
If ((Cp_1 = 0#) Or (Cp_2 = 0#)) Then
Hpbar = hp_1 + hp_2
ElseIf (Abs(hp_1 - hp_2) <= 180#) Then
Hpbar = (hp_1 + hp_2) / 2#
ElseIf ((Abs(hp_1 - hp_2) > 180#) And (hp_1 + hp_2 < 360)) Then
Hpbar = (hp_1 + hp_2 + 360#) / 2#
Else
Hpbar = (hp_1 + hp_2 - 360#) / 2#
End If
T = 1# - 0.17 * Cos(WorksheetFunction.Radians(Hpbar - 30#)) _
+ 0.24 * Cos(WorksheetFunction.Radians(2# * Hpbar)) _
+ 0.32 * Cos(WorksheetFunction.Radians(3# * Hpbar + 6#)) _
- 0.2 * Cos(WorksheetFunction.Radians(4# * Hpbar - 63#))
S_L = 1# + (0.015 * (Lbar - 50#) ^ 2#) / (Sqr(20# + (Lbar - 50#) ^ 2#))
S_C = 1# + 0.045 * Cpbar
S_H = 1# + 0.015 * Cpbar * T
R_C = 2# * Sqr(Cpbar ^ 7# / (Cpbar ^ 7# + 25# ^ 7#))
dtheta = 30# * Exp(-1# * ((Hpbar - 275#) / 25#) ^ 2#)
R_T = -1# * Sin(WorksheetFunction.Radians(2 * dtheta)) * R_C
term1 = (dLp / (kL * S_L)) ^ 2#
term2 = (dCp / (kC * S_C)) ^ 2#
term3 = (dHp / (kH * S_H)) ^ 2#
term4 = R_T * (dCp / (kC * S_C)) * (dHp / (kH * S_H))
deltaE00 = Sqr(term1 + term2 + term3 + term4)
End Function
Function lab2XYZ(L, a, b)
Dim result(2, 0) As Double
Xn = D65wp(0)
Yn = D65wp(1)
Zn = D65wp(2)
fy = (L + 16#) / 116#
fx = fy + (a / 500#)
fz = fy - (b / 200#)
If fx > (6# / 29#) Then
X = fx ^ 3# * Xn
Else
X = (3# / 29#) ^ 3# * (116# * fx - 16#) * Xn
End If
If fy > (6# / 29#) Then
Y = fy ^ 3# * Yn
Else
Y = (3# / 29#) ^ 3# * (116# * fy - 16#) * Yn
End If
If fz > (6# / 29#) Then
Z = fz ^ 3# * Zn
Else
Z = (3# / 29#) ^ 3# * (116# * fz - 16#) * Zn
End If
result(0, 0) = X
result(1, 0) = Y
result(2, 0) = Z
lab2XYZ = result
End Function
Function XYZ2RGB(X, Y, Z)
Dim result(2, 0) As Double
M_inv00 = 3.24097
M_inv01 = -1.537383
M_inv02 = -0.498611
M_inv10 = -0.969244
M_inv11 = 1.875968
M_inv12 = 0.041555
M_inv20 = 0.05563
M_inv21 = -0.203977
M_inv22 = 1.056972
R_tmp = M_inv00 * X / 100# + M_inv01 * Y / 100# + M_inv02 * Z / 100#
G_tmp = M_inv10 * X / 100# + M_inv11 * Y / 100# + M_inv12 * Z / 100#
B_tmp = M_inv20 * X / 100# + M_inv21 * Y / 100# + M_inv22 * Z / 100#
If R_tmp <= 0.0031308 Then
r = (12.92 * R_tmp) * 255#
Else
r = (1.055 * R_tmp ^ (1# / 2.4) - 0.055) * 255#
End If
If G_tmp <= 0.0031308 Then
g = (12.92 * G_tmp) * 255#
Else
g = (1.055 * G_tmp ^ (1# / 2.4) - 0.055) * 255#
End If
If B_tmp <= 0.0031308 Then
b = (12.92 * B_tmp) * 255#
Else
b = (1.055 * B_tmp ^ (1# / 2.4) - 0.055) * 255#
End If
result(0, 0) = r
result(1, 0) = g
result(2, 0) = b
XYZ2RGB = result
End Function
Sub Coloring_Cell()
i = 9
Do While Cells(4, i) <> ""
L = Cells(4, i)
a = Cells(5, i)
b = Cells(6, i)
X = lab2XYZ(L, a, b)(0, 0)
Y = lab2XYZ(L, a, b)(1, 0)
Z = lab2XYZ(L, a, b)(2, 0)
r = XYZ2RGB(X, Y, Z)(0, 0)
g = XYZ2RGB(X, Y, Z)(1, 0)
b = XYZ2RGB(X, Y, Z)(2, 0)
Cells(9, i).Interior.Color = RGB(r, g, b)
i = i + 1
Loop
End Sub