Floyd-Steinberg dithering (C)

From LiteratePrograms
Jump to: navigation, search
Other implementations: C | Java

This is an implementation of Floyd-Steinberg dithering in C for color images. As input we receive an image in 24-bit RGB format and a palette of at most 256 colors, and as output we produce a palettized image where each pixel value is an index into the palette:

<<type declarations>>=
typedef struct {
    unsigned char R, G, B;
} RGBTriple;

typedef struct {
    int size;
    RGBTriple* table;
} RGBPalette;

typedef struct {
    int width, height;
    RGBTriple* pixels;
} RGBImage;

typedef struct {
    int width, height;
    unsigned char* pixels;
} PalettizedImage;

<<FloydSteinbergDither declaration>>=
PalettizedImage FloydSteinbergDither(RGBImage image, RGBPalette palette)

Note that we reserve a full 8 bits for each pixel in the output even when the palette is small — we leave compacting of bits to the caller.

On some platforms, the RGBTriple structure may be padded with extra bits, resulting in wasted space. On such platforms an alternative implementation is to directly use arrays of unsigned characters for pixel data, indexed using macros.

Next we implement the core algorithm. Floyd-Steinberg is based on the concept of error diffusion, where we choose the nearest palette color that we can to the current pixel, and then compute the difference of that color from the original color in each RGB channel. Pieces of this difference are dispersed throughout several adjacent pixels not yet visited:

<<FloydSteinbergDither definition>>=
FloydSteinbergDither declaration
{
    PalettizedImage result;
    set up result image

    loop over all pixels
        RGBTriple* currentPixel = &(image.pixels[x + y*image.width]);
        find nearest color to current pixel
        compute and disperse error
    close loop over all pixels
    return result;
}

There are many ways to find the nearest palette color with varying levels of efficiency and quality. Here we use a trivial algorithm that finds the color with minimum straight-line distance in the color cube from the given color:

<<find nearest color to current pixel>>=
unsigned char index = FindNearestColor(*currentPixel, palette);
result.pixels[x + y*result.width] = index;

<<FindNearestColor definition>>=
unsigned char FindNearestColor(RGBTriple color, RGBPalette palette) {
    int i, distanceSquared, minDistanceSquared, bestIndex = 0;
    minDistanceSquared = 255*255 + 255*255 + 255*255 + 1;
    for (i=0; i<palette.size; i++) {
        int Rdiff = ((int)color.R) - palette.table[i].R;
        int Gdiff = ((int)color.G) - palette.table[i].G;
        int Bdiff = ((int)color.B) - palette.table[i].B;
        distanceSquared = Rdiff*Rdiff + Gdiff*Gdiff + Bdiff*Bdiff;
        if (distanceSquared < minDistanceSquared) {
            minDistanceSquared = distanceSquared;
            bestIndex = i;
        }
    }
    return bestIndex;
}

Note that although the actual straight-line distance is the square root of the quantity measured above, the square of the distance (which is much easier to compute) is minimized by the same color as the distance itself because x2 is strictly increasing on the positive real numbers.

Now, we do the actual dithering. The Floyd-Steinberg dithering matrix looks like this:

X 7/16
3/16 5/16 1/16

For each channel, we determine the difference between the original and palettized value of the pixel labelled "X". We add 7/16 of this difference to the pixel to its right, 5/16 to the pixel below it, and so on. The following macro implements this:

<<compute and disperse error macro simple>>=
#define compute_disperse(channel) \
error = ((int)(currentPixel->channel)) - palette.table[index].channel; \
image.pixels[(x+1) + (y+0)*image.width].channel += (error*7) >> 4; \
image.pixels[(x-1) + (y+1)*image.width].channel += (error*3) >> 4; \
image.pixels[(x+0) + (y+1)*image.width].channel += (error*5) >> 4; \
image.pixels[(x+1) + (y+1)*image.width].channel += (error*1) >> 4;

However, it's possible that enough error could accumulate in a pixel to cause overflow or underflow. Because "wrap around" semantics cause undesirable artifacts, we instead truncate the value to 0 or 255 — this is the most expensive part of the algorithm:

<<plus truncate macro>>=
#define plus_truncate_uchar(a, b) \
    if (((int)(a)) + (b) < 0) \
        (a) = 0; \
    else if (((int)(a)) + (b) > 255) \
        (a) = 255; \
    else \
        (a) += (b);

<<compute and disperse error macro>>=
#define compute_disperse(channel) \
error = ((int)(currentPixel->channel)) - palette.table[index].channel; \
if (x + 1 < image.width) { \
    plus_truncate_uchar(image.pixels[(x+1) + (y+0)*image.width].channel, (error*7) >> 4); \
} \
if (y + 1 < image.height) { \
    if (x - 1 > 0) { \
        plus_truncate_uchar(image.pixels[(x-1) + (y+1)*image.width].channel, (error*3) >> 4); \
    } \
    plus_truncate_uchar(image.pixels[(x+0) + (y+1)*image.width].channel, (error*5) >> 4); \
    if (x + 1 < image.width) { \
        plus_truncate_uchar(image.pixels[(x+1) + (y+1)*image.width].channel, (error*1) >> 4); \
    } \
}

Additionally, we must ensure that we ignore pixels outside of the image area. This might be made more efficient by placing the image inside a buffer with a border around it, but here we just test for it.

This implementation actually modifies the source image in place for efficiency and simplicity, but we could avoid this by instead computing the value passed to FindNearestColor based on the error of adjacent previously visited pixels.

Now we simply invoke the macro once for each channel:

<<compute and disperse error>>=
{
    int error;
    compute_disperse(R);
    compute_disperse(G);
    compute_disperse(B);
}

To loop over the pixels, we proceed in reading order, left to right, top to bottom:

<<loop over all pixels>>=
{
int x, y;
for(y = 0; y < image.height; y++) {
    for(x = 0; x < image.width; x++) {
<<close loop over all pixels>>=
    }
}
}

Another common order used in Floyd-Steinberg dithering alternates between traversing rows left-to-right and right-to-left, which helps avoid some kinds of artifacts. The dithering matrix must be reflected when traversing right-to-left. For simplicity we don't do this here.

Finally, we initialize the result image with the metadata from the original image and a sufficiently large buffer to hold all the pixels' indexes:

<<set up result image>>=
result.width = image.width;
result.height = image.height;
result.pixels = malloc(sizeof(unsigned char) * result.width * result.height);

We use sizeof(unsigned char), which is always 1, as a reminder that if we used a different word size for indexes that a multiplier would be required. The malloc() call requires stdlib.h:

<<header files>>=
#include <stdlib.h>

Finally, we put together the implementation in a source and header file for reuse:

<<floyd_steinberg_dither.h>>=
#ifndef _FLOYD_STEINBERG_DITHER_H_

type declarations
FloydSteinbergDither declaration;

#endif /* #ifndef _FLOYD_STEINBERG_DITHER_H_ */
<<floyd_steinberg_dither.c>>=
header files
#include "floyd_steinberg_dither.h"

plus truncate macro

static
FindNearestColor definition

compute and disperse error macro

FloydSteinbergDither definition

This completes the implementation.

[edit] Sample code

The following sample code reads in an RGB image in a raw image file format, dithers it to a specific optimized 16-color palette found using Photoshop, and writes the output to a raw image file. The original image and dithered result are shown below. The raw input image file for this sample can be downloaded at Image:Amorphophallus titanum USBG small.raw.

Amorphophallus titanum USBG small.png Amorphophallus titanum USBG small dithered.png Amorphophallus titanum palette.png
<<floyd_steinberg_dither_sample.c>>=
#include <stdio.h>
#include <stdlib.h>
#include "floyd_steinberg_dither.h"

int main(int argc, char* argv[]) {
    int x, y;
    RGBImage image;
    RGBPalette palette;
    PalettizedImage result;
    FILE * raw_in, * raw_out;

    image.width = 100;
    image.height = 145;
    image.pixels = malloc(sizeof(RGBTriple) *
                   image.width * image.height);

    raw_in = fopen(argv[1], "rb");
    for(y = 0; y < image.height; y++) {
        for(x = 0; x < image.width; x++) {
            fread(&image.pixels[x + y*image.width],
                  3, 1, raw_in);
        }
    }
    fclose(raw_in);

    palette.size = 16;
    palette.table = malloc(sizeof(RGBTriple) * 16);
    palette.table[0].R = 149;
    palette.table[0].G = 91;
    palette.table[0].B = 110;
    palette.table[1].R = 176;
    palette.table[1].G = 116;
    palette.table[1].B = 137;
    palette.table[2].R = 17;
    palette.table[2].G = 11;
    palette.table[2].B = 15;
    palette.table[3].R = 63;
    palette.table[3].G = 47;
    palette.table[3].B = 69;
    palette.table[4].R = 93;
    palette.table[4].G = 75;
    palette.table[4].B = 112;
    palette.table[5].R = 47;
    palette.table[5].G = 62;
    palette.table[5].B = 24;
    palette.table[6].R = 76;
    palette.table[6].G = 90;
    palette.table[6].B = 55;
    palette.table[7].R = 190;
    palette.table[7].G = 212;
    palette.table[7].B = 115;
    palette.table[8].R = 160;
    palette.table[8].G = 176;
    palette.table[8].B = 87;
    palette.table[9].R = 116;
    palette.table[9].G = 120;
    palette.table[9].B = 87;
    palette.table[10].R = 245;
    palette.table[10].G = 246;
    palette.table[10].B = 225;
    palette.table[11].R = 148;
    palette.table[11].G = 146;
    palette.table[11].B = 130;
    palette.table[12].R = 200;
    palette.table[12].G = 195;
    palette.table[12].B = 180;
    palette.table[13].R = 36;
    palette.table[13].G = 32;
    palette.table[13].B = 27;
    palette.table[14].R = 87;
    palette.table[14].G = 54;
    palette.table[14].B = 45;
    palette.table[15].R = 121;
    palette.table[15].G = 72;
    palette.table[15].B = 72;

    result = FloydSteinbergDither(image, palette);

    raw_out = fopen(argv[2], "wb");
    for(y = 0; y < result.height; y++) {
        for(x = 0; x < result.width; x++) {
            fwrite(&palette.table[result.pixels[x + y*image.width]],
                   3, 1, raw_out);
        }
    }
    fclose(raw_out);

    return 0;
}
Download code
hijacker
hijacker
hijacker
hijacker