Improve turbulence renderer

(bzr r9508.1.35)

3 files changed, 322 insertions, 284 deletions

diff --git a/src/display/cairo-templates.h b/src/display/cairo-templates.h
index e0b1bfd98..88f3bf41a 100644
--- a/src/display/cairo-templates.h
+++ b/src/display/cairo-templates.h
@@ -644,9 +644,11 @@ private:
 };*/
 
 // Some helpers for pixel manipulation
-
 G_GNUC_CONST inline gint32
 pxclamp(gint32 v, gint32 low, gint32 high) {
+    // NOTE: it is possible to write a "branchless" clamping operation.
+    // However, it will be slower than this function, because the code below
+    // is compiled to conditional moves.
     if (v < low) return low;
     if (v > high) return high;
     return v;
diff --git a/src/display/nr-filter-turbulence.cpp b/src/display/nr-filter-turbulence.cpp
index de46b1565..d60c4f617 100644
--- a/src/display/nr-filter-turbulence.cpp
+++ b/src/display/nr-filter-turbulence.cpp
@@ -30,45 +30,282 @@
 namespace Inkscape {
 namespace Filters{
 
-/* Produces results in the range [1, 2**31 - 2].
-Algorithm is: r = (a * r) mod m
-where a = 16807 and m = 2**31 - 1 = 2147483647
-See [Park & Miller], CACM vol. 31 no. 10 p. 1195, Oct. 1988
-To test: the algorithm should produce the result 1043618065
-as the 10,000th generated number if the original seed is 1.
-*/
-#define RAND_m 2147483647 /* 2**31 - 1 */
-#define RAND_a 16807 /* 7**5; primitive root of m */
-#define RAND_q 127773 /* m / a */
-#define RAND_r 2836 /* m % a */
-//#define BSize 0x100 // defined in the header
-#define BM 0xff
-#define PerlinN 0x1000
-#define NP 12 /* 2^PerlinN */
-#define NM 0xfff
-#define s_curve(t) ( t * t * (3. - 2. * t) )
-#define turb_lerp(t, a, b) ( a + t * (b - a) )
-
-struct StitchInfo
-{
-  int nWidth; // How much to subtract to wrap for stitching.
-  int nHeight;
-  int nWrapX; // Minimum value to wrap.
-  int nWrapY;
+class TurbulenceGenerator {
+public:
+    TurbulenceGenerator()
+        : _wrapx(0)
+        , _wrapy(0)
+        , _wrapw(0)
+        , _wraph(0)
+        , _inited(false)
+    {}
+
+    void init(long seed, Geom::Rect const &tile, Geom::Point const &freq, bool stitch,
+        bool fractalnoise, int octaves)
+    {
+        // setup random number generator
+        _setupSeed(seed);
+
+        // set values
+        _tile = tile;
+        _baseFreq = freq;
+        _stitchTiles = stitch;
+        _fractalnoise = fractalnoise;
+        _octaves = octaves;
+
+        int i;
+        for (int k = 0; k < 4; ++k) {
+            for (i = 0; i < BSize; ++i) {
+                _latticeSelector[i] = i;
+
+                _gradient[i][k][0] = static_cast<double>(_random() % (BSize*2) - BSize) / BSize;
+                _gradient[i][k][1] = static_cast<double>(_random() % (BSize*2) - BSize) / BSize;
+
+                // normalize gradient
+                double s = hypot(_gradient[i][k][0], _gradient[i][k][1]);
+                _gradient[i][k][0] /= s;
+                _gradient[i][k][1] /= s;
+            }
+        }
+        while (--i) {
+            // shuffle lattice selectors
+            int j = _random() % BSize;
+            std::swap(_latticeSelector[i], _latticeSelector[j]);
+        }
+
+        // fill out the remaining part of the gradient
+        for (i = 0; i < BSize + 2; ++i)
+        {
+            _latticeSelector[BSize + i] = _latticeSelector[i];
+
+            for(int k = 0; k < 4; ++k) {
+                _gradient[BSize + i][k][0] = _gradient[i][k][0];
+                _gradient[BSize + i][k][1] = _gradient[i][k][1];
+            }
+        }
+
+        // When stitching tiled turbulence, the frequencies must be adjusted
+        // so that the tile borders will be continuous.
+        if (_stitchTiles) {
+            if (_baseFreq[Geom::X] != 0.0)
+            {
+                double freq = _baseFreq[Geom::X];
+                double lo = floor(_tile.width() * freq) / _tile.width();
+                double hi = ceil(_tile.width() * freq) / _tile.width();
+                _baseFreq[Geom::X] = freq / lo < hi / freq ? lo : hi;
+            }
+            if (_baseFreq[Geom::Y] != 0.0)
+            {
+                double freq = _baseFreq[Geom::Y];
+                double lo = floor(_tile.height() * freq) / _tile.height();
+                double hi = ceil(_tile.height() * freq) / _tile.height();
+                _baseFreq[Geom::Y] = freq / lo < hi / freq ? lo : hi;
+            }
+
+            _wrapw = _tile.width() * _baseFreq[Geom::X] + 0.5;
+            _wraph = _tile.height() * _baseFreq[Geom::Y] + 0.5;
+            _wrapx = _tile.left() * _baseFreq[Geom::X] + PerlinOffset + _wrapw;
+            _wrapy = _tile.top() * _baseFreq[Geom::Y] + PerlinOffset + _wraph;
+        }
+        _inited = true;
+    }
+
+    G_GNUC_PURE
+    guint32 turbulencePixel(Geom::Point const &p) const {
+        int wrapx = _wrapx, wrapy = _wrapy, wrapw = _wrapw, wraph = _wraph;
+
+        double pixel[4];
+        double x = p[Geom::X] * _baseFreq[Geom::X];
+        double y = p[Geom::Y] * _baseFreq[Geom::Y];
+        double ratio = 1.0;
+
+        for (int k = 0; k < 4; ++k)
+            pixel[k] = 0.0;
+
+        for(int octave = 0; octave < _octaves; ++octave)
+        {
+            double tx = x + PerlinOffset;
+            double bx = floor(tx);
+            double rx0 = tx - bx, rx1 = rx0 - 1.0;
+            int bx0 = bx, bx1 = bx0 + 1;
+
+            double ty = y + PerlinOffset;
+            double by = floor(ty);
+            double ry0 = ty - by, ry1 = ry0 - 1.0;
+            int by0 = by, by1 = by0 + 1;
+
+            if (_stitchTiles) {
+                if (bx0 >= wrapx) bx0 -= wrapw;
+                if (bx1 >= wrapx) bx1 -= wrapw;
+                if (by0 >= wrapy) by0 -= wraph;
+                if (by1 >= wrapy) by1 -= wraph;
+            }
+            bx0 &= BMask;
+            bx1 &= BMask;
+            by0 &= BMask;
+            by1 &= BMask;
+
+            int i = _latticeSelector[bx0];
+            int j = _latticeSelector[bx1];
+            int b00 = _latticeSelector[i + by0];
+            int b01 = _latticeSelector[i + by1];
+            int b10 = _latticeSelector[j + by0];
+            int b11 = _latticeSelector[j + by1];
+
+            double sx = _scurve(rx0);
+            double sy = _scurve(ry0);
+
+            double result[4];
+            // channel numbering: R=0, G=1, B=2, A=3
+            for (int k = 0; k < 4; ++k) {
+                double const *qxa = _gradient[b00][k];
+                double const *qxb = _gradient[b10][k];
+                double a = _lerp(sx, rx0 * qxa[0] + ry0 * qxa[1],
+                                     rx1 * qxb[0] + ry0 * qxb[1]);
+                double const *qya = _gradient[b01][k];
+                double const *qyb = _gradient[b11][k];
+                double b = _lerp(sx, rx0 * qya[0] + ry1 * qya[1],
+                                     rx1 * qyb[0] + ry1 * qyb[1]);
+                result[k] = _lerp(sy, a, b);
+            }
+
+            if (_fractalnoise) {
+                for (int k = 0; k < 4; ++k)
+                    pixel[k] += result[k] / ratio;
+            } else {
+                for (int k = 0; k < 4; ++k)
+                    pixel[k] += fabs(result[k]) / ratio;
+            }
+
+            x *= 2;
+            y *= 2;
+            ratio *= 2;
+
+            if(_stitchTiles)
+            {
+                // Update stitch values. Subtracting PerlinOffset before the multiplication and
+                // adding it afterward simplifies to subtracting it once.
+                wrapw *= 2;
+                wraph *= 2;
+                wrapx = wrapx*2 - PerlinOffset;
+                wrapy = wrapy*2 - PerlinOffset;
+            }
+        }
+
+        if (_fractalnoise) {
+            guint32 r = CLAMP_D_TO_U8((pixel[0]*255.0 + 255.0) / 2);
+            guint32 g = CLAMP_D_TO_U8((pixel[1]*255.0 + 255.0) / 2);
+            guint32 b = CLAMP_D_TO_U8((pixel[2]*255.0 + 255.0) / 2);
+            guint32 a = CLAMP_D_TO_U8((pixel[3]*255.0 + 255.0) / 2);
+            r = premul_alpha(r, a);
+            g = premul_alpha(g, a);
+            b = premul_alpha(b, a);
+            ASSEMBLE_ARGB32(pxout, a,r,g,b);
+            return pxout;
+        } else {
+            guint32 r = CLAMP_D_TO_U8(pixel[0]*255.0);
+            guint32 g = CLAMP_D_TO_U8(pixel[1]*255.0);
+            guint32 b = CLAMP_D_TO_U8(pixel[2]*255.0);
+            guint32 a = CLAMP_D_TO_U8(pixel[3]*255.0);
+            r = premul_alpha(r, a);
+            g = premul_alpha(g, a);
+            b = premul_alpha(b, a);
+            ASSEMBLE_ARGB32(pxout, a,r,g,b);
+            return pxout;
+        }
+    }
+
+    //G_GNUC_PURE
+    /*guint32 turbulencePixel(Geom::Point const &p) const {
+        if (!_fractalnoise) {
+            guint32 r = CLAMP_D_TO_U8(turbulence(0, p)*255.0);
+            guint32 g = CLAMP_D_TO_U8(turbulence(1, p)*255.0);
+            guint32 b = CLAMP_D_TO_U8(turbulence(2, p)*255.0);
+            guint32 a = CLAMP_D_TO_U8(turbulence(3, p)*255.0);
+            r = premul_alpha(r, a);
+            g = premul_alpha(g, a);
+            b = premul_alpha(b, a);
+            ASSEMBLE_ARGB32(pxout, a,r,g,b);
+            return pxout;
+        } else {
+            guint32 r = CLAMP_D_TO_U8((turbulence(0, p)*255.0 + 255.0) / 2);
+            guint32 g = CLAMP_D_TO_U8((turbulence(1, p)*255.0 + 255.0) / 2);
+            guint32 b = CLAMP_D_TO_U8((turbulence(2, p)*255.0 + 255.0) / 2);
+            guint32 a = CLAMP_D_TO_U8((turbulence(3, p)*255.0 + 255.0) / 2);
+            r = premul_alpha(r, a);
+            g = premul_alpha(g, a);
+            b = premul_alpha(b, a);
+            ASSEMBLE_ARGB32(pxout, a,r,g,b);
+            return pxout;
+        }
+    }*/
+
+    bool ready() const { return _inited; }
+    void dirty() { _inited = false; }
+
+private:
+    void _setupSeed(long seed) {
+        _seed = seed;
+        if (_seed <= 0) _seed = -(_seed % (RAND_m - 1)) + 1;
+        if (_seed > RAND_m - 1) _seed = RAND_m - 1;
+    }
+    long _random() {
+        /* Produces results in the range [1, 2**31 - 2].
+         * Algorithm is: r = (a * r) mod m
+         * where a = 16807 and m = 2**31 - 1 = 2147483647
+         * See [Park & Miller], CACM vol. 31 no. 10 p. 1195, Oct. 1988
+         * To test: the algorithm should produce the result 1043618065
+         * as the 10,000th generated number if the original seed is 1. */
+        _seed = RAND_a * (_seed % RAND_q) - RAND_r * (_seed / RAND_q);
+        if (_seed <= 0) _seed += RAND_m;
+        return _seed;
+    }
+    static inline double _scurve(double t) {
+        return t * t * (3.0 - 2.0*t);
+    }
+    static inline double _lerp(double t, double a, double b) {
+        return a + t * (b-a);
+    }
+
+    // random number generator constants
+    static long const
+        RAND_m = 2147483647, // 2**31 - 1
+        RAND_a = 16807, // 7**5; primitive root of m
+        RAND_q = 127773, // m / a
+        RAND_r = 2836; // m % a
+
+    // other constants
+    static int const
+        BSize = 0x100,
+        BMask = 0xff;
+    static double const
+        PerlinOffset = 4096.0;
+
+    Geom::Rect _tile;
+    Geom::Point _baseFreq;
+    int _latticeSelector[2*BSize + 2];
+    double _gradient[2*BSize + 2][4][2];
+    long _seed;
+    int _octaves;
+    bool _stitchTiles;
+    int _wrapx, _wrapy, _wrapw, _wraph;
+    bool _inited;
+    bool _fractalnoise;
 };
 
 FilterTurbulence::FilterTurbulence()
-: XbaseFrequency(0),
-  YbaseFrequency(0),
-  numOctaves(1),
-  seed(0),
-  updated(false),
-  updated_area(NR::IPoint(), NR::IPoint()),
-  pix(NULL),
-  fTileWidth(10), //guessed
-  fTileHeight(10), //guessed
-  fTileX(1), //guessed
-  fTileY(1) //guessed
+    : gen(new TurbulenceGenerator())
+    , XbaseFrequency(0)
+    , YbaseFrequency(0)
+    , numOctaves(1)
+    , seed(0)
+    , updated(false)
+    , updated_area(NR::IPoint(), NR::IPoint())
+    , pix(NULL)
+    , fTileWidth(10) //guessed
+    , fTileHeight(10) //guessed
+    , fTileX(1) //guessed
+    , fTileY(1) //guessed
 {
 }
 
@@ -78,6 +315,8 @@ FilterPrimitive * FilterTurbulence::create() {
 
 FilterTurbulence::~FilterTurbulence()
 {
+    delete gen;
+
     if (pix) {
         nr_pixblock_release(pix);
         delete pix;
@@ -87,26 +326,30 @@ FilterTurbulence::~FilterTurbulence()
 void FilterTurbulence::set_baseFrequency(int axis, double freq){
     if (axis==0) XbaseFrequency=freq;
     if (axis==1) YbaseFrequency=freq;
+    gen->dirty();
 }
 
 void FilterTurbulence::set_numOctaves(int num){
-    numOctaves=num;
+    numOctaves = num;
+    gen->dirty();
 }
 
 void FilterTurbulence::set_seed(double s){
-    seed=s;
+    seed = s;
+    gen->dirty();
 }
 
 void FilterTurbulence::set_stitchTiles(bool st){
-    stitchTiles=st;
+    stitchTiles = st;
+    gen->dirty();
 }
 
 void FilterTurbulence::set_type(FilterTurbulenceType t){
-    type=t;
+    type = t;
+    gen->dirty();
 }
 
 void FilterTurbulence::set_updated(bool u){
-    updated=u;
 }
 
 void FilterTurbulence::render_area(NRPixBlock *pix, NR::IRect &full_area, FilterUnits const &units) {
@@ -160,7 +403,7 @@ void FilterTurbulence::update_pixbuffer(NR::IRect &area, FilterUnits const &unit
     int bbox_x1 = area.max()[NR::X];
     int bbox_y1 = area.max()[NR::Y];
 
-    TurbulenceInit((long)seed);
+    //TurbulenceInit((long)seed);
 
     if (!pix){
         pix = new NRPixBlock;
@@ -192,67 +435,41 @@ void FilterTurbulence::update_pixbuffer(NR::IRect &area, FilterUnits const &unit
     updated_area = area;
 }
 
+struct Turbulence {
+    Turbulence(TurbulenceGenerator const &gen, Geom::Matrix const &trans, int x0, int y0)
+        : _gen(gen)
+        , _trans(trans)
+        , _x0(x0), _y0(y0)
+    {}
+    guint32 operator()(int x, int y) {
+        Geom::Point point(x + _x0, y + _y0);
+        point *= _trans;
+        return _gen.turbulencePixel(point);
+    }
+private:
+    TurbulenceGenerator const &_gen;
+    Geom::Matrix _trans;
+    int _x0, _y0;
+};
+
 void FilterTurbulence::render_cairo(FilterSlot &slot)
 {
     cairo_surface_t *input = slot.getcairo(_input);
     cairo_surface_t *out = ink_cairo_surface_create_same_size(input, CAIRO_CONTENT_COLOR_ALPHA);
 
-    if (!updated) {
-        TurbulenceInit((long)seed);
-        updated = true;
+    if (!gen->ready()) {
+        Geom::Point ta(fTileX, fTileY);
+        Geom::Point tb(fTileX + fTileWidth, fTileY + fTileHeight);
+        gen->init(seed, Geom::Rect(ta, tb),
+            Geom::Point(XbaseFrequency, YbaseFrequency), stitchTiles,
+            type == TURBULENCE_FRACTALNOISE, numOctaves);
     }
 
     // TODO: convert this to ink_cairo_surface_synthesize
     Geom::Matrix unit_trans = slot.get_units().get_matrix_primitiveunits2pb().inverse();
     NRRectL const &slot_area = slot.get_slot_area();
 
-    int w = cairo_image_surface_get_width(out);
-    int h = cairo_image_surface_get_height(out);
-    int stride = cairo_image_surface_get_stride(out);
-    unsigned char *data = cairo_image_surface_get_data(out);
-
-    if (type == TURBULENCE_TURBULENCE) {
-        for (int i = 0; i < h; ++i) {
-            guint32 *out_p = reinterpret_cast<guint32*>(data + i*stride);
-            for (int j = 0; j < w; ++j) {
-                Geom::Point pt(slot_area.x0 + j, slot_area.y0 + i);
-                pt *= unit_trans;
-
-                guint32 r = CLAMP_D_TO_U8(turbulence(0, pt)*255);
-                guint32 g = CLAMP_D_TO_U8(turbulence(1, pt)*255);
-                guint32 b = CLAMP_D_TO_U8(turbulence(2, pt)*255);
-                guint32 a = CLAMP_D_TO_U8(turbulence(3, pt)*255);
-
-                r = premul_alpha(r, a);
-                g = premul_alpha(g, a);
-                b = premul_alpha(b, a);
-
-                ASSEMBLE_ARGB32(result, a,r,g,b)
-                *out_p++ = result;
-            }
-        }
-    } else {
-        // TURBULENCE_FRACTALNOISE
-        for (int i = 0; i < h; ++i) {
-            guint32 *out_p = reinterpret_cast<guint32*>(data + i*stride);
-            for (int j = 0; j < w; ++j) {
-                Geom::Point pt(slot_area.x0 + j, slot_area.y0 + i);
-                pt *= unit_trans;
-
-                guint32 r = CLAMP_D_TO_U8((turbulence(0, pt)*255 + 255)/2);
-                guint32 g = CLAMP_D_TO_U8((turbulence(1, pt)*255 + 255)/2);
-                guint32 b = CLAMP_D_TO_U8((turbulence(2, pt)*255 + 255)/2);
-                guint32 a = CLAMP_D_TO_U8((turbulence(3, pt)*255 + 255)/2);
-
-                r = premul_alpha(r, a);
-                g = premul_alpha(g, a);
-                b = premul_alpha(b, a);
-
-                ASSEMBLE_ARGB32(result, a,r,g,b)
-                *out_p++ = result;
-            }
-        }
-    }
+    ink_cairo_surface_synthesize(out, Turbulence(*gen, unit_trans, slot_area.x0, slot_area.y0));
 
     cairo_surface_mark_dirty(out);
 
@@ -260,193 +477,6 @@ void FilterTurbulence::render_cairo(FilterSlot &slot)
     cairo_surface_destroy(out);
 }
 
-#if 0
-int FilterTurbulence::render(FilterSlot &slot, FilterUnits const &units) {
-    NR::IRect area = units.get_pixblock_filterarea_paraller();
-    // TODO: could be faster - updated_area only has to be same size as area
-    if (!updated || updated_area != area) update_pixbuffer(area, units);
-
-    NRPixBlock *in = slot.get(_input);
-    if (!in) {
-        g_warning("Missing source image for feTurbulence (in=%d)", _input);
-        return 1;
-    }
-
-    NRPixBlock *out = new NRPixBlock;
-    int x0 = in->area.x0, y0 = in->area.y0;
-    int x1 = in->area.x1, y1 = in->area.y1;
-    nr_pixblock_setup_fast(out, NR_PIXBLOCK_MODE_R8G8B8A8N, x0, y0, x1, y1, true);
-
-    if (pix_data) {
-        /* If pre-rendered output of whole filter area exists, just copy it. */
-        nr_blit_pixblock_pixblock(out, pix);
-    } else {
-        /* No pre-rendered output, render the required area here. */
-        render_area(out, area, units);
-    }
-
-    out->empty = FALSE;
-    slot.set(_output, out);
-    return 0;
-}
-#endif
-
-long FilterTurbulence::Turbulence_setup_seed(long lSeed)
-{
-  if (lSeed <= 0) lSeed = -(lSeed % (RAND_m - 1)) + 1;
-  if (lSeed > RAND_m - 1) lSeed = RAND_m - 1;
-  return lSeed;
-}
-
-long FilterTurbulence::TurbulenceRandom(long lSeed)
-{
-  long result;
-  result = RAND_a * (lSeed % RAND_q) - RAND_r * (lSeed / RAND_q);
-  if (result <= 0) result += RAND_m;
-  return result;
-}
-
-void FilterTurbulence::TurbulenceInit(long lSeed)
-{
-  double s;
-  int i, j, k;
-  lSeed = Turbulence_setup_seed(lSeed);
-  for(k = 0; k < 4; k++)
-  {
-    for(i = 0; i < BSize; i++)
-    {
-      uLatticeSelector[i] = i;
-      for (j = 0; j < 2; j++)
-        fGradient[k][i][j] = (double)(((lSeed = TurbulenceRandom(lSeed)) % (BSize + BSize)) - BSize) / BSize;
-      s = double(sqrt(fGradient[k][i][0] * fGradient[k][i][0] + fGradient[k][i][1] * fGradient[k][i][1]));
-      fGradient[k][i][0] /= s;
-      fGradient[k][i][1] /= s;
-    }
-  }
-  while(--i)
-  {
-    k = uLatticeSelector[i];
-    uLatticeSelector[i] = uLatticeSelector[j = (lSeed = TurbulenceRandom(lSeed)) % BSize];
-    uLatticeSelector[j] = k;
-  }
-  for(i = 0; i < BSize + 2; i++)
-  {
-    uLatticeSelector[BSize + i] = uLatticeSelector[i];
-    for(k = 0; k < 4; k++)
-      for(j = 0; j < 2; j++)
-        fGradient[k][BSize + i][j] = fGradient[k][i][j];
-  }
-}
-
-double FilterTurbulence::TurbulenceNoise2(int nColorChannel, double vec[2], StitchInfo *pStitchInfo)
-{
-  int bx0, bx1, by0, by1, b00, b10, b01, b11;
-  double rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v;
-  int i, j;
-  t = vec[0] + PerlinN;
-  bx0 = (int)t;
-  bx1 = bx0+1;
-  rx0 = t - (int)t;
-  rx1 = rx0 - 1.0f;
-  t = vec[1] + PerlinN;
-  by0 = (int)t;
-  by1 = by0+1;
-  ry0 = t - (int)t;
-  ry1 = ry0 - 1.0f;
-  // If stitching, adjust lattice points accordingly.
-  if(pStitchInfo != NULL)
-  {
-    if(bx0 >= pStitchInfo->nWrapX)
-      bx0 -= pStitchInfo->nWidth;
-    if(bx1 >= pStitchInfo->nWrapX)
-      bx1 -= pStitchInfo->nWidth;
-    if(by0 >= pStitchInfo->nWrapY)
-      by0 -= pStitchInfo->nHeight;
-    if(by1 >= pStitchInfo->nWrapY)
-      by1 -= pStitchInfo->nHeight;
-  }
-  bx0 &= BM;
-  bx1 &= BM;
-  by0 &= BM;
-  by1 &= BM;
-  i = uLatticeSelector[bx0];
-  j = uLatticeSelector[bx1];
-  b00 = uLatticeSelector[i + by0];
-  b10 = uLatticeSelector[j + by0];
-  b01 = uLatticeSelector[i + by1];
-  b11 = uLatticeSelector[j + by1];
-  sx = double(s_curve(rx0));
-  sy = double(s_curve(ry0));
-  q = fGradient[nColorChannel][b00]; u = rx0 * q[0] + ry0 * q[1];
-  q = fGradient[nColorChannel][b10]; v = rx1 * q[0] + ry0 * q[1];
-  a = turb_lerp(sx, u, v);
-  q = fGradient[nColorChannel][b01]; u = rx0 * q[0] + ry1 * q[1];
-  q = fGradient[nColorChannel][b11]; v = rx1 * q[0] + ry1 * q[1];
-  b = turb_lerp(sx, u, v);
-  return turb_lerp(sy, a, b);
-}
-
-double FilterTurbulence::turbulence(int nColorChannel, Geom::Point const &point)
-{
-  StitchInfo stitch;
-  StitchInfo *pStitchInfo = NULL; // Not stitching when NULL.
-  // Adjust the base frequencies if necessary for stitching.
-  if(stitchTiles)
-  {
-    // When stitching tiled turbulence, the frequencies must be adjusted
-    // so that the tile borders will be continuous.
-    if(XbaseFrequency != 0.0)
-    {
-      double fLoFreq = double(floor(fTileWidth * XbaseFrequency)) / fTileWidth;
-      double fHiFreq = double(ceil(fTileWidth * XbaseFrequency)) / fTileWidth;
-      if(XbaseFrequency / fLoFreq < fHiFreq / XbaseFrequency)
-        XbaseFrequency = fLoFreq;
-      else
-        XbaseFrequency = fHiFreq;
-    }
-    if(YbaseFrequency != 0.0)
-    {
-      double fLoFreq = double(floor(fTileHeight * YbaseFrequency)) / fTileHeight;
-      double fHiFreq = double(ceil(fTileHeight * YbaseFrequency)) / fTileHeight;
-      if(YbaseFrequency / fLoFreq < fHiFreq / YbaseFrequency)
-        YbaseFrequency = fLoFreq;
-      else
-        YbaseFrequency = fHiFreq;
-    }
-    // Set up TurbulenceInitial stitch values.
-    pStitchInfo = &stitch;
-    stitch.nWidth = int(fTileWidth * XbaseFrequency + 0.5f);
-    stitch.nWrapX = int(fTileX * XbaseFrequency + PerlinN + stitch.nWidth);
-    stitch.nHeight = int(fTileHeight * YbaseFrequency + 0.5f);
-    stitch.nWrapY = int(fTileY * YbaseFrequency + PerlinN + stitch.nHeight);
-  }
-  double fSum = 0.0f;
-  double vec[2];
-  vec[0] = point[Geom::X] * XbaseFrequency;
-  vec[1] = point[Geom::Y] * YbaseFrequency;
-  double ratio = 1;
-  for(int nOctave = 0; nOctave < numOctaves; nOctave++)
-  {
-    if(type==TURBULENCE_FRACTALNOISE)
-      fSum += double(TurbulenceNoise2(nColorChannel, vec, pStitchInfo) / ratio);
-    else
-      fSum += double(fabs(TurbulenceNoise2(nColorChannel, vec, pStitchInfo)) / ratio);
-    vec[0] *= 2;
-    vec[1] *= 2;
-    ratio *= 2;
-    if(pStitchInfo != NULL)
-    {
-      // Update stitch values. Subtracting PerlinN before the multiplication and
-      // adding it afterward simplifies to subtracting it once.
-      stitch.nWidth *= 2;
-      stitch.nWrapX = 2 * stitch.nWrapX - PerlinN;
-      stitch.nHeight *= 2;
-      stitch.nWrapY = 2 * stitch.nWrapY - PerlinN;
-    }
-  }
-  return fSum;
-}
-
 } /* namespace Filters */
 } /* namespace Inkscape */
 
diff --git a/src/display/nr-filter-turbulence.h b/src/display/nr-filter-turbulence.h
index 53ea5dd39..b2bc3a185 100644
--- a/src/display/nr-filter-turbulence.h
+++ b/src/display/nr-filter-turbulence.h
@@ -38,7 +38,8 @@ enum FilterTurbulenceType {
 
 struct StitchInfo;
 
-#define BSize 0x100
+//#define BSize 0x100
+class TurbulenceGenerator;
 
 class FilterTurbulence : public FilterPrimitive {
 public:
@@ -58,12 +59,16 @@ public:
     void set_updated(bool u);
 private:
 
+    TurbulenceGenerator *gen;
+
+    void turbulenceInit(long seed);
+/*
     long Turbulence_setup_seed(long lSeed);
     long TurbulenceRandom(long lSeed);
     void TurbulenceInit(long lSeed);
     double TurbulenceNoise2(int nColorChannel, double vec[2], StitchInfo *pStitchInfo);
     double turbulence(int nColorChannel, Geom::Point const &point);
-
+*/
     double XbaseFrequency, YbaseFrequency;
     int numOctaves;
     double seed;
@@ -74,14 +79,15 @@ private:
     NRPixBlock *pix;
     unsigned char *pix_data;
 
-    int uLatticeSelector[BSize + BSize + 2];
-    double fGradient[4][BSize + BSize + 2][2];
+    //int uLatticeSelector[BSize + BSize + 2];
+    //double fGradient[4][BSize + BSize + 2][2];
 
     double fTileWidth;
     double fTileHeight;
 
     double fTileX;
     double fTileY;
+
 };
 
 } /* namespace Filters */