Skip to content

Instantly share code, notes, and snippets.

@BlockoS

BlockoS/CMakeLists.txt

Last active July 31, 2020 20:59
Show Gist options
  • Save BlockoS/7de1abdb8a979f632b9d805041f29192 to your computer and use it in GitHub Desktop.
Save BlockoS/7de1abdb8a979f632b9d805041f29192 to your computer and use it in GitHub Desktop.
Download ZIP
PCE 16 colors conversion using kmeans++ and CIE Lab
Raw
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)
Raw
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;
}
Raw
stb_image.c
#define STB_IMAGE_IMPLEMENTATION
#include "stb_image.h"
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image_write.h"
@BlockoS
Copy link
Author

BlockoS commented Jul 31, 2020

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment