diff --git a/rtengine/color.cc b/rtengine/color.cc
index 68fd467e2..3be412b6d 100644
--- a/rtengine/color.cc
+++ b/rtengine/color.cc
@@ -34,6 +34,7 @@ extern const Settings* settings;
 
 cmsToneCurve* Color::linearGammaTRC;
 LUTf Color::cachef;
+LUTf Color::cachefy;
 LUTf Color::gamma2curve;
 
 LUTf Color::gammatab;
@@ -55,21 +56,6 @@ LUTf Color::igammatab_115_2;
 LUTf Color::gammatab_145_3;
 LUTf Color::igammatab_145_3;
 
-// 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; //(MAXVALF* 216.0f/24389.0);
-const double Color::eps = 216.0f / 24389.0; //0.008856
-const double Color::kappa = 24389.0 / 27.0; //903.29630;
-
-const float Color::D50x = 0.9642f; //0.96422;
-const float Color::D50z = 0.8249f; //0.82521;
-const double Color::u0 = 4.0 * D50x / (D50x + 15 + 3 * D50z);
-const double Color::v0 = 9.0 / (D50x + 15 + 3 * D50z);
-const double Color::epskap = 8.0;
 /*
  * Munsell Lch correction
  * Copyright (c) 2011  Jacques Desmis <jdesmis@gmail.com>
@@ -137,6 +123,7 @@ void Color::init ()
     constexpr auto maxindex = 65536;
 
     cachef(maxindex, LUT_CLIP_BELOW);
+    cachefy(maxindex, LUT_CLIP_BELOW);
     gammatab(maxindex, 0);
     gammatabThumb(maxindex, 0);
 
@@ -180,6 +167,23 @@ void Color::init ()
         }
 #ifdef _OPENMP
         #pragma omp section
+#endif
+        {
+            int i = 0;
+            int epsmaxint = eps_max;
+
+            for (; i <= epsmaxint; i++)
+            {
+                cachefy[i] = 327.68 * (kappa * i / MAXVALF);
+            }
+
+            for(; i < maxindex; i++)
+            {
+                cachefy[i] = 327.68 * (116.0 * std::cbrt((double)i / MAXVALF) - 16.0);
+            }
+        }
+#ifdef _OPENMP
+        #pragma omp section
 #endif
         {
             for (int i = 0; i < maxindex; i++)
@@ -360,150 +364,141 @@ void Color::cleanup ()
     }
 }
 
-void Color::rgb2lab (Glib::ustring profile, Glib::ustring profileW, int r, int g, int b, float &LAB_l, float &LAB_a, float &LAB_b, bool workingSpace)
-{
-    double xyz_rgb[3][3];
-    const double ep = 216.0 / 24389.0;
-    const double ka = 24389.0 / 27.0;
+void Color::rgb2lab01 (const Glib::ustring &profile, const Glib::ustring &profileW, float r, float g, float b, float &LAB_l, float &LAB_a, float &LAB_b, bool workingSpace)
+{ // do not use this function in a loop. It really eats processing time caused by Glib::ustring comparisons
 
-    double var_R = r / 65535.0;
-    double var_G = g / 65535.0;
-    double var_B = b / 65535.0;
+    Glib::ustring profileCalc = "sRGB"; //default
 
-    Glib::ustring profileCalc;
-    profileCalc = "sRGB"; //default
-
-    if (workingSpace) {
-        profileCalc = profileW;    //display working
-    }
-
-    else {// if you want display = output space
-        if (profile == "RT_sRGB" || profile == "RT_sRGB_gBT709" || profile == "RT_sRGB_g10") {
-            profileCalc = "sRGB";
-        }
-
-        if (profile == "ProPhoto" || profile == "RT_Large_gBT709" || profile == "RT_Large_g10"  || profile == "RT_Large_gsRGB") {
-            profileCalc = "ProPhoto";
-        }
-
-        if (profile == "AdobeRGB1998" || profile == "RT_Medium_gsRGB") {
-            profileCalc = "Adobe RGB";
-        }
-
-        if (profile == "WideGamutRGB") {
-            profileCalc = "WideGamut";
-        }
-    }
-
-    if (workingSpace) {//display working
+    if (workingSpace) {//display working profile
+        profileCalc = profileW;
         if (profileW == "sRGB") { //apply sRGB inverse gamma
 
-            if ( var_R > 0.04045 ) {
-                var_R = pow ( ( ( var_R + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
+            if (r > 0.04045f) {
+                r = pow_F(((r + 0.055f) / 1.055f), rtengine::Color::sRGBGammaCurve);
             } else {
-                var_R = var_R / 12.92;
+                r /= 12.92f;
             }
 
-            if ( var_G > 0.04045 ) {
-                var_G = pow ( ( ( var_G + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
+            if (g > 0.04045f) {
+                g = pow_F(((g + 0.055f) / 1.055f), rtengine::Color::sRGBGammaCurve);
             } else {
-                var_G = var_G / 12.92;
+                g /= 12.92f;
             }
 
-            if ( var_B > 0.04045 ) {
-                var_B = pow ( ( ( var_B + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
+            if (b > 0.04045f) {
+                b = pow_F(((b + 0.055f) / 1.055f), rtengine::Color::sRGBGammaCurve);
             } else {
-                var_B = var_B / 12.92;
+                b /= 12.92f;
             }
         } else if (profileW == "ProPhoto") { // apply inverse gamma 1.8
-            var_R = pow ( var_R, 1.8);
-            var_G = pow ( var_G, 1.8);
-            var_B = pow ( var_B, 1.8);
+            r = pow_F(r, 1.8f);
+            g = pow_F(g, 1.8f);
+            b = pow_F(b, 1.8f);
+        } else if (profile == "Rec2020") {
+            if (r > 0.0795f) {
+                r = pow_F(((r + 0.0954f) / 1.0954f), 2.2f);
+            } else {
+                r /= 4.5f;
+            }
+
+            if (g > 0.0795f) {
+                g = pow_F(((g + 0.0954f) / 1.0954f), 2.2f);
+            } else {
+                g /= 4.5f;
+            }
+
+            if (b > 0.0795f) {
+                b = pow_F(((b + 0.0954f) / 1.0954f), 2.2f);
+            } else {
+                b /= 4.5f;
+            }
         } else { // apply inverse gamma 2.2
-            var_R = pow ( var_R, 2.2);
-            var_G = pow ( var_G, 2.2);
-            var_B = pow ( var_B, 2.2);
+            r = pow_F(r, 2.2f);
+            g = pow_F(g, 2.2f);
+            b = pow_F(b, 2.2f);
+        }
+    } else { //display output profile
+        if (profile == "RT_sRGB" || profile == "RT_sRGB_gBT709" || profile == "RT_sRGB_g10") {
+            // use default "sRGB"
+        } else if (profile == "ProPhoto" || profile == "RT_Large_gBT709" || profile == "RT_Large_g10"  || profile == "RT_Large_gsRGB") {
+            profileCalc = "ProPhoto";
+        } else if (profile == "AdobeRGB1998" || profile == "RT_Medium_gsRGB") {
+            profileCalc = "Adobe RGB";
+        } else if (profile == "WideGamutRGB") {
+            profileCalc = "WideGamut";
         }
-    } else { //display outout profile
 
         if (profile == "RT_sRGB" || profile == "RT_Large_gsRGB"  || profile == "RT_Medium_gsRGB") { //apply sRGB inverse gamma
-            if ( var_R > 0.04045 ) {
-                var_R = pow ( ( ( var_R + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
+            if (r > 0.04045f) {
+                r = pow_F(((r + 0.055f) / 1.055f), rtengine::Color::sRGBGammaCurve);
             } else {
-                var_R = var_R / 12.92;
+                r /= 12.92f;
             }
 
-            if ( var_G > 0.04045 ) {
-                var_G = pow ( ( ( var_G + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
+            if (g > 0.04045f) {
+                g = pow_F(((g + 0.055f) / 1.055f), rtengine::Color::sRGBGammaCurve);
             } else {
-                var_G = var_G / 12.92;
+                g /= 12.92f;
             }
 
-            if ( var_B > 0.04045 ) {
-                var_B = pow ( ( ( var_B + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
+            if (b > 0.04045f) {
+                b = pow_F(((b + 0.055f) / 1.055f), rtengine::Color::sRGBGammaCurve);
             } else {
-                var_B = var_B / 12.92;
+                b /= 12.92f;
             }
-        }
-
-        else if (profile == "RT_sRGB_gBT709"  || profile == "RT_Large_gBT709") { //
-            if ( var_R > 0.0795 ) {
-                var_R = pow ( ( ( var_R + 0.0954 ) / 1.0954 ), 2.2);
+        } else if (profile == "RT_sRGB_gBT709"  || profile == "RT_Large_gBT709" || profile == "Rec2020") {
+            if (r > 0.0795f) {
+                r = pow_F(((r + 0.0954f) / 1.0954f), 2.2f);
             } else {
-                var_R = var_R / 4.5;
+                r /= 4.5f;
             }
 
-            if ( var_G > 0.0795 ) {
-                var_G = pow ( ( ( var_G + 0.0954 ) / 1.0954 ), 2.2);
+            if (g > 0.0795f) {
+                g = pow_F(((g + 0.0954f) / 1.0954f), 2.2f);
             } else {
-                var_G = var_G / 4.5;
+                g /= 4.5f;
             }
 
-            if ( var_B > 0.0795 ) {
-                var_B = pow ( ( ( var_B + 0.0954 ) / 1.0954 ), 2.2);
+            if (b > 0.0795f) {
+                b = pow_F(((b + 0.0954f) / 1.0954f), 2.2f);
             } else {
-                var_B = var_B / 4.5;
+                b /= 4.5f;
             }
 
         } else if (profile == "ProPhoto") { // apply inverse gamma 1.8
 
-            var_R = pow ( var_R, 1.8);
-            var_G = pow ( var_G, 1.8);
-            var_B = pow ( var_B, 1.8);
-        } else if (profile == "RT_sRGB_g10"  || profile == "RT_Large_g10") { // apply inverse gamma 1.8
+            r = pow_F(r, 1.8f);
+            g = pow_F(g, 1.8f);
+            b = pow_F(b, 1.8f);
+        } else if (profile == "RT_sRGB_g10"  || profile == "RT_Large_g10") {
+            // gamma 1.0, do nothing
 
-            var_R = pow ( var_R, 1.);
-            var_G = pow ( var_G, 1.);
-            var_B = pow ( var_B, 1.);
-        }
+        } else {// apply inverse gamma 2.2
 
-        else {// apply inverse gamma 2.2
-            var_R = pow ( var_R, 2.2);
-            var_G = pow ( var_G, 2.2);
-            var_B = pow ( var_B, 2.2);
+            r = pow_F(r, 2.2f);
+            g = pow_F(g, 2.2f);
+            b = pow_F(b, 2.2f);
         }
     }
 
-    // TMatrix wprof = rtengine::ICCStore::getInstance()->workingSpaceMatrix (profileW);
+    const TMatrix wprof = rtengine::ICCStore::getInstance()->workingSpaceMatrix(profileCalc);
 
-    TMatrix wprof = rtengine::ICCStore::getInstance()->workingSpaceMatrix (profileCalc);
+    const float xyz_rgb[3][3] = { {static_cast<float>(wprof[0][0]), static_cast<float>(wprof[0][1]), static_cast<float>(wprof[0][2])},
+                                  {static_cast<float>(wprof[1][0]), static_cast<float>(wprof[1][1]), static_cast<float>(wprof[1][2])},
+                                  {static_cast<float>(wprof[2][0]), static_cast<float>(wprof[2][1]), static_cast<float>(wprof[2][2])}
+                                };
 
-    for (int m = 0; m < 3; m++)
-        for (int n = 0; n < 3; n++) {
-            xyz_rgb[m][n] = wprof[m][n];
-        }
+    const float var_X = (xyz_rgb[0][0] * r + xyz_rgb[0][1] * g + xyz_rgb[0][2] * b) / Color::D50x;
+    const float var_Y = (xyz_rgb[1][0] * r + xyz_rgb[1][1] * g + xyz_rgb[1][2] * b);
+    const float var_Z = (xyz_rgb[2][0] * r + xyz_rgb[2][1] * g + xyz_rgb[2][2] * b) / Color::D50z;
 
-    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 ) / 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 ) / Color::D50z;
+    const float varxx = var_X > epsf ? xcbrtf(var_X) : (kappaf * var_X  + 16.f) / 116.f ;
+    const float varyy = var_Y > epsf ? xcbrtf(var_Y) : (kappaf * var_Y  + 16.f) / 116.f ;
+    const float varzz = var_Z > epsf ? xcbrtf(var_Z) : (kappaf * var_Z  + 16.f) / 116.f ;
 
-    varxx = var_X > ep ? cbrt(var_X) : ( ka * var_X  +  16.0) / 116.0 ;
-    varyy = var_Y > ep ? cbrt(var_Y) : ( ka * var_Y  +  16.0) / 116.0 ;
-    varzz = var_Z > ep ? cbrt(var_Z) : ( ka * var_Z  +  16.0) / 116.0 ;
-    LAB_l = ( 116 * varyy ) - 16;
-    LAB_a = 500 * ( varxx - varyy );
-    LAB_b = 200 * ( varyy - varzz );
+    LAB_l = var_Y > epsf ? (xcbrtf(var_Y) * 116.f) - 16.f : kappaf * var_Y;
+    LAB_a = 500.f * (varxx - varyy);
+    LAB_b = 200.f * (varyy - varzz);
 
 }
 
@@ -746,6 +741,36 @@ void Color::rgb2hsv(float r, float g, float b, float &h, float &s, float &v)
     }
 }
 
+void Color::rgb2hsv01(float r, float g, float b, float &h, float &s, float &v)
+{
+
+    const float minVal = min(r, g, b);
+    v = max(r, g, b);
+    const float delta = v - minVal;
+
+    h = 0.f;
+
+    if (delta < 0.00001f) {
+        s = 0.f;
+    } else {
+        s = delta / (v == 0.f ? 1.f : v);
+
+        if (r == v) {
+            h = (g - b) / delta;
+        } else if (g == v) {
+            h = 2.f + (b - r) / delta;
+        } else if (b == v) {
+            h = 4.f + (r - g) / delta;
+        }
+
+        h /= 6.f;
+
+        if (h < 0.f) {
+            h += 1.f;
+        }
+    }
+}
+
 void Color::hsv2rgb (float h, float s, float v, float &r, float &g, float &b)
 {
 
@@ -1683,7 +1708,7 @@ void Color::Lab2XYZ(float L, float a, float b, float &x, float &y, float &z)
     float LL = L / 327.68f;
     float aa = a / 327.68f;
     float bb = b / 327.68f;
-    float fy = (0.00862069f * LL) + 0.137932f; // (L+16)/116
+    float fy = (c1By116 * LL) + c16By116; // (L+16)/116
     float fx = (0.002f * aa) + fy;
     float fz = fy - (0.005f * bb);
     x = 65535.0f * f2xyz(fx) * D50x;
@@ -1694,7 +1719,7 @@ void Color::Lab2XYZ(float L, float a, float b, float &x, float &y, float &z)
 void Color::L2XYZ(float L, float &x, float &y, float &z) // for black & white
 {
     float LL = L / 327.68f;
-    float fy = (0.00862069f * LL) + 0.137932f; // (L+16)/116
+    float fy = (c1By116 * LL) + c16By116; // (L+16)/116
     float fxz = 65535.f * f2xyz(fy);
     x = fxz * D50x;
     z = fxz * D50z;
@@ -1709,7 +1734,7 @@ void Color::Lab2XYZ(vfloat L, vfloat a, vfloat b, vfloat &x, vfloat &y, vfloat &
     L /= c327d68;
     a /= c327d68;
     b /= c327d68;
-    vfloat fy = F2V(0.00862069f) * L + F2V(0.137932f);
+    vfloat fy = F2V(c1By116) * L + F2V(c16By116);
     vfloat fx = F2V(0.002f) * a + fy;
     vfloat fz = fy - (F2V(0.005f) * b);
     vfloat c65535 = F2V(65535.f);
@@ -1727,8 +1752,6 @@ void Color::RGB2Lab(float *R, float *G, float *B, float *L, float *a, float *b,
 
 #ifdef __SSE2__
     vfloat maxvalfv = F2V(MAXVALF);
-    vfloat c116v = F2V(116.f);
-    vfloat c5242d88v = F2V(5242.88f);
     vfloat c500v = F2V(500.f);
     vfloat c200v = F2V(200.f);
 #endif
@@ -1753,7 +1776,7 @@ void Color::RGB2Lab(float *R, float *G, float *B, float *L, float *a, float *b,
                 float fy = (y <= 65535.f ? cachef[y] : (327.68f * xcbrtf(y / MAXVALF)));
                 float fz = (z <= 65535.f ? cachef[z] : (327.68f * xcbrtf(z / MAXVALF)));
 
-                L[i + k] = (116.f *  fy - 5242.88f); //5242.88=16.0*327.68;
+                L[i + k] = (y <= 65535.0f ? cachefy[y] : 327.68f * (116.f * xcbrtf(y / MAXVALF) - 16.f));
                 a[i + k] = (500.f * (fx - fy) );
                 b[i + k] = (200.f * (fy - fz) );
             }
@@ -1762,7 +1785,7 @@ void Color::RGB2Lab(float *R, float *G, float *B, float *L, float *a, float *b,
             const vfloat fy = cachef[yv];
             const vfloat fz = cachef[zv];
 
-            STVFU(L[i], c116v *  fy - c5242d88v); //5242.88=16.0*327.68;
+            STVFU(L[i], cachefy[yv]);
             STVFU(a[i], c500v * (fx - fy));
             STVFU(b[i], c200v * (fy - fz));
         }
@@ -1781,7 +1804,7 @@ void Color::RGB2Lab(float *R, float *G, float *B, float *L, float *a, float *b,
         fy = (y <= 65535.0f ? cachef[y] : (327.68f * xcbrtf(y / MAXVALF)));
         fz = (z <= 65535.0f ? cachef[z] : (327.68f * xcbrtf(z / MAXVALF)));
 
-        L[i] = 116.0f *  fy - 5242.88f; //5242.88=16.0*327.68;
+        L[i] = (y <= 65535.0f ? cachefy[y] : 327.68f * (116.f * xcbrtf(y / MAXVALF) - 16.f));
         a[i] = 500.0f * (fx - fy);
         b[i] = 200.0f * (fy - fz);
     }
@@ -1799,14 +1822,14 @@ void Color::XYZ2Lab(float X, float Y, float Z, float &L, float &a, float &b)
     fy = (y <= 65535.0f ? cachef[y] : (327.68f * xcbrtf(y / MAXVALF)));
     fz = (z <= 65535.0f ? cachef[z] : (327.68f * xcbrtf(z / MAXVALF)));
 
-    L = (116.0f *  fy - 5242.88f); //5242.88=16.0*327.68;
+    L = (y <= 65535.0f ? cachefy[y] : 327.68f * (116.f * xcbrtf(y / MAXVALF) - 16.f));
     a = (500.0f * (fx - fy) );
     b = (200.0f * (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 fy = (c1By116 * L / 327.68) + c16By116; // (L+16)/116
     float fx = (0.002 * a / 327.68) + fy;
     float fz = fy - (0.005 * b / 327.68);
     float LL = L / 327.68;
@@ -1835,7 +1858,7 @@ void Color::Yuv2Lab(float Yin, float u, float v, float &L, float &a, float &b, c
     float fy = (Y <= 65535.0 ? cachef[Y] : (327.68 * std::cbrt(Y / MAXVALF)));
     float fz = (Z <= 65535.0 ? cachef[Z] : (327.68 * std::cbrt(Z / MAXVALF)));
 
-    L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68;
+    L = (Y <= 65535.0f ? cachefy[Y] : 327.68f * (116.f * xcbrtf(Y / MAXVALF) - 16.f));
     a = (500.0 * (fx - fy) );
     b = (200.0 * (fy - fz) );
 }
@@ -2304,7 +2327,7 @@ void Color::gamutLchonly (float HH, float &Lprov1, float &Chprov1, float &R, flo
         float bprov1 = Chprov1 * sincosval.x;
 
         //conversion Lab RGB to limit Lab values - this conversion is useful before Munsell correction
-        float fy = (0.00862069f * Lprov1 ) + 0.137932f;
+        float fy = (c1By116 * Lprov1 ) + c16By116;
         float fx = (0.002f * aprov1) + fy;
         float fz = fy - (0.005f * bprov1);
 
@@ -2424,7 +2447,7 @@ void Color::gamutLchonly (float HH, float2 sincosval, float &Lprov1, float &Chpr
         float bprov1 = Chprov1 * sincosval.x;
 
         //conversion Lab RGB to limit Lab values - this conversion is useful before Munsell correction
-        float fy = (0.00862069f * Lprov1 ) + 0.137932f;
+        float fy = (c1By116 * Lprov1 ) + c16By116;
         float fx = (0.002f * aprov1) + fy;
         float fz = fy - (0.005f * bprov1);
 
@@ -2533,7 +2556,7 @@ void Color::gamutLchonly (float2 sincosval, float &Lprov1, float &Chprov1, const
         float bprov1 = Chprov1 * sincosval.x;
 
         //conversion Lab RGB to limit Lab values - this conversion is useful before Munsell correction
-        float fy = (0.00862069f * Lprov1 ) + 0.137932f;
+        float fy = (c1By116 * Lprov1 ) + c16By116;
         float fx = (0.002f * aprov1) + fy;
         float fz = fy - (0.005f * bprov1);
 
diff --git a/rtengine/color.h b/rtengine/color.h
index d56b30e52..4b865f1d9 100644
--- a/rtengine/color.h
+++ b/rtengine/color.h
@@ -120,15 +120,36 @@ public:
         ID_DOWN      /// Interpolate color by decreasing the hue value, crossing the lower limit
     } eInterpolationDirection;
 
-    const static double sRGBGamma;        // standard average gamma
-    const static double sRGBGammaCurve;   // 2.4 in the curve
-    const static double eps, eps_max, kappa, epskap;
-    const static float D50x, D50z;
-    const static double u0, v0;
+    // 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.
+    constexpr static double sRGBGamma = 2.2;
+    constexpr static double sRGBGammaCurve = 2.4;
+
+    constexpr static double eps = 216.0 / 24389.0; //0.008856
+    constexpr static double eps_max = MAXVALF * eps; //580.40756;
+    constexpr static double kappa = 24389.0 / 27.0; //903.29630;
+    constexpr static double kappaInv = 27.0 / 24389.0;
+    constexpr static double epsilonExpInv3 = 6.0 / 29.0;
+
+    constexpr static float epsf = eps;
+    constexpr static float kappaf = kappa;
+    constexpr static float kappaInvf = kappaInv;
+    constexpr static float epsilonExpInv3f = epsilonExpInv3;
+
+    constexpr static float D50x = 0.9642f; //0.96422;
+    constexpr static float D50z = 0.8249f; //0.82521;
+    constexpr static double u0 = 4.0 * D50x / (D50x + 15 + 3 * D50z);
+    constexpr static double v0 = 9.0 / (D50x + 15 + 3 * D50z);
+    constexpr static double epskap = 8.0;
+
+    constexpr static float c1By116 = 1.0 / 116.0;
+    constexpr static float c16By116 = 16.0 / 116.0;
 
     static cmsToneCurve* linearGammaTRC;
 
     static LUTf cachef;
+    static LUTf cachefy;
     static LUTf gamma2curve;
 
     // look-up tables for the standard srgb gamma and its inverse (filled by init())
@@ -184,16 +205,16 @@ public:
     * @brief Convert red/green/blue to hue/saturation/luminance
     * @param profile output profile name
     * @param profileW working profile name
-    * @param r red channel [0 ; 65535]
-    * @param g green channel [0 ; 65535]
-    * @param b blue channel [0 ; 65535]
+    * @param r red channel [0 ; 1]
+    * @param g green channel [0 ; 1]
+    * @param b blue channel [0 ; 1]
     * @param L Lab L channel [0 ; 1] (return value)
-    * @param a Lab  a channel [0 ; 1] (return value)
+    * @param a Lab a channel [0 ; 1] (return value)
     * @param b Lab b channel [0; 1] (return value)
     * @param workingSpace true: compute the Lab value using the Working color space ; false: use the Output color space
     */
-    static void rgb2lab (Glib::ustring profile, Glib::ustring profileW, int r, int g, int b, float &LAB_l, float &LAB_a, float &LAB_b, bool workingSpace);
-
+    // do not use this function in a loop. It really eats processing time caused by Glib::ustring comparisons
+    static void rgb2lab01 (const Glib::ustring &profile, const Glib::ustring &profileW, float r, float g, float b, float &LAB_l, float &LAB_a, float &LAB_b, bool workingSpace);
 
     /**
     * @brief Convert red/green/blue to hue/saturation/luminance
@@ -327,6 +348,17 @@ public:
     */
     static void rgb2hsv (float r, float g, float b, float &h, float &s, float &v);
 
+    /**
+    * @brief Convert red green blue to hue saturation value
+    * @param r red channel [0 ; 1]
+    * @param g green channel [0 ; 1]
+    * @param b blue channel [0 ; 1]
+    * @param h hue channel [0 ; 1] (return value)
+    * @param s saturation channel [0 ; 1] (return value)
+    * @param v value channel [0 ; 1] (return value)
+    */
+    static void rgb2hsv01 (float r, float g, float b, float &h, float &s, float &v);
+
     static inline float rgb2s(float r, float g, float b) // fast version if only saturation is needed
     {
         float var_Min = min(r, g, b);
@@ -675,27 +707,21 @@ public:
     */
     static inline double 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;
 
     }
     static inline float f2xyz(float f)
     {
-        const float epsilonExpInv3 = 0.20689655f; // 6.0f/29.0f;
-        const float kappaInv = 0.0011070565f; // 27.0f/24389.0f;  // inverse of kappa
-
-        return (f > epsilonExpInv3) ? f * f * f : (116.f * f - 16.f) * kappaInv;
+        return (f > epsilonExpInv3f) ? f * f * f : (116.f * f - 16.f) * kappaInvf;
     }
 #ifdef __SSE2__
     static inline vfloat f2xyz(vfloat f)
     {
-        const vfloat epsilonExpInv3 = F2V(0.20689655f); // 6.0f/29.0f;
-        const vfloat kappaInv = F2V(0.0011070565f); // 27.0f/24389.0f;  // inverse of kappa
+        const vfloat epsilonExpInv3v = F2V(epsilonExpInv3f);
+        const vfloat kappaInvv = F2V(kappaInvf);
         vfloat res1 = f * f * f;
-        vfloat res2 = (F2V(116.f) * f - F2V(16.f)) * kappaInv;
-        return vself(vmaskf_gt(f, epsilonExpInv3), res1, res2);
+        vfloat res2 = (F2V(116.f) * f - F2V(16.f)) * kappaInvv;
+        return vself(vmaskf_gt(f, epsilonExpInv3v), res1, res2);
     }
 #endif
 
diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc
index 90a202a1c..1bbea204c 100644
--- a/rtengine/demosaic_algos.cc
+++ b/rtengine/demosaic_algos.cc
@@ -3861,7 +3861,7 @@ void RawImageSource::cielab (const float (*rgb)[3], float* l, float* a, float *b
         if(!cbrtinit) {
             for (int i = 0; i < 0x14000; i++) {
                 double r = i / 65535.0;
-                cbrt[i] = r > 0.008856f ? std::cbrt(r) : 7.787f * r + 16.f / 116.f;
+                cbrt[i] = r > Color::eps ? std::cbrt(r) : (Color::kappa * r + 16.0) / 116.0;
             }
 
             cbrtinit = true;
@@ -3871,7 +3871,6 @@ void RawImageSource::cielab (const float (*rgb)[3], float* l, float* a, float *b
     }
 
 #ifdef __SSE2__
-    vfloat zd5v = F2V(0.5f);
     vfloat c116v = F2V(116.f);
     vfloat c16v = F2V(16.f);
     vfloat c500v = F2V(500.f);
@@ -3892,12 +3891,12 @@ void RawImageSource::cielab (const float (*rgb)[3], float* l, float* a, float *b
         for(; j < labWidth - 3; j += 4) {
             vfloat redv, greenv, bluev;
             vconvertrgbrgbrgbrgb2rrrrggggbbbb(rgb[i * width + j], redv, greenv, bluev);
-            vfloat xyz0v = zd5v + redv * xyz_camv[0][0] + greenv * xyz_camv[0][1] + bluev * xyz_camv[0][2];
-            vfloat xyz1v = zd5v + redv * xyz_camv[1][0] + greenv * xyz_camv[1][1] + bluev * xyz_camv[1][2];
-            vfloat xyz2v = zd5v + redv * xyz_camv[2][0] + greenv * xyz_camv[2][1] + bluev * xyz_camv[2][2];
-            xyz0v = cbrt[_mm_cvttps_epi32(xyz0v)];
-            xyz1v = cbrt[_mm_cvttps_epi32(xyz1v)];
-            xyz2v = cbrt[_mm_cvttps_epi32(xyz2v)];
+            vfloat xyz0v = redv * xyz_camv[0][0] + greenv * xyz_camv[0][1] + bluev * xyz_camv[0][2];
+            vfloat xyz1v = redv * xyz_camv[1][0] + greenv * xyz_camv[1][1] + bluev * xyz_camv[1][2];
+            vfloat xyz2v = redv * xyz_camv[2][0] + greenv * xyz_camv[2][1] + bluev * xyz_camv[2][2];
+            xyz0v = cbrt[_mm_cvtps_epi32(xyz0v)];
+            xyz1v = cbrt[_mm_cvtps_epi32(xyz1v)];
+            xyz2v = cbrt[_mm_cvtps_epi32(xyz2v)];
 
             STVFU(l[i * labWidth + j], c116v * xyz1v - c16v);
             STVFU(a[i * labWidth + j], c500v * (xyz0v - xyz1v));
diff --git a/rtengine/improcfun.cc b/rtengine/improcfun.cc
index a76d4cd84..1a545bdf5 100644
--- a/rtengine/improcfun.cc
+++ b/rtengine/improcfun.cc
@@ -3877,8 +3877,7 @@ void ImProcFunctions::rgbProc (Imagefloat* working, LabImage* lab, PipetteBuffer
                                 // Luminosity after
                                 // only Luminance in Lab
                                 float newy = toxyz[1][0] * r + toxyz[1][1] * g + toxyz[1][2] * b;
-                                float newfy = newy < MAXVALF ? Color::cachef[newy] : 327.68f * std::cbrt (newy / MAXVALF);
-                                float L_2 = 116.0f * newfy - 5242.88f;
+                                float L_2 = newy <= MAXVALF ? Color::cachefy[newy] : 327.68f * (116.f * xcbrtf(newy / MAXVALF) - 16.f);
 
                                 //gamut control
                                 if (settings->rgbcurveslumamode_gamut) {
@@ -5775,7 +5774,7 @@ void ImProcFunctions::chromiLuminanceCurve (PipetteBuffer *pipetteBuffer, int pW
                         float aprov1 = Chprov2 * sincosval.y;
                         float bprov1 = Chprov2 * sincosval.x;
 
-                        float fy = (0.00862069f * Lprov1 ) + 0.137932f;
+                        float fy = (Color::c1By116 * Lprov1 ) + Color::c1By116;
                         float fx = (0.002f * aprov1) + fy;
                         float fz = fy - (0.005f * bprov1);
 
diff --git a/rtengine/iplab2rgb.cc b/rtengine/iplab2rgb.cc
index 27d9982fb..6979224a0 100644
--- a/rtengine/iplab2rgb.cc
+++ b/rtengine/iplab2rgb.cc
@@ -318,7 +318,7 @@ Imagefloat* ImProcFunctions::lab2rgbOut (LabImage* lab, int cx, int cy, int cw,
 
             for (int j = cx; j < cx + cw; j++) {
 
-                float fy = (0.0086206897f * rL[j]) / 327.68f + 0.1379310345f; // (L+16)/116
+                float fy = (Color::c1By116 * rL[j]) / 327.68f + Color::c16By116; // (L+16)/116
                 float fx = (0.002f * ra[j]) / 327.68f + fy;
                 float fz = fy - (0.005f * rb[j]) / 327.68f;
                 float LL = rL[j] / 327.68f;
diff --git a/rtengine/ipvibrance.cc b/rtengine/ipvibrance.cc
index f17d5d9b1..d7fb0da97 100644
--- a/rtengine/ipvibrance.cc
+++ b/rtengine/ipvibrance.cc
@@ -700,7 +700,7 @@ void ImProcFunctions::vibrance (LabImage* lab)
                     aprovn = Chprov * sincosval.y;
                     bprovn = Chprov * sincosval.x;
 
-                    float fyy = (0.00862069f * Lprov ) + 0.137932f;
+                    float fyy = (Color::c1By116 * Lprov ) + Color::c16By116;
                     float fxx = (0.002f * aprovn) + fyy;
                     float fzz = fyy - (0.005f * bprovn);
                     float xx_ = 65535.f * Color::f2xyz (fxx) * Color::D50x;
diff --git a/rtgui/histogrampanel.cc b/rtgui/histogrampanel.cc
index 51e288951..c713b51a8 100644
--- a/rtgui/histogrampanel.cc
+++ b/rtgui/histogrampanel.cc
@@ -476,7 +476,7 @@ void HistogramRGBArea::updateFreeze (bool f)
     return;
 }
 
-void HistogramRGBArea::updateBackBuffer (int r, int g, int b, Glib::ustring profile, Glib::ustring profileW)
+void HistogramRGBArea::updateBackBuffer (int r, int g, int b, const Glib::ustring &profile, const Glib::ustring &profileW)
 {
     if (!get_realized () || frozen || !showMode) {
         return;
@@ -531,7 +531,7 @@ void HistogramRGBArea::updateBackBuffer (int r, int g, int b, Glib::ustring prof
 
             if(needLuma || needChroma) {
                 float Lab_L, Lab_a, Lab_b;
-                rgb2lab( profile, profileW, r, g, b, Lab_L, Lab_a, Lab_b);
+                rtengine::Color::rgb2lab01(profile, profileW, r / 255.f, g / 255.f, b / 255.f, Lab_L, Lab_a, Lab_b, options.rtSettings.HistogramWorking);
 
                 if (needLuma) {
                     // Luma
@@ -557,153 +557,6 @@ void HistogramRGBArea::updateBackBuffer (int r, int g, int b, Glib::ustring prof
     setDirty(false);
 }
 
-void HistogramRGBArea::rgb2lab (Glib::ustring profile, Glib::ustring profileW, int r, int g, int b, float &LAB_l, float &LAB_a, float &LAB_b)
-{
-    double xyz_rgb[3][3];
-    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;
-    double var_B = b / 255.0;
-
-    Glib::ustring profileCalc;
-    profileCalc = "sRGB"; //default
-
-    if(options.rtSettings.HistogramWorking) {
-        profileCalc = profileW;    //display working
-    }
-
-    else {// if you want display = output space
-        if (profile == "RT_sRGB" || profile == "RT_sRGB_gBT709" || profile == "RT_sRGB_g10") {
-            profileCalc = "sRGB";
-        }
-
-        if (profile == "ProPhoto" || profile == "RT_Large_gBT709" || profile == "RT_Large_g10"  || profile == "RT_Large_gsRGB") {
-            profileCalc = "ProPhoto";
-        }
-
-        if (profile == "AdobeRGB1998" || profile == "RT_Medium_gsRGB") {
-            profileCalc = "Adobe RGB";
-        }
-
-        if (profile == "WideGamutRGB") {
-            profileCalc = "WideGamut";
-        }
-    }
-
-    if(options.rtSettings.HistogramWorking) {//display working
-        if (profileW == "sRGB") { //apply sRGB inverse gamma
-
-            if ( var_R > 0.04045 ) {
-                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::Color::sRGBGammaCurve);
-            } else {
-                var_G = var_G / 12.92;
-            }
-
-            if ( var_B > 0.04045 ) {
-                var_B = pow ( ( ( var_B + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
-            } else {
-                var_B = var_B / 12.92;
-            }
-        } else if (profileW == "ProPhoto") { // apply inverse gamma 1.8
-            var_R = pow ( var_R, 1.8);
-            var_G = pow ( var_G, 1.8);
-            var_B = pow ( var_B, 1.8);
-        } else { // apply inverse gamma 2.2
-            var_R = pow ( var_R, 2.2);
-            var_G = pow ( var_G, 2.2);
-            var_B = pow ( var_B, 2.2);
-        }
-    } else { //display outout profile
-
-        if (profile == "RT_sRGB" || profile == "RT_Large_gsRGB"  || profile == "RT_Medium_gsRGB") { //apply sRGB inverse gamma
-            if ( var_R > 0.04045 ) {
-                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::Color::sRGBGammaCurve);
-            } else {
-                var_G = var_G / 12.92;
-            }
-
-            if ( var_B > 0.04045 ) {
-                var_B = pow ( ( ( var_B + 0.055 ) / 1.055 ), rtengine::Color::sRGBGammaCurve);
-            } else {
-                var_B = var_B / 12.92;
-            }
-        }
-
-        else if (profile == "RT_sRGB_gBT709"  || profile == "RT_Large_gBT709") { //
-            if ( var_R > 0.0795 ) {
-                var_R = pow ( ( ( var_R + 0.0954 ) / 1.0954 ), 2.2);
-            } else {
-                var_R = var_R / 4.5;
-            }
-
-            if ( var_G > 0.0795 ) {
-                var_G = pow ( ( ( var_G + 0.0954 ) / 1.0954 ), 2.2);
-            } else {
-                var_G = var_G / 4.5;
-            }
-
-            if ( var_B > 0.0795 ) {
-                var_B = pow ( ( ( var_B + 0.0954 ) / 1.0954 ), 2.2);
-            } else {
-                var_B = var_B / 4.5;
-            }
-
-        } else if (profile == "ProPhoto") { // apply inverse gamma 1.8
-
-            var_R = pow ( var_R, 1.8);
-            var_G = pow ( var_G, 1.8);
-            var_B = pow ( var_B, 1.8);
-        } else if (profile == "RT_sRGB_g10"  || profile == "RT_Large_g10") { // apply inverse gamma 1.8
-
-            var_R = pow ( var_R, 1.);
-            var_G = pow ( var_G, 1.);
-            var_B = pow ( var_B, 1.);
-        }
-
-        else {// apply inverse gamma 2.2
-            var_R = pow ( var_R, 2.2);
-            var_G = pow ( var_G, 2.2);
-            var_B = pow ( var_B, 2.2);
-        }
-    }
-
-    // TMatrix wprof = rtengine::ICCStore::getInstance()->workingSpaceMatrix (profileW);
-
-    TMatrix wprof = rtengine::ICCStore::getInstance()->workingSpaceMatrix (profileCalc);
-
-    for (int m = 0; m < 3; m++)
-        for (int n = 0; n < 3; n++) {
-            xyz_rgb[m][n] = wprof[m][n];
-        }
-
-    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 ) / 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 ) / Color::D50z;
-
-    varxx = var_X > ep ? cbrt(var_X) : ( ka * var_X  +  16.0) / 116.0 ;
-    varyy = var_Y > ep ? cbrt(var_Y) : ( ka * var_Y  +  16.0) / 116.0 ;
-    varzz = var_Z > ep ? cbrt(var_Z) : ( ka * var_Z  +  16.0) / 116.0 ;
-    LAB_l = ( 116 * varyy ) - 16;
-    LAB_a = 500 * ( varxx - varyy );
-    LAB_b = 200 * ( varyy - varzz );
-
-}
-
 void HistogramRGBArea::update (int valh, int rh, int  gh, int bh)
 {
 
diff --git a/rtgui/histogrampanel.h b/rtgui/histogrampanel.h
index 25fa15c23..0ce62956c 100644
--- a/rtgui/histogrampanel.h
+++ b/rtgui/histogrampanel.h
@@ -76,7 +76,7 @@ public:
     HistogramRGBArea();
     ~HistogramRGBArea();
 
-    void updateBackBuffer (int r, int g, int b, Glib::ustring profile = "", Glib::ustring profileW = "");
+    void updateBackBuffer (int r, int g, int b, const Glib::ustring &profile = "", const Glib::ustring &profileW = "");
     void updateFreeze (bool f);
     bool getFreeze ();
     bool getShow ();
@@ -92,7 +92,6 @@ public:
     bool on_draw(const ::Cairo::RefPtr< Cairo::Context> &cr);
     bool on_button_press_event (GdkEventButton* event);
 private:
-    void rgb2lab (Glib::ustring profile, Glib::ustring profileW, int r, int g, int b, float &LAB_l, float &LAB_a, float &LAB_b);
     Gtk::SizeRequestMode get_request_mode_vfunc () const;
     void get_preferred_height_vfunc (int& minimum_height, int& natural_height) const;
     void get_preferred_width_vfunc (int &minimum_width, int &natural_width) const;
diff --git a/rtgui/lockablecolorpicker.cc b/rtgui/lockablecolorpicker.cc
index d10dfb438..859057e2c 100644
--- a/rtgui/lockablecolorpicker.cc
+++ b/rtgui/lockablecolorpicker.cc
@@ -277,8 +277,8 @@ void LockableColorPicker::setRGB (const float R, const float G, const float B, c
     gpreview = previewG;
     bpreview = previewB;
 
-    rtengine::Color::rgb2hsv(r*65535.f, g*65535.f, b*65535.f, hue, sat, val);
-    rtengine::Color::rgb2lab (*outputProfile, *workingProfile, r * 65535.f, g * 65535.f, b * 65535.f, L, a, bb, options.rtSettings.HistogramWorking);  // TODO: Really sure this function works?
+    rtengine::Color::rgb2hsv01(r, g, b, hue, sat, val);
+    rtengine::Color::rgb2lab01(*outputProfile, *workingProfile, r, g, b, L, a, bb, options.rtSettings.HistogramWorking);  // TODO: Really sure this function works?
 
     if (validity != Validity::OUTSIDE) {
         setDirty(true);
diff --git a/rtgui/navigator.cc b/rtgui/navigator.cc
index 7222436b6..7f159861d 100644
--- a/rtgui/navigator.cc
+++ b/rtgui/navigator.cc
@@ -290,13 +290,13 @@ void Navigator::pointerMoved (bool validPos, Glib::ustring profile, Glib::ustrin
         G->set_text (s2);
         B->set_text (s3);
 
-        Color::rgb2hsv (r * 0xffff / 0xff, g * 0xffff / 0xff, b * 0xffff / 0xff, h, s, v);
+        Color::rgb2hsv01(r / 255.f, g / 255.f, b / 255.f, h, s, v);
         getHSVText (h, s, v, s1, s2, s3);
         H->set_text (s1);
         S->set_text (s2);
         V->set_text (s3);
 
-        Color::rgb2lab (profile, profileW, r * 0xffff / 0xff, g * 0xffff / 0xff, b * 0xffff / 0xff, LAB_l, LAB_a, LAB_b, options.rtSettings.HistogramWorking);  // TODO: Really sure this function works?
+        Color::rgb2lab01(profile, profileW, r / 255.f, g / 255.f, b / 255.f, LAB_l, LAB_a, LAB_b, options.rtSettings.HistogramWorking);  // TODO: Really sure this function works?
         getLABText (LAB_l, LAB_a, LAB_b, s1, s2, s3);
         LAB_L->set_text (s1);
         LAB_A->set_text (s2);