color difference

VB source code for color difference calculation.

Public Function DeltaOfColors(ByVal mode As Integer, ByVal r1 As Integer, ByVal g1 As Integer, ByVal b1 As Integer, ByVal r2 As Integer, ByVal g2 As Integer, ByVal b2 As Integer) As Double

Dim LAB1 As HSB

Dim LAB2 As HSB

Dim rez As Double

LAB1 = RGBtoLAB(r1, g1, b1)

LAB2 = RGBtoLAB(r2, g2, b2)

rez = ColorDelta(mode, LAB1.Hue, LAB1.Saturation, LAB1.Brightness, LAB2.Hue, LAB2.Saturation, LAB2.Brightness)

Return rez

End Function

Public Function ColorDelta(ByVal mode As Integer, ByVal L1 As Double, ByVal a1 As Double, ByVal b1 As Double, ByVal L2 As Double, ByVal a2 As Double, ByVal b2 As Double) As Double

Dim rez, xDE As Double

Dim xC1, xC2, xDL, xDC, xDH, xSC, xSH, xCX, xGX, xNN, xH1, xH2 As Double

Dim xLX, xCY, xHX, xTX, xPH, xRC, xSL, xRT, xff, xTT As Double

Dim C1 As Double

Dim C2 As Double

Dim CIE_1_a_squared As Double '= a1 * a1

Dim CIE_1_b_squared As Double '= b1 * b1

Dim CIE_2_a_squared As Double '= a2 * a2

Dim CIE_2_b_squared As Double '= b2 * b2

Dim delta_a As Double

Dim delta_a_squared As Double

Dim delta_b As Double

Dim delta_b_squared As Double

Dim delta_C_ab As Double

Dim delta_C_ab_divisor As Double

Dim delta_C_ab_squared As Double

Dim delta_E_Lab As Double

Dim delta_H_ab As Double

Dim delta_H_ab_divisor As Double

Dim delta_L As Double

Dim delta_L_squared As Double

Dim K_1 As Double

Dim K_2 As Double

Dim dummi As Double

Dim DC As Double

Dim DH As Double

Dim DL As Double

Dim Da As Double

Dim Db As Double

Dim SL As Double

Dim SC As Double

Dim SH As Double

Dim T As Double

Dim F As Double

Dim H1 As Double

'Dim t As Double

Dim weighting_factor_C As Double = 1

Dim weighting_factor_H As Double = 1

Dim weighting_factor_L As Double = 1

Select Case mode

Case 0 ', 7 'delta_E_2000

Dim c As Double = Math.Pow(25, 7)

CIE_1_a_squared = a1 * a1

CIE_1_b_squared = b1 * b1

CIE_2_a_squared = a2 * a2

CIE_2_b_squared = b2 * b2

xC1 = Math.Sqrt((CIE_1_a_squared + CIE_1_b_squared))

xC2 = Math.Sqrt((CIE_2_a_squared + CIE_2_b_squared))

xCX = ((xC1 + xC2) / 2.0)

T = Math.Pow(xCX, 7)

xGX = (0.5 * (1 - Math.Sqrt(T / (T + c))))

xNN = ((1 + xGX) * a1)

xC1 = Math.Sqrt((xNN * xNN + CIE_1_b_squared))

xH1 = CieLab2Hue(xNN, b1)

xNN = ((1 + xGX) * a2)

xC2 = Math.Sqrt(((xNN * xNN) + CIE_2_b_squared))

xH2 = CieLab2Hue(xNN, b2)

xDL = (L2 - L1)

xDC = (xC2 - xC1)

If (xC1 * xC2) = 0 Then

xDH = 0.0

Else

T = xH2 - xH1

xNN = Math.Round(T, 12)

If Math.Abs(xNN) <= 180 Then

xDH = T

Else

If xNN > 180 Then

xDH = T - 360.0

Else

xDH = T + 360.0

End If

xDH = 2.0 * Math.Sqrt(xC1 * xC2) * Math.Sin(deg2rad(xDH / 2.0))

xLX = (L1 - L2) / 2.0

xCY = (xC1 + xC2) / 2.0

T = xH1 + xH2

If (xC1 * xC2) = 0 Then

xHX = T

Else

xNN = Math.Abs(Math.Round((xH1 - xH2), 12))

If xNN > 180 Then

If T < 360.0 Then

xHX = T + 360.0

Else

xHX = T - 360.0

End If

Else

xHX = T

End If

xHX /= 2

End If

xTX = 1.0 - 0.17 * Math.Cos(deg2rad(xHX - 30.0)) + 0.24 * Math.Cos(deg2rad(2.0 * xHX)) + 0.32 * Math.Cos(deg2rad(3.0 * xHX + 6.0)) - 0.2 * Math.Cos(deg2rad(4.0 * xHX - 63.0))

T = (xHX - 275.0) / 25.0

xPH = 30.0 * Math.Exp(-(T * T))

T = Math.Pow(xCY, 7)

xRC = 2.0 * Math.Sqrt(T / (T + c))

T = xLX - 50.0

xSL = 1.0 + (0.015 * (T * T)) / Math.Sqrt(20.0 + (T * T))

xSC = 1.0 + 0.045 * xCY

xSH = 1.0 + 0.015 * xCY * xTX

xRT = -Math.Sin(deg2rad(2.0 * xPH)) * xRC

xDL /= (weighting_factor_L * xSL)

xDC /= (weighting_factor_C * xSC)

xDH /= (weighting_factor_H * xSH)

rez = Math.Sqrt((xDL * xDL) + (xDC * xDC) + (xDH * xDH) + (xRT * xDC * xDH))

Case 1 'DElta 1994

CIE_1_a_squared = a1 * a1

CIE_1_b_squared = b1 * b1

CIE_2_a_squared = a2 * a2

CIE_2_b_squared = b2 * b2

delta_L = L1 - L2

delta_L_squared = delta_L * delta_L

delta_a = a1 - a2

delta_a_squared = delta_a * delta_a

delta_b = b1 - b2

delta_b_squared = delta_b * delta_b

delta_E_Lab = Math.Sqrt(delta_L_squared + delta_a_squared + delta_b_squared)

C1 = Math.Sqrt(CIE_1_a_squared + CIE_1_b_squared)

C2 = Math.Sqrt(CIE_2_a_squared + CIE_2_b_squared)

delta_C_ab = C1 - C2

delta_C_ab_squared = delta_C_ab * delta_C_ab

If (delta_a_squared + delta_b_squared) > delta_C_ab_squared Then

delta_H_ab = Math.Sqrt(delta_a_squared + delta_b_squared - delta_C_ab_squared)

Else

delta_H_ab = 0.0

End If

K_1 = 0.045

K_2 = 0.015

delta_C_ab_divisor = 1.0 + (K_1 * C1)

delta_H_ab_divisor = 1.0 + (K_2 * C1)

delta_C_ab /= delta_C_ab_divisor

delta_H_ab /= delta_H_ab_divisor

rez = (Math.Sqrt(delta_L_squared + delta_C_ab * delta_C_ab + delta_H_ab * delta_H_ab))

Case 2 ', 6 'CIE 1976 DeltaE

rez = Math.Sqrt((L2 - L1) * (L2 - L1) + (a2 - a1) * (a2 - a1) + (b2 - b1) * (b2 - b1))

Case 3 '13 'DeltaE (CMC)

C1 = Math.Sqrt(a1 * a1 + b1 * b1)

C2 = Math.Sqrt(a2 * a2 + b2 * b2)

DC = C1 - C2

DL = L1 - L2

Da = a1 - a2

Db = b1 - b2

dummi = Da * Da + Db * Db - DC * DC

If dummi < 0 Then

DH = 0

Else

DH = Math.Sqrt(dummi)

End If

If L1 < 16 Then

SL = 0.511

Else

SL = (0.040975 * L1) / (1 + 0.01765 * L1)

End If

' H1 = 1 / Math.Tan(b1 / a1)

H1 = rad2deg(Math.Atan2(a1, b1))

If H1 < 0 Then

H1 = H1 + 360

ElseIf H1 > 360 Then

H1 = H1 - 360

End If

SC = (0.0638 * C1) / (1 + 0.0131 * C1) + 0.638

If 164 <= H1 <= 345 Then

T = 0.56 + Math.Abs(0.2 * Math.Cos(H1 + 168.0))

Else

T = 0.36 + Math.Abs(0.4 * Math.Cos(H1 + 35.0))

End If

F = Math.Sqrt(C1 * C1 * C1 * C1 / (C1 * C1 * C1 * C1 + 1900.0))

SH = SC * (F * T + 1 - F)

'rez = Math.Sqrt((DL / SL) * (DL / SL) + (DC / SC) * (DC / SC) + (DH / SH) * (DH / SH))

rez = Math.Sqrt((DL / 2.0 * SL) * (DL / 2.0 * SL) + (DC / SC) * (DC / SC) + (DH / SH) * (DH / SH))

Case 4 'delta_C

rez = Math.Abs(Math.Sqrt(a2 * a2 + b2 * b2) - Math.Sqrt(a1 * a1 + b1 * b1))

Case 5 '6 'delta_H

xDE = Math.Sqrt(a2 * a2 + b2 * b2) - Math.Sqrt(a1 * a1 + b1 * b1)

' Try

dummi = (a2 - a1) * (a2 - a1) + (b2 - b1) * (b2 - b1) - xDE * xDE

If dummi < 0 Then

rez = 0

Else

rez = Math.Sqrt(dummi)

End If

Case 6 'DElta 1994

DL = L1

Da = a1

Db = b1

L1 = L2

a1 = a2

b1 = b2

L2 = DL

a2 = Da

b2 = Db

CIE_1_a_squared = a1 * a1

CIE_1_b_squared = b1 * b1

CIE_2_a_squared = a2 * a2

CIE_2_b_squared = b2 * b2

delta_L = L1 - L2

delta_L_squared = delta_L * delta_L

delta_a = a1 - a2

delta_a_squared = delta_a * delta_a

delta_b = b1 - b2

delta_b_squared = delta_b * delta_b

delta_E_Lab = Math.Sqrt(delta_L_squared + delta_a_squared + delta_b_squared)

C1 = Math.Sqrt(CIE_1_a_squared + CIE_1_b_squared)

C2 = Math.Sqrt(CIE_2_a_squared + CIE_2_b_squared)

delta_C_ab = C1 - C2

delta_C_ab_squared = delta_C_ab * delta_C_ab

If (delta_a_squared + delta_b_squared) > delta_C_ab_squared Then

delta_H_ab = Math.Sqrt(delta_a_squared + delta_b_squared - delta_C_ab_squared)

Else

delta_H_ab = 0.0

End If

K_1 = 0.045

K_2 = 0.015

delta_C_ab_divisor = 1.0 + (K_1 * C1)

delta_H_ab_divisor = 1.0 + (K_2 * C1)

delta_C_ab /= delta_C_ab_divisor

delta_H_ab /= delta_H_ab_divisor

rez = (Math.Sqrt(delta_L_squared + delta_C_ab * delta_C_ab + delta_H_ab * delta_H_ab))

Case 7 'delta_E_2000

DL = L1

Da = a1

Db = b1

L1 = L2

a1 = a2

b1 = b2

L2 = DL

a2 = Da

b2 = Db

Dim c As Double = Math.Pow(25, 7)

CIE_1_a_squared = a1 * a1

CIE_1_b_squared = b1 * b1

CIE_2_a_squared = a2 * a2

CIE_2_b_squared = b2 * b2

xC1 = Math.Sqrt((CIE_1_a_squared + CIE_1_b_squared))

xC2 = Math.Sqrt((CIE_2_a_squared + CIE_2_b_squared))

xCX = ((xC1 + xC2) / 2.0)

T = Math.Pow(xCX, 7)

xGX = (0.5 * (1 - Math.Sqrt(T / (T + c))))

xNN = ((1 + xGX) * a1)

xC1 = Math.Sqrt((xNN * xNN + CIE_1_b_squared))

xH1 = CieLab2Hue(xNN, b1)

xNN = ((1 + xGX) * a2)

xC2 = Math.Sqrt(((xNN * xNN) + CIE_2_b_squared))

xH2 = CieLab2Hue(xNN, b2)

xDL = (L2 - L1)

xDC = (xC2 - xC1)

If (xC1 * xC2) = 0 Then

xDH = 0.0

Else

T = xH2 - xH1

xNN = Math.Round(T, 12)

If Math.Abs(xNN) <= 180 Then

xDH = T

Else

If xNN > 180 Then

xDH = T - 360.0

Else

xDH = T + 360.0

End If

xDH = 2.0 * Math.Sqrt(xC1 * xC2) * Math.Sin(deg2rad(xDH / 2.0))

xLX = (L1 - L2) / 2.0

xCY = (xC1 + xC2) / 2.0

T = xH1 + xH2

If (xC1 * xC2) = 0 Then

xHX = T

Else

xNN = Math.Abs(Math.Round((xH1 - xH2), 12))

If xNN > 180 Then

If T < 360.0 Then

xHX = T + 360.0

Else

xHX = T - 360.0

End If

Else

xHX = T

End If

xHX /= 2

End If

xTX = 1.0 - 0.17 * Math.Cos(deg2rad(xHX - 30.0)) + 0.24 * Math.Cos(deg2rad(2.0 * xHX)) + 0.32 * Math.Cos(deg2rad(3.0 * xHX + 6.0)) - 0.2 * Math.Cos(deg2rad(4.0 * xHX - 63.0))

T = (xHX - 275.0) / 25.0

xPH = 30.0 * Math.Exp(-(T * T))

T = Math.Pow(xCY, 7)

xRC = 2.0 * Math.Sqrt(T / (T + c))

T = xLX - 50.0

xSL = 1.0 + (0.015 * (T * T)) / Math.Sqrt(20.0 + (T * T))

xSC = 1.0 + 0.045 * xCY

xSH = 1.0 + 0.015 * xCY * xTX

xRT = -Math.Sin(deg2rad(2.0 * xPH)) * xRC

xDL /= (weighting_factor_L * xSL)

xDC /= (weighting_factor_C * xSC)

xDH /= (weighting_factor_H * xSH)

rez = Math.Sqrt((xDL * xDL) + (xDC * xDC) + (xDH * xDH) + (xRT * xDC * xDH))

Case 15 '6 '6 'delta_H

C1 = Math.Sqrt(a1 * a1 + b1 * b1)

C2 = Math.Sqrt(a2 * a2 + b2 * b2)

DC = C1 - C2

' xDE = Math.Sqrt(a2 * a2 + b2 * b2) - Math.Sqrt(a1 * a1 + b1 * b1)

DH = Math.Sqrt((a2 - a1) * (a2 - a1) + (b2 - b1) * (b2 - b1) - xDE * xDE)

DL = L1 - L2

Da = a1 - a2

Db = b1 - b2

rez = Math.Sqrt(DL * DL + (DC / (1 + 0.045 * C1)) ^ 2 + (DH / (1 + 0.015 * C1)) ^ 2)

Case 46 '5 'delta_E_CMC

xC1 = Math.Sqrt((a1 ^ 2) + (b1 ^ 2))

xC2 = Math.Sqrt((a2 ^ 2) + (b2 ^ 2))

xff = Math.Sqrt((xC1 ^ 4) / ((xC1 ^ 4) + 1900.0))

xH1 = CieLab2Hue(a1, b1)

If (xH1 < 164 Or xH1 > 345) Then

xTT = 0.36 + Math.Abs(0.4 * Math.Cos(DtoR(35 + xH1)))

Else

xTT = 0.56 + Math.Abs(0.2 * Math.Cos(DtoR(168 + xH1)))

End If

If (L1 < 16) Then

xSL = 0.511

Else

xSL = (0.040975 * L1) / (1 + (0.01765 * L1))

End If

xSC = ((0.0638 * xC1) / (1 + (0.0131 * xC1))) + 0.638

xSH = ((xff * xTT) + 1 - xff) * xSC

xDH = Math.Sqrt((a2 - a1) ^ 2 + (b2 - b1) ^ 2 - (xC2 - xC1) ^ 2)

xSL = (L2 - L1) * xSL

xSC = (xC2 - xC1) * xSC

xSH = xDH / xSH

rez = Math.Sqrt(xSL ^ 2 + xSC ^ 2 + xSH ^ 2)

End Select

Return rez

End Function

Private Function CieLab2Hue(ByVal a As Double, ByVal b As Double) As Double

Dim bias As Double = 0.0

If (a >= 0.0) AndAlso (b = 0.0) Then Return 0.0

If (a < 0.0) AndAlso (b = 0.0) Then Return 180.0

If (a = 0.0) AndAlso (b > 0.0) Then Return 90.0

If (a = 0.0) AndAlso (b < 0.0) Then Return 270.0

If (a > 0.0) AndAlso (b > 0.0) Then bias = 0.0

If a < 0.0 Then bias = 180.0

If (a > 0.0) AndAlso (b < 0.0) Then bias = 360.0

Return (rad2deg(Math.Atan(b / a)) + bias)

End Function

Private Function deg2rad(ByVal deg As Double) As Double

Return deg * Math.PI / 180.0

End Function

Private Function rad2deg(ByVal rad As Double) As Double

Return rad / Math.PI * 180.0

End Function

Function DtoR(ByVal Degree As Integer) As Double

DtoR = Degree * 3.14159265358979 / 180

End Function