low_discrepancy.py

import numpy

gamma = 2.2
dark_cutoff_fraction = 0.2
non_grey = 0.4

# Blues don't appear all that bright
blue_dist_multiplier = 0.4

# Raw and gamma-corrected colors
r, g, b = numpy.mgrid[:256, :256, :256]
r_gamma = (r / 255) ** gamma
g_gamma = (g / 255) ** gamma
b_gamma = (b / 255) ** gamma

distances = numpy.full(shape=(256, 256, 256), fill_value=float("inf"))

# Break ties in this random order
resolve_priority = numpy.random.random(distances.shape)

# Remove, or penalize, some unwanted colors
resolve_mask = numpy.ones_like(resolve_priority)

# Skip darks
sum_rgb = r_gamma + g_gamma + b_gamma
dark_cutoff = numpy.percentile(sum_rgb, 100 * dark_cutoff_fraction)
resolve_mask[sum_rgb <= dark_cutoff] = 0

# Discourage greys
max_rgb = numpy.maximum.reduce([r_gamma, g_gamma, b_gamma])
min_rgb = numpy.minimum.reduce([r_gamma, g_gamma, b_gamma])
resolve_mask[max_rgb - min_rgb < non_grey] = 0.5

def calc_distances(new_r, new_g, new_b):
    return (
        (r_gamma - (new_r / 255)**gamma) ** 2 +
        (g_gamma - (new_g / 255)**gamma) ** 2 +
        ((b_gamma - (new_b / 255)**gamma) * blue_dist_multiplier) ** 2
    ) ** 0.5

print('<svg width="512" height="512" xmlns="http://www.w3.org/2000/svg" xmlns:svg="http://www.w3.org/2000/svg">')

# Start on grey
new_r, new_g, new_b = 127, 127, 127

for y in range(16):
    for x in range(16):
        print(
            f'<rect '
            f' width="32" height="32" x="{x * 32:3}" y="{y * 32:3}"'
            f' style="fill:rgb({new_r:3},{new_g:3},{new_b:3})"'
            f' />'
        )

        new_distances = calc_distances(new_r, new_g, new_b)
        distances = numpy.minimum(distances, new_distances)
        distances *= resolve_mask
        new_rgb = ((distances == distances.max()) * resolve_priority).argmax()
        new_r, new_g, new_b = numpy.unravel_index(new_rgb, shape=distances.shape)

print('</svg>')

output.svg