diff --git a/rtengine/CMakeLists.txt b/rtengine/CMakeLists.txt
index 9ebcce5ea..fad3c3cbc 100644
--- a/rtengine/CMakeLists.txt
+++ b/rtengine/CMakeLists.txt
@@ -5,7 +5,7 @@ link_directories ("${PROJECT_SOURCE_DIR}/rtexif" ${EXTRA_LIBDIR} ${GTHREAD_LIBRA
     ${GOBJECT_LIBRARY_DIRS} ${GLIB2_LIBRARY_DIRS} ${GLIBMM_LIBRARY_DIRS}
     ${IPTCDATA_LIBRARY_DIRS} ${LCMS_LIBRARY_DIRS} ${EXPAT_LIBRARY_DIRS})
 
-set (RTENGINESOURCEFILES safegtk.cc colortemp.cc curves.cc flatcurves.cc diagonalcurves.cc dcraw.cc iccstore.cc
+set (RTENGINESOURCEFILES safegtk.cc colortemp.cc curves.cc flatcurves.cc diagonalcurves.cc dcraw.cc iccstore.cc color.cc
     dfmanager.cc ffmanager.cc rawimage.cc image8.cc image16.cc imagefloat.cc imagedata.cc imageio.cc improcfun.cc init.cc dcrop.cc
     loadinitial.cc procparams.cc rawimagesource.cc demosaic_algos.cc shmap.cc simpleprocess.cc refreshmap.cc
     stdimagesource.cc myfile.cc iccjpeg.cc hlmultipliers.cc improccoordinator.cc
diff --git a/rtengine/color.cc b/rtengine/color.cc
index a7b8be1b4..4c95ee5b2 100644
--- a/rtengine/color.cc
+++ b/rtengine/color.cc
@@ -1,21 +1,21 @@
 /*
- *  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/>.
- */
+*  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/>.
+*/
 
 #include "rt_math.h"
 #include "color.h"
@@ -25,330 +25,362 @@ using namespace std;
 
 namespace rtengine {
 
-#define eps_max 580.40756 //(MAXVAL* 216.0f/24389.0);
-#define kappa	903.29630 //24389.0/27.0;
+    LUTf Color::cachef;
+    LUTf Color::gamma2curve = 0;
 
-LUTf Color::cachef ;
-LUTf Color::gamma2curve = 0;
+    LUTf Color::gammatab;
+    LUTf Color::igammatab_srgb;
+    LUTf Color::gammatab_srgb;
 
-LUTf Color::gammatab;
-LUTf Color::igammatab_srgb;
-LUTf Color::gammatab_srgb;
+    // Wikipedia sRGB: Unlike most other RGB color spaces, the sRGB gamma cannot be expressed as a single numerical value.
+    // The overall gamma is approximately 2.2, consisting of a linear (gamma 1.0) section near black, and a non-linear section elsewhere involving a 2.4 exponent
+    // and a gamma (slope of log output versus log input) changing from 1.0 through about 2.3.
+    const double Color::sRGBGamma = 2.2;
+    const double Color::sRGBGammaCurve = 2.4;
 
-// Wikipedia sRGB: Unlike most other RGB color spaces, the sRGB gamma cannot be expressed as a single numerical value.
-// The overall gamma is approximately 2.2, consisting of a linear (gamma 1.0) section near black, and a non-linear section elsewhere involving a 2.4 exponent
-// and a gamma (slope of log output versus log input) changing from 1.0 through about 2.3.
-const double Color::sRGBGamma = 2.2;
-const double Color::sRGBGammaCurve = 2.4;
+    const double Color::eps_max=580.40756;  //(MAXVAL* 216.0f/24389.0);
+    const double Color::kappa=903.29630;  //24389.0/27.0;
 
-void Color::init () {
+    const float Color::D50x=0.96422;
+    const float Color::D50z=0.82521;
+    const double Color::u0=4.0*D50x/(D50x+15+3*D50z);
+    const double Color::v0=9.0/(D50x+15+3*D50z);
 
-	int maxindex = 65536;
-	cachef(maxindex,0/*LUT_CLIP_BELOW*/);
+    void Color::init () {
 
-	gamma2curve(maxindex,0);
+        int maxindex = 65536;
+        cachef(maxindex,0/*LUT_CLIP_BELOW*/);
 
-	for (int i=0; i<maxindex; i++) {
-		if (i>eps_max) {
-			cachef[i] = 327.68*( exp(1.0/3.0 * log((double)i / MAXVAL) ));
-		}
-		else {
-			cachef[i] = 327.68*((kappa*i/MAXVAL+16.0)/116.0);
-		}
-	}
+        gamma2curve(maxindex,0);
 
-	for (int i=0; i<maxindex; i++) {
-		gamma2curve[i] = (gamma2(i/65535.0) * 65535.0);
-	}
+        for (int i=0; i<maxindex; i++) {
+            if (i>eps_max) {
+                cachef[i] = 327.68*( exp(1.0/3.0 * log((double)i / MAXVAL) ));
+            }
+            else {
+                cachef[i] = 327.68*((kappa*i/MAXVAL+16.0)/116.0);
+            }
+        }
 
-	/*******************************************/
+        for (int i=0; i<maxindex; i++) {
+            gamma2curve[i] = (gamma2(i/65535.0) * 65535.0);
+        }
 
-	gammatab(65536,0);
-	igammatab_srgb(65536,0);
-	gammatab_srgb(65536,0);
+        /*******************************************/
 
-	for (int i=0; i<65536; i++)
-	gammatab_srgb[i] = (65535.0 * gamma2 (i/65535.0));
-	for (int i=0; i<65536; i++)
-	igammatab_srgb[i] = (65535.0 * igamma2 (i/65535.0));
-	for (int i=0; i<65536; i++)
-	gammatab[i] = (65535.0 * pow (i/65535.0, 0.454545));
+        gammatab(65536,0);
+        igammatab_srgb(65536,0);
+        gammatab_srgb(65536,0);
 
-	/*FILE* f = fopen ("c.txt", "wt");
-	for (int i=0; i<256; i++)
-		fprintf (f, "%g %g\n", i/255.0, clower (i/255.0, 2.0, 1.0));
-	fclose (f);*/
-
-}
-
-void Color::cleanup () {
-
-}
-
-void Color::rgb2hsv (float r, float g, float b, float &h, float &s, float &v) {
-
-	double var_R = r / 65535.0;
-	double var_G = g / 65535.0;
-	double var_B = b / 65535.0;
-
-	double var_Min = min(var_R,var_G,var_B);
-	double var_Max = max(var_R,var_G,var_B);
-	double del_Max = var_Max - var_Min;
-	v = var_Max;
-	if (fabs(del_Max)<0.00001) {
-		h = 0;
-		s = 0;
-	}
-	else {
-		s = del_Max/var_Max;
-
-		if      ( var_R == var_Max ) h = (var_G - var_B)/del_Max;
-		else if ( var_G == var_Max ) h = 2.0 + (var_B - var_R)/del_Max;
-		else if ( var_B == var_Max ) h = 4.0 + (var_R - var_G)/del_Max;
-		h /= 6.0;
-
-		if ( h < 0 )  h += 1;
-		if ( h > 1 )  h -= 1;
-	}
-}
-
-void Color::rgb2hsv (int r, int g, int b, float &h, float &s, float &v) {
-
-	double var_R = r / 65535.0;
-	double var_G = g / 65535.0;
-	double var_B = b / 65535.0;
-
-	double var_Min = min(var_R,var_G,var_B);
-	double var_Max = max(var_R,var_G,var_B);
-	double del_Max = var_Max - var_Min;
-	v = var_Max;
-	if (fabs(del_Max)<0.00001) {
-		h = 0;
-		s = 0;
-	}
-	else {
-		s = del_Max/var_Max;
-
-		if      ( var_R == var_Max ) h = (var_G - var_B)/del_Max;
-		else if ( var_G == var_Max ) h = 2.0 + (var_B - var_R)/del_Max;
-		else if ( var_B == var_Max ) h = 4.0 + (var_R - var_G)/del_Max;
-		h /= 6.0;
-
-		if ( h < 0 )  h += 1;
-		if ( h > 1 )  h -= 1;
-	}
-}
-
-void Color::hsv2rgb (float h, float s, float v, float &r, float &g, float &b) {
-
-	float h1 = h*6; // sector 0 to 5
-	int i = floor( h1 );
-	float f = h1 - i; // fractional part of h
-
-	float p = v * ( 1 - s );
-	float q = v * ( 1 - s * f );
-	float t = v * ( 1 - s * ( 1 - f ) );
-
-	float r1,g1,b1;
-
-	if      (i==0) {r1 = v;  g1 = t;  b1 = p;}
-	else if (i==1) {r1 = q;  g1 = v;  b1 = p;}
-	else if (i==2) {r1 = p;  g1 = v;  b1 = t;}
-	else if (i==3) {r1 = p;  g1 = q;  b1 = v;}
-	else if (i==4) {r1 = t;  g1 = p;  b1 = v;}
-	else if (i==5) {r1 = v;  g1 = p;  b1 = q;}
-
-	r = ((r1)*65535.0);
-	g = ((g1)*65535.0);
-	b = ((b1)*65535.0);
-}
-
-
-// The same function but set float values instead if int
-// Function copied for speed concerns
-/* Not exactly the same as above ; this one return a result in the [0.0 ; 1.0] range
-void Color::hsv2rgb (float h, float s, float v, float &r, float &g, float &b) {
-
-	float h1 = h*6; // sector 0 to 5
-	int i = floor( h1 );
-	float f = h1 - i; // fractional part of h
-
-	float p = v * ( 1 - s );
-	float q = v * ( 1 - s * f );
-	float t = v * ( 1 - s * ( 1 - f ) );
-
-	if      (i==0) {r = v;  g = t;  b = p;}
-	else if (i==1) {r = q;  g = v;  b = p;}
-	else if (i==2) {r = p;  g = v;  b = t;}
-	else if (i==3) {r = p;  g = q;  b = v;}
-	else if (i==4) {r = t;  g = p;  b = v;}
-	else if (i==5) {r = v;  g = p;  b = q;}
-}
-*/
-
-void Color::hsv2rgb (float h, float s, float v, int &r, int &g, int &b) {
-
-	float h1 = h*6; // sector 0 to 5
-	int i = floor( h1 );
-	float f = h1 - i; // fractional part of h
-
-	float p = v * ( 1 - s );
-	float q = v * ( 1 - s * f );
-	float t = v * ( 1 - s * ( 1 - f ) );
-
-	float r1,g1,b1;
-
-	if (i==0) {r1 = v;  g1 = t;  b1 = p;}
-	else if (i==1) {r1 = q;  g1 = v;  b1 = p;}
-	else if (i==2) {r1 = p;  g1 = v;  b1 = t;}
-	else if (i==3) {r1 = p;  g1 = q;  b1 = v;}
-	else if (i==4) {r1 = t;  g1 = p;  b1 = v;}
-	else if (i==5) {r1 = v;  g1 = p;  b1 = q;}
-
-	r = (int)( r1 * 65535);
-	g = (int)( g1 * 65535);
-	b = (int)( b1 * 65535);
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-void Color::xyz2srgb (float x, float y, float z, float &r, float &g, float &b) {
-
-	//Transform to output color.  Standard sRGB is D65, but internal representation is D50
-	//Note that it is only at this point that we should have need of clipping color data
-
-	/*float x65 = d65_d50[0][0]*x + d65_d50[0][1]*y + d65_d50[0][2]*z ;
-	 float y65 = d65_d50[1][0]*x + d65_d50[1][1]*y + d65_d50[1][2]*z ;
-	 float z65 = d65_d50[2][0]*x + d65_d50[2][1]*y + d65_d50[2][2]*z ;
-
-	 r = sRGB_xyz[0][0]*x65 + sRGB_xyz[0][1]*y65 + sRGB_xyz[0][2]*z65;
-	 g = sRGB_xyz[1][0]*x65 + sRGB_xyz[1][1]*y65 + sRGB_xyz[1][2]*z65;
-	 b = sRGB_xyz[2][0]*x65 + sRGB_xyz[2][1]*y65 + sRGB_xyz[2][2]*z65;*/
-
-	/*r = sRGBd65_xyz[0][0]*x + sRGBd65_xyz[0][1]*y + sRGBd65_xyz[0][2]*z ;
-	 g = sRGBd65_xyz[1][0]*x + sRGBd65_xyz[1][1]*y + sRGBd65_xyz[1][2]*z ;
-	 b = sRGBd65_xyz[2][0]*x + sRGBd65_xyz[2][1]*y + sRGBd65_xyz[2][2]*z ;*/
-
-	r = ((sRGB_xyz[0][0]*x + sRGB_xyz[0][1]*y + sRGB_xyz[0][2]*z)) ;
-	g = ((sRGB_xyz[1][0]*x + sRGB_xyz[1][1]*y + sRGB_xyz[1][2]*z)) ;
-	b = ((sRGB_xyz[2][0]*x + sRGB_xyz[2][1]*y + sRGB_xyz[2][2]*z)) ;
-
-}
-
-
-void Color::xyz2rgb (float x, float y, float z, float &r, float &g, float &b, float rgb_xyz[3][3]) {
-
-	//Transform to output color.  Standard sRGB is D65, but internal representation is D50
-	//Note that it is only at this point that we should have need of clipping color data
-
-	/*float x65 = d65_d50[0][0]*x + d65_d50[0][1]*y + d65_d50[0][2]*z ;
-	 float y65 = d65_d50[1][0]*x + d65_d50[1][1]*y + d65_d50[1][2]*z ;
-	 float z65 = d65_d50[2][0]*x + d65_d50[2][1]*y + d65_d50[2][2]*z ;
-
-	 r = sRGB_xyz[0][0]*x65 + sRGB_xyz[0][1]*y65 + sRGB_xyz[0][2]*z65;
-	 g = sRGB_xyz[1][0]*x65 + sRGB_xyz[1][1]*y65 + sRGB_xyz[1][2]*z65;
-	 b = sRGB_xyz[2][0]*x65 + sRGB_xyz[2][1]*y65 + sRGB_xyz[2][2]*z65;*/
-
-	/*r = sRGBd65_xyz[0][0]*x + sRGBd65_xyz[0][1]*y + sRGBd65_xyz[0][2]*z ;
-	 g = sRGBd65_xyz[1][0]*x + sRGBd65_xyz[1][1]*y + sRGBd65_xyz[1][2]*z ;
-	 b = sRGBd65_xyz[2][0]*x + sRGBd65_xyz[2][1]*y + sRGBd65_xyz[2][2]*z ;*/
-
-
-
-	r = ((rgb_xyz[0][0]*x + rgb_xyz[0][1]*y + rgb_xyz[0][2]*z)) ;
-	g = ((rgb_xyz[1][0]*x + rgb_xyz[1][1]*y + rgb_xyz[1][2]*z)) ;
-	b = ((rgb_xyz[2][0]*x + rgb_xyz[2][1]*y + rgb_xyz[2][2]*z)) ;
-
-}
-
-void Color::calcGamma (double pwr, double ts, int mode, int imax, double &gamma0, double &gamma1, double &gamma2, double &gamma3, double &gamma4, double &gamma5) {
-//from Dcraw (D.Coffin)
-	int i;
-	double g[6], bnd[2]={0,0};
-
-	g[0] = pwr;
-	g[1] = ts;
-	g[2] = g[3] = g[4] = 0;
-	bnd[g[1] >= 1] = 1;
-	if (g[1] && (g[1]-1)*(g[0]-1) <= 0) {
-		for (i=0; i < 48; i++) {
-			g[2] = (bnd[0] + bnd[1])/2;
-			if (g[0])
-				bnd[(pow(g[2]/g[1],-g[0]) - 1)/g[0] - 1/g[2] > -1] = g[2];
-			else
-				bnd[g[2]/exp(1-1/g[2]) < g[1]] = g[2];
-		}
-		g[3] = g[2] / g[1];
-		if (g[0]) g[4] = g[2] * (1/g[0] - 1);
-	}
-	if (g[0])
-		g[5] = 1 / (g[1]*SQR(g[3])/2 - g[4]*(1 - g[3]) + (1 - pow(g[3],1+g[0]))*(1 + g[4])/(1 + g[0])) - 1;
-	else
-		g[5] = 1 / (g[1]*SQR(g[3])/2 + 1 - g[2] - g[3] - g[2]*g[3]*(log(g[3]) - 1)) - 1;
-	if (!mode--) {
-		gamma0=g[0];gamma1=g[1];gamma2=g[2];gamma3=g[3];gamma4=g[4];gamma5=g[5];
-		return;
-	}
-}
-
-void Color::Lab2XYZ(float L, float a, float b, float &x, float &y, float &z) {
-	float fy = (0.00862069 * L) + 0.137932; // (L+16)/116
-	float fx = (0.002 * a) + fy;
-	float fz = fy - (0.005 * b);
-
-	x = 65535.0*f2xyz(fx)*D50x;
-	y = 65535.0*f2xyz(fy);
-	z = 65535.0*f2xyz(fz)*D50z;
-}
-
-void Color::XYZ2Lab(float X, float Y, float Z, float &L, float &a, float &b) {
-
-	float X1 = X/D50x;
-	float Z1 = Z/D50z;
-
-	float fx = (X1<65535.0 ? cachef[X1] : (327.68*exp(log(X1/MAXVAL)/3.0 )));
-	float fy = (Y<65535.0 ? cachef[Y] : (327.68*exp(log(Y/MAXVAL)/3.0 )));
-	float fz = (Z1<65535.0 ? cachef[Z1] : (327.68*exp(log(Z1/MAXVAL)/3.0 )));
-
-	L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68;
-	a = (500.0 * (fx - fy) );
-	b = (200.0 * (fy - fz) );
-
-}
-
-void Color::Lab2Yuv(float L, float a, float b, float &Y, float &u, float &v) {
-	float fy = (0.00862069 * L/327.68) + 0.137932; // (L+16)/116
-	float fx = (0.002 * a/327.68) + fy;
-	float fz = fy - (0.005 * b/327.68);
-
-	float X = 65535.0*f2xyz(fx)*D50x;
-	Y = 65535.0*f2xyz(fy);
-	float Z = 65535.0*f2xyz(fz)*D50z;
-
-	u = 4.0*X/(X+15*Y+3*Z)-u0;
-	v = 9.0*Y/(X+15*Y+3*Z)-v0;
-
-}
-
-void Color::Yuv2Lab(float Yin, float u, float v, float &L, float &a, float &b, double wp[3][3]) {
-
-	float u1 = u + u0;
-	float v1 = v + v0;
-
-	float Y = Yin;
-	float X = (9*u1*Y)/(4*v1*D50x);
-	float Z = (12 - 3*u1 - 20*v1)*Y/(4*v1*D50z);
-
-	gamutmap(X,Y,Z,wp);
-
-	float fx = (X<65535.0 ? cachef[X] : (327.68*exp(log(X/MAXVAL)/3.0 )));
-	float fy = (Y<65535.0 ? cachef[Y] : (327.68*exp(log(Y/MAXVAL)/3.0 )));
-	float fz = (Z<65535.0 ? cachef[Z] : (327.68*exp(log(Z/MAXVAL)/3.0 )));
-
-	L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68;
-	a = (500.0 * (fx - fy) );
-	b = (200.0 * (fy - fz) );
-
-}
+        for (int i=0; i<65536; i++)
+            gammatab_srgb[i] = (65535.0 * gamma2 (i/65535.0));
+        for (int i=0; i<65536; i++)
+            igammatab_srgb[i] = (65535.0 * igamma2 (i/65535.0));
+        for (int i=0; i<65536; i++)
+            gammatab[i] = (65535.0 * pow (i/65535.0, 0.454545));
 
+        /*FILE* f = fopen ("c.txt", "wt");
+        for (int i=0; i<256; i++)
+        fprintf (f, "%g %g\n", i/255.0, clower (i/255.0, 2.0, 1.0));
+        fclose (f);*/
+
+    }
+
+    void Color::cleanup () {
+
+    }
+
+    void Color::rgb2hsv(float r, float g, float b, float &h, float &s, float &v) {
+        double var_R = r / 65535.0;
+        double var_G = g / 65535.0;
+        double var_B = b / 65535.0;
+
+        double var_Min = min(var_R,var_G,var_B);
+        double var_Max = max(var_R,var_G,var_B);
+        double del_Max = var_Max - var_Min;
+        v = var_Max;
+        if (del_Max<0.00001 && del_Max>-0.00001) {  // no fabs, slow!
+            h = 0;
+            s = 0;
+        }
+        else {
+            s = del_Max/var_Max;
+
+            if      ( var_R == var_Max ) h = (var_G - var_B)/del_Max; 
+            else if ( var_G == var_Max ) h = 2.0 + (var_B - var_R)/del_Max; 
+            else if ( var_B == var_Max ) h = 4.0 + (var_R - var_G)/del_Max; 
+            h /= 6.0;
+
+            if ( h < 0 )  h += 1;
+            if ( h > 1 )  h -= 1;
+        }
+    }
+
+    void Color::hsv2rgb (float h, float s, float v, float &r, float &g, float &b) {
+
+        float h1 = h*6; // sector 0 to 5
+        int i = (int)h1;  // floor() is very slow, and h1 is always >0
+        float f = h1 - i; // fractional part of h
+
+        float p = v * ( 1 - s );
+        float q = v * ( 1 - s * f );
+        float t = v * ( 1 - s * ( 1 - f ) );
+
+        float r1,g1,b1;
+
+        if      (i==1)    {r1 = q;  g1 = v;  b1 = p;}
+        else if (i==2)    {r1 = p;  g1 = v;  b1 = t;}
+        else if (i==3)    {r1 = p;  g1 = q;  b1 = v;}
+        else if (i==4)    {r1 = t;  g1 = p;  b1 = v;}
+        else if (i==5)    {r1 = v;  g1 = p;  b1 = q;}
+        else /*i==(0|6)*/ {r1 = v;  g1 = t;  b1 = p;}
+
+        r = ((r1)*65535.0);
+        g = ((g1)*65535.0);
+        b = ((b1)*65535.0);
+    }
+
+    // Function copied for speed concerns
+    // Not exactly the same as above ; this one return a result in the [0.0 ; 1.0] range
+    void Color::hsv2rgb01 (float h, float s, float v, float &r, float &g, float &b) {
+        float h1 = h*6; // sector 0 to 5
+        int i = int(h1);
+        float f = h1 - i; // fractional part of h
+
+        float p = v * ( 1 - s );
+        float q = v * ( 1 - s * f );
+        float t = v * ( 1 - s * ( 1 - f ) );
+
+        if      (i==1)    {r = q;  g = v;  b = p;}
+        else if (i==2)    {r = p;  g = v;  b = t;}
+        else if (i==3)    {r = p;  g = q;  b = v;}
+        else if (i==4)    {r = t;  g = p;  b = v;}
+        else if (i==5)    {r = v;  g = p;  b = q;}
+        else /*(i==0|6)*/ {r = v;  g = t;  b = p;}
+    }
+
+    void Color::hsv2rgb (float h, float s, float v, int &r, int &g, int &b) {
+
+        float h1 = h*6; // sector 0 to 5
+        int i = floor( h1 );
+        float f = h1 - i; // fractional part of h
+
+        float p = v * ( 1 - s );
+        float q = v * ( 1 - s * f );
+        float t = v * ( 1 - s * ( 1 - f ) );
+
+        float r1,g1,b1;
+
+        if (i==0) {r1 = v;  g1 = t;  b1 = p;}
+        else if (i==1) {r1 = q;  g1 = v;  b1 = p;}
+        else if (i==2) {r1 = p;  g1 = v;  b1 = t;}
+        else if (i==3) {r1 = p;  g1 = q;  b1 = v;}
+        else if (i==4) {r1 = t;  g1 = p;  b1 = v;}
+        else if (i==5) {r1 = v;  g1 = p;  b1 = q;}
+
+        r = (int)( r1 * 65535);
+        g = (int)( g1 * 65535);
+        b = (int)( b1 * 65535);
+    }
+
+    //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    void Color::xyz2srgb (float x, float y, float z, float &r, float &g, float &b) {
+
+        //Transform to output color.  Standard sRGB is D65, but internal representation is D50
+        //Note that it is only at this point that we should have need of clipping color data
+
+        /*float x65 = d65_d50[0][0]*x + d65_d50[0][1]*y + d65_d50[0][2]*z ;
+        float y65 = d65_d50[1][0]*x + d65_d50[1][1]*y + d65_d50[1][2]*z ;
+        float z65 = d65_d50[2][0]*x + d65_d50[2][1]*y + d65_d50[2][2]*z ;
+
+        r = sRGB_xyz[0][0]*x65 + sRGB_xyz[0][1]*y65 + sRGB_xyz[0][2]*z65;
+        g = sRGB_xyz[1][0]*x65 + sRGB_xyz[1][1]*y65 + sRGB_xyz[1][2]*z65;
+        b = sRGB_xyz[2][0]*x65 + sRGB_xyz[2][1]*y65 + sRGB_xyz[2][2]*z65;*/
+
+        /*r = sRGBd65_xyz[0][0]*x + sRGBd65_xyz[0][1]*y + sRGBd65_xyz[0][2]*z ;
+        g = sRGBd65_xyz[1][0]*x + sRGBd65_xyz[1][1]*y + sRGBd65_xyz[1][2]*z ;
+        b = sRGBd65_xyz[2][0]*x + sRGBd65_xyz[2][1]*y + sRGBd65_xyz[2][2]*z ;*/
+
+        r = ((sRGB_xyz[0][0]*x + sRGB_xyz[0][1]*y + sRGB_xyz[0][2]*z)) ;
+        g = ((sRGB_xyz[1][0]*x + sRGB_xyz[1][1]*y + sRGB_xyz[1][2]*z)) ;
+        b = ((sRGB_xyz[2][0]*x + sRGB_xyz[2][1]*y + sRGB_xyz[2][2]*z)) ;
+
+    }
+
+    void Color::xyz2rgb (float x, float y, float z, float &r, float &g, float &b, double rgb_xyz[3][3]) {
+        //Transform to output color.  Standard sRGB is D65, but internal representation is D50
+        //Note that it is only at this point that we should have need of clipping color data
+
+        /*float x65 = d65_d50[0][0]*x + d65_d50[0][1]*y + d65_d50[0][2]*z ;
+        float y65 = d65_d50[1][0]*x + d65_d50[1][1]*y + d65_d50[1][2]*z ;
+        float z65 = d65_d50[2][0]*x + d65_d50[2][1]*y + d65_d50[2][2]*z ;
+
+        r = sRGB_xyz[0][0]*x65 + sRGB_xyz[0][1]*y65 + sRGB_xyz[0][2]*z65;
+        g = sRGB_xyz[1][0]*x65 + sRGB_xyz[1][1]*y65 + sRGB_xyz[1][2]*z65;
+        b = sRGB_xyz[2][0]*x65 + sRGB_xyz[2][1]*y65 + sRGB_xyz[2][2]*z65;*/
+
+        /*r = sRGBd65_xyz[0][0]*x + sRGBd65_xyz[0][1]*y + sRGBd65_xyz[0][2]*z ;
+        g = sRGBd65_xyz[1][0]*x + sRGBd65_xyz[1][1]*y + sRGBd65_xyz[1][2]*z ;
+        b = sRGBd65_xyz[2][0]*x + sRGBd65_xyz[2][1]*y + sRGBd65_xyz[2][2]*z ;*/
+
+        r = ((rgb_xyz[0][0]*x + rgb_xyz[0][1]*y + rgb_xyz[0][2]*z)) ;
+        g = ((rgb_xyz[1][0]*x + rgb_xyz[1][1]*y + rgb_xyz[1][2]*z)) ;
+        b = ((rgb_xyz[2][0]*x + rgb_xyz[2][1]*y + rgb_xyz[2][2]*z)) ;
+    }
+
+    // same for float
+    void Color::xyz2rgb (float x, float y, float z, float &r, float &g, float &b, float rgb_xyz[3][3]) {
+        r = ((rgb_xyz[0][0]*x + rgb_xyz[0][1]*y + rgb_xyz[0][2]*z)) ;
+        g = ((rgb_xyz[1][0]*x + rgb_xyz[1][1]*y + rgb_xyz[1][2]*z)) ;
+        b = ((rgb_xyz[2][0]*x + rgb_xyz[2][1]*y + rgb_xyz[2][2]*z)) ;
+    }
+
+    void Color::calcGamma (double pwr, double ts, int mode, int imax, double &gamma0, double &gamma1, double &gamma2, double &gamma3, double &gamma4, double &gamma5) {
+        //from Dcraw (D.Coffin)
+        int i;
+        double g[6], bnd[2]={0,0};
+
+        g[0] = pwr;
+        g[1] = ts;
+        g[2] = g[3] = g[4] = 0;
+        bnd[g[1] >= 1] = 1;
+        if (g[1] && (g[1]-1)*(g[0]-1) <= 0) {
+            for (i=0; i < 48; i++) {
+                g[2] = (bnd[0] + bnd[1])/2;
+                if (g[0])
+                    bnd[(pow(g[2]/g[1],-g[0]) - 1)/g[0] - 1/g[2] > -1] = g[2];
+                else
+                    bnd[g[2]/exp(1-1/g[2]) < g[1]] = g[2];
+            }
+            g[3] = g[2] / g[1];
+            if (g[0]) g[4] = g[2] * (1/g[0] - 1);
+        }
+        if (g[0])
+            g[5] = 1 / (g[1]*SQR(g[3])/2 - g[4]*(1 - g[3]) + (1 - pow(g[3],1+g[0]))*(1 + g[4])/(1 + g[0])) - 1;
+        else
+            g[5] = 1 / (g[1]*SQR(g[3])/2 + 1 - g[2] - g[3] - g[2]*g[3]*(log(g[3]) - 1)) - 1;
+        if (!mode--) {
+            gamma0=g[0];gamma1=g[1];gamma2=g[2];gamma3=g[3];gamma4=g[4];gamma5=g[5];
+            return;
+        }
+    }
+
+    void Color::Lab2XYZ(float L, float a, float b, float &x, float &y, float &z) {
+        float fy = (0.00862069 * L) + 0.137932; // (L+16)/116
+        float fx = (0.002 * a) + fy;
+        float fz = fy - (0.005 * b);
+
+        x = 65535.0*f2xyz(fx)*D50x;
+        y = 65535.0*f2xyz(fy);
+        z = 65535.0*f2xyz(fz)*D50z;
+    }
+
+    void Color::XYZ2Lab(float X, float Y, float Z, float &L, float &a, float &b) {
+
+        float X1 = X/D50x;
+        float Z1 = Z/D50z;
+
+        float fx = (X1<65535.0 ? cachef[X1] : (327.68*exp(log(X1/MAXVAL)/3.0 )));
+        float fy = (Y<65535.0 ? cachef[Y] : (327.68*exp(log(Y/MAXVAL)/3.0 )));
+        float fz = (Z1<65535.0 ? cachef[Z1] : (327.68*exp(log(Z1/MAXVAL)/3.0 )));
+
+        L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68;
+        a = (500.0 * (fx - fy) );
+        b = (200.0 * (fy - fz) );
+
+    }
+
+    void Color::Lab2Yuv(float L, float a, float b, float &Y, float &u, float &v) {
+        float fy = (0.00862069 * L/327.68) + 0.137932; // (L+16)/116
+        float fx = (0.002 * a/327.68) + fy;
+        float fz = fy - (0.005 * b/327.68);
+
+        float X = 65535.0*f2xyz(fx)*D50x;
+        Y = 65535.0*f2xyz(fy);
+        float Z = 65535.0*f2xyz(fz)*D50z;
+
+        u = 4.0*X/(X+15*Y+3*Z)-u0;
+        v = 9.0*Y/(X+15*Y+3*Z)-v0;
+    }
+
+    void Color::Yuv2Lab(float Yin, float u, float v, float &L, float &a, float &b, double wp[3][3]) {
+        float u1 = u + u0;
+        float v1 = v + v0;
+
+        float Y = Yin;
+        float X = (9*u1*Y)/(4*v1*D50x);
+        float Z = (12 - 3*u1 - 20*v1)*Y/(4*v1*D50z);
+
+        gamutmap(X,Y,Z,wp);
+
+        float fx = (X<65535.0 ? cachef[X] : (327.68*exp(log(X/MAXVAL)/3.0 )));
+        float fy = (Y<65535.0 ? cachef[Y] : (327.68*exp(log(Y/MAXVAL)/3.0 )));
+        float fz = (Z<65535.0 ? cachef[Z] : (327.68*exp(log(Z/MAXVAL)/3.0 )));
+
+        L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68;
+        a = (500.0 * (fx - fy) );
+        b = (200.0 * (fy - fz) );
+    }
+
+    double Color::f2xyz(double f) {
+        const double epsilonExpInv3 = 6.0/29.0;
+        const double kappaInv = 27.0/24389.0;  // inverse of kappa
+
+        return (f > epsilonExpInv3) ? f*f*f : (116 * f - 16) * kappaInv;
+    }
+
+/* Gamut mapping algorithm
+ copyright (c) 2010-2011  Emil Martinec <ejmartin@uchicago.edu>
+ solutions to scaling u and v to XYZ paralleliped boundaries
+ some equations:
+ fu(X,Y,Z) = 4 X/(X + 15 Y + 3 Z);
+ fv(X,Y,Z) = 9 Y/(X + 15 Y + 3 Z);
+ 
+ take the plane spanned by X=a*Xr+b*Xg+c*Xb etc with one of a,b,c equal to 0 or 1,
+ and itersect with the line u0+lam*u, or in other words solve
+ 
+ u0+lam*u=fu(X,Y,Z)
+ v0+lam*v=fv(X,Y,Z)
+ 
+ the value of lam is the scale factor that takes the color to the gamut boundary
+
+ columns of the matrix p=xyz_rgb are RGB tristimulus primaries in XYZ
+ c is the color fixed on the boundary; and m=0 for c=0, m=1 for c=255
+*/
+
+    void Color::gamutmap(float &X, float &Y, float &Z, const double p[3][3]) 
+    {
+        float u = 4*X/(X+15*Y+3*Z)-u0;
+        float v = 9*Y/(X+15*Y+3*Z)-v0;
+
+        float lam[3][2];
+        float lam_min = 1.0;
+
+        for (int c=0; c<3; c++)
+            for (int m=0; m<2; m++) {
+
+                int c1=(c+1)%3;
+                int c2=(c+2)%3;
+
+                lam[c][m] = (-(p[0][c1]*p[1][c]*((-12 + 3*u0 + 20*v0)*Y + 4*m*65535*v0*p[2][c2])) +
+                    p[0][c]*p[1][c1]*((-12 + 3*u0 + 20*v0)*Y + 4*m*65535*v0*p[2][c2]) -
+                    4*v0*p[0][c1]*(Y - m*65535*p[1][c2])*p[2][c] + 4*v0*p[0][c]*(Y - m*65535*p[1][c2])*p[2][c1] -
+                    (4*m*65535*v0*p[0][c2] - 9*u0*Y)*(p[1][c1]*p[2][c] - p[1][c]*p[2][c1]));
+
+                lam[c][m] /= (3*u*Y*(p[0][c1]*p[1][c] - p[1][c1]*(p[0][c] + 3*p[2][c]) + 3*p[1][c]*p[2][c1]) +
+                    4*v*(p[0][c1]*(5*Y*p[1][c] + m*65535*p[1][c]*p[2][c2] + Y*p[2][c] - m*65535*p[1][c2]*p[2][c]) -
+                    p[0][c]*(5*Y*p[1][c1] + m*65535*p[1][c1]*p[2][c2] + Y*p[2][c1] - m*65535*p[1][c2]*p[2][c1]) +
+                    m*65535*p[0][c2]*(p[1][c1]*p[2][c] - p[1][c]*p[2][c1])));
+
+                if (lam[c][m]<lam_min && lam[c][m]>0) lam_min=lam[c][m];
+
+            }
+
+            u = u*lam_min + u0;
+            v = v*lam_min + v0;
+
+            X = (9*u*Y)/(4*v); 
+            Z = (12 - 3*u - 20*v)*Y/(4*v);
+    }
 }
diff --git a/rtengine/color.h b/rtengine/color.h
new file mode 100644
index 000000000..7b6a3b57a
--- /dev/null
+++ b/rtengine/color.h
@@ -0,0 +1,104 @@
+/*
+ *  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 _COLOR_H_
+#define _COLOR_H_
+
+#include <math.h>
+#include "LUT.h"
+
+namespace rtengine {
+
+class Color {
+
+public:
+	const static double sRGBGamma;        // standard average gamma
+    const static double sRGBGammaCurve;   // 2.4 in the curve
+    const static double eps_max, kappa;
+    const static float D50x, D50z;
+    const static double u0, v0;
+
+    static LUTf cachef;
+	static LUTf gamma2curve;
+
+    // look-up tables for the standard srgb gamma and its inverse (filled by init())
+    static LUTf igammatab_srgb;
+    static LUTf gammatab_srgb;
+    // look-up tables for the simple exponential gamma
+    static LUTf gammatab;
+
+
+	static void init ();
+	static void cleanup ();
+
+	static void rgb2hsv (float r, float g, float b, float &h, float &s, float &v);
+	static void hsv2rgb (float h, float s, float v, float &r, float &g, float &b);
+	static void hsv2rgb (float h, float s, float v, int &r, int &g, int &b);
+    static void hsv2rgb01 (float h, float s, float v, float &r, float &g, float &b);
+	static void xyz2srgb (float x, float y, float z, float &r, float &g, float &b);
+	static void xyz2rgb (float x, float y, float z, float &r, float &g, float &b, double rgb_xyz[3][3]);
+    static void xyz2rgb (float x, float y, float z, float &r, float &g, float &b, float rgb_xyz[3][3]);
+	static void Lab2XYZ(float L, float a, float b, float &x, float &y, float &z);
+	static void XYZ2Lab(float X, float Y, float Z, float &L, float &a, float &b);
+	static void Lab2Yuv(float L, float a, float b, float &Y, float &u, float &v);
+	static void Yuv2Lab(float Y, float u, float v, float &L, float &a, float &b, double wp[3][3]);
+    static double f2xyz(double f);
+	static void calcGamma (double pwr, double ts, int mode, int imax, double &gamma0, double &gamma1, double &gamma2, double &gamma3, double &gamma4,double &gamma5);
+
+	// standard srgb gamma and its inverse
+	static inline double gamma2     (double x) {
+											return x <= 0.00304 ? x*12.92 : 1.055*exp(log(x)/sRGBGammaCurve)-0.055;
+									}
+	static inline double igamma2    (double x) {
+										return x <= 0.03928 ? x/12.92 : exp(log((x+0.055)/1.055)*sRGBGammaCurve);
+									}
+	// gamma function with adjustable parameters
+	static inline double gamma      (double x, double gamma, double start, double slope, double mul, double add){
+										return (x <= start ? x*slope : exp(log(x)/gamma)*mul-add);
+									}
+	static inline double igamma     (double x, double gamma, double start, double slope, double mul, double add){
+										return (x <= start*slope ? x/slope : exp(log((x+add)/mul)*gamma) );
+									}
+
+	// gamma functions on [0,65535] based on look-up tables
+	static inline float  gamma_srgb       (int x) { return gammatab_srgb[x]; }
+	static inline float  gamma            (int x) { return gammatab[x]; }
+	static inline float  igamma_srgb      (int x) { return igammatab_srgb[x]; }
+	static inline float  gamma_srgb       (float x) { return gammatab_srgb[x]; }
+	static inline float  gamma            (float x) { return gammatab[x]; }
+	static inline float  igamma_srgb      (float x) { return igammatab_srgb[x]; }
+	//static inline float  gamma_srgb       (double x) { return gammatab_srgb[x]; }
+	//static inline float  gamma            (double x) { return gammatab[x]; }
+	//static inline float  igamma_srgb      (double x) { return igammatab_srgb[x]; }
+
+	//void gamutmap(LabImage* );
+	static void gamutmap(float &X, float &Y, float &Z, const double p[3][3]);
+
+	static inline float f2xyz(float f) {
+		const float epsilonExpInv3 = 6.0/29.0;
+		const float kappaInv = 27.0/24389.0;  // inverse of kappa
+
+		return (f > epsilonExpInv3) ? f*f*f : (116 * f - 16) * kappaInv;
+	}
+
+};
+
+}
+
+#endif
diff --git a/rtengine/dcp.cc b/rtengine/dcp.cc
index 60acca23f..a490b0090 100644
--- a/rtengine/dcp.cc
+++ b/rtengine/dcp.cc
@@ -271,12 +271,12 @@ void DCPProfile::Apply(Imagefloat *pImg, DCPLightType preferredProfile, Glib::us
 
                 // if point is in negative area, just the matrix, but not the LUT
                 if (newr>=0 && newg>=0 && newb>=0) {
-                ImProcFunctions::rgb2hsv(newr, newg, newb, h , s, v);
+                Color::rgb2hsv(newr, newg, newb, h , s, v);
                 h*=6.f;  // RT calculates in [0,1]
 
                 if (useRawWhite) {
                     // Retro-calculate what the point was like before RAW white came in
-                    ImProcFunctions::rgb2hsv(newr/rawWhiteFac, newg/rawWhiteFac, newb/rawWhiteFac, hs, ss, vs);
+                    Color::rgb2hsv(newr/rawWhiteFac, newg/rawWhiteFac, newb/rawWhiteFac, hs, ss, vs);
                     hs*=6.f;  // RT calculates in [0,1]
                 } else {
                     hs=h; ss=s; vs=v;
@@ -422,7 +422,7 @@ void DCPProfile::Apply(Imagefloat *pImg, DCPLightType preferredProfile, Glib::us
                 if (h < 0.0f) h += 6.0f;
                 if (h >= 6.0f) h -= 6.0f;
                 h/=6.f;  
-                ImProcFunctions::hsv2rgb( h, s, v, newr, newg, newb);
+                Color::hsv2rgb( h, s, v, newr, newg, newb);
                 }
 
                 pImg->r[y][x] = m2Work[0][0]*newr + m2Work[0][1]*newg + m2Work[0][2]*newb;
diff --git a/rtengine/dcrop.cc b/rtengine/dcrop.cc
index 84ea38cc8..6b52463e3 100644
--- a/rtengine/dcrop.cc
+++ b/rtengine/dcrop.cc
@@ -393,7 +393,7 @@ void Crop::fullUpdate () {
     parent->updaterThreadStart.lock ();
     if (parent->updaterRunning && parent->thread) {
 		// Do NOT reset changes here, since in a long chain of events it will lead to chroma_scale not being updated,
-		// causing ImProcFunctions::lab2rgb to return a black image on some opens
+		// causing Color::lab2rgb to return a black image on some opens
         //parent->changeSinceLast = 0;
         parent->thread->join ();
     }
diff --git a/rtengine/gamutbdy.cc b/rtengine/gamutbdy.cc
deleted file mode 100644
index 1b57b1cfb..000000000
--- a/rtengine/gamutbdy.cc
+++ /dev/null
@@ -1,97 +0,0 @@
-////////////////////////////////////////////////////////////////
-//
-//		Gamut mapping algorithm
-//
-//		copyright (c) 2010-2011  Emil Martinec <ejmartin@uchicago.edu>
-//
-//
-// code dated: February 2, 2011
-//
-//	sRGBgamutbdy.cc 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.
-//
-//	This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
-//
-////////////////////////////////////////////////////////////////
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-	
-#include "rt_math.h"
-
-#define u0 4.0*D50x/(D50x+15+3*D50z)
-#define v0 9.0/(D50x+15+3*D50z)
-
-#define sgn(x) (x<0 ? -1 : 1)
-#define sqrt2 (sqrt(2))
-#define sqrt3 (sqrt(3))
-#define sqrt6 (sqrt(6))
-#define K 20
-#define rmax	(3*K + 1)
-#define rctr	(2*K)
-#define cctr	(2*K + 1)
-#define maxindx	(2*(3*K+1) + SQR(3*K+1))
-
-
-
-// solutions to scaling u and v to XYZ paralleliped boundaries
-/* some equations:
- fu(X,Y,Z) = 4 X/(X + 15 Y + 3 Z);
- fv(X,Y,Z) = 9 Y/(X + 15 Y + 3 Z);
- 
- take the plane spanned by X=a*Xr+b*Xg+c*Xb etc with one of a,b,c equal to 0 or 1,
- and itersect with the line u0+lam*u, or in other words solve
- 
- u0+lam*u=fu(X,Y,Z)
- v0+lam*v=fv(X,Y,Z)
- 
- the value of lam is the scale factor that takes the color to the gamut boundary
-*/
-// columns of the matrix p=xyz_rgb are RGB tristimulus primaries in XYZ
-// c is the color fixed on the boundary; and m=0 for c=0, m=1 for c=255
-
-void ImProcFunctions::gamutmap(float &X, float &Y, float &Z, const double p[3][3]) 
-{
-	float u = 4*X/(X+15*Y+3*Z)-u0;
-	float v = 9*Y/(X+15*Y+3*Z)-v0;
-	
-	float lam[3][2];
-	float lam_min = 1.0;
-	
-	for (int c=0; c<3; c++)
-		for (int m=0; m<2; m++) {
-			
-			int c1=(c+1)%3;
-			int c2=(c+2)%3;
-			
-			lam[c][m] = (-(p[0][c1]*p[1][c]*((-12 + 3*u0 + 20*v0)*Y + 4*m*65535*v0*p[2][c2])) +
-						 p[0][c]*p[1][c1]*((-12 + 3*u0 + 20*v0)*Y + 4*m*65535*v0*p[2][c2]) -
-						 4*v0*p[0][c1]*(Y - m*65535*p[1][c2])*p[2][c] + 4*v0*p[0][c]*(Y - m*65535*p[1][c2])*p[2][c1] -
-						 (4*m*65535*v0*p[0][c2] - 9*u0*Y)*(p[1][c1]*p[2][c] - p[1][c]*p[2][c1]));
-			
-			lam[c][m] /= (3*u*Y*(p[0][c1]*p[1][c] - p[1][c1]*(p[0][c] + 3*p[2][c]) + 3*p[1][c]*p[2][c1]) +
-					4*v*(p[0][c1]*(5*Y*p[1][c] + m*65535*p[1][c]*p[2][c2] + Y*p[2][c] - m*65535*p[1][c2]*p[2][c]) -
-						 p[0][c]*(5*Y*p[1][c1] + m*65535*p[1][c1]*p[2][c2] + Y*p[2][c1] - m*65535*p[1][c2]*p[2][c1]) +
-						 m*65535*p[0][c2]*(p[1][c1]*p[2][c] - p[1][c]*p[2][c1])));
-			
-			if (lam[c][m]<lam_min && lam[c][m]>0) lam_min=lam[c][m];
-
-		}
-	
-	u = u*lam_min + u0;
-	v = v*lam_min + v0;
-	
-	X = (9*u*Y)/(4*v); 
-	Z = (12 - 3*u - 20*v)*Y/(4*v);
-	
-	
-}
-
-	
-//};//namespace
diff --git a/rtengine/improcfun.cc b/rtengine/improcfun.cc
index 5b4d13664..c9cfc4a98 100644
--- a/rtengine/improcfun.cc
+++ b/rtengine/improcfun.cc
@@ -33,6 +33,7 @@
 #include "rtthumbnail.h"
 #include "utils.h"
 #include "iccmatrices.h"
+#include "color.h"
 #include "calc_distort.h"
 #include "rt_math.h"
 
@@ -54,12 +55,6 @@ namespace rtengine {
 #define CLIP2(a) ((a)<MAXVAL ? a : MAXVAL )
 #define FCLIP(a) ((a)>0.0?((a)<65535.5?(a):65535.5):0.0)
 	
-#define u0 4.0*D50x/(D50x+15+3*D50z)
-#define v0 9.0/(D50x+15+3*D50z)
-	
-#define eps_max 580.40756 //(MAXVAL* 216.0f/24389.0);
-#define kappa	903.29630 //24389.0/27.0;
-	
 	
 extern const Settings* settings;
 
@@ -74,11 +69,11 @@ void ImProcFunctions::initCache () {
 	gamma2curve(maxindex,0);
 
     for (int i=0; i<maxindex; i++) {
-        if (i>eps_max) {
+        if (i>Color::eps_max) {
 			cachef[i] = 327.68*( exp(1.0/3.0 * log((double)i / MAXVAL) ));
         }
         else {
-			cachef[i] = 327.68*((kappa*i/MAXVAL+16.0)/116.0);
+			cachef[i] = 327.68*((Color::kappa*i/MAXVAL+16.0)/116.0);
         }
 	}
 
@@ -209,17 +204,17 @@ void ImProcFunctions::rgbProc (Imagefloat* working, LabImage* lab, LUTf & hltone
 
     double toxyz[3][3] = {
         {
-        	( wprof[0][0] / D50x),
-        	( wprof[0][1] / D50x),
-        	( wprof[0][2] / D50x)
+        	( wprof[0][0] / Color::D50x),
+        	( wprof[0][1] / Color::D50x),
+        	( wprof[0][2] / Color::D50x)
         },{
 			( wprof[1][0]		),
 			( wprof[1][1]		),
 			( wprof[1][2]		)
         },{
-			( wprof[2][0] / D50z),
-			( wprof[2][1] / D50z),
-			( wprof[2][2] / D50z)
+			( wprof[2][0] / Color::D50z),
+			( wprof[2][1] / Color::D50z),
+			( wprof[2][2] / Color::D50z)
         }
     };
 
@@ -330,7 +325,7 @@ void ImProcFunctions::rgbProc (Imagefloat* working, LabImage* lab, LUTf & hltone
 
 			if (sat!=0 || hCurveEnabled || sCurveEnabled || vCurveEnabled) {
 				float h,s,v;
-				rgb2hsv(r,g,b,h,s,v);
+				Color::rgb2hsv(r,g,b,h,s,v);
 				if (sat > 0.5) {
 					s = (1-(float)sat/100)*s+(float)sat/100*(1-SQR(SQR(1-min(s,1.0f))));
 					if (s<0) s=0;
@@ -373,7 +368,7 @@ void ImProcFunctions::rgbProc (Imagefloat* working, LabImage* lab, LUTf & hltone
 					}
 					
 				}
-				hsv2rgb(h,s,v,r,g,b);
+				Color::hsv2rgb(h,s,v,r,g,b);
 			}
 			 
             float x = toxyz[0][0] * r + toxyz[0][1] * g + toxyz[0][2] * b;
@@ -409,9 +404,9 @@ void ImProcFunctions::rgbProc (Imagefloat* working, LabImage* lab, LUTf & hltone
 			float fx = (0.002 * lab->a[i][j])/327.68 + fy;
 			float fz = fy - (0.005 * lab->b[i][j])/327.68;
 			
-			float x_ = 65535*Lab2xyz(fx)*D50x;
+			float x_ = 65535*Lab2xyz(fx)*Color::D50x;
 			float y_ = 65535*Lab2xyz(fy);
-			float z_ = 65535*Lab2xyz(fz)*D50z;
+			float z_ = 65535*Lab2xyz(fz)*Color::D50z;
 			
 			int R,G,B;
 			xyz2srgb(x_,y_,z_,R,G,B);
@@ -482,9 +477,9 @@ void ImProcFunctions::chrominanceCurve (LabImage* lold, LabImage* lnew, LUTf & a
             if (params->labCurve.avoidclip) {
 				//Luv limiter
 				float Y,u,v;
-				Lab2Yuv(lnew->L[i][j],atmp,btmp,Y,u,v);
+				Color::Lab2Yuv(lnew->L[i][j],atmp,btmp,Y,u,v);
 				//Yuv2Lab includes gamut restriction map
-				Yuv2Lab(Y,u,v,lnew->L[i][j],lnew->a[i][j],lnew->b[i][j], wp);
+				Color::Yuv2Lab(Y,u,v,lnew->L[i][j],lnew->a[i][j],lnew->b[i][j], wp);
 				
             } else {
 				//Luv limiter only
@@ -929,229 +924,7 @@ fclose(f);*/
 		else
 			return 0.0;
 	}
-	 
-	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-	
-	void ImProcFunctions::rgb2hsv (float r, float g, float b, float &h, float &s, float &v) {
-		
-		double var_R = r / 65535.0;
-		double var_G = g / 65535.0;
-		double var_B = b / 65535.0;
-		
-		double var_Min = min(var_R,var_G,var_B);
-		double var_Max = max(var_R,var_G,var_B);
-		double del_Max = var_Max - var_Min;
-		v = var_Max;
-		if (del_Max<0.00001 && del_Max>-0.00001) {  // no fabs, slow!
-			h = 0;
-			s = 0;
-		}
-		else {
-			s = del_Max/var_Max;
-			
-			if      ( var_R == var_Max ) h = (var_G - var_B)/del_Max; 
-			else if ( var_G == var_Max ) h = 2.0 + (var_B - var_R)/del_Max; 
-			else if ( var_B == var_Max ) h = 4.0 + (var_R - var_G)/del_Max; 
-			h /= 6.0;
-			
-			if ( h < 0 )  h += 1;
-			if ( h > 1 )  h -= 1;
-		}
-	}
 	
 	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-	//
-	// WARNING: rtengine/color.cc have similar functions. I guess that color.cc was mistakenly committed from the denoise branch
-	//          in changeset #2493dd00b9f3. Code cleanup about those methods will have to be made when the denoise branch will be merged
-	//
-	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-	
-	void ImProcFunctions::hsv2rgb (float h, float s, float v, float &r, float &g, float &b) {
-		
-		float h1 = h*6; // sector 0 to 5
-		int i = (int)h1;  // floor() is very slow, and h1 is always >0
-		float f = h1 - i; // fractional part of h
-		
-		float p = v * ( 1 - s );
-		float q = v * ( 1 - s * f );
-		float t = v * ( 1 - s * ( 1 - f ) );
-		
-		float r1,g1,b1;
 
-		if      (i==1)    {r1 = q;  g1 = v;  b1 = p;}
-		else if (i==2)    {r1 = p;  g1 = v;  b1 = t;}
-		else if (i==3)    {r1 = p;  g1 = q;  b1 = v;}
-		else if (i==4)    {r1 = t;  g1 = p;  b1 = v;}
-		else if (i==5)    {r1 = v;  g1 = p;  b1 = q;}
-		else /*i==(0|6)*/ {r1 = v;  g1 = t;  b1 = p;}
-		
-		r = ((r1)*65535.0);
-		g = ((g1)*65535.0);
-		b = ((b1)*65535.0);
-	}
-	
-	// Function copied for speed concerns
-	// Not exactly the same as above ; this one return a result in the [0.0 ; 1.0] range
-	void ImProcFunctions::hsv2rgb01 (float h, float s, float v, float &r, float &g, float &b) {
-
-		float h1 = h*6; // sector 0 to 5
-		int i = int(h1);
-		float f = h1 - i; // fractional part of h
-
-		float p = v * ( 1 - s );
-		float q = v * ( 1 - s * f );
-		float t = v * ( 1 - s * ( 1 - f ) );
-
-		if      (i==1)    {r = q;  g = v;  b = p;}
-		else if (i==2)    {r = p;  g = v;  b = t;}
-		else if (i==3)    {r = p;  g = q;  b = v;}
-		else if (i==4)    {r = t;  g = p;  b = v;}
-		else if (i==5)    {r = v;  g = p;  b = q;}
-		else /*(i==0|6)*/ {r = v;  g = t;  b = p;}
-	}
-
-	void ImProcFunctions::xyz2srgb (float x, float y, float z, float &r, float &g, float &b) {
-		
-		//Transform to output color.  Standard sRGB is D65, but internal representation is D50
-		//Note that it is only at this point that we should have need of clipping color data
-		
-		/*float x65 = d65_d50[0][0]*x + d65_d50[0][1]*y + d65_d50[0][2]*z ;
-		 float y65 = d65_d50[1][0]*x + d65_d50[1][1]*y + d65_d50[1][2]*z ;
-		 float z65 = d65_d50[2][0]*x + d65_d50[2][1]*y + d65_d50[2][2]*z ;
-		 
-		 r = sRGB_xyz[0][0]*x65 + sRGB_xyz[0][1]*y65 + sRGB_xyz[0][2]*z65;
-		 g = sRGB_xyz[1][0]*x65 + sRGB_xyz[1][1]*y65 + sRGB_xyz[1][2]*z65;
-		 b = sRGB_xyz[2][0]*x65 + sRGB_xyz[2][1]*y65 + sRGB_xyz[2][2]*z65;*/
-		
-		/*r = sRGBd65_xyz[0][0]*x + sRGBd65_xyz[0][1]*y + sRGBd65_xyz[0][2]*z ;
-		 g = sRGBd65_xyz[1][0]*x + sRGBd65_xyz[1][1]*y + sRGBd65_xyz[1][2]*z ;
-		 b = sRGBd65_xyz[2][0]*x + sRGBd65_xyz[2][1]*y + sRGBd65_xyz[2][2]*z ;*/
-		
-		r = ((sRGB_xyz[0][0]*x + sRGB_xyz[0][1]*y + sRGB_xyz[0][2]*z)) ;
-		g = ((sRGB_xyz[1][0]*x + sRGB_xyz[1][1]*y + sRGB_xyz[1][2]*z)) ;
-		b = ((sRGB_xyz[2][0]*x + sRGB_xyz[2][1]*y + sRGB_xyz[2][2]*z)) ;
-		
-	}
-	
-	
-	void ImProcFunctions::xyz2rgb (float x, float y, float z, float &r, float &g, float &b, double rgb_xyz[3][3]) {
-		
-		//Transform to output color.  Standard sRGB is D65, but internal representation is D50
-		//Note that it is only at this point that we should have need of clipping color data
-		
-		/*float x65 = d65_d50[0][0]*x + d65_d50[0][1]*y + d65_d50[0][2]*z ;
-		 float y65 = d65_d50[1][0]*x + d65_d50[1][1]*y + d65_d50[1][2]*z ;
-		 float z65 = d65_d50[2][0]*x + d65_d50[2][1]*y + d65_d50[2][2]*z ;
-		 
-		 r = sRGB_xyz[0][0]*x65 + sRGB_xyz[0][1]*y65 + sRGB_xyz[0][2]*z65;
-		 g = sRGB_xyz[1][0]*x65 + sRGB_xyz[1][1]*y65 + sRGB_xyz[1][2]*z65;
-		 b = sRGB_xyz[2][0]*x65 + sRGB_xyz[2][1]*y65 + sRGB_xyz[2][2]*z65;*/
-		
-		/*r = sRGBd65_xyz[0][0]*x + sRGBd65_xyz[0][1]*y + sRGBd65_xyz[0][2]*z ;
-		 g = sRGBd65_xyz[1][0]*x + sRGBd65_xyz[1][1]*y + sRGBd65_xyz[1][2]*z ;
-		 b = sRGBd65_xyz[2][0]*x + sRGBd65_xyz[2][1]*y + sRGBd65_xyz[2][2]*z ;*/
-		
-		
-		
-		r = ((rgb_xyz[0][0]*x + rgb_xyz[0][1]*y + rgb_xyz[0][2]*z)) ;
-		g = ((rgb_xyz[1][0]*x + rgb_xyz[1][1]*y + rgb_xyz[1][2]*z)) ;
-		b = ((rgb_xyz[2][0]*x + rgb_xyz[2][1]*y + rgb_xyz[2][2]*z)) ;
-		
-	}
- 	
-void ImProcFunctions::calcGamma (double pwr, double ts, int mode, int imax, double &gamma0, double &gamma1, double &gamma2, double &gamma3, double &gamma4, double &gamma5) {
-{//from Dcraw (D.Coffin)
-  int i;
-  double g[6], bnd[2]={0,0};
-
-  g[0] = pwr;
-  g[1] = ts;
-  g[2] = g[3] = g[4] = 0;
-  bnd[g[1] >= 1] = 1;
-  if (g[1] && (g[1]-1)*(g[0]-1) <= 0) {
-    for (i=0; i < 48; i++) {
-      g[2] = (bnd[0] + bnd[1])/2;
-      if (g[0]) bnd[(pow(g[2]/g[1],-g[0]) - 1)/g[0] - 1/g[2] > -1] = g[2];
-      else	bnd[g[2]/exp(1-1/g[2]) < g[1]] = g[2];
-    }
-    g[3] = g[2] / g[1];
-    if (g[0]) g[4] = g[2] * (1/g[0] - 1);
-  }
-  if (g[0]) g[5] = 1 / (g[1]*SQR(g[3])/2 - g[4]*(1 - g[3]) +
-		(1 - pow(g[3],1+g[0]))*(1 + g[4])/(1 + g[0])) - 1;
-  else      g[5] = 1 / (g[1]*SQR(g[3])/2 + 1
-		- g[2] - g[3] -	g[2]*g[3]*(log(g[3]) - 1)) - 1;
-  if (!mode--) {
-   gamma0=g[0];gamma1=g[1];gamma2=g[2];gamma3=g[3];gamma4=g[4];gamma5=g[5];
-    return;
-  }
-}
-}
-	
-	void ImProcFunctions::Lab2XYZ(float L, float a, float b, float &x, float &y, float &z) {
-		float fy = (0.00862069 * L) + 0.137932; // (L+16)/116
-		float fx = (0.002 * a) + fy;
-		float fz = fy - (0.005 * b);
-		
-		x = 65535*f2xyz(fx)*D50x;
-		y = 65535*f2xyz(fy);
-		z = 65535*f2xyz(fz)*D50z;
-	}
-	
-	void ImProcFunctions::XYZ2Lab(float X, float Y, float Z, float &L, float &a, float &b) {
-		
-		float fx = (X<65535.0 ? cachef[X] : (327.68*exp(log(X/MAXVAL)/3.0 )));
-		float fy = (Y<65535.0 ? cachef[Y] : (327.68*exp(log(Y/MAXVAL)/3.0 )));
-		float fz = (Z<65535.0 ? cachef[Z] : (327.68*exp(log(Z/MAXVAL)/3.0 )));
-		
-		L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68;
-		a = (500.0 * (fx - fy) );
-		b = (200.0 * (fy - fz) );
-		
-	}
-	
-	void ImProcFunctions::Lab2Yuv(float L, float a, float b, float &Y, float &u, float &v) {
-		float fy = (0.00862069 * L/327.68) + 0.137932; // (L+16)/116
-		float fx = (0.002 * a/327.68) + fy;
-		float fz = fy - (0.005 * b/327.68);
-		
-		float X = 65535.0*f2xyz(fx)*D50x;
-		Y = 65535.0*f2xyz(fy);
-		float Z = 65535.0*f2xyz(fz)*D50z;
-		
-		u = 4.0*X/(X+15*Y+3*Z)-u0;
-		v = 9.0*Y/(X+15*Y+3*Z)-v0;
-		
-		/*float u0 = 4*D50x/(D50x+15+3*D50z);
-		 float v0 = 9/(D50x+15+3*D50z);
-		 u -= u0;
-		 v -= v0;*/
-		
-	}
-	
-	void ImProcFunctions::Yuv2Lab(float Yin, float u, float v, float &L, float &a, float &b, double wp[3][3]) {
-		
-		float u1 = u + u0;
-		float v1 = v + v0;
-		
-		float Y = Yin;
-		float X = (9*u1*Y)/(4*v1*D50x); 
-		float Z = (12 - 3*u1 - 20*v1)*Y/(4*v1*D50z);
-		
-		gamutmap(X,Y,Z,wp);
-		
-		float fx = (X<65535.0 ? cachef[X] : (327.68*exp(log(X/MAXVAL)/3.0 )));
-		float fy = (Y<65535.0 ? cachef[Y] : (327.68*exp(log(Y/MAXVAL)/3.0 )));
-		float fz = (Z<65535.0 ? cachef[Z] : (327.68*exp(log(Z/MAXVAL)/3.0 )));
-		
-		L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68;
-		a = (500.0 * (fx - fy) );
-		b = (200.0 * (fy - fz) );
-		
-	}
-	
-#include "gamutbdy.cc"
-}
-
-#undef eps_max
-#undef kappa
+}
\ No newline at end of file
diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h
index 205f97243..4b22a3038 100644
--- a/rtengine/improcfun.h
+++ b/rtengine/improcfun.h
@@ -25,6 +25,7 @@
 #include "procparams.h"
 #include "shmap.h"
 #include "coord2d.h"
+#include "color.h"
 #include "labimage.h"
 #include "LUT.h"
 #include "lcp.h"
@@ -166,27 +167,6 @@ class ImProcFunctions {
 		void getAutoExp       (LUTu & histogram, int histcompr, double defgain, double clip, double& expcomp, int& bright, int& contr, int& black, int& hlcompr, int& hlcomprthresh);
 		static double getAutoDistor  (const Glib::ustring& fname, int thumb_size);	
 		double getTransformAutoFill (int oW, int oH, const LCPMapper *pLCPMap=NULL);
-
-		static void rgb2hsv (float r, float g, float b, float &h, float &s, float &v);
-		static void hsv2rgb (float h, float s, float v, float &r, float &g, float &b);
-		static void hsv2rgb01 (float h, float s, float v, float &r, float &g, float &b);
-		void xyz2srgb (float x, float y, float z, float &r, float &g, float &b);
-		void xyz2rgb (float x, float y, float z, float &r, float &g, float &b, double rgb_xyz[3][3]);
-		void Lab2XYZ(float L, float a, float b, float &x, float &y, float &z);
-		void XYZ2Lab(float X, float Y, float Z, float &L, float &a, float &b);
-		void Lab2Yuv(float L, float a, float b, float &Y, float &u, float &v);
-		void Yuv2Lab(float Y, float u, float v, float &L, float &a, float &b, double wp[3][3]);
-		void calcGamma (double pwr, double ts, int mode, int imax, double &gamma0, double &gamma1, double &gamma2, double &gamma3, double &gamma4,double &gamma5);
-
-		//void gamutmap(LabImage* );
-		void gamutmap(float &X, float &Y, float &Z, const double p[3][3]);
-
-		static inline float f2xyz(float f) {
-			const float epsilonExpInv3 = 6.0/29.0;
-			const float kappaInv = 27.0/24389.0;  // inverse of kappa
-
-			return (f > epsilonExpInv3) ? f*f*f : (116 * f - 16) * kappaInv;
-		}
 };
 }
 #endif
diff --git a/rtengine/iplab2rgb.cc b/rtengine/iplab2rgb.cc
index fbdad662b..c77d1b6a0 100644
--- a/rtengine/iplab2rgb.cc
+++ b/rtengine/iplab2rgb.cc
@@ -74,9 +74,9 @@ void ImProcFunctions::lab2monitorRgb (LabImage* lab, Image8* image) {
 				fx = (0.002 * ra[j]) / 327.68 + fy;
 				fz = fy - (0.005 * rb[j]) / 327.68;
 				
-				x_ = f2xyz(fx)*D50x;
-				y_ = f2xyz(fy);
-				z_ = f2xyz(fz)*D50z;
+				x_ = Color::f2xyz(fx)*Color::D50x;
+				y_ = Color::f2xyz(fy);
+				z_ = Color::f2xyz(fz)*Color::D50z;
 
                 buffer[iy++] = (unsigned short)CLIP(x_* MAXVAL+0.5);
                 buffer[iy++] = (unsigned short)CLIP(y_* MAXVAL+0.5);
@@ -108,11 +108,11 @@ void ImProcFunctions::lab2monitorRgb (LabImage* lab, Image8* image) {
 				fx = (0.002 * ra[j]) / 327.68 + fy;
 				fz = fy - (0.005 * rb[j]) / 327.68;
 				
-				x_ = 65535.0 * f2xyz(fx)*D50x;
-				y_ = 65535.0 * f2xyz(fy);
-				z_ = 65535.0 * f2xyz(fz)*D50z;
+				x_ = 65535.0 * Color::f2xyz(fx)*Color::D50x;
+				y_ = 65535.0 * Color::f2xyz(fy);
+				z_ = 65535.0 * Color::f2xyz(fz)*Color::D50z;
 
-				xyz2srgb(x_,y_,z_,R,G,B);
+				Color::xyz2srgb(x_,y_,z_,R,G,B);
 				
 				/* copy RGB */
 				//int R1=((int)gamma2curve[(R)]) 
@@ -165,9 +165,9 @@ Image8* ImProcFunctions::lab2rgb (LabImage* lab, int cx, int cy, int cw, int ch,
 				float fx = (0.002 * ra[j])/327.68 + fy;
 				float fz = fy - (0.005 * rb[j])/327.68;
 				
-				float x_ = 65535.0 * f2xyz(fx)*D50x;
-				float y_ = 65535.0 * f2xyz(fy);
-				float z_ = 65535.0 * f2xyz(fz)*D50z;
+				float x_ = 65535.0 * Color::f2xyz(fx)*Color::D50x;
+				float y_ = 65535.0 * Color::f2xyz(fy);
+				float z_ = 65535.0 * Color::f2xyz(fz)*Color::D50z;
 
                 buffer[iy++] = CLIP((int)(x_+0.5));
                 buffer[iy++] = CLIP((int)(y_+0.5));
@@ -205,11 +205,11 @@ Image8* ImProcFunctions::lab2rgb (LabImage* lab, int cx, int cy, int cw, int ch,
 				float fx = (0.002 * ra[j])/327.68 + fy;
 				float fz = fy - (0.005 * rb[j])/327.68;
 				
-				float x_ = 65535.0 * f2xyz(fx)*D50x;
-				float y_ = 65535.0 * f2xyz(fy);
-				float z_ = 65535.0 * f2xyz(fz)*D50z;
+				float x_ = 65535.0 * Color::f2xyz(fx)*Color::D50x;
+				float y_ = 65535.0 * Color::f2xyz(fy);
+				float z_ = 65535.0 * Color::f2xyz(fz)*Color::D50z;
 
-				xyz2rgb(x_,y_,z_,R,G,B,rgb_xyz);
+				Color::xyz2rgb(x_,y_,z_,R,G,B,rgb_xyz);
 
                 image->data[ix++] = (int)gamma2curve[CLIP(R)] >> 8;
                 image->data[ix++] = (int)gamma2curve[CLIP(G)] >> 8;
@@ -249,9 +249,9 @@ Image16* ImProcFunctions::lab2rgb16 (LabImage* lab, int cx, int cy, int cw, int
 				float fx = (0.002 * ra[j])/327.68 + fy;
 				float fz = fy - (0.005 * rb[j])/327.68;
 				
-				float x_ = 65535.0 * f2xyz(fx)*D50x;
-				float y_ = 65535.0 * f2xyz(fy);
-				float z_ = 65535.0 * f2xyz(fz)*D50z;
+				float x_ = 65535.0 * Color::f2xyz(fx)*Color::D50x;
+				float y_ = 65535.0 * Color::f2xyz(fy);
+				float z_ = 65535.0 * Color::f2xyz(fz)*Color::D50z;
 
 				xa[j-cx] = CLIP((int)(x_+0.5));
 				ya[j-cx] = CLIP((int)(y_+0.5));
@@ -280,11 +280,11 @@ Image16* ImProcFunctions::lab2rgb16 (LabImage* lab, int cx, int cy, int cw, int
 				float fx = (0.002 * ra[j])/327.68 + fy;
 				float fz = fy - (0.005 * rb[j])/327.68;
 				
-				float x_ = 65535.0 * f2xyz(fx)*D50x;
-				float y_ = 65535.0 * f2xyz(fy);
-				float z_ = 65535.0 * f2xyz(fz)*D50z;
+				float x_ = 65535.0 * Color::f2xyz(fx)*Color::D50x;
+				float y_ = 65535.0 * Color::f2xyz(fy);
+				float z_ = 65535.0 * Color::f2xyz(fz)*Color::D50z;
 
-				xyz2srgb(x_,y_,z_,R,G,B);
+				Color::xyz2srgb(x_,y_,z_,R,G,B);
 
 				image->r[i-cy][j-cx] = (int)gamma2curve[CLIP(R)];
 				image->g[i-cy][j-cx] = (int)gamma2curve[CLIP(G)];
@@ -342,7 +342,7 @@ Image16* ImProcFunctions::lab2rgb16b (LabImage* lab, int cx, int cy, int cw, int
 	else //free gamma selected
 	{
 	if(slpos==0) slpos=eps;
-	calcGamma(pwr, ts, mode, imax,g_a0,g_a1,g_a2,g_a3,g_a4,g_a5);// call to calcGamma with selected gamma and slope : return parameters for LCMS2
+	Color::calcGamma(pwr, ts, mode, imax,g_a0,g_a1,g_a2,g_a3,g_a4,g_a5);// call to calcGamma with selected gamma and slope : return parameters for LCMS2
 	ga4=g_a3*ts;
 	//printf("g_a0=%f g_a1=%f g_a2=%f g_a3=%f g_a4=%f\n", g_a0,g_a1,g_a2,g_a3,g_a4);
 	ga0=gampos;ga1=1./(1.0+g_a4);ga2=g_a4/(1.0 + g_a4);ga3=1./slpos;ga5=0.0;
@@ -389,9 +389,9 @@ Image16* ImProcFunctions::lab2rgb16b (LabImage* lab, int cx, int cy, int cw, int
 				float fx = (0.002 * ra[j])/327.68 + fy;
 				float fz = fy - (0.005 * rb[j])/327.68;
 				
-				float x_ = 65535.0 * f2xyz(fx)*D50x;
-				float y_ = 65535.0 * f2xyz(fy);
-				float z_ = 65535.0 * f2xyz(fz)*D50z;
+				float x_ = 65535.0 * Color::f2xyz(fx)*Color::D50x;
+				float y_ = 65535.0 * Color::f2xyz(fy);
+				float z_ = 65535.0 * Color::f2xyz(fz)*Color::D50z;
 
 				xa[j-cx] = CLIP((int)x_);
 				ya[j-cx] = CLIP((int)y_);
@@ -420,11 +420,11 @@ Image16* ImProcFunctions::lab2rgb16b (LabImage* lab, int cx, int cy, int cw, int
 				float fx = (0.002 * ra[j])/327.68 + fy;
 				float fz = fy - (0.005 * rb[j])/327.68;
 				
-				float x_ = 65535.0 * f2xyz(fx)*D50x;
-				float y_ = 65535.0 * f2xyz(fy);
-				float z_ = 65535.0 * f2xyz(fz)*D50z;
+				float x_ = 65535.0 * Color::f2xyz(fx)*Color::D50x;
+				float y_ = 65535.0 * Color::f2xyz(fy);
+				float z_ = 65535.0 * Color::f2xyz(fz)*Color::D50z;
 
-				xyz2srgb(x_,y_,z_,R,G,B);
+				Color::xyz2srgb(x_,y_,z_,R,G,B);
 
 				image->r[i-cy][j-cx] = (int)gamma2curve[CLIP(R)];
 				image->g[i-cy][j-cx] = (int)gamma2curve[CLIP(G)];
diff --git a/rtengine/ipvibrance.cc b/rtengine/ipvibrance.cc
index c3e773884..af06a172a 100644
--- a/rtengine/ipvibrance.cc
+++ b/rtengine/ipvibrance.cc
@@ -2145,10 +2145,10 @@ void ImProcFunctions::vibrance (LabImage* lab) {
 				fx = (0.002 * aprov1) + fy;
 				fz = fy - (0.005 * bprov1);
 
-				x_ = 65535.0 * f2xyz(fx)*D50x;
-				y_ = 65535.0 * f2xyz(fy);
-				z_ = 65535.0 * f2xyz(fz)*D50z;
-				xyz2rgb(x_,y_,z_,R,G,B,wip);
+				x_ = 65535.0 * Color::f2xyz(fx)*Color::D50x;
+				y_ = 65535.0 * Color::f2xyz(fy);
+				z_ = 65535.0 * Color::f2xyz(fz)*Color::D50z;
+				Color::xyz2rgb(x_,y_,z_,R,G,B,wip);
 
 				// gamut control before saturation to put Lab values in future gamut, but not RGB
 				if (allwaysingamut) {        // gamut control
@@ -2323,10 +2323,10 @@ void ImProcFunctions::vibrance (LabImage* lab) {
 					fxx = (0.002 * aprovn) + fyy;
 					fzz = fyy - (0.005 * bprovn);
 
-					xx_ = 65535.0 * f2xyz(fxx)*D50x;
-					yy_ = 65535.0 * f2xyz(fyy);
-					zz_ = 65535.0 * f2xyz(fzz)*D50z;
-					xyz2rgb(xx_,yy_,zz_,RR,GG,BB,wip);
+					xx_ = 65535.0 * Color::f2xyz(fxx)*Color::D50x;
+					yy_ = 65535.0 * Color::f2xyz(fyy);
+					zz_ = 65535.0 * Color::f2xyz(fzz)*Color::D50z;
+					Color::xyz2rgb(xx_,yy_,zz_,RR,GG,BB,wip);
 
 					if(RR<0.0 || GG < 0.0 || BB < 0.0) {
 #ifdef _DEBUG
diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc
index 16908fce4..2bf944b0e 100644
--- a/rtengine/rawimagesource.cc
+++ b/rtengine/rawimagesource.cc
@@ -2214,8 +2214,8 @@ void RawImageSource::HLRecovery_CIELab (float* rin, float* gin, float* bin, floa
             double fz = fy - y + z;
             double fx = fy + x - y;
 
-			double zr = ImProcFunctions::f2xyz(fz);
-            double xr = ImProcFunctions::f2xyz(fx);
+			double zr = Color::f2xyz(fz);
+            double xr = Color::f2xyz(fx);
 
             x = xr*65535.0 ;
             y = yy;
diff --git a/rtengine/rt_math.h b/rtengine/rt_math.h
index 1a1c339d9..c1049e279 100644
--- a/rtengine/rt_math.h
+++ b/rtengine/rt_math.h
@@ -5,9 +5,6 @@
 
 #include "rtengine.h"
 
-#define D50x 0.96422
-#define D50z 0.82521
-
 namespace rtengine {
 	static const int MAXVAL = 0xffff;
 
diff --git a/rtgui/hsvequalizer.cc b/rtgui/hsvequalizer.cc
index bc247c449..b5e77f1b1 100644
--- a/rtgui/hsvequalizer.cc
+++ b/rtgui/hsvequalizer.cc
@@ -18,7 +18,7 @@
  */
 
 #include "hsvequalizer.h"
-#include "../rtengine/improcfun.h"
+#include "../rtengine/color.h"
 
 using namespace rtengine;
 using namespace rtengine::procparams;
@@ -195,19 +195,19 @@ void HSVEqualizer::colorForValue (double valX, double valY) {
 			h -= 1.0f;
 		else if (h < 0.0f)
 			h += 1.0f;
-		ImProcFunctions::hsv2rgb01(h, 0.5f, 0.5f, r, g, b);
+		Color::hsv2rgb01(h, 0.5f, 0.5f, r, g, b);
 		red = double(r);
 		green = double(g);
 		blue = double(b);
 	}
 	else if (ce == sshape) {   // Saturation = f(Hue)
-		ImProcFunctions::hsv2rgb01(float(valX), float(valY), 0.5f, r, g, b);
+		Color::hsv2rgb01(float(valX), float(valY), 0.5f, r, g, b);
 		red = double(r);
 		green = double(g);
 		blue = double(b);
 	}
 	else if (ce == vshape) {   // Value = f(Hue)
-		ImProcFunctions::hsv2rgb01(float(valX), 0.5f, float(valY), r, g, b);
+		Color::hsv2rgb01(float(valX), 0.5f, float(valY), r, g, b);
 		red = double(r);
 		green = double(g);
 		blue = double(b);
diff --git a/rtgui/navigator.cc b/rtgui/navigator.cc
index 8252f959c..591fa2302 100644
--- a/rtgui/navigator.cc
+++ b/rtgui/navigator.cc
@@ -21,8 +21,11 @@
 #include "../rtengine/iccmatrices.h"
 #include "../rtengine/iccstore.h"
 #include "../rtengine/curves.h"
+#include "../rtengine/color.h"
 #include "../rtengine/rt_math.h"
 
+using namespace rtengine;
+
 Navigator::Navigator () {
 
 	set_label (M("MAIN_MSG_NAVIGATOR"));
@@ -88,11 +91,11 @@ void Navigator::pointerMoved (bool validPos, Glib::ustring profile, int x, int y
 		R->set_text (Glib::ustring::compose (M("NAVIGATOR_R_VALUE"), r));
 		G->set_text (Glib::ustring::compose (M("NAVIGATOR_G_VALUE"), g));
 		B->set_text (Glib::ustring::compose (M("NAVIGATOR_B_VALUE"), b));
-		int h, s, v;
-		rgb2hsv (r, g, b, h, s, v);
-		H->set_text (Glib::ustring::compose (M("NAVIGATOR_H_VALUE"), h));
-		S->set_text (Glib::ustring::compose (M("NAVIGATOR_S_VALUE"), s));
-		V->set_text (Glib::ustring::compose (M("NAVIGATOR_V_VALUE"), v));
+		float h, s, v;
+		Color::rgb2hsv (r*0xffff/0xff, g*0xffff/0xff, b*0xffff/0xff, h, s, v);
+		H->set_text (Glib::ustring::compose (M("NAVIGATOR_H_VALUE"), (int)(h*(float)0xff)));
+		S->set_text (Glib::ustring::compose (M("NAVIGATOR_S_VALUE"), (int)(s*(float)0xff)));
+		V->set_text (Glib::ustring::compose (M("NAVIGATOR_V_VALUE"), (int)(v*(float)0xff)));
 		int LAB_a, LAB_b, LAB_l;
 		//rgb2lab (r, g, b, LAB_l, LAB_a, LAB_b);
 		rgb2lab (profile, r, g, b, LAB_l, LAB_a, LAB_b);
@@ -104,67 +107,12 @@ void Navigator::pointerMoved (bool validPos, Glib::ustring profile, int x, int y
 	}
 }
 
-void Navigator::rgb2hsv (int r, int g, int b, int &h, int &s, int &v) {
-
-	volatile double var_R = r / 255.0;
-	volatile double var_G = g / 255.0;
-	volatile double var_B = b / 255.0;
-
-	/*volatile double var_Min = min(var_R,var_G,var_B);
-	volatile double var_Max = max(var_R,var_G,var_B);
-	double del_Max = var_Max - var_Min;
-    double	V = (var_Max + var_Min) / 2;
-    double H, S;
-	if (fabs(del_Max)<0.0000001) {
-		H = 0;
-		S = 0;
-	}
-	else {
-		if (V < 0.5) S = del_Max / (var_Max + var_Min);
-		else         S = del_Max / (2 - var_Max - var_Min);
-
-		double del_R = ( ( ( var_Max - var_R ) / 6.0 ) + ( del_Max / 2.0 ) ) / del_Max;
-		double del_G = ( ( ( var_Max - var_G ) / 6.0 ) + ( del_Max / 2.0 ) ) / del_Max;
-		double del_B = ( ( ( var_Max - var_B ) / 6.0 ) + ( del_Max / 2.0 ) ) / del_Max;
-		if      ( var_R == var_Max ) H = del_B - del_G; 
-		else if ( var_G == var_Max ) H = (1.0 / 3.0) + del_R - del_B; 
-		else if ( var_B == var_Max ) H = (2.0 / 3.0) + del_G - del_R; 
-
-		if ( H < 0 )  H += 1;
-		if ( H > 1 )  H -= 1;
-	}*/
-	
-	double var_Min = rtengine::min(var_R,var_G,var_B);
-	double var_Max = rtengine::max(var_R,var_G,var_B);
-	double del_Max = var_Max - var_Min;
-	double V = var_Max;
-	double H, S;
-	if (fabs(del_Max)<0.001) {
-		H = 0;
-		S = 0;
-	}
-	else {
-		S = del_Max/var_Max;
-		
-		if      ( var_R == var_Max ) H = (var_G - var_B)/del_Max; 
-		else if ( var_G == var_Max ) H = 2.0 + (var_B - var_R)/del_Max; 
-		else if ( var_B == var_Max ) H = 4.0 + (var_R - var_G)/del_Max; 
-		H /= 6.0;
-		
-		if ( H < 0 )  H += 1;
-		if ( H > 1 )  H -= 1;
-	}
-    
-    h = (int)(H*255.0);
-    s = (int)(S*255.0);
-    v = (int)(V*255.0);
-}
 
 void Navigator::rgb2lab (Glib::ustring profile, int r, int g, int b, int &LAB_l, int &LAB_a, int &LAB_b) {
 	
 	double xyz_rgb[3][3];
-	double ep=216.0/24389.0;
-	double ka=24389.0/27.0;
+	const double ep=216.0/24389.0;
+	const double ka=24389.0/27.0;
 
 	double var_R = r / 255.0;
 	double var_G = g / 255.0;
@@ -177,15 +125,15 @@ void Navigator::rgb2lab (Glib::ustring profile, int r, int g, int b, int &LAB_l,
 // today as the gamma output can not be configured
 // it is better that the user has the gamma of the output space
 	if ( var_R > 0.04045 ) 
-		var_R = pow ( ( ( var_R + 0.055 ) / 1.055 ), rtengine::CurveFactory::sRGBGammaCurve);
+		var_R = pow ( ( ( var_R + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
 	else    
 		var_R = var_R / 12.92;
 	if ( var_G > 0.04045 ) 
-		var_G = pow ( ( ( var_G + 0.055 ) / 1.055 ), rtengine::CurveFactory::sRGBGammaCurve);
+		var_G = pow ( ( ( var_G + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
 	else                   
 		var_G = var_G / 12.92;
 	if ( var_B > 0.04045 ) 
-		var_B = pow ( ( ( var_B + 0.055 ) / 1.055 ), rtengine::CurveFactory::sRGBGammaCurve);
+		var_B = pow ( ( ( var_B + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
 	else    
 		var_B = var_B / 12.92;
 	} 
@@ -220,9 +168,9 @@ void Navigator::rgb2lab (Glib::ustring profile, int r, int g, int b, int &LAB_l,
 		}
 
 	double varxx,varyy,varzz;
-	double var_X = ( xyz_rgb[0][0]*var_R + xyz_rgb[0][1]*var_G + xyz_rgb[0][2]*var_B ) / D50x;
+	double var_X = ( xyz_rgb[0][0]*var_R + xyz_rgb[0][1]*var_G + xyz_rgb[0][2]*var_B ) / Color::D50x;
 	double var_Y = ( xyz_rgb[1][0]*var_R + xyz_rgb[1][1]*var_G + xyz_rgb[1][2]*var_B ) ;
-	double var_Z = ( xyz_rgb[2][0]*var_R + xyz_rgb[2][1]*var_G + xyz_rgb[2][2]*var_B ) / D50z;
+	double var_Z = ( xyz_rgb[2][0]*var_R + xyz_rgb[2][1]*var_G + xyz_rgb[2][2]*var_B ) / Color::D50z;
 
 	varxx = var_X>ep?pow (var_X, 1.0/3.0):( ka * var_X  +  16.0) / 116.0 ;
 	varyy = var_Y>ep?pow (var_Y, 1.0/3.0):( ka * var_Y  +  16.0) / 116.0 ;
diff --git a/rtgui/navigator.h b/rtgui/navigator.h
index 3d8a0c35c..207a9daf3 100644
--- a/rtgui/navigator.h
+++ b/rtgui/navigator.h
@@ -34,8 +34,6 @@ class Navigator : public Gtk::Frame, public PointerMotionListener {
 		Gtk::Label *H, *S, *V;
 		Gtk::Label *LAB_A, *LAB_B, *LAB_L;
 
-		void rgb2hsv 	(int r, int g, int b, int &h, int &s, int &v);
-		//void rgb2lab (int r, int g, int b, int &LAB_l, int &LAB_a, int &LAB_b);
 		void rgb2lab (Glib::ustring profile, int r, int g, int b, int &LAB_l, int &LAB_a, int &LAB_b);
 		
 		void setInvalid (int fullWidth=-1, int fullHeight=-1);