Color quantization/C: Difference between revisions
Content added Content deleted
(Created page with "This is a complete program that takes a PPM P6 image and a number, then writes out the image reduced to the number of colors to out.ppm. <lang c>#include <stdio.h> #include <stdl...") |
m (Replace <lang> by <syntaxhighlight> block so it's better formatted) |
||
(3 intermediate revisions by 2 users not shown) | |||
Line 1: | Line 1: | ||
This is a complete program that takes a PPM P6 image and a number, then writes out the image reduced to the number of colors to out.ppm. |
[[file:Quantum_frog_dithered.png|200px|thumb]]This is a complete program that takes a PPM P6 image and a number, then writes out the image reduced to the number of colors to out.ppm. There is optional dithering, too, which doesn't make a whole lot of difference with say 64 colors or more. And with low colors, the quantization did such a good job of picking average colors that it actually hurts the dithering process. |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include <fcntl.h> |
#include <fcntl.h> |
||
Line 7: | Line 7: | ||
#include <string.h> |
#include <string.h> |
||
#include <stdint.h> |
#include <stdint.h> |
||
#include <math.h> |
|||
typedef struct { |
typedef struct { |
||
Line 15: | Line 16: | ||
int write_ppm(image im, char *fn) |
int write_ppm(image im, char *fn) |
||
{ |
{ |
||
FILE *fp = fopen(fn, " |
FILE *fp = fopen(fn, "wb"); |
||
if (!fp) return 0; |
if (!fp) return 0; |
||
fprintf(fp, "P6\n%d %d\n255\n", im->w, im->h); |
fprintf(fp, "P6\n%d %d\n255\n", im->w, im->h); |
||
Line 46: | Line 47: | ||
image read_ppm(char *fn) |
image read_ppm(char *fn) |
||
{ |
{ |
||
FILE *fp = fopen(fn, " |
FILE *fp = fopen(fn, "rb"); |
||
int w, h, maxval; |
int w, h, maxval; |
||
image im = 0; |
image im = 0; |
||
Line 68: | Line 69: | ||
#define ON_INHEAP 1 |
#define ON_INHEAP 1 |
||
typedef struct oct_node_t oct_node_t, *oct_node; |
typedef struct oct_node_t oct_node_t, *oct_node; |
||
struct oct_node_t{ |
struct oct_node_t{ |
||
/* sum of all |
int64_t r, g, b; /* sum of all child node colors */ |
||
uint64_t r, g, b; |
|||
int count, heap_idx; |
int count, heap_idx; |
||
oct_node kids[8], parent; |
|||
unsigned char n_kids, kid_idx, flags, depth; |
unsigned char n_kids, kid_idx, flags, depth; |
||
oct_node kids[8], parent; |
|||
}; |
}; |
||
Line 82: | Line 83: | ||
} node_heap; |
} node_heap; |
||
/* cmp function that decides the ordering in the heap. This is how we determine |
|||
which octree node to fold next, the heart of the algorithm. */ |
|||
inline int cmp_node(oct_node a, oct_node b) |
inline int cmp_node(oct_node a, oct_node b) |
||
{ |
{ |
||
Line 89: | Line 88: | ||
if (a->n_kids > b->n_kids) return 1; |
if (a->n_kids > b->n_kids) return 1; |
||
int ac = a->count |
int ac = a->count >> a->depth; |
||
int bc = b->count |
int bc = b->count >> b->depth; |
||
return ac < bc ? -1 : ac > bc; |
return ac < bc ? -1 : ac > bc; |
||
} |
} |
||
Line 195: | Line 194: | ||
oct_node node_insert(oct_node root, unsigned char *pix) |
oct_node node_insert(oct_node root, unsigned char *pix) |
||
{ |
{ |
||
# define OCT_DEPTH 8 |
|||
/* 8: number of significant bits used for tree. It's probably good enough |
|||
for most images to use a value of 5. This affects how many nodes eventually |
|||
end up in the tree and heap, thus smaller values helps with both speed |
|||
and memory. */ |
|||
unsigned char i, bit, depth = 0; |
unsigned char i, bit, depth = 0; |
||
for (bit = 1 << 7; ++depth < OCT_DEPTH; bit >>= 1) { |
|||
for (bit = 1 << 7; ++depth < 8; bit >>= 1) { |
|||
i = !!(pix[1] & bit) * 4 + !!(pix[0] & bit) * 2 + !!(pix[2] & bit); |
i = !!(pix[1] & bit) * 4 + !!(pix[0] & bit) * 2 + !!(pix[2] & bit); |
||
if (!root->kids[i]) |
if (!root->kids[i]) |
||
Line 231: | Line 225: | ||
} |
} |
||
/* traverse the octree just like construction, but this time we replace the pixel |
|||
color with color stored in the tree node */ |
|||
void color_replace(oct_node root, unsigned char *pix) |
void color_replace(oct_node root, unsigned char *pix) |
||
{ |
{ |
||
Line 248: | Line 240: | ||
} |
} |
||
void |
void error_diffuse(image im, node_heap *h) |
||
{ |
|||
oct_node nearest_color(int *v) { |
|||
int i; |
|||
int diff, max = 100000000; |
|||
oct_node o = 0; |
|||
for (i = 1; i < h->n; i++) { |
|||
diff = 3 * abs(h->buf[i]->r - v[0]) |
|||
+ 5 * abs(h->buf[i]->g - v[1]) |
|||
+ 2 * abs(h->buf[i]->b - v[2]); |
|||
if (diff < max) { |
|||
max = diff; |
|||
o = h->buf[i]; |
|||
} |
|||
} |
|||
return o; |
|||
} |
|||
# define POS(i, j) (3 * ((i) * im->w + (j))) |
|||
int i, j; |
|||
int *npx = calloc(sizeof(int), im->h * im->w * 3), *px; |
|||
int v[3]; |
|||
unsigned char *pix = im->pix; |
|||
oct_node nd; |
|||
#define C10 7 |
|||
#define C01 5 |
|||
#define C11 2 |
|||
#define C00 1 |
|||
#define CTOTAL (C00 + C11 + C10 + C01) |
|||
for (px = npx, i = 0; i < im->h; i++) { |
|||
for (j = 0; j < im->w; j++, pix += 3, px += 3) { |
|||
px[0] = (int)pix[0] * CTOTAL; |
|||
px[1] = (int)pix[1] * CTOTAL; |
|||
px[2] = (int)pix[2] * CTOTAL; |
|||
} |
|||
} |
|||
#define clamp(x, i) if (x[i] > 255) x[i] = 255; if (x[i] < 0) x[i] = 0 |
|||
pix = im->pix; |
|||
for (px = npx, i = 0; i < im->h; i++) { |
|||
for (j = 0; j < im->w; j++, pix += 3, px += 3) { |
|||
px[0] /= CTOTAL; |
|||
px[1] /= CTOTAL; |
|||
px[2] /= CTOTAL; |
|||
clamp(px, 0); clamp(px, 1); clamp(px, 2); |
|||
nd = nearest_color(px); |
|||
v[0] = px[0] - nd->r; |
|||
v[1] = px[1] - nd->g; |
|||
v[2] = px[2] - nd->b; |
|||
pix[0] = nd->r; pix[1] = nd->g; pix[2] = nd->b; |
|||
if (j < im->w - 1) { |
|||
npx[POS(i, j+1) + 0] += v[0] * C10; |
|||
npx[POS(i, j+1) + 1] += v[1] * C10; |
|||
npx[POS(i, j+1) + 2] += v[2] * C10; |
|||
} |
|||
if (i >= im->h - 1) continue; |
|||
npx[POS(i+1, j) + 0] += v[0] * C01; |
|||
npx[POS(i+1, j) + 1] += v[1] * C01; |
|||
npx[POS(i+1, j) + 2] += v[2] * C01; |
|||
if (j < im->w - 1) { |
|||
npx[POS(i+1, j+1) + 0] += v[0] * C11; |
|||
npx[POS(i+1, j+1) + 1] += v[1] * C11; |
|||
npx[POS(i+1, j+1) + 2] += v[2] * C11; |
|||
} |
|||
if (j) { |
|||
npx[POS(i+1, j-1) + 0] += v[0] * C00; |
|||
npx[POS(i+1, j-1) + 1] += v[1] * C00; |
|||
npx[POS(i+1, j-1) + 2] += v[2] * C00; |
|||
} |
|||
} |
|||
} |
|||
free(npx); |
|||
} |
|||
void color_quant(image im, int n_colors, int dither) |
|||
{ |
{ |
||
int i; |
int i; |
||
Line 268: | Line 340: | ||
got->g = got->g / c + .5; |
got->g = got->g / c + .5; |
||
got->b = got->b / c + .5; |
got->b = got->b / c + .5; |
||
printf("%2d | %3llu %3llu %3llu (%d pixels)\n", |
|||
i, got->r, got->g, got->b, got->count); |
|||
} |
} |
||
if (dither) error_diffuse(im, &heap); |
|||
for (i = 0, pix = im->pix; i < im->w * im->h; i++, pix += 3) |
|||
else |
|||
color_replace(root, pix); |
|||
for (i = 0, pix = im->pix; i < im->w * im->h; i++, pix += 3) |
|||
color_replace(root, pix); |
|||
node_free(); |
node_free(); |
||
Line 289: | Line 361: | ||
image im = read_ppm(v[1]); |
image im = read_ppm(v[1]); |
||
color_quant(im, c); |
color_quant(im, c, 0); |
||
write_ppm(im, "out.pnm"); |
write_ppm(im, "out.pnm"); |
||
free(im); |
free(im); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
Latest revision as of 11:49, 27 October 2022
This is a complete program that takes a PPM P6 image and a number, then writes out the image reduced to the number of colors to out.ppm. There is optional dithering, too, which doesn't make a whole lot of difference with say 64 colors or more. And with low colors, the quantization did such a good job of picking average colors that it actually hurts the dithering process.
#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <unistd.h>
#include <ctype.h>
#include <string.h>
#include <stdint.h>
#include <math.h>
typedef struct {
int w, h;
unsigned char *pix;
} image_t, *image;
int write_ppm(image im, char *fn)
{
FILE *fp = fopen(fn, "wb");
if (!fp) return 0;
fprintf(fp, "P6\n%d %d\n255\n", im->w, im->h);
fwrite(im->pix, 1, 3 * im->w * im->h, fp);
fclose(fp);
return 1;
}
image img_new(int w, int h)
{
image im = malloc(sizeof(image_t) + h * w * 3);
im->w = w; im->h = h;
im->pix = (unsigned char *)(im + 1);
return im;
}
int read_num(FILE *f)
{
int n;
while (!fscanf(f, "%d ", &n)) {
if ((n = fgetc(f)) == '#') {
while ((n = fgetc(f)) != '\n')
if (n == EOF) return 0;
} else
return 0;
}
return n;
}
image read_ppm(char *fn)
{
FILE *fp = fopen(fn, "rb");
int w, h, maxval;
image im = 0;
if (!fp) return 0;
if (fgetc(fp) != 'P' || fgetc(fp) != '6' || !isspace(fgetc(fp)))
goto bail;
w = read_num(fp);
h = read_num(fp);
maxval = read_num(fp);
if (!w || !h || !maxval) goto bail;
im = img_new(w, h);
fread(im->pix, 1, 3 * w * h, fp);
bail:
if (fp) fclose(fp);
return im;
}
#define ON_INHEAP 1
typedef struct oct_node_t oct_node_t, *oct_node;
struct oct_node_t{
int64_t r, g, b; /* sum of all child node colors */
int count, heap_idx;
unsigned char n_kids, kid_idx, flags, depth;
oct_node kids[8], parent;
};
typedef struct {
int alloc, n;
oct_node* buf;
} node_heap;
inline int cmp_node(oct_node a, oct_node b)
{
if (a->n_kids < b->n_kids) return -1;
if (a->n_kids > b->n_kids) return 1;
int ac = a->count >> a->depth;
int bc = b->count >> b->depth;
return ac < bc ? -1 : ac > bc;
}
void down_heap(node_heap *h, oct_node p)
{
int n = p->heap_idx, m;
while (1) {
m = n * 2;
if (m >= h->n) break;
if (m + 1 < h->n && cmp_node(h->buf[m], h->buf[m + 1]) > 0) m++;
if (cmp_node(p, h->buf[m]) <= 0) break;
h->buf[n] = h->buf[m];
h->buf[n]->heap_idx = n;
n = m;
}
h->buf[n] = p;
p->heap_idx = n;
}
void up_heap(node_heap *h, oct_node p)
{
int n = p->heap_idx;
oct_node prev;
while (n > 1) {
prev = h->buf[n / 2];
if (cmp_node(p, prev) >= 0) break;
h->buf[n] = prev;
prev->heap_idx = n;
n /= 2;
}
h->buf[n] = p;
p->heap_idx = n;
}
void heap_add(node_heap *h, oct_node p)
{
if ((p->flags & ON_INHEAP)) {
down_heap(h, p);
up_heap(h, p);
return;
}
p->flags |= ON_INHEAP;
if (!h->n) h->n = 1;
if (h->n >= h->alloc) {
while (h->n >= h->alloc) h->alloc += 1024;
h->buf = realloc(h->buf, sizeof(oct_node) * h->alloc);
}
p->heap_idx = h->n;
h->buf[h->n++] = p;
up_heap(h, p);
}
oct_node pop_heap(node_heap *h)
{
if (h->n <= 1) return 0;
oct_node ret = h->buf[1];
h->buf[1] = h->buf[--h->n];
h->buf[h->n] = 0;
h->buf[1]->heap_idx = 1;
down_heap(h, h->buf[1]);
return ret;
}
static oct_node pool = 0;
oct_node node_new(unsigned char idx, unsigned char depth, oct_node p)
{
static int len = 0;
if (len <= 1) {
oct_node p = calloc(sizeof(oct_node_t), 2048);
p->parent = pool;
pool = p;
len = 2047;
}
oct_node x = pool + len--;
x->kid_idx = idx;
x->depth = depth;
x->parent = p;
if (p) p->n_kids++;
return x;
}
void node_free()
{
oct_node p;
while (pool) {
p = pool->parent;
free(pool);
pool = p;
}
}
oct_node node_insert(oct_node root, unsigned char *pix)
{
unsigned char i, bit, depth = 0;
for (bit = 1 << 7; ++depth < 8; bit >>= 1) {
i = !!(pix[1] & bit) * 4 + !!(pix[0] & bit) * 2 + !!(pix[2] & bit);
if (!root->kids[i])
root->kids[i] = node_new(i, depth, root);
root = root->kids[i];
}
root->r += pix[0];
root->g += pix[1];
root->b += pix[2];
root->count++;
return root;
}
oct_node node_fold(oct_node p)
{
if (p->n_kids) abort();
oct_node q = p->parent;
q->count += p->count;
q->r += p->r;
q->g += p->g;
q->b += p->b;
q->n_kids --;
q->kids[p->kid_idx] = 0;
return q;
}
void color_replace(oct_node root, unsigned char *pix)
{
unsigned char i, bit;
for (bit = 1 << 7; bit; bit >>= 1) {
i = !!(pix[1] & bit) * 4 + !!(pix[0] & bit) * 2 + !!(pix[2] & bit);
if (!root->kids[i]) break;
root = root->kids[i];
}
pix[0] = root->r;
pix[1] = root->g;
pix[2] = root->b;
}
void error_diffuse(image im, node_heap *h)
{
oct_node nearest_color(int *v) {
int i;
int diff, max = 100000000;
oct_node o = 0;
for (i = 1; i < h->n; i++) {
diff = 3 * abs(h->buf[i]->r - v[0])
+ 5 * abs(h->buf[i]->g - v[1])
+ 2 * abs(h->buf[i]->b - v[2]);
if (diff < max) {
max = diff;
o = h->buf[i];
}
}
return o;
}
# define POS(i, j) (3 * ((i) * im->w + (j)))
int i, j;
int *npx = calloc(sizeof(int), im->h * im->w * 3), *px;
int v[3];
unsigned char *pix = im->pix;
oct_node nd;
#define C10 7
#define C01 5
#define C11 2
#define C00 1
#define CTOTAL (C00 + C11 + C10 + C01)
for (px = npx, i = 0; i < im->h; i++) {
for (j = 0; j < im->w; j++, pix += 3, px += 3) {
px[0] = (int)pix[0] * CTOTAL;
px[1] = (int)pix[1] * CTOTAL;
px[2] = (int)pix[2] * CTOTAL;
}
}
#define clamp(x, i) if (x[i] > 255) x[i] = 255; if (x[i] < 0) x[i] = 0
pix = im->pix;
for (px = npx, i = 0; i < im->h; i++) {
for (j = 0; j < im->w; j++, pix += 3, px += 3) {
px[0] /= CTOTAL;
px[1] /= CTOTAL;
px[2] /= CTOTAL;
clamp(px, 0); clamp(px, 1); clamp(px, 2);
nd = nearest_color(px);
v[0] = px[0] - nd->r;
v[1] = px[1] - nd->g;
v[2] = px[2] - nd->b;
pix[0] = nd->r; pix[1] = nd->g; pix[2] = nd->b;
if (j < im->w - 1) {
npx[POS(i, j+1) + 0] += v[0] * C10;
npx[POS(i, j+1) + 1] += v[1] * C10;
npx[POS(i, j+1) + 2] += v[2] * C10;
}
if (i >= im->h - 1) continue;
npx[POS(i+1, j) + 0] += v[0] * C01;
npx[POS(i+1, j) + 1] += v[1] * C01;
npx[POS(i+1, j) + 2] += v[2] * C01;
if (j < im->w - 1) {
npx[POS(i+1, j+1) + 0] += v[0] * C11;
npx[POS(i+1, j+1) + 1] += v[1] * C11;
npx[POS(i+1, j+1) + 2] += v[2] * C11;
}
if (j) {
npx[POS(i+1, j-1) + 0] += v[0] * C00;
npx[POS(i+1, j-1) + 1] += v[1] * C00;
npx[POS(i+1, j-1) + 2] += v[2] * C00;
}
}
}
free(npx);
}
void color_quant(image im, int n_colors, int dither)
{
int i;
unsigned char *pix = im->pix;
node_heap heap = { 0, 0, 0 };
oct_node root = node_new(0, 0, 0), got;
for (i = 0; i < im->w * im->h; i++, pix += 3)
heap_add(&heap, node_insert(root, pix));
while (heap.n > n_colors + 1)
heap_add(&heap, node_fold(pop_heap(&heap)));
double c;
for (i = 1; i < heap.n; i++) {
got = heap.buf[i];
c = got->count;
got->r = got->r / c + .5;
got->g = got->g / c + .5;
got->b = got->b / c + .5;
}
if (dither) error_diffuse(im, &heap);
else
for (i = 0, pix = im->pix; i < im->w * im->h; i++, pix += 3)
color_replace(root, pix);
node_free();
free(heap.buf);
}
int main(int c, char *v[])
{
if (c < 3) {
fprintf(stderr, "usage: %s ppm_file n_colors\n", v[0]);
return 0;
}
c = atoi(v[2]) ? : 16; /* GCC extension */
image im = read_ppm(v[1]);
color_quant(im, c, 0);
write_ppm(im, "out.pnm");
free(im);
return 0;
}