de2000の計算コード

de2000を計算する式は書いてあってもコードは置いてないことが多いので

from math import sin, cos, atan2, exp, radians, degrees
import numpy as np
import cv2
import matplotlib.pyplot as plt

#input
L1 = 50
a1 = 20
b1 = 20
L2 = 50
a2 = 50
b2 = 10
k_L = 1
k_C = 1
k_H = 1

tmp1 = [[[L1/255*100,a1-126,b1-126] for _ in range(50)] for _ in range(50)]
tmp2 = [[[L2/255*100,a2-126,b2-126] for _ in range(50)] for _ in range(50)]
tmp1 = np.array(tmp1).reshape(50,50,3).astype('uint8')
tmp2 = np.array(tmp2).reshape(50,50,3).astype('uint8')
plt.imshow(cv2.cvtColor(tmp1,cv2.COLOR_Lab2RGB))
plt.show()
plt.imshow(cv2.cvtColor(tmp2,cv2.COLOR_Lab2RGB))
plt.show()

#equation
def calcDelta_h(h1,h2):
  if abs(h1-h2) <= 180:
    return h2-h1
  elif (abs(h1-h2) > 180 and h2 <= h1):
    return h2-h1+360
  elif (abs(h1-h2) > 180 and h2 > h1):
    return h2-h1-360
  else:
    return None

def calcH_bar_dash(h1,h2,c1,c2):
  if abs(h1-h2) <= 180:
    if c1 ==0 or c2 == 0:
      return (h1+h2)
    else:
      return (h1+h2)/2
  elif (abs(h1-h2) > 180 and h1 + h2 < 360):
    if c1 ==0 or c2 == 0:
      return (h1+h2+360)
    else:
      return (h1+h2+360)/2
  elif (abs(h1-h2) > 180 and h1+h2 >= 360):
    if c1 ==0 or c2 == 0:
      return (h1+h2-360)
    else:
      return (h1+h2-360)/2
  else:
    return None

def calc_h_dash(a,b):
  if a == 0 or b == 0:
    return 0
  elif atan2(b,a)<0:
    return degrees(atan2(b,a))+360
  else:
    return degrees(atan2(b,a))

delta_L_dash = L2-L1
L_bar = (L1+L2)/2

C1 = (a1**2+b1**2)**0.5
C2 = (a2**2+b2**2)**0.5
C_bar = (C1+C2)/2

a1_dash = a1+a1/2*(1-(C_bar**7/(C_bar**7+25**7))**0.5)
a2_dash = a2+a2/2*(1-(C_bar**7/(C_bar**7+25**7))**0.5)

c1_dash = (a1_dash**2+b1**2)**0.5
c2_dash = (a2_dash**2+b2**2)**0.5

c_bar_dash = (c1_dash+c2_dash)/2
delta_C_dash = c2_dash - c1_dash

h1_dash = calc_h_dash(a1_dash, b1)
h2_dash = calc_h_dash(a2_dash, b2)

if c1_dash == 0 or c2_dash == 0:
  delta_h = 0
else:
  delta_h = calcDelta_h(h1_dash,h2_dash)

delta_H_dash = 2*(c1_dash*c2_dash)**0.5 * sin(radians(delta_h/2))

H_bar_dash = calcH_bar_dash(h1_dash, h2_dash, c1_dash, c2_dash)

T = 1-0.17*cos(radians(H_bar_dash-30))+0.24*cos(radians(2*H_bar_dash))+0.32*cos(radians(3*H_bar_dash+6))-0.20*cos(radians(4*H_bar_dash-63))

S_L = 1+0.015*(L_bar-50)**2/(20+(L_bar-50)**2)**0.5
S_C = 1+0.045*c_bar_dash
S_H = 1+0.015*c_bar_dash*T

R_T = -2*(c_bar_dash**7/(c_bar_dash**7+25**7))**0.5*sin(radians(60*exp(-1*((H_bar_dash-275)/25)**2)))

delta_E00 = ((delta_L_dash/(k_L*S_L))**2+(delta_C_dash/(k_C*S_C))**2+(delta_H_dash/(k_H*S_H))**2+R_T*delta_C_dash*delta_H_dash/(k_C*S_C*k_H*S_H))**0.5

print(delta_E00)