diff --git a/rtengine/array2D.h b/rtengine/array2D.h
index 25d644c85..48a789bf8 100644
--- a/rtengine/array2D.h
+++ b/rtengine/array2D.h
@@ -267,11 +267,11 @@ public:
         ar_realloc(w, h);
         memcpy(data, copy, w * h * sizeof(T));
     }
-    int width()
+    int width() const
     {
         return x;
     }
-    int height()
+    int height() const
     {
         return y;
     }
diff --git a/rtengine/rescale.h b/rtengine/rescale.h
new file mode 100644
index 000000000..ba9a01c99
--- /dev/null
+++ b/rtengine/rescale.h
@@ -0,0 +1,99 @@
+/* -*- 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 {
+
+inline float getBilinearValue(const array2D<float> &src, float x, float y)
+{
+    const int W = src.width();
+    const int H = src.height();
+    
+    // Get integer and fractional parts of numbers
+    int xi = x;
+    int yi = y;
+    float xf = x - xi;
+    float yf = y - yi;
+    int xi1 = std::min(xi + 1, W - 1);
+    int yi1 = std::min(yi + 1, H - 1);
+
+    float bl = src[yi][xi];
+    float br = src[yi][xi1];
+    float tl = src[yi1][xi];
+    float tr = src[yi1][xi1];
+
+    // interpolate
+    float b = xf * br + (1.f - xf) * bl;
+    float t = xf * tr + (1.f - xf) * tl;
+    float pxf = yf * t + (1.f - yf) * b;
+    return pxf;
+}
+
+
+inline void rescaleBilinear(const array2D<float> &src, array2D<float> &dst, bool multithread)
+{
+    const int Ws = src.width();
+    const int Hs = src.height();
+    const int Wd = dst.width();
+    const int Hd = dst.height();
+    
+    float col_scale = float (Ws) / float (Wd);
+    float row_scale = float (Hs) / float (Hd);
+
+#ifdef _OPENMP
+    #pragma omp parallel for if (multithread)
+#endif
+
+    for (int y = 0; y < Hd; ++y) {
+        float ymrs = y * row_scale;
+
+        for (int x = 0; x < Wd; ++x) {
+            dst[y][x] = getBilinearValue(src, x * col_scale, ymrs);
+        }
+    }
+}
+
+
+inline void rescaleNearest(const array2D<float> &src, array2D<float> &dst, bool multithread)
+{
+    const int width = src.width();
+    const int height = src.height();
+    const int nw = dst.width();
+    const int nh = dst.height();
+
+#ifdef _OPENMP
+    #pragma omp parallel for if (multithread)
+#endif
+
+    for (int y = 0; y < nh; ++y) {
+        int sy = y * height / nh;
+
+        for (int x = 0; x < nw; ++x) {
+            int sx = x * width / nw;
+            dst[y][x] = src[sy][sx];
+        }
+    }
+}
+
+
+} // namespace rtengine
diff --git a/rtengine/tmo_fattal02.cc b/rtengine/tmo_fattal02.cc
index 359b33575..dc7826501 100644
--- a/rtengine/tmo_fattal02.cc
+++ b/rtengine/tmo_fattal02.cc
@@ -73,6 +73,7 @@
 #include "sleef.c"
 #include "opthelper.h"
 #include "rt_algo.h"
+#include "rescale.h"
 
 namespace rtengine
 {
@@ -938,66 +939,14 @@ void solve_pde_fft (Array2Df *F, Array2Df *U, Array2Df *buf, bool multithread)/*
  * RT code from here on
  *****************************************************************************/
 
-inline float get_bilinear_value (const Array2Df &src, float x, float y)
+inline void rescale_bilinear (const Array2Df &src, Array2Df &dst, bool multithread)
 {
-    // Get integer and fractional parts of numbers
-    int xi = x;
-    int yi = y;
-    float xf = x - xi;
-    float yf = y - yi;
-    int xi1 = std::min (xi + 1, src.getCols() - 1);
-    int yi1 = std::min (yi + 1, src.getRows() - 1);
-
-    float bl = src (xi, yi);
-    float br = src (xi1, yi);
-    float tl = src (xi, yi1);
-    float tr = src (xi1, yi1);
-
-    // interpolate
-    float b = xf * br + (1.f - xf) * bl;
-    float t = xf * tr + (1.f - xf) * tl;
-    float pxf = yf * t + (1.f - yf) * b;
-    return pxf;
+    rescaleBilinear(src, dst, multithread);
 }
 
-
-void rescale_bilinear (const Array2Df &src, Array2Df &dst, bool multithread)
+inline void rescale_nearest (const Array2Df &src, Array2Df &dst, bool multithread)
 {
-    float col_scale = float (src.getCols()) / float (dst.getCols());
-    float row_scale = float (src.getRows()) / float (dst.getRows());
-
-#ifdef _OPENMP
-    #pragma omp parallel for if (multithread)
-#endif
-
-    for (int y = 0; y < dst.getRows(); ++y) {
-        float ymrs = y * row_scale;
-
-        for (int x = 0; x < dst.getCols(); ++x) {
-            dst (x, y) = get_bilinear_value (src, x * col_scale, ymrs);
-        }
-    }
-}
-
-void rescale_nearest (const Array2Df &src, Array2Df &dst, bool multithread)
-{
-    const int width = src.getCols();
-    const int height = src.getRows();
-    const int nw = dst.getCols();
-    const int nh = dst.getRows();
-
-#ifdef _OPENMP
-    #pragma omp parallel for if (multithread)
-#endif
-
-    for (int y = 0; y < nh; ++y) {
-        int sy = y * height / nh;
-
-        for (int x = 0; x < nw; ++x) {
-            int sx = x * width / nw;
-            dst (x, y) = src (sx, sy);
-        }
-    }
+    rescaleNearest(src, dst, multithread);
 }