arbitrary-palette positional dithering opencl implementation

dither.cl

// Dithering matrices.
// See https://en.wikipedia.org/wiki/Ordered_dithering

__constant float indexMatrix2x2[] = {
    0.25, 0.75,
    1.00, 0.50
};
__constant float indexMatrix4x4[] = {
    0.0625, 0.5625, 0.1875, 0.6875,
    0.8125, 0.3125, 0.9375, 0.4375,
    0.25,   0.75,   0.125,  0.625,
    1.0,    0.5,    0.875,  0.375
};
__constant float indexMatrix8x8[] = {
    0.015625, 0.515625, 0.140625, 0.640625, 0.046875, 0.546875, 0.171875, 0.671875,
    0.765625, 0.265625, 0.890625, 0.390625, 0.796875, 0.296875, 0.921875, 0.421875,
    0.203125, 0.703125, 0.078125, 0.578125, 0.234375, 0.734375, 0.109375, 0.609375,
    0.953125, 0.453125, 0.828125, 0.328125, 0.984375, 0.484375, 0.859375, 0.359375,
    0.0625,   0.5625,   0.1875,   0.6875,   0.03125,  0.53125,  0.15625,  0.65625,
    0.8125,   0.3125,   0.9375,   0.4375,   0.78125,  0.28125,  0.90625,  0.40625,
    0.25,     0.75,     0.125,    0.625,    0.21875,  0.71875,  0.09375,  0.59375,
    1.0,      0.5,      0.875,    0.375,    0.96875,  0.46875,  0.84375,  0.34375
};

// Typically does not need to be larger than 8
#define indexMatrix indexMatrix8x8
#define indexMatrixSize 8

int index(uint2 pos) {
    uint2 ipos = pos % indexMatrixSize;
    return ipos.y * indexMatrixSize + ipos.x;
}
float indexValue(uint2 pos) {
    return indexMatrix[index(pos)];
}

float square(float x) {
    return x * x;
}
// Distance starting at 0 for equal colors, some high value for opposite colors
// No defined max range
float colorDistance(float3 a, float3 b) {
    // See https://en.wikipedia.org/wiki/Luma_%28video%29
    float3 d = b - a;
    // luma corrected distance
    float luma1 = (a.x*0.299 + a.y*0.587 + a.z*0.114);
    float luma2 = (b.x*0.299 + b.y*0.587 + b.z*0.114);
    float dluma = luma2 - luma1;
    return (square(d.x)*0.299 + square(d.y)*0.587 + square(d.z)*0.114)*0.75
         + square(dluma);
}
// Find closest color in palette
float3 closestColor(float3 color, __global const float3 *palette, uint palette_size) {
    float3 result = palette[0];
    for (int i=1; i < palette_size; i++) {
        float3 testColor = palette[i];
        if (colorDistance(testColor, color) < colorDistance(result, color)) {
            result = testColor;
        }
    }
    return result;
}

__kernel void dither(__global float3 *out, __global const float3 *in, uint2 dimensions, __global const float3 *palette, uint palette_size) {
    int id = get_global_id(0);
    uint2 coord = (uint2)(id % dimensions.x, id / dimensions.x);
    float3 color = in[id];
    float3 closest = closestColor(color, palette, palette_size);
    if (closest.x != color.x || closest.y != color.y || closest.z != color.z) {
        float threshold = 1.0 / sqrt((float)palette_size);
        color.x += (indexValue(coord) - 0.5) * threshold;
        color.y += (indexValue(coord) - 0.5) * threshold;
        color.z += (indexValue(coord) - 0.5) * threshold;
        color = clamp(color, (float)0.0, (float)1.0);
        color = closestColor(color, palette, palette_size);
    }
    out[id] = color;
}

@kravohi This is OpenCL, for use on GPUs. You would have to rewrite it in Script-Fu.

By the way, GIMP does have built-in (though I would argue inferior) palette dithering already (Image > Mode > Indexed...).