From fa92226ba3e24645fd3fa682a4389fd75642aaf8 Mon Sep 17 00:00:00 2001 From: Desmis Date: Sat, 29 Jun 2019 15:53:26 +0200 Subject: [PATCH] Fixed bad behavior sigma FFTW --- rtengine/improcfun.h | 4 ++-- rtengine/iplocalcontrast.cc | 4 ++-- rtengine/iplocallab.cc | 38 ++++++++++++++++++++++++------------- rtengine/ipretinex.cc | 6 +++--- 4 files changed, 32 insertions(+), 20 deletions(-) diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h index 3f7a9e455..b121bd20d 100644 --- a/rtengine/improcfun.h +++ b/rtengine/improcfun.h @@ -300,8 +300,8 @@ public: //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, float *dE); - void fftw_convol_blur(float *input, float *output, int bfw, int bfh, float radius, int fftkern); - void fftw_convol_blur2(float **input2, float **output2, int bfw, int bfh, float radius, int fftkern); + void fftw_convol_blur(float *input, float *output, int bfw, int bfh, float radius, int fftkern, int algo); + void fftw_convol_blur2(float **input2, float **output2, int bfw, int bfh, float radius, int fftkern, int algo); void fftw_tile_blur(int GW, int GH, int tilssize , int max_numblox_W, int min_numblox_W, float **tmp1, int numThreads, double radius); void MSRLocal(int sp, bool fftw, 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, diff --git a/rtengine/iplocalcontrast.cc b/rtengine/iplocalcontrast.cc index fdc5fb61f..b35c3d3da 100644 --- a/rtengine/iplocalcontrast.cc +++ b/rtengine/iplocalcontrast.cc @@ -55,8 +55,8 @@ void ImProcFunctions::localContrast(LabImage *lab, float **destination, const Lo #endif gaussianBlur(lab->L, buf, width, height, sigma); } else { - sigma *= 50.f;//approximation to convert sigma "gaussianBlur" to FFTW blur - ImProcFunctions::fftw_convol_blur2(lab->L, buf, width, height, sigma, 0); + sigma *= 1.f; + ImProcFunctions::fftw_convol_blur2(lab->L, buf, width, height, sigma, 0, 0); } #ifdef _OPENMP #pragma omp parallel for if(multiThread) diff --git a/rtengine/iplocallab.cc b/rtengine/iplocallab.cc index 4b9ee24f2..a85eefbde 100644 --- a/rtengine/iplocallab.cc +++ b/rtengine/iplocallab.cc @@ -3994,7 +3994,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 fftkern) +void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int bfh, float radius, int fftkern, int algo) { /* ** Jacques Desmis june 2019 - inspired by Copyright 2013 IPOL Image Processing On Line http://www.ipol.im/ @@ -4026,7 +4026,9 @@ void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int fftwf_plan p; fftwf_plan pkern;//plan for FFT int image_size, image_sizechange; - float n_x, n_y;//relative coordonates for kernel Gauss + float n_x = 1.f; + float n_y = 1.f;//relative coordonates for kernel Gauss + float radsig = 1.f; out = (float*) fftwf_malloc(sizeof(float) * (bfw * bfh));//allocate real datas for FFT @@ -4043,14 +4045,17 @@ void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int fftwf_execute(p); fftwf_destroy_plan(p); /*define the gaussian constants for the convolution kernel*/ - - // n_x = rtengine::RT_PI/(double) bfw;//ipol - // n_y = rtengine::RT_PI/(double) bfh; - n_x = 1.f/(float) bfw;//gauss - n_y = 1.f/(float) bfh; + if(algo == 0) { + n_x = rtengine::RT_PI/(double) bfw;//ipol + n_y = rtengine::RT_PI/(double) bfh; + } else if(algo == 1) { + n_x = 1.f/(float) bfw;//gauss + n_y = 1.f/(float) bfh; + radsig = 1.f / (2.f * rtengine::RT_PI * radius * radius);//gauss + } n_x = n_x * n_x; n_y = n_y * n_y; - float radsig = 1.f / (2.f * rtengine::RT_PI * radius * radius);//gauss + image_size = bfw * bfh; image_sizechange = 4 * image_size; @@ -4061,8 +4066,11 @@ void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int for(int j = 0; j < bfh; j++){ int index = j * bfw; for(int i = 0; i < bfw; i++) - // kern[ i+ index] = exp((float)(-radius)*(n_x * i * i + n_y * j * j));//calculate Gauss kernel Ipol formula - kern[ i+ index] = radsig * exp((float)(-(n_x * i * i + n_y * j * j)/ (2.f * radius * radius)));//calculate Gauss kernel with Gauss formula + if(algo == 0) { + kern[ i+ index] = exp((float)(-radius)*(n_x * i * i + n_y * j * j));//calculate Gauss kernel Ipol formula + } else if(algo == 1) { + kern[ i+ index] = radsig * exp((float)(-(n_x * i * i + n_y * j * j)/ (2.f * radius * radius)));//calculate Gauss kernel with Gauss formula + } } /*compute the Fourier transform of the kernel data*/ pkern = fftwf_plan_r2r_2d(bfh, bfw, kern, outkern, FFTW_REDFT10, FFTW_REDFT10,FFTW_ESTIMATE);//FFT 2 dimensions forward @@ -4088,7 +4096,11 @@ void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int for(int j = 0; j < bfh; j++){ int index = j * bfw; for(int i = 0; i < bfw; i++) - out[i + index] *= exp((float)(-radius)*(n_x * i * i + n_y * j * j));//apply Gauss kernel whithout FFT + if(algo == 0) { + out[i + index] *= exp((float)(-radius)*(n_x * i * i + n_y * j * j));//apply Gauss kernel whithout FFT + } else if(algo == 1) { + out[i + index] *= radsig * exp((float)(-(n_x * i * i + n_y * j * j)/ (2.f * radius * radius)));//calculate Gauss kernel with Gauss formula + } } } @@ -4107,7 +4119,7 @@ void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int } } -void ImProcFunctions::fftw_convol_blur2(float **input2, float **output2, int bfw, int bfh, float radius, int fftkern) +void ImProcFunctions::fftw_convol_blur2(float **input2, float **output2, int bfw, int bfh, float radius, int fftkern, int algo) { MyMutex::MyLock lock(*fftwMutex); @@ -4131,7 +4143,7 @@ void ImProcFunctions::fftw_convol_blur2(float **input2, float **output2, int bfw input[y * bfw + x] = input2[y][x]; } } - ImProcFunctions::fftw_convol_blur(input, output, bfw, bfh, radius, fftkern); + ImProcFunctions::fftw_convol_blur(input, output, bfw, bfh, radius, fftkern, algo); #ifdef _OPENMP #pragma omp parallel for schedule(dynamic,16) diff --git a/rtengine/ipretinex.cc b/rtengine/ipretinex.cc index eae962483..f15f99d6e 100644 --- a/rtengine/ipretinex.cc +++ b/rtengine/ipretinex.cc @@ -961,7 +961,7 @@ void ImProcFunctions::MSRLocal(int sp, bool fftw, int lum, LabImage * bufreti, L } float *buffer = new float[W_L * H_L]; - float mulradiusfftw = 40.f; + float mulradiusfftw = 1.f; for (int scale = scal - 1; scale >= 0; scale--) { // printf("retscale=%f scale=%i \n", mulradiusfftw * RetinexScales[scale], scale); if(!fftw) { @@ -982,11 +982,11 @@ void ImProcFunctions::MSRLocal(int sp, bool fftw, int lum, LabImage * bufreti, L } else { if (scale == scal - 1) { - ImProcFunctions::fftw_convol_blur2(src, out, W_L, H_L, mulradiusfftw * RetinexScales[scale], 0); + ImProcFunctions::fftw_convol_blur2(src, out, W_L, H_L, mulradiusfftw * RetinexScales[scale], 0, 0); } else // reuse result of last iteration { // out was modified in last iteration => restore it - ImProcFunctions::fftw_convol_blur2(out, out, W_L, H_L,sqrtf(SQR(mulradiusfftw * RetinexScales[scale]) - SQR(mulradiusfftw * RetinexScales[scale + 1])), 0); + ImProcFunctions::fftw_convol_blur2(out, out, W_L, H_L,sqrtf(SQR(mulradiusfftw * RetinexScales[scale]) - SQR(mulradiusfftw * RetinexScales[scale + 1])), 0, 0); } }