From dc4788e00b86d3c40796198e0a1941b6b12ae275 Mon Sep 17 00:00:00 2001
From: Desmis <jdesmis@gmail.com>
Date: Mon, 1 Jul 2019 10:17:07 +0200
Subject: [PATCH] another similar improvment

---
 rtengine/iplocallab.cc |  5 +++--
 rtengine/ipretinex.cc  | 29 +++++++++++++++++++++--------
 2 files changed, 24 insertions(+), 10 deletions(-)

diff --git a/rtengine/iplocallab.cc b/rtengine/iplocallab.cc
index 0195e48f8..3674613ec 100644
--- a/rtengine/iplocallab.cc
+++ b/rtengine/iplocallab.cc
@@ -4009,8 +4009,9 @@ void ImProcFunctions::fftw_convol_blur(float *input, float *output, int bfw, int
     ** radius = sigma for kernel
     ** n_x n_y relative width and high for kernel
     ** Gaussian blur is given by G(x,y) = (1/2*PI*sigma) * exp(-(x2 + y2) / 2* sigma2)
-    ** its traduction in Fourier transform is G(x,y) =  exp((-sigma)*(PI * x2 + PI * y2))
-    ** after several test the only result that works very well is with fftkern = 0 and algo = 0
+    ** its traduction in Fourier transform is G(x,y) =  exp((-sigma)*(PI * x2 + PI * y2)), for some authors it is not sigma but sigma^2..I have tried...huge diffrences with Gaussianblur
+    ** after several test the only result that works very well is with fftkern = 0 and algo = 0, and as there is big differences with Gaussianblur, I put an empirical correction in Ipretinex and Iplocalcontrast
+    ** in fact no importance....if it is this function (for sigma) or another... we are not in research :)
 */
     BENCHFUN
   
diff --git a/rtengine/ipretinex.cc b/rtengine/ipretinex.cc
index 356c58d66..cfd03b2b5 100644
--- a/rtengine/ipretinex.cc
+++ b/rtengine/ipretinex.cc
@@ -962,19 +962,32 @@ void ImProcFunctions::MSRLocal(int sp, bool fftw, int lum, LabImage * bufreti, L
         }
 
         float *buffer = new float[W_L * H_L];
-        float kr = 1.f;
+        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
-        //above 400 ==> 1.f
+            //emprical adjustement between FFTW radius and Gaussainblur
+            //under 50 ==> 10.f
+            // 400 ==> 1.f
             float sigm = RetinexScales[scale];
-            float mulradiusfftw = 1.f;
             float ak = -9.f / 350.f;
             float bk = 10.f - 50.f * ak;
             kr = ak * sigm + bk;
             if(sigm < 50.f) kr = 10.f;
-            if(sigm > 400.f) kr = 1.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.5f / 4600.f;
+                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;
+                }
+                
+            }
            
         if(!fftw) {
 #ifdef _OPENMP
@@ -984,11 +997,11 @@ void ImProcFunctions::MSRLocal(int sp, bool fftw, int lum, LabImage * bufreti, L
 
                 if (scale == scal - 1)
                 {
-                    gaussianBlur(src, out, W_L, H_L, RetinexScales[scale], buffer);
+                    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(RetinexScales[scale]) - SQR(RetinexScales[scale + 1])), buffer);
+                    gaussianBlur(out, out, W_L, H_L, sqrtf(SQR(kg * RetinexScales[scale]) - SQR(kg * RetinexScales[scale + 1])), buffer);
                 }
             }
         } else {