Process color distances to find max color distances between all different colors.

colorclusters.py

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

# 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" % (r, g, 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 nearest_pair(colors):
    min_distance = float('inf')
    min_pair = None
    for j, q1 in enumerate(colors):
        for k, q2 in enumerate(colors):
            if j == k:
                continue
            current_distance = color_distance(*q1, *q2)
            if current_distance < min_distance:
                min_distance = current_distance
                min_pair = (j,k)
    return min_pair


colors = [[0, 0, 0]]


skip = 1  # Setting this higher makes it go faster, but less brute force.
print("Phase 1: Find our colors.")

for values in range(4):
    max_distance = -float("inf")
    max_value = None
    colors.append([0, 0, 0])
    for r in range(0, 256, skip):
        colors[-1][0] = r
        for g in range(0, 256, skip):
            colors[-1][1] = g
            for b in range(0, 256, skip):
                colors[-1][2] = b
                current_distance = maximized_distance(colors)
                if current_distance > max_distance:
                    max_distance = current_distance
                    max_value = [r, g, b]
    print(max_distance)
    if max_value is not None:
        colors[-1] = max_value
    print_colors(colors)
    # print(colors)

print("Phase 1 complete.")
print_colors(colors)

print("Phase 2 adjust our colors")


def attempt_improvement_brute(colors, level):
    distance = maximized_distance(colors)
    attempt = deepcopy(colors)
    for manipulation in range(3 ** (3 * len(attempt))):
        for b in range(3):
            for c in range(len(attempt)):
                change = manipulation % 3
                manipulation //= 3
                if change == 1 and attempt[c][b] != 0:
                    attempt[c][b] -= 1
                if change == 2 and attempt[c][b] != 255:
                    attempt[c][b] += 1
        attempt_distance = maximized_distance(attempt)
        if attempt_distance > distance:
            print_colors(attempt)
            attempt = deepcopy(attempt)
            return attempt
    return None


def attempt_improvement_pair(colors, level):
    improved = False
    distance = maximized_distance(colors)
    attempt = deepcopy(colors)
    best = None
    for manipulation in range(3 ** 6):
        pair = nearest_pair(attempt)
        j, k = pair
        for b in range(3):
            change = manipulation % 3
            manipulation //= 3
            if change == 1 and attempt[j][b] != 0:
                attempt[j][b] -= 1
            if change == 2 and attempt[j][b] != 255:
                attempt[j][b] += 1
        for b in range(3):
            change = manipulation % 3
            manipulation //= 3
            if change == 1 and attempt[k][b] != 0:
                attempt[k][b] -= 1
            if change == 2 and attempt[k][b] != 255:
                attempt[k][b] += 1
        attempt_distance = maximized_distance(attempt)
        if attempt_distance > distance:
            print_colors(attempt)
            distance = attempt_distance
            attempt = deepcopy(attempt)
            best = deepcopy(attempt)
            improved = True
    if improved:
        return best
    return None


level = 1
while True:
    print_colors(colors)
    c = attempt_improvement_pair(colors, level)
    if c is None:
        c = attempt_improvement_brute(colors, level)
        if c is None:
            break
    print_colors(c)
    colors = c
    print("Phase2-%d Improved." % level)
    print_colors(colors)
    level += 1

print("Phase 3 making image.")
print_colors(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")

A naive early result.

Code! Superb! I'm almost certainly going to use this in the future.

You almost did the entire work of my immediate needs for me — these are backgrounds for black text so the near-black is out. In fact the original first 6 "furthest distance from previous" colors seem so perfect that it's pure curiosity on my part whether such excellence could be exceeded, because my brain says no, but the math says... theoretically yes? But maybe at such a small sample the natural imperfection of any perceptual approximation, like CIEDE2000, makes the results a wash.

I suppose for the equidistant you would compute 6 with one of them locked in as #000000 as opposed to discarding the darkest? Hmm... 🤔

Forgive me I haven't poured over the code yet but it sounds like an optimization might be simulated annealing?
Not that we need to bother. 😅