Floyd Steinberg Dithering Matlab - What am I doing wrong?

Question

I am trying to implement Floyd Steinberg Dithering in MATLAB, using the pseudocode on the Wikipedia page https://en.wikipedia.org/wiki/Floyd%E2%80%93Steinberg_dithering

My code is below

image = double(imread("dithering.jpg")) ./ 255;
levels = 2;

image_quantised = round(image .* (levels - 1)) ./ (levels - 1);
error = image - image_quantised;

height = size(image(:, :, 1), 1);
width = size(image(:, :, 1), 2);

image_dithered = image_quantised;

for y = 1:height - 1
    for x = 2:width - 1
                
        image_dithered(y    , x + 1, :) = image_dithered(y    , x + 1, :) + error(y, x, :) .* 7 / 16;
        image_dithered(y + 1, x - 1, :) = image_dithered(y + 1, x - 1, :) + error(y, x, :) .* 3 / 16;
        image_dithered(y + 1, x    , :) = image_dithered(y + 1, x    , :) + error(y, x, :) .* 5 / 16;
        image_dithered(y + 1, x + 1, :) = image_dithered(y + 1, x + 1, :) + error(y, x, :) .* 1 / 16;
        
    end
end

imshow(image_dithered)  % Image 1
imshow(dither(mean(image, 3)))  % Image 2

Image 1

Image 2

I am expecting the result in Image 2, but I am getting Image 1. It looks as though the algorithm isn't doing anything. Any ideas? :)

Edit: I have tried initialising image_dithered with different values; all zeros, the quantised image, and the original image. None of them work correctly

Edit 2: I'm getting closer by now calculating the error and quantisation within the loop. Still not spot on however.

for y = 1:height - 1
    for x = 2:width - 1
        new_pixel = round(image_dithered(y, x, :) .* (levels - 1)) ./ (levels - 1);
        image_dithered(y, x, :) = new_pixel;

        error = image(y, x, :) - new_pixel;

        image_dithered(y    , x + 1, :) = image_dithered(y    , x + 1, :) + error .* 7 / 16;
        image_dithered(y + 1, x - 1, :) = image_dithered(y + 1, x - 1, :) + error .* 3 / 16;
        image_dithered(y + 1, x    , :) = image_dithered(y + 1, x    , :) + error .* 5 / 16;
        image_dithered(y + 1, x + 1, :) = image_dithered(y + 1, x + 1, :) + error .* 1 / 16;
    end
end

Edit 3: Thanks for the @saastn and @Cris Luengo, both of the answers helped me work out where I was going wrong and it appears to be working as expected now!

The fixed code is below for completeness.

height = size(image(:, :, 1), 1);
width = size(image(:, :, 1), 2);

image_dithered = image;

for y = 1:height - 1
    for x = 2:width - 1
        old_pixel = image_dithered(y, x, :);
        new_pixel = round(image_dithered(y, x, :) .* (levels - 1)) ./ (levels - 1);
        
        image_dithered(y, x, :) = new_pixel;
        
        error = old_pixel - new_pixel;
        
        image_dithered(y    , x + 1, :) = image_dithered(y    , x + 1, :) + error .* 7 / 16;
        image_dithered(y + 1, x - 1, :) = image_dithered(y + 1, x - 1, :) + error .* 3 / 16;
        image_dithered(y + 1, x    , :) = image_dithered(y + 1, x    , :) + error .* 5 / 16;
        image_dithered(y + 1, x + 1, :) = image_dithered(y + 1, x + 1, :) + error .* 1 / 16;
        
    end
end

imshow(image_dithered)
imshow(dither(mean(image, 3)))

Answer 1

The problem is that you tried to improve the pseudocode and remove oldpixel! Note that the algorithm does not calculate an error between the quantized pixel and its corresponding value in the original image, but rather the error between the quantized pixel and the previous value in the dithered image, which may have been updated while scanning the previous pixels. Bring oldpixel back and review the whole algorithm one more time.

But even after the modifications, you can not expect the results to match the MATLAB output, which could be the result of differences in the details of the two implementations.

Answer 2

You’re propagating errors in image_quantized instead of in the original image. Remove this quantized image, it is not part of the algorithm.

You need to quantize one pixel, then find the difference with the original value, and propagate that error into future pixels.

Note that Wikipedia pseudocode does this in-place, there is only one copy of the image that works as both input and output.

On this old blog post of mine has MATLAB code for Floyd-Steinberg dithering.

Answer 3

Assuming you first follow saastn's answer (and replying to what it says about the differences between those images) : I'd say that by just looking at them in saastn's image, Matlab's patterns look like they're a sideways version of Floyd-Steinberg here, either by rotation, or more probably by transposition (swap x with y, which is a reflexion across the x+y diagonal axis). This can't be imitated by just changing coefficients, the image has to be processed by columns ("for x" outside of "for y", etc).