Calculate color lists for various numbers of equidistant colors.

calc_colors.py

import math
from copy import deepcopy
from functools import lru_cache
from PIL import ImageDraw, Image

from random import Random
r = Random()


# Color conversion formula borrowed from:
# http://www.easyrgb.com/index.php?X=MATH&H=02#text2

ref_X = 95.047  # ref_X =  95.047   Observer= 2°, Illuminant= D65
ref_Y = 100.000  # ref_Y = 100.000
ref_Z = 108.883  # ref_Z = 108.883
TwentyFivePowerSeven = 25 * 25 * 25 * 25 * 25 * 25 * 25


@lru_cache(maxsize=0xFFF)
def convertRGBLab(var_R: float, var_G: float, var_B: float):
    """
    Convert RGB to Lab colors.
    :param var_R:
    :param var_G:
    :param var_B:
    :return:
    """
    x, y, z = convertRGBXYZ(var_R, var_G, var_B)
    return convertXYZLab(x, y, z)


def convertRGBXYZ(var_R: float, var_G: float, var_B: float):
    """
    Convert RGB to XYZ colors

    :param var_R:
    :param var_G:
    :param var_B:
    :return:
    """
    if var_R > 0.04045:
        var_R = math.pow(((var_R + 0.055) / 1.055), 2.4)

    else:
        var_R = var_R / 12.92
    if var_G > 0.04045:
        var_G = math.pow(((var_G + 0.055) / 1.055), 2.4)

    else:
        var_G = var_G / 12.92
    if var_B > 0.04045:
        var_B = math.pow(((var_B + 0.055) / 1.055), 2.4)
    else:
        var_B = var_B / 12.92
    var_R = var_R * 100
    var_G = var_G * 100
    var_B = var_B * 100  # Observer. = 2°, Illuminant = D65
    X = var_R * 0.4124 + var_G * 0.3576 + var_B * 0.1805
    Y = var_R * 0.2126 + var_G * 0.7152 + var_B * 0.0722
    Z = var_R * 0.0193 + var_G * 0.1192 + var_B * 0.9505
    return X, Y, Z


def convertXYZLab(X: float, Y: float, Z: float):
    """
    Convert xyz to lab colors.
    :param X:
    :param Y:
    :param Z:
    :return:
    """
    var_X = X / ref_X
    var_Y = Y / ref_Y
    var_Z = Z / ref_Z

    if var_X > 0.008856:
        var_X = math.pow(var_X, 1.0 / 3.0)
    else:
        var_X = (7.787 * var_X) + (16.0 / 116.0)
    if var_Y > 0.008856:
        var_Y = math.pow(var_Y, 1.0 / 3.0)
    else:
        var_Y = (7.787 * var_Y) + (16.0 / 116.0)
    if var_Z > 0.008856:
        var_Z = math.pow(var_Z, 1.0 / 3.0)
    else:
        var_Z = (7.787 * var_Z) + (16.0 / 116.0)
    CIE_L = (116 * var_Y) - 16
    CIE_a = 500 * (var_X - var_Y)
    CIE_b = 200 * (var_Y - var_Z)
    return CIE_L, CIE_a, CIE_b


def polarAngle(var_a: float, var_b: float) -> float:  # Function returns CIE-H° value
    return math.degrees(math.atan2(var_a, var_b))


def DeltaE00(
    LabL1: float, Laba1: float, Labb1: float, LabL2: float, Laba2: float, Labb2: float
) -> float:
    """
    Apply the CIEDeltaE00 color distance algorithm.

    :param LabL1:
    :param Laba1:
    :param Labb1:
    :param LabL2:
    :param Laba2:
    :param Labb2:
    :return:
    """
    # CIE-L*1, CIE-a*1, CIE-b*1          //Color #1 CIE-L*ab values
    # CIE-L*2, CIE-a*2, CIE-b*2          //Color #2 CIE-L*ab values
    WHTL = 1.0
    WHTC = 1.0
    WHTH = 1.0  # Weight factor

    xC1 = math.sqrt(Laba1 * Laba1 + Labb1 * Labb1)
    xC2 = math.sqrt(Laba2 * Laba2 + Labb2 * Labb2)

    xCX = (xC1 + xC2) / 2
    xCXSeven = xCX * xCX * xCX * xCX * xCX * xCX * xCX

    xGX = 0.5 * (1 - math.sqrt((xCXSeven) / ((xCXSeven) + (TwentyFivePowerSeven))))
    xNN = (1 + xGX) * Laba1
    xC1 = math.sqrt(xNN * xNN + Labb1 * Labb1)
    xH1 = polarAngle(xNN, Labb1)
    xNN = (1 + xGX) * Laba2
    xC2 = math.sqrt(xNN * xNN + Labb2 * Labb2)
    xH2 = polarAngle(xNN, Labb2)
    xDL = LabL2 - LabL1
    xDC = xC2 - xC1
    xDH = 0.0
    if (xC1 * xC2) == 0:
        xDH = 0
    else:
        xNN = xH2 - xH1  # round(xH2 - xH1, 12);
        if abs(xNN) <= 180:
            xDH = xH2 - xH1
        else:
            if xNN > 180:
                xDH = xH2 - xH1 - 360
            else:
                xDH = xH2 - xH1 + 360
    xDH = 2 * math.sqrt(xC1 * xC2) * math.sin(math.radians(xDH / 2.0))
    xLX = (LabL1 + LabL2) / 2.0
    xCY = (xC1 + xC2) / 2.0
    xHX = 0.0
    if (xC1 * xC2) == 0:
        xHX = xH1 + xH2
    else:
        xNN = abs(xH1 - xH2)  # Math.round(xH1 - xH2, 12)
        if xNN > 180:
            if (xH2 + xH1) < 360:
                xHX = xH1 + xH2 + 360
            else:
                xHX = xH1 + xH2 - 360
        else:
            xHX = xH1 + xH2
        xHX /= 2.0
    xTX = (
        1
        - 0.17 * math.cos(math.radians(xHX - 30))
        + 0.24 * math.cos(math.radians(2 * xHX))
        + 0.32 * math.cos(math.radians(3 * xHX + 6))
        - 0.20 * math.cos(math.radians(4 * xHX - 63))
    )
    xPH = 30 * math.exp(-((xHX - 275) / 25) * ((xHX - 275) / 25))
    xCYSeven = xCY * xCY * xCY * xCY * xCY * xCY * xCY
    xRC = 2 * math.sqrt((xCYSeven) / ((xCYSeven) + (TwentyFivePowerSeven)))
    xSL = 1 + (
        (0.015 * ((xLX - 50) * (xLX - 50))) / math.sqrt(20 + ((xLX - 50) * (xLX - 50)))
    )
    xSC = 1 + 0.045 * xCY
    xSH = 1 + 0.015 * xCY * xTX
    xRT = -math.sin(math.radians(2 * xPH)) * xRC
    xDL = xDL / (WHTL * xSL)
    xDC = xDC / (WHTC * xSC)
    xDH = xDH / (WHTH * xSH)
    DeltaE00 = math.sqrt(xDL * xDL + xDC * xDC + xDH * xDH + xRT * xDC * xDH)
    return DeltaE00


@lru_cache(maxsize=0xFFF)
def color_distance(r1, g1, b1, r2, g2, b2):
    """
    Find our color distance, using CIEDeltaE00

    :param r1:
    :param g1:
    :param b1:
    :param r2:
    :param g2:
    :param b2:
    :return:
    """
    Lab1 = convertRGBLab(r1 / 255.0, g1 / 255.0, b1 / 255.0)
    Lab2 = convertRGBLab(r2 / 255.0, g2 / 255.0, b2 / 255.0)
    return DeltaE00(*Lab1, *Lab2)


def hexcolor(r, g, b):
    return "#%02X%02X%02X" % (int(round(r)), int(round(g)), int(round(b)))


def hexcolors(colors):
    c = []
    for color in colors:
        c.append(hexcolor(*color))
    return c


def print_colors(colors):
    dist = maximized_distance(colors)
    print(str(dist) + ": " + " ".join(hexcolors(colors)))


def maximized_distance(colors):
    distances = []
    for j, q1 in enumerate(colors):
        for k, q2 in enumerate(colors):
            if j == k:
                continue
            distances.append(color_distance(*q1, *q2))
    if len(distances):
        return min(distances)
    else:
        return -float("inf")


def print_distance(colors):
    distances = []
    for j, q1 in enumerate(colors):
        for k, q2 in enumerate(colors):
            if j == k:
                continue
            distances.append(color_distance(*q1, *q2))
    print(distances)


def swatch(colors):
    img = Image.new("RGB", (len(colors) * 20, 20))
    draw = ImageDraw.Draw(img)
    hcolors = hexcolors(colors)
    for i, c in enumerate(hcolors):
        draw.rectangle((i * 20, 0, (i + 1) * 20, 20), fill=c)
    img.save("swatch.png")



print("Phase 1: Find our colors.")

# CHANGE DATA HERE
num_of_colors = 10

random = Random()
colors = []
for i in range(num_of_colors):
    colors.append([0.0, 0.0, 0.0])

for cc in colors:
    for i in range(len(cc)):
        cc[i] = random.randint(0, 256)

max_value = deepcopy(colors)
max_distance = -float("inf")


while True:
    step = random.random() * 10
    for j, cc in enumerate(max_value):
        for k in range(len(cc)):
            m = random.random() * step
            c = cc[k]
            c += m - (step / 2.0)
            if c > 255:
                c = 255
            if c < 0:
                c = 0
            colors[j][k] = c
    current_distance = maximized_distance(colors)
    if current_distance > max_distance:
        print((current_distance, max_distance))
        print(colors)
        max_distance = current_distance
        max_value = deepcopy(colors)
        print_colors(colors)
        swatch(max_value)