/*
*  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 <https://www.gnu.org/licenses/>.

    *   adaptation to RawTherapee
    *   2015 Jacques Desmis <jdesmis@gmail.com>
    *   2015 Ingo Weyrich <heckflosse67@gmx.de>

    * D. J. Jobson, Z. Rahman, and G. A. Woodell. A multi-scale
    * Retinex for bridging the gap between color images and the
    * human observation of scenes. IEEE Transactions on Image Processing,
    * 1997, 6(7): 965-976

    * Fan Guo Zixing Cai Bin Xie Jin Tang
    * School of Information Science and Engineering, Central South University Changsha, China

    * Weixing Wang and Lian Xu
    * College of Physics and Information Engineering, Fuzhou University, Fuzhou, China

    * inspired from 2003 Fabien Pelisson <Fabien.Pelisson@inrialpes.fr>
    * some ideas taken (use of mask) Russell Cottrell - The Retinex .8bf Plugin

*/

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>

#include "gauss.h"
#include "improcfun.h"
#include "median.h"
#include "opthelper.h"
#include "procparams.h"
#include "rawimagesource.h"
#include "rtengine.h"

#define BENCHMARK
#include "StopWatch.h"
#include "guidedfilter.h"

#define clipretinex( val, minv, maxv )    (( val = (val < minv ? minv : val ) ) > maxv ? maxv : val )
#define CLIPLOC(x) LIM(x,0.f,32767.f)
#define CLIPC(a) LIM(a, -42000.f, 42000.f)  // limit a and b  to 130 probably enough ?
#define CLIPMAX(x) LIM(x,0.f,500000.f)

namespace
{

void calcGammaLut(double gamma, double ts, LUTf &gammaLut)
{
    double pwr = 1.0 / gamma;
    double gamm = gamma;
    const double gamm2 = gamma;
    rtengine::GammaValues g_a;

    if (gamm2 < 1.0) {
        std::swap(pwr, gamm);
    }

    rtengine::Color::calcGamma(pwr, ts, 0, g_a); // call to calcGamma with selected gamma and slope

    const double start = gamm2 < 1. ? g_a[2] : g_a[3];
    const double add = g_a[4];
    const double mul = 1.0 + g_a[4];

    if (gamm2 < 1.) {
        #pragma omp parallel for schedule(dynamic, 1024)

        for (int i = 0; i < 65536; i++) {
            const double x = rtengine::Color::igammareti(i / 65535.0, gamm, start, ts, mul, add);
            gammaLut[i] = 0.5 * rtengine::CLIP(x * 65535.0);  // CLIP avoid in some case extra values
        }
    } else {
        #pragma omp parallel for schedule(dynamic, 1024)

        for (int i = 0; i < 65536; i++) {
            const double x = rtengine::Color::gammareti(i / 65535.0, gamm, start, ts, mul, add);
            gammaLut[i] = 0.5 * rtengine::CLIP(x * 65535.0);  // CLIP avoid in some case extra values
        }
    }
}

void retinex_scales(float* scales, int nscales, int mode, int s, float high)
{
    if (s < 3) {
        s = 3;    //to avoid crash in MSRlocal if nei small
    }

    if (nscales == 1) {
        scales[0] = (float)s / 2.f;
    } else if (nscales == 2) {
        scales[1] = (float) s / 2.f;
        scales[0] = (float) s;
    } else {
        float size_step = (float) s / (float) nscales;

        if (mode == 0) {
            for (int i = 0; i < nscales; ++i) {
                scales[nscales - i - 1] = 2.0f + (float)i * size_step;
            }
        } else if (mode == 1) {
            size_step = (float)log(s - 2.0f) / (float) nscales;

            for (int i = 0; i < nscales; ++i) {
                scales[nscales - i - 1] = 2.0f + (float)pow(10.f, (i * size_step) / log(10.f));
            }
        } else if (mode == 2) {
            size_step = (float) log(s - 2.0f) / (float) nscales;

            for (int i = 0; i < nscales; ++i) {
                scales[i] = s - (float)pow(10.f, (i * size_step) / log(10.f));
            }
        } else if (mode == 3) {
            size_step = (float) log(s - 2.0f) / (float) nscales;

            for (int i = 0; i < nscales; ++i) {
                scales[i] = high * s - (float)pow(10.f, (i * size_step) / log(10.f));
            }
        }
    }
}

void mean_stddv2(float **dst, float &mean, float &stddv, int W_L, int H_L, float &maxtr, float &mintr)
{
    // summation using double precision to avoid too large summation error for large pictures
    double vsquared = 0.f;
    double sum = 0.f;
    maxtr = -999999.f;
    mintr = 999999.f;
#ifdef _OPENMP
    #pragma omp parallel
#endif
    {
        float lmax = -999999.f, lmin = 999999.f;
#ifdef _OPENMP
        #pragma omp for reduction(+:sum,vsquared) nowait // this leads to differences, but parallel summation is more accurate
#endif

        for (int i = 0; i < H_L; i++)
            for (int j = 0; j < W_L; j++) {
                sum += dst[i][j];
                vsquared += (dst[i][j] * dst[i][j]);

                lmax = dst[i][j] > lmax ? dst[i][j] : lmax;
                lmin = dst[i][j] < lmin ? dst[i][j] : lmin;
            }

#ifdef _OPENMP
        #pragma omp critical
#endif
        {
            maxtr = maxtr > lmax ? maxtr : lmax;
            mintr = mintr < lmin ? mintr : lmin;
        }

    }
    mean = sum / (double)(W_L * H_L);
    vsquared /= (double) W_L * H_L;
    stddv = (vsquared - (mean * mean));
    stddv = (float)sqrt(stddv);
}

}


namespace rtengine
{

extern const Settings* settings;

void RawImageSource::MSR(float** luminance, float** originalLuminance, float **exLuminance,  LUTf & mapcurve, bool &mapcontlutili, int width, int height, const RetinexParams &deh, const RetinextransmissionCurve & dehatransmissionCurve, const RetinexgaintransmissionCurve & dehagaintransmissionCurve, float &minCD, float &maxCD, float &mini, float &maxi, float &Tmean, float &Tsigma, float &Tmin, float &Tmax)
{

    if (deh.enabled) {//enabled
        float maxtr, mintr;
        constexpr float eps = 2.f;
        bool useHsl = deh.retinexcolorspace == "HSLLOG";
        bool useHslLin = deh.retinexcolorspace == "HSLLIN";
        float offse = (float) deh.offs; //def = 0  not use
        int iter = deh.iter;
        int gradient = deh.scal;
        int scal = 3;//disabled scal
        int nei = (int)(2.8f * deh.neigh);  //def = 220
        float vart = (float)deh.vart / 100.f;//variance
        float gradvart = (float)deh.grad;
        float gradstr = (float)deh.grads;
        float strength = (float) deh.str / 100.f; // Blend with original L channel data
        float limD = (float) deh.limd;
        limD = pow(limD, 1.7f);  //about 2500 enough
        limD *= useHslLin ? 10.f : 1.f;
        float ilimD = 1.f / limD;
        float hig = ((float) deh.highl) / 100.f;
        scal = deh.skal;

        int H_L = height;
        int W_L = width;

        float *tran[H_L] ALIGNED16;
        float *tranBuffer = nullptr;

        constexpr float elogt = 2.71828f;
        bool lhutili = false;

        FlatCurve* shcurve = new FlatCurve(deh.lhcurve);  //curve L=f(H)

        if (!shcurve || shcurve->isIdentity()) {
            if (shcurve) {
                delete shcurve;
                shcurve = nullptr;
            }
        } else {
            lhutili = true;
        }


        bool higplus = false ;
        int moderetinex = 2; // default to 2 ( deh.retinexMethod == "high" )

        if (deh.retinexMethod == "highliplus") {
            higplus = true;
            moderetinex = 3;
        } else if (deh.retinexMethod == "uni") {
            moderetinex = 0;
        } else if (deh.retinexMethod == "low") {
            moderetinex = 1;
        } else { /*if (deh.retinexMethod == "highli") */
            moderetinex = 3;
        }

        constexpr float aahi = 49.f / 99.f; ////reduce sensibility 50%
        constexpr float bbhi = 1.f - aahi;

        for (int it = 1; it < iter + 1; it++) { //iter nb max of iterations
            float high = bbhi + aahi * (float) deh.highl;

            float grad = 1.f;
            float sc = scal;

            if (gradient == 0) {
                grad = 1.f;
                sc = 3.f;
            } else if (gradient == 1) {
                grad = 0.25f * it + 0.75f;
                sc = -0.5f * it + 4.5f;
            } else if (gradient == 2) {
                grad = 0.5f * it + 0.5f;
                sc = -0.75f * it + 5.75f;
            } else if (gradient == 3) {
                grad = 0.666f * it + 0.333f;
                sc = -0.75f * it + 5.75f;
            } else if (gradient == 4) {
                grad = 0.8f * it + 0.2f;
                sc = -0.75f * it + 5.75f;
            } else if (gradient == 5) {
                if (moderetinex != 3) {
                    grad = 2.5f * it - 1.5f;
                } else {
                    float aa = (11.f * high - 1.f) / 4.f;
                    float bb = 1.f - aa;
                    grad = aa * it + bb;
                }

                sc = -0.75f * it + 5.75f;
            } else if (gradient == 6) {
                if (moderetinex != 3) {
                    grad = 5.f * it - 4.f;
                } else {
                    float aa = (21.f * high - 1.f) / 4.f;
                    float bb = 1.f - aa;
                    grad = aa * it + bb;
                }

                sc = -0.75f * it + 5.75f;
            }

            else if (gradient == -1) {
                grad = -0.125f * it + 1.125f;
                sc = 3.f;
            }

            if (iter == 1) {
                sc = scal;
            } else {
                //adjust sc in function of choice of scale by user if iterations
                if (scal < 3) {
                    sc -= 1;

                    if (sc < 1.f) { //avoid 0
                        sc = 1.f;
                    }
                }

                if (scal > 4) {
                    sc += 1;
                }
            }

            float varx = vart;
            float limdx = limD;
            float ilimdx = ilimD;

            if (gradvart != 0) {
                if (gradvart == 1) {
                    varx = vart * (-0.125f * it + 1.125f);
                    limdx = limD * (-0.125f * it + 1.125f);
                    ilimdx = 1.f / limdx;
                } else if (gradvart == 2) {
                    varx = vart * (-0.2f * it + 1.2f);
                    limdx = limD * (-0.2f * it + 1.2f);
                    ilimdx = 1.f / limdx;
                } else if (gradvart == -1) {
                    varx = vart * (0.125f * it + 0.875f);
                    limdx = limD * (0.125f * it + 0.875f);
                    ilimdx = 1.f / limdx;
                } else if (gradvart == -2) {
                    varx = vart * (0.4f * it + 0.6f);
                    limdx = limD * (0.4f * it + 0.6f);
                    ilimdx = 1.f / limdx;
                }
            }

            scal = round(sc);
            float ks = 1.f;

            if (gradstr != 0) {
                if (gradstr == 1) {
                    if (it <= 3) {
                        ks = -0.3f * it + 1.6f;
                    } else {
                        ks = 0.5f;
                    }
                } else if (gradstr == 2) {
                    if (it <= 3) {
                        ks = -0.6f * it + 2.2f;
                    } else {
                        ks = 0.3f;
                    }
                } else if (gradstr == -1) {
                    if (it <= 3) {
                        ks = 0.2f * it + 0.6f;
                    } else {
                        ks = 1.2f;
                    }
                } else if (gradstr == -2) {
                    if (it <= 3) {
                        ks = 0.4f * it + 0.2f;
                    } else {
                        ks = 1.5f;
                    }
                }
            }

            float strengthx = ks * strength;

            constexpr auto maxRetinexScales = 8;
            float RetinexScales[maxRetinexScales];

            retinex_scales(RetinexScales, scal, moderetinex, nei / grad, high);

            float *src[H_L] ALIGNED16;
            float *srcBuffer = new float[H_L * W_L];

            for (int i = 0; i < H_L; i++) {
                src[i] = &srcBuffer[i * W_L];
            }

            int h_th = 0, s_th = 0;

            int shHighlights = deh.highlights;
            int shShadows = deh.shadows;

            int mapmet = 0;

            if (deh.mapMethod == "map") {
                mapmet = 2;
            } else if (deh.mapMethod == "mapT") {
                mapmet = 3;
            } else if (deh.mapMethod == "gaus") {
                mapmet = 4;
            }

            const double shradius = mapmet == 4 ? (double) deh.radius : 40.;

            int viewmet = 0;

            if (deh.viewMethod == "mask") {
                viewmet = 1;
            } else if (deh.viewMethod == "tran") {
                viewmet = 2;
            } else if (deh.viewMethod == "tran2") {
                viewmet = 3;
            } else if (deh.viewMethod == "unsharp") {
                viewmet = 4;
            }

#ifdef _OPENMP
            #pragma omp parallel for
#endif

            for (int i = 0; i < H_L; i++)
                for (int j = 0; j < W_L; j++) {
                    src[i][j] = luminance[i][j] + eps;
                    luminance[i][j] = 0.f;
                }

            float *out[H_L] ALIGNED16;
            float *outBuffer = new float[H_L * W_L];

            for (int i = 0; i < H_L; i++) {
                out[i] = &outBuffer[i * W_L];
            }

            if (viewmet == 3  || viewmet == 2) {
                tranBuffer = new float[H_L * W_L];

                for (int i = 0; i < H_L; i++) {
                    tran[i] = &tranBuffer[i * W_L];
                }
            }

            const float logBetaGain = xlogf(16384.f);
            float pond = logBetaGain / (float) scal;

            if (!useHslLin) {
                pond /= log(elogt);
            }

            auto shmap = ((mapmet == 2 || mapmet == 3 || mapmet == 4) && it == 1) ? new SHMap(W_L, H_L) : nullptr;

            float *buffer = new float[W_L * H_L];;

            for (int scale = scal - 1; scale >= 0; scale--) {
#ifdef _OPENMP
                #pragma omp parallel
#endif
                {
                    if (scale == scal - 1)
                    {
                        gaussianBlur(src, out, W_L, H_L, RetinexScales[scale], buffer);
                    } else   // reuse result of last iteration
                    {
                        // out was modified in last iteration => restore it
                        if ((((mapmet == 2 && scale > 1) || mapmet == 3 || mapmet == 4) || (mapmet > 0 && mapcontlutili)) && it == 1) {
#ifdef _OPENMP
                            #pragma omp for
#endif

                            for (int i = 0; i < H_L; i++) {
                                for (int j = 0; j < W_L; j++) {
                                    out[i][j] = buffer[i * W_L + j];
                                }
                            }
                        }

                        gaussianBlur(out, out, W_L, H_L, sqrtf(SQR(RetinexScales[scale]) - SQR(RetinexScales[scale + 1])), buffer);
                    }

                    if ((((mapmet == 2 && scale > 2) || mapmet == 3 || mapmet == 4) || (mapmet > 0 && mapcontlutili)) && it == 1 && scale > 0)
                    {
                        // out will be modified => store it for use in next iteration. We even don't need a new buffer because 'buffer' is free after gaussianBlur :)
#ifdef _OPENMP
                        #pragma omp for
#endif

                        for (int i = 0; i < H_L; i++) {
                            for (int j = 0; j < W_L; j++) {
                                buffer[i * W_L + j] = out[i][j];
                            }
                        }
                    }
                }

                if (((mapmet == 2 && scale > 2) || mapmet == 3 || mapmet == 4) && it == 1) {
                    shmap->updateL(out, shradius, true, 1);
                    h_th = shmap->max_f - deh.htonalwidth * (shmap->max_f - shmap->avg) / 100;
                    s_th = deh.stonalwidth * (shmap->avg - shmap->min_f) / 100;
                }

#ifdef __SSE2__
                vfloat pondv = F2V(pond);
                vfloat limMinv = F2V(ilimdx);
                vfloat limMaxv = F2V(limdx);

#endif

                if (mapmet > 0 && mapcontlutili && it == 1) {
                    // TODO: When rgbcurvespeedup branch is merged into master we can simplify the code by
                    // 1) in rawimagesource.retinexPrepareCurves() insert
                    //    mapcurve *= 0.5f;
                    //    after
                    //    CurveFactory::mapcurve (mapcontlutili, retinexParams.mapcurve, mapcurve, 1, lhist16RETI, histLRETI);
                    // 2) remove the division by 2.f from the code 7 lines below this line
#ifdef _OPENMP
                    #pragma omp parallel for
#endif

                    for (int i = 0; i < H_L; i++) {
                        for (int j = 0; j < W_L; j++) {
                            out[i][j] = mapcurve[2.f * out[i][j]] / 2.f;
                        }
                    }

                }

                if (((mapmet == 2 && scale > 2) || mapmet == 3 || mapmet == 4) && it == 1) {
                    float hWeight = (100.f - shHighlights) / 100.f;
                    float sWeight = (100.f - shShadows) / 100.f;
#ifdef _OPENMP
                    #pragma omp parallel for schedule(dynamic,16)
#endif

                    for (int i = 0; i < H_L; i++) {
                        for (int j = 0; j < W_L; j++) {
                            float mapval = 1.f + shmap->map[i][j];
                            float factor = 1.f;

                            if (mapval > h_th) {
                                factor = (h_th + hWeight * (mapval - h_th)) / mapval;
                            } else if (mapval < s_th) {
                                factor = (s_th - sWeight * (s_th - mapval)) / mapval;
                            }

                            out[i][j] *= factor;

                        }
                    }
                }

#ifdef _OPENMP
                #pragma omp parallel for
#endif

                for (int i = 0; i < H_L; i++) {
                    int j = 0;

#ifdef __SSE2__

                    if (useHslLin) {
                        for (; j < W_L - 3; j += 4) {
                            _mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * (vclampf(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv)));
                        }
                    } else {
                        for (; j < W_L - 3; j += 4) {
                            _mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv *  xlogf(vclampf(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv)));
                        }
                    }

#endif

                    if (useHslLin) {
                        for (; j < W_L; j++) {
                            luminance[i][j] +=  pond * (LIM(src[i][j] / out[i][j], ilimdx, limdx));
                        }
                    } else {
                        for (; j < W_L; j++) {
                            luminance[i][j] +=  pond * xlogf(LIM(src[i][j] / out[i][j], ilimdx, limdx));   //  /logt ?
                        }
                    }
                }
            }

            if (mapmet > 1) {
                if (shmap) {
                    delete shmap;
                }
            }

            shmap = nullptr;

            delete [] buffer;
            delete [] srcBuffer;

            float mean = 0.f;
            float stddv = 0.f;
            // I call mean_stddv2 instead of mean_stddv ==> logBetaGain

            mean_stddv2(luminance, mean, stddv, W_L, H_L, maxtr, mintr);
            //printf("mean=%f std=%f delta=%f maxtr=%f mintr=%f\n", mean, stddv, delta, maxtr, mintr);

            //mean_stddv( luminance, mean, stddv, W_L, H_L, logBetaGain, maxtr, mintr);
            if (dehatransmissionCurve && mean != 0.f && stddv != 0.f) { //if curve
                float asig = 0.166666f / stddv;
                float bsig = 0.5f - asig * mean;
                float amax = 0.333333f / (maxtr - mean - stddv);
                float bmax = 1.f - amax * maxtr;
                float amin = 0.333333f / (mean - stddv - mintr);
                float bmin = -amin * mintr;

                asig *= 500.f;
                bsig *= 500.f;
                amax *= 500.f;
                bmax *= 500.f;
                amin *= 500.f;
                bmin *= 500.f;

#ifdef _OPENMP
                #pragma omp parallel
#endif
                {
                    float absciss;
#ifdef _OPENMP
                    #pragma omp for schedule(dynamic,16)
#endif

                    for (int i = 0; i < H_L; i++)
                        for (int j = 0; j < W_L; j++) { //for mintr to maxtr evalate absciss in function of original transmission
                            if (LIKELY(fabsf(luminance[i][j] - mean) < stddv)) {
                                absciss = asig * luminance[i][j] + bsig;
                            } else if (luminance[i][j] >= mean) {
                                absciss = amax * luminance[i][j] + bmax;
                            } else { /*if(luminance[i][j] <= mean - stddv)*/
                                absciss = amin * luminance[i][j] + bmin;
                            }

                            //TODO : move multiplication by 4.f and subtraction of 1.f inside the curve
                            luminance[i][j] *= (-1.f + 4.f * dehatransmissionCurve[absciss]); //new transmission

                            if (viewmet == 3 || viewmet == 2) {
                                tran[i][j] = luminance[i][j];
                            }
                        }
                }

                // median filter on transmission  ==> reduce artifacts
                if (deh.medianmap && it == 1) { //only one time
                    int wid = W_L;
                    int hei = H_L;
                    float *tmL[hei] ALIGNED16;
                    float *tmLBuffer = new float[wid * hei];
                    int borderL = 1;

                    for (int i = 0; i < hei; i++) {
                        tmL[i] = &tmLBuffer[i * wid];
                    }

#ifdef _OPENMP
                    #pragma omp parallel for
#endif

                    for (int i = borderL; i < hei - borderL; i++) {
                        for (int j = borderL; j < wid - borderL; j++) {
                            tmL[i][j] = median(luminance[i][j], luminance[i - 1][j], luminance[i + 1][j], luminance[i][j + 1], luminance[i][j - 1], luminance[i - 1][j - 1], luminance[i - 1][j + 1], luminance[i + 1][j - 1], luminance[i + 1][j + 1]);  //3x3
                        }
                    }

#ifdef _OPENMP
                    #pragma omp parallel for
#endif

                    for (int i = borderL; i < hei - borderL; i++) {
                        for (int j = borderL; j < wid - borderL; j++) {
                            luminance[i][j] = tmL[i][j];
                        }
                    }

                    delete [] tmLBuffer;

                }

                // I call mean_stddv2 instead of mean_stddv ==> logBetaGain
                //mean_stddv( luminance, mean, stddv, W_L, H_L, 1.f, maxtr, mintr);
                mean_stddv2(luminance, mean, stddv, W_L, H_L, maxtr, mintr);

            }

            float epsil = 0.1f;

            mini = mean - varx * stddv;

            if (mini < mintr) {
                mini = mintr + epsil;
            }

            maxi = mean + varx * stddv;

            if (maxi > maxtr) {
                maxi = maxtr - epsil;
            }

            float delta = maxi - mini;
            //printf("maxi=%f mini=%f mean=%f std=%f delta=%f maxtr=%f mintr=%f\n", maxi, mini, mean, stddv, delta, maxtr, mintr);

            if (!delta) {
                delta = 1.0f;
            }

            float cdfactor = 32768.f / delta;
            maxCD = -9999999.f;
            minCD = 9999999.f;
            // coeff for auto    "transmission" with 2 sigma #95% data
            float aza = 16300.f / (2.f * stddv);
            float azb = -aza * (mean - 2.f * stddv);
            float bza = 16300.f / (2.f * stddv);
            float bzb = 16300.f - bza * (mean);

//prepare work for curve gain
#ifdef _OPENMP
            #pragma omp parallel for
#endif

            for (int i = 0; i < H_L; i++) {
                for (int j = 0; j < W_L; j++) {
                    luminance[i][j] = luminance[i][j] - mini;
                }
            }

            mean = 0.f;
            stddv = 0.f;
            // I call mean_stddv2 instead of mean_stddv ==> logBetaGain

            mean_stddv2(luminance, mean, stddv, W_L, H_L, maxtr, mintr);
            float asig = 0.f, bsig = 0.f, amax = 0.f, bmax = 0.f, amin = 0.f, bmin = 0.f;

            if (dehagaintransmissionCurve && mean != 0.f && stddv != 0.f) { //if curve
                asig = 0.166666f / stddv;
                bsig = 0.5f - asig * mean;
                amax = 0.333333f / (maxtr - mean - stddv);
                bmax = 1.f - amax * maxtr;
                amin = 0.333333f / (mean - stddv - mintr);
                bmin = -amin * mintr;

                asig *= 500.f;
                bsig *= 500.f;
                amax *= 500.f;
                bmax *= 500.f;
                amin *= 500.f;
                bmin *= 500.f;
            }

#ifdef _OPENMP
            #pragma omp parallel
#endif
            {
                float cdmax = -999999.f, cdmin = 999999.f;
#ifdef _OPENMP
                #pragma omp for schedule(dynamic,16) nowait
#endif

                for (int i = 0; i < H_L; i ++)
                    for (int j = 0; j < W_L; j++) {
                        float gan;

                        if (dehagaintransmissionCurve && mean != 0.f && stddv != 0.f) {
                            float absciss;

                            if (LIKELY(fabsf(luminance[i][j] - mean) < stddv)) {
                                absciss = asig * luminance[i][j] + bsig;
                            } else if (luminance[i][j] >= mean) {
                                absciss = amax * luminance[i][j] + bmax;
                            } else { /*if(luminance[i][j] <= mean - stddv)*/
                                absciss = amin * luminance[i][j] + bmin;
                            }


                            //    float cd = cdfactor * ( luminance[i][j]  - mini ) + offse;
                            // TODO : move multiplication by 2.f inside the curve
                            gan = 2.f * (dehagaintransmissionCurve[absciss]); //new gain function transmission
                        } else {
                            gan = 0.5f;
                        }

                        float cd = gan * cdfactor * (luminance[i][j]) + offse;

                        cdmax = cd > cdmax ? cd : cdmax;
                        cdmin = cd < cdmin ? cd : cdmin;

                        float str = strengthx;

                        if (lhutili && it == 1) { // S=f(H)
                            {
                                float HH = exLuminance[i][j];
                                float valparam;

                                if (useHsl || useHslLin) {
                                    valparam = float ((shcurve->getVal(HH) - 0.5f));
                                } else {
                                    valparam = float ((shcurve->getVal(Color::huelab_to_huehsv2(HH)) - 0.5f));
                                }

                                str *= (1.f + 2.f * valparam);
                            }
                        }

                        if (higplus && exLuminance[i][j] > 65535.f * hig) {
                            str *= hig;
                        }

                        if (viewmet == 0) {
                            luminance[i][j] = intp(str, clipretinex(cd, 0.f, 32768.f), originalLuminance[i][j]);
                        } else if (viewmet == 1) {
                            luminance[i][j] = out[i][j];
                        } else if (viewmet == 4) {
                            luminance[i][j] = originalLuminance[i][j] + str * (originalLuminance[i][j] - out[i][j]);//unsharp
                        } else if (viewmet == 2) {
                            if (tran[i][j] <= mean) {
                                luminance[i][j] = azb + aza * tran[i][j];    //auto values
                            } else {
                                luminance[i][j] = bzb + bza * tran[i][j];
                            }
                        } else { /*if(viewmet == 3) */
                            luminance[i][j] = 1000.f + tran[i][j] * 700.f;    //arbitrary values to help display log values which are between -20 to + 30 - usage values -4 + 5
                        }

                    }

#ifdef _OPENMP
                #pragma omp critical
#endif
                {
                    maxCD = maxCD > cdmax ? maxCD : cdmax;
                    minCD = minCD < cdmin ? minCD : cdmin;
                }

            }

            delete [] outBuffer;
            outBuffer = nullptr;
            //printf("cdmin=%f cdmax=%f\n",minCD, maxCD);
            Tmean = mean;
            Tsigma = stddv;
            Tmin = mintr;
            Tmax = maxtr;

            if (shcurve) {
                delete shcurve;
                shcurve = nullptr;
            }
        }

        if (tranBuffer) {
            delete [] tranBuffer;
        }

    }
}


void ImProcFunctions::maskforretinex(int sp, int before, float ** luminance, float ** out, int W_L, int H_L, int skip,
                                     const LocCCmaskCurve & locccmasretiCurve, bool &lcmasretiutili, const  LocLLmaskCurve & locllmasretiCurve, bool &llmasretiutili, const  LocHHmaskCurve & lochhmasretiCurve, bool & lhmasretiutili,
                                     int llretiMask, bool retiMasktmap, bool retiMask, float rad, float lap, bool pde, float gamm, float slop, float chro, float blend,
                                     LUTf & lmaskretilocalcurve, bool & localmaskretiutili,
                                     LabImage * bufreti, LabImage * bufmask, LabImage * buforig, LabImage * buforigmas, bool multiThread,
                                     bool delt, const float hueref, const float chromaref, const float lumaref,
                                     float maxdE, float mindE, float maxdElim,  float mindElim, float iterat, float limscope, int scope, float balance)
{
    array2D<float> loctemp(W_L, H_L);
    array2D<float> ble(W_L, H_L);
    array2D<float> guid(W_L, H_L);
    std::unique_ptr<LabImage> bufmaskblurreti;
    bufmaskblurreti.reset(new LabImage(W_L, H_L));
    std::unique_ptr<LabImage> bufmaskorigreti;
    bufmaskorigreti.reset(new LabImage(W_L, H_L));
    std::unique_ptr<LabImage> bufprov;
    bufprov.reset(new LabImage(W_L, H_L));

#ifdef _OPENMP
    #pragma omp parallel for schedule(dynamic,16)
#endif

    for (int y = 0; y < H_L; y++) {
        for (int x = 0; x < W_L; x++) {
            if (before == 1 && retiMasktmap) {
                loctemp[y][x] = LIM(luminance[y][x], 0.f, 32768.f);
            } else if (before == 0 && retiMasktmap) {
                loctemp[y][x] = out[y][x];
            } else {
                loctemp[y][x] = bufreti->L[y][x];
            }
        }
    }


    float fab = 4000.f;//value must be good in most cases

#ifdef _OPENMP
    #pragma omp parallel for schedule(dynamic,16)
#endif

    for (int ir = 0; ir < H_L; ir++) {
        for (int jr = 0; jr < W_L; jr++) {
            float kmaskLexp = 0;
            float kmaskCH = 0;

            if (locllmasretiCurve  && llmasretiutili) {
                float ligh = loctemp[ir][jr] / 32768.f;
                kmaskLexp = 32768.f * LIM01(1.f - locllmasretiCurve[500.f * ligh]);

            }


            if (locllmasretiCurve  && llmasretiutili && retiMasktmap) {
            }

            if (llretiMask != 4) {
                if (locccmasretiCurve && lcmasretiutili) {
                    float chromask = 0.0001f + sqrt(SQR((bufreti->a[ir][jr]) / fab) + SQR((bufreti->b[ir][jr]) / fab));
                    kmaskCH = LIM01(1.f - locccmasretiCurve[500.f *  chromask]);
                }
            }

            if (lochhmasretiCurve && lhmasretiutili) {
                float huema = xatan2f(bufreti->b[ir][jr], bufreti->a[ir][jr]);
                float h = Color::huelab_to_huehsv2(huema);
                h += 1.f / 6.f;

                if (h > 1.f) {
                    h -= 1.f;
                }

                float valHH = LIM01(1.f - lochhmasretiCurve[500.f *  h]);

                if (llretiMask != 4) {
                    kmaskCH += valHH;
                }

                kmaskLexp += 32768.f * valHH;
            }

            bufmaskblurreti->L[ir][jr] = kmaskLexp;
            bufmaskblurreti->a[ir][jr] = kmaskCH;
            bufmaskblurreti->b[ir][jr] = kmaskCH;
            ble[ir][jr] = bufmaskblurreti->L[ir][jr] / 32768.f;
            guid[ir][jr] = bufreti->L[ir][jr] / 32768.f;
            bufprov->L[ir][jr] = bufmaskblurreti->L[ir][jr];
            bufprov->a[ir][jr] = bufmaskblurreti->a[ir][jr];
            bufprov->b[ir][jr] = bufmaskblurreti->b[ir][jr];

        }
    }

    if (rad > 0.f) {
        guidedFilter(guid, ble, ble, rad * 10.f / skip, 0.001, multiThread, 4);
    }

    LUTf lutTonemaskreti(65536);
    calcGammaLut(gamm, slop, lutTonemaskreti);
    float radiusb = 1.f / skip;

#ifdef _OPENMP
    #pragma omp parallel for schedule(dynamic,16)
#endif

    for (int ir = 0; ir < H_L; ir++)
        for (int jr = 0; jr < W_L; jr++) {
            float L_;
            bufmaskblurreti->L[ir][jr] = LIM01(ble[ir][jr]) * 32768.f;
            L_ = 2.f * bufmaskblurreti->L[ir][jr];
            bufmaskblurreti->L[ir][jr] = lutTonemaskreti[L_];
        }

    if (lmaskretilocalcurve && localmaskretiutili) {
#ifdef _OPENMP
        #pragma omp parallel for schedule(dynamic,16)
#endif

        for (int ir = 0; ir < H_L; ir++)
            for (int jr = 0; jr < W_L; jr++) {
                bufmaskblurreti->L[ir][jr] = 0.5f * lmaskretilocalcurve[2.f * bufmaskblurreti->L[ir][jr]];
            }
    }


    if (delt) {
        float *rdE[H_L] ALIGNED16;
        float *rdEBuffer = new float[H_L * W_L];

        for (int i = 0; i < H_L; i++) {
            rdE[i] = &rdEBuffer[i * W_L];
        }

        deltaEforMask(rdE, W_L, H_L, bufreti, hueref, chromaref, lumaref, maxdE, mindE, maxdElim, mindElim, iterat, limscope, scope, balance);
        // printf("rde1=%f rde2=%f\n", rdE[1][1], rdE[100][100]);
        std::unique_ptr<LabImage> delta(new LabImage(W_L, H_L));
#ifdef _OPENMP
        #pragma omp parallel for schedule(dynamic,16)
#endif

        for (int ir = 0; ir < H_L; ir++)
            for (int jr = 0; jr < W_L; jr++) {
                delta->L[ir][jr] = bufmaskblurreti->L[ir][jr] - bufprov->L[ir][jr];
                delta->a[ir][jr] = bufmaskblurreti->a[ir][jr] - bufprov->a[ir][jr];
                delta->b[ir][jr] = bufmaskblurreti->b[ir][jr] - bufprov->b[ir][jr];

                bufmaskblurreti->L[ir][jr] = bufprov->L[ir][jr] + rdE[ir][jr] * delta->L[ir][jr];
                bufmaskblurreti->a[ir][jr] = bufprov->a[ir][jr] + rdE[ir][jr] * delta->a[ir][jr];
                bufmaskblurreti->b[ir][jr] = bufprov->b[ir][jr] + rdE[ir][jr] * delta->b[ir][jr];
            }

        delete [] rdEBuffer;

    }


    if (lap > 0.f) {
        float *datain = new float[H_L * W_L];
        float *data_tmp = new float[H_L * W_L];

#ifdef _OPENMP
        #pragma omp parallel for
#endif

        for (int y = 0; y < H_L; y++) {
            for (int x = 0; x < W_L; x++) {
                datain[y * W_L + x] =  bufmaskblurreti->L[y][x];
            }
        }

        if (!pde) {
            ImProcFunctions::discrete_laplacian_threshold(data_tmp, datain, W_L, H_L, 200.f * lap);
        } else {
            ImProcFunctions::retinex_pde(datain, data_tmp, W_L, H_L, 12.f * lap, 1.f, nullptr, 0, 0, 1);
        }

#ifdef _OPENMP
        #pragma omp parallel for
#endif

        for (int y = 0; y < H_L; y++) {
            for (int x = 0; x < W_L; x++) {
                bufmaskblurreti->L[y][x] = data_tmp[y * W_L + x];
            }
        }

        delete [] datain;
        delete [] data_tmp;

    }

//blend
#ifdef _OPENMP
    #pragma omp parallel
#endif
    {
        gaussianBlur(bufmaskblurreti->L, bufmaskorigreti->L, W_L, H_L, radiusb);
        gaussianBlur(bufmaskblurreti->a, bufmaskorigreti->a, W_L, H_L, 1.f + (0.5f * rad) / skip);
        gaussianBlur(bufmaskblurreti->b, bufmaskorigreti->b, W_L, H_L, 1.f + (0.5f * rad) / skip);
    }

    float modr = 0.01f * (float) blend;

    if (llretiMask != 3 && retiMask) {
#ifdef _OPENMP
        #pragma omp parallel for schedule(dynamic,16)
#endif

        for (int y = 0; y < H_L; y++) {
            for (int x = 0; x < W_L; x++) {
                if (before == 0 && retiMasktmap) {
                    out[y][x] += fabs(modr) * bufmaskorigreti->L[y][x];
                    out[y][x] = LIM(out[y][x], 0.f, 100000.f);
                } else {
                    bufreti->L[y][x] +=  bufmaskorigreti->L[y][x] * modr;
                    bufreti->L[y][x] = CLIPLOC(bufreti->L[y][x]);

                }

                bufreti->a[y][x] *= (1.f + bufmaskorigreti->a[y][x] * modr * (1.f + 0.01f * chro));
                bufreti->b[y][x] *= (1.f + bufmaskorigreti->b[y][x] * modr * (1.f + 0.01f * chro));
                bufreti->a[y][x] = CLIPC(bufreti->a[y][x]);
                bufreti->b[y][x] = CLIPC(bufreti->b[y][x]);
            }
        }


    }

    if (!retiMasktmap  && retiMask) { //new original blur mask for deltaE
#ifdef _OPENMP
        #pragma omp parallel for schedule(dynamic,16)
#endif

        for (int y = 0; y < H_L; y++) {
            for (int x = 0; x < W_L; x++) {

                buforig->L[y][x] += (modr * bufmaskorigreti->L[y][x]);
                buforig->a[y][x] *= (1.f + modr * bufmaskorigreti->a[y][x]);
                buforig->b[y][x] *= (1.f + modr * bufmaskorigreti->b[y][x]);

                buforig->L[y][x] = CLIP(buforig->L[y][x]);
                buforig->a[y][x] = CLIPC(buforig->a[y][x]);
                buforig->b[y][x] = CLIPC(buforig->b[y][x]);

                buforig->L[y][x] = CLIP(buforig->L[y][x] - bufmaskorigreti->L[y][x]);
                buforig->a[y][x] = CLIPC(buforig->a[y][x] * (1.f - bufmaskorigreti->a[y][x]));
                buforig->b[y][x] = CLIPC(buforig->b[y][x] * (1.f - bufmaskorigreti->b[y][x]));
            }
        }

        float radius = 3.f / skip;

#ifdef _OPENMP
        #pragma omp parallel if (multiThread)
#endif
        {
            gaussianBlur(buforig->L, buforigmas->L, W_L, H_L, radius);
            gaussianBlur(buforig->a, buforigmas->a, W_L, H_L, radius);
            gaussianBlur(buforig->b, buforigmas->b, W_L, H_L, radius);
        }

    }


    if (llretiMask == 3) {

#ifdef _OPENMP
        #pragma omp parallel for schedule(dynamic,16)
#endif

        for (int y = 0; y < H_L; y++) {
            for (int x = 0; x < W_L; x++) {
                bufmask->L[y][x] = 6000.f + CLIPLOC(bufmaskorigreti->L[y][x]);
                bufmask->a[y][x] = CLIPC(bufreti->a[y][x] * bufmaskorigreti->a[y][x]);
                bufmask->b[y][x] = CLIPC(bufreti->b[y][x] * bufmaskorigreti->b[y][x]);
            }
        }
    }

}



void ImProcFunctions::MSRLocal(int sp, bool fftw, int lum, float** reducDE, LabImage * bufreti, LabImage * bufmask, LabImage * buforig, LabImage * buforigmas, float** luminance, float** templ, const float* const *originalLuminance,
                               const int width, const int height, int bfwr, int bfhr, const procparams::LocallabParams &loc, const int skip, const LocretigainCurve &locRETgainCcurve, const LocretitransCurve &locRETtransCcurve,
                               const int chrome, const int scall, const float krad, float &minCD, float &maxCD, float &mini, float &maxi, float &Tmean, float &Tsigma, float &Tmin, float &Tmax,
                               const LocCCmaskCurve & locccmasretiCurve, bool &lcmasretiutili, const  LocLLmaskCurve & locllmasretiCurve, bool &llmasretiutili, const  LocHHmaskCurve & lochhmasretiCurve, bool & lhmasretiutili, int llretiMask,
                               LUTf & lmaskretilocalcurve, bool & localmaskretiutili,
                               LabImage * transformed, bool retiMasktmap, bool retiMask,
                               bool delt, const float hueref, const float chromaref, const float lumaref,
                               float maxdE, float mindE, float maxdElim,  float mindElim, float iterat, float limscope, int scope, float balance)

{
    BENCHFUN
    bool py = true;

    if (py) {//enabled
        float mean, stddv, maxtr, mintr;
        float delta;
        constexpr float eps = 2.f;
        bool useHslLin = false;
        const float offse = loc.spots.at(sp).offs;
        const float chrT = (float)(loc.spots.at(sp).chrrt) / 100.f;
        const int scal = (loc.spots.at(sp).scalereti);
        float vart = loc.spots.at(sp).vart / 100.f;//variance
        const float strength = loc.spots.at(sp).str / 100.f; // Blend with original L channel data
        const float dar = loc.spots.at(sp).darkness;
        const float lig = loc.spots.at(sp).lightnessreti;
        float value = pow(strength, 0.4f);
        float value_1 = pow(strength, 0.3f);
        bool logli = loc.spots.at(sp).loglin;
        float limD = loc.spots.at(sp).limd;//10.f
        limD = pow(limD, 1.7f);  //about 2500 enough
        float ilimD = 1.f / limD;
        float threslum = loc.spots.at(sp).limd;
        const float elogt = 2.71828f;

        if (!logli) {
            useHslLin = true;
        }

        //empirical skip evaluation : very difficult  because quasi all parameters interfere
        //to test on several images
        int nei = (int)(krad * loc.spots.at(sp).neigh);
        // printf("neigh=%i\n", nei);
        //several test to find good values ???!!!
        //very difficult to do because 4 factor are correlate with skip and cannot been solved easily
        // size of spots
        // radius - neigh
        // scal
        // variance vart
        //not too bad proposition
        float divsca = 1.f;

        if (scal >= 3) {
            divsca = sqrt(scal / 3.f);
        }

        if (skip >= 4) {
            //nei = (int)(0.1f * nei + 2.f);     //not too bad
            nei = (int)(nei / (1.5f * skip)) / divsca;
            vart *= sqrt(skip);
        } else if (skip > 1 && skip < 4) {
            //nei = (int)(0.3f * nei + 2.f);
            nei = (int)(nei / skip) / divsca;
            vart *= sqrt(skip);
        }

        int moderetinex = 0;

        if (loc.spots.at(sp).retinexMethod == "uni") {
            moderetinex = 0;
        } else if (loc.spots.at(sp).retinexMethod == "low") {
            moderetinex = 1;
        } else if (loc.spots.at(sp).retinexMethod == "high") {
            moderetinex = 2;
        }

        const float high = 0.f; // Dummy to pass to retinex_scales(...)

        constexpr auto maxRetinexScales = 10;
        float RetinexScales[maxRetinexScales];

        retinex_scales(RetinexScales, scal, moderetinex, nei, high);


        const int H_L = height;
        const int W_L = width;

        float *src[H_L] ALIGNED16;
        float *srcBuffer = new float[H_L * W_L];

        for (int i = 0; i < H_L; i++) {
            src[i] = &srcBuffer[i * W_L];
        }


#ifdef _OPENMP
        #pragma omp parallel for
#endif

        for (int i = 0; i < H_L; i++)
            for (int j = 0; j < W_L; j++) {
                src[i][j] = luminance[i][j] + eps;
                luminance[i][j] = 0.f;
            }

        float *out[H_L] ALIGNED16;
        float *outBuffer = new float[H_L * W_L];

        for (int i = 0; i < H_L; i++) {
            out[i] = &outBuffer[i * W_L];
        }

        float clipt = loc.spots.at(sp).cliptm;

        const float logBetaGain = xlogf(16384.f);
        float pond = logBetaGain / (float) scal;

        if (!useHslLin) {
            pond /= log(elogt);
        }

        float *buffer = new float[W_L * H_L];
        float kr = 1.f;//on FFTW
        float kg = 1.f;//on Gaussianblur

        for (int scale = scal - 1; scale >= 0; scale--) {
            //    printf("retscale=%f scale=%i \n", mulradiusfftw * RetinexScales[scale], scale);
            //emprical adjustement between FFTW radius and Gaussainblur
            //under 50 ==> 10.f
            // 400 ==> 1.f
            float sigm = 1.f;

            if (settings->fftwsigma == false) { //empirical formula
                sigm = RetinexScales[scale];
                float ak = -9.f / 350.f;
                float bk = 10.f - 50.f * ak;
                kr = ak * sigm + bk;

                if (sigm < 50.f) {
                    kr = 10.f;
                }

                //above 400 at 5000 ==> 20.f
                if (sigm > 400.f) { //increase ==> 5000
                    float ka = 19.f / 4600.f;
                    float kb = 1.f - 400 * ka;
                    kr = ka * sigm + kb;
                    float kga = -0.14f / 4600.f;//decrease
                    float kgb = 1.f - 400.f * kga;
                    kg = kga * sigm + kgb;

                    if (sigm > 5000.f) {
                        kr = ka * 5000.f + kb;
                        kg = kga * 5000.f + kgb;
                    }

                }
            } else {//sigma *= sigma
                kg = 1.f;
                kr = sigm;
            }

            if (!fftw) {
#ifdef _OPENMP
                #pragma omp parallel  //disabled with FFTW
#endif
                {

                    if (scale == scal - 1)
                    {
                        gaussianBlur(src, out, W_L, H_L, kg * RetinexScales[scale], buffer);
                    } else   // reuse result of last iteration
                    {
                        // out was modified in last iteration => restore it
                        gaussianBlur(out, out, W_L, H_L, sqrtf(SQR(kg * RetinexScales[scale]) - SQR(kg * RetinexScales[scale + 1])), buffer);
                    }
                }
            } else {
                if (scale == scal - 1) {
                    if (settings->fftwsigma == false) { //empirical formula
                        ImProcFunctions::fftw_convol_blur2(src, out, bfwr, bfhr, (kr * RetinexScales[scale]), 0, 0);
                    } else {
                        ImProcFunctions::fftw_convol_blur2(src, out, bfwr, bfhr, (SQR(RetinexScales[scale])), 0, 0);
                    }
                } else { // reuse result of last iteration
                    // out was modified in last iteration => restore it
                    if (settings->fftwsigma == false) { //empirical formula
                        ImProcFunctions::fftw_convol_blur2(out, out, bfwr, bfhr, sqrtf(SQR(kr * RetinexScales[scale]) - SQR(kr * RetinexScales[scale + 1])), 0, 0);
                    } else {
                        ImProcFunctions::fftw_convol_blur2(out, out, bfwr, bfhr, (SQR(RetinexScales[scale]) - SQR(RetinexScales[scale + 1])), 0, 0);
                    }
                }
            }

            if (scale == 1) { //equalize last scale with darkness and lightness of course acts on TM!

                if (dar != 1.f || lig != 1.f) {

#ifdef _OPENMP
                    #pragma omp parallel for
#endif

                    for (int y = 0; y < H_L; ++y) {
                        for (int x = 0; x < W_L; ++x) {
                            float buf = (src[y][x] - out[y][x]) * value;
                            buf *= (buf > 0.f) ? lig : dar;
                            out[y][x] = LIM(out[y][x] + buf, -100000.f, 100000.f);
                        }
                    }
                }
            }

            /*
                        if (lum == 1 && scale == 1 && (llretiMask == 3 || llretiMask == 0 || llretiMask == 2 || llretiMask == 4)) { //only mask with luminance on last scale
                            int before = 0;
                            maskforretinex(sp, before, luminance, out, W_L, H_L, skip, locccmasretiCurve, lcmasretiutili, locllmasretiCurve, llmasretiutili, lochhmasretiCurve, lhmasretiutili, llretiMask, retiMasktmap, retiMask,
                                           loc, bufreti, bufmask, buforig, buforigmas, multiThread);
                        }
            */

#ifdef __SSE2__
            vfloat pondv = F2V(pond);
            vfloat limMinv = F2V(ilimD);
            vfloat limMaxv = F2V(limD);
#endif
#ifdef _OPENMP
            #pragma omp parallel for
#endif

            for (int i = 0; i < H_L; i++) {
                int j = 0;

#ifdef __SSE2__

                if (useHslLin) {
                    for (; j < W_L - 3; j += 4) {
                        _mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * (vclampf(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv)));
                    }
                } else {
                    for (; j < W_L - 3; j += 4) {
                        _mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv *  xlogf(vclampf(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv)));
                    }
                }

#endif

                if (useHslLin) {
                    for (; j < W_L; j++) {
                        luminance[i][j] +=  pond * (LIM(src[i][j] / out[i][j], ilimD, limD));
                    }
                } else {
                    for (; j < W_L; j++) {
                        luminance[i][j] +=  pond * xlogf(LIM(src[i][j] / out[i][j], ilimD, limD));
                    }
                }
            }

        }


        if (scal == 1) {//only if user select scal = 1
            float kval = 1.f;
#ifdef _OPENMP
            #pragma omp parallel for
#endif

            for (int y = 0; y < H_L; ++y) {
                for (int x = 0; x < W_L; ++x) {
                    float threslow = threslum * 163.f;

                    if (src[y][x] < threslow) {
                        kval = src[y][x] / threslow;
                    }
                }
            }


#ifdef _OPENMP
            #pragma omp parallel for
#endif

            for (int y = 0; y < H_L; ++y) {
                for (int x = 0; x < W_L; ++x) {
                    float buf = (src[y][x] - out[y][x]) * value_1;
                    buf *= (buf > 0.f) ? lig : dar;
                    luminance[y][x] = LIM(src[y][x] + (1.f + kval) * buf, -32768.f, 32768.f);
                }
            }

            double avg = 0.f;
            int ng = 0;

#ifdef _OPENMP
            #pragma omp parallel for
#endif

            for (int i = 0; i < H_L; i++) {
                for (int j = 0; j < W_L; j++) {
                    avg += luminance[i][j];
                    ng++;
                }
            }

            avg /= ng;
            avg /= 32768.f;
            avg = LIM01(avg);
            float contreal = 0.5f * vart;
            DiagonalCurve reti_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
#endif

            for (int i = 0; i < H_L; i++)
                for (int j = 0; j < W_L; j++) {
                    float buf = LIM01(luminance[i][j] / 32768.f);
                    buf = reti_contrast.getVal(buf);
                    buf *= 32768.f;
                    luminance[i][j] = buf;
                }
        }


        delete [] buffer;
        delete [] outBuffer;
        outBuffer = nullptr;
        delete [] srcBuffer;
        srcBuffer = nullptr;

        float str = strength * (chrome == 0 ? 1.f : 0.8f * (chrT - 0.4f));
        const float maxclip = (chrome == 0 ? 32768.f : 50000.f);

        if (scal != 1) {
            mean = 0.f;
            stddv = 0.f;

            mean_stddv2(luminance, mean, stddv, W_L, H_L, maxtr, mintr);
            // printf("mean=%f std=%f delta=%f maxtr=%f mintr=%f\n", mean, stddv, delta, maxtr, mintr);

            if (locRETtransCcurve && mean != 0.f && stddv != 0.f) { //if curve
                float asig = 0.166666f / stddv;
                float bsig = 0.5f - asig * mean;
                float amax = 0.333333f / (maxtr - mean - stddv);
                float bmax = 1.f - amax * maxtr;
                float amin = 0.333333f / (mean - stddv - mintr);
                float bmin = -amin * mintr;

                asig *= 500.f;
                bsig *= 500.f;
                amax *= 500.f;
                bmax *= 500.f;
                amin *= 500.f;
                bmin *= 500.f;

#ifdef _OPENMP
                #pragma omp parallel
#endif
                {
                    float absciss;
#ifdef _OPENMP
                    #pragma omp for schedule(dynamic,16)
#endif

                    for (int i = 0; i < H_L; i++)
                        for (int j = 0; j < W_L; j++) { //for mintr to maxtr evalate absciss in function of original transmission
                            if (LIKELY(fabsf(luminance[i][j] - mean) < stddv)) {
                                absciss = asig * luminance[i][j] + bsig;
                            } else if (luminance[i][j] >= mean) {
                                absciss = amax * luminance[i][j] + bmax;
                            } else { /*if(luminance[i][j] <= mean - stddv)*/
                                absciss = amin * luminance[i][j] + bmin;
                            }

                            //TODO : move multiplication by 4.f and subtraction of 1.f inside the curve
                            luminance[i][j] *= (-1.f + 4.f * locRETtransCcurve[absciss]); //new transmission

                        }
                }
            }

            mean = 0.f;
            stddv = 0.f;
            mean_stddv2(luminance, mean, stddv, W_L, H_L, maxtr, mintr);//new calculation of mean...

            float epsil = 0.1f;

            mini = mean - vart * stddv;

            if (mini < mintr) {
                mini = mintr + epsil;
            }

            maxi = mean + vart * stddv;

            if (maxi > maxtr) {
                maxi = maxtr - epsil;
            }

            delta = maxi - mini;
            //   printf("maxi=%f mini=%f mean=%f std=%f delta=%f maxtr=%f mintr=%f\n", maxi, mini, mean, stddv, delta, maxtr, mintr);

            if (!delta) {
                delta = 1.0f;
            }


            float *copylum[H_L] ALIGNED16;
            float *copylumBuffer = new float[H_L * W_L];

            for (int i = 0; i < H_L; i++) {
                copylum[i] = &copylumBuffer[i * W_L];
            }

            float cdfactor = (clipt * 32768.f) / delta;
            maxCD = -9999999.f;
            minCD = 9999999.f;
            //prepare work for curve gain
#ifdef _OPENMP
            #pragma omp parallel for
#endif

            for (int i = 0; i < H_L; i++) {
                for (int j = 0; j < W_L; j++) {
                    luminance[i][j] = luminance[i][j] - mini;
                }
            }

            mean = 0.f;
            stddv = 0.f;

            mean_stddv2(luminance, mean, stddv, W_L, H_L, maxtr, mintr);
//         printf("meanun=%f stdun=%f maxtr=%f mintr=%f\n", mean, stddv, maxtr, mintr);

            float asig = 0.f, bsig = 0.f, amax = 0.f, bmax = 0.f, amin = 0.f, bmin = 0.f;
            const bool hasRetGainCurve =  locRETgainCcurve && mean != 0.f && stddv != 0.f;

            if (hasRetGainCurve) { //if curve
                asig = 0.166666f / stddv;
                bsig = 0.5f - asig * mean;
                amax = 0.333333f / (maxtr - mean - stddv);
                bmax = 1.f - amax * maxtr;
                amin = 0.333333f / (mean - stddv - mintr);
                bmin = -amin * mintr;

                asig *= 500.f;
                bsig *= 500.f;
                amax *= 500.f;
                bmax *= 500.f;
                amin *= 500.f;
                bmin *= 500.f;
                cdfactor *= 2.f;
            }


#ifdef _OPENMP
            #pragma omp parallel
#endif
            {
                //            float absciss;
                float cdmax = -999999.f, cdmin = 999999.f;
                float gan = 0.5f;

#ifdef _OPENMP
                #pragma omp for schedule(dynamic,16)
#endif

                for (int i = 0; i < H_L; i ++)
                    for (int j = 0; j < W_L; j++) {

                        if (hasRetGainCurve) {
                            float absciss;

                            if (LIKELY(fabsf(luminance[i][j] - mean) < stddv)) {
                                absciss = asig * luminance[i][j] + bsig;
                            } else if (luminance[i][j] >= mean) {
                                absciss = amax * luminance[i][j] + bmax;
                            } else {
                                absciss = amin * luminance[i][j] + bmin;
                            }

                            gan = locRETgainCcurve[absciss]; //new gain function transmission
                        }

                        //but we don't update mean stddv for display only...
                        copylum[i][j] = gan * luminance[i][j];//update datas for display
                        float cd = gan * cdfactor * luminance[i][j] + offse;

                        cdmax = cd > cdmax ? cd : cdmax;
                        cdmin = cd < cdmin ? cd : cdmin;
                        luminance[i][j] = intp(str * reducDE[i][j], clipretinex(cd, 0.f, maxclip), originalLuminance[i][j]);
                    }



#ifdef _OPENMP
                #pragma omp critical
#endif
                {
                    maxCD = maxCD > cdmax ? maxCD : cdmax;
                    minCD = minCD < cdmin ? minCD : cdmin;
                }
            }
            mean = 0.f;
            stddv = 0.f;

            mean_stddv2(copylum, mean, stddv, W_L, H_L, maxtr, mintr);
            delete [] copylumBuffer;
            copylumBuffer = nullptr;

//         printf("mean=%f std=%f maxtr=%f mintr=%f\n", mean, stddv, maxtr, mintr);

        } else {
#ifdef _OPENMP
            #pragma omp for schedule(dynamic,16)
#endif

            for (int i = 0; i < H_L; i ++)
                for (int j = 0; j < W_L; j++) {
                    luminance[i][j] =   LIM(luminance[i][j], 0.f, maxclip) * str + (1.f - str) * originalLuminance[i][j];

                }

        }

        float rad = loc.spots.at(sp).radmaskreti;
        float slop = loc.spots.at(sp).slomaskreti;
        float gamm = loc.spots.at(sp).gammaskreti;
        float blend = loc.spots.at(sp).blendmaskreti;
        float chro = loc.spots.at(sp).chromaskreti;
        float lap = loc.spots.at(sp).lapmaskreti;
        bool pde = params->locallab.spots.at(sp).laplac;

        if (lum == 1  && (llretiMask == 3 || llretiMask == 0 || llretiMask == 2 || llretiMask == 4)) { //only mask with luminance on last scale
            int before = 1;
            maskforretinex(sp, before, luminance, nullptr, W_L, H_L, skip,
                           locccmasretiCurve, lcmasretiutili, locllmasretiCurve, llmasretiutili, lochhmasretiCurve, lhmasretiutili, llretiMask, retiMasktmap, retiMask,
                           rad, lap, pde, gamm, slop, chro, blend,
                           lmaskretilocalcurve, localmaskretiutili,
                           bufreti, bufmask, buforig, buforigmas, multiThread,
                           delt, hueref, chromaref, lumaref,
                           maxdE, mindE, maxdElim, mindElim, iterat, limscope, scope, balance
                          );
        }

        //mask does not interfered with datas displayed

        Tmean = mean;
        Tsigma = stddv;
        Tmin = mintr;
        Tmax = maxtr;
    }
}
}