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 outAnd 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.
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] = idxThis 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. 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) | |
|---|---|
|
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.
How to find t
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.
An analytic t solution
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.
The same as code
c_k = palette_float[k] # B
delta = c_k - mean # B - A
t = 0.0
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^2Limits to t
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. Let’s go one step further.
Fixed t
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.
Three variants and some observations
Since I’ve modified the original Yliluoma-2 approach, I dubbed the three variants above as “EMA” color selection algorithms:
- EMA-Sweep: This is Yliluoma-2 reimplemented in Python. Uses a 4x4 threshold matrix and tests 8 values of t in [0.2, 1.0). The C++ code tests 1–16 values in (0.0, 1.0) based on the number of candidates chosen so far.
- EMA-Exact: Calculates an exact solution for t \in [0.2, 1.0). Does less work than the original.
- EMA-Constant: Sets a fixed t=0.3. Simplest of the three.
Time complexities of the new variants
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 Yliluoma-2. EMA-Exact and EMA-Constant test only a single t value per candidate, so they are faster but in the same O(NK) complexity class.
Results and discussion
In these images, EMA and Knoll’s algorithms use Euclidean sRGB distances, while Yliluoma-2 uses its luma-weighted distance.
First of all, EMA-Sweep and EMA-Exact produce close to identical results. You can see the two images below are indistinguishable. The t=0.3 constant chosen for EMA-Constant creates slightly weaker dithering. See the supplement for a demonstration.
The exact ts chosen by EMA-Exact produce results similar to Yliluoma-2’s parameter sweep. Since the solved t is clamped to [0.2, 1.0), the resulting dithering is also a tad weaker than Yliluoma-2’s.
Below is an EMA-Exact vs Yliluoma-2 comparison with 16 colors, N=16 candidate iterations, and a 4x4 threshold matrix in both. They look mostly the same.
The number of candidates N has a large influence on the noisiness of the result. This applies to all N-candidate methods. The comparison below demonstrates with EMA-Exact how much worse 16 candidates look than 32. See for example the back of the yellow-turquoise macaw.
At high enough candidate counts, the EMA-Exact and EMA-Constant variants are comparable to Knoll’s algorithm at 20% and 10% dither strengths, respectively. At least that’s the case in the Chrono Cross screenshot used as an example in Yliluoma’s 2011 article.
Below is a 16-color comparison with the fixed palette used in the article. I used N=32 candidate iterations, so Yliluoma-2 was left out of this as its C++ code does only N=16. See the supplement for a direct N=16 comparison against Yliluoma-2.
Conclusion
Yliluoma’s second color selection algorithm does not need a perceptual color difference function to work. The algorithm calculates repeated weighted averages, and the chosen color candidate set’s quality is determined by the mixing factor t along with the number of candidates N. The optimal value for t can be solved exactly when using an Euclidean color difference formula. A suitable N for a K-color image seems to be around N=2K.
The simplified “EMA” variants yield comparable quality to the original despite being faster. While these variants occasionally overtake Knoll’s algorithm in quality, this isn’t the case most of the time. Also, the variants inherit their predecessor’s structure and thus have to test all K palette colors on each candidate search iteration. Knoll’s algorithm, on the other hand, can use a faster O(\log K) closest point lookup.
For cases when only light dithering is needed, the EMA-Constant variant presents a new, simple alternative. But since it’s still slower than Knoll’s algorithm, it’s difficult to recommend it.
Material
The article supplement includes more result images and an artificial 2D comparison at the end.
A standalone script that implements the Knoll, EMA-Sweep, EMA-Exact, EMA-Constant, and Offset ordered dithers with a 4x4 threshold matrix:
Usage: uv run color_selection.py bigbuckbunny_bird.png 16 32 --algo=ema-exact
My doctored C++ version of the original Yliluoma-2 algorithm:
This article is based on research I did for a book I’m writing. It’s about color reduction algorithms. Sign up here if you’re interested.