18feb: Halftoning & Dithering

This week I experimented with different techniques in halftone printing using MATLAB. First I describe a naive method and then I will go into Floyd-Steinberg dithering

In newspaper grayscale printing, usually only 2 shades (black and white) available. But the lack of shades is made of for by high spatial resolution (which means very tiny and closely spaced dots can be printed). The idea is to print densely packed black dots to denote dark shades of gray, and print sparsely spaced dots to denote lighter shades. The resolution of the eye is not enough to detect each dot separately, hence it detects only an averaged and blurred out grayscale image. This helps produce a better looking image than simple thresholding, which produces strong contouring effects.

Activity 1: Naive halftoning

First I experiment with a simple halftoning idea. Given an image, map each pixel to say an 8x8 grid. The grid is filled with interspersed black and white dots. The density of the dots in the 8x8 grid is inversely proportional to the colour value of the original pixel. For example a black pixel would map to an 8x8 grid with all black pixel while a gray pixel maps to an 8x8 grid that looks like a chess board.

First lets see what the original looks like and what happens when we apply a simple threshold to it.

As we can see the results are not so good for simple thresholding. Now I map each pixel of the original image to 8x8 grid and create a new image with the following code:

function halftone = getHalftone(img, factor) %each pixel is mapped to a factor x factor grid

sz = size(img);

halftone = zeros(sz(1) * factor, sz(2) * factor);

pixelExpansion = zeros(factor, factor, 256);

for i = 1:256

%if colour is of intensity i, then it is i/256th of max. so it can have

%i/256th of pixels white. since there are a total of "factor" pixels in

%eaxh expanded pixel row, the number of pixel in each row is (i/256)*factor

val = ceil(((i - 1)/256)*factor);

temp = zeros(factor, 1);

for j = 1:val

if (round(factor/val) * j <= factor)

temp(round(factor/val) * j) = 255;

end

end

%create the map for each colour.

for x = 1:factor

pixelExpansion(:, x, i) = circshift(temp, x);

end

end

for y = 1:sz(1)

for x = 1:sz(2)

%Fill factor x factor blocks in the resulting image using the map created earlier

halftone((factor*(x-1) + 1):factor*(x-1) + factor, (factor*(y-1) + 1):factor*(y-1)+factor) = circshift(pixelExpansion(:, :, img(x, y) + 1), round(rand(1)*7));

end

end

end

Of course this image looks very bad! But wait, if we remember printers have great spatial resolution therefore what actually will happen is that the above halftone image will actually be displayed in a smaller area and then our eyes will not be able to distinguish individual dots. Only then will we see it in gray colours as was originally intended.

The above image shows the halftone image and the simple thresholding image, when I clicked them on my laptop screen from a distance. Now the halftone image looks better than the thresholded image as some shades of grey are visible. Note that the halftone image is of a higher resolution as explained earlier, but both images use only black(0) and white(255).

Activity 2: Floyd-Steinberg Dithering

Though the above halftone result was better than simple thresholding, it is still not satisfactory. So now let me describe a simple but effective technique that yields true images that have same resolution as the original image and look good too! It is the Floyd-Steinberg technique, whose details can be found here. The algorithm uses error diffusion. It scans the image in raster order. At the current pixel it quantises it to the nearest available gray level, and then distributes the error to the pixel to its right and on the row below it (ie those neighnbours that haven't yet been encountered in the raster scan). Below is the MATLAB code for a generalised dithering where we can choose over a given palette of available gray shades instead of only black and white. Of course the only black and white case is a special case of this function.

function [img, hardcolourmapped] = getHalftone(img, colourmap, mode)

sz = size(img);

hardcolourmapped = zeros(sz);

%map 0-255 colours to the given colormap.

lencolourmap = length(colourmap);

mappedcolours = zeros(255,1);

for i = 0:255

%subtract i from all elements of given colourmap

%find the one which gives minimum absolute difference.

%that is the value to which i will be mapped to.

[val, loc] = min(abs(repmat(i, 1, lencolourmap) - colourmap));

mappedcolours(i + 1) = colourmap(loc);

end

%generating the hard colour mapped image (without error spreading) ie simple thresholding

for y = 1:sz(1) %y coordinate

for x = 1:sz(2) %x coordinate

hardcolourmapped(x,y) = mappedcolours(img(x,y) + 1);

end

end

if (mode == 1)

figure; imshow(uint8(hardcolourmapped));

end

%begin dithering operation;

for y = 1:sz(1) %y coordinate (rows)

for x = 1:sz(2) %x coordinate (cols)

oldval = img(y, x);

if (img(y, x) > 255)

img(y,x) = mappedcolours(255 + 1);

else

img(y,x) = mappedcolours(img(y, x) + 1);

end

error = oldval - img(y, x);

%spread the error to the pixels that haven't been printed yet.

if (x + 1 <= sz(2))

img(y, x + 1) = img(y,x + 1) + round(error*7/16); %Pixel to the right

end

if (x - 1 >= 1 && y + 1 <= sz(1))

img(y + 1, x - 1) = img(y + 1, x - 1) + round(error*3/16); %Pixel on the bottom left

end

if (y + 1 <= sz(1))

img(y + 1, x) = img(y + 1, x) + round(error*5/16); %Pixel below

end

if (x + 1 <= sz(2) && y + 1 <= sz(1))

img(y + 1, x + 1) = img(y + 1, x + 1) + round(error*1/16); %Pixel on the bottom right

end

end

end

end

Lets see the results on a Lena image.

Clearly we can see that dithering halftone gives superior results. Notice that as the number of available gray shades increase, the simple threshold and the dithering techniques produce similar results, which is as expected. Let us study the psnr between the original image and the halftone and the thresholded images.

We observe that as the number of shades available for halftoning increases the PSNR decreases quite sharply. That tells us that if we are given between 10 to 20 shades only, we can do a pretty good job with this algorithm. That is evident if you look at the halftone image produced by 26 shades shown in an earlier table. Another observation is that the PSNR for halftone image (dithering) is slightly higher than that of the simple thresholded image. But visually the dithered image looks better. This is an example where PSNR fails to capture the "goodness" of an image.

Activity 3: Floyd-Steinberg dithering for aesthetic purposes

If you look at the images produced by the dithering technique, they look pretty nice, so I thought, I could use them for some image "enhancement" and produce aesthetically pleasing images. here are my results

The first image is the result of using the halftoning algorithm with 5 shades of grey. In the second I add a coloured background. The colour of the background changes based on the row number of the image linearly, much like my activity 2 on the 11th of feb page. In the third image, I sharpen the dark regions by superposing on the second image, a mask that I generate by a halftone image of 2 shades only.

The relevant code for the second image is:

halftoneimage1 = getHalftone(img, [0 255], 1); %halftone image with only 2 shades (will be used to as a mask to generate the 3rd image from the 2nd)

halftoneimage = getHalftone(img, [0 64 128 192 255], 1); %generate the 1st image

rgbImage = repmat(halftoneimage,[1 1 3]); %convert grayscale to RGB

s = size(rgbImage);

col1 = [255, 160, 0]; %These coulors will be used for the background

col2 = [255, 0, 0];

%generate the 2nd image

for i = 1:s(1)

for j = 1:s(2)

rgbImage(i,j,1) = uint8(0.6*rgbImage(i,j,1) + 0.4*(col1(1) * (i/s(1)) + col2(1) * (1 - (i/s(1)))));

rgbImage(i,j,2) = uint8(0.6*rgbImage(i,j,2) + 0.4*(col1(2) * (i/s(1)) + col2(2) * (1 - (i/s(1)))));

rgbImage(i,j,3) = uint8(0.6*rgbImage(i,j,3) + 0.4*(col1(3) * (i/s(1)) + col2(3) * (1 - (i/s(1)))));

end

end

%generate the 3rd image

for i = 1:s(1)

for j = 1:s(2)

if (halftoneimage1(i, j, 1) == 0) %use as a mask for sharpening

rgbImage(i,j,1) = 0;

rgbImage(i,j,2) = 0;

rgbImage(i,j,3) = 0;

end

end

end