/*----------------------------------------------------------------------------
Color reduction (quantization) demo code
Chris Giese	<ChrisGiese@SisAndHappy.com>
http://SisAndHappy.com/ChrisGiese

This code converts a .BMP file that uses 24-bit color
into a .BMP file that uses 8-bit color.

I, the copyright holder of this work, hereby release it into the
public domain. This applies worldwide. If this is not legally possible:
I grant any entity the right to use this work for any purpose,
without any conditions, unless such conditions are required by law.

Feb 2014:
- Heavily rewritten; seems to work well now.
- Don't need OPEN()/CLOSE()/READ()/etc. any more
  so they've been replaced with STDIO.H I/O functions

28 Sep 2013:
- Added OPEN()/CLOSE()/READ()/WRITE() macros/functions so SMALL memory
  model works with 16-bit compilers
- Replaced read_leXX()/write_leXX() macros with functions that are
  independent of CPU bit order

6 Mar 2004: initial release
----------------------------------------------------------------------------*/
#include <string.h> /* NULL, memset() */
#include <stdlib.h> /* calloc(), free() */
/* FILE, SEEK_..., EOF, L_tmpnam, printf(), tmpnam(), remove(), rename()
fopen(), fseek(), fread(), fwrite(), fputc(), fclose() */
#include <stdio.h>
#if 0
/* C99 fixed-width types */
#include <stdint.h>
#else
typedef unsigned char	uint8_t;
typedef unsigned short	uint16_t;
typedef unsigned long	uint32_t;
#endif

#if 0
#define	DEBUG(X)	X
#else
#define	DEBUG(X)	/* nothing */
#endif

#if defined(__TURBOC__)
#if __TURBOC__==0x401
#error Sorry, 'huge' is broken in Turbo C++ 3.0. Use TC++ 1.0 or BC++ 3.1
#endif
#include <alloc.h> /* farcalloc(), farfree() */
#define	HUGE		huge
#define	CALLOC(N)	farcalloc(1, N)
#define	FREE(N)		farfree(N)

#elif (defined(__WATCOMC__) && !defined(__386__))
#include <malloc.h> /* _fcalloc(), _ffree() */
#define	HUGE		huge
#define	CALLOC(N)	_fcalloc(1, N)
#define	FREE(N)		_ffree(N)

#elif (defined(__WATCOMC__) && defined(__386__))
#define	HUGE		/* nothing */
#define	CALLOC(N)	calloc(1, N)
#define	FREE(N)		free(N)

#elif defined(__GNUC__)
#define	HUGE		/* ...tracts of land */
#define	CALLOC(N)	calloc(1, N)
#define	FREE(N)		free(N)

#else
#error Sorry, unsupported compiler
#endif

/* 16-bit color value (High Color) */
typedef union
{
	struct
	{
		unsigned b : 5;
		unsigned g : 6;
		unsigned r : 5;
	} c;
	uint16_t i;
} ci_t;

/* Histogram entry.
Only 5 bytes, I hope -- I allocate 65536 of these =320 Kbytes. */
#pragma pack(1)
typedef struct
{
/* number of pixels with this color in input image */
	uint32_t pop;
/* palette index for this color in output image */
	uint8_t pal_index;
} hist_t;

/* A rectangular prism in the RGB colorspace.
The corner closest  to   the origin is RGB=(min_r, min_g, min_b).
The corner farthest from the origin is RGB=(max_r-1, max_g-1, max_b-1). */
typedef struct
{
	short min_r, min_g, min_b, max_r, max_g, max_b;
} box_t;

typedef struct
{
	unsigned long pop;
	unsigned long vol;
	box_t box;
} state_t;
/****************************************************************************/
static void dump_box(box_t *b)
{
	printf("RGB=(%3u,%3u,%3u) to RGB=(%3u,%3u,%3u)",
		b->min_r, b->min_g, b->min_b,
		b->max_r, b->max_g, b->max_b);
}
/****************************************************************************/
static void normalize_box(state_t *s, hist_t HUGE *hist)
{
	unsigned long pop = 0, vol;
	box_t *box = &s->box;
	short r, g, b;
	ci_t ci;

/* shrink box_t so the outermost pixels lie on the edges of the box */
	for(r = box->min_r; r < box->max_r; r++)
		for(g = box->min_g; g < box->max_g; g++)
			for(b = box->min_b; b < box->max_b; b++)
			{
				ci.c.r = r;
				ci.c.g = g;
				ci.c.b = b;
				if(hist[ci.i].pop != 0)
				{
					box->min_r = r;
					goto MIN_R;
				}
			}
MIN_R:
	for(r = box->max_r - 1; r >= box->min_r; r--)
		for(g = box->max_g - 1; g >= box->min_g; g--)
			for(b = box->max_b - 1; b >= box->min_b; b--)
			{
				ci.c.r = r;
				ci.c.g = g;
				ci.c.b = b;
				if(hist[ci.i].pop != 0)
				{
					box->max_r = r + 1;
					goto MAX_R;
				}
			}
MAX_R:
	for(g = box->min_g; g < box->max_g; g++)
		for(r = box->min_r; r < box->max_r; r++)
			for(b = box->min_b; b < box->max_b; b++)
			{
				ci.c.r = r;
				ci.c.g = g;
				ci.c.b = b;
				if(hist[ci.i].pop != 0)
				{
					box->min_g = g;
					goto MIN_G;
				}
			}
MIN_G:
	for(g = box->max_g - 1; g >= box->min_g; g--)
		for(r = box->max_r - 1; r >= box->min_r; r--)
			for(b = box->max_b - 1; b >= box->min_b; b--)
			{
				ci.c.r = r;
				ci.c.g = g;
				ci.c.b = b;
				if(hist[ci.i].pop != 0)
				{
					box->max_g = g + 1;
					goto MAX_G;
				}
			}
MAX_G:
	for(b = box->min_b; b < box->max_b; b++)
		for(g = box->min_g; g < box->max_g; g++)
			for(r = box->min_r; r < box->max_r; r++)
			{
				ci.c.r = r;
				ci.c.g = g;
				ci.c.b = b;
				if(hist[ci.i].pop != 0)
				{
					box->min_b = b;
					goto MIN_B;
				}
			}
MIN_B:
	for(b = box->max_b - 1; b >= box->min_b; b--)
		for(g = box->max_g - 1; g >= box->min_g; g--)
			for(r = box->max_r - 1; r >= box->min_r; r--)
			{
				ci.c.r = r;
				ci.c.g = g;
				ci.c.b = b;
				if(hist[ci.i].pop != 0)
				{
					box->max_b = b + 1;
					goto MAX_B;
				}
			}
MAX_B: ;
/* calculate population */
	for(r = box->min_r; r < box->max_r; r++)
		for(g = box->min_g; g < box->max_g; g++)
			for(b = box->min_b; b < box->max_b; b++)
			{
				ci.c.r = r;
				ci.c.g = g;
				ci.c.b = b;
				if(hist[ci.i].pop != 0)
					pop++;
			}
	s->pop = pop;
/* calculate volume */
	r = box->max_r - box->min_r;
	g = box->max_g - box->min_g;
	b = box->max_b - box->min_b;
#if 1
	vol  = (unsigned long)r * g * b;
#else
/* Not really volume, but I think this is what libjpeg uses.
Output doesn't seem to be any different than using real volume. */
	vol = r * r + g * g + b * b;
#endif
	s->vol =  vol;
}
/****************************************************************************/
#define	NUM_COLORS	255

static void reduce(state_t *s, hist_t HUGE *hist)
{
	unsigned i, j, max, r, g, b;
	box_t *box1, *box2;
	ci_t ci;

/* create a box containing all possible colors */
	s[0].box.min_r = 0;	s[0].box.max_r = 32;
	s[0].box.min_g = 0;	s[0].box.max_g = 64;
	s[0].box.min_b = 0;	s[0].box.max_b = 32;
/* trim it and get volume & population */
	normalize_box(&s[0], hist);
/* reduce to NUM_COLORS colors */
	for(i = 1; i < NUM_COLORS; i++)
	{
		max = 0;
		for(j = 0; j < i; j++)
		{
#if 1
/* find box with greatest volume */
			if(s[max].vol < s[j].vol)
				max = j;
		}
		if(s[max].vol <= 1)
			break;
#else
/* Find box with greatest number of unique colors.
Output doesn't seem to be any different vs. using volume. */
			if(s[max].pop < s[j].pop)
				max = j;
		}
		if(s[max].pop == 0)
{ printf("Trying to split empty box...\n");
			break;
		}
#endif
	DEBUG(	printf("Splitting box with palette index #%u; "
			"creating box #%u\n", max, i);)
		box1 = &s[max].box;
		box2 = &s[i].box;
/* copy to new box (structure copy) */
		*box2 = *box1;
/* chose axis and location of split:
find longest axis of box */
#if 0
		r = (box1->max_r - box1->min_r) * 8;
		g = (box1->max_g - box1->min_g) * 4;
		b = (box1->max_b - box1->min_b) * 8;
#else
/* libjpeg weighting.
This makes the gravel in the test image (AK24.BMP) look better. */
		r = (box1->max_r - box1->min_r) * 8 * 2;
		g = (box1->max_g - box1->min_g) * 4 * 3;
		b = (box1->max_b - box1->min_b) * 8 * 1;
#endif
/* libjpeg: break ties in favor of green, then red, blue last */
		if(g >= r)
			if(g >= b)
			{
/* do split */
				j = (box1->max_g + box1->min_g) / 2;
				box2->min_g = j;
				box1->max_g = j;
			}
			else
BLUE:			{
				j = (box1->max_b + box1->min_b) / 2;
				box2->min_b = j;
				box1->max_b = j;
			}
		else
			if(r >= b)
			{
				j = (box1->max_r + box1->min_r) / 2;
				box2->min_r = j;
				box1->max_r = j;
			}
			else
				goto BLUE;
/* trim the newly-created boxes and get volume and population of each */
		normalize_box(&s[max], hist);
		normalize_box(&s[i], hist);
	}
#if 0
for(i = 0; i < NUM_COLORS; i++)
{
  box1 = &s[i].box;
  for(j = 0; j < NUM_COLORS; j++)
  {
    char mr, mg, mb, min, xr, xg, xb, max;

    if(j == i)
      continue;
    box2 = &s[j].box;
    min = max = 0;

    mr = (box2->min_r >= box1->min_r && box2->min_r < box1->max_r);
    mg = (box2->min_g >= box1->min_g && box2->min_g < box1->max_g);
    mb = (box2->min_b >= box1->min_b && box2->min_b < box1->max_b);
    if(mr && mg && mb)
	min = 1;
/*  xr = (box2->max_r - 1 >= box1->min_r && box2->max_r - 1 < box1->max_r); */
    xr = (box2->max_r >= 1 + box1->min_r && box2->max_r < 1 + box1->max_r);
    xg = (box2->max_g >= 1 + box1->min_g && box2->max_g < 1 + box1->max_g);
    xb = (box2->max_b >= 1 + box1->min_b && box2->max_b < 1 + box1->max_b);
    if(xr && xg && xb)
	max = 1;
    if(min || max)
    {
      printf("Error: color with palette index %u (", j);
      dump_box(box2);
      printf(") overlaps color with palette index %u (", i);
      dump_box(box1);
      printf(") mr=%u, mg=%u, mb=%u, min=%u;\txr=%u, xg=%u, xb=%u, max=%u\n",
	mr, mg, mb, min, xr, xg, xb, max);
    }
  }
}
#endif
	for(i = 0; i < NUM_COLORS; i++)
	{
		unsigned long tr, tg, tb, total;

/* store palette index 'i' in all hist[] for each color in this box
xxx - can this be done faster? */
		box1 = &s[i].box;
	DEBUG(	j = 0;)
		tr = tg = tb = total = 0;
		for(r = box1->min_r; r < box1->max_r; r++)
			for(g = box1->min_g; g < box1->max_g; g++)
				for(b = box1->min_b; b < box1->max_b; b++)
				{
					ci.c.r = r;
					ci.c.g = g;
					ci.c.b = b;
					hist[ci.i].pal_index = i;
				DEBUG(	j++;)
/* calculate sums for weighted average */
					total += hist[ci.i].pop;
					tr += hist[ci.i].pop * r;
					tg += hist[ci.i].pop * g;
					tb += hist[ci.i].pop * b;
				}
	DEBUG(	printf("Palette entry #%3u used for %4u color(s); "
			"from ", i, j);
		dump_box(box1);
		printf("\n");)
/* store weighted average color */
		box1->min_r = tr / total;
		box1->min_g = tg / total;
		box1->min_b = tb / total;
	}
}
/****************************************************************************/
static unsigned long read_le32(const void *buf0)
{
	const unsigned char *buf = buf0;
	unsigned long rv;

	rv = buf[3]; /* MSB */
	rv <<= 8;
	rv |= buf[2];
	rv <<= 8;
	rv |= buf[1];
	rv <<= 8;
	rv |= buf[0]; /* LSB */
	return rv;
}
/****************************************************************************/
static void write_le16(void *buf0, unsigned val)
{
	unsigned char *buf = buf0;

	buf[0] = val;
	val >>= 8;
	buf[1] = val;
}
/****************************************************************************/
static void write_le32(void *buf0, unsigned long val)
{
	unsigned char *buf = buf0;

	buf[0] = val;
	val >>= 8;
	buf[1] = val;
	val >>= 8;
	buf[2] = val;
	val >>= 8;
	buf[3] = val;
}
/****************************************************************************/
#define	BMP_HDR_LEN	54

int main(int arg_c, char *arg_v[])
{
static state_t state[256];
	unsigned wd, ht, ipad, opad, y, x;
	unsigned char buf[BMP_HDR_LEN];
	unsigned long raster_start;
	char out_name[L_tmpnam];
	FILE *ifile, *ofile;
	hist_t HUGE *hist;
	ci_t ci;

	ci.c.r = 31;
	ci.c.g = 0;
	ci.c.b = 7;
	if(ci.i != 0xF807)
	{
		printf("Software error: bit fields in ci_t structure are "
			"not\nwhere they should be\n");
		return 1;
	}
/* check args */
	if(arg_c != 2)
	{
		printf("Converts .BMP files with 24-bit color to 8-bit color\n");
		return 2;
	}
/* open input file */
	if((ifile = fopen(arg_v[1], "rb")) == NULL)
	{
		printf("Error: can't open input file '%s'\n", arg_v[1]);
		return 3;
	}
/* read .BMP file header */
	if(fread(buf, 1, BMP_HDR_LEN, ifile) != BMP_HDR_LEN)
NOT:	{
		printf("Error: file '%s' is not a 24-bit .BMP file\n",
			arg_v[1]);
		fclose(ifile);
		return 4;
	}
/* validate */
	if(buf[0] != 'B' || buf[1] != 'M'	/* signature bytes */
		|| buf[28] != 24		/* color depth */
		|| read_le32(&buf[30]) != 0)	/* compression (0=none) */
			goto NOT;
/* get info */
	wd = (unsigned)read_le32(&buf[18]);
	ht = (unsigned)read_le32(&buf[22]);
	raster_start = read_le32(&buf[10]);
	ipad =  wd & 3; /* <- this works for 24-bit .BMP only */
	opad = -wd & 3; /* <- this works for  8-bit .BMP only */
/* allocate memory for histogram */
	if((hist = CALLOC(65536L * sizeof(hist_t))) == NULL)
	{
		printf("Error: out of memory\n");
		fclose(ifile);
		return 5;
	}
/* Seek to and read input raster.
It's upside-down but that doesn't matter here. */
	fseek(ifile, raster_start, SEEK_SET);
	for(y = 0; y < ht; y++)
	{
		for(x = 0; x < wd; x++)
		{
/* read one 24-bit pixel (note BGR order) */
			if(fread(buf, 1, 3, ifile) != 3)
			{
				printf("Error reading input file '%s'\n",
					arg_v[1]);
				FREE(hist);
				fclose(ifile);
				return 6;
			}
/* reduce color to 16-bit and enter value into histogram */
			ci.c.b = buf[0] >> 3;
			ci.c.g = buf[1] >> 2;
			ci.c.r = buf[2] >> 3;
			hist[ci.i].pop++;
		}
/* skip pad bytes at end of row */
		fseek(ifile, ipad, SEEK_CUR);
	}
/* DO COLOR REDUCTION */
	reduce(state, hist);
/* open temporary output file */
	tmpnam(out_name);
	if((ofile = fopen(out_name, "wb")) == NULL)
	{
		printf("Error: can't open temporary output file '%s'\n",
			out_name);
		FREE(hist);
		fclose(ifile);
		return 3;
	}
/* write header */
	memset(&buf, 0, sizeof(buf));
	buf[0] = 'B';
	buf[1] = 'M';
	write_le32(&buf[2], BMP_HDR_LEN + 1024
		+ (unsigned long)wd * ht); /* file size */
	write_le32(&buf[10], BMP_HDR_LEN + 1024); /* raster offset */
	write_le32(&buf[14], 40);	/* size of hdr2 */
	write_le32(&buf[18], wd);
	write_le32(&buf[22], ht);
	write_le16(&buf[26], 1);	/* number of planes */
	write_le16(&buf[28], 8);	/* depth (bits/pixel) */
	if(fwrite(buf, 1, BMP_HDR_LEN, ofile) != BMP_HDR_LEN)
ERR:	{
		printf("Error writing output file '%s' "
			"(disk full?)\n", out_name);
		fclose(ofile);
		remove(out_name);
		FREE(hist);
		fclose(ifile);
		return 6;
	}
/* write palette */
	buf[3] = 0;
	for(x = 0; x < 256; x++)
	{
		buf[2] = state[x].box.min_r << 3;
		buf[1] = state[x].box.min_g << 2;
		buf[0] = state[x].box.min_b << 3;
		if(fwrite(buf, 1, 4, ofile) != 4)
			goto ERR;
	}
/* go back to start of raster in input file */
	fseek(ifile, raster_start, SEEK_SET);
	for(y = 0; y < ht; y++)
	{
		for(x = 0; x < wd; x++)
		{
/* read one 24-bit pixel (note BGR order) */
			if(fread(buf, 1, 3, ifile) != 3)
			{
				printf("Error reading input file '%s'\n",
					arg_v[1]);
				fclose(ofile);
				remove(out_name);
				FREE(hist);
				fclose(ifile);
				return 6;
			}
/* reduce color to 16-bit */
			ci.c.b = buf[0] >> 3;
			ci.c.g = buf[1] >> 2;
			ci.c.r = buf[2] >> 3;
/* lookup 8-bit palette entry and write to output file */
			if(fputc(hist[ci.i].pal_index, ofile) == EOF)
				goto ERR;
		}
/* skip pad bytes at end of row */
		fseek(ifile, ipad, SEEK_CUR);
		fseek(ofile, opad, SEEK_CUR);
	}
/* all done */
	fclose(ofile);
	FREE(hist);
	fclose(ifile);
/* replace input file with output file */
	if(remove(arg_v[1]))
	{
		printf("Error: input file '%s' is read-only\n", arg_v[1]);
		remove(out_name);
		return 7;
	}
	else
	{
		if(rename(out_name, arg_v[1]))
		{
			printf("Error: can't rename file '%s' to '%s'\n",
				out_name, arg_v[1]);
			return 7;
		}
	}
	return 0;
}