Image dithering in Matlab

Dithering

Dithering is a method used to reduce the color palette (color quantization) used for displaying an image, while the viewer still perceives almost the same image. The technique is based on approximating colors that are not on the palette, by diffusing colors which are available. By juxtaposing pixels of two colors, an illusion that a third color (the mixture of the two diffused colors) is present is created. Seven common diffusion algorithms are presented (and are available to download as Matlab code below). We apply those dithering methods on two images and compare the results. The two images used for the experiments are shown below (Lena and Wool). However, I present the results of the dithering algorithms only on Lena. Download the package to see the results on the other image (or any other picture you want).

Download

The package contains the following files:

main.m: Applies the different methods on the two sample images and reports the mean squared error

Fixed threshold

The first dithering method used compares each pixel with a fixed threshold. If its intensity is less or equal to the threshold, the pixel is assigned to black. If its intensity is higher than this threshold, the pixel is assigned to white.

Random threshold

By adding discrete uniform noise to the image prior to performing the threshold one can remedy some of the drawbacks of thresholding. The generated dithered is free of any artifacts. In practice this method compares each pixel with a random threshold. This type of dithering works particularly well on low-frequency images due to the absence of artifacts on the resulting picture.

Ordered threshold

In the ordered threshold method instead of using a fixed threshold for the whole image (or random values for each pixel), a threshold matrix (of usually small size, less than 10) of a predefined pattern is applied. The matrix paved over the image plane and each pixel of the image is thresholded by the corresponding value of the matrix. An example of such a matrix is shown below:
To apply the ordered threshold method, we first create a filter of the same size as the image by replicating the threshold matrix S on both directions as many times needed, depending on the same of the image. Specifically we replicate ts = (tsx,tsy) times, where is the number of times replicated along x axis and the number of times replicated along y axis. ts = ⌈size(image) / size(S)⌉ where size(A) is the (height, width) of A (division is element-by-element). Then we discard the (potential) extra few rows/columns of filter Simg by keeping just the rectangle:(1, size(image,x)) * (1,size(image,y))
So in the end we have a filter Simg of the same size as original image. The next step is to quantize the image to N gray levels, where N-1 is the number of thresholds in the matrix S (N = 37 in this case). The maximum value in S is N-2 and the minimum is 0. So there are N-1 threshold values. Obviously, the number of levels is 1 more than the amount of threshold values. There should be one gray level above the max threshold value and one below the min threshold value. Since the min value is 0, we shift the values of S by 0.5 (S = S + 0.5), so that we can correctly threshold the image. By doing this we can efficiently quantize the image without needing use floor/round/ceil operations. The quantization step is then simply defined as D = 255 / (N-1) and the quantized image is Q = I/D where is the original image and the division is applied element-by-element.
Finally we threshold each pixel position of the original image with the corresponding value of Simg.

A number of different matrices can be used to produce the dithered image. Below we present 5 different matrices and their corresponding results:

Clustered dots

Central white point

Balanced centered point

Diagonal ordered matrix with balanced centered points

In this case the filter used is rhombus-shaped as shown below.
This can be generated by combining these two matrices:
like this:

Dispersed dots / Bayer

This last ordered threshold technique was proposed by Bayer and is based on a sparse pattern. The filter matrix is defined as:
where D can be (more variations of different sizes have also been proposed):
and U is a matrix of the same size with all elements equal to 1.

Error diffusion

The error diffusion method works as follows. Each pixel in the original image I is thresholded (at 50% for black and white dithering) and the difference between the original image value and the thresholded value is calculated. This is the error value err(i,j) = I(i,j) – T(i,j). This error value is distributed to the neighboring pixels that haven’t been processed yet according to the diffusion filter F. The proportion of the error value that corresponds to those adjucent pixels is added to the original value of the pixel and this modified value is used for processing the pixel as we scan the image line by line. Two different diffusion filters are given in the library:
The Floyd & Steinberg:
and the Stucki:

As can be seen on the results below, this family of dithering methods (error diffusion) gives the best results in terms of visual quality.

Evaluation

Complexity

For all complexities stated below n is the number of pixels in image (width*height).Fixed Threshold
Complexity: O(n), since we threshold each pixel with a value.
Space: O(n), we only need to store the image (the output image can as well use the same matrix, not done in our implementation for more code readability).Ordered Threshold (All methods)
First, we create a matrix with the same dimensions as the image, by replicating the filter, using repmat(). The cost is O(n/k), where k=size(image)/size(filter). Then we compare each pixel of the two images which takes O(n) time.
Complexity: O(n+n/k) = O(n), since we threshold each pixel with the value of the same pixel in the big filter.
Space: O(2n) = O(n). We need to store both the image and the big filter.Error Diffusion
Complexity: O(n*m), where m is the size of the filter.
Space: O(n+m) since we need to store the original image and the filter. In our implementation we use a different matrix to store the output values, but that is not necessary. Since (usually) n >> m, space complexity is O(n).
The table below summarizes the complexities.

Time complexity

Space complexity

Fixed Threshold

O(n)

O(n)

Ordered Threshold

O(n)

O(n)

Error Diffusion

O(n*m)

O(n+m)

Mean Squared Errors and Visual Quality

The mean squared errors of all methods, for both images and their average, can be found below.
The interesting thing to notice is that even though the MSE of the fixed threshold less than that of all other methods, the images actually look much “worse” than the rest. All the other methods have comparable MSEs.