Revisiting Yliluoma’s ordered dither algorithm

# Revisiting Yliluoma’s ordered dither algorithm Source: [https://30fps.net/pages/revisiting-yliluoma-2/](https://30fps.net/pages/revisiting-yliluoma-2/) **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 ditheringToday 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 ofN*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*](https://matejlou.blog/2023/12/06/ordered-dithering-for-arbitrary-or-irregular-palettes/)that explains everything you need to know\. Please read at least its[“The Probability Matrix”](https://matejlou.blog/2023/12/06/ordered-dithering-for-arbitrary-or-irregular-palettes/#n-candidate-methods)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 candidatesw\_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\}\_1to the palette color closest to the input color\\mathbf\{p\}\. This palette color\\mathbf\{r\}\_1is the first candidate color\. Then the algorithm measures how much error there is between\\mathbf\{r\}\_1and\\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\}\_2is 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\. AfterNrounds 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 loop1. Find the palette color\\mathbf\{r\}\_iclosest to the goal point\\mathbf\{x\}\_i\. 2. Increase\\mathbf\{r\}\_i’s weight by 1\. 3. Move the goal\\mathbf\{x\}\_\{i\+1\} = \\mathbf\{x\}\_i \+ \(\\mathbf\{p\} \- \\mathbf\{r\}\_i\)\. 4. RepeatNtimes\. 5. Normalize weights byN\. 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](https://matejlou.blog/2023/12/06/ordered-dithering-for-arbitrary-or-irregular-palettes/#appendix-i-candidate-sorting)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](https://www.shadertoy.com/view/dlcGzN)\. This should be enough context to understand the alternative solution presented next\. ## Yliluoma’s algorithms [](https://30fps.net/pages/revisiting-yliluoma-2/yliluoma2-full.jpg)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*](https://bisqwit.iki.fi/story/howto/dither/jy/)\. 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](https://bisqwit.iki.fi/story/howto/dither/jy/#YliluomaSOrderedDitheringAlgorithm%202)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](https://30fps.net/pages/better-srgb-to-greyscale/)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](https://pngquant.org/lib/)multiplies RGB colors by\(0\.5, 1, 0\.45\)and raises each channel to the power of0\.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\)*](https://en.wikipedia.org/wiki/Exponential_smoothing)of the candidate colors chosen so far\. In other words, on each search iteration, the routine updates the moving average\\mathbf\{x\}\_iwith a new candidate\\mathbf\{r\}\_ivia\\mathbf\{x\}\_\{i\+1\} = \(1 \- t\)\\mathbf\{x\}\_i \+ t\\mathbf\{r\}\_i, wheret \\in \[0,1\]is a mixing weight that varies on every iteration\. \(Perhapstshould be calledt\_iinstead?\) Okay, so the candidate is mixed in to the running average via**lerp\(\)**\. But how is the candidate chosen and where does the mixing factortcome from? The candidate point is chosen by testing line segments between every palette color\\mathbf\{c\}\_kand 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 oftis 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 is1\. The rest of the iterations proceed like this: Yliluoma\-2’s exponential moving average loop \(simplified\)1. Consider every palette color\\mathbf\{c\}\_kas a possible candidate\.- Find a mixing factortthat 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:- Updatet\_\\text\{best\} = t - Update\\mathbf\{c\}\_\\text\{best\} = \\mathbf\{c\}\_k\. - Updatek\_\\text\{best\} = k\. 2. Choose the candidate\\mathbf\{r\}\_i=\\mathbf\{c\}\_\\text\{best\}\. 3. Move the mean point towards the candidate via the best mixing factor:\\mathbf\{x\}\_\{i\+1\} = \(1 \- t\_\\text\{best\}\)\\mathbf\{x\}\_i \+ t\_\\text\{best\}\\mathbf\{r\}\_\{i\} 4. Addt\_\\text\{best\}to candidate\\mathbf\{r\}\_i’s weight\. 5. RepeatNtimes\. Normalize weights via`weights /= sum\(weights\)`\. 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\}\_iand 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 oft \\in \[0,1\]and keep the one that was closest to the input point\. Here’s a diagram\. We are trying to find the pointDthat’s closest toCon the line segmentAB\. Below, the red dots stand for individual values oftwe 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 fortand solve the geometry problem directly\. ### An analytic*t*solution AssumingA \\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 betweenAandB, 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^2 ``` ### Limits to*t* As the exact solution can be anywhere in\[0,1\], it can pick really weird points\. For examplet=0equals 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 found0\.2to 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 choset = 0\.3and 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 oftin\[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 fort \\in \[0\.2, 1\.0\)\. Does less work than the original\. - **EMA\-Constant**: Sets a fixedt=0\.3\. Simplest of the three\. ### Time complexities of the new variants AssumeKcolors andNcandidates\. Per\-pixel work of Knoll’s algorithm isO\(N \\log K\)assuming closest point lookups accelerated with a kd\-tree\. Yliluoma\-2 has to doO\(N^2 K\)work because the number of distinct values ofttested seems to depend onN\(not totally clear from the code\)\. EMA\-Sweep tests 8 points on a segment for every candidate, so it’sO\(NK\*8\) = O\(NK\)but basically the same speed as Yliluoma\-2\. EMA\-Exact and EMA\-Constant test only a singletvalue per candidate, so they are faster but in the sameO\(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\. Thet=0\.3constant chosen for EMA\-Constant creates slightly weaker dithering\. See the supplement for[a demonstration](https://30fps.net/pages/revisiting-yliluoma-2/supplement.html#effect-of-ema-constants-fixed-t-value)\. The exactts chosen by EMA\-Exact produce results similar to Yliluoma\-2’s parameter sweep\. Since the solvedtis 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 candidatesNhas 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](https://30fps.net/pages/revisiting-yliluoma-2/supplement.html#chronocross)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 candidatesN\. The optimal value fortcan be solved exactly when using an Euclidean color difference formula\. A suitableNfor aK\-color image seems to be aroundN=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 allKpalette colors on each candidate search iteration\. Knoll’s algorithm, on the other hand, can use a fasterO\(\\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](https://30fps.net/pages/revisiting-yliluoma-2/supplement.html)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: - [color\_selection\.py](https://30fps.net/pages/revisiting-yliluoma-2/color_selection.py) Usage:`uv run color\_selection\.py bigbuckbunny\_bird\.png 16 32 \-\-algo=ema\-exact` My doctored C\+\+ version of the original Yliluoma\-2 algorithm: - [yliluoma2\.cpp](https://30fps.net/pages/revisiting-yliluoma-2/yliluoma2.cpp) *This article is based on research I did for a book I’m writing\. It’s about color reduction algorithms\.[Sign up here](https://paletteprogramming.com/)if you’re interested\.* ### Further reading - K\. Lemström, J\. Tarhio & T\. Takala:[*Color Dithering with n\-Best Algorithm*](https://www.researchgate.net/publication/2945114_Color_Dithering_with_n-Best_Algorithm), 1996\. - K\. Lemström & P\. Fränti:[*N\-Candidate methods for location invariant dithering of color images*](https://cs.joensuu.fi/sipu/pub/N-candidateDithering-IVC.pdf), 2000\. - Joel Yliluoma,[*Arbitrary\-palette positional dithering algorithm*](https://bisqwit.iki.fi/story/howto/dither/jy/), 2011\. - mateljou,[*Ordered Dithering with Arbitrary or Irregular Colour Palettes*](https://matejlou.blog/2023/12/06/ordered-dithering-for-arbitrary-or-irregular-palettes/), 2023\. - mateljou,[Thomas Knoll’s Pattern Dithering](https://www.shadertoy.com/view/dlcGzN)shader, 2023\.