Introduction

Dithering is a fascinating and seemingly counter-intuitive technique. In order to demonstrate the technique, take a look at the following image. How many different shades of gray do you see in the picture?

How many did you guess? The answer is 2. Taking a closer look, we begin to see the pixels.

The original image is quantized to just two colors, fully black and white. How can an image with just black and white pixels look so good? The pixels are distributed just right, in such a way that, when you squint your eyes, you don’t really notice that the image actually consists of very few colors. The pixel distribution approximates the gray value of the original image. We have counter-intuitively improved an image with just two colors by intentionally introducing noise to it. It’s both an intellectually interesting effect, as well as a cool visual gimmick.

This is an old technique with which images used to be rendered on display devices which were severely limited in the number of different colors they could reproduce. Some early displays had only binary pixels, and printers to this day have to approximate different shades of colors using only a very small number of colored inks.

In order to learn the dithering process, we need to start from the beginning – by quantizing an image to just black and white, without any additional processing. Then we will make incremental improvements to the image using three different dithering methods: Bayer matrices, white noise, and blue noise.

While this post isn’t intended to be a programming language tutorial, the algorithm was implemented in a tiny Rust program. In order to keep the post general, the Rust implementation notes are included in the final section.

Quantizing the image

We would like to quantize an image to just black and white. The quantization algorithm is simple. We will convert each pixel value to black or white by comparing the original pixel to a threshold value. The threshold value is what determines whether the pixel will be black or white.

To start, we need the original image.

Next, we need a signal which we will quantize. We can use the image’s luma, \(l\), which is a weighed sum of the image’s \(R\), \(G\), \(B\) components. This signal is typically associated with the “black and whiteness” of an image. Even though the original image is already black and white, the pixels still contain all three RGB values. Using luma in our code enables us to also use color images as input for our dithering algorithm, and still produce sensible output images.

Finally, we need a threshold value, a number which chooses whether a luma value is converted into black (0.0) or white (1.0). Choosing 0.5, we obtain the following step function for each pixel \(p_{x,y}\).

Running our image through this function with using gives us the following result.

Not a very good looking result. Luckily, we already know that we can do better than this!

Dithering

The dithering algorithm is based on the quantization algorithm above and requires just one modification to it: it adds noise to the threshold values. This has the effect of distributing the harsh boundaries between black and white in the quantized image. The classical way of distributing noise is by using Bayer matrices.

Bayer matrices

A Bayer matrix, \(M\), looks like this.

The matrix is square and is characterized by its side length \(N\) which in this case is 4. The matrix consists of elements \((0,;1,;\ldots,;N^2-1/N^2)\) which are distributed in such a way that the distance between elements of the sequence is as large and as even as possible across the matrix. We can generate a Bayer matrix for side length 2, 8, or 16 as well. The last section of this post contains notes on how you can generate your own Bayer matrix.

For any given pixel coordinate \((x,,y)\) in our image, we can map the pixel coordinate to the Bayer matrix element \((i,,j)\) using the modulo operator: \((i,,j)=(y,\text{mod},N,;x,\text{mod},N)\). The idea is to use the matching Bayer matrix element as the threshold value while doing quantization, instead of just a single number like 0.5. Inserting the Bayer matrix element into the previous step function, we obtain

Running our image through this function using \(N=4\), we obtain the following result.

The details of the original image are now clearly visible, although we do have a visible grid-like pattern. We can dither with Bayer matrices using \(N=2,8,16\) as well.

Looking at the image dithered using \(N=2\), we can still see detail resembling the original images in places of finer detail, such as the snow banks. But in places where finer gradients exist, such as the sky, the result is closer to the original naive quantization. There is almost no perceptible difference between the images dithered using \(N=8\) and \(N=16\). Comparing the sky using \(N=4\) and \(N=8\), you can see that there is clearly more detail using \(N=8\). That is where the sweet-spot between detail and matrix size seems to lie.

Although the results using larger Bayer matrices look good, there is just no way of getting rid of the slight crosshatch artifact, arising from the deterministic way the values are distributed. Our next two dithering methods will attempt to get rid of the pattern, by introducing randomness.

Random noise

The simplest way to add some randomness to the threshold value is to replace the Bayer matrix lookup with a call to a uniform random number generator.

While we no longer have any visible crosshatch artifacts, the image has even less detail than in the previous section. Unlike the Bayer matrix, where neighboring values are guaranteed to be distributed evenly about the matrix, we can obtain multiple similarly valued uniform random numbers in sequence. As a result, our image loses detail due to the clustering of black and white. Nowhere is this more visible than when comparing the white plume from the smokestack with the images dithered with a Bayer matrix.

Blue noise

Blue noise is going to come to the rescue and restore the detail that we lost with the random noise. Blue noise values, like the random numbers used in the previous section, are distributed uniformly. But the frequency distribution of blue noise is different. Blue noise has less low-frequency components than the random numbers used in the previous section. Subsequent random numbers are more evenly distributed from each other.

This is what blue noise values look like in an image.

If we were to blur this image, it would look uniformly gray. This is the reason images look better when dithered with blue noise. We can verify that the pixel value distribution in blue noise is uniform by comparing it against uniform random numbers in a histogram.

Blue noise image masks, like the one above, are very expensive to compute. I used a free texture obtained from http://momentsingraphics.de/BlueNoise.html. These textures can be tiled, meaning that patterns due to repetition should be impercetable.

You can use the blue noise image mask exactly like a Bayer matrix. Load the pixels into a large matrix, and map pixel coordinates \(x\), \(y\) to image coordinates \(i\), \(j\) using the modulo operator. The resulting dithering function is exactly the same as the Bayer dithering function. Doing so, we obtain the final image, which was also the first image shown in this post.

Looking at the sky and snow detail, we are able to resolve far more detail than before, by simply adjusting the spatial distribution of our black and white pixels.

Implementation notes

Using the https://crates.io/crates/image crate, it’s remarkably easy to load, transform and output images. You can obtain a black and white image with two lines of code.

1
2

let user_img = image::open(opts.input_image).unwrap();
let user_img = user_img.grayscale();

With just a few more lines of code, it’s possible to transform, and output a result. This is a really fun way to get familiar with Rust.

Dithering function

Having loaded an image, it’s time to run its pixels through a dithering function. The dithering function uses a ThresholdMatrix struct and applies it to each pixel using the step function for Bayer dithering seen earlier. All dithering methods in this post are implemented using matrix look ups.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

fn dithered_rgb_image(
    threshold_matrix: &ThresholdMatrix,
    img: &image::DynamicImage,
) -> image::RgbImage {
    let (width, height) = img.dimensions();
    let img_luma = img.as_luma8().unwrap();
    let mut img_dithered = image::RgbImage::new(width, height);

    for (x, y) in (0..width).cartesian_product(0..height) {
        let luma = img_luma.get_pixel(x, y)[0];

        let dithered_pxl = if luma > threshold_matrix.look_up(x, y) {
            image::Rgb([255, 255, 255])
        } else {
            image::Rgb([0, 0, 0])
        };
        img_dithered.put_pixel(x, y, dithered_pxl);
    }

    img_dithered
}

The ThresholdMatrix is a matrix containing luma offset values, with dimensions nx and ny.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

struct ThresholdMatrix {
    nx: u32,
    ny: u32,
    matrix: Vec<u8>,
}

impl ThresholdMatrix {
       fn look_up(&self, x: u32, y: u32) -> u8 {
        let j = x % self.nx;
        let i = y % self.ny;
        let idx: usize = (i * self.nx + j)
            .try_into()
            .expect("i * side_length + j does not fit into usize");

        self.matrix[idx]
    }
 }

In order to implement the different dithering methods, all we need to do is to initialize the matrix with values in three different ways. The next sections adds associated functions to ThresholdMatrix to do just that.

Bayer matrix

Generating Bayer matrices is perhaps the most interesting part of the implementation. Starting from Wikipedia, we could just hard code the matrices presented in the article. But what if we want larger matrices than \(N=8\)? In that case, we need to generate our own matrices of arbitrary size.

Wikipedia mentions two ways of generating your own Bayer matrix: a recursive algorithm, and a matrix elementwise bit manipulation algorithm.

On the surface, the elementwise algorithm seems simpler to implement in Rust. However, the Wikipedia article seems to have misunderstood the meaning of the source. It says that you can bit reverse the bit interleaved values x ^ y and y. However, in the original source, Joel Yliluoma’s arbitrary-palette positional dithering algorithm, the bits are interleaved in reverse:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

// The original C algorithm

const unsigned dim = 1 << M;
for (unsigned y = 0; y < dim; ++y) {
    for (unsigned x = 0; x < dim; ++x) {
        unsigned v = 0, mask = M - 1, xc = x ^ y, yc = y;
        for (unsigned bit = 0; bit < 2 * M; --mask) {
            v |= ((yc >> mask) & 1) << bit++;
            v |= ((xc >> mask) & 1) << bit++;
        }
        printf("%4d\n", v);
    }
}

This is not quite the same as bit-interleaving xc, y, and reversing the bits of that. Doing so would give you really large integer values, closer to the maximum value of the integer than to \(N\).

The Rust version of this algorithm isn’t totally straightforward, since there are a few C things which can’t just be copy-pasted over.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

// Rust version of the algorithm

impl ThresholdMatrix {
    fn bayer_matrix(power_of_two: u32) -> Self {
        let side = 2_u32.pow(power_of_two);
        let num_elements = side * side;
        let norm_factor = 255_f32 / (num_elements as f32);
        let matrix = (0..side)
            .cartesian_product(0..side)
            .map(|(x, y)| {
                let xc = x ^ y;
                let yc = y;

                let mut v = 0;
                for p in (0..power_of_two).rev() {
                    let bit_idx = 2 * (power_of_two - p - 1);
                    v |= ((yc >> p) & 1) << bit_idx;
                    v |= ((xc >> p) & 1) << (bit_idx + 1);
                }
                (v as f32 * norm_factor) as u8
            })
            .collect();

        Self {
            nx: side,
            ny: side,
            matrix,
        }
    }

// ...

Random noise

Random noise is implemented by generating \(x\times y\) random numbers and storing them in the matrix. \(x\) and \(y\) need to be the image dimension in order for this to work.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

fn white_noise(xdim: u32, ydim: u32) -> Self {
    let mut matrix = Vec::new();
    for _ in 0..(xdim * ydim) {
        matrix.push(rand::thread_rng().gen::<u8>());
    }

    Self {
        nx: xdim,
        ny: ydim,
        matrix,
    }
}

Blue noise

Finally, blue noise is “generated” by loading the blue noise image into the matrix.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

fn blue_noise() -> Self {
    let texture = image::load(
        std::io::Cursor::new(include_bytes!("128x128_blue.png")),
        image::ImageFormat::Png,
    )
    .unwrap()
    .grayscale();
    let dim = texture.dimensions();

    let matrix = (0..dim.0)
        .cartesian_product(0..dim.1)
        .map(|(x, y)| texture.get_pixel(x, y)[0])
        .collect();

    Self {
        nx: dim.0,
        ny: dim.1,
        matrix,
    }
}

References

https://en.wikipedia.org/wiki/Ordered_dithering

Joel Yliluoma’s arbitrary-palette positional dithering algorithm

Free blue noise textures

Image quantization and plots were done using Python.

Image quantization notebook

Noise distribution plot notebook