Added an implementation of the Fast Guided Filter

derived directly from the pseudo-code of the paper:

 Fast Guided Filter
 by Kaiming He, Jian Sun

available at https://arxiv.org/abs/1505.00996

rtengine/CMakeLists.txt

@@ -127,6 +127,7 @@ set(RTENGINESOURCEFILES
     xtrans_demosaic.cc
     vng4_demosaic_RT.cc
     ipsoftlight.cc
+    guidedfilter.cc
     )
 
 if(LENSFUN_HAS_LOAD_DIRECTORY)

rtengine/boxblur.h

@@ -35,6 +35,8 @@ namespace rtengine
 template<class T, class A> void boxblur (T** src, A** dst, int radx, int rady, int W, int H)
 {
     //box blur image; box range = (radx,rady)
+    radx = min(radx, W-1);
+    rady = min(rady, H-1);
 
     AlignedBuffer<float>* buffer = new AlignedBuffer<float> (W * H);
     float* temp = buffer->data;

rtengine/guidedfilter.cc

@@ -0,0 +1,194 @@
+/* -*- C++ -*-
+ *
+ *  This file is part of RawTherapee.
+ *
+ *  Copyright (c) 2018 Alberto Griggio <alberto.griggio@gmail.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 <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * This is a Fast Guided Filter implementation, derived directly from the
+ * pseudo-code of the paper:
+ *
+ * Fast Guided Filter
+ * by Kaiming He, Jian Sun
+ *
+ * available at https://arxiv.org/abs/1505.00996
+ */
+
+#include "guidedfilter.h"
+#include "boxblur.h"
+#include "rescale.h"
+#include "imagefloat.h"
+
+namespace rtengine {
+
+#if 0
+#  define DEBUG_DUMP(arr)                                                 \
+    do {                                                                \
+        Imagefloat im(arr.width(), arr.height());                      \
+        const char *out = "/tmp/" #arr ".tif";                     \
+        for (int y = 0; y < im.getHeight(); ++y) {                      \
+            for (int x = 0; x < im.getWidth(); ++x) {                   \
+                im.r(y, x) = im.g(y, x) = im.b(y, x) = arr[y][x] * 65535.f; \
+            }                                                           \
+        }                                                               \
+        im.saveTIFF(out, 16);                                           \
+    } while (false)
+#else
+#  define DEBUG_DUMP(arr)
+#endif
+
+
+void guidedFilter(const array2D<float> &guide, const array2D<float> &src, array2D<float> &dst, int r, float epsilon, bool multithread, int subsampling)
+{
+    const int W = src.width();
+    const int H = src.height();
+
+    enum Op { MUL, DIVEPSILON, ADD, SUB, ADDMUL, SUBMUL };
+
+    const auto apply =
+        [=](Op op, array2D<float> &res, const array2D<float> &a, const array2D<float> &b, const array2D<float> &c=array2D<float>()) -> void
+        {
+            const int w = res.width();
+            const int h = res.height();
+            
+#ifdef _OPENMP
+            #pragma omp parallel for if (multithread)
+#endif
+            for (int y = 0; y < h; ++y) {
+                for (int x = 0; x < w; ++x) {
+                    float r;
+                    float aa = a[y][x];
+                    float bb = b[y][x];
+                    switch (op) {
+                    case MUL:
+                        r = aa * bb;
+                        break;
+                    case DIVEPSILON:
+                        r = aa / (bb + epsilon);
+                        break;
+                    case ADD:
+                        r = aa + bb;
+                        break;
+                    case SUB:
+                        r = aa - bb;
+                        break;
+                    case ADDMUL:
+                        r = aa * bb + c[y][x];
+                        break;
+                    case SUBMUL:
+                        r = c[y][x] - (aa * bb);
+                        break;
+                    default:
+                        assert(false);
+                        r = 0;
+                        break;
+                    }
+                    res[y][x] = r;
+                }
+            }
+        };
+
+    // use the terminology of the paper (Algorithm 2)
+    const array2D<float> &I = guide;
+    const array2D<float> &p = src;
+    array2D<float> &q = dst;
+
+    const auto f_mean =
+        [](array2D<float> &d, array2D<float> &s, int rad) -> void
+        {
+            boxblur<float, float>(s, d, rad, rad, s.width(), s.height());
+        };
+
+    const auto f_subsample =
+        [=](array2D<float> &d, const array2D<float> &s) -> void
+        {
+            rescaleBilinear(s, d, multithread);
+        };
+
+    const auto f_upsample = f_subsample;
+    
+    const int w = W / subsampling;
+    const int h = H / subsampling;
+
+    array2D<float> I1(w, h);
+    array2D<float> p1(w, h);
+
+    f_subsample(I1, I);
+    f_subsample(p1, p);
+
+    DEBUG_DUMP(I);
+    DEBUG_DUMP(p);
+    DEBUG_DUMP(I1);
+    DEBUG_DUMP(p1);
+
+    float r1 = float(r) / subsampling;
+
+    array2D<float> meanI(w, h);
+    f_mean(meanI, I1, r1);
+    DEBUG_DUMP(meanI);
+
+    array2D<float> meanp(w, h);
+    f_mean(meanp, p1, r1);
+    DEBUG_DUMP(meanp);
+
+    array2D<float> &corrIp = p1;
+    apply(MUL, corrIp, I1, p1);
+    f_mean(corrIp, corrIp, r1);
+    DEBUG_DUMP(corrIp);
+
+    array2D<float> &corrI = I1;
+    apply(MUL, corrI, I1, I1);
+    f_mean(corrI, corrI, r1);
+    DEBUG_DUMP(corrI);
+
+    array2D<float> &varI = corrI;
+    apply(SUBMUL, varI, meanI, meanI, corrI);
+    DEBUG_DUMP(varI);
+
+    array2D<float> &covIp = corrIp;
+    apply(SUBMUL, covIp, meanI, meanp, corrIp);
+    DEBUG_DUMP(covIp);
+
+    array2D<float> &a = varI;
+    apply(DIVEPSILON, a, covIp, varI);
+    DEBUG_DUMP(a);
+
+    array2D<float> &b = covIp;
+    apply(SUBMUL, b, a, meanI, meanp);
+    DEBUG_DUMP(b);
+
+    array2D<float> &meana = a;
+    f_mean(meana, a, r1);
+    DEBUG_DUMP(meana);
+
+    array2D<float> &meanb = b;
+    f_mean(meanb, b, r1);
+    DEBUG_DUMP(meanb);
+
+    array2D<float> meanA(W, H);
+    f_upsample(meanA, meana);
+    DEBUG_DUMP(meanA);
+
+    array2D<float> &meanB = q;
+    f_upsample(meanB, meanb);
+    DEBUG_DUMP(meanB);
+
+    apply(ADDMUL, q, meanA, I, meanB);
+    DEBUG_DUMP(q);
+}
+
+} // namespace rtengine

rtengine/guidedfilter.h

@@ -0,0 +1,29 @@
+/* -*- C++ -*-
+ *
+ *  This file is part of RawTherapee.
+ *
+ *  Copyright (c) 2018 Alberto Griggio <alberto.griggio@gmail.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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "array2D.h"
+
+namespace rtengine {
+
+void guidedFilter(const array2D<float> &guide, const array2D<float> &src, array2D<float> &dst, int r, float epsilon, bool multithread, int subsampling=4);
+
+} // namespace rtengine