diff --git a/rtengine/ipgrain.cc b/rtengine/ipgrain.cc
new file mode 100644
index 000000000..75c9218b7
--- /dev/null
+++ b/rtengine/ipgrain.cc
@@ -0,0 +1,374 @@
+/* -*- C++ -*-
+ *
+ *  This file is part of RawTherapee.
+ *
+ *  Copyright (c) 2018 Alberto Griggio <alberto.griggio@gmail.com>
+ *  Small adaptation to Rawtherapee Locallab October 2019
+ *  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/>.
+ */
+
+/* film grain emulation.
+ * Ported from darktable (src/iop/grain.c). Original copyright/license follows
+ */
+/*
+    This file is part of darktable,
+    copyright (c) 2010-2012 Henrik Andersson.
+
+    darktable 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.
+
+    darktable 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 darktable.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#include "improcfun.h"
+#include "rt_math.h"
+
+
+namespace rtengine {
+
+namespace {
+
+constexpr float GRAIN_LIGHTNESS_STRENGTH_SCALE = 0.15f;
+constexpr float GRAIN_SCALE_FACTOR = 213.2f;
+
+constexpr int GRAIN_LUT_SIZE = 128;
+constexpr float GRAIN_LUT_DELTA_MAX = 2.0f;
+constexpr float GRAIN_LUT_DELTA_MIN = 0.0001f;
+constexpr float GRAIN_LUT_PAPER_GAMMA = 1.0f;
+
+
+const int grad3[12][3] = { { 1, 1, 0 },
+                           { -1, 1, 0 },
+                           { 1, -1, 0 },
+                           { -1, -1, 0 },
+                           { 1, 0, 1 },
+                           { -1, 0, 1 },
+                           { 1, 0, -1 },
+                           { -1, 0, -1 },
+                           { 0, 1, 1 },
+                           { 0, -1, 1 },
+                           { 0, 1, -1 },
+                           { 0, -1, -1 } };
+
+const int permutation[]
+    = { 151, 160, 137, 91,  90,  15,  131, 13,  201, 95,  96,  53,  194, 233, 7,   225, 140, 36,  103, 30,
+        69,  142, 8,   99,  37,  240, 21,  10,  23,  190, 6,   148, 247, 120, 234, 75,  0,   26,  197, 62,
+        94,  252, 219, 203, 117, 35,  11,  32,  57,  177, 33,  88,  237, 149, 56,  87,  174, 20,  125, 136,
+        171, 168, 68,  175, 74,  165, 71,  134, 139, 48,  27,  166, 77,  146, 158, 231, 83,  111, 229, 122,
+        60,  211, 133, 230, 220, 105, 92,  41,  55,  46,  245, 40,  244, 102, 143, 54,  65,  25,  63,  161,
+        1,   216, 80,  73,  209, 76,  132, 187, 208, 89,  18,  169, 200, 196, 135, 130, 116, 188, 159, 86,
+        164, 100, 109, 198, 173, 186, 3,   64,  52,  217, 226, 250, 124, 123, 5,   202, 38,  147, 118, 126,
+        255, 82,  85,  212, 207, 206, 59,  227, 47,  16,  58,  17,  182, 189, 28,  42,  223, 183, 170, 213,
+        119, 248, 152, 2,   44,  154, 163, 70,  221, 153, 101, 155, 167, 43,  172, 9,   129, 22,  39,  253,
+        19,  98,  108, 110, 79,  113, 224, 232, 178, 185, 112, 104, 218, 246, 97,  228, 251, 34,  242, 193,
+        238, 210, 144, 12,  191, 179, 162, 241, 81,  51,  145, 235, 249, 14,  239, 107, 49,  192, 214, 31,
+        181, 199, 106, 157, 184, 84,  204, 176, 115, 121, 50,  45,  127, 4,   150, 254, 138, 236, 205, 93,
+        222, 114, 67,  29,  24,  72,  243, 141, 128, 195, 78,  66,  215, 61,  156, 180 };
+
+
+class GrainEvaluator {
+public:
+    GrainEvaluator(int offset_x, int offset_y, int full_width, int full_height, double scale):
+        ox(offset_x),
+        oy(offset_y),
+        fw(full_width),
+        fh(full_height),
+        scale(scale)
+    {
+        simplex_noise_init();
+        constexpr float mb = 100.f;
+        evaluate_grain_lut(mb);
+    }
+    
+    void operator()(int isogr, int strengr, int scalegr, Imagefloat *lab, bool multithread)
+    {
+        const double strength = (strengr / 100.0);
+        const double octaves = 3;
+        const double wd = std::min(fw, fh);
+        const double zoom = (1.0 + 8 * (double(isogr) / GRAIN_SCALE_FACTOR) / 100.0) / 800.0;
+        const double s = std::max(scale / 3.0, 1.0) / (double(std::max(scalegr, 1)) / 100.0);
+
+        const int W = lab->getWidth();
+        const int H = lab->getHeight();
+        float **lab_L = lab->g.ptrs;
+
+#ifdef _OPENMP
+#       pragma omp parallel for if (multithread)
+#endif
+        for (int j = 0; j < H; ++j) {
+            double wy = oy + j;
+            double y = wy / wd;
+            for (int i = 0; i < W; ++i) {
+                double wx = ox + i;
+                double x = wx / wd;
+                double noise = simplex_2d_noise(x, y, octaves, 1.0, zoom) / s;
+                lab_L[j][i] += lut_lookup(noise * strength * GRAIN_LIGHTNESS_STRENGTH_SCALE, lab_L[j][i] / 32768.f);
+            }
+        }
+    }
+
+private:
+    void simplex_noise_init()
+    {
+        for(int i = 0; i < 512; i++) perm[i] = permutation[i & 255];
+    }
+
+    double dot(const int *g, double x, double y, double z)
+    {
+        return g[0] * x + g[1] * y + g[2] * z;
+    }
+
+    float FASTFLOOR(float x)
+    {
+        return (x > 0 ? (int)(x) : (int)(x)-1);
+    }
+
+    double simplex_noise(double xin, double yin, double zin)
+    {
+        double n0, n1, n2, n3; // Noise contributions from the four corners
+        // Skew the input space to determine which simplex cell we're in
+        const double F3 = 1.0 / 3.0;
+        const double s = (xin + yin + zin) * F3; // Very nice and simple skew factor for 3D
+        const int i = FASTFLOOR(xin + s);
+        const int j = FASTFLOOR(yin + s);
+        const int k = FASTFLOOR(zin + s);
+        const double G3 = 1.0 / 6.0; // Very nice and simple unskew factor, too
+        const double t = (i + j + k) * G3;
+        const double X0 = i - t; // Unskew the cell origin back to (x,y,z) space
+        const double Y0 = j - t;
+        const double Z0 = k - t;
+        const double x0 = xin - X0; // The x,y,z distances from the cell origin
+        const double y0 = yin - Y0;
+        const double z0 = zin - Z0;
+        // For the 3D case, the simplex shape is a slightly irregular tetrahedron.
+        // Determine which simplex we are in.
+        int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords
+        int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords
+        if(x0 >= y0)
+        {
+            if(y0 >= z0)
+            {
+                i1 = 1; // X Y Z order
+                j1 = 0;
+                k1 = 0;
+                i2 = 1;
+                j2 = 1;
+                k2 = 0;
+            }
+            else if(x0 >= z0)
+            {
+                i1 = 1; // X Z Y order
+                j1 = 0;
+                k1 = 0;
+                i2 = 1;
+                j2 = 0;
+                k2 = 1;
+            }
+            else
+            {
+                i1 = 0; // Z X Y order
+                j1 = 0;
+                k1 = 1;
+                i2 = 1;
+                j2 = 0;
+                k2 = 1;
+            }
+        }
+        else // x0<y0
+        {
+            if(y0 < z0)
+            {
+                i1 = 0; // Z Y X order
+                j1 = 0;
+                k1 = 1;
+                i2 = 0;
+                j2 = 1;
+                k2 = 1;
+            }
+            else if(x0 < z0)
+            {
+                i1 = 0; // Y Z X order
+                j1 = 1;
+                k1 = 0;
+                i2 = 0;
+                j2 = 1;
+                k2 = 1;
+            }
+            else
+            {
+                i1 = 0; // Y X Z order
+                j1 = 1;
+                k1 = 0;
+                i2 = 1;
+                j2 = 1;
+                k2 = 0;
+            }
+        }
+        //  A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z),
+        //  a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and
+        //  a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where
+        //  c = 1/6.
+        const double x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z) coords
+        const double y1 = y0 - j1 + G3;
+        const double z1 = z0 - k1 + G3;
+        const double x2 = x0 - i2 + 2.0 * G3; // Offsets for third corner in (x,y,z) coords
+        const double y2 = y0 - j2 + 2.0 * G3;
+        const double z2 = z0 - k2 + 2.0 * G3;
+        const double x3 = x0 - 1.0 + 3.0 * G3; // Offsets for last corner in (x,y,z) coords
+        const double y3 = y0 - 1.0 + 3.0 * G3;
+        const double z3 = z0 - 1.0 + 3.0 * G3;
+        // Work out the hashed gradient indices of the four simplex corners
+        const int ii = i & 255;
+        const int jj = j & 255;
+        const int kk = k & 255;
+        const int gi0 = perm[ii + perm[jj + perm[kk]]] % 12;
+        const int gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]] % 12;
+        const int gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]] % 12;
+        const int gi3 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]] % 12;
+        // Calculate the contribution from the four corners
+        double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0;
+        if(t0 < 0)
+            n0 = 0.0;
+        else
+        {
+            t0 *= t0;
+            n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0);
+        }
+        double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1;
+        if(t1 < 0)
+            n1 = 0.0;
+        else
+        {
+            t1 *= t1;
+            n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1);
+        }
+        double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2;
+        if(t2 < 0)
+            n2 = 0.0;
+        else
+        {
+            t2 *= t2;
+            n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2);
+        }
+        double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3;
+        if(t3 < 0)
+            n3 = 0.0;
+        else
+        {
+            t3 *= t3;
+            n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3);
+        }
+        // Add contributions from each corner to get the final noise value.
+        // The result is scaled to stay just inside [-1,1]
+        return 32.0 * (n0 + n1 + n2 + n3);
+    }
+
+    double simplex_2d_noise(double x, double y, uint32_t octaves, double persistance, double z)
+    {
+        double total = 0;
+
+        // parametrization of octaves to match power spectrum of real grain scans
+        static double f[] = {0.4910, 0.9441, 1.7280};
+        static double a[] = {0.2340, 0.7850, 1.2150};
+
+        for(uint32_t o = 0; o < octaves; o++)
+        {
+            total += (simplex_noise(x * f[o] / z, y * f[o] / z, o) * a[o]);
+        }
+        return total;
+    }
+
+    float paper_resp(float exposure, float mb, float gp)
+    {
+        const float delta = GRAIN_LUT_DELTA_MAX * expf((mb / 100.0f) * logf(GRAIN_LUT_DELTA_MIN));
+        const float density = (1.0f + 2.0f * delta) / (1.0f + expf( (4.0f * gp * (0.5f - exposure)) / (1.0f + 2.0f * delta) )) - delta;
+        return density;
+    }
+
+    float paper_resp_inverse(float density, float mb, float gp)
+    {
+        const float delta = GRAIN_LUT_DELTA_MAX * expf((mb / 100.0f) * logf(GRAIN_LUT_DELTA_MIN));
+        const float exposure = -logf((1.0f + 2.0f * delta) / (density + delta) - 1.0f) * (1.0f + 2.0f * delta) / (4.0f * gp) + 0.5f;
+        return exposure;
+    }
+
+    void evaluate_grain_lut(const float mb)
+    {
+        for(int i = 0; i < GRAIN_LUT_SIZE; i++)
+        {
+            for(int j = 0; j < GRAIN_LUT_SIZE; j++)
+            {
+                const float gu = (float)i / (GRAIN_LUT_SIZE - 1) - 0.5;
+                const float l = (float)j / (GRAIN_LUT_SIZE - 1);
+                grain_lut[j * GRAIN_LUT_SIZE + i] = 32768.f * (paper_resp(gu + paper_resp_inverse(l, mb, GRAIN_LUT_PAPER_GAMMA), mb, GRAIN_LUT_PAPER_GAMMA) - l);
+            }
+        }
+    }
+
+    float lut_lookup(const float x, const float y)
+    {
+        const float _x = LIM((x + 0.5f) * (GRAIN_LUT_SIZE - 1), 0.f, float(GRAIN_LUT_SIZE - 1));
+        const float _y = LIM(y * (GRAIN_LUT_SIZE - 1), 0.f, float(GRAIN_LUT_SIZE - 1));
+
+        const int _x0 = _x < GRAIN_LUT_SIZE - 2 ? _x : GRAIN_LUT_SIZE - 2;
+        const int _y0 = _y < GRAIN_LUT_SIZE - 2 ? _y : GRAIN_LUT_SIZE - 2;
+
+        const int _x1 = _x0 + 1;
+        const int _y1 = _y0 + 1;
+
+        const float x_diff = _x - _x0;
+        const float y_diff = _y - _y0;
+
+        const float l00 = grain_lut[_y0 * GRAIN_LUT_SIZE + _x0];
+        const float l01 = grain_lut[_y0 * GRAIN_LUT_SIZE + _x1];
+        const float l10 = grain_lut[_y1 * GRAIN_LUT_SIZE + _x0];
+        const float l11 = grain_lut[_y1 * GRAIN_LUT_SIZE + _x1];
+
+        const float xy0 = (1.0 - y_diff) * l00 + l10 * y_diff;
+        const float xy1 = (1.0 - y_diff) * l01 + l11 * y_diff;
+        return xy0 * (1.0f - x_diff) + xy1 * x_diff;
+    }
+
+    int ox;
+    int oy;
+    int fw;
+    int fh;
+    double scale;
+    int perm[512];
+    float grain_lut[GRAIN_LUT_SIZE * GRAIN_LUT_SIZE];
+};
+
+} // namespace
+
+
+void ImProcFunctions::filmGrain(Imagefloat *rgb, int isogr, int strengr, int scalegr, int bfw, int bfh)
+{
+
+    GrainEvaluator ge(0, 0, bfw, bfh, scale);
+    ge(isogr, strengr, scalegr, rgb, multiThread);
+}
+
+} // namespace rtengine