/*
	STColor
	Color conversions and value->color maps
	
	Stilianos Louca
	September 2011
*/

#ifndef ST_COLOR_DEF
#define ST_COLOR_DEF

#include "STColor.h"
#include <iostream>

namespace ST{


#pragma mark Constants
#pragma mark -

	const double STCOLOR_RGB_to_LMS[3][3] = {	{17.8824,   43.5161,  4.11935},
												{3.45565,   27.1554,  3.86714},
												{0.0299566, 0.184309, 1.46709} };
	const double STCOLOR_LMS_to_RGB[3][3] = {	{0.0809444479,		-0.130504409,	0.116721066},
												{-0.0102485335,    	0.0540193266,  	-0.113614708},
												{-0.000365296938, 	-0.00412161469, 0.693511405} };
  	const double STCOLOR_RR_to_MOD[3][3]	=  {{0.0, 0.0, 0.0},
												{0.7, 1.0, 0.0},
												{0.7, 0.0, 1.0} };
  	const double STCOLOR_Deuteranope[3][3] = {	{1.0, 		0.0, 	0.0},
												{0.494207, 	0.0, 	1.24827},
												{0.0,      	0.0, 	1.0} };
   	const double STCOLOR_Protanope[3][3] = 	{	{0.0, 	2.02344, 	-2.52581},
												{0.0, 	1.0,      	0.0},
												{0.0, 	0.0,      	1.0} };
  	const double STCOLOR_Tritanope[3][3] = 	{	{1.0,   	0.0,    	0.0},
												{0.0,       1.0,      	0.0},
												{-0.395913, 0.801109, 	0.0} };
	const double STCOLOR_ERR_to_MOD[3][3] = {	{0.0, 0.0, 0.0},
												{0.7, 1.0, 0.0},
												{0.7, 0.0, 1.0} };
												
	const long N_STCOLOR_DISCRETE_COLOR_PALETTE = 34;
	char const * const STCOLOR_DISCRETE_COLOR_PALETTE[N_STCOLOR_DISCRETE_COLOR_PALETTE] = {"E30505", "3C41F0", "48BD5A", "ED9C02", "E81ECA", "2BD6CD", "9BB0BA", "35B842", "D69449", "e31a1c", "1f78b4", "0C064B","481800","0C064B","5A2D00","008B8B","550000","311EF1","D8D800","651265","0093F0D","4A2816","E00BE0","D76C00","8E008E","004141","616102","600592","9C1B1B","7E7E2D","1608A2","250239","002F03","370037"};
	
	const long N_STCOLOR_DISCRETE_COLOR_PALETTE_DARK = 29;
	char const * const STCOLOR_DISCRETE_COLOR_PALETTE_DARK[N_STCOLOR_DISCRETE_COLOR_PALETTE_DARK]= {"C70000", "292EC4", "308C3E", "916004", "AD1096", "1AA19A", "697F8A", "39573C", "855724", "e31a1c", "6a3d9a",  "1f78b4", "0C064B","481800","0C064B","5A2D00","550000","651265","0093F0D","4A2816","D76C00","004141","616102","600592","9C1B1B","1608A2","250239","002F03","370037"};

		
	
	
	
	
	
	
	

#pragma mark -
#pragma mark Color conversion
#pragma mark -


unsigned int RGB2ColorNumber(double r, double g, double b){
	return int(r * 255) * 65536 + int(g * 255) * 256 + int(b * 255) * 1;
}

void colorNumber2RGB(unsigned int colorNumber, double &r, double &g, double &b){
	colorNumber = max(0u, min((unsigned int)(256 * 256 * 256- 1), colorNumber));
	r	= floor(colorNumber / 65536)/255.0;
	g	= floor((colorNumber % 65536) / 256) / 255.0;
	b	= (colorNumber % 256) / 255.0;
}


void RGB2HSV(double r, double g, double b, double &h, double &s, double &v){
	double M = max(r,max(g,b)), m = min(r, min(g,b));
	double C = M - m;
	double EPS = 0.0001;
	
	//hue
	if(C <= EPS){
		h = 0.0;
	}else if((M-r) <= EPS){
		h = ((g-b)/C + 0.0)/6.0;
	}else if((M-g) <= EPS){
		h = ((b-r)/C + 2.0)/6.0;
	}else if((M-b) <= EPS){
		h = ((r-g)/C + 4.0)/6.0;
	}
	h = fmod(1.0 + h, 1.0);
	
	//saturation
	if(M <= EPS){
		s = 0.0;
	}else{
		s = C/M;
	}
	
	//value
	v = M;
}



//implementation taken from: http://de.wikipedia.org/wiki/HSV-Farbraum, 04.09.2011
void HSV2RGB(double h, double s, double v, double &r, double &g, double &b){
	//h is modulo 1, so project to interval [0,1]
	if(h >= 0.0){
		h = h - floor(h);
	}else{
		h = h + ceil(abs(h));
	}
	unsigned int hi = (unsigned int) floor(h * 6.0);
	double 	f = (h * 6.0 - hi);
	double 	p = v * (1-s), 
			q = v * (1-s*f), 
			t = v * (1-s*(1-f));
	switch(hi){
	case 0: r=v; g=t; b=p; break;
	case 1: r=q; g=v; b=p; break;
	case 2: r=p; g=v; b=t; break;
	case 3: r=p; g=q; b=v; break;
	case 4: r=t; g=p; b=v; break;
	case 5: r=v; g=p; b=q; break;
	case 6: r=v; g=t; b=p; break;
	}
}



inline void STCOLOR_left_matmul3(double vec[3], const double mat[3][3], double scratch[3]){
	int i;
	for(i=0; i<3; ++i){
		scratch[i] = vec[0]*mat[0][i] + vec[1]*mat[1][i] + vec[2]*mat[2][i];
	}
	for(i=0; i<3; ++i){
		vec[i] = scratch[i];
	}
}

inline void STCOLOR_right_matmul3(const double mat[3][3], double vec[3], double scratch[3]){
	int i;
	for(i=0; i<3; ++i){
		scratch[i] = vec[0]*mat[i][0] + vec[1]*mat[i][1] + vec[2]*mat[i][2];
	}
	for(i=0; i<3; ++i){
		vec[i] = scratch[i];
	}
}



#pragma mark -
#pragma mark Complex numbers to colors
#pragma mark -




void complex2HSV(double realPart, double imaginaryPart, double &h, double &s, double &v){
	complex2HSV(complex<double>(realPart,imaginaryPart), 0.0, 1.0, h, s, v);
}


void complex2HSV(double realPart, double imaginaryPart, double zeroRadius, double inftyRadius, double &h, double &s, double &v){
	/*
	h = atan2(imaginaryPart, realPart) / (2.0 * PI);
	h = (h < 0.0 ? h + 1.0 : h);
	
	double norm = sqrt(pow(realPart, 2.0) + pow(imaginaryPart, 2.0));
	double a,b;
	a = inftyRadius; b = 0.05 / (EPSILON + inftyRadius);
	s = 0.5 - 0.5 * (b*norm - b*a*a/norm)/(1.0 + abs(b*norm - b*a*a/norm));
	
	a = zeroRadius; b = 0.33 / (EPSILON + zeroRadius);
	v = 0.5 + 0.5 * (b*norm - b*a*a/norm)/(1.0 + abs(b*norm - b*a*a/norm));
	*/
	complex2HSV(complex<double>(realPart, imaginaryPart), zeroRadius, inftyRadius, h, s, v);
}


void complex2HSV(complex<double> z, double zeroRadius, double inftyRadius, double &h, double &s, double &v){
	complexPolar2HSV(arg(z), abs(z), zeroRadius, inftyRadius, h, s, v);
}


void complexPolar2HSV(double argument, double norm, double zeroRadius, double inftyRadius, double &h, double &s, double &v){
	h = argument / (2.0 * PI);
	if(h >= 0.0){
		h = h - floor(h);
	}else{
		h = h + ceil(abs(h));
	}
	
	double a,b;
	a = inftyRadius; b = 0.05 / (EPSILON + inftyRadius);
	s = 0.5 - 0.5 * (b*norm - b*a*a/norm)/(1.0 + abs(b*norm - b*a*a/norm));
	
	a = zeroRadius; b = 0.33 / (EPSILON + zeroRadius);
	v = 0.5 + 0.5 * (b*norm - b*a*a/norm)/(1.0 + abs(b*norm - b*a*a/norm));
}
	


void complex2RGB(double realPart, double imaginaryPart, double &r, double &g, double &b){
	complex2HSV(realPart, imaginaryPart, r, g, b);
	HSV2RGB(r, g, b, r, g, b);
}


void complex2RGB(double realPart, double imaginaryPart, double zeroRadius, double inftyRadius, double &r, double &g, double &b){
	complex2HSV(complex<double>(realPart,imaginaryPart), zeroRadius, inftyRadius, r, g, b);
	HSV2RGB(r, g, b, r, g, b);
}

void complex2RGB(complex<double> z, double zeroRadius, double inftyRadius, double &r, double &g, double &b){
	complexPolar2HSV(arg(z), abs(z), zeroRadius, inftyRadius, r, g, b);
	HSV2RGB(r, g, b, r, g, b);
}

void complexPolar2RGB(double argument, double norm, double zeroRadius, double inftyRadius, double &r, double &g, double &b){
	complexPolar2HSV(argument, norm, zeroRadius, inftyRadius, r, g, b);
	HSV2RGB(r, g, b, r, g, b);
}




#pragma mark Vectors to colors
#pragma mark -

Color vector2RGBColor(double vx, double vy, double vz, double norm, VectorColoring coloring){
	double r,g,b;
	switch(coloring){
	case VectorColoringNone:
		r = g = b = 0;
		break;
	case VectorColoringNormGrey:
		r = g = b = 0.8 * norm;
		break;
	case VectorColoringXGrey:
		r = g = b = 0.8 * vx;
		break;
	case VectorColoringYGrey:
		r = g = b = 0.8 * vy;
		break;
	case VectorColoringZGrey:
		r = g = b = 0.8 * vz;
		break;
	case VectorColoringNormHeat:
		value2RGBViaPalette(0.8 * norm, ColorPaletteHeat, r, g, b);
		break;
	case VectorColoringXHeat:
		value2RGBViaPalette(0.8 * vx, ColorPaletteHeat, r, g, b);
		break;
	case VectorColoringYHeat:
		value2RGBViaPalette(0.8 * vy, ColorPaletteHeat, r, g, b);
		break;
	case VectorColoringZHeat:
		value2RGBViaPalette(0.8 * vz, ColorPaletteHeat, r, g, b);
		break;
	case VectorColoringNormOcean:
		value2RGBViaPalette(0.8 * norm, ColorPaletteOcean, r, g, b);
		break;
	case VectorColoringXOcean:
		value2RGBViaPalette(0.8 * vx, ColorPaletteOcean, r, g, b);
		break;
	case VectorColoringYOcean:
		value2RGBViaPalette(0.8 * vy, ColorPaletteOcean, r, g, b);
		break;
	case VectorColoringZOcean:
		value2RGBViaPalette(0.8 * vz, ColorPaletteOcean, r, g, b);
		break;
	case VectorColoringXYRedBlue:
		r	= vx;
		g 	= 0.0;
		b	= vy;
		break;
	case VectorColoringXZRedBlue:
		r 	= vx; 
		g 	= 0.0;
		b	= vz;
		break;
	case VectorColoringYZRedBlue:
		r	= vy;
		g	= 0.0;
		b	= vz;
		break;
	case VectorColoringXYHeat:
		r	= min(1.0, vx + 2 * vy);
		g	= vy;
		b	= 0.0;
		break;
	case VectorColoringXZHeat:
		r	= min(1.0, vx + 2 * vz);
		g	= vz;
		b	= 0.0;
		break;
	case VectorColoringYZHeat:
		r	= min(1.0, vy + 2 * vz);
		g	= vz;
		b	= 0.0;
		break;
	case VectorColoringXYZRGB:
		r	= 0.8 * vx;
		g	= 0.8 * vy;
		b	= 0.8 * vz;
		break;
	default:
		r = g = b = 0;
	}
	return Color(r,g,b);
}



void generateDiscreteColors(long count, 
							ColorPalette palette, 
							bool random, 
							bool preferDark, 
							vector<Color> &colors){
	colors.clear(); colors.reserve(count);
	if(count==0) return;
	long i,j;
	double minV = 0, maxV = 1, r,g,b, dV;
	if(preferDark){
		int bd = paletteBrighteningDirection(palette);
		if(bd>0){
			maxV = 0.5;
		}else if(bd<0){
			minV = 0.5;
		}
	}
	dV = (maxV - minV)/count;
	
	if(random){
		vector<double> values(count); 
		for(i=0; i<count; ++i){ values[i] = i; }
		std::random_shuffle(values.begin(), values.end());
		for(i=0; i<count; ++i){ 
			value2RGBViaPalette(minV + dV*values[i], palette, r, g, b);
			colors.push_back(Color(r,g,b));
		}
	}else{
		for(i=0; i<count; ++i){ 
			value2RGBViaPalette(minV + dV*i, palette, r, g, b);
			colors.push_back(Color(r,g,b));
		}
	}
}


void generateDiscreteColors(long count,
							bool preferDark,
							vector<Color> &colors,
							ColorBlindness correctForColorBlindness){
	if(count<=(preferDark ? N_STCOLOR_DISCRETE_COLOR_PALETTE_DARK : N_STCOLOR_DISCRETE_COLOR_PALETTE)){
		colors.clear(); colors.reserve(count);
		for(long n=0; n<count; ++n){
			colors.push_back(Color(preferDark ? STCOLOR_DISCRETE_COLOR_PALETTE_DARK[n] : STCOLOR_DISCRETE_COLOR_PALETTE[n]));
		}
	}else{
		generateDiscreteColors(count, ColorPaletteRainbow, true, preferDark, colors);
	}	
	
	if(correctForColorBlindness!=ColorBlindnessNone){
		//correct colors for color blindness
		for(long n=0; n<count; ++n){
			colors[n].daltonize(correctForColorBlindness);
		}
	}
}





#pragma mark -
#pragma mark Color palettes
#pragma mark -



void colorPalette2PaletteCodes(ColorPalette palette, int &redCode, int &greenCode, int &blueCode){
	switch(palette){
	case ColorPaletteGrey:
		redCode = 		greenCode = 		blueCode = 3;
		break;
	case ColorPaletteIGrey:
		redCode = 		greenCode = 		blueCode = -3;
		break;
	case ColorPaletteHeat:
		redCode = 21; 	greenCode = 22; 	blueCode = 23;
		break;
	case ColorPaletteIHeat:
		redCode = -21; 	greenCode = -22; 	blueCode = -23;
		break;
	case ColorPaletteOcean:
		redCode = 23; 	greenCode = 28; 	blueCode = 3;
		break;
	case ColorPaletteIOcean:
		redCode = -23; 	greenCode = -28; 	blueCode = -3;
		break;
	case ColorPaletteRainbow:
		redCode = 33; 	greenCode = 13; 	blueCode = 10;
		break;
	case ColorPaletteIRainbow:
		redCode = -33; 	greenCode = -13; 	blueCode = -10;
		break;
	case ColorPaletteColorPrintableOnGrey:
		redCode = 30; 	greenCode = 31; 	blueCode = 32;
		break;
	case ColorPaletteIColorPrintableOnGrey:
		redCode = -30; 	greenCode = -31; 	blueCode = -32;
		break;
	case ColorPaletteParrot:
		redCode = 7; 	greenCode = 5; 		blueCode = 15;
		break;
	case ColorPaletteIParrot:
		redCode = -7; 	greenCode = -5; 	blueCode = -15;
		break;
	case ColorPalettePassionFruit:
		redCode = 3; 	greenCode = 11; 	blueCode = 6;
		break;
	case ColorPaletteIPassionFruit:
		redCode = -3; 	greenCode = -11; 	blueCode = -6;
		break;
	default:
		redCode = greenCode = blueCode = 3;
	}
}



	  
double paletteMapping(int mappingCode, double x){
	if(mappingCode < 0){
		return 1.0 - paletteMapping(- mappingCode, x);
	}
	
	switch(mappingCode){
	case 0: return 0;
	case 1: return 0.5;
	case 2: return 1;
	case 3: return x;
	case 4: return pow(x,2);
	case 5: return pow(x,3);
	case 6: return pow(x,4);
	case 7: return sqrt(x);
	case 8: return pow(x,0.25);
	case 9: return sin(PI * x / 2.0);
	case 10: return cos(PI * x / 2.0);
	case 11: return abs(x - 0.5);
	case 12: return pow(2 * x - 1.0, 2);
	case 13: return sin(PI * x);
	case 14: return abs(cos(PI * x));
	case 15: return max(0.0, sin(2 * PI * x));
	case 16: return max(0.0, cos(2 * PI * x));
	case 17: return abs(sin(2 * PI * x));
	case 18: return abs(cos(2 * PI * x));
	case 19: return abs(sin(4 * PI * x));
	case 20: return abs(cos(4 * PI * x));
	case 21: return min(1.0, 3 * x);
	case 22: return min(1.0, max(0.0, 3 * x - 1.0));
	case 23: return min(1.0, max(0.0, 3 * x - 2.0));
	case 24: return min(1.0, abs(3 * x - 1.0));
	case 25: return min(1.0, abs(3 * x - 2.0));
	case 26: return max(0.0, (3 * x - 1.0)/2.0);
	case 27: return max(0.0, (3 * x - 2.0)/2.0);
	case 28: return abs((3 * x - 1.0)/2.0);
	case 29: return abs((3 * x - 2.0)/2.0);
	case 30: return min(1.0, max(0.0, x/0.32 - 0.78125));
	case 31: return min(1.0, max(0.0, 2 * x - 0.84));
	case 32: return (x <= 0.25 ? 4 * x : (x <= 0.42 ? 1 : (x <= 0.92 ? -2 * x + 1.84 : (x/0.08 - 11.5))));
	case 33: return min(1.0, abs(2 * x - 0.5));
	case 34: return min(1.0, 2 * x);
	case 35: return min(1.0, max(0.0, 2 * x - 0.5));
	case 36: return max(0.0, 2 * x - 1.0);
	default: return 0;
	}
}



void value2RGBViaPalette(double value, ColorPalette palette, double &r, double &g, double &b){
	int redCode,greenCode,blueCode;
	colorPalette2PaletteCodes(palette, redCode, greenCode, blueCode);
	value 	= max(0.0, min(1.0, value));
	r		= paletteMapping(redCode, value);
	g		= paletteMapping(greenCode, value);
	b		= paletteMapping(blueCode, value);
}


unsigned int value2ColorNumberViaPalette(double value, ColorPalette palette){
	double r,g,b;
	value2RGBViaPalette(value, palette, r, g, b);
	return RGB2ColorNumber(r, g, b);
}


int	paletteBrighteningDirection(ColorPalette palette){
	switch(palette){
	case ColorPaletteGrey:
	case ColorPaletteHeat:
	case ColorPaletteOcean:
		return 1;
	case ColorPaletteIGrey:
	case ColorPaletteIHeat:
	case ColorPaletteIOcean:
		return -1;
	case ColorPaletteRainbow:
	case ColorPaletteColorPrintableOnGrey:
	case ColorPaletteParrot:
	case ColorPalettePassionFruit:
		return 0;
	default:
		return 1;
	}
}



#pragma mark Color class
#pragma mark -

Color::Color(double r, double g, double b){
	red = max(0.0, min(1.0, r));
	green = max(0.0, min(1.0, g));
	blue = max(0.0, min(1.0, b));
}



Color::Color(const string &HexRGB){
	if(!setHex(HexRGB)) red=green=blue=0;
}



bool Color::setRGB(double r, double g, double b){
	if((r<0) || (r>1.0 + EPSILON) || (g<0) || (g>1.0+EPSILON) || (b<0) || (b>1.0+EPSILON)) return false;
	red = r;
	green = g;
	blue = b;
	return true;
}

bool Color::setHSV(double H, double S, double V){
	if((H<0) || (H>1.0 + EPSILON) || (S<0) || (S>1.0+EPSILON) || (V<0) || (V>1.0+EPSILON)) return false;
	HSV2RGB(H, S, V, red, green, blue);
	return true;
}

bool Color::setHex(const string &hex){
	if(hex.size() != 6) return false;
	for(unsigned int d = 0; d<hex.size(); ++d){
		if(!isHexChar(hex[d])) return false;
	}
	red = (hexchar2int(hex[0])*16.0 + hexchar2int(hex[1]))/255.0;
	green = (hexchar2int(hex[2])*16.0 + hexchar2int(hex[3]))/255.0;
	blue = (hexchar2int(hex[4])*16.0 + hexchar2int(hex[5]))/255.0;
	return true;
}


void Color::getRGB(double &r, double &g, double &b) const{
	r = red;
	g = green;
	b = blue;
}


void Color::getHSV(double &H, double &S, double &V) const{
	RGB2HSV(red, green, blue, H, S, V);
}


string 	Color::decimal2hex(unsigned int d) const{
	unsigned int d1, d2;
	d2 = d%16;
	d1 = (d-d2)/16;
	char hexMap[] = {'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'};
	string s = "00";
	s[0] = hexMap[d1];
	s[1] = hexMap[d2];
	return s;
}


int Color::hexchar2int(char hex) const{
	if((hex >= 48) && (hex <= 57)){
		return hex-48;
	}
	if((hex >= 65) && (hex <= 70)){
		return 10 + hex - 65;
	}
	if((hex >= 97) && (hex <= 102)){
		return 10 + hex - 97;
	}
	return -1;
}

bool	Color::isHexChar(char c){
	return (((c>=48) && (c<=57)) || ((c>=65) && (c<=70)) || ((c>=97) && (c<=102)));
}


string Color::getHex() const{
	return decimal2hex(red*255.0) + decimal2hex(green*255.0) + decimal2hex(blue*255.0);
}


bool Color::isWhite() const{
	return ((red > 1.0 - EPSILON) && (green > 1.0 - EPSILON) && (blue > 1.0 - EPSILON));
}


bool Color::isBlack() const{
	return ((red < EPSILON) && (green < EPSILON) && (blue < EPSILON));
}

Color Color::white(){
	return Color(1,1,1);
}

Color Color::black(){
	return Color(0,0,0);
}


void Color::saturate(double factor){
	double H,S,V;
	getHSV(H,S,V);
	setHSV(H,max(0.0,min(1.0,S*factor)),V);
}

		
void Color::turnGray(){
	// use 'luminosity' formula
	red = green = blue = 0.21*red + 0.72*green + 0.07*blue;
}


void Color::daltonize(ColorBlindness CB){
	if(CB == ColorBlindnessNone) return; //nothing to be done
	double z[3], scratch[3], RGB[3];
	int i;
	getRGB(RGB[0], RGB[1], RGB[2]);
	for(i=0; i<3; ++i){ z[i] = RGB[i]; }
	
	// transform into LMS colorspace
	STCOLOR_right_matmul3(STCOLOR_RGB_to_LMS, z, scratch);
	
	// transform into a color blind's perspective
	switch(CB){
	case ColorBlindnessDeuteranope:
		STCOLOR_right_matmul3(STCOLOR_Deuteranope, z, scratch);
		break;
	case ColorBlindnessProtanope:
		STCOLOR_right_matmul3(STCOLOR_Protanope, z, scratch);
		break;
	case ColorBlindnessTritanope:
		STCOLOR_right_matmul3(STCOLOR_Tritanope, z, scratch);
		break;
	}
	
	// transform back into RGB colorspace
	STCOLOR_right_matmul3(STCOLOR_LMS_to_RGB, z, scratch); 
	
	// calculate deviance from normal
	for(i=0; i<3; ++i){ z[i] = RGB[i] - z[i]; }		
	STCOLOR_right_matmul3(STCOLOR_ERR_to_MOD, z, scratch); 
	
	// apply compensation
	for(i=0; i<3; ++i){ z[i] = max(0.0,min(1.0,z[i] + RGB[i])); }
	
	// clamp
	for(i=0; i<3; ++i){ z[i] = min(1.0, max(0.0, z[i])); }

	setRGB(z[0], z[1], z[2]);
}


void Color::printRGB(ostream &s, bool bracket) const{
	if(bracket){
		s << "(" << red << ", " << green << ", " << blue << ")";
	}else{
		s << red << " " << green << " " << blue;
	}
}


Color Color::getNicePlotColor(unsigned long number){
	return Color(STCOLOR_DISCRETE_COLOR_PALETTE_DARK[number % N_STCOLOR_DISCRETE_COLOR_PALETTE_DARK]);
}




}//end of namespace ST

#endif