PCE 16 colors conversion using kmeans++ and CIE Lab

CMakeLists.txt

cmake_minimum_required (VERSION 3.12)

project(palette C)

add_library(stb_image STATIC stb_image.c stb_image.h stb_image_write.h)

add_executable(kmeans kmeans.c)
target_compile_options(kmeans PUBLIC -Wall -Wshadow -Wextra)
target_link_libraries(kmeans stb_image m)

kmeans.c

#include <errno.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <float.h>

#include "stb_image.h"
#include "stb_image_write.h"

static uint8_t pce_palette[512*3] = {
	0, 0, 0,
    0, 0, 22, 
	0, 0, 52, 
	0, 0, 72, 
	0, 0, 103, 
	0, 0, 128, 
	0, 0, 167, 
	0, 0, 207, 
	27, 0, 3, 
	31, 3, 31, 
	27, 0, 53, 
	31, 1, 82, 
	27, 0, 105, 
	31, 0, 136, 
	25, 0, 163, 
	27, 0, 207, 
	62, 8, 11, 
	57, 3, 34, 
	61, 6, 64, 
	57, 1, 85, 
	61, 4, 114, 
	57, 0, 136, 
	60, 2, 170, 
	54, 0, 201, 
	87, 8, 13, 
	83, 3, 36, 
	87, 6, 64, 
	83, 1, 87, 
	87, 5, 115, 
	83, 0, 138, 
	86, 3, 170, 
	81, 0, 197, 
	114, 7, 15, 
	118, 11, 44, 
	114, 6, 67, 
	118, 9, 98, 
	114, 4, 119, 
	118, 7, 147, 
	114, 2, 169, 
	117, 5, 203, 
	152, 13, 27, 
	148, 8, 48, 
	152, 11, 78, 
	148, 7, 99, 
	152, 10, 129, 
	148, 5, 150, 
	151, 8, 181, 
	147, 3, 204, 
	185, 8, 31, 
	181, 3, 51, 
	185, 7, 81, 
	181, 2, 103, 
	185, 5, 134, 
	180, 0, 155, 
	185, 4, 183, 
	180, 0, 206, 
	225, 0, 35, 
	229, 2, 66, 
	225, 0, 87, 
	229, 0, 117, 
	225, 0, 138, 
	229, 0, 168, 
	225, 0, 189, 
	229, 0, 220, 
	10, 38, 10, 
	6, 34, 30, 
	10, 37, 61, 
	6, 32, 81, 
	10, 35, 112, 
	5, 30, 135, 
	8, 32, 171, 
	1, 24, 206, 
	37, 38, 10, 
	33, 33, 33, 
	37, 36, 64, 
	33, 31, 84, 
	37, 35, 113, 
	33, 30, 135, 
	36, 32, 169, 
	30, 25, 200, 
	64, 37, 13, 
	68, 40, 44, 
	64, 35, 64, 
	68, 39, 95, 
	64, 34, 115, 
	68, 37, 146, 
	64, 32, 169, 
	67, 34, 205, 
	99, 45, 23, 
	95, 41, 44, 
	98, 44, 75, 
	94, 39, 97, 
	98, 42, 126, 
	94, 37, 147, 
	98, 40, 177, 
	93, 35, 203, 
	124, 45, 26, 
	120, 41, 46, 
	124, 44, 77, 
	120, 39, 97, 
	124, 42, 128, 
	120, 38, 148, 
	124, 41, 179, 
	119, 35, 202, 
	151, 45, 28, 
	155, 48, 56, 
	151, 44, 77, 
	155, 47, 108, 
	151, 41, 130, 
	155, 45, 160, 
	151, 40, 180, 
	155, 43, 211, 
	189, 51, 37, 
	185, 46, 60, 
	189, 49, 89, 
	185, 44, 111, 
	189, 48, 140, 
	185, 43, 162, 
	189, 46, 191, 
	184, 41, 214, 
	222, 46, 42, 
	218, 41, 64, 
	222, 44, 93, 
	218, 40, 113, 
	222, 43, 144, 
	218, 38, 167, 
	222, 41, 196, 
	218, 36, 216, 
	13, 68, 10, 
	17, 71, 40, 
	13, 66, 63, 
	17, 69, 92, 
	13, 65, 113, 
	17, 68, 143, 
	12, 62, 169, 
	14, 63, 207, 
	47, 76, 20, 
	43, 71, 43, 
	47, 75, 71, 
	43, 70, 94, 
	47, 73, 122, 
	43, 68, 145, 
	38, 63, 168, 
	41, 65, 204, 
	74, 75, 22, 
	70, 71, 43, 
	74, 74, 74, 
	70, 69, 96, 
	74, 72, 126, 
	70, 67, 146, 
	74, 71, 176, 
	69, 65, 202, 
	101, 74, 25, 
	106, 78, 54, 
	101, 73, 76, 
	106, 76, 105, 
	101, 71, 127, 
	106, 75, 156, 
	101, 70, 179, 
	105, 73, 210, 
	136, 83, 35, 
	132, 78, 56, 
	136, 82, 85, 
	131, 76, 107, 
	135, 79, 138, 
	131, 75, 159, 
	127, 70, 180, 
	131, 73, 210, 
	162, 83, 36, 
	157, 78, 59, 
	161, 82, 87, 
	157, 77, 110, 
	161, 80, 138, 
	157, 75, 161, 
	161, 79, 189, 
	157, 73, 212, 
	188, 83, 38, 
	192, 86, 69, 
	188, 81, 89, 
	193, 84, 118, 
	188, 79, 141, 
	192, 82, 172, 
	188, 78, 192, 
	192, 81, 221, 
	227, 89, 47, 
	222, 84, 70, 
	226, 87, 101, 
	222, 82, 122, 
	218, 77, 144, 
	222, 80, 173, 
	218, 75, 195, 
	222, 79, 224, 
	25, 105, 20, 
	20, 100, 42, 
	24, 104, 71, 
	20, 99, 93, 
	24, 102, 122, 
	20, 97, 145, 
	24, 100, 176, 
	18, 94, 203, 
	50, 106, 22, 
	54, 109, 51, 
	50, 104, 73, 
	54, 107, 104, 
	50, 102, 125, 
	54, 106, 154, 
	50, 101, 176, 
	53, 103, 210, 
	76, 105, 25, 
	80, 109, 53, 
	76, 104, 75, 
	80, 107, 104, 
	76, 102, 126, 
	80, 106, 155, 
	76, 101, 178, 
	79, 104, 209, 
	111, 113, 35, 
	107, 108, 55, 
	112, 112, 84, 
	107, 107, 107, 
	111, 110, 137, 
	107, 105, 158, 
	111, 108, 187, 
	107, 103, 209, 
	139, 112, 36, 
	143, 115, 67, 
	139, 111, 87, 
	143, 114, 118, 
	139, 109, 138, 
	143, 112, 168, 
	139, 108, 189, 
	142, 111, 220, 
	165, 112, 38, 
	169, 116, 68, 
	165, 111, 89, 
	169, 114, 118, 
	164, 109, 140, 
	168, 112, 171, 
	164, 108, 192, 
	168, 111, 221, 
	199, 121, 46, 
	195, 116, 69, 
	198, 119, 100, 
	194, 114, 120, 
	198, 117, 151, 
	194, 113, 171, 
	198, 116, 202, 
	194, 111, 222, 
	226, 121, 45, 
	230, 124, 79, 
	225, 118, 102, 
	230, 122, 130, 
	226, 117, 151, 
	229, 120, 182, 
	225, 115, 204, 
	229, 119, 233, 
	27, 135, 21, 
	31, 139, 50, 
	26, 134, 73, 
	30, 137, 103, 
	26, 132, 124, 
	30, 136, 153, 
	26, 130, 175, 
	29, 133, 209, 
	62, 143, 32, 
	58, 138, 53, 
	61, 141, 83, 
	57, 137, 104, 
	61, 140, 134, 
	57, 135, 155, 
	61, 138, 186, 
	57, 133, 209, 
	87, 143, 34, 
	92, 147, 63, 
	87, 142, 84, 
	91, 145, 114, 
	87, 140, 137, 
	91, 143, 166, 
	87, 139, 186, 
	91, 142, 217, 
	113, 143, 35, 
	117, 146, 66, 
	113, 142, 86, 
	117, 145, 117, 
	113, 140, 137, 
	117, 143, 167, 
	113, 139, 188, 
	117, 142, 219, 
	149, 151, 45, 
	144, 146, 67, 
	149, 149, 96, 
	145, 144, 117, 
	149, 147, 148, 
	144, 142, 170, 
	148, 146, 199, 
	144, 141, 220, 
	176, 150, 45, 
	180, 153, 77, 
	176, 148, 99, 
	180, 152, 128, 
	176, 147, 150, 
	180, 150, 179, 
	176, 145, 201, 
	180, 149, 230, 
	203, 151, 45, 
	206, 153, 79, 
	202, 148, 101, 
	206, 152, 130, 
	202, 147, 151, 
	206, 150, 181, 
	201, 145, 204, 
	206, 148, 233, 
	237, 160, 51, 
	232, 154, 78, 
	236, 157, 110, 
	231, 152, 133, 
	236, 155, 161, 
	231, 150, 184, 
	236, 154, 212, 
	231, 149, 234, 
	34, 175, 31, 
	38, 179, 60, 
	34, 174, 82, 
	38, 177, 111, 
	34, 172, 133, 
	38, 176, 162, 
	34, 170, 184, 
	29, 165, 208, 
	64, 173, 33, 
	68, 176, 62, 
	64, 172, 83, 
	68, 175, 114, 
	63, 170, 136, 
	67, 173, 165, 
	63, 168, 186, 
	67, 171, 216, 
	99, 181, 42, 
	95, 175, 65, 
	99, 179, 94, 
	94, 174, 116, 
	99, 177, 145, 
	94, 172, 167, 
	99, 176, 196, 
	94, 171, 219, 
	125, 181, 45, 
	129, 184, 75, 
	125, 179, 96, 
	129, 183, 125, 
	124, 178, 147, 
	120, 173, 170, 
	124, 176, 199, 
	120, 171, 219, 
	151, 181, 44, 
	155, 184, 76, 
	150, 179, 98, 
	155, 183, 127, 
	150, 178, 149, 
	154, 181, 178, 
	150, 176, 200, 
	154, 180, 229, 
	187, 189, 52, 
	182, 183, 78, 
	186, 187, 109, 
	182, 182, 129, 
	186, 185, 158, 
	182, 180, 181, 
	185, 183, 211, 
	181, 179, 232, 
	215, 189, 51, 
	218, 191, 86, 
	213, 186, 109, 
	217, 189, 140, 
	213, 184, 161, 
	209, 179, 183, 
	213, 183, 212, 
	209, 178, 234, 
	242, 191, 47, 
	244, 192, 86, 
	239, 186, 112, 
	243, 189, 142, 
	239, 185, 163, 
	243, 188, 192, 
	239, 183, 214, 
	243, 186, 245, 
	38, 218, 38, 
	34, 213, 59, 
	38, 216, 88, 
	34, 211, 110, 
	38, 214, 141, 
	34, 210, 162, 
	38, 213, 191, 
	34, 208, 213, 
	71, 213, 41, 
	67, 208, 64, 
	71, 211, 92, 
	67, 206, 115, 
	71, 210, 143, 
	67, 205, 166, 
	71, 208, 194, 
	67, 203, 217, 
	101, 211, 44, 
	105, 214, 74, 
	101, 209, 95, 
	105, 213, 124, 
	101, 207, 147, 
	105, 211, 177, 
	100, 206, 198, 
	105, 209, 227, 
	136, 219, 52, 
	132, 213, 75, 
	136, 216, 106, 
	132, 212, 127, 
	136, 215, 157, 
	132, 210, 178, 
	136, 213, 208, 
	132, 209, 229, 
	163, 219, 52, 
	158, 214, 77, 
	162, 217, 108, 
	158, 212, 128, 
	162, 216, 158, 
	157, 210, 180, 
	161, 213, 211, 
	157, 209, 232, 
	189, 220, 50, 
	192, 222, 85, 
	188, 217, 108, 
	192, 220, 139, 
	188, 215, 160, 
	191, 219, 190, 
	187, 214, 211, 
	191, 217, 241, 
	226, 229, 55, 
	220, 222, 85, 
	223, 224, 119, 
	219, 219, 141, 
	223, 223, 170, 
	219, 218, 191, 
	223, 221, 222, 
	218, 216, 244, 
	255, 230, 49, 
	248, 222, 84, 
	251, 224, 119, 
	246, 219, 142, 
	250, 222, 173, 
	246, 217, 194, 
	250, 220, 224, 
	246, 216, 245, 
	27, 255, 35, 
	31, 255, 66, 
	27, 255, 86, 
	31, 255, 117, 
	27, 254, 137, 
	31, 255, 168, 
	27, 252, 188, 
	30, 255, 219, 
	75, 255, 49, 
	71, 250, 71, 
	75, 254, 100, 
	71, 249, 121, 
	75, 252, 152, 
	71, 247, 174, 
	75, 251, 203, 
	71, 246, 224, 
	109, 251, 51, 
	104, 246, 74, 
	108, 249, 105, 
	104, 244, 125, 
	108, 247, 156, 
	104, 243, 176, 
	108, 246, 207, 
	104, 241, 227, 
	139, 249, 51, 
	142, 252, 85, 
	138, 247, 107, 
	142, 250, 136, 
	138, 245, 157, 
	142, 248, 188, 
	138, 243, 210, 
	142, 247, 239, 
	175, 255, 58, 
	170, 251, 85, 
	173, 254, 116, 
	169, 249, 139, 
	173, 253, 168, 
	169, 248, 190, 
	173, 251, 219, 
	169, 246, 241, 
	202, 255, 54, 
	196, 252, 85, 
	199, 255, 119, 
	195, 250, 141, 
	199, 253, 169, 
	195, 248, 190, 
	199, 251, 221, 
	194, 246, 244, 
	229, 255, 48, 
	231, 255, 91, 
	225, 255, 118, 
	229, 255, 150, 
	225, 253, 172, 
	229, 255, 201, 
	224, 251, 223, 
	229, 255, 252, 
	255, 255, 47, 
	255, 255, 88, 
	255, 255, 126, 
	255, 255, 152, 
	255, 255, 182, 
	255, 255, 203, 
	255, 255, 232, 
	255, 254, 255
};

void rgb2xyz(uint8_t rgb[3], float xyz[3]) {
    float r = rgb[0] / 255.f;
    r = ((r > 0.04045f) ? pow((r + 0.055f) / 1.055f, 2.4f) : (r / 12.92f)) * 100.f;

    float g = rgb[1] / 255.f;
    g = ((g > 0.04045f) ? pow((g + 0.055f) / 1.055f, 2.4f) : (g / 12.92f)) * 100.f;

    float b = rgb[2] / 255.f;
    b = ((b > 0.04045f) ? pow((b + 0.055f) / 1.055f, 2.4f) : (b / 12.92f)) * 100.f;

    xyz[0] = 0.4124f * r + 0.3576f * g + 0.1805f * b;
    xyz[1] = 0.2126f * r + 0.7152f * g + 0.0722f * b;
    xyz[2] = 0.0193f * r + 0.1192f * g + 0.9505f * b;
}

void xyz2rgb(float xyz[3], uint8_t rgb[3]) {
    float x = xyz[0] / 100.f;
    float y = xyz[1] / 100.f;
    float z = xyz[2] / 100.f;
    
    float r =  3.2406f * x - 1.5372f * y - 0.4986f * z;
    float g = -0.9689f * x + 1.8758f * y + 0.0415f * z;
    float b =  0.0557f * x - 0.2040f * y + 1.0570f * z;

    rgb[0] = ((r > 0.0031308f) ? (1.055f * pow(r, 1.f/2.4f) - 0.055f) : (12.92f * r)) * 255.f;
    rgb[1] = ((g > 0.0031308f) ? (1.055f * pow(g, 1.f/2.4f) - 0.055f) : (12.92f * g)) * 255.f;
    rgb[2] = ((b > 0.0031308f) ? (1.055f * pow(b, 1.f/2.4f) - 0.055f) : (12.92f * b)) * 255.f;
}

static const float xyz_whitepoint[3] = {
    96.720f, 100.f, 81.427f
};

void xyz2lab(float xyz[3], float lab[3]) {
    float x = xyz[0] / xyz_whitepoint[0];
    float y = xyz[1] / xyz_whitepoint[1];
    float z = xyz[2] / xyz_whitepoint[2];

    x = (x > 0.008856f) ? cbrt(x) : ((7.787f * x) + 16.f/116.f);
    y = (y > 0.008856f) ? cbrt(y) : ((7.787f * y) + 16.f/116.f);
    z = (z > 0.008856f) ? cbrt(z) : ((7.787f * z) + 16.f/116.f);

    lab[0] = (116.f * y) - 16.f;
    lab[1] = 500.f * (x - y);
    lab[2] = 200.f * (y - z);
}

void lab2xyz(float lab[3], float xyz[3]) {
    float y = (lab[0] + 16.f) / 116.f;
    float x = (lab[1] / 500.f) + y;
    float z = y - (lab[2] / 200.f);

    float x3 = x*x*x;
    xyz[0] = ((x3 > 0.008856f) ? x3 : ((x - 16.f/116.f) / 7.787)) * xyz_whitepoint[0];
    float y3 = y*y*y;
    xyz[1] = ((y3 > 0.008856f) ? y3 : ((y - 16.f/116.f) / 7.787)) * xyz_whitepoint[1];
    float z3 = z*z*z;
    xyz[2] = ((z3 > 0.008856f) ? z3 : ((z - 16.f/116.f) / 7.787)) * xyz_whitepoint[2];
}

void rgb2lab(uint8_t rgb[3], float lab[3]) {
    float xyz[3];
    rgb2xyz(rgb, xyz);
    xyz2lab(xyz, lab);
}

void lab2rgb(float lab[3], uint8_t rgb[3]) {
    float xyz[3];
    lab2xyz(lab, xyz);
    xyz2rgb(xyz, rgb);
}

uint16_t rand_range(uint16_t n) {
    int32_t limit = RAND_MAX - (RAND_MAX % n);
    int32_t r;
    do {
        r = random();
    } while(r >= limit);
    return r % n;
}

float distance(float *u, float *v) {
    return ((u[0]-v[0])*(u[0]-v[0])) + ((u[1]-v[1])*(u[1]-v[1])) + ((u[2]-v[2])*(u[2]-v[2]));
}

void kmeans_init(float *lut, uint16_t *index, int index_count, uint16_t *centroid, int centroid_count) {
    int i, j, k;

    centroid[0] = index[rand_range(index_count)];

    for(j=1; j<centroid_count; j++) {
        float d_max = 0.f;
        int farthest = -1;

        for(i=0; i<index_count; i++) {
            float d_min = FLT_MAX;
            int p0 = index[i];
            for(k=0; k<j; k++) {
                float d;
                int p1 = centroid[k];
                d = distance(&lut[p0*3], &lut[p1*3]);
                if(d < d_min) {
                    d_min = d;
                }
            }
            if(d_min > d_max){
                d_max = d_min;
                farthest = p0;
            }
        }
        centroid[j] = farthest;
    }
}

void kmeans_update(float *lut, uint16_t *index, int index_count, uint16_t *centroid, int centroid_count) {
    float tmp[16*3];
    float sse[2], diff;
    int32_t count[16];
    int i, j, iteration;

    for(iteration=0; iteration<16; iteration++) {
        for(i=0; i<centroid_count; i++) {
            tmp[3*i+0] = tmp[3*i+1] = tmp[3*i+2] = 0.f;
            count[i] = 0;
        }
        for(i=0; i<index_count; i++) {
            float d_min = FLT_MAX;
            int k = 0;
            uint16_t p0 = index[i];
            for(j=0; j<16; j++) {
                uint16_t p1 = centroid[j];
                float d = distance(&lut[p0*3], &lut[p1*3]);
                if(d < d_min) {
                    k = j;
                    d_min = d;
                }
            }
            
            float l = count[k] * tmp[k*3+0];
            float a = count[k] * tmp[k*3+1];
            float b = count[k] * tmp[k*3+2];
            count[k]++;
            tmp[3*k+0] = (lut[p0*3+0] + l) / count[k];
            tmp[3*k+1] = (lut[p0*3+1] + a) / count[k];
            tmp[3*k+2] = (lut[p0*3+2] + b) / count[k];
        }
        for(i=0; i<16; i++) {
            float d_min = FLT_MAX;
            int k = 0;
            for(j=0; j<512; j++) {
                float d = distance(&tmp[3*i], &lut[3*j]);
                if(d < d_min) {
                    d_min = d;
                    k = j;
                }
            }
            centroid[i] = k;
        }

        sse[1] = sse[0];
        sse[0] = 0.f;
        for(i=0; i<index_count; i++) {
            float d_min = FLT_MAX;
            for(j=0; j<16; j++) {
                float d = distance(&lut[centroid[j]*3], &lut[index[i]*3]);
                if(d < d_min) {
                    d_min = d;
                }
            }
            sse[0] += d_min / (100.f*100.f);
        }
        
        diff = sse[0] - sse[1];
        if((diff <= 0.f) && (diff > -0.1f)) {
            return;
        }
    }
}

int main(int argc, char **argv) {
    int i, j, width, height;
    uint8_t *input, *output;
    uint16_t *index;
    uint16_t centroid[16];
    float lab_lut[512*3], *lab;
    
    int ret = EXIT_FAILURE;

    input = output = NULL;
    index = NULL;
    lab = NULL;
    
    if(argc != 3) {
		fprintf(stderr, "Usage: palette <input> <output>\n input - 24bpp RGB image\n output - output map\n");
		goto err;
	}

	input = stbi_load(argv[1], &width, &height, NULL, 3);
	if(input == NULL) {
		fprintf(stderr, "Failed to read %s\n", argv[1]);
        goto err;
	}

    index = (uint16_t*)malloc(width * height * sizeof(uint16_t));
    if(index == NULL) {
        fprintf(stderr, "failed to allocate index buffer: %s\n", strerror(errno));
        goto err;
    }

    output = (uint8_t*)malloc(3 * width * height);
    if(output == NULL) {
        fprintf(stderr, "failed to allocate output image: %s\n", strerror(errno));
        goto err;
    }

    lab = (float*)malloc(3 * width * height * sizeof(float));
    if(lab == NULL) {
        fprintf(stderr, "failed to allocate lab buffer: %s\n", strerror(errno));
        goto err;
    }

    for(i=0; i<512; i++) {
        rgb2lab(&pce_palette[i*3], &lab_lut[i*3]);
    }

    for(i=0; i<(width*height); i++) {
       rgb2lab(&input[i*3], &lab[i*3]);
    }

    for(i=0; i<(width*height); i++) {
        int k = 0;
        float d_min = FLT_MAX;
        for(j=0; j<512; j++) {
            float d = distance(&lab_lut[j*3], &lab[i*3]);
            if(d < d_min) {
                d_min = d;
                k = j;
            }
        }
        index[i] = k;
    }

    kmeans_init(lab_lut, index, width*height, centroid, 16);
    kmeans_update(lab_lut, index, width*height, centroid, 16);

    for(i=0; i<(width*height); i++) {
        int k = 0;
        float d_min = FLT_MAX;
        for(j=0; j<16; j++) {
            float d = distance(&lab_lut[centroid[j]*3], &lab_lut[index[i]*3]);
            if(d < d_min) {
                d_min = d;
                k = centroid[j];
            }
        }
        output[3*i+0] = pce_palette[3*k+0];
        output[3*i+1] = pce_palette[3*k+1];
        output[3*i+2] = pce_palette[3*k+2];
    }

    stbi_write_png(argv[2], width, height, 3, output, 0);

    ret = EXIT_SUCCESS;
err:
    if(lab) {
        free(lab);
    }
    if(output) {
        free(output);
    }
    if(index) {
        free(index);
    }
    if(input) {
    	stbi_image_free(input);
    }
    return ret;
}

stb_image.c

#define STB_IMAGE_IMPLEMENTATION
#include "stb_image.h"

#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image_write.h"

• https://github.com/nothings/stb/blob/master/stb_image_write.h
• https://github.com/nothings/stb/blob/master/stb_image.h