Code

- new: Print Colors Preview Dialog and rendermode
[inkscape.git] / src / display / nr-filter-gaussian.cpp
index be0c9c5be34c7d5a99e1708c33aba02500930a26..d80782d2009bb7687a11ac811f02558feeee1e7a 100644 (file)
@@ -8,26 +8,38 @@
  *   bulia byak
  *   Jasper van de Gronde <th.v.d.gronde@hccnet.nl>
  *
- * Copyright (C) 2006 authors
+ * Copyright (C) 2006-2008 authors
  *
  * Released under GNU GPL, read the file 'COPYING' for more information
  */
 
+#include "config.h" // Needed for HAVE_OPENMP
+
 #include <algorithm>
 #include <cmath>
 #include <complex>
+#include <cstdlib>
 #include <glib.h>
 #include <limits>
+#if HAVE_OPENMP
+#include <omp.h>
+#endif //HAVE_OPENMP
 
-using std::isnormal;
+#include "2geom/isnan.h"
 
 #include "display/nr-filter-primitive.h"
 #include "display/nr-filter-gaussian.h"
 #include "display/nr-filter-types.h"
+#include "display/nr-filter-units.h"
+#include "libnr/nr-blit.h"
 #include "libnr/nr-pixblock.h"
-#include "libnr/nr-matrix.h"
+#include <2geom/matrix.h>
 #include "util/fixed_point.h"
-#include "prefs-utils.h"
+#include "preferences.h"
+
+#ifndef INK_UNUSED
+#define INK_UNUSED(x) ((void)(x))
+#endif
 
 // IIR filtering method based on:
 // L.J. van Vliet, I.T. Young, and P.W. Verbeek, Recursive Gaussian Derivative Filters,
@@ -38,21 +50,21 @@ using std::isnormal;
 // Using the backwards-pass initialization procedure from:
 // Boundary Conditions for Young - van Vliet Recursive Filtering
 // Bill Triggs, Michael Sdika
-// IEEE Transactions on Signal Processing, Volume 54, Number 5 - may 2006 
+// IEEE Transactions on Signal Processing, Volume 54, Number 5 - may 2006
 
 // Number of IIR filter coefficients used. Currently only 3 is supported.
 // "Recursive Gaussian Derivative Filters" says this is enough though (and
 // some testing indeed shows that the quality doesn't improve much if larger
 // filters are used).
-const size_t N = 3;
+static size_t const N = 3;
 
 template<typename InIt, typename OutIt, typename Size>
-void copy_n(InIt beg_in, Size N, OutIt beg_out) {
+inline void copy_n(InIt beg_in, Size N, OutIt beg_out) {
     std::copy(beg_in, beg_in+N, beg_out);
 }
 
 // Type used for IIR filter coefficients (can be 10.21 signed fixed point, see Anisotropic Gaussian Filtering Using Fixed Point Arithmetic, Christoph H. Lampert & Oliver Wirjadi, 2006)
-typedef double IIRValue; 
+typedef double IIRValue;
 
 // Type used for FIR filter coefficients (can be 16.16 unsigned fixed point, should have 8 or more bits in the fractional part, the integer part should be capable of storing approximately 20*255)
 typedef Inkscape::Util::FixedPoint<unsigned int,16> FIRValue;
@@ -67,23 +79,23 @@ template<typename T> static inline T clip(T const& v, T const& a, T const& b) {
 
 template<typename Tt, typename Ts>
 static inline Tt round_cast(Ts const& v) {
-       static Ts const rndoffset(.5);
-       return static_cast<Tt>(v+rndoffset);
+    static Ts const rndoffset(.5);
+    return static_cast<Tt>(v+rndoffset);
 }
 
 template<typename Tt, typename Ts>
-static inline Tt clip_round_cast(Ts const& v) {
-       static Ts const rndoffset(.5);
-    if ( v < std::numeric_limits<Tt>::min() ) return std::numeric_limits<Tt>::min();
-    if ( v > std::numeric_limits<Tt>::max() ) return std::numeric_limits<Tt>::max();
-       return static_cast<Tt>(v+rndoffset);
+static inline Tt clip_round_cast(Ts const& v, Tt const minval=std::numeric_limits<Tt>::min(), Tt const maxval=std::numeric_limits<Tt>::max()) {
+    if ( v < minval ) return minval;
+    if ( v > maxval ) return maxval;
+    return round_cast<Tt>(v);
 }
 
-namespace NR {
+namespace Inkscape {
+namespace Filters {
 
 FilterGaussian::FilterGaussian()
 {
-    _deviation_x = _deviation_y = prefs_get_double_attribute("options.filtertest", "value", 1.0);
+    _deviation_x = _deviation_y = 0.0;
 }
 
 FilterPrimitive *FilterGaussian::create()
@@ -96,12 +108,14 @@ FilterGaussian::~FilterGaussian()
     // Nothing to do here
 }
 
-int _effect_area_scr(double deviation)
+static int
+_effect_area_scr(double const deviation)
 {
     return (int)std::ceil(deviation * 3.0);
 }
 
-void _make_kernel(FIRValue *kernel, double deviation)
+static void
+_make_kernel(FIRValue *const kernel, double const deviation)
 {
     int const scr_len = _effect_area_scr(deviation);
     double const d_sq = sqr(deviation) * 2;
@@ -109,18 +123,23 @@ void _make_kernel(FIRValue *kernel, double deviation)
 
     // Compute kernel and sum of coefficients
     // Note that actually only half the kernel is computed, as it is symmetric
-    FIRValue sum = 0;
-    for ( int i = 0; i <= scr_len ; i++ ) {
+    double sum = 0;
+    for ( int i = scr_len; i >= 0 ; i-- ) {
         k[i] = std::exp(-sqr(i) / d_sq);
-        sum += static_cast<FIRValue>(k[i]);
+        if ( i > 0 ) sum += k[i];
     }
-    // the sum of the complete kernel is twice as large (minus the center element to avoid counting it twice)
-    sum = 2*sum - static_cast<FIRValue>(k[0]);
-
-    // Normalize kernel
-    for ( int i = 0; i <= scr_len ; i++ ) {
-        kernel[i] = k[i]/static_cast<double>(sum);
+    // the sum of the complete kernel is twice as large (plus the center element which we skipped above to prevent counting it twice)
+    sum = 2*sum + k[0];
+
+    // Normalize kernel (making sure the sum is exactly 1)
+    double ksum = 0;
+    FIRValue kernelsum = 0;
+    for ( int i = scr_len; i >= 1 ; i-- ) {
+        ksum += k[i]/sum;
+        kernel[i] = ksum-static_cast<double>(kernelsum);
+        kernelsum += kernel[i];
     }
+    kernel[0] = FIRValue(1)-2*kernelsum;
 }
 
 // Return value (v) should satisfy:
@@ -130,7 +149,8 @@ void _make_kernel(FIRValue *kernel, double deviation)
 //  8<=32-2*v
 //  2*v<=24
 //  v<=12
-int _effect_subsample_step_log2(double deviation, int quality)
+static int
+_effect_subsample_step_log2(double const deviation, int const quality)
 {
     // To make sure FIR will always be used (unless the kernel is VERY big):
     //  deviation/step <= 3
@@ -177,7 +197,8 @@ int _effect_subsample_step_log2(double deviation, int quality)
  * Catches reading and writing outside the pixblock area.
  * When enabled, decreases filter rendering speed massively.
  */
-inline void _check_index(NRPixBlock const * const pb, int const location, int const line)
+static inline void
+_check_index(NRPixBlock const * const pb, int const location, int const line)
 {
     if (false) {
         int max_loc = pb->rs * (pb->area.y1 - pb->area.y0);
@@ -199,12 +220,6 @@ static void calcFilter(double const sigma, double b[N]) {
         // Compute scaled filter coefficients
         std::complex<double> const d1 = pow(d1_org, 1.0/q);
         double const d3 = pow(d3_org, 1.0/q);
-        double const absd1sqr = std::norm(d1);
-        double const re2d1 = 2*d1.real();
-        double const bscale = 1.0/(absd1sqr*d3);
-        b[2] = -bscale;
-        b[1] =  bscale*(d3+re2d1);
-        b[0] = -bscale*(absd1sqr+d3*re2d1);
         // Compute actual sigma^2
         double const ssqr = 2*(2*(d1/sqr(d1-1.)).real()+d3/sqr(d3-1.));
         if ( ssqr < sigmasqr ) {
@@ -214,6 +229,16 @@ static void calcFilter(double const sigma, double b[N]) {
         }
         s = sqrt(ssqr);
     } while(qend-qbeg>(sigma/(1<<30)));
+    // Compute filter coefficients
+    double const q = (qbeg+qend)/2;
+    std::complex<double> const d1 = pow(d1_org, 1.0/q);
+    double const d3 = pow(d3_org, 1.0/q);
+    double const absd1sqr = std::norm(d1); // d1*d2 = d1*conj(d1) = |d1|^2 = std::norm(d1)
+    double const re2d1 = 2*d1.real(); // d1+d2 = d1+conj(d1) = 2*real(d1)
+    double const bscale = 1.0/(absd1sqr*d3);
+    b[2] = -bscale;
+    b[1] =  bscale*(d3+re2d1);
+    b[0] = -bscale*(absd1sqr+d3*re2d1);
 }
 
 static void calcTriggsSdikaM(double const b[N], double M[N*N]) {
@@ -252,9 +277,24 @@ static void calcTriggsSdikaInitialization(double const M[N*N], IIRValue const uo
 }
 
 // Filters over 1st dimension
-template<typename PT, unsigned int PC>
-void filter2D_IIR(PT *dest, int dstr1, int dstr2, PT const *src, int sstr1, int sstr2, int n1, int n2, IIRValue const b[N+1], double const M[N*N], IIRValue *const tmpdata) {
+template<typename PT, unsigned int PC, bool PREMULTIPLIED_ALPHA>
+static void
+filter2D_IIR(PT *const dest, int const dstr1, int const dstr2,
+             PT const *const src, int const sstr1, int const sstr2,
+             int const n1, int const n2, IIRValue const b[N+1], double const M[N*N],
+             IIRValue *const tmpdata[], int const num_threads)
+{
+#if HAVE_OPENMP
+#pragma omp parallel for num_threads(num_threads)
+#else
+    INK_UNUSED(num_threads);
+#endif // HAVE_OPENMP
     for ( int c2 = 0 ; c2 < n2 ; c2++ ) {
+#if HAVE_OPENMP
+        unsigned int tid = omp_get_thread_num();
+#else
+        unsigned int tid = 0;
+#endif // HAVE_OPENMP
         // corresponding line in the source and output buffer
         PT const * srcimg = src  + c2*sstr2;
         PT       * dstimg = dest + c2*dstr2 + n1*dstr1;
@@ -272,23 +312,33 @@ void filter2D_IIR(PT *dest, int dstr1, int dstr2, PT const *src, int sstr1, int
             for(unsigned int i=1; i<N+1; i++) {
                 for(unsigned int c=0; c<PC; c++) u[0][c] += u[i][c]*b[i];
             }
-            copy_n(u[0], PC, tmpdata+c1*PC);
+            copy_n(u[0], PC, tmpdata[tid]+c1*PC);
         }
         // Backward pass
         IIRValue v[N+1][PC];
-        calcTriggsSdikaInitialization(M, u, iplus, iplus, b[0], v);
+        calcTriggsSdikaInitialization<PC>(M, u, iplus, iplus, b[0], v);
         dstimg -= dstr1;
-        for(unsigned int c=0; c<PC; c++) dstimg[c] = clip_round_cast<PT>(v[0][c]);
+        if ( PREMULTIPLIED_ALPHA ) {
+            dstimg[PC-1] = clip_round_cast<PT>(v[0][PC-1]);
+            for(unsigned int c=0; c<PC-1; c++) dstimg[c] = clip_round_cast<PT>(v[0][c], std::numeric_limits<PT>::min(), dstimg[PC-1]);
+        } else {
+            for(unsigned int c=0; c<PC; c++) dstimg[c] = clip_round_cast<PT>(v[0][c]);
+        }
         int c1=n1-1;
         while(c1-->0) {
             for(unsigned int i=N; i>0; i--) copy_n(v[i-1], PC, v[i]);
-            copy_n(tmpdata+c1*PC, PC, v[0]);
+            copy_n(tmpdata[tid]+c1*PC, PC, v[0]);
             for(unsigned int c=0; c<PC; c++) v[0][c] *= b[0];
             for(unsigned int i=1; i<N+1; i++) {
                 for(unsigned int c=0; c<PC; c++) v[0][c] += v[i][c]*b[i];
             }
             dstimg -= dstr1;
-            for(unsigned int c=0; c<PC; c++) dstimg[c] = clip_round_cast<PT>(v[0][c]);
+            if ( PREMULTIPLIED_ALPHA ) {
+                dstimg[PC-1] = clip_round_cast<PT>(v[0][PC-1]);
+                for(unsigned int c=0; c<PC-1; c++) dstimg[c] = clip_round_cast<PT>(v[0][c], std::numeric_limits<PT>::min(), dstimg[PC-1]);
+            } else {
+                for(unsigned int c=0; c<PC; c++) dstimg[c] = clip_round_cast<PT>(v[0][c]);
+            }
         }
     }
 }
@@ -297,10 +347,19 @@ void filter2D_IIR(PT *dest, int dstr1, int dstr2, PT const *src, int sstr1, int
 // Assumes kernel is symmetric
 // scr_len should be size of kernel - 1
 template<typename PT, unsigned int PC>
-void filter2D_FIR(PT *dst, int dstr1, int dstr2, PT const *src, int sstr1, int sstr2, int n1, int n2, FIRValue const *const kernel, int scr_len) {
+static void
+filter2D_FIR(PT *const dst, int const dstr1, int const dstr2,
+             PT const *const src, int const sstr1, int const sstr2,
+             int const n1, int const n2, FIRValue const *const kernel, int const scr_len, int const num_threads)
+{
     // Past pixels seen (to enable in-place operation)
     PT history[scr_len+1][PC];
 
+#if HAVE_OPENMP
+#pragma omp parallel for num_threads(num_threads) private(history)
+#else
+    INK_UNUSED(num_threads);
+#endif // HAVE_OPENMP
     for ( int c2 = 0 ; c2 < n2 ; c2++ ) {
 
         // corresponding line in the source buffer
@@ -371,8 +430,8 @@ void filter2D_FIR(PT *dst, int dstr1, int dstr2, PT const *src, int sstr1, int s
                 // store the result in bufx
                 dst[dst_disp + byte] = round_cast<PT>(sum);
 
-                // optimization: if there was no variation within this point's neighborhood, 
-                // skip ahead while we keep seeing the same last_in byte: 
+                // optimization: if there was no variation within this point's neighborhood,
+                // skip ahead while we keep seeing the same last_in byte:
                 // blurring flat color would not change it anyway
                 if (different_count <= 1) {
                     int pos = c1 + 1;
@@ -392,7 +451,11 @@ void filter2D_FIR(PT *dst, int dstr1, int dstr2, PT const *src, int sstr1, int s
 }
 
 template<typename PT, unsigned int PC>
-void downsample(PT *dst, int dstr1, int dstr2, int dn1, int dn2, PT const *src, int sstr1, int sstr2, int sn1, int sn2, int step1_l2, int step2_l2) {
+static void
+downsample(PT *const dst, int const dstr1, int const dstr2, int const dn1, int const dn2,
+           PT const *const src, int const sstr1, int const sstr2, int const sn1, int const sn2,
+           int const step1_l2, int const step2_l2)
+{
     unsigned int const divisor_l2 = step1_l2+step2_l2; // step1*step2=2^(step1_l2+step2_l2)
     unsigned int const round_offset = (1<<divisor_l2)/2;
     int const step1 = 1<<step1_l2;
@@ -422,7 +485,11 @@ void downsample(PT *dst, int dstr1, int dstr2, int dn1, int dn2, PT const *src,
 }
 
 template<typename PT, unsigned int PC>
-void upsample(PT *dst, int dstr1, int dstr2, unsigned int dn1, unsigned int dn2, PT const *src, int sstr1, int sstr2, unsigned int sn1, unsigned int sn2, unsigned int step1_l2, unsigned int step2_l2) {
+static void
+upsample(PT *const dst, int const dstr1, int const dstr2, unsigned int const dn1, unsigned int const dn2,
+         PT const *const src, int const sstr1, int const sstr2, unsigned int const sn1, unsigned int const sn2,
+         unsigned int const step1_l2, unsigned int const step2_l2)
+{
     assert(((sn1-1)<<step1_l2)>=dn1 && ((sn2-1)<<step2_l2)>=dn2); // The last pixel of the source image should fall outside the destination image
     unsigned int const divisor_l2 = step1_l2+step2_l2; // step1*step2=2^(step1_l2+step2_l2)
     unsigned int const round_offset = (1<<divisor_l2)/2;
@@ -470,38 +537,66 @@ void upsample(PT *dst, int dstr1, int dstr2, unsigned int dn1, unsigned int dn2,
     }
 }
 
-int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
+int FilterGaussian::render(FilterSlot &slot, FilterUnits const &units)
 {
+    // TODO: Meaningful return values? (If they're checked at all.)
+
     /* in holds the input pixblock */
-    NRPixBlock *in = slot.get(_input);
+    NRPixBlock *original_in = slot.get(_input);
 
     /* If to either direction, the standard deviation is zero or
      * input image is not defined,
      * a transparent black image should be returned. */
-    if (_deviation_x <= 0 || _deviation_y <= 0 || in == NULL) {
-        NRPixBlock *out = new NRPixBlock;
-        if (in == NULL) {
+    if (_deviation_x <= 0 || _deviation_y <= 0 || original_in == NULL) {
+        NRPixBlock *src = original_in;
+        if (src == NULL) {
+            g_warning("Missing source image for feGaussianBlur (in=%d)", _input);
             // A bit guessing here, but source graphic is likely to be of
             // right size
-            in = slot.get(NR_FILTER_SOURCEGRAPHIC);
+            src = slot.get(NR_FILTER_SOURCEGRAPHIC);
         }
-        nr_pixblock_setup_fast(out, in->mode, in->area.x0, in->area.y0,
-                               in->area.x1, in->area.y1, true);
-        if (out->data.px != NULL) {
+        NRPixBlock *out = new NRPixBlock;
+        nr_pixblock_setup_fast(out, src->mode, src->area.x0, src->area.y0,
+                               src->area.x1, src->area.y1, true);
+        if (out->size != NR_PIXBLOCK_SIZE_TINY && out->data.px != NULL) {
             out->empty = false;
             slot.set(_output, out);
         }
         return 0;
     }
 
+    // Gaussian blur is defined to operate on non-premultiplied color values.
+    //   So, convert the input first it uses non-premultiplied color values.
+    //   And please note that this should not be done AFTER resampling, as resampling a non-premultiplied image
+    //   does not simply yield a non-premultiplied version of the resampled premultiplied image!!!
+    NRPixBlock *in = original_in;
+    if (in->mode == NR_PIXBLOCK_MODE_R8G8B8A8N) {
+        in = nr_pixblock_new_fast(NR_PIXBLOCK_MODE_R8G8B8A8P,
+                                  original_in->area.x0, original_in->area.y0,
+                                  original_in->area.x1, original_in->area.y1,
+                                  false);
+        if (!in) {
+            // ran out of memory
+            return 0;
+        }
+        nr_blit_pixblock_pixblock(in, original_in);
+    }
+
+    Geom::Matrix trans = units.get_matrix_primitiveunits2pb();
+
     // Some common constants
     int const width_org = in->area.x1-in->area.x0, height_org = in->area.y1-in->area.y0;
     double const deviation_x_org = _deviation_x * trans.expansionX();
     double const deviation_y_org = _deviation_y * trans.expansionY();
     int const PC = NR_PIXBLOCK_BPP(in);
+#if HAVE_OPENMP
+    int const NTHREADS = std::max(1,std::min(8, Inkscape::Preferences::get()->getInt("/options/threading/numthreads", omp_get_num_procs())));
+#else
+    int const NTHREADS = 1;
+#endif // HAVE_OPENMP
 
     // Subsampling constants
-    int const quality = prefs_get_int_attribute("options.blurquality", "value", 0);\r
+    int const quality = slot.get_blurquality();
     int const x_step_l2 = _effect_subsample_step_log2(deviation_x_org, quality);
     int const y_step_l2 = _effect_subsample_step_log2(deviation_y_org, quality);
     int const x_step = 1<<x_step_l2;
@@ -528,17 +623,26 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
                                           in->area.x0/x_step+width, in->area.y0/y_step+height, true);
     if (out->size != NR_PIXBLOCK_SIZE_TINY && out->data.px == NULL) {
         // alas, we've accomplished a lot, but ran out of memory - so abort
+        if (in != original_in) nr_pixblock_free(in);
+        delete out;
         return 0;
     }
     // Temporary storage for IIR filter
     // NOTE: This can be eliminated, but it reduces the precision a bit
-    IIRValue * tmpdata = 0;
+    IIRValue * tmpdata[NTHREADS];
+    std::fill_n(tmpdata, NTHREADS, (IIRValue*)0);
     if ( use_IIR_x || use_IIR_y ) {
-        tmpdata = new IIRValue[std::max(width,height)*PC];
-        if (tmpdata == NULL) {
-            nr_pixblock_release(out);
-            delete out;
-            return 0;
+        for(int i=0; i<NTHREADS; i++) {
+            tmpdata[i] = new IIRValue[std::max(width,height)*PC];
+            if (tmpdata[i] == NULL) {
+                if (in != original_in) nr_pixblock_free(in);
+                nr_pixblock_release(out);
+                while(i-->0) {
+                    delete[] tmpdata[i];
+                }
+                delete out;
+                return 0;
+            }
         }
     }
     NRPixBlock *ssin = in;
@@ -552,10 +656,10 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
             downsample<unsigned char,3>(NR_PIXBLOCK_PX(out), 3, out->rs, width, height, NR_PIXBLOCK_PX(in), 3, in->rs, width_org, height_org, x_step_l2, y_step_l2);
             break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
-            downsample<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, width, height, NR_PIXBLOCK_PX(in), 4, in->rs, width_org, height_org, x_step_l2, y_step_l2);
-            break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
+        //case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
+        //    downsample<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, width, height, NR_PIXBLOCK_PX(in), 4, in->rs, width_org, height_org, x_step_l2, y_step_l2);
+        //    break;
+        case NR_PIXBLOCK_MODE_R8G8B8A8P: ///< Premultiplied 8 bit RGBA
             downsample<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, width, height, NR_PIXBLOCK_PX(in), 4, in->rs, width_org, height_org, x_step_l2, y_step_l2);
             break;
         default:
@@ -584,21 +688,21 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
         // Filter (x)
         switch(in->mode) {
         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
-            filter2D_IIR<unsigned char,1>(NR_PIXBLOCK_PX(out), 1, out->rs, NR_PIXBLOCK_PX(ssin), 1, ssin->rs, width, height, b, M, tmpdata);
+            filter2D_IIR<unsigned char,1,false>(NR_PIXBLOCK_PX(out), 1, out->rs, NR_PIXBLOCK_PX(ssin), 1, ssin->rs, width, height, b, M, tmpdata, NTHREADS);
             break;
         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
-            filter2D_IIR<unsigned char,3>(NR_PIXBLOCK_PX(out), 3, out->rs, NR_PIXBLOCK_PX(ssin), 3, ssin->rs, width, height, b, M, tmpdata);
-            break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
-            filter2D_IIR<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, NR_PIXBLOCK_PX(ssin), 4, ssin->rs, width, height, b, M, tmpdata);
+            filter2D_IIR<unsigned char,3,false>(NR_PIXBLOCK_PX(out), 3, out->rs, NR_PIXBLOCK_PX(ssin), 3, ssin->rs, width, height, b, M, tmpdata, NTHREADS);
             break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
-            filter2D_IIR<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, NR_PIXBLOCK_PX(ssin), 4, ssin->rs, width, height, b, M, tmpdata);
+        //case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
+        //    filter2D_IIR<unsigned char,4,false>(NR_PIXBLOCK_PX(out), 4, out->rs, NR_PIXBLOCK_PX(ssin), 4, ssin->rs, width, height, b, M, tmpdata, NTHREADS);
+        //    break;
+        case NR_PIXBLOCK_MODE_R8G8B8A8P: ///< Premultiplied 8 bit RGBA
+            filter2D_IIR<unsigned char,4,true >(NR_PIXBLOCK_PX(out), 4, out->rs, NR_PIXBLOCK_PX(ssin), 4, ssin->rs, width, height, b, M, tmpdata, NTHREADS);
             break;
         default:
             assert(false);
         };
-    } else { // !use_IIR_x
+    } else if ( scr_len_x > 1 ) { // !use_IIR_x
         // Filter kernel for x direction
         FIRValue kernel[scr_len_x];
         _make_kernel(kernel, deviation_x);
@@ -606,20 +710,22 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
         // Filter (x)
         switch(in->mode) {
         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
-            filter2D_FIR<unsigned char,1>(NR_PIXBLOCK_PX(out), 1, out->rs, NR_PIXBLOCK_PX(ssin), 1, ssin->rs, width, height, kernel, scr_len_x);
+            filter2D_FIR<unsigned char,1>(NR_PIXBLOCK_PX(out), 1, out->rs, NR_PIXBLOCK_PX(ssin), 1, ssin->rs, width, height, kernel, scr_len_x, NTHREADS);
             break;
         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
-            filter2D_FIR<unsigned char,3>(NR_PIXBLOCK_PX(out), 3, out->rs, NR_PIXBLOCK_PX(ssin), 3, ssin->rs, width, height, kernel, scr_len_x);
+            filter2D_FIR<unsigned char,3>(NR_PIXBLOCK_PX(out), 3, out->rs, NR_PIXBLOCK_PX(ssin), 3, ssin->rs, width, height, kernel, scr_len_x, NTHREADS);
             break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
-            filter2D_FIR<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, NR_PIXBLOCK_PX(ssin), 4, ssin->rs, width, height, kernel, scr_len_x);
-            break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
-            filter2D_FIR<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, NR_PIXBLOCK_PX(ssin), 4, ssin->rs, width, height, kernel, scr_len_x);
+        //case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
+        //    filter2D_FIR<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, NR_PIXBLOCK_PX(ssin), 4, ssin->rs, width, height, kernel, scr_len_x, NTHREADS);
+        //    break;
+        case NR_PIXBLOCK_MODE_R8G8B8A8P: ///< Premultiplied 8 bit RGBA
+            filter2D_FIR<unsigned char,4>(NR_PIXBLOCK_PX(out), 4, out->rs, NR_PIXBLOCK_PX(ssin), 4, ssin->rs, width, height, kernel, scr_len_x, NTHREADS);
             break;
         default:
             assert(false);
         };
+    } else if ( out != ssin ) { // out can be equal to ssin if resampling is used
+        nr_blit_pixblock_pixblock(out, ssin);
     }
 
     if (use_IIR_y) {
@@ -643,21 +749,21 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
         // Filter (y)
         switch(in->mode) {
         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
-            filter2D_IIR<unsigned char,1>(NR_PIXBLOCK_PX(out), out->rs, 1, NR_PIXBLOCK_PX(out), out->rs, 1, height, width, b, M, tmpdata);
+            filter2D_IIR<unsigned char,1,false>(NR_PIXBLOCK_PX(out), out->rs, 1, NR_PIXBLOCK_PX(out), out->rs, 1, height, width, b, M, tmpdata, NTHREADS);
             break;
         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
-            filter2D_IIR<unsigned char,3>(NR_PIXBLOCK_PX(out), out->rs, 3, NR_PIXBLOCK_PX(out), out->rs, 3, height, width, b, M, tmpdata);
-            break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
-            filter2D_IIR<unsigned char,4>(NR_PIXBLOCK_PX(out), out->rs, 4, NR_PIXBLOCK_PX(out), out->rs, 4, height, width, b, M, tmpdata);
+            filter2D_IIR<unsigned char,3,false>(NR_PIXBLOCK_PX(out), out->rs, 3, NR_PIXBLOCK_PX(out), out->rs, 3, height, width, b, M, tmpdata, NTHREADS);
             break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
-            filter2D_IIR<unsigned char,4>(NR_PIXBLOCK_PX(out), out->rs, 4, NR_PIXBLOCK_PX(out), out->rs, 4, height, width, b, M, tmpdata);
+        //case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
+        //    filter2D_IIR<unsigned char,4,false>(NR_PIXBLOCK_PX(out), out->rs, 4, NR_PIXBLOCK_PX(out), out->rs, 4, height, width, b, M, tmpdata, NTHREADS);
+        //    break;
+        case NR_PIXBLOCK_MODE_R8G8B8A8P: ///< Premultiplied 8 bit RGBA
+            filter2D_IIR<unsigned char,4,true >(NR_PIXBLOCK_PX(out), out->rs, 4, NR_PIXBLOCK_PX(out), out->rs, 4, height, width, b, M, tmpdata, NTHREADS);
             break;
         default:
             assert(false);
         };
-    } else { // !use_IIR_y
+    } else if ( scr_len_y > 1 ) { // !use_IIR_y
         // Filter kernel for y direction
         FIRValue kernel[scr_len_y];
         _make_kernel(kernel, deviation_y);
@@ -665,23 +771,25 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
         // Filter (y)
         switch(in->mode) {
         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
-            filter2D_FIR<unsigned char,1>(NR_PIXBLOCK_PX(out), out->rs, 1, NR_PIXBLOCK_PX(out), out->rs, 1, height, width, kernel, scr_len_y);
+            filter2D_FIR<unsigned char,1>(NR_PIXBLOCK_PX(out), out->rs, 1, NR_PIXBLOCK_PX(out), out->rs, 1, height, width, kernel, scr_len_y, NTHREADS);
             break;
         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
-            filter2D_FIR<unsigned char,3>(NR_PIXBLOCK_PX(out), out->rs, 3, NR_PIXBLOCK_PX(out), out->rs, 3, height, width, kernel, scr_len_y);
-            break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
-            filter2D_FIR<unsigned char,4>(NR_PIXBLOCK_PX(out), out->rs, 4, NR_PIXBLOCK_PX(out), out->rs, 4, height, width, kernel, scr_len_y);
+            filter2D_FIR<unsigned char,3>(NR_PIXBLOCK_PX(out), out->rs, 3, NR_PIXBLOCK_PX(out), out->rs, 3, height, width, kernel, scr_len_y, NTHREADS);
             break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
-            filter2D_FIR<unsigned char,4>(NR_PIXBLOCK_PX(out), out->rs, 4, NR_PIXBLOCK_PX(out), out->rs, 4, height, width, kernel, scr_len_y);
+        //case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
+        //    filter2D_FIR<unsigned char,4>(NR_PIXBLOCK_PX(out), out->rs, 4, NR_PIXBLOCK_PX(out), out->rs, 4, height, width, kernel, scr_len_y, NTHREADS);
+        //    break;
+        case NR_PIXBLOCK_MODE_R8G8B8A8P: ///< Premultiplied 8 bit RGBA
+            filter2D_FIR<unsigned char,4>(NR_PIXBLOCK_PX(out), out->rs, 4, NR_PIXBLOCK_PX(out), out->rs, 4, height, width, kernel, scr_len_y, NTHREADS);
             break;
         default:
             assert(false);
         };
     }
 
-    delete[] tmpdata; // deleting a nullptr has no effect, so this is save
+    for(int i=0; i<NTHREADS; i++) {
+        delete[] tmpdata[i]; // deleting a nullptr has no effect, so this is safe
+    }
 
     if ( !resampling ) {
         // No upsampling needed
@@ -694,6 +802,7 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
                                                    in->area.x1, in->area.y1, true);
         if (finalout->size != NR_PIXBLOCK_SIZE_TINY && finalout->data.px == NULL) {
             // alas, we've accomplished a lot, but ran out of memory - so abort
+            if (in != original_in) nr_pixblock_free(in);
             nr_pixblock_release(out);
             delete out;
             return 0;
@@ -707,10 +816,10 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
             upsample<unsigned char,3>(NR_PIXBLOCK_PX(finalout), 3, finalout->rs, width_org, height_org, NR_PIXBLOCK_PX(out), 3, out->rs, width, height, x_step_l2, y_step_l2);
             break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
-            upsample<unsigned char,4>(NR_PIXBLOCK_PX(finalout), 4, finalout->rs, width_org, height_org, NR_PIXBLOCK_PX(out), 4, out->rs, width, height, x_step_l2, y_step_l2);
-            break;
-        case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
+        //case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
+        //    upsample<unsigned char,4>(NR_PIXBLOCK_PX(finalout), 4, finalout->rs, width_org, height_org, NR_PIXBLOCK_PX(out), 4, out->rs, width, height, x_step_l2, y_step_l2);
+        //    break;
+        case NR_PIXBLOCK_MODE_R8G8B8A8P: ///< Premultiplied 8 bit RGBA
             upsample<unsigned char,4>(NR_PIXBLOCK_PX(finalout), 4, finalout->rs, width_org, height_org, NR_PIXBLOCK_PX(out), 4, out->rs, width, height, x_step_l2, y_step_l2);
             break;
         default:
@@ -726,32 +835,45 @@ int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
         slot.set(_output, finalout);
     }
 
+    if (in != original_in) nr_pixblock_free(in);
+
     return 0;
 }
 
-int FilterGaussian::get_enlarge(Matrix const &trans)
+void FilterGaussian::area_enlarge(NRRectL &area, Geom::Matrix const &trans)
 {
     int area_x = _effect_area_scr(_deviation_x * trans.expansionX());
     int area_y = _effect_area_scr(_deviation_y * trans.expansionY());
-    return std::max(area_x, area_y);
+    // maximum is used because rotations can mix up these directions
+    // TODO: calculate a more tight-fitting rendering area
+    int area_max = std::max(area_x, area_y);
+    area.x0 -= area_max;
+    area.x1 += area_max;
+    area.y0 -= area_max;
+    area.y1 += area_max;
+}
+
+FilterTraits FilterGaussian::get_input_traits() {
+    return TRAIT_PARALLER;
 }
 
 void FilterGaussian::set_deviation(double deviation)
 {
-    if(isnormal(deviation) && deviation >= 0) {
+    if(IS_FINITE(deviation) && deviation >= 0) {
         _deviation_x = _deviation_y = deviation;
     }
 }
 
 void FilterGaussian::set_deviation(double x, double y)
 {
-    if(isnormal(x) && x >= 0 && isnormal(y) && y >= 0) {
+    if(IS_FINITE(x) && x >= 0 && IS_FINITE(y) && y >= 0) {
         _deviation_x = x;
         _deviation_y = y;
     }
 }
 
-} /* namespace NR */
+} /* namespace Filters */
+} /* namespace Inkscape */
 
 /*
   Local Variables: