/*
 *  This file is part of RawTherapee.
 *
 *  Copyright (c) 2004-2010 Gabor Horvath <hgabor@rawtherapee.com>
 *
 *  RawTherapee is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 * 
 *  RawTherapee is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with RawTherapee.  If not, see <http://www.gnu.org/licenses/>.
 */
#ifndef _COLORTEMP_
#define _COLORTEMP_

#include <gtkmm.h>
#include <cmath>
#include "colortemp.h"
#include "imagefloat.h"
#include "labimage.h"
#define MAXR(a,b) ((a) > (b) ? (a) : (b))
#define pow_F(a,b) (exp(b*log(a)))


namespace rtengine {

#define MINTEMP 2000
#define MAXTEMP 25000
#define MINGREEN 0.02
#define MAXGREEN 5.0
#define INITIALBLACKBODY 4000

class ColorTemp {

	private:
		double temp;
		double green;
		Glib::ustring method;

		static void clip (double &temp, double &green);

	public:

		ColorTemp () : temp(-1), green(-1), method("Custom") {}
		ColorTemp (double t, double g, Glib::ustring m);
		ColorTemp (double mulr, double mulg, double mulb);

		inline double getTemp ()    { return temp;  }
		inline double getGreen ()   { return green; }

		void   getMultipliers (double &mulr, double &mulg, double &mulb) { temp2mul (temp, green, mulr, mulg, mulb); }

		void mul2temp (double rmul, double gmul, double bmul, double& temp, double& green);
		void temp2mul (double temp, double green, double& rmul, double& gmul, double& bmul);
		static void temp2mulxyz (double tem, double gree, Glib::ustring method, double &Xxyz, double &Zxyz);

		int XYZtoCorColorTemp(double x0,double y0 ,double z0, double &temp);
		static void cieCAT02(double Xw, double Yw, double Zw,double &CAM02BB00,double &CAM02BB01,double &CAM02BB02, double &CAM02BB10,double &CAM02BB11,double &CAM02BB12,double &CAM02BB20,double &CAM02BB21,double &CAM02BB22, double adap );
		//static	void CAT02 (Imagefloat* baseImg, const ProcParams* params);
		//static void ciecam_02 (LabImage* lab, const ProcParams* params);

		static double d_factor( double f, double la ) {
			return f * (1.0 - ((1.0 / 3.6) * exp((-la - 42.0) / 92.0)));
		}
		static float d_factorfloat( float f, float la ) {
			return f * (1.0f - ((1.0f / 3.6f) * exp((-la - 42.0f) / 92.0f)));
		}

		static double calculate_fl_from_la_ciecam02( double la ) {
			double la5 = la * 5.0;
			double k = 1.0 / (la5 + 1.0);

			/* Calculate k^4. */
			k = k * k;
			k = k * k;

			return (0.2 * k * la5) + (0.1 * (1.0 - k) * (1.0 - k) * pow(la5, 1.0 / 3.0));
		}
		static float calculate_fl_from_la_ciecam02float( float la ) {
			float la5 = la * 5.0f;
			float k = 1.0f / (la5 + 1.0f);

			/* Calculate k^4. */
			k = k * k;
			k = k * k;

			return (0.2f * k * la5) + (0.1f * (1.0f - k) * (1.0f - k) * pow_F(la5, 1.0f / 3.0f));
		}

		static double achromatic_response_to_white( double x, double y, double z, double d, double fl, double nbb, int gamu ) {
			double r, g, b;
			double rc, gc, bc;
			double rp, gp, bp;
			double rpa, gpa, bpa;
			gamu=1;
			xyz_to_cat02( r, g, b, x, y, z, gamu );

			rc = r * (((y * d) / r) + (1.0 - d));
			gc = g * (((y * d) / g) + (1.0 - d));
			bc = b * (((y * d) / b) + (1.0 - d));

			cat02_to_hpe( rp, gp, bp, rc, gc, bc, gamu );
			if(gamu==1){//gamut correction M.H.Brill S.Susstrunk
				rp=MAXR(rp,0.0);
				gp=MAXR(gp,0.0);
				bp=MAXR(bp,0.0);
			}

			rpa = nonlinear_adaptation( rp, fl );
			gpa = nonlinear_adaptation( gp, fl );
			bpa = nonlinear_adaptation( bp, fl );

			return ((2.0 * rpa) + gpa + ((1.0 / 20.0) * bpa) - 0.305) * nbb;
		}
		
		static float achromatic_response_to_whitefloat( float x, float y, float z, float d, float fl, float nbb, int gamu ) {
			float r, g, b;
			float rc, gc, bc;
			float rp, gp, bp;
			float rpa, gpa, bpa;
			gamu=1;
			xyz_to_cat02float( r, g, b, x, y, z, gamu );

			rc = r * (((y * d) / r) + (1.0f - d));
			gc = g * (((y * d) / g) + (1.0f - d));
			bc = b * (((y * d) / b) + (1.0f - d));

			cat02_to_hpefloat( rp, gp, bp, rc, gc, bc, gamu );
			if(gamu==1){//gamut correction M.H.Brill S.Susstrunk
				rp=MAXR(rp,0.0f);
				gp=MAXR(gp,0.0f);
				bp=MAXR(bp,0.0f);
			}

			rpa = nonlinear_adaptationfloat( rp, fl );
			gpa = nonlinear_adaptationfloat( gp, fl );
			bpa = nonlinear_adaptationfloat( bp, fl );

			return ((2.0f * rpa) + gpa + ((1.0f / 20.0f) * bpa) - 0.305f) * nbb;
		}

		static void xyz_to_cat02 ( double &r,  double &g,  double &b,  double x, double y, double z, int gamu );
		static void cat02_to_hpe ( double &rh, double &gh, double &bh, double r, double g, double b, int gamu );
		static void cat02_to_xyz ( double &x,  double &y,  double &z,  double r, double g, double b, int gamu );
		static void hpe_to_xyz   ( double &x,  double &y,  double &z,  double r, double g, double b );

		static void xyz_to_cat02float ( float &r,  float &g,  float &b,  float x, float y, float z, int gamu );
		static void cat02_to_hpefloat ( float &rh, float &gh, float &bh, float r, float g, float b, int gamu );
		static void cat02_to_xyzfloat ( float &x,  float &y,  float &z,  float r, float g, float b, int gamu );
		static void hpe_to_xyzfloat   ( float &x,  float &y,  float &z,  float r, float g, float b );

		static void Aab_to_rgb( double &r, double &g, double &b, double A, double aa, double bb, double nbb );
		static void Aab_to_rgbfloat( float &r, float &g, float &b, float A, float aa, float bb, float nbb );
		static void calculate_ab( double &aa, double &bb, double h, double e, double t, double nbb, double a );
		static void calculate_abfloat( float &aa, float &bb, float h, float e, float t, float nbb, float a );

		
		static double nonlinear_adaptation( double c, double fl ) {
		double p;
			if(c<0.0){ p = pow( (-1.0*fl * c) / 100.0, 0.42 );return ((-1.0*400.0 * p) / (27.13 + p)) + 0.1;}
			else {p = pow( (fl * c) / 100.0, 0.42 ); return ((400.0 * p) / (27.13 + p)) + 0.1;}
		}
		static float nonlinear_adaptationfloat( float c, float fl ) {
		float p;
			if(c<0.0f){ p = pow_F( (-1.0f*fl * c) / 100.0f, 0.42f );return ((-1.0f*400.0f * p) / (27.13f + p)) + 0.1f;}
			else {p = pow_F( (fl * c) / 100.0f, 0.42f ); return ((400.0f * p) / (27.13f + p)) + 0.1f;}
		}

		static double inverse_nonlinear_adaptation( double c, double fl ) {
		    int c1;
		    if(c-0.1 < 0.0) c1=-1; else c1=1;
			return c1*(100.0 / fl) * pow( (27.13 * fabs( c - 0.1 )) / (400.0 - fabs( c - 0.1 )), 1.0 / 0.42 );
		}
		static float inverse_nonlinear_adaptationfloat( float c, float fl ) {
		    int c1;
		    if(c-0.1f < 0.0f) c1=-1; else c1=1;
			return c1*(100.0f / fl) * pow_F( (27.13f * fabs( c - 0.1f )) / (400.0f - fabs( c - 0.1f )), 1.0f / 0.42f );
		}
		
		
		static void curvecolor(double satind, double satval, double &sres, double parsat); 				
		static void curvecolorfloat(float satind, float satval, float &sres, float parsat); 				
		static void curveJ (double br, double contr, int db, LUTf & outCurve , LUTu & histogram ) ;
		static void curveJfloat (float br, float contr, int db, LUTf & outCurve , LUTu & histogram ) ;

		bool operator== (const ColorTemp& other) { return fabs(temp-other.temp)<1e-10 && fabs(green-other.green)<1e-10; }
		bool operator!= (const ColorTemp& other) { return !(*this==other); }

		static double blackbody_spect (double wavelength, double m1, double m2, double temp);
		static double daylight_spect  (double wavelength, double m1, double m2, double temp);
		static const double Cloudy6200_spect[97];
		static const double Daylight5300_spect[97];
		static const double Shade7600_spect[97];
		static const double A2856_spect[97];
		static const double FluoF1_spect[97];
		static const double FluoF2_spect[97];
		static const double FluoF3_spect[97];
		static const double FluoF4_spect[97];
		static const double FluoF5_spect[97];
		static const double FluoF6_spect[97];
		static const double FluoF7_spect[97];
		static const double FluoF8_spect[97];
		static const double FluoF9_spect[97];
		static const double FluoF10_spect[97];
		static const double FluoF11_spect[97];
		static const double FluoF12_spect[97];
		static const double HMI_spect[97];
		static const double GTI_spect[97];
		static const double JudgeIII_spect[97];
		static const double Solux3500_spect[97];
		static const double Solux4100_spect[97];
		static const double Solux4700_spect[97];
		static const double NG_Solux4700_spect[97];
		static const double NG_LEDLSI2040_spect[97];
		static const double NG_CRSSP12WWMR16_spect[97];
		static const double Flash5500_spect[97];
		static const double Flash6000_spect[97];
		static const double Flash6500_spect[97];

		//spectral data 8 color Colorchecker24
		static double get_spectral_color (double wavelength, const double* array) {
			int wlm = (int) ((wavelength -350.)/5.);
			return (array[wlm]);
		}

		static const double ColorchechredC3_spect[97];
		static const double ColorchechOraA2_spect[97];
		static const double ColorchechYelD3_spect[97];
		static const double ColorchechGreE2_spect[97];
		static const double ColorchechGreB3_spect[97];
		static const double ColorchechCyaF3_spect[97];
		static const double ColorchechPurD2_spect[97];
		static const double ColorchechMagE3_spect[97];
		static const double ColorchechSkiA138_13_14_spect[97];
		static const double ColorchechGraC4_67_spect[97];
		static const double ColorchechSkiB166_18_18_spect[97];
		static const double ColorchechBluC150_m5_m22_spect[97];
		static const double ColorchechDCBluN881_m7_m14_spect[97];//ColorChecker DC  N8
		static const double ColorchechSGSkiF763_14_26_spect[97];//ColorChecker SG  F7
		static const double ColorchechSGSkiK285_11_17_spect[97];//ColorChecker SG  K2
		static const double ColorchechWhiA496_spect[97];
		static const double ColorchechGreD1_spect[97];
		static const double ColorchechSGBlaN3_6_spect[97];//ColorChecker SG  N3
		static const double JDC468_GraK14_44_spect[97];//468  K14
		static const double JDC468_BluH10_spect[97];//468  H10
		static const double ColabSkin35_15_17_spect[97];//Skin L 35
		static const double ColabSkin57_22_18_spect[97];//Skin L 57
		static const double ColabSkin40_17_17_spect[97];//Skin L 40
		static const double ColabSkin91_4_14_spect[97];//Skin L 91
		static const double ColabSkin87_8_8_spect[97];//Skin L 87
		static const double ColabSkin89_8_21_spect[97];//Skin L 89
		static const double ColabSkin75_8_4_spect[97];//Skin L 75
		static const double ColabSkin75_10_33_spect[97];//Skin L 75
		static const double ColabSkin65_33_11_spect[97];//Skin L 65
		static const double ColabSkin65_7_24_spect[97];//Skin L 65
		static const double ColabSkin57_19_6_spect[97];//Skin L 57
		static const double ColabSkin57_4_19_spect[97];//Skin L 57
		static const double ColabSkin57_10_28_spect[97];//Skin L 57
		static const double ColabSkin40_7_19_spect[97];//Skin L 57
		static const double ColabSkin40_17_6_spect[97];//Skin L 40
		static const double ColabSkin40_4_11_spect[97];//Skin L 40
		static const double ColabSkin33_6_15_spect[97];//Skin L 33
		static const double ColabSkin33_15_5_spect[97];//Skin L 33
		static const double ColabSkin33_10_15_spect[97];//Skin L 33
		static const double ColabSkin24_5_6_spect[97];//Skin L 24
		static const double ColabSkin26_18_18_spect[97];//Skin L 26
		static const double ColabSkin24_7_5_spect[97];//Skin L 24
		static const double ColabSkin20_4_2_spect[97];//Skin L 20
		static const double ColabSkin98_m2_10_spect[97];//Skin L 98
		static const double ColabSkin90_m1_20_spect[97];//Skin L 90
		static const double ColabSkin95_0_4_spect[97];//Skin L 95
		static const double ColabSkin81_2_14_spect[97];//Skin L 81
		static const double ColabSkin87_3_10_spect[97];//Skin L 87
		static const double ColabSkin77_12_21_spect[97];//Skin L 77
		static const double ColabSkin70_7_32_spect[97];//Skin L 77
		static const double ColabSky60_0_m31_spect[97];//Sky L=60
		static const double ColabSky42_0_m24_spect[97];//Sky L=42

		static void spectrum_to_xyz_daylight  (double _m1, double _m2, double _temp, double &x, double &y, double &z);
		static void spectrum_to_xyz_blackbody (double _m1, double _m2, double _temp, double &x, double &y, double &z);
		static void spectrum_to_xyz_preset    (const double* spec_intens, double &x, double &y, double &z);

		static void spectrum_to_color_xyz_daylight  (const double* spec_color, double _m1, double _m2, double _temp, double &xx, double &yy, double &zz);
		static void spectrum_to_color_xyz_blackbody (const double* spec_color, double _m1, double _m2, double _temp, double &xx, double &yy, double &zz);
		static void spectrum_to_color_xyz_preset    (const double* spec_color, const double* spec_intens, double &xx, double &yy, double &zz);



/**
 * Inverse transform from CIECAM02 JCh to XYZ.
 */
		static void jch2xyz_ciecam02( double &x, double &y, double &z,
		                              double J, double C, double h,
		                              double xw, double yw, double zw,
		                              double yb, double la,
		                              double f, double c, double nc, bool doneinit2, int gamu );
									  
		static void jch2xyz_ciecam02float( float &x, float &y, float &z,
		                              float J, float C, float h,
		                              float xw, float yw, float zw,
		                              float yb, float la,
		                              float f, float c, float nc, bool doneinit2, int gamu );
									  
/**
 * Forward transform from XYZ to CIECAM02 JCh.
 */
		static void xyz2jchqms_ciecam02( double &J, double &C, double &h,
		                                 double &Q, double &M, double &s,double &aw, double &fl, double &wh,
		                                 double x, double y, double z,
		                                 double xw, double yw, double zw,
		                                 double yb, double la,
		                                 double f, double c, double nc,  double pilotd, bool doneinit1, int gamu  );

		static void xyz2jchqms_ciecam02float( float &J, float &C, float &h,
		                                 float &Q, float &M, float &s,float &aw, float &fl, float &wh,
		                                 float x, float y, float z,
		                                 float xw, float yw, float zw,
		                                 float yb, float la,
		                                 float f, float c, float nc,  float pilotd, bool doneinit1, int gamu  );
										 
										 
};
}
#endif