Rerunning the Floyd Steinberg algorithm on a image

Question

I implemented a version of the Dithering algorithm in Python (from Image Processing). The algorithm that I used is the Floyd Steinberg algorithm.

I was wondering how the images would change if I will rerun my algorithm on the same image over and over, repeatedly. I noticed that it didn't change at all:

First iteration:

10th iteration:

First of all, is it the correct behavior or something is up with my implementation? It it's correct, I was wondering why after one iteration, it does not do change to the image at all? Is there some math-explanation behind it?

Answer 1

Dithering by itself won't affect anything if you already start with an image that's in the target bit-depth: for error-diffusion dithering there's nothing to quantize further, and for ordered dithering the shifts are always less than 0.5 in magnitude if input and output depths are the same.

However, there is a subtle point that if you take that lowered-bit depth image, upconvert to a higher-bit depth image (as you usually do when you create an 8bpc PNG to display on the web), then run dithering again, your image will be modified.

This is because as explained here increasing bit-depth isn't just a simple multiply by power of two (bit shift left). If you want to properly use the full dynamic range you have to scale it by (outlevels - 1)/(inlevels - 1), e.g. mapping in=(0, 63) -> out=(0, 255). This is usually approximated as a bit shift left and replicate of higher order bits. It's easy to see then that this actually adds back some error that can be quantized again.

Concretely, if you implement quantization from 8 to 3-bits as round((x / 255.0) * 7.0) then 36 gets mapped to 1 (you can also bit shift right in this case, since we don't care about lower order bits anyway). And conversely upconversion from 3 to 8-bits is round((x / 7.0) * 255) and 1 gets mapped to 36. But now notice that starting with the 3bpc sample of value 1, converting to 8pc as 36, and calculating theoretical quantization error from conversion back to 3bpc gives 1 - (36/255)*7 = 0.98 which is non-zero.

This fact can be really counterintuitive since it seems like upconversion should be lossless, but if you think you realize it can't be: Take 2bpc, which means you have 4 different levels. We usually like to fix 0 as the lowest and 1 as the highest representable level, this gives us light levels of {0, 1/3, 2/3, 1}. Since the entire space of representable light is assigned one of these 4 levels, it splits the space into fourths [^1]. Likewise with 8bpc we split the space into 256, at gradients of 1/255. A single sample thus represents a "range" of actual light values that cannot be distinguished at the given precision.

[^1] The assignment of ranges to levels may not be even, depending on the rounding scheme you use). E.g. you could choose to do normal rounding, which means [0, 0.5) is assiged to level 0 while [0.5, 1.5) is assigned to 1. In such a case, fewer light levels are mapped to 0 compared to 1. Or you could be uniform and assign [0, 0.25) to first level, [0.25, 0,5) to second level and so on. With this scheme note that [0, 0.5) gets mapped to a value lower than the interval's average while [0.75,1] gets mapped to a value higher than its average, which is a necessary fact that we try to fix the endpoints while maintaining equal spacing. Gamma encoding could be seen as a way to assign buckets that in a way that's optimally perceptually uniform, and gamma must also be taken into account when dithering to lower bit depths (it seems to only significantly matter at < 6bpc though).

Concretely, the 3-bit value 1 maps to 8-bit values from 19 to 54, so the value 1 in 3-bit is thus representative of this whole range. Actually it's a bit more nuanced since the range of light values represented by a 3bpc level of 1 is [0.5, 1.5), which corresponds to a 8bpc level of [~18.2, 54.6). So in terms of discrete actual 8bpc samples it's 19-54, but you can see that the boundaries do not overlap exactly. Moreover in terms of light levels, a 3bpc value of 1 represents the midpoint of the interval (0.07, 0.21) at 1/7 light output, which cannot be represented exactly at a gradient of 1/255 as it would map to 36.4/255. This is where the error in upconversion comes from, simply put the discrete steps in each bit depth (1/7 vs 1/255) don't perfectly align.

The takeaway is that if you try to implement error-diffusion dithering with a get_closest_color of

def get_closest_color(v, nc = 8):
    """
    input and output are normalized to range [0, 1]. Return closest quantized value
    for a bit depth that results in nc unique samples (i.e. range [0, nc -1])
    """
    return round(v * (nc - 1)) / (nc - 1)

then the calculated errors are only correct assuming the values are kept at a bit depth of nc samples. If you know for sure that you will be converting the samples back to a higher bit-depth later, then you should take that into account when calculating error:

def get_closest_color(v, nc = 8):
    """
    input and output are normalized to range [0, 1]. Return closest quantized value
    for a bit depth that results in nc unique samples (i.e. range [0, nc -1]),
    assuming the resulting quantized sample will later be upconverted to 8bpc [0, 255].
    """
    quantized = np.round(old_val * (nc - 1)) / (nc - 1)
    return np.round(quantized * 255.0)/255.0

Some examples of floyd-steinberg dithering directly use a hardcoded palette with values expressed in standard 8bpc, and that's an analogous way to get around the issue.

Answer 2

Yes, it is the correct behavior.

And the reason why repeating the same algorithm does not change the image is:

This algorithm simplify the image so it try to change the color of the pixel to the nearest color from a small colors set. So if the pixel color becomes equal to a Cor from the small colors set, then it is the nearest color to itself and will remain unchanged.