From 52d2f4e6e49d0d92f0758135d2af952b299d4ad7 Mon Sep 17 00:00:00 2001 From: Desmis Date: Fri, 7 Jun 2019 15:17:05 +0200 Subject: [PATCH] Improve retinex_pde --- rtdata/languages/default | 2 +- rtengine/improcfun.h | 4 +- rtengine/iplocallab.cc | 99 +++++++++++++++++++++++++++++++++++----- rtgui/locallab.cc | 4 +- 4 files changed, 93 insertions(+), 16 deletions(-) diff --git a/rtdata/languages/default b/rtdata/languages/default index 79cd54438..41d85d171 100644 --- a/rtdata/languages/default +++ b/rtdata/languages/default @@ -2030,7 +2030,7 @@ TP_LOCALLAB_DEHAZ;Dehaze TP_LOCALLAB_GRIDONE;Color Toning TP_LOCALLAB_GRIDTWO;Direct TP_LOCALLAB_LUM;Curves LC -TP_LOCALLAB_LAPLACE;Laplacian threshold +TP_LOCALLAB_LAPLACE;Laplacian threshold deltaE TP_LOCALLAB_HLH;Curves H TP_LOCALLAB_CHROMACBDL;Chroma TP_LOCALLAB_CHROMACB_TOOLTIP;Acts as an amplifier-reducer action compare to sliders of luminance.\nUnder 100 reduce, above 100 amplifie diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h index 947e88eef..acbb778d0 100644 --- a/rtengine/improcfun.h +++ b/rtengine/improcfun.h @@ -299,8 +299,8 @@ public: int pitch, int scale, const int luma, const int chroma/*, LUTf & Lcurve, LUTf & abcurve*/); //locallab void normalize_mean_dt(float *data, const float *ref, size_t size); - void retinex_pde(float *datain, float * dataout, int bfw, int bfh, float thresh, float multy, int numThreads); - void fftw_convol_blur(float *input, float *output, int bfw, int bfh, float radius, int sigmafft); + void retinex_pde(float *datain, float * dataout, int bfw, int bfh, float thresh, float multy, float *dE); + void fftw_convol_blur(float *input, float *output, int bfw, int bfh, float radius); void MSRLocal(int sp, int lum, LabImage * bufreti, LabImage * bufmask, LabImage * buforig, LabImage * buforigmas, float** luminance, float** templ, const float* const *originalLuminance, const int width, const int height, const procparams::LocallabParams &loc, const int skip, const LocretigainCurve &locRETgainCcurve, 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 LocCCmaskretiCurve & locccmasretiCurve, bool &lcmasretiutili, const LocLLmaskretiCurve & locllmasretiCurve, bool &llmasretiutili, const LocHHmaskretiCurve & lochhmasretiCurve, bool & lhmasretiutili, int llretiMask, LabImage * transformed, bool retiMasktmap, bool retiMask); diff --git a/rtengine/iplocallab.cc b/rtengine/iplocallab.cc index 5dab6a77d..714ad60b9 100644 --- a/rtengine/iplocallab.cc +++ b/rtengine/iplocallab.cc @@ -1955,6 +1955,54 @@ static void blendmask(const local_params& lp, int xstart, int ystart, int cx, in } } + +static void deltaEforLaplace (float *dE, const local_params& lp, int bfw, int bfh, LabImage* bufexporig, const float hueref, const float chromaref, const float lumaref) +{ + + const float refa = chromaref * cos(hueref); + const float refb = chromaref * sin(hueref); + const float refL = lumaref; + float maxdE = 5.f + MAXSCOPE * lp.lap; + float *dEforLaplace = new float [bfw * bfh]; + float maxC = sqrt((SQR(refa - bufexporig->a[0][0] / 327.68f) + SQR(refb - bufexporig->b[0][0]/327.68f)) + SQR(refL - bufexporig->L[0][0]/327.68f)); + // float sumde = 0.f; +#ifdef _OPENMP + #pragma omp parallel for reduction(max:maxC) // reduction(+:sumde) +#endif + for (int y = 0; y < bfh; y++) { + for (int x = 0; x < bfw; x++) { + dEforLaplace[y * bfw + x] = sqrt((SQR(refa - bufexporig->a[y][x] / 327.68f) + SQR(refb - bufexporig->b[y][x]/327.68f)) + SQR(refL - bufexporig->L[y][x]/327.68f)); + maxC = rtengine::max(maxC, dEforLaplace[y * bfw + x]); + // sumde += dEforLaplace[y * bfw + x]; + } + } + // float mxde = sumde /(bfh * bfw); + // maxC = 0.5f * (mxde + maxC); + if(maxdE > maxC) { + maxdE = maxC - 1.f; + } + float ade = 1.f / (maxdE - maxC); + float bde = -ade * maxC; + +#ifdef _OPENMP + #pragma omp parallel for schedule(dynamic,16) +#endif + for (int y = 0; y < bfh; y++) { + for (int x = 0; x < bfw; x++) { + + float reducdEforLap = 1.f; + if(dEforLaplace[y * bfw + x] < maxdE) { + reducdEforLap = 1.f; + } else { + reducdEforLap = ade * dEforLaplace[y * bfw + x] + bde; + } + dE[y * bfw + x] = reducdEforLap; + } + } + delete [] dEforLaplace; +} + + static void showmask(const local_params& lp, int xstart, int ystart, int cx, int cy, int bfw, int bfh, LabImage* bufexporig, LabImage* transformed, LabImage* bufmaskorigSH) { #ifdef _OPENMP @@ -3761,6 +3809,7 @@ static float *discrete_laplacian_threshold(float *data_out, const float *data_in if (fabs(diff) > t) *ptr_out += diff; } + ptr_in++; ptr_in_xm1++; ptr_in_xp1++; @@ -3773,7 +3822,7 @@ static float *discrete_laplacian_threshold(float *data_out, const float *data_in return data_out; } -void ImProcFunctions::retinex_pde(float *datain, float * dataout, int bfw, int bfh, float thresh, float multy, int numThreads) +void ImProcFunctions::retinex_pde(float *datain, float * dataout, int bfw, int bfh, float thresh, float multy, float *dE) { /* * Copyright 2009-2011 IPOL Image Processing On Line http://www.ipol.im/ @@ -3783,30 +3832,57 @@ void ImProcFunctions::retinex_pde(float *datain, float * dataout, int bfw, int b * @brief laplacian, DFT and Poisson routines * * @author Nicolas Limare + * adapted by Jacques Desmis 2019 */ - fftwf_plan dct_fw, dct_bw; - float *data_fft, *data_tmp, *data; - + fftwf_plan dct_fw, dct_fw04, dct_bw; + float *data_fft, *data_fft04, *data_tmp, *data, *data_tmp04; if (NULL == (data_tmp = (float *) fftwf_malloc(sizeof(float) * bfw * bfh))) { fprintf(stderr, "allocation error\n"); abort(); } - + if (NULL == (data_tmp04 = (float *) fftwf_malloc(sizeof(float) * bfw * bfh))) { + fprintf(stderr, "allocation error\n"); + abort(); + } + //first call to laplacian with plein strength (void) discrete_laplacian_threshold(data_tmp, datain, bfw, bfh, thresh); if (NULL == (data_fft = (float *) fftwf_malloc(sizeof(float) * bfw * bfh))) { fprintf(stderr, "allocation error\n"); abort(); } + //second call to laplacian with 40% strength ==> reduce effect if we are far from ref (deltaE) + (void) discrete_laplacian_threshold(data_tmp04, datain, bfw, bfh, 0.4f * thresh); + if (NULL == (data_fft04 = (float *) fftwf_malloc(sizeof(float) * bfw * bfh))) { + fprintf(stderr, "allocation error\n"); + abort(); + } + if (NULL == (data = (float *) fftwf_malloc(sizeof(float) * bfw * bfh))) { fprintf(stderr, "allocation error\n"); abort(); } - + //execute first dct_fw = fftwf_plan_r2r_2d(bfh, bfw, data_tmp, data_fft, FFTW_REDFT10, FFTW_REDFT10, FFTW_ESTIMATE | FFTW_DESTROY_INPUT); fftwf_execute(dct_fw); - fftwf_free(data_tmp); + //execute second + dct_fw04 = fftwf_plan_r2r_2d(bfh, bfw, data_tmp04, data_fft04, FFTW_REDFT10, FFTW_REDFT10, FFTW_ESTIMATE | FFTW_DESTROY_INPUT); + fftwf_execute(dct_fw04); +#ifdef _OPENMP + #pragma omp parallel for +#endif + for (int y = 0; y < bfh ; y++) {//mix two fftw Laplacian : plein if dE near ref + for (int x = 0; x < bfw; x++) { + float prov = pow(dE[y * bfw + x], 4.5f); + data_fft[y * bfw + x] = prov * data_fft[y * bfw + x] + (1.f - prov) * data_fft04[y * bfw + x]; + } + } + fftwf_free(data_fft04); + fftwf_free(data_tmp); + fftwf_free(data_tmp04); + fftwf_destroy_plan(dct_fw04); + /* solve the Poisson PDE in Fourier space */ /* 1. / (float) (bfw * bfh)) is the DCT normalisation term, see libfftw */ (void) retinex_poisson_dct(data_fft, bfw, bfh, 1./(double) (bfw * bfh)); @@ -3829,7 +3905,7 @@ void ImProcFunctions::retinex_pde(float *datain, float * dataout, int bfw, int b } } -void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int bfh, float radius, int sigmafft) +void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int bfh, float radius) { float *out; // *outfft; fftwf_plan p;// ps; @@ -5932,19 +6008,19 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o } else if(lp.softmet == 1) { float *datain = new float[bfw*bfh]; float *dataout = new float[bfw*bfh]; + float *dE = new float[bfw*bfh]; + deltaEforLaplace (dE, lp, bfw, bfh, bufexpfin.get(), hueref, chromaref, lumaref); #ifdef _OPENMP #pragma omp parallel for schedule(dynamic,16) #endif for (int y = 0; y < bfh; y++) { for (int x = 0; x < bfw; x++) { - // datain[y * bfw + x] = bufexpfin->L[y][x] / (10.9f * lp.lap); datain[y * bfw + x] = bufexpfin->L[y][x]; } } - // ImProcFunctions::retinex_pde(datain, dataout, bfw, bfh, 5.f * lp.strng, 10.9f * lp.lap, 1); - ImProcFunctions::retinex_pde(datain, dataout, bfw, bfh, 8.f * lp.strng, 1.f, 1); + ImProcFunctions::retinex_pde(datain, dataout, bfw, bfh, 8.f * lp.strng, 1.f, dE); #ifdef _OPENMP #pragma omp parallel for schedule(dynamic,16) #endif @@ -5955,6 +6031,7 @@ void ImProcFunctions::Lab_Local(int call, int sp, float** shbuffer, LabImage * o } delete [] datain; delete [] dataout; + delete [] dE; } #ifdef _OPENMP #pragma omp parallel for schedule(dynamic,16) diff --git a/rtgui/locallab.cc b/rtgui/locallab.cc index 3b4992e97..abae8e982 100644 --- a/rtgui/locallab.cc +++ b/rtgui/locallab.cc @@ -133,7 +133,7 @@ Locallab::Locallab(): sensiv(Gtk::manage(new Adjuster(M("TP_LOCALLAB_SENSI"), 0, 100, 1, 15))), //Soft Light streng(Gtk::manage(new Adjuster(M("TP_LOCALLAB_STRENG"), 1, 100, 1, 1))), - laplace(Gtk::manage(new Adjuster(M("TP_LOCALLAB_LAPLACE"), 0.5, 60., 0.1, 20.))), + laplace(Gtk::manage(new Adjuster(M("TP_LOCALLAB_LAPLACE"), 0., 100., 0.5, 20.))), sensisf(Gtk::manage(new Adjuster(M("TP_LOCALLAB_SENSI"), 1, 100, 1, 15))), // Blur & Noise radius(Gtk::manage(new Adjuster(M("TP_LOCALLAB_RADIUS"), 1.0, 100.0, 0.1, 1.0))), @@ -743,7 +743,7 @@ Locallab::Locallab(): ToolParamBlock* const softBox = Gtk::manage(new ToolParamBlock()); softBox->pack_start(*softMethod); softBox->pack_start(*streng); -// softBox->pack_start(*laplace); + softBox->pack_start(*laplace); softBox->pack_start(*sensisf); expsoft->add(*softBox, false); expsoft->setLevel(2);