Revisiting Yliluoma’s ordered dither algorithm

Summary: This article discusses Joel Yliluoma’s 2011 ordered dither algorithm (left), explains its internals in greater detail than other treatments, presents new simplified variants (middle), and compares the results to a state-of-the-art algorithm (right). Source code is included at the end. A 1-bit introduction Let’s lay some foundations first. I assume you know what ordered dithering looks like (if not, see above). In black and white it’s easy to code. First you acquire a threshold matrix from somewhere. For example, a 4x4 Bayer matrix like this:

\begin{bmatrix} 0 & 8 & 2 & 10 \\ 12 & 4 & 14 & 6 \\ 3 & 11 & 1 & 9 \\ 15 & 7 & 13 & 5 \end{bmatrix}

To apply the matrix to a grayscale image, you go over each pixel, find out which matrix element it corresponds to (conceptually, the matrix is tiled over the image), and read a threshold from the matrix. If the input pixel’s gray value is higher than the threshold, you output a white pixel. The code could look like this: bayer_4x4 = np.array([ [ 0, 8, 2, 10], [12, 4, 14, 6], [ 3, 11, 1, 9], [15, 7, 13, 5]])

def dither_bayer4x4_1bit(img: np.ndarray): img_float = img / 255 # process in [0,1] range H, W, _ = img.shape

# The floating point output image out = np.zeros((H, W, 1), dtype=float)

for y in range(H): for x in range(W): color = img_float[y, x] threshold = (bayer_4x4[y%4,x%4] + 0.5) / 16.0 if color > threshold: out[y, x] = 1 else: out[y, x] = 0

return out And below you see what it produces.

For 1-bit output, an ordered matrix like this creates pretty grainy results. If I had to choose, something like Atkinson error diffusion dithering would work better in black and white. But this was a mere warmup exercise. The real question is how to do the same for color images. To create something like this:

16-color indexed image with ordered dithering

Today we are discussing how exactly it’s done. Ordered dithering in color is both harder and easier than you think Doing ordered dithering well is surprisingly tricky. Getting started is easy though. We can reinterpret the threshold matrix as structured noise, repeat it over the whole image and add it to the original colors, and then find the closest color for each pixel. The 1-bit example code turns into something like this: # Assume "offset" is the average distance between two palette colors.

for y in range(H): for x in range(W): # input color color = img_float[y, x]

# (0, 1) range threshold threshold = (bayer_4x4[y%4,x%4] + 0.5) / 16.0

# bias threshold to (-0.5, 0.5) range, then add moved = color + offset * (threshold - 0.5)

# find closest palette color to 'moved' idx = find_index_of_closest(moved) inds[y, x] = idx This actually works better than it should, at least with the palettes used for my test images. The only question is the noise magnitude, the value of offset above. There’s no right answer, but in my experiments this magic formula worked OK: offset = dither_strength * 0.5 * median(pairwise_color_distances) Here’s how this “greyscale offset” method looks:

Not bad! The black dots in the sky are a bit rough but those could be cleaned up by reducing dither strength. For arbitrary palettes, this approach tends to produce desaturated results, unfortunately. But when a palette is designed for the image like here, it works OK.

N-candidate methods

One family of more advanced solutions to the color selection problem involves collecting a set of N candidate palette colors for each pixel, assigning a probability to each, and then picking one either at random or via a threshold matrix. These are called N-candidate methods. I won’t go into detail here, but will point you to a 2023 blog post titled Ordered Dithering with Arbitrary or Irregular Colour Palettes that explains everything you need to know. Please read at least its “The Probability Matrix” section. A key property these algorithms try to satisfy is local mean reproduction: when the dithered result is seen from afar (or blurred), it should look like the original. This translates to minimizing the distance between the input pixel color \mathbf{p} and the weighted sum of all chosen candidates w_1 \mathbf{r}_1 + ... + w_N \mathbf{r}_N. For more details, I suggest reading the blog post linked above. Knoll’s algorithm

One algorithm that minimizes the above distance astonishingly well is Thomas Knoll’s dither algorithm, famously used in Adobe Photoshop (he’s its inventor after all). Knoll’s algorithm first sets a goal color \mathbf{x}_1 to the palette color closest to the input color \mathbf{p}. This palette color \mathbf{r}_1 is the first candidate color. Then the algorithm measures how much error there is between \mathbf{r}_1 and \mathbf{p}, and moves the goal in the opposite direction:

\mathbf{x}_{2} = \mathbf{x}_1 + (\mathbf{p} - \mathbf{r}_1).

Now the process repeats and a new closest point to the goal \mathbf{x}_2 is found. This will be the second candidate color. It compensates for the error of the earlier candidate; if the first palette color found was too blue, then this one will be yellowish. After N rounds of this, the iteration has visited different colors around \mathbf{p}. Each color could’ve been visited multiple times, and the final candidate probabilities are relative to the frequency of the visits.

The process can be summarized as follows:

Knoll’s algorithm’s error compensation loop

Find the palette color \mathbf{r}_i closest to the goal point \mathbf{x}_i. Increase \mathbf{r}_i’s weight by 1.

Move the goal \mathbf{x}_{i+1} = \mathbf{x}_i + (\mathbf{p} - \mathbf{r}_i). Repeat N times. Normalize weights by N.

When the process repeats for \mathbf{x}_{i+1}, the closest-point-query is done on the opposite side of \mathbf{p}. Selecting successive points this way approximates a convex hull.

I know the procedure may still seem abstract, but the point is that the algorithm is both simple and effective. There are some subtleties like sorting by brightness to keep the color choices consistent across pixels, but in Python it boils down to this: # Work arrays allocated outside the main loop weights = np.zeros(K, dtype=float) error = np.zeros(3, dtype=float)

for y in range(H): for x in range(W): color = img_float[y, x]

# Clear work arrays weights[:] = 0.0 # candidate weights weight_sum = 0.0 error[:] = 0.0 # accumulated error compensation

# Find candidates, can be less than N unique colors for _ in range(N): idx = find_closest(color + error * dither_strength) weights[idx] += 1 weight_sum += 1 error += color - palette_float[idx]

weights /= weight_sum

threshold = (bayer_4x4[y%4,x%4] + 0.5) / 16.0

# Output the candidate index on which cumulative weight sum # first crosses the threshold read from the 4x4 matrix # If a palette index wasn't used, its weight will be zero # and it won't be selected by this loop. cumulative_sum = 0.0 for idx in luma_order: w_ko = weights[idx] cumulative_sum += w_ko if cumulative_sum > threshold: break inds[y, x] = idx Notice how dither strength can be adjusted by modulating the scale of the error compensation. You can also study mateljou’s shader version. This should be enough context to understand the alternative solution presented next.

Yliluoma’s algorithms

In 2011, Joel Yliluoma presented a series of dithering algorithms as alternatives to Knoll’s. They were described on his website in an article titled Arbitrary-palette positional dithering algorithm. As far as I know, they haven’t gained much traction. I studied them carefully and found one particularly interesting. Yliluoma-2 simplified

I will focus on the second algorithm presented in the article. The variant in its C++ implementation, to be exact. Let’s call it Yliluoma-2. On a high level, Yliluoma-2 is similar to other N-candidate algorithms like Knoll’s: on each pixel, select weighted candidate colors from the palette, then output one based on a threshold matrix (or noise). The difference is how the candidates are selected, which in this case is done by testing every palette color using a unique error formula. No need for complex color difference formulas There’s a common belief that you need complex color difference calculations for high quality dithering. I think this is mistaken, and what you really want are (a) more weight for the green channel and (b) emphasis on differences of bright colors. The article proposes a luma-weighted (previously on this site) color difference formula, which I implemented. You can see its result in the first image from the left below. A simpler way to achieve the same thing is to desaturate your image and palette before dithering. The final indexed image still uses the original, colorful palette. For example libimagequant multiplies RGB colors by (0.5, 1, 0.45) and raises each channel to the power of 0.8. I did that in the second picture below, and it works just as well.

Of course, if you want a high-quality reconstruction, as in “looks like the original when squinting”, then do dithering in linear space. I’d also argue that if your output resolution is low and pixels big, you can do whatever you like. You’re constructing something very different-looking that you hope the viewer will interpret the same way as the original. Candidates are chosen via closest-point-to-line-segment tests Back to the main topic. I spent a pleasant afternoon with the code and arrived at a geometric interpretation of it.

Recall how N-candidate methods both find candidate palette colors and assign weights (probabilities) to them. The color selection loop in Yliluoma-2 keeps track of an exponential moving average (EMA) of the candidate colors chosen so far. In other words, on each search iteration, the routine updates the moving average \mathbf{x}_i with a new candidate \mathbf{r}_i via \mathbf{x}_{i+1} = (1 - t)\mathbf{x}_i + t\mathbf{r}_i, where t \in [0,1] is a mixing weight that varies on every iteration. (Perhaps t should be called t_i instead?) Okay, so the candidate is mixed in to the running average via lerp(). But how is the candidate chosen and where does the mixing factor t come from? The candidate point is chosen by testing line segments between every palette color \mathbf{c}_k and the current mean \mathbf{x}_i. The line that passes closest to the input color \mathbf{p} tells which palette color to choose as the candidate \mathbf{r}_i. The value of t is then the mixing factor that produced the closest point on the segment that passed nearest to the input. It’s hard to explain, but visually things should make more sense:

To summarize: of all palette colors, the chosen candidate best approximates the input when linearly combined with the current mean. Exponential moving average loop in detail Here’s the algorithm in more detail. Studying it is not necessary to understand the developments below, but I thought it prudent to document it. On the first iteration, the candidate is chosen as the closest palette color to the input: \mathbf{x}_1 = \text{closest}(\mathbf{p}). Its weight is 1. The rest of the iterations proceed like this:

Yliluoma-2’s exponential moving average loop (simplified)

Consider every palette color \mathbf{c}_k as a possible candidate.

Find a mixing factor t that minimizes the squared distance ||((1-t)\mathbf{x}_i + t \mathbf{c}_k) - \mathbf{p}||^2. If the distance was the shortest found so far, then:

Update t_\text{best} = t Update \mathbf{c}_\text{best} = \mathbf{c}_k.