/* dither.c -- convert 16- or 24-bit bmp to 8 or fewer bits (Floyd-Steinberg) */
/* see end-of-file for (c) */
/* This was written for Microsoft 32-bit-integer compilers--it
 * should compile with
 *	cl gif2bmp.c
 * to produce the executable gif2bmp.exe
 *
 * The Microsoft oddities are:
 *   The funny include files -- they may be omitted on compilers
 *   that can not find them (stdio.h better be there, however)
 * #pragma pack(1) -- This forces the BITMAPFILEHEADER structure
 *   to be the same length in memory as in files (otherwise, the structure
 *   is padded out to be a multiple of 8? bytes long
 * Note that Microsoft uses all-caps for typedefs.  Ugh.
 */
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <string.h>
#include <math.h>
typedef unsigned char byte;
typedef unsigned short WORD;
typedef unsigned long DWORD;

#define SIZE_BITMAPFILEHEADER 14
#define SIZE_BITMAPINFOHEADER 40
#define BI_BITFIELDS 3

double gamma= 1.85;
int outmode= 0;
#define OUT16	1	// 16-color output
#define OUT2    2	// 2-color output
int reversing= 1;
int (*nearrtn)(int,int,int);
int bpp= 8;
int threshhold;
int startnoisy= 1;

unsigned char winmap[]= {
  0,  0,  0,0,   0,  0,128,0,   0,128,  0,0,   0,128,128,0,
128,  0,  0,0, 128,  0,128,0, 128,128,  0,0, 192,192,192,0,
128,128,128,0,   0,  0,255,0,   0,255,  0,0,   0,255,255,0,
255,  0,  0,0, 255,  0,255,0, 255,255,  0,0, 255,255,255,0};
unsigned char monomap[]= {0,0,0,0, 255,255,255,0};

#pragma pack(1)
typedef struct tagRGBQUAD {
	byte	rgbBlue;
	byte	rgbGreen;
	byte	rgbRed;
	byte	rgbReserved;
} RGBQUAD;
typedef struct tagBITMAPFILEHEADER {
	WORD	bfType;
	DWORD	bfSize;
        WORD    bfReserved1;
        WORD    bfReserved2;
	DWORD	bfOffBits;
} BITMAPFILEHEADER;
typedef BITMAPFILEHEADER *PBITMAPFILEHEADER;

typedef struct tagBITMAPINFOHEADER{
  	DWORD	   biSize;
  	DWORD	   biWidth;
  	DWORD	   biHeight;
  	WORD	   biPlanes;
  	WORD	   biBitCount;

	DWORD	   biCompression;
	DWORD	   biSizeImage;
	DWORD	   biXPelsPerMeter;
	DWORD	   biYPelsPerMeter;
	DWORD	   biClrUsed;
	DWORD	   biClrImportant;
} BITMAPINFOHEADER;

BITMAPFILEHEADER hdr;
BITMAPINFOHEADER bi;	// Source image header
BITMAPINFOHEADER bo;	// Output image header

long wib;
long dwib;

RGBQUAD srccmap[256];
unsigned char cmap[256*4];	/* Bytes of blue, green, red, reserved */
unsigned char *scanline;
unsigned char used[256];

FILE *fpin;	/* Input file fp */

static int dib_wib(int bitcount, int wi){
	switch (bitcount){
	case 1: wi= (wi+31) >> 3; break;
	case 4:	wi= (wi+7)  >> 1; break;
	case 8:	wi=  wi+3; break;
	case 16:wi= (wi*2)+3; break;
        case 24:wi= (wi*3)+3; break;
	case 32:return wi*4;
	}
	return wi & ~3;
}
void err(char *msg){
	fprintf(stderr, "%s\n", msg);
	exit(5);
}
/* These are used for 16- and 32-bit DIBS */		
static DWORD rgbmask[3]= {0x7c00,0x03e0,0x001f};
static int rgbshft[3]= {7, 2, -3};
void gulp_header(char *dibname){
	int v, y;

	if((fpin= fopen (dibname, "rb"))==0){
		perror(dibname);
		exit(0);
	}
	fread(&hdr,1,sizeof(hdr),fpin);
	if(hdr.bfType!= 0x4d42)
		err("Image not a .bmp");

	fread(&bi, 1, sizeof(bi), fpin);
	wib= dib_wib((int)bi.biBitCount, (int)bi.biWidth);
	bi.biSizeImage= bi.biHeight * wib;
	if(bi.biBitCount<=8){
		v= (int)((hdr.bfOffBits - bi.biSize - sizeof(BITMAPFILEHEADER))/sizeof(RGBQUAD));
		bi.biClrUsed= v;
		if(v>0)
			fread(srccmap, sizeof(RGBQUAD), v, fpin);
	}else if((bi.biBitCount==16 || bi.biBitCount==32)
	&& bi.biCompression==BI_BITFIELDS){
		/* Read 3 masks */
		fread(rgbmask, sizeof(DWORD), 3, fpin);
		for(y= 0; y<3; y++){
			for(v= sizeof(DWORD); --v>=0 && 0==(rgbmask[y]&(1<<v)); v++);
			rgbshft[y]= v-8;
		}
	}
	v= (int)wib;
	if(v<(int)(3*bi.biWidth))
		v= (int)3*bi.biWidth;
	scanline= malloc(v);
	if(scanline==0)
		err("Not enough room for the scanline");
}
void readnextscanline(){
	int k= 0;
	int x;
	DWORD *dwp;
	WORD *wp;
	fread(scanline, 1, (int)wib, fpin);
	switch((int)bi.biBitCount){
	case 32:
		dwp= (DWORD *)scanline;
		for(x= 0; x<(int)bi.biWidth; x++){
			if(bi.biCompression==BI_BITFIELDS){
				scanline[k++]= (byte)((*dwp&rgbmask[2]) >> rgbshft[2]);
				scanline[k++]= (byte)((*dwp&rgbmask[1]) >> rgbshft[1]);
				scanline[k++]= (byte)((*dwp&rgbmask[0]) >> rgbshft[0]);
				++dwp;
			}else{
				scanline[k++]= scanline[4*x];
				scanline[k++]= scanline[4*x+1];
				scanline[k++]= scanline[4*x+2];
			}
		} break;
	case 24:
		break;
	case 16:
		k= (int)(3*bi.biWidth);
		wp= (WORD *)(&scanline[wib]);
		for(x= (int)bi.biWidth; --x>=0; ){
			--wp;
			scanline[--k]= (byte)((*wp&rgbmask[0]) >> rgbshft[0]);
			scanline[--k]= (byte)((*wp&rgbmask[1]) >> rgbshft[1]);
			scanline[--k]= (byte)((*wp&rgbmask[2]) >> rgbshft[2]);
		} break;
	case 8: case 4: case 1:
		k= (int)(3*bi.biWidth);
		for(x= (int)bi.biWidth; --x>=0; ){
			int pv;
			if(bi.biBitCount==8)
				pv= scanline[x];
			else if(bi.biBitCount==4)
				pv= (x&1? scanline[x>>1] : scanline[x>>1] >> 4) & 15;
			else
				pv= scanline[x>>3] & (1 << (x&7)) ? 255 : 0;
			scanline[--k]= srccmap[pv].rgbRed;
			scanline[--k]= srccmap[pv].rgbGreen;
			scanline[--k]= srccmap[pv].rgbBlue;
		} break;
	}
}
/* Output the image */
int spew_header(FILE *fp, unsigned char *cmap){
	int n= 0;

	dwib= dib_wib((int)bo.biBitCount, (int)bo.biWidth);
	if(bo.biBitCount==1)n= 2;
	if(bo.biBitCount==4)n= 16;
	if(bo.biBitCount==8)n= 256;
	bo.biSize= SIZE_BITMAPINFOHEADER;
	bo.biPlanes= 1;
	bo.biCompression= 0;
	bo.biSizeImage= bo.biHeight * dwib;
	bo.biClrUsed= n;
	bo.biClrImportant= n;
	bo.biXPelsPerMeter= 0;
	bo.biYPelsPerMeter= 0;

	hdr.bfType= 0x4d42;	/* BM */
	hdr.bfReserved1= 0;
	hdr.bfReserved2= 0;
	hdr.bfOffBits = SIZE_BITMAPFILEHEADER + bo.biSize + n*sizeof(RGBQUAD);
	hdr.bfSize= hdr.bfOffBits + bo.biSizeImage;

	fwrite(&hdr,1,sizeof(hdr),fp);
	fwrite(&bo, 1, sizeof(bo), fp);
	if(n>0)
		fwrite(cmap, sizeof(RGBQUAD), n, fp);
	return 1;
}
char *note[]= {
	"Usage: [Options] source.bmp dest.bmp",
	"Dither a 24-bit .bmp file to 8 bits (or 4 bits Windows colors)",
	"Options:",
	"-n### Dither to # red, # green, # blue values [default is -n666]",
	"-w    Dither to the 16 standard windows colors (4 bits/pixel output)",
	"-m    Dither to black and white (1 bit/pixel output)",
	"-g #  Set image gamma [default 1.85]",
	"-r    Send errors to the right [default alternates right/left]",
	"-i    Start with no noise (when background color is an output color)",
	0};

int tab7[1024], tab5[1024], tab3[1024], tab1[1024];
unsigned char neartable[16*16*16];
int *rtab, *gtab, *btab;
int gammatab[256];
unsigned char degamma[4*512];
int gammacmap[4*256];

/* The 'get the pvalue of a color near this color' routines: */

int pvnear(int r, int g, int b){	/* reds x greens x blues colormap */
	return rtab[degamma[512+r]] + gtab[degamma[512+g]] + btab[degamma[512+b]];
}
int pvtab_near(int r, int g, int b){	/* disordered colormap */
	r >>= 5; g >>= 5; b >>= 5;
	if((r|g|b)&-16){	// Color off cube, then force to a cube corner!
		r= r>7? 15:0;
		g= g>7? 15:0;
		b= b>7? 15:0;
	}
	return neartable[b+16*(g+16*r)];
}
int mono_near(int r, int g, int b){	/* one-bit output */
	return b + 2*(r+2*g) < threshhold? 0 : 1;
}

void neartable_make(int n){	/* Set up the disordered table from gammacmap */
	int r,g,b;
	int j= 0;
	for(r= 0; r<512; r+= 34)
	for(g= 0; g<512; g+= 34)
	for(b= 0; b<512; b+= 34){
		int pv= 0;
		int k;
		long dmin= 0x7fffffff;
		for(k= 0; k<n; k++){
			int k4= 4*k;
			long dr= r-gammacmap[k4+2];
			long dg= g-gammacmap[k4+1];
			long db= b-gammacmap[k4+0];
			long d= dr*dr+db*db+dg*dg;
			if(d<dmin){
				dmin= d;
				pv= k;
			}
		}
		neartable[j++]= pv;
	}
}
void createmap(unsigned char *map, int nr, int ng, int nb){
	int r,g,b,k;
	if(nr*ng*nb>256)
		return;

	for(k= 0; k<512; k++){
		tab7[512+k]= (7*k+8)/16;
		tab7[512-k]= -tab7[512+k];
		tab5[512+k]= (5*k+8)/16;
		tab5[512-k]= -tab5[512+k];
		tab3[512+k]= (3*k+8)/16;
		tab3[512-k]= -tab3[512+k];
		tab1[512+k]= (1*k+8)/16;
		tab1[512-k]= -tab1[512+k];
	}
	for(k= 0; k<256; k++)
		gammatab[k]= (int)(511.0*pow(k/255.0, gamma) + 0.5);
	for(k= 0; k<512; k++){
		degamma[k]= 0;
		degamma[k+512]= (int)(255.0*pow(k/511.0, 1/gamma) + 0.5);
		degamma[k+1024]= 255;
		degamma[k+1536]= 255;
	}
	if(outmode==OUT2){			// One-bit colormap
		threshhold= (int)(7*511*pow(0.5, 1/gamma) + 0.5);
		for(k= 0; k<sizeof(monomap); k++)
			map[k]= monomap[k];
		for(k= 0; k<4*2; k++)
			gammacmap[k]= map[k]? 511 : 0;
	}else if(outmode==OUT16){		// Windows colormap
		for(k= 0; k<sizeof(winmap); k++)
			map[k]= winmap[k];
		for(k= 0; k<4*16; k++)
			gammacmap[k]= gammatab[map[k]];
		neartable_make(16);
	}else{					// Ordered colormap
		for(k= r= 0; r<nr; r++){
			for(g= 0; g<ng; g++){
				for(b= 0; b<nb; b++){
					map[k++]= b*255/(nb-1);
					map[k++]= g*255/(ng-1);
					map[k++]= r*255/(nr-1);
					k++;
				}
			}
		}
		rtab= (int *)(&neartable[0]);
		gtab= rtab+256;
		btab= gtab+256;
		for(k= 0; k<256; k++){
			rtab[k]= ng * nb * ( k * nr/256 );
			gtab[k]= nb * ( k * ng/256 );
			btab[k]= k * nb/256;
		}
		for(k= 0; k<4*256; k++)
			gammacmap[k]= gammatab[map[k]];
	}
}
void fsdither(FILE *fo){
	int y, r, g, b, pv;
	int ru, gu, bu;
	unsigned char *outbuf= (char *)malloc((int)bo.biWidth);
	int x= 3*bi.biWidth+6;
	int *scanerr= (int *)malloc(sizeof(int)*x);	// Scanline of errors

	if(startnoisy){
		int mask;
		switch(bo.biBitCount){
		case 1:	mask= 511; break;
		case 4: mask= 255; break;
		case 8: mask= 128; break;
		}
		srand(12345);	// Always use the same noise
		while(--x>=0)scanerr[x]= 2*(rand()&mask)-mask;
	}else
		while(--x>=0)scanerr[x]= 0;	
	ru= 0; gu= 0; bu= 0;
	for(y= 0; y<(int)bi.biHeight; y++){
		unsigned char *pc;
		int *pe;
		readnextscanline();
		if(y&reversing){
			while(--x>=0){
				pc -= 3;
				pe -= 3;
				r= gammatab[pc[2]] + tab7[ru] + pe[2];
				g= gammatab[pc[1]] + tab7[gu] + pe[1];
				b= gammatab[pc[0]] + tab7[bu] + pe[0];
				pv= (*nearrtn)(r,g,b);
				outbuf[x]= pv;
				pv *= 4;
				r -= gammacmap[pv+2]-512; if(r>1023)r= 1023;
				g -= gammacmap[pv+1]-512; if(g>1023)g= 1023;
				b -= gammacmap[pv+0]-512; if(b>1023)b= 1023;
				pe[5] += tab3[r];
				pe[4] += tab3[g];
				pe[3] += tab3[b];
				pe[2]= tab5[r] + tab1[ru];
				pe[1]= tab5[g] + tab1[gu];
				pe[0]= tab5[b] + tab1[bu];
				ru= r; gu= g; bu= b;
			}
		}else{
			pc= scanline;
			pe= scanerr+3;
			for(x= 0; x<(int)bi.biWidth; x++){
				r= gammatab[pc[2]] + tab7[ru] + pe[2];
				g= gammatab[pc[1]] + tab7[gu] + pe[1];
				b= gammatab[pc[0]] + tab7[bu] + pe[0];
				pv= (*nearrtn)(r,g,b);
				outbuf[x]= pv;
				pv *= 4;
				r -= gammacmap[pv+2]-512; if(r>1023)r= 1023;
				g -= gammacmap[pv+1]-512; if(g>1023)g= 1023;
				b -= gammacmap[pv+0]-512; if(b>1023)b= 1023;
				pe[-1] += tab3[r];
				pe[-2] += tab3[g];
				pe[-3] += tab3[b];
				pe[2]= tab5[r] + tab1[ru];
				pe[1]= tab5[g] + tab1[gu];
				pe[0]= tab5[b] + tab1[bu];
				ru= r; gu= g; bu= b;
				pc += 3;
				pe += 3;
			}
		}
		if(bpp==1){
			b= 64;
			r= 0;
			if(outbuf[0])
				outbuf[0]= 128;
			for(x= 1; x<(int)bi.biWidth; x++){
				if(outbuf[x])
					outbuf[r] |= b;
				if(0==(b >>= 1)){
					b= 128;
					outbuf[++r]= 0;
				}
			}
		}else if(bpp==4){
			for(r= x= 0; x<(int)bi.biWidth; x+= 2)
				outbuf[r++]= (unsigned char)((outbuf[x]<<4) | outbuf[x+1]);
		}
		fwrite(outbuf, 1, (int)dwib, fo);
	}
}
void main(int ac, char **av){
	FILE *fp;
	char *cp;
	char *srcname= 0;
	char *dstname= 0;
	int nr= 6, ng= 6, nb= 6;

	bpp= 8;
	nearrtn= pvnear;

	while(cp= *++av)if(*cp++=='-')switch(*cp++){
	case 'i':	startnoisy= 0; break;
	case 'r':	reversing= 0; break;
	case 'w':	outmode= OUT16; nearrtn= pvtab_near; bpp= 4; break;
	case 'm':	outmode= OUT2;  nearrtn= mono_near; bpp= 1; break;
	case 'g':	gamma= atof(*++av); break;
	case 'n':	if(*cp==0)cp= *++av;
			nr= *cp++ - '0'; ng= *cp++ - '0'; nb= *cp++ - '0';
			break;
	default:
	erra:		for(av= note; *av; ++av)
				fprintf(stderr, "%s\n", *av);
			exit(1);
	}else if(srcname==0)
		srcname= *av;
	else if(dstname==0)
		dstname= *av;
	else
		goto erra;
	if(dstname==0)
		goto erra;

	/* Open the source bmp */
	gulp_header(srcname);

	/* Create the gif file */
	fp= fopen(dstname, "wb");
	if(fp==0){
		perror(dstname);
		exit(1);
	}
	bo= bi;	// The output image resembles the source image
	bo.biBitCount= bpp;	// With fewer bits per pixel
	createmap(cmap, nr, ng, nb);
	spew_header(fp, cmap);
	fsdither(fp);
}
/*(C)1996 Ephraim Cohen--
** permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notices appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
* 
*/