Dithering

First of all: Happy New Year!

Over the holidays I’ve been learning about dithering, the process of creating the illusion of many grey levels using only black and white dots. This is used when displaying an image on a device with fewer than the 64 or so grey levels that we can distinguish, such as an ink-jet printer (which prints small, solid dots), and also when quantizing an image to use a colour map (remember the EGA and VGA computer displays?). It turns out that this is still an active research field. I ran into the paper Structure-aware error diffusion, ACM Transactions on Graphics 28(5), 2009, which improves upon a method presented a year earlier, which in turn improved on the state of the art by placing dots to optimize the appearance of thin lines. This got me interested in the basic algorithms, which I had never studied before. Hopefully this post will give an understanding of dithering and its history.

(As always, the code in this post is for MATLAB, and some of it requires DIPimage.)

Random dithering

The reason one needs dithering is quite obvious. If your display can only display black and white pixels, a grey-value image like this:

img = [newim(150,64)+1;newim(150,64)+64;newim(150,64)+85;newim(150,64)+127];
img = [readim('trui'),img,xx('corner')]

Grey-value image

will end up looking like this:

out = 255*(img > 127)

Thresholded image

Obviously, more is lost than necessary. What we want to do is put 50% of pixels to on in a middle-grey (127) region, 25% of pixels to on in a region of grey value 63, etc. The simplest way to accomplish this is using a random process, in which a pixel of value 127 has a 50% chance of becoming white, a pixel of value 63 has a 25% chance, etc. This two lines of code accomplishes this:

t = noise(newim(img),'uniform',0,255);
out = 255*(img > t)

Dithering with white noise

We simply replaced the constant threshold value of 127 with a random one, in the range 0 to 255. The problem with this output is that it looks noisy! We basically added white noise to the image, which has equal power at high and low frequencies. The goal is to avoid changing the low frequencies, so that, when viewed from an appropriate distance, the image looks identical to the original (the further away from the display you are, the more high-frequency components get lost). Thus, we want a high-frequency noise with no power at lower frequencies. Such noise is called blue noise, and can be generated simply by removing the low-frequency components of the noise image t:

t = hist_equalize(t-gaussf(t,2))*(254/255);
out = 255*(img > t)

Dithering with blue noise

Note that I used hist_equalize to make sure that the values in the image t are uniformly distributed. The result is still noisy, of course, but much more detail of the original image is preserved. This is because the low-frequency components of the image are not disturbed as much. I urge you to look at the Fourier transform of these various images and compare what happens at low and high frequencies in each.

Halftoning

Artists have been using techniques similar to dithering for many centuries: wood, leather or metal plates were etched, then inked and pressed against paper to produce a black and white impression. By creating small structures, the artists could create the impression of different levels of grey. With the advent of photography, an optical technique to create such an etching from a photograph was developed. This technique, halftoning, allowed photographs to be printed with ink on paper, for use in books and newspapers. In brief, the method involved projecting the image through a screen composed of tiny dots (such as a silk cloth). Brighter parts of the image will produce larger dots on the screen than darker parts. A photochemical process etches the dots into the metal plate. Nowadays this is done digitally, and at much higher resolution, but newspapers still use variable-sized dots to reproduce photographs.

Halftoning can be implemented in a quick-and-dirty way thus:

s = 2;
t = (cos((pi/3/s)*(xx(img)+0.2))*cos((pi/3/s)*(yy(img)+0.3))+1)*127;
out = 255*(gaussf(img,s) > t)

But please note, this is not the way it should be done!

Halftoning, s=2

I created a mask for thresholding that, similarly to earlier masks, is uniformly distributed in the range 0 to 255. However, the points now are ordered, so that the “on” pixels in a grey area are clustered together. The parameter s determines the size of the dots, and I’ve smoothed the image before thresholding so that the dots look more like dots (it’s not really necessary to do the smoothing). Because we have an ordered mask now, there is a limited number of grey values. Our dots grow with whole pixel increments. When reducing the size of the dots, the number of reproducible grey values is also reduced:

s = 1;
t = (cos((pi/3/s)*(xx(img)+0.2))*cos((pi/3/s)*(yy(img)+0.3))+1)*127;
out = 255*(gaussf(img,s) > t)

Halftoning, s=2

Ordered dithering

The quick-and-dirty halftoning I implemented above is a form of ordered dithering. The randomness of the noise mask is replaced by a fixed, repeated pattern. However, clustering the “on” pixels in this way creates strong artefacts (the dots!) that are not necessary for electronic displays. A printer benefits from the clustering because it becomes easier to control the amount of ink on the paper (ink tends to diffuse into the paper, for example). But for your computer’s screen a more pleasant effect is achieved when the “on” pixels are distributed as evenly as possible. Bayer (1973) descried a mask to be used for this purpose:

m = [
 1 , 49 , 13 , 61 ,  4 , 52 , 16 , 64
33 , 17 , 45 , 29 , 36 , 20 , 48 , 32
 9 , 57 ,  5 , 53 , 12 , 60 ,  8 , 56
41 , 25 , 37 , 21 , 44 , 28 , 40 , 24
 3 , 51 , 15 , 63 ,  2 , 50 , 14 , 62
35 , 19 , 47 , 31 , 34 , 18 , 46 , 30
11 , 59 ,  7 , 55 , 10 , 58 ,  6 , 54
43 , 27 , 39 , 23 , 42 , 26 , 38 , 22 ] * (255/65);

The mask contains 64 different values, thus it is possible to threshold the image at 64 different grey levels, producing a different pattern for each of 65 different grey levels. It is possible to use smaller subsets of this mask, such as [1,4;3,2], for 5 different grey levels. This pattern is repeated to create a mask of the same size as the image, then we threshold as before:

sz = size(img);
t = repmat(dip_image(m),ceil(sz/n));
t = t(0:sz(1)-1,0:sz(2)-1);
out = 255*(img > t)

Bayer dithering

This method is most commonly used in computer displays with limited palettes, because it is simple, fast, and deterministic (meaning it can be used in animations). Also, this is the best method to render large areas of uniform intensity, such as in GUIs. But for a photograph it leaves quite a bit to be desired. This is where error diffusion comes in.

Error diffusion dithering

In 1975, Floyd and Steinberg introduced error diffusion, a process where the error made at one pixel, when it is assigned a 0 or a 255, is passed on to its neighbours. Much like a real diffusion process, the error is spread so that it can be compensated for later on. The image is scanned row by row once. Each pixel visited is assigned a value, and the difference between that value and the input value is split among the neighbours not yet visited:

out = double(img);
sz = size(out);
for ii=1:sz(1)
   for jj=1:sz(2)
      old = out(ii,jj);
      new = 255*(old >= 128);
      out(ii,jj) = new;
      err = new-old;
         if jj<sz(2)
            % right
            out(ii  ,jj+1) = out(ii  ,jj+1)-err*(7/16);
         end
      if ii<sz(1)
         if jj<sz(2)
            % right-down
            out(ii+1,jj+1) = out(ii+1,jj+1)-err*(1/16);
         end
            % down
            out(ii+1,jj  ) = out(ii+1,jj  )-err*(5/16);
         if jj>1
            % left-down
            out(ii+1,jj-1) = out(ii+1,jj-1)-err*(3/16);
         end
      end
   end
end
out = dip_image(out)

Floyd-Steinberg dithering

This scheme produces a pseudo-random pattern, somewhere in between the fully random and fully ordered patterns we have seen so far.

It is a very, very simple algorithm, the code is only slightly complicated by the checking for out-of-bounds writes, and could be simplified by adding a dummy border to the image. Note, however, several structured patterns appearing at specific grey values (best visible in the ramp on the right side of the image). The four uniform blocks also show problems. These specific grey values are known to cause these problems, and is why I included them in the test image. A regular pattern is produced, but broken at some locations. These breaks produce visual artefacts that were not present in the input image. A third problem is that, sometimes, dots will be placed in little diagonal lines. These are called worms and significantly affect the visual quality of the output. Again, look at the Fourier transform of this image; you’ll notice some strong peaks at middle and high frequencies, caused by these artefacts.

Since this algorithm was published, many publications suggested improvements to the basic scheme: larger windows over which to diffuse the error, different weights to the diffusion kernel, changing the traversal direction on even and odd lines, etc. etc. The most simple and effective one is to add a little bit of randomness to the process, though this is unsuitable for some purposes. In the code above, change the threshold code to read:

      % new = 255*(old >= 128);
      new = 255*(old >= 128+(rand-0.5)*100);

Floyd-Steinberg dithering with noise

In 2001, Ostromoukhov proposed to change the diffusion coefficients depending on the input grey level. They determined “difficult” grey levels experimentally, and chose weights for them that optimize certain criterion in the output of a uniform image with that grey level (for those interested: they aimed at making the spectrum of the error image most similar to blue noise as possible). Weights for other grey levels were determined through linear interpolation. Thus, this algorithm is as fast as the original error diffusion, but with much better results. Zhou and Fang extended the Ostromoukhov algorithm with threshold modulation. The threshold they used also depends on the input grey level, and is determined experimentally for a set of key grey values and interpolated linearly for other grey values. The code below implements this Zhou-Fang dithering, which includes also the alternating left-right, right-left processing of rows. Again, if we were to add a dummy border to the image, the code would be a bit simpler.

coef = [ % [ index, right, left-down, down ]
  0,      13,      0,      5
  1, 1300249,      0, 499250
  2,  214114,    287,  99357
  3,  351854,      0, 199965
  4,  801100,      0, 490999
 10,  704075, 297466, 303694
 22,   46613,  31917,  21469
 32,   47482,  30617,  21900
 44,   43024,  42131,  14826
 64,   36411,  43219,  20369
 72,   38477,  53843,   7678
 77,   40503,  51547,   7948
 85,   35865,  34108,  30026
 95,   34117,  36899,  28983
102,   35464,  35049,  29485
107,   16477,  18810,  14712
112,   33360,  37954,  28685
127,   35269,  36066,  28664];
x = coef(:,1);
y = coef(:,2:4);
y = bsxfun(@rdivide,y,sum(y,2));
coef = [interp1(x,y(:,1),0:127)',interp1(x,y(:,2),0:127)',interp1(x,y(:,3),0:127)'];
coef = [coef;flipud(coef)];
strg = [
0,   0.00
44,  0.34
64,  0.50
85,  1.00
95,  0.17
102, 0.50
107, 0.70
112, 0.79
127, 1.00];
strg = interp1(strg(:,1),strg(:,2),0:127)';
strg = [strg;flipud(strg)];
img = double(img);
out = img;
sz = size(out);
for ii=1:2:sz(1)
   for jj=1:sz(2)
      old = out(ii,jj);
      new = 255*(old >= 128+(rand*128)*strg(img(ii,jj)+1));
      out(ii,jj) = new;
      err = new-old;
      weights = coef(img(ii,jj)+1,:); % +1 to account for MATLAB indexing into table
         if jj<sz(2)
            % right
            out(ii  ,jj+1) = out(ii  ,jj+1)-err*weights(1);
         end
      if ii<sz(1)
            % down
            out(ii+1,jj  ) = out(ii+1,jj  )-err*weights(3);
         if jj>1
            % left-down
            out(ii+1,jj-1) = out(ii+1,jj-1)-err*weights(2);
         end
      end
   end
   ii = ii+1; % even rows go right to left
   for jj=sz(2):-1:1
      old = out(ii,jj);
      new = 255*(old >= 128+(rand*128)*strg(img(ii,jj)+1));
      out(ii,jj) = new;
      err = new-old;
      weights = coef(img(ii,jj)+1,:); % +1 to account for MATLAB indexing into table
         if jj>1
            % right
            out(ii  ,jj-1) = out(ii  ,jj-1)-err*weights(1);
         end
      if ii<sz(1)
            % down
            out(ii+1,jj  ) = out(ii+1,jj  )-err*weights(3);
         if jj<sz(2)
            % left-down
            out(ii+1,jj+1) = out(ii+1,jj+1)-err*weights(2);
         end
      end
   end
end
out = dip_image(out)

Zhou-Fang dithering

Iterative dithering

I think this is as good as simple error diffusion can possibly be. More complex algorithms use an iterative approach, where dots are placed to minimize some error criterion. These algorithms are typically several orders of magnitude slower than error diffusion, and, of course, the most difficult thing about them is to design a correct error measure. As I mentioned at the beginning of this post, a recent paper took the improvements of an iterative algorithm, designed to optimize dot placement for preserving thin lines (“structure-aware dithering”), and developed an error-diffusion scheme with similar results. It uses even larger lookup tables than the Zhou-Fang scheme (about 1.3 Mb of them), together with frequency and orientation analysis of local neighbourhoods, to determine error diffusion weights and threshold modulation weights. The paper shows some fantastic results, but this method is way too complex for one blog post! The paper comes with the source code in the supplementary material, which I compiled and ran on my test image:

Structure-aware error diffusion dithering

Doing all the analysis takes several seconds on this small image, meaning that the whole process takes several seconds instead of being instantaneous. But then again, a few years down the road this will be instantaneous as well.

5 Responses to “Dithering”

  1. On March 7th, 2011, at 11:15, Filip Malmberg said:

    Ah, the nostalgia! This brings back memories of my old Macintosh SE.

    On the same note, Apple just released the source code for the classic MacPaint drawing program:

    http://www.computerhistory.org/highlights/macpaint/

  2. On August 28th, 2012, at 13:28, Matteo said:

    Hi Cris

    I’ve been following your blog for two years now but somehow I missed this post. Better late than never because it is a great post. For me dithering and half toning are among those topics that are there in a corner of my mind which I’d like to investigate, but never do. This post does a great job at explaining the basics with clarity and readers can see you had fun pushing it further. So, thank you. I was tickled by your encouragement to look at the Fourier transform of the first few images and had all the intention of trying, but by the time I finished reading and I had Matlab up and running I ended up playing around in a different direction. I thought it’d be interesting to see how well Canny and Sobel did at detecting edges on the last few results compared to the performance on the original. Another way to evaluate the results that I often use. I uploaded my results here:

    http://mycarta.wordpress.com/834-2/1465/edge-detection-on-dithering-results/

    The structure-aware dithering does a much better job at preserving the edges in the original and the Canny and Sobel pick up on this (as do our eyes). As you say, the process is far from instantaneous but it’s good to be aware of the advantages.

  3. On August 28th, 2012, at 15:46, Cris Luengo said:

    Hi Matteo,

    That’s a good point, and interesting way of evaluating visual fidelity. Thanks!

  4. On August 29th, 2012, at 12:24, Matteo said:

    Hi Cris,

    Thanks. I use this method often to evaluate effectiveness of filtering in the Fourier domain.

    I will soon show on my blog some results of a fun project using King Tut cat scan. I will keep you posted.

    Matteo

  5. On February 2nd, 2013, at 23:11, edge detection as image fidelity test | MyCarta said:

    […] post is a quick follow-up to Dithering, a very interesting post by Cris Luengo, developer of DIPimage, a free Matlab image analysis […]

Leave a Reply

You can use these HTML tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Note: I moderate all comments. Comments without a clear relation to the text above will not be published.