My take on a Numpy version of step sorting from "The incredibly challenging task of sorting colours" by Alan Zucconi (https://www.alanzucconi.com/2015/09/30/colour-sorting/)

generate_sorted_color_image.py

import itertools
import warnings

import numpy as np
import cv2
from skimage.color import deltaE_ciede2000, rgb2lab, lab2rgb, rgb2hsv
from skimage import img_as_ubyte, img_as_float32
from PIL import Image

h, w = (20, 5)

# # produce all possible colors
all_colors = [color for color in itertools.product(range(256), repeat=3)]
# # equivalent but ~8% slower on my computer
# all_colors = [(r, g, b) for r in range(256) for g in range(256) for b in range(256)]

all_colors = img_as_float32(np.expand_dims(np.array(all_colors, dtype=np.uint8), axis=0))
# # equivalent but ~5% slower on my computer
# all_colors = np.expand_dims(np.array(all_colors, dtype=np.float32) / 255, axis=0)

all_colors = cv2.cvtColor(all_colors, cv2.COLOR_RGB2Lab)[0]
# # equivalent but ~2400% slower on my computer
# all_colors = rgb2lab(all_colors)[0]

# select only sufficiently dissimilar colors
final_colors = []
while all_colors.size > 0:
    color = all_colors[np.random.choice(all_colors.shape[0])]
    final_colors.append(color)
    similarity = deltaE_ciede2000(color, all_colors)
    all_colors = all_colors[similarity > 10]

# # deltaE_ciede2000 distance:
# # <= 1.0 	Not perceptible by human eyes.
# # 1 - 2 	Perceptible through close observation.
# # 2 - 10 	Perceptible at a glance.
# # 11 - 49 	Colors are more similar than opposite
# # 100 	Colors are exact opposite

# # If you want all color instead, but beware of memory errors
# final_colors = all_colors

def transform_to_sorting_color_space(hsv_arr, hue_nb):
    arr = hsv_arr.copy()
    sat = arr[:, :, 1]
    grey_threshold = 0.15

    arr[:, :, 0] = (arr[:, :, 0] * (hue_nb - 2)).round() + 1

    # # personal remix, to regroup color that are really a shade of grey to the left
    indexes_greys = (sat < grey_threshold)
    arr[:, :, 0][indexes_greys] = 0

    #  invert the luminosity of every other segment
    indexes_odd_hue = (arr[:, :, 0] % 2 == 1)
    arr[:, :, 1:][indexes_odd_hue] = 1 - arr[:, :, 1:][indexes_odd_hue]

    arr[:, :, 1:] = arr[:, :, [2, 1]]

    return arr


final_colors = np.expand_dims(np.array(final_colors, dtype=np.float32), axis=0)
final_colors = lab2rgb(final_colors)
hsv_colors = rgb2hsv(final_colors)

# Step sorting
sorting_colors = transform_to_sorting_color_space(hsv_colors, 12)[0]
sort_indices = np.lexsort(sorting_colors[:, ::-1].T)
final_colors[0] = final_colors[0, sort_indices]

# # If you prefer Hue sort
# final_colors[0] = final_colors[0][np.lexsort(hsv_colors[0][:, ::-1].T)]

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    final_colors = img_as_ubyte(final_colors)[0]

print(f'nb_colors : {final_colors.shape[0]}')

# produce image
color_img = Image.new('RGB', (w * final_colors.shape[0], h))
final_colors = sum([[tuple(color)] * w for color in final_colors] * h, [])
color_img.putdata(final_colors)
color_img.save('./all_colors.png')