kL <- 1
kC <- 1
kH <- 1
K1 <- 0.045
K2 <- 0.015

# x and y should be two different colors in L*a*b* space
CIEDE2000 <- function(x, y) {
  L1 <- x$L
  L2 <- y$L
  a1 <- x$a
  a2 <- y$a
  b1 <- x$b
  b2 <- y$b

  delta_L <- L2 - L1
  L_bar <- mean(c(L1, L2))

  C1 <- sqrt(a1^2 + b1^2)
  C2 <- sqrt(a2^2 + b2^2)
  C_bar <- mean(C1, C2)
  delta_C <- C1 - C2

  sqrt_factor <- sqrt(C_bar^7 / (C_bar^7 + 25^7))
  a_prime1 <- a1 + a1 / 2 * (1 - sqrt_factor)
  a_prime2 <- a2 + a2 / 2 * (1 - sqrt_factor)

  C_prime1 <- sqrt(a_prime1^2 + b1^2)
  C_prime2 <- sqrt(a_prime2^2 + b2^2)
  C_prime_bar <- mean(C_prime1, C_prime2)
  delta_C_prime <- C_prime2 - C_prime1

  h_prime1 <- (atan2(b1, a_prime1) * 180 / pi) %% 360
  h_prime2 <- (atan2(b2, a_prime2) * 180 / pi) %% 360

  if (C_prime1 == 0 || C_prime2 == 0) {
    delta_h_prime <- 0
  } else if (abs(h_prime1 - h_prime2) <= 180) {
    delta_h_prime <- h_prime2 - h_prime1
  } else if (h_prime2 <= h_prime1) {
    delta_h_prime <- h_prime2 - h_prime1 + 360
  } else {
    delta_h_prime <- h_prime2 - h_prime1 - 360
  }

  delta_H_prime <- 2 * sqrt(C_prime1 * C_prime2) * sin((delta_h_prime * pi / 180) / 2)

  if (C_prime1 == 0 || C_prime2 == 0) {
    H_prime_bar <- h_prime1 + h_prime2
  } else if (abs(h_prime1 - h_prime2) <= 180) {
    H_prime_bar <- mean(h_prime1, h_prime2)
  } else if (h_prime1 + h_prime2 < 360) {
    H_prime_bar <- (h_prime1 + h_prime2 + 360) / 2
  } else {
    H_prime_bar <- (h_prime1 + h_prime2 - 360) / 2
  }

  T <- 1 - 0.17 * cos((H_prime_bar - 30) * pi / 180) +
    0.24 * cos(2 * H_prime_bar * pi / 180) +
    0.32 * cos((3 * H_prime_bar + 6) * pi / 180) -
    0.2 * cos((4 * H_prime_bar - 63) * pi / 180)

  SL <- 1 + (0.015 * (L_bar - 50)^2) / sqrt(20 + (L_bar - 50)^2)
  SC <- 1 + 0.045 * C_prime_bar
  SH <- 1 + 0.015 * C_prime_bar * T

  RT <- -2 * (C_prime_bar^7 / (C_prime_bar^7 + 25^7)) * sin(60 * exp(-((H_prime_bar - 275) / 25)^2) * pi / 180)

  term1 <- (delta_L / (kL * SL))^2
  term2 <- (delta_C_prime / (kC * SC))^2
  term3 <- (delta_H_prime / (kH * SH))^2
  term4 <- RT * delta_C_prime * delta_H_prime / (kC * SC * kH * SH)

  return(sqrt(term1 + term2 + term3 + term4))
}