// converts images to PCM audio so you can compress them with mp3
// this works about as well as you'd think

#include <stdint.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <math.h>
#include <stdio.h>

// store x*x as a sample instead of x
// this is smaller and makes it sound bit more like "real" sound
// probably worse for lossless audio codecs though
#define SQR

static int16_t bswap16(int16_t a) {return (a >> 8) | ((a & 0xff) << 8);}
static int clip(int a, int min, int max) {if (a > max) a = max; if (a < min) a = min; return a;}

// Hilbert fractal curve
// http://forums.xkcd.com/viewtopic.php?f=3&t=866&p=19105#p19105
static int transform_table[4][4]={{0,1,2,3},{0,2,1,3},{3,2,1,0},{3,1,2,0}};
static int locations[4]={0,1,3,2};

void hilbert(int n,int *xp,int *yp)
{
	static int transforms[4]={1,0,0,3};
	int x=0,y=0;
	int trans=0;
	for(int i=30;i>=0;i-=2)
	{
		int m=(n>>i)&3;
		int bits=transform_table[trans][locations[m]];
		
		x=(x<<1)|((bits>>1)&1);
		y=(y<<1)|(bits&1);
		
		trans^=transforms[m];
	}
	*xp=x; *yp=y;
}

// unlike libjpeg yuv<>rgb, this is actually precise when inverted, since i don't suck at rounding
// hm, if it's going to end up in 16-bit anyway, maybe i shouldn't bother reducing to 8-bit yuv
static void yuv601torgb(unsigned char y, unsigned char cb, unsigned char cr, unsigned char *r, unsigned char *g, unsigned char *b)
{
    int _r,_g,_b;
    _r = abs(y * 9539 + cr * 13075 + 4096 - 1826224);
    _g = abs(y * 9539 + cb * -3209 + cr * -6660 + 4096 + 1110633);
    _b = abs(y * 9539 + cb * 16525 + 4096 - 2267824);
    *r = (_r<=2088960?_r:2088960)>>13;
    *g = (_g<=2088960?_g:2088960)>>13;
    *b = (_b<=2088960?_b:2088960)>>13;
}

static void rgbtoyuv601(unsigned char r, unsigned char g, unsigned char b, unsigned char *y, unsigned char *cb, unsigned char *cr)
{
    int _y,_cb,_cr;
    _y = abs(r * 2104 + g * 4130 + b * 802 + 4096 + 131072);
    _cb = abs(r * -1214 + g * -2384 + b * 3598 + 4096 + 1048576);
    _cr = abs(r * 3598 + g * -3013 + b * -585 + 4096 + 1048576);
    *y = clip((_y<=2088960?_y:2088960)>>13, 0, 235);
    *cb = clip((_cb<=2088960?_cb:2088960)>>13, 0, 240);
    *cr = clip((_cr<=2088960?_cr:2088960)>>13, 0, 240);
}

// YCbCr to stereo sound
// Y in center, "hue" (Cb - Cr) in stereo balance
// loses "saturation" (value of Cr)
static void forward_stereo(uint8_t y, uint8_t u, uint8_t v, float *l_, float *r_)
{
	float l, r;
	char sy, su, sv;
	
	sy = y - 128; // center of Y range is actually 125.5, maybe?
	su = u - 128;
	sv = v - 128;
	
	// we can safely scale by 128 and 32768 in this since yuv is always less than 255/256
	float fy = sy/128., fu = su/128., fv = sv/128.;
	
#ifdef SQR
	fy = copysignf(fy*fy, fy);
	fu = copysignf(fu*fu, fu);
	fv = copysignf(fv*fv, fv);
#endif SQR
	
	float hue = (fu-fv)/2.;
	
	l = fy + hue;
	r = fy - hue;
	
	*l_ = l;
	*r_ = r;
}

FILE *fout;
uint8_t *img;
int w, h;

// 16-bit samples are written native endian
static void encode_sample(uint8_t r, uint8_t g, uint8_t b)
{
	uint8_t y,u,v;
	float ls, rs;
	short pl, pr;
	
	rgbtoyuv601(r,g,b,&y,&u,&v);
	forward_stereo(y,u,v,&ls,&rs);
	
	pl = lrintf(ls * 32768.);
	pr = lrintf(rs * 32768.);
	
#ifdef __BIG_ENDIAN__
	pl = bswap16(pl);
	pr = bswap16(pr);
#endif
	
	fwrite(&pl, 2, 1, fout);
	fwrite(&pr, 2, 1, fout);
}

static void decode_sample(int16_t pl, int16_t pr, int x, int y)
{
	float ls = pl / 32768., rs = pr / 32768.;
	float hue = ls - rs;
	
	float fy = (ls+rs)/2., fu = hue, fv = -hue;
	
#ifdef SQR
	fy = copysignf(sqrtf(fabs(fy)), fy);
	fu = copysignf(sqrtf(fabs(fu)), fu);
	fv = copysignf(sqrtf(fabs(fv)), fv);
#endif
	
	int y_ = lrintf(fy * 128.) + 128, u = lrintf(fu * 128.) + 128, v = lrintf(fv * 128.) + 128;
	
	uint8_t *im = &img[(y*w + x)*3];
	
	yuv601torgb(clip(y_, 16, 235), clip(u, 16, 240), clip(v, 16, 240), &im[0], &im[1], &im[2]);
}

static size_t filesize(char *fn)
{
	struct stat st;
	stat(fn, &st);
	
	return st.st_size;
}

static int find_delay(int16_t *snd)
{
	//int d = 0;
	//while (abs(snd[d*2])<4096 && abs(snd[d*2+1])<4096) {d++;}
	//return d;
	if (!snd[0] && !snd[1]) return 1681;
	return 0;
}

int main(int argc, char *argv[])
{
	if (argc < 6) {
		printf("usage: c|d <input> <output> <width> <height>\n"
			   "image format: raw 8-bit RGB\n"
			   "sound format: 16-bit pcm at native endian, no wave header\n"
			   "tries to handle MP3 sample delay but not very well\n");
		return 0;
	}
	
	FILE *fin = fopen(argv[2], "rb");
	w = atoi(argv[4]); h = atoi(argv[5]);
	fout = fopen(argv[3], "wb");
	img = malloc(w*h*3);
	
	if (argv[1][0] == 'c') {
		fread(img, 1, w*h*3, fin);
		
		int i = 0, hi = 0;
		
		while (i < w*h) {
			int x, y;
			hilbert(hi++, &x, &y);
			
			if (x >= w || y >= h) continue;
			uint8_t *im = &img[(y * w + x)*3];
			
			encode_sample(im[0], im[1], im[2]);
			i++;
		}
	} else {
		size_t scount = filesize(argv[2]);
		int16_t *snd = malloc(scount); 
		fread(snd, 1, scount, fin);
		int delay = find_delay(snd);
		printf("estimated %d samples delay\n",delay);
		
		int i = 0, hi = 0;
		
		while (i < w*h) {
			int x, y;
			hilbert(hi++, &x, &y);
			
			if (x >= w || y >= h) continue;
			int16_t *s = &snd[(i+delay)*2];
			
			decode_sample(s[0], s[1], x, y);
			i++;
		}
		
		fwrite(img, 1, w*h*3, fout);
	}
	
	free(img);
	return 0;
}