Fix blur radius when using the Fourier method

Fix the formula, adapt old pp3s, and change the maximum blur radius
where applicable.

rtengine/iplocalcontrast.cc

@@ -52,21 +52,8 @@ void ImProcFunctions::localContrast(LabImage *lab, float **destination, const rt
 #endif
         gaussianBlur(lab->L, buf, width, height, sigma);
     } else {
-        float kr = 1.f;
-        //emprical adjustment between FFTW radius and Gaussainblur
-        //under 50 ==> 10.f
-        //above 400 ==> 1.f
-        if(settings->fftwsigma == false) {//empirical formula
-            float ak = -9.f / 350.f;
-            float bk = 10.f - 50.f * ak;
-            kr = ak * sigma + bk;
-            if(sigma < 50.f) kr = 10.f;
-            if(sigma > 400.f) kr = 1.f;
-        } else {//sigma *= sigma
-            kr = sigma;
-        }
         //OPENMP disabled
-        ImProcFunctions::fftw_convol_blur2(lab->L, buf, width, height, kr * sigma, 0, 0);
+        ImProcFunctions::fftw_convol_blur2(lab->L, buf, width, height, sigma, 0, 0);
     }
 #ifdef _OPENMP
     #pragma omp parallel for if(multiThread)

rtengine/iplocallab.cc

@@ -9523,8 +9523,8 @@ void ImProcFunctions::fftw_convol_blur(float * input, float * output, int bfw, i
 
     /*define the gaussian constants for the convolution kernel*/
     if (algo == 0) {
-        n_x = rtengine::RT_PI / (double) bfw; //ipol
+        n_x = rtengine::RT_PI / (double) bfw / std::sqrt(2.); //ipol
-        n_y = rtengine::RT_PI / (double) bfh;
+        n_y = rtengine::RT_PI / (double) bfh / std::sqrt(2.);
     } else if (algo == 1) {
         n_x = 1.f / bfw; //gauss
         n_y = 1.f / bfh;
@@ -9547,7 +9547,7 @@ void ImProcFunctions::fftw_convol_blur(float * input, float * output, int bfw, i
 
             for (int i = 0; i < bfw; i++)
                 if (algo == 0) {
-                    kern[ i + index] = exp((float)(-radius) * (n_x * i * i + n_y * j * j)); //calculate Gauss kernel Ipol formula
+                    kern[ i + index] = exp((float)(-radius * 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
                 }
@@ -9585,7 +9585,7 @@ void ImProcFunctions::fftw_convol_blur(float * input, float * output, int bfw, i
                 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 without FFT - some authors says radius*radius but differences with Gaussianblur
+                    out[i + index] *= exp((float)(-radius * radius) * (n_x * i * i + n_y * j * j));    //apply Gauss kernel without FFT - some authors says radius*radius but differences with Gaussianblur
                 }
             }
         } else if (algo == 1) {

rtengine/ipretinex.cc

@@ -1314,14 +1314,26 @@ void ImProcFunctions::MSRLocal(int call, int sp, bool fftw, int lum, float** red
                 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);
+                    // FFT blur radius fixed in 5.12, resulting in different
+                    // blur amount. Here we preserve the original behavior by
+                    // multiplying the original sigma SQR(RetinexScales[scale])
+                    // with 2 then taking the square root.
+                    const auto gaussianSigma = std::sqrt(2.f) * RetinexScales[scale];
+                    ImProcFunctions::fftw_convol_blur2(src, out, bfwr, bfhr, gaussianSigma, 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);
+                    // FFT blur radius fixed in 5.12, resulting in different
+                    // blur amount. Here we preserve the original behavior by
+                    // multiplying the original sigma
+                    // SQR(RetinexScales[scale]) - SQR(RetinexScales[scale + 1])
+                    // with 2 then taking the square root.
+                    const auto gaussianSigma =
+                        std::sqrt(2 * (SQR(RetinexScales[scale]) - SQR(RetinexScales[scale + 1])));
+                    ImProcFunctions::fftw_convol_blur2(out, out, bfwr, bfhr, gaussianSigma, 0, 0);
                 }
             }
         }

rtengine/procparams.cc

@@ -10084,6 +10084,40 @@ int ProcParams::load(const Glib::ustring& fname, ParamsEdited* pedited)
                 assignFromKeyfile(keyFile, "Locallab", "shadmaskcie_" + index_str, spot.shadmaskcie, spotEdited.shadmaskcie);
                 assignFromKeyfile(keyFile, "Locallab", "LLmaskcieCurvewav_" + index_str, spot.LLmaskciecurvewav, spotEdited.LLmaskciecurvewav);
 
+                if (ppVersion < 352) {
+                    // FFT Gaussian blur fixed in 5.12. Converts old radii to
+                    // the equivalent new radii if FFT mode is enabled. The
+                    // fixed blur radius has a factor of radius / 2 compared to
+                    // the original, so the base conversion is sqrt(2 * radius).
+                    // Derivation: old_radius = new_radius * new_radius / 2
+                    //             2 * old_radius = new_radius * new_radius
+                    //             sqrt(2 * old_radius) = new_radius
+                    if (spot.fftColorMask) {
+                        spot.blurcol = std::sqrt(2 * spot.blurcol);
+                        spotEdited.blurcol = true;
+                    }
+                    if (spot.fftwbl || spot.radius > 30.) {
+                        // Internally, the radius is divided by 1.4, so the
+                        // factor is actually radius / 1.4 / 2.
+                        spot.radius = std::sqrt(2.8 * spot.radius);
+                        spotEdited.radius = true;
+                    }
+                    if (spot.fftwlc) {
+                        // Internally, the radius was squared, so replace
+                        // old_radius with old_radius^2 before solving.
+                        spot.lcradius = std::sqrt(2.) * spot.lcradius;
+                        spotEdited.lcradius = true;
+                    }
+                    if (spot.fftmask) {
+                        spot.blurmask = std::sqrt(2 * spot.blurmask);
+                        spotEdited.blurmask = true;
+                    }
+                    if (spot.fftcieMask) {
+                        spot.blurcie = std::sqrt(2 * spot.blurcie);
+                        spotEdited.blurcie = true;
+                    }
+                }
+
                 if (keyFile.has_key("Locallab", "CSThresholdcie_" + index_str)) {
                     const std::vector<int> thresh = keyFile.get_integer_list("Locallab", "CSThresholdcie_" + index_str);
 

rtgui/locallabtools.cc

@@ -28,7 +28,7 @@ fft *
 #include "../rtengine/color.h"
 
 #define MINRAD 1.5
-#define MAXRAD 10000
+#define MAXRAD 1000
 #define CENTERRAD 100
 #define MINCHRO 0.
 #define MAXCHRO 150.

rtgui/locallabtools2.cc

@@ -4660,7 +4660,7 @@ void LocallabContrast::updateContrastGUI3()
     const double temp = lcradius->getValue();
 
     if (fftwlc->get_active()) {
-        lcradius->setLimits(20, 1000, 1, 80);
+        lcradius->setLimits(20, 1500, 1, 80);
     } else {
         lcradius->setLimits(20, 100, 1, 80);
     }

rtgui/ppversion.h

@@ -1,11 +1,13 @@
 #pragma once
 
 // This number has to be incremented whenever the PP3 file format is modified or the behaviour of a tool changes
-#define PPVERSION 351
+#define PPVERSION 352
 #define PPVERSION_AEXP 301 //value of PPVERSION when auto exposure algorithm was modified
 
 /*
   Log of version changes
+   352  2024-10-27
+        FFT Gaussian blur radius fixed
    351  2024-06-19
         take into account Global in selective editing
    350  2023-03-05