Revisiting Yliluoma's ordered dither algorithm

Revisiting Yliluoma's ordered dither algorithm

Comments

June 29th, 2026

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.

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:

Today we are discussing how exactly it's done. 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 = (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:

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.

Yliluoma-2 Algorithm

This should be enough context to understand the alternative solution presented next. 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. 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.

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.

Geometric Interpretation

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.

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:

Normalize weights via

Like Knoll's algorithm, Yliluoma-2 also ends up approximating a convex hull. On each iteration, the goal $\mathbf{x}_{i+1}$ is moved to the closest point on the segment that was found earlier. Because this point is in between the old goal $\mathbf{x}_i$ and the input color $\mathbf{p}$, it must be on the opposite side (loosely speaking) of $\mathbf{p}$. That's why the goal, the running average, moves around $\mathbf{p}$.

In the above explanation, I tried to make the procedure easy to understand by choosing a formulation with explicit averaging steps. The original C++ code uses a running average instead, presumably for speed.

Finding the Closest Point

Above, I left it open as to how exactly the closest point on a line is found. The simplest way is to do a parameter sweep: test different values of $t \in [0,1]$ and keep the one that was closest to the input point. Here's a diagram. We are trying to find the point D that's closest to C on the line segment AB. Below, the red dots stand for individual values of $t$ we are testing:

Unfortunately, a parameter sweep is unlikely to hit the exact optimal $t$ value. It also takes time, but if you have a complex color difference formula, it might be the only thing that works. But remember how we established that luma weighting can be done via preprocessing? That's why we can ignore the perceptual stuff when searching for $t$ and solve the geometry problem directly.

Assuming $A \neq B$, we get

$$t_\text{closest} = \frac{(C - A) \cdot (B - A)}{||(B - A)||^2}$$

where "$\cdot$" denotes a dot product. Since we want the point in between A and B, we have to clamp the range to $[0,1]$:

$$t = \max(0, \min(1, t_\text{closest}))$$

and like earlier, $D = (1-t)A + tB$.

c_k = palette_float[k]  # B
delta = c_k - mean  # B - A
delta2 = delta.dot(delta)  # ||B - A||^2

# Avoid division by zero
if delta2 >= 1e-9:
    t = (p - mean).dot(delta) / delta2
    # Limit the choice to the given range
    t = max(min_t, min(max_t, t))
    mix = (1-t) * mean + t * c_k  # D
    diff = p - mix
    cand_dist = diff.dot(diff)  # d^2

As the exact solution can be anywhere in $[0,1]$, it can pick really weird points. For example $t=0$ equals not moving the mean at all, which is probably not what you want. That's why I experimented with various minimum values for $t$, and found 0.2 to be a decent compromise. So in practice, my code for the closest point calculates this instead:

$$t = \max(\textcolor{red}{0.2}, \min(1, t_\text{closest})$$

We have already strayed far from the original Yliluoma-2 algorithm that has no such limit.

EMA Variants

Let's go one step further. In my experiments, I noticed that the optimal $t$ was near zero 97% of the time. This makes sense if you consider a situation where the mean color has drifted near the input color; keeping it still minimizes error. This motivated an alternative method where instead of solving for $t$, it's just always a constant. I chose $t = 0.3$ and in most cases the results are indistinguishable from the exact method.

This approach can be thought of as choosing candidates that are in the direction of the input color, but quite far. In practice its output quality seems on par with the exact method, perhaps slightly worse.

Since I've modified the original Yliluoma-2 approach, I dubbed the three variants above as "EMA" color selection algorithms:

Assume K colors and N candidates. Per-pixel work of Knoll's algorithm is $O(N \log K)$ assuming closest point lookups accelerated with a kd-tree. Yliluoma-2 has to do $O(N^2 K)$ work because the number of distinct values of $t$ tested seems to depend on N (not totally clear from the code). EMA-Sweep tests 8 points on a segment for every candidate, so it's $O(NK*8) = O(NK)$ but basically the same speed as Yliluo