diff --git a/rtdata/languages/default b/rtdata/languages/default
index cba87d663..c76d5da54 100644
--- a/rtdata/languages/default
+++ b/rtdata/languages/default
@@ -1018,6 +1018,7 @@ HISTORY_MSG_777;Local - Blur Mask local contrast curve
 HISTORY_MSG_778;Local - Blur Mask Wavelet level
 HISTORY_MSG_779;Local - Color Mask local contrast curve
 HISTORY_MSG_780;Local - Color Mask Wavelet level
+HISTORY_MSG_781;Local - Local contrast Wavelet level
 HISTORY_MSG_CLAMPOOG;Clip out-of-gamut colors
 HISTORY_MSG_COLORTONING_LABGRID_VALUE;CT - Color correction
 HISTORY_MSG_COLORTONING_LABREGION_AB;CT - Color correction
@@ -2332,6 +2333,7 @@ TP_LOCALLAB_LOCCONT;Unsharp Mask
 TP_LOCALLAB_WAV;Levels local contrast
 TP_LOCALLAB_WAVMASK;Mask Levels local contrast
 TP_LOCALLAB_LEVELWAV;Wavelets Levels
+TP_LOCALLAB_CSTHRESHOLD;Wavelets Levels
 TP_LOCALLAB_RESIDCONT;Residual image Contrast
 TP_LOCALLAB_RESIDCHRO;Residual image Chroma
 TP_LOCALLAB_CURVEEDITOR_TONES_LABEL;Tone curve
diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h
index df445c739..058f5d30e 100644
--- a/rtengine/improcfun.h
+++ b/rtengine/improcfun.h
@@ -256,6 +256,7 @@ public:
     static void strcurv_data(std::string retistr, int *s_datc, int &siz);
     void blendstruc(int bfw, int bfh, LabImage* bufcolorig, float radius, float stru, array2D<float> & blend2, int sk, bool multiThread);
     void wavcontrast(float ** tmp, float contrast, int bfw, int bfh, int wavelet_level, int sk, bool numThreads, const LocwavCurve & locwavCurve, bool & locwavutili, int & maxlvl);
+    void wavcontrast4(float ** tmp, float contrast, int bfw, int bfh, int level_bl, int level_hl, int level_br, int level_hr, int sk, bool numThreads, const LocwavCurve & locwavCurve, bool & locwavutili, int & maxlvl);
 
     void transit_shapedetect_retinex(int call, int senstype, LabImage * bufexporig, LabImage * bufmask, LabImage * buforigmas, float **buflight, float **bufchro, const float hueref, const float chromaref,  const float lumaref, const struct local_params & lp, LabImage * original, LabImage * transformed, int cx, int cy, int sk);
     void transit_shapedetect(int senstype, const LabImage *bufexporig, LabImage * originalmask, float **buflight, float **bufchro, float **buf_a_cat, float ** buf_b_cat, float ** bufhh, bool HHutili, const float hueref, const float chromaref, const float lumaref, float sobelref, float meansobel, float ** blend2, const struct local_params & lp, LabImage * original, LabImage * transformed, int cx, int cy, int sk);
diff --git a/rtengine/iplocallab.cc b/rtengine/iplocallab.cc
index 01241ef5f..d937bbed4 100644
--- a/rtengine/iplocallab.cc
+++ b/rtengine/iplocallab.cc
@@ -2984,36 +2984,38 @@ void ImProcFunctions::maskcalccol(bool invmask, bool pde, int bfw, int bfh, int
                 }
         }
 
-            int wavelet_level = wavlevel;
+        int wavelet_level = wavlevel;
 
-            int minwin = min(bfw, bfh);
-            int maxlevelspot = 9;
-            while ((1 << maxlevelspot) >= (minwin * sk) && maxlevelspot  > 1) {
-                --maxlevelspot ;
-            }
-            wavelet_level = min(wavelet_level, maxlevelspot);
-            int maxlvl;
-            float contrast = 0.f;
-            bool wavcurvemask = false;
+        int minwin = min(bfw, bfh);
+        int maxlevelspot = 9;
 
-            if (loclmasCurvecolwav && lmasutilicolwav) {
-                for (int i = 0; i < 500; i++) {
-                    if (loclmasCurvecolwav[i] != 0.5) {
-                        wavcurvemask = true;
-                    }
+        while ((1 << maxlevelspot) >= (minwin * sk) && maxlevelspot  > 1) {
+            --maxlevelspot ;
+        }
+
+        wavelet_level = min(wavelet_level, maxlevelspot);
+        int maxlvl;
+        float contrast = 0.f;
+        bool wavcurvemask = false;
+
+        if (loclmasCurvecolwav && lmasutilicolwav) {
+            for (int i = 0; i < 500; i++) {
+                if (loclmasCurvecolwav[i] != 0.5) {
+                    wavcurvemask = true;
                 }
             }
-            
-            if(wavcurvemask) {
+        }
+
+        if (wavcurvemask) {
 #ifdef _OPENMP
-                    const int numThreads = omp_get_max_threads();
+            const int numThreads = omp_get_max_threads();
 #else
-                    const int numThreads = 1;
+            const int numThreads = 1;
 
 #endif
-                
-                wavcontrast(bufmaskblurcol->L, contrast, bfw, bfh, wavelet_level, sk, numThreads, loclmasCurvecolwav, lmasutilicolwav, maxlvl);
-            }
+
+            wavcontrast(bufmaskblurcol->L, contrast, bfw, bfh, wavelet_level, sk, numThreads, loclmasCurvecolwav, lmasutilicolwav, maxlvl);
+        }
 
 
         if (lap > 0.f) {
@@ -5295,119 +5297,272 @@ void ImProcFunctions::fftw_tile_blur(int GW, int GH, int tilssize, int max_numbl
     fftwf_cleanup();
 }
 
-void ImProcFunctions::wavcontrast(float ** tmp, float contrast, int bfw, int bfh, int wavelet_level, int sk, bool numThreads, const LocwavCurve & locwavCurve, bool & locwavutili, int & maxlvl)
-{                    
-                    wavelet_decomposition *wdspot = new wavelet_decomposition(tmp[0], bfw, bfh, wavelet_level, 1, sk, numThreads, 6);
+void ImProcFunctions::wavcontrast4(float ** tmp, float contrast, int bfw, int bfh, int level_bl, int level_hl, int level_br, int level_hr, int sk, bool numThreads, const LocwavCurve & locwavCurve, bool & locwavutili, int & maxlvl)
+{
+    wavelet_decomposition *wdspot = new wavelet_decomposition(tmp[0], bfw, bfh, level_br, 1, sk, numThreads, 6);
 
-                    if (wdspot->memoryAllocationFailed) {
-                        return;
-                    }
+    if (wdspot->memoryAllocationFailed) {
+        return;
+    }
 
-                    maxlvl = wdspot->maxlevel();
+    maxlvl = wdspot->maxlevel();
 
-                    if (contrast != 0) {
-                        int W_L = wdspot->level_W(0);
-                        int H_L = wdspot->level_H(0);
-                        float *wav_L0 = wdspot->coeff0;
+    if (contrast != 0) {
+        int W_L = wdspot->level_W(0);
+        int H_L = wdspot->level_H(0);
+        float *wav_L0 = wdspot->coeff0;
 
-                        double avedbl = 0.0; // use double precision for large summations
+        double avedbl = 0.0; // use double precision for large summations
 
 #ifdef _OPENMP
-                        #pragma omp parallel for reduction(+:avedbl) if (multiThread)
+        #pragma omp parallel for reduction(+:avedbl) if (multiThread)
 #endif
 
-                        for (int i = 0; i < W_L * H_L; i++) {
-                            avedbl += wav_L0[i];
-                        }
+        for (int i = 0; i < W_L * H_L; i++) {
+            avedbl += wav_L0[i];
+        }
 
-                        float ave = avedbl / double(W_L * H_L);
+        float ave = avedbl / double(W_L * H_L);
 
-                        float avg = ave / 32768.f;
-                        avg = LIM01(avg);
-                        double contreal = 0.6 * contrast;
-                        DiagonalCurve resid_contrast({
-                            DCT_NURBS,
-                            0, 0,
-                            avg - avg * (0.6 - contreal / 250.0), avg - avg * (0.6 + contreal / 250.0),
-                            avg + (1. - avg) * (0.6 - contreal / 250.0), avg + (1. - avg) * (0.6 + contreal / 250.0),
-                            1, 1
-                        });
+        float avg = ave / 32768.f;
+        avg = LIM01(avg);
+        double contreal = 0.6 * contrast;
+        DiagonalCurve resid_contrast({
+            DCT_NURBS,
+            0, 0,
+            avg - avg * (0.6 - contreal / 250.0), avg - avg * (0.6 + contreal / 250.0),
+            avg + (1. - avg) * (0.6 - contreal / 250.0), avg + (1. - avg) * (0.6 + contreal / 250.0),
+            1, 1
+        });
 #ifdef _OPENMP
-                        #pragma omp parallel for if (multiThread)
+        #pragma omp parallel for if (multiThread)
 #endif
 
-                        for (int i = 0; i < W_L * H_L; i++) {
-                            float buf = LIM01(wav_L0[i] / 32768.f);
-                            buf = resid_contrast.getVal(buf);
-                            buf *= 32768.f;
-                            wav_L0[i] = buf;
-                        }
+        for (int i = 0; i < W_L * H_L; i++) {
+            float buf = LIM01(wav_L0[i] / 32768.f);
+            buf = resid_contrast.getVal(buf);
+            buf *= 32768.f;
+            wav_L0[i] = buf;
+        }
 
-                    }
+    }
 
 
-                    float mean[10];
-                    float meanN[10];
-                    float sigma[10];
-                    float sigmaN[10];
-                    float MaxP[10];
-                    float MaxN[10];
-                    Evaluate2(*wdspot, mean, meanN, sigma, sigmaN, MaxP, MaxN);
+    float mean[10];
+    float meanN[10];
+    float sigma[10];
+    float sigmaN[10];
+    float MaxP[10];
+    float MaxN[10];
+    Evaluate2(*wdspot, mean, meanN, sigma, sigmaN, MaxP, MaxN);
+    float alow = 1.f;
+    float blow = 0.f;
 
-                    if (locwavCurve && locwavutili) {
-                        for (int dir = 1; dir < 4; dir++) {
-                            for (int level = 0; level < maxlvl; ++level) {
-                                int W_L = wdspot->level_W(level);
-                                int H_L = wdspot->level_H(level);
-                                float **wav_L = wdspot->level_coeffs(level);
+    if (level_hl != level_bl) {
+    //    alow = 0.5f / (level_hl - level_bl);//to test with 0.5
+    //    blow = 0.5f -alow * level_bl;
+        alow = 1.f / (level_hl - level_bl);
+        blow = -alow * level_bl;
+    }
 
-                                if (MaxP[level] > 0.f && mean[level] != 0.f && sigma[level] != 0.f) {
-                                    float insigma = 0.666f; //SD
-                                    float logmax = log(MaxP[level]); //log Max
-                                    float rapX = (mean[level] + sigma[level]) / MaxP[level]; //rapport between sD / max
-                                    float inx = log(insigma);
-                                    float iny = log(rapX);
-                                    float rap = inx / iny; //koef
-                                    float asig = 0.166f / sigma[level];
-                                    float bsig = 0.5f - asig * mean[level];
-                                    float amean = 0.5f / mean[level];
+    float ahigh = 1.f;
+    float bhigh = 0.f;
+
+    if (level_hr != level_br) {
+     //   ahigh = 0.5f / (level_hr - level_br);//to test with 0.5
+     //   bhigh = 0.5f -ahigh * level_br;
+        ahigh = 1.f / (level_hr - level_br);
+        bhigh =  -ahigh * level_br;
+    }
+
+    if (locwavCurve && locwavutili) {
+        for (int dir = 1; dir < 4; dir++) {
+            for (int level = level_bl; level < maxlvl; ++level) {
+                int W_L = wdspot->level_W(level);
+                int H_L = wdspot->level_H(level);
+                float **wav_L = wdspot->level_coeffs(level);
+
+                if (MaxP[level] > 0.f && mean[level] != 0.f && sigma[level] != 0.f) {
+                    float insigma = 0.666f; //SD
+                    float logmax = log(MaxP[level]); //log Max
+                    float rapX = (mean[level] + sigma[level]) / MaxP[level]; //rapport between sD / max
+                    float inx = log(insigma);
+                    float iny = log(rapX);
+                    float rap = inx / iny; //koef
+                    float asig = 0.166f / sigma[level];
+                    float bsig = 0.5f - asig * mean[level];
+                    float amean = 0.5f / mean[level];
 
 #ifdef _OPENMP
-                                    #pragma omp parallel for if (multiThread)
+                    #pragma omp parallel for if (multiThread)
 #endif
 
-                                    for (int i = 0; i < W_L * H_L; i++) {
-                                        if (locwavCurve && locwavutili) {
-                                            float absciss;
-                                            float &val = wav_L[dir][i];
+                    for (int i = 0; i < W_L * H_L; i++) {
+                        if (locwavCurve && locwavutili) {
+                            float absciss;
+                            float &val = wav_L[dir][i];
 
-                                            if (fabsf(val) >= (mean[level] + sigma[level])) { //for max
-                                                float valcour = xlogf(fabsf(val));
-                                                float valc = valcour - logmax;
-                                                float vald = valc * rap;
-                                                absciss = xexpf(vald);
-                                            } else if (fabsf(val) >= mean[level]) {
-                                                absciss = asig * fabsf(val) + bsig;
-                                            } else {
-                                                absciss = amean * fabsf(val);
-                                            }
+                            if (fabsf(val) >= (mean[level] + sigma[level])) { //for max
+                                float valcour = xlogf(fabsf(val));
+                                float valc = valcour - logmax;
+                                float vald = valc * rap;
+                                absciss = xexpf(vald);
+                            } else if (fabsf(val) >= mean[level]) {
+                                absciss = asig * fabsf(val) + bsig;
+                            } else {
+                                absciss = amean * fabsf(val);
+                            }
 
-                                            float kc = locwavCurve[absciss * 500.f] - 0.5f;
-                                            float reduceeffect = kc <= 0.f ? 1.f : 1.5f;
+                            float klev = 1.f;
 
-                                            float kinterm = 1.f + reduceeffect * kc;
-                                            kinterm = kinterm <= 0.f ? 0.01f : kinterm;
+                            if (level >= level_hl && level <= level_hr) {
+                                klev = 1.f;
+                            }
 
-                                            val *=  kinterm;
-                                        }
-                                    }
+                            if (level_hl != level_bl) {
+                                if (level >= level_bl && level < level_hl) {
+                                    klev = alow * level + blow;
                                 }
                             }
+
+                            if (level_hr != level_br) {
+                                if (level > level_hr && level <= level_br) {
+                                    klev = ahigh * level + bhigh;
+                                }
+                            }
+
+                            float kc = klev * (locwavCurve[absciss * 500.f] - 0.5f);
+                            float reduceeffect = kc <= 0.f ? 1.f : 1.5f;
+
+                            float kinterm = 1.f + reduceeffect * kc;
+                            kinterm = kinterm <= 0.f ? 0.01f : kinterm;
+
+                            val *=  kinterm;
                         }
                     }
+                }
+            }
+        }
+    }
 
-                    wdspot->reconstruct(tmp[0], 1.f);
-                    delete wdspot;
+    wdspot->reconstruct(tmp[0], 1.f);
+    delete wdspot;
+}
+
+
+void ImProcFunctions::wavcontrast(float ** tmp, float contrast, int bfw, int bfh, int wavelet_level, int sk, bool numThreads, const LocwavCurve & locwavCurve, bool & locwavutili, int & maxlvl)
+{
+    wavelet_decomposition *wdspot = new wavelet_decomposition(tmp[0], bfw, bfh, wavelet_level, 1, sk, numThreads, 6);
+
+    if (wdspot->memoryAllocationFailed) {
+        return;
+    }
+
+    maxlvl = wdspot->maxlevel();
+
+    if (contrast != 0) {
+        int W_L = wdspot->level_W(0);
+        int H_L = wdspot->level_H(0);
+        float *wav_L0 = wdspot->coeff0;
+
+        double avedbl = 0.0; // use double precision for large summations
+
+#ifdef _OPENMP
+        #pragma omp parallel for reduction(+:avedbl) if (multiThread)
+#endif
+
+        for (int i = 0; i < W_L * H_L; i++) {
+            avedbl += wav_L0[i];
+        }
+
+        float ave = avedbl / double(W_L * H_L);
+
+        float avg = ave / 32768.f;
+        avg = LIM01(avg);
+        double contreal = 0.6 * contrast;
+        DiagonalCurve resid_contrast({
+            DCT_NURBS,
+            0, 0,
+            avg - avg * (0.6 - contreal / 250.0), avg - avg * (0.6 + contreal / 250.0),
+            avg + (1. - avg) * (0.6 - contreal / 250.0), avg + (1. - avg) * (0.6 + contreal / 250.0),
+            1, 1
+        });
+#ifdef _OPENMP
+        #pragma omp parallel for if (multiThread)
+#endif
+
+        for (int i = 0; i < W_L * H_L; i++) {
+            float buf = LIM01(wav_L0[i] / 32768.f);
+            buf = resid_contrast.getVal(buf);
+            buf *= 32768.f;
+            wav_L0[i] = buf;
+        }
+
+    }
+
+
+    float mean[10];
+    float meanN[10];
+    float sigma[10];
+    float sigmaN[10];
+    float MaxP[10];
+    float MaxN[10];
+    Evaluate2(*wdspot, mean, meanN, sigma, sigmaN, MaxP, MaxN);
+
+    if (locwavCurve && locwavutili) {
+        for (int dir = 1; dir < 4; dir++) {
+            for (int level = 0; level < maxlvl; ++level) {
+                int W_L = wdspot->level_W(level);
+                int H_L = wdspot->level_H(level);
+                float **wav_L = wdspot->level_coeffs(level);
+
+                if (MaxP[level] > 0.f && mean[level] != 0.f && sigma[level] != 0.f) {
+                    float insigma = 0.666f; //SD
+                    float logmax = log(MaxP[level]); //log Max
+                    float rapX = (mean[level] + sigma[level]) / MaxP[level]; //rapport between sD / max
+                    float inx = log(insigma);
+                    float iny = log(rapX);
+                    float rap = inx / iny; //koef
+                    float asig = 0.166f / sigma[level];
+                    float bsig = 0.5f - asig * mean[level];
+                    float amean = 0.5f / mean[level];
+
+#ifdef _OPENMP
+                    #pragma omp parallel for if (multiThread)
+#endif
+
+                    for (int i = 0; i < W_L * H_L; i++) {
+                        if (locwavCurve && locwavutili) {
+                            float absciss;
+                            float &val = wav_L[dir][i];
+
+                            if (fabsf(val) >= (mean[level] + sigma[level])) { //for max
+                                float valcour = xlogf(fabsf(val));
+                                float valc = valcour - logmax;
+                                float vald = valc * rap;
+                                absciss = xexpf(vald);
+                            } else if (fabsf(val) >= mean[level]) {
+                                absciss = asig * fabsf(val) + bsig;
+                            } else {
+                                absciss = amean * fabsf(val);
+                            }
+
+                            float kc = locwavCurve[absciss * 500.f] - 0.5f;
+                            float reduceeffect = kc <= 0.f ? 1.f : 1.5f;
+
+                            float kinterm = 1.f + reduceeffect * kc;
+                            kinterm = kinterm <= 0.f ? 0.01f : kinterm;
+
+                            val *=  kinterm;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    wdspot->reconstruct(tmp[0], 1.f);
+    delete wdspot;
 }
 
 
@@ -7122,23 +7277,26 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o
             }
 
 
-        if (lmaskbllocalcurve && localmaskblutili) {
+            if (lmaskbllocalcurve && localmaskblutili) {
 #ifdef _OPENMP
-            #pragma omp parallel for schedule(dynamic,16)
+                #pragma omp parallel for schedule(dynamic,16)
 #endif
 
-            for (int ir = 0; ir < GH; ir++)
-                for (int jr = 0; jr < GW; jr++) {
-                    bufmaskblurbl->L[ir][jr] = 0.5f * lmaskbllocalcurve[2.f * bufmaskblurbl->L[ir][jr]];
-                }
-        }
+                for (int ir = 0; ir < GH; ir++)
+                    for (int jr = 0; jr < GW; jr++) {
+                        bufmaskblurbl->L[ir][jr] = 0.5f * lmaskbllocalcurve[2.f * bufmaskblurbl->L[ir][jr]];
+                    }
+            }
+
             int wavelet_level = params->locallab.spots.at(sp).wavmaskbl;
 
             int minwin = min(GW, GH);
             int maxlevelspot = 9;
+
             while ((1 << maxlevelspot) >= (minwin * sk) && maxlevelspot  > 1) {
                 --maxlevelspot ;
             }
+
             wavelet_level = min(wavelet_level, maxlevelspot);
             int maxlvl;
             float contrast = 0.f;
@@ -7151,15 +7309,15 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o
                     }
                 }
             }
-            
-            if(wavcurvemask) {
+
+            if (wavcurvemask) {
 #ifdef _OPENMP
-                    const int numThreads = omp_get_max_threads();
+                const int numThreads = omp_get_max_threads();
 #else
-                    const int numThreads = 1;
+                const int numThreads = 1;
 
 #endif
-                
+
                 wavcontrast(bufmaskblurbl->L, contrast, GW, GH, wavelet_level, sk, numThreads, loclmasCurveblwav, lmasutiliblwav, maxlvl);
             }
 
@@ -7898,6 +8056,7 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o
                     vibranceParams.pastsattog = params->locallab.spots.at(sp).pastsattog;
                     vibranceParams.skintonescurve = params->locallab.spots.at(sp).skintonescurve;
 
+
                     bufexpfin->CopyFrom(bufexporig.get());
                     ImProcFunctions::vibrance(bufexpfin.get(), vibranceParams, params->toneCurve.hrenabled, params->icm.workingProfile);
 
@@ -8683,6 +8842,7 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o
                     // adap maximum level wavelet to size of RT-spot
                     int minwin = min(bfw, bfh);
                     int maxlevelspot = 9;
+
                     while ((1 << maxlevelspot) >= (minwin * sk) && maxlevelspot  > 1) {
                         --maxlevelspot ;
                     }
@@ -8814,10 +8974,16 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o
                         wdspotresidb->reconstruct(tmpresid->b[0], 1.f);
                         delete wdspotresidb;
                     }
+
                     int maxlvl;
                     const float contrast = params->locallab.spots.at(sp).residcont;
+                    int level_bl = params->locallab.spots.at(sp).csthreshold.getBottomLeft();
+                    int level_hl = params->locallab.spots.at(sp).csthreshold.getTopLeft();
+                    int level_br = params->locallab.spots.at(sp).csthreshold.getBottomRight();
+                    int level_hr = params->locallab.spots.at(sp).csthreshold.getTopRight();
 
-                    wavcontrast(tmp1->L, contrast, tmp1->W, tmp1->H, wavelet_level, sk, numThreads, locwavCurve, locwavutili, maxlvl);
+                    //     wavcontrast(tmp1->L, contrast, tmp1->W, tmp1->H, wavelet_level, sk, numThreads, locwavCurve, locwavutili, maxlvl);
+                    wavcontrast4(tmp1->L, contrast, tmp1->W, tmp1->H, level_bl, level_hl, level_br, level_hr, sk, numThreads, locwavCurve, locwavutili, maxlvl);
 
                     const float satur = params->locallab.spots.at(sp).residchro;
 
@@ -9309,7 +9475,7 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o
                 //for Retinex Mask are incorporated in MSR
                 ImProcFunctions::MSRLocal(sp, fftw, 1, reducDE, bufreti, bufmask, buforig, buforigmas, orig, tmpl->L, orig1,
                                           Wd, Hd, Wd, Hd, params->locallab, sk, locRETgainCcurve, locRETtransCcurve, 0, 4, 1.f, minCD, maxCD, mini, maxi, Tmean, Tsigma, Tmin, Tmax,
-                                          locccmasretiCurve, lcmasretiutili, locllmasretiCurve, llmasretiutili, lochhmasretiCurve, lhmasretiutili, llretiMask, 
+                                          locccmasretiCurve, lcmasretiutili, locllmasretiCurve, llmasretiutili, lochhmasretiCurve, lhmasretiutili, llretiMask,
                                           lmaskretilocalcurve, localmaskretiutili,
                                           transformed, lp.enaretiMasktmap, lp.enaretiMask);
 #ifdef _OPENMP
@@ -9438,11 +9604,11 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o
                     bool fftw = false;
 
                     if (params->locallab.spots.at(sp).chrrt > 40.f) { //second step active Retinex Chroma
-                        ImProcFunctions::MSRLocal(sp, fftw, 0, nullptr, bufreti, bufmask, buforig, buforigmas, orig, tmpl->L, orig1, 
-                                          Wd, Hd, Wd, Hd, params->locallab, sk, locRETgainCcurve, locRETtransCcurve, 1, 4, 0.8f, minCD, maxCD, mini, maxi, Tmean, Tsigma, Tmin, Tmax,
-                                          locccmasretiCurve, lcmasretiutili, locllmasretiCurve, llmasretiutili, lochhmasretiCurve, lhmasretiutili, llretiMask,
-                                          lmaskretilocalcurve, localmaskretiutili,
-                                          transformed, lp.enaretiMasktmap, lp.enaretiMask);
+                        ImProcFunctions::MSRLocal(sp, fftw, 0, nullptr, bufreti, bufmask, buforig, buforigmas, orig, tmpl->L, orig1,
+                                                  Wd, Hd, Wd, Hd, params->locallab, sk, locRETgainCcurve, locRETtransCcurve, 1, 4, 0.8f, minCD, maxCD, mini, maxi, Tmean, Tsigma, Tmin, Tmax,
+                                                  locccmasretiCurve, lcmasretiutili, locllmasretiCurve, llmasretiutili, lochhmasretiCurve, lhmasretiutili, llretiMask,
+                                                  lmaskretilocalcurve, localmaskretiutili,
+                                                  transformed, lp.enaretiMasktmap, lp.enaretiMask);
                     }
 
                     if (!lp.invret && call == 1) {
@@ -9963,11 +10129,11 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o
                     bool fftw = false;
 
                     if (params->locallab.spots.at(sp).chrrt > 40.f) { //second step active Retinex Chroma
-                        ImProcFunctions::MSRLocal(sp, fftw, 0, nullptr, bufreti, bufmask, buforig, buforigmas, orig, tmpl->L, orig1, 
-                                         Wd, Hd, Wd, Hd, params->locallab, sk, locRETgainCcurve, locRETtransCcurve, 1, 4, 0.8f, minCD, maxCD, mini, maxi, Tmean, Tsigma, Tmin, Tmax,
-                                         locccmasretiCurve, lcmasretiutili, locllmasretiCurve, llmasretiutili, lochhmasretiCurve, lhmasretiutili, llretiMask, 
-                                         lmaskretilocalcurve, localmaskretiutili,
-                                         transformed, lp.enaretiMasktmap, lp.enaretiMask);
+                        ImProcFunctions::MSRLocal(sp, fftw, 0, nullptr, bufreti, bufmask, buforig, buforigmas, orig, tmpl->L, orig1,
+                                                  Wd, Hd, Wd, Hd, params->locallab, sk, locRETgainCcurve, locRETtransCcurve, 1, 4, 0.8f, minCD, maxCD, mini, maxi, Tmean, Tsigma, Tmin, Tmax,
+                                                  locccmasretiCurve, lcmasretiutili, locllmasretiCurve, llmasretiutili, lochhmasretiCurve, lhmasretiutili, llretiMask,
+                                                  lmaskretilocalcurve, localmaskretiutili,
+                                                  transformed, lp.enaretiMasktmap, lp.enaretiMask);
                     }
 
                     if (!lp.invret && call == 2) {
diff --git a/rtengine/procevents.h b/rtengine/procevents.h
index a25ca3aea..94e2f0f14 100644
--- a/rtengine/procevents.h
+++ b/rtengine/procevents.h
@@ -807,6 +807,7 @@ enum ProcEventCode {
     Evlocallabwavmaskbl = 777,
     EvlocallabLLmaskcolshapewav = 778,
     Evlocallabwavmaskcol = 779,
+    EvlocallabcsThreshold = 780,
     NUMOFEVENTS
 };
 
diff --git a/rtengine/procparams.cc b/rtengine/procparams.cc
index f298b372d..538479862 100644
--- a/rtengine/procparams.cc
+++ b/rtengine/procparams.cc
@@ -2658,6 +2658,7 @@ LocallabParams::LocallabSpot::LocallabSpot() :
     fftwlc(false),
     localcontMethod("loc"),
     locwavcurve{(double)FCT_MinMaxCPoints, 0.0, 0.5, 0.35, 0.35, 1., 0.5, 0.35, 0.35},
+    csthreshold(1, 1, 5, 5, false),
     // Contrast by detail levels
     expcbdl(false),
     mult{1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
@@ -2951,6 +2952,7 @@ bool LocallabParams::LocallabSpot::operator ==(const LocallabSpot& other) const
         && fftwlc == other.fftwlc
         && localcontMethod == other.localcontMethod
         && locwavcurve == other.locwavcurve
+        && csthreshold == other.csthreshold
         // Constrast by detail levels
         && expcbdl == other.expcbdl
         && [this, &other]()->bool {
@@ -4230,6 +4232,7 @@ int ProcParams::save(const Glib::ustring& fname, const Glib::ustring& fname2, bo
                 saveToKeyfile(!pedited || pedited->locallab.spots.at(i).fftwlc, "Locallab", "Fftwlc_" + std::to_string(i), spot.fftwlc, keyFile);
                 saveToKeyfile(!pedited || pedited->locallab.spots.at(i).localcontMethod, "Locallab", "localcontMethod_" + std::to_string(i), spot.localcontMethod, keyFile);
                 saveToKeyfile(!pedited || pedited->locallab.spots.at(i).locwavcurve, "Locallab", "LocwavCurve_" + std::to_string(i), spot.locwavcurve, keyFile);
+                saveToKeyfile(!pedited || pedited->locallab.spots.at(i).csthreshold, "Locallab", "CSThreshold_" + std::to_string(i), spot.csthreshold.toVector(), keyFile);
                 // Contrast by detail levels
                 saveToKeyfile(!pedited || pedited->locallab.spots.at(i).expcbdl, "Locallab", "Expcbdl_" + std::to_string(i), spot.expcbdl, keyFile);
 
@@ -5636,6 +5639,19 @@ int ProcParams::load(const Glib::ustring& fname, ParamsEdited* pedited)
                 assignFromKeyfile(keyFile, "Locallab", "Fftwlc_" + std::to_string(i), pedited, spot.fftwlc, spotEdited.fftwlc);
                 assignFromKeyfile(keyFile, "Locallab", "localcontMethod_" + std::to_string(i), pedited, spot.localcontMethod, spotEdited.localcontMethod);
                 assignFromKeyfile(keyFile, "Locallab", "LocwavCurve_" + std::to_string(i), pedited, spot.locwavcurve, spotEdited.locwavcurve);
+                if (keyFile.has_key("Locallab", "CSThreshold_" + std::to_string(i))) {
+
+                    const std::vector<int> thresh = keyFile.get_integer_list("Locallab", "CSThreshold_" + std::to_string(i));
+
+                    if (thresh.size() >= 4) {
+                        spot.csthreshold.setValues(thresh[0], thresh[1], min(thresh[2], 10), min(thresh[3], 10));
+                    }
+
+                    if (pedited) {
+                        spotEdited.csthreshold = true;
+                    }
+                }
+
                 // Contrast by detail levels
                 assignFromKeyfile(keyFile, "Locallab", "Expcbdl_" + std::to_string(i), pedited, spot.expcbdl, spotEdited.expcbdl);
 
diff --git a/rtengine/procparams.h b/rtengine/procparams.h
index addf5e272..ed6a196e8 100644
--- a/rtengine/procparams.h
+++ b/rtengine/procparams.h
@@ -1207,6 +1207,8 @@ struct LocallabParams {
         bool fftwlc;
         Glib::ustring localcontMethod;
         std::vector<double> locwavcurve;
+        Threshold<int> csthreshold;
+
         // Contrast by detail levels
         bool expcbdl;
         double mult[6];
diff --git a/rtengine/refreshmap.cc b/rtengine/refreshmap.cc
index 34a460d3a..c8570ed0a 100644
--- a/rtengine/refreshmap.cc
+++ b/rtengine/refreshmap.cc
@@ -806,7 +806,8 @@ int refreshmap[rtengine::NUMOFEVENTS] = {
     LUMINANCECURVE,          //EvlocallabLLmaskblshapewav
     LUMINANCECURVE,          //Evlocallabwavmaskbl
     LUMINANCECURVE,          //EvlocallabLLmaskcolshapewav
-    LUMINANCECURVE          //Evlocallabwavmaskcol
+    LUMINANCECURVE,          //Evlocallabwavmaskcol
+    LUMINANCECURVE          //EvlocallabcsThreshold
 };
 
 namespace rtengine
diff --git a/rtgui/locallab.cc b/rtgui/locallab.cc
index 055f0c3d0..41f6fa9f8 100644
--- a/rtgui/locallab.cc
+++ b/rtgui/locallab.cc
@@ -424,7 +424,8 @@ Locallab::Locallab():
     showmaskretiMethod(Gtk::manage(new MyComboBoxText())),
     //Local contrast
     localcontMethod(Gtk::manage(new MyComboBoxText())),
-    //CBDL
+    csThreshold(Gtk::manage(new ThresholdAdjuster(M("TP_LOCALLAB_CSTHRESHOLD"), 1, 9, 1, 1, 5, 5, 0, false))),
+   //CBDL
     showmaskcbMethod(Gtk::manage(new MyComboBoxText())),
 
     // ThresholdAdjuster widgets
@@ -1789,6 +1790,7 @@ Locallab::Locallab():
     expcontrast->signal_button_release_event().connect_notify(sigc::bind(sigc::mem_fun(this, &Locallab::foldAllButMe), expcontrast));
     enablecontrastConn = expcontrast->signal_enabled_toggled().connect(sigc::bind(sigc::mem_fun(this, &Locallab::enableToggled), expcontrast));
     fftwlcConn  = fftwlc->signal_toggled().connect(sigc::mem_fun(*this, &Locallab::fftwlcChanged));
+    csThreshold->setAdjusterListener(this);
 
     if (showtooltip) {
         fftwlc->set_tooltip_text(M("TP_LOCALLAB_LC_FFTW_TOOLTIP"));
@@ -1856,7 +1858,8 @@ Locallab::Locallab():
     contrastBox->pack_start(*lcdarkness);
     contrastBox->pack_start(*lclightness);
     contrastBox->pack_start(*LocalcurveEditorwav, Gtk::PACK_SHRINK, 4);
-    contrastBox->pack_start(*levelwav);
+//    contrastBox->pack_start(*levelwav);
+    contrastBox->pack_start(*csThreshold);
     contrastBox->pack_start(*residcont);
     contrastBox->pack_start(*residchro);
     contrastBox->pack_start(*clariFrame);
@@ -3552,6 +3555,7 @@ void Locallab::write(ProcParams* pp, ParamsEdited* pedited)
                     pp->locallab.spots.at(pp->locallab.selspot).sensilc = sensilc->getIntValue();
                     pp->locallab.spots.at(pp->locallab.selspot).fftwlc = fftwlc->get_active();
                     pp->locallab.spots.at(pp->locallab.selspot).locwavcurve = wavshape->getCurve();
+                    pp->locallab.spots.at(pp->locallab.selspot).csthreshold = csThreshold->getValue<int>();
 
                     if (localcontMethod->get_active_row_number() == 0) {
                         pp->locallab.spots.at(pp->locallab.selspot).localcontMethod = "loc";
@@ -3854,6 +3858,7 @@ void Locallab::write(ProcParams* pp, ParamsEdited* pedited)
                         pe->locallab.spots.at(pp->locallab.selspot).fftwlc = pe->locallab.spots.at(pp->locallab.selspot).fftwlc || !fftwlc->get_inconsistent();
                         pe->locallab.spots.at(pp->locallab.selspot).localcontMethod = pe->locallab.spots.at(pp->locallab.selspot).localcontMethod || localcontMethod->get_active_text() != M("GENERAL_UNCHANGED");
                         pe->locallab.spots.at(pp->locallab.selspot).locwavcurve = pe->locallab.spots.at(pp->locallab.selspot).locwavcurve || !wavshape->isUnChanged();
+                        pe->locallab.spots.at(pp->locallab.selspot).csthreshold = pe->locallab.spots.at(pp->locallab.selspot).csthreshold || csThreshold->getEditedState();
                         // Contrast by detail levels
                         pe->locallab.spots.at(pp->locallab.selspot).expcbdl = pe->locallab.spots.at(pp->locallab.selspot).expcbdl || !expcbdl->get_inconsistent();
 
@@ -4152,6 +4157,7 @@ void Locallab::write(ProcParams* pp, ParamsEdited* pedited)
                         pedited->locallab.spots.at(pp->locallab.selspot).sensilc = pedited->locallab.spots.at(pp->locallab.selspot).sensilc || sensilc->getEditedState();
                         pedited->locallab.spots.at(pp->locallab.selspot).fftwlc = pedited->locallab.spots.at(pp->locallab.selspot).fftwlc || !fftwlc->get_inconsistent();
                         pedited->locallab.spots.at(pp->locallab.selspot).locwavcurve = pedited->locallab.spots.at(pp->locallab.selspot).locwavcurve || !wavshape->isUnChanged();
+                        pedited->locallab.spots.at(pp->locallab.selspot).csthreshold = pedited->locallab.spots.at(pp->locallab.selspot).csthreshold || csThreshold->getEditedState();
                         // Contrast by detail levels
                         pedited->locallab.spots.at(pp->locallab.selspot).expcbdl = pedited->locallab.spots.at(pp->locallab.selspot).expcbdl || !expcbdl->get_inconsistent();
 
@@ -6056,6 +6062,7 @@ void Locallab::setDefaults(const ProcParams * defParams, const ParamsEdited * pe
     claricres->setDefault(defSpot->claricres);
     clarisoft->setDefault(defSpot->clarisoft);
     sensilc->setDefault((double)defSpot->sensilc);
+    csThreshold->setDefault<int>(defSpot->csthreshold);
 
     // Contrast by detail levels
     for (int i = 0; i < 6; i++) {
@@ -6234,6 +6241,7 @@ void Locallab::setDefaults(const ProcParams * defParams, const ParamsEdited * pe
         claricres->setDefaultEditedState(Irrelevant);
         clarisoft->setDefaultEditedState(Irrelevant);
         sensilc->setDefaultEditedState(Irrelevant);
+        csThreshold->setDefaultEditedState(Irrelevant);
 
         // Contrast by detail levels
         for (int i = 0; i < 6; i++) {
@@ -6416,6 +6424,7 @@ void Locallab::setDefaults(const ProcParams * defParams, const ParamsEdited * pe
         claricres->setDefaultEditedState(defSpotState->claricres ? Edited : UnEdited);
         clarisoft->setDefaultEditedState(defSpotState->clarisoft ? Edited : UnEdited);
         sensilc->setDefaultEditedState(defSpotState->sensilc ? Edited : UnEdited);
+        csThreshold->setDefaultEditedState(defSpotState->csthreshold ? Edited : UnEdited);
 
         // Contrast by detail levels
         for (int i = 0; i < 6; i++) {
@@ -6475,7 +6484,12 @@ void Locallab::adjusterChanged(ThresholdAdjuster* a, int newBottomLeft, int newT
 
 void Locallab::adjusterChanged2(ThresholdAdjuster* a, int newBottomL, int newTopL, int newBottomR, int newTopR)
 {
-    // Not used
+    if (getEnabled() && expcontrast->getEnabled()) {
+        if (listener) {
+            listener->panelChanged(EvlocallabcsThreshold, csThreshold->getHistoryString());
+        }
+    }
+
 }
 
 void Locallab::adjusterChanged(ThresholdAdjuster* a, int newBottom, int newTop)
@@ -7696,6 +7710,7 @@ void Locallab::setBatchMode(bool batchMode)
     clarilres->showEditedCB();
     claricres->showEditedCB();
     clarisoft->showEditedCB();
+    csThreshold->showEditedCB();
     sensilc->showEditedCB();
 
     // Contrast by detail levels
@@ -8368,6 +8383,7 @@ void Locallab::updateLocallabGUI(const rtengine::procparams::ProcParams* pp, con
         clarisoft->setValue(pp->locallab.spots.at(index).clarisoft);
         sensilc->setValue(pp->locallab.spots.at(index).sensilc);
         fftwlc->set_active(pp->locallab.spots.at(index).fftwlc);
+        csThreshold->setValue<int>(pp->locallab.spots.at(index).csthreshold);
 
         if (pp->locallab.spots.at(index).localcontMethod == "loc") {
             localcontMethod->set_active(0);
@@ -8720,6 +8736,7 @@ void Locallab::updateLocallabGUI(const rtengine::procparams::ProcParams* pp, con
                 sensilc->setEditedState(spotState->sensilc ? Edited : UnEdited);
                 fftwlc->set_inconsistent(multiImage && !spotState->fftwlc);
                 wavshape->setUnChanged(!spotState->locwavcurve);
+                csThreshold->setEditedState(spotState->csthreshold ? Edited : UnEdited);
 
                 if (!spotState->retinexMethod) {
                     localcontMethod->set_active_text(M("GENERAL_UNCHANGED"));
diff --git a/rtgui/locallab.h b/rtgui/locallab.h
index 18c7141b4..cf7c009af 100644
--- a/rtgui/locallab.h
+++ b/rtgui/locallab.h
@@ -429,6 +429,7 @@ private:
     //local contrast
     MyComboBoxText* const localcontMethod;
     sigc::connection localcontMethodConn;
+    ThresholdAdjuster* const csThreshold;
     //CBDL
     MyComboBoxText* const showmaskcbMethod;
     sigc::connection showmaskcbMethodConn;
diff --git a/rtgui/paramsedited.cc b/rtgui/paramsedited.cc
index 1528b2d89..f5d9ebd35 100644
--- a/rtgui/paramsedited.cc
+++ b/rtgui/paramsedited.cc
@@ -1198,6 +1198,7 @@ void ParamsEdited::initFrom(const std::vector<rtengine::procparams::ProcParams>&
                 locallab.spots.at(j).fftwlc = locallab.spots.at(j).fftwlc && pSpot.fftwlc == otherSpot.fftwlc;
                 locallab.spots.at(j).localcontMethod = locallab.spots.at(j).localcontMethod && pSpot.localcontMethod == otherSpot.localcontMethod;
                 locallab.spots.at(j).locwavcurve = locallab.spots.at(j).locwavcurve && pSpot.locwavcurve == otherSpot.locwavcurve;
+                locallab.spots.at(j).csthreshold = locallab.spots.at(j).csthreshold && pSpot.csthreshold == otherSpot.csthreshold;
                 // Contrast by detail levels
                 locallab.spots.at(j).expcbdl = locallab.spots.at(j).expcbdl && pSpot.expcbdl == otherSpot.expcbdl;
 
@@ -3664,6 +3665,10 @@ void ParamsEdited::combine(rtengine::procparams::ProcParams& toEdit, const rteng
             toEdit.locallab.spots.at(i).locwavcurve = mods.locallab.spots.at(i).locwavcurve;
         }
 
+        if (locallab.spots.at(i).csthreshold) {
+            toEdit.locallab.spots.at(i).csthreshold = mods.locallab.spots.at(i).csthreshold;
+        }
+
         // Contrast by detail levels
         if (locallab.spots.at(i).expcbdl) {
             toEdit.locallab.spots.at(i).expcbdl = mods.locallab.spots.at(i).expcbdl;
@@ -5021,6 +5026,7 @@ LocallabParamsEdited::LocallabSpotEdited::LocallabSpotEdited(bool v) :
     fftwlc(v),
     localcontMethod(v),
     locwavcurve(v),
+    csthreshold(v),
     // Contrast by detail levels
     expcbdl(v),
     mult{v, v, v, v, v, v},
@@ -5311,6 +5317,7 @@ void LocallabParamsEdited::LocallabSpotEdited::set(bool v)
     fftwlc = v;
     localcontMethod = v;
     locwavcurve = v;
+    csthreshold = v;
     // Contrast by detail levels
     expcbdl = v;
 
diff --git a/rtgui/paramsedited.h b/rtgui/paramsedited.h
index ebbd1213f..4229ea49c 100644
--- a/rtgui/paramsedited.h
+++ b/rtgui/paramsedited.h
@@ -617,6 +617,7 @@ public:
         bool fftwlc;
         bool localcontMethod;
         bool locwavcurve;
+        bool csthreshold;
         // Contrast by detail levels
         bool expcbdl;
         bool mult[6];