path: root/src/display/nr-filter-gaussian.cpp

// SPDX-License-Identifier: GPL-2.0-or-later
/*
 * Gaussian blur renderer
 *
 * Authors:
 *   Niko Kiirala <niko@kiirala.com>
 *   bulia byak
 *   Jasper van de Gronde <th.v.d.gronde@hccnet.nl>
 *
 * Copyright (C) 2006-2008 authors
 *
 * Released under GNU GPL v2+, read the file 'COPYING' for more information.
 */

#ifdef HAVE_CONFIG_H
# include "config.h"  // only include where actually required!
#endif

#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdlib>
#include <glib.h>
#include <limits>
#if HAVE_OPENMP
#include <omp.h>
#endif //HAVE_OPENMP

#include "display/cairo-utils.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 "display/nr-filter-slot.h"
#include <2geom/affine.h>
#include "util/fixed_point.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,
// in: A.K. Jain, S. Venkatesh, B.C. Lovell (eds.),
// ICPR'98, Proc. 14th Int. Conference on Pattern Recognition (Brisbane, Aug. 16-20),
// IEEE Computer Society Press, Los Alamitos, 1998, 509-514.
//
// 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

// 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).
static size_t const N = 3;

template<typename InIt, typename OutIt, typename Size>
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;

// 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;

template<typename T> static inline T sqr(T const& v) { return v*v; }

template<typename T> static inline T clip(T const& v, T const& a, T const& b) {
    if ( v < a ) return a;
    if ( v > b ) return b;
    return v;
}

template<typename Tt, typename Ts>
static inline Tt round_cast(Ts v) {
    static Ts const rndoffset(.5);
    return static_cast<Tt>(v+rndoffset);
}
/*
template<>
inline unsigned char round_cast(double v) {
    // This (fast) rounding method is based on:
    // http://stereopsis.com/sree/fpu2006.html
#if G_BYTE_ORDER==G_LITTLE_ENDIAN
    double const dmr = 6755399441055744.0;
    v = v + dmr;
    return ((unsigned char*)&v)[0];
#elif G_BYTE_ORDER==G_BIG_ENDIAN
    double const dmr = 6755399441055744.0;
    v = v + dmr;
    return ((unsigned char*)&v)[7];
#else
    static double const rndoffset(.5);
    return static_cast<unsigned char>(v+rndoffset);
#endif
}*/

template<typename Tt, typename Ts>
static inline Tt clip_round_cast(Ts const v) {
    Ts const minval = std::numeric_limits<Tt>::min();
    Ts const maxval = std::numeric_limits<Tt>::max();
    Tt const minval_rounded = std::numeric_limits<Tt>::min();
    Tt const maxval_rounded = std::numeric_limits<Tt>::max();
    if ( v < minval ) return minval_rounded;
    if ( v > maxval ) return maxval_rounded;
    return round_cast<Tt>(v);
}

template<typename Tt, typename Ts>
static inline Tt clip_round_cast_varmax(Ts const v, Tt const maxval_rounded) {
    Ts const minval = std::numeric_limits<Tt>::min();
    Tt const maxval = maxval_rounded;
    Tt const minval_rounded = std::numeric_limits<Tt>::min();
    if ( v < minval ) return minval_rounded;
    if ( v > maxval ) return maxval_rounded;
    return round_cast<Tt>(v);
}

namespace Inkscape {
namespace Filters {

FilterGaussian::FilterGaussian()
{
    _deviation_x = _deviation_y = 0.0;
}

FilterPrimitive *FilterGaussian::create()
{
    return new FilterGaussian();
}

FilterGaussian::~FilterGaussian()
{
    // Nothing to do here
}

static int
_effect_area_scr(double const deviation)
{
    return (int)std::ceil(std::fabs(deviation) * 3.0);
}

static void
_make_kernel(FIRValue *const kernel, double const deviation)
{
    int const scr_len = _effect_area_scr(deviation);
    g_assert(scr_len >= 0);
    double const d_sq = sqr(deviation) * 2;
    double k[scr_len+1]; // This is only called for small kernel sizes (above approximately 10 coefficients the IIR filter is used)

    // Compute kernel and sum of coefficients
    // Note that actually only half the kernel is computed, as it is symmetric
    double sum = 0;
    for ( int i = scr_len; i >= 0 ; i-- ) {
        k[i] = std::exp(-sqr(i) / d_sq);
        if ( i > 0 ) sum += k[i];
    }
    // 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:
//  2^(2*v)*255<2^32
//  255<2^(32-2*v)
//  2^8<=2^(32-2*v)
//  8<=32-2*v
//  2*v<=24
//  v<=12
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
    //  deviation/3 <= step
    //  log(deviation/3) <= log(step)
    // So when x below is >= 1/3 FIR will almost always be used.
    // This means IIR is almost only used with the modes BETTER or BEST.
    int stepsize_l2;
    switch (quality) {
        case BLUR_QUALITY_WORST:
            // 2 == log(x*8/3))
            // 2^2 == x*2^3/3
            // x == 3/2
            stepsize_l2 = clip(static_cast<int>(log(deviation*(3./2.))/log(2.)), 0, 12);
            break;
        case BLUR_QUALITY_WORSE:
            // 2 == log(x*16/3))
            // 2^2 == x*2^4/3
            // x == 3/2^2
            stepsize_l2 = clip(static_cast<int>(log(deviation*(3./4.))/log(2.)), 0, 12);
            break;
        case BLUR_QUALITY_BETTER:
            // 2 == log(x*32/3))
            // 2 == x*2^5/3
            // x == 3/2^4
            stepsize_l2 = clip(static_cast<int>(log(deviation*(3./16.))/log(2.)), 0, 12);
            break;
        case BLUR_QUALITY_BEST:
            stepsize_l2 = 0; // no subsampling at all
            break;
        case BLUR_QUALITY_NORMAL:
        default:
            // 2 == log(x*16/3))
            // 2 == x*2^4/3
            // x == 3/2^3
            stepsize_l2 = clip(static_cast<int>(log(deviation*(3./8.))/log(2.)), 0, 12);
            break;
    }
    return stepsize_l2;
}

static void calcFilter(double const sigma, double b[N]) {
    assert(N==3);
    std::complex<double> const d1_org(1.40098,  1.00236);
    double const d3_org = 1.85132;
    double qbeg = 1; // Don't go lower than sigma==2 (we'd probably want a normal convolution in that case anyway)
    double qend = 2*sigma;
    double const sigmasqr = sqr(sigma);
    do { // Binary search for right q (a linear interpolation scheme is suggested, but this should work fine as well)
        double const q = (qbeg+qend)/2;
        // Compute scaled filter coefficients
        std::complex<double> const d1 = pow(d1_org, 1.0/q);
        double const d3 = pow(d3_org, 1.0/q);
        // Compute actual sigma^2
        double const ssqr = 2*(2*(d1/sqr(d1-1.)).real()+d3/sqr(d3-1.));
        if ( ssqr < sigmasqr ) {
            qbeg = q;
        } else {
            qend = q;
        }
    } 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]) {
    assert(N==3);
    double a1=b[0], a2=b[1], a3=b[2];
    double const Mscale = 1.0/((1+a1-a2+a3)*(1-a1-a2-a3)*(1+a2+(a1-a3)*a3));
    M[0] = 1-a2-a1*a3-sqr(a3);
    M[1] = (a1+a3)*(a2+a1*a3);
    M[2] = a3*(a1+a2*a3);
    M[3] = a1+a2*a3;
    M[4] = (1-a2)*(a2+a1*a3);
    M[5] = a3*(1-a2-a1*a3-sqr(a3));
    M[6] = a1*(a1+a3)+a2*(1-a2);
    M[7] = a1*(a2-sqr(a3))+a3*(1+a2*(a2-1)-sqr(a3));
    M[8] = a3*(a1+a2*a3);
    for(unsigned int i=0; i<9; i++) M[i] *= Mscale;
}

template<unsigned int SIZE>
static void calcTriggsSdikaInitialization(double const M[N*N], IIRValue const uold[N][SIZE], IIRValue const uplus[SIZE], IIRValue const vplus[SIZE], IIRValue const alpha, IIRValue vold[N][SIZE]) {
    for(unsigned int c=0; c<SIZE; c++) {
        double uminp[N];
        for(unsigned int i=0; i<N; i++) uminp[i] = uold[i][c] - uplus[c];
        for(unsigned int i=0; i<N; i++) {
            double voldf = 0;
            for(unsigned int j=0; j<N; j++) {
                voldf += uminp[j]*M[i*N+j];
            }
            // Properly takes care of the scaling coefficient alpha and vplus (which is already appropriately scaled)
            // This was arrived at by starting from a version of the blur filter that ignored the scaling coefficient
            // (and scaled the final output by alpha^2) and then gradually reintroducing the scaling coefficient.
            vold[i][c] = voldf*alpha;
            vold[i][c] += vplus[c];
        }
    }
}

// Filters over 1st dimension
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)
{
    assert(src && dest);

#if G_BYTE_ORDER == G_LITTLE_ENDIAN
    static unsigned int const alpha_PC = PC-1;
    #define PREMUL_ALPHA_LOOP for(unsigned int c=0; c<PC-1; ++c)
#else
    static unsigned int const alpha_PC = 0;
    #define PREMUL_ALPHA_LOOP for(unsigned int c=1; c<PC; ++c)
#endif

INK_UNUSED(num_threads); // to suppress unused argument compiler warning
#if HAVE_OPENMP
#pragma omp parallel for num_threads(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;
        // Border constants
        IIRValue imin[PC];  copy_n(srcimg + (0)*sstr1, PC, imin);
        IIRValue iplus[PC]; copy_n(srcimg + (n1-1)*sstr1, PC, iplus);
        // Forward pass
        IIRValue u[N+1][PC];
        for(unsigned int i=0; i<N; i++) copy_n(imin, PC, u[i]);
        for ( int c1 = 0 ; c1 < n1 ; c1++ ) {
            for(unsigned int i=N; i>0; i--) copy_n(u[i-1], PC, u[i]);
            copy_n(srcimg, PC, u[0]);
            srcimg += sstr1;
            for(unsigned int c=0; c<PC; c++) u[0][c] *= b[0];
            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[tid]+c1*PC);
        }
        // Backward pass
        IIRValue v[N+1][PC];
        calcTriggsSdikaInitialization<PC>(M, u, iplus, iplus, b[0], v);
        dstimg -= dstr1;
        if ( PREMULTIPLIED_ALPHA ) {
            dstimg[alpha_PC] = clip_round_cast<PT>(v[0][alpha_PC]);
            PREMUL_ALPHA_LOOP dstimg[c] = clip_round_cast_varmax<PT>(v[0][c], dstimg[alpha_PC]);
        } 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[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;
            if ( PREMULTIPLIED_ALPHA ) {
                dstimg[alpha_PC] = clip_round_cast<PT>(v[0][alpha_PC]);
                PREMUL_ALPHA_LOOP dstimg[c] = clip_round_cast_varmax<PT>(v[0][c], dstimg[alpha_PC]);
            } else {
                for(unsigned int c=0; c<PC; c++) dstimg[c] = clip_round_cast<PT>(v[0][c]);
            }
        }
    }
}

// Filters over 1st dimension
// Assumes kernel is symmetric
// Kernel should have scr_len+1 elements
template<typename PT, unsigned int PC>
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)
{
    assert(src && dst);

    // Past pixels seen (to enable in-place operation)
    PT history[scr_len+1][PC];

INK_UNUSED(num_threads); // suppresses unused argument compiler warning
#if HAVE_OPENMP
#pragma omp parallel for num_threads(num_threads) private(history)
#endif // HAVE_OPENMP
    for ( int c2 = 0 ; c2 < n2 ; c2++ ) {

        // corresponding line in the source buffer
        int const src_line = c2 * sstr2;

        // current line in the output buffer
        int const dst_line = c2 * dstr2;

        int skipbuf[4] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN};

        // history initialization
        PT imin[PC]; copy_n(src + src_line, PC, imin);
        for(int i=0; i<scr_len; i++) copy_n(imin, PC, history[i]);

        for ( int c1 = 0 ; c1 < n1 ; c1++ ) {

            int const src_disp = src_line + c1 * sstr1;
            int const dst_disp = dst_line + c1 * dstr1;

            // update history
            for(int i=scr_len; i>0; i--) copy_n(history[i-1], PC, history[i]);
            copy_n(src + src_disp, PC, history[0]);

            // for all bytes of the pixel
            for ( unsigned int byte = 0 ; byte < PC ; byte++) {

                if(skipbuf[byte] > c1) continue;

                FIRValue sum = 0;
                int last_in = -1;
                int different_count = 0;

                // go over our point's neighbours in the history
                for ( int i = 0 ; i <= scr_len ; i++ ) {
                    // value at the pixel
                    PT in_byte = history[i][byte];

                    // is it the same as last one we saw?
                    if(in_byte != last_in) different_count++;
                    last_in = in_byte;

                    // sum pixels weighted by the kernel
                    sum += in_byte * kernel[i];
                }

                // go over our point's neighborhood on x axis in the in buffer
                int nb_src_disp = src_disp + byte;
                for ( int i = 1 ; i <= scr_len ; i++ ) {
                    // the pixel we're looking at
                    int c1_in = c1 + i;
                    if (c1_in >= n1) {
                        c1_in = n1 - 1;
                    } else {
                        nb_src_disp += sstr1;
                    }

                    // value at the pixel
                    PT in_byte = src[nb_src_disp];

                    // is it the same as last one we saw?
                    if(in_byte != last_in) different_count++;
                    last_in = in_byte;

                    // sum pixels weighted by the kernel
                    sum += in_byte * kernel[i];
                }

                // 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:
                // blurring flat color would not change it anyway
                if (different_count <= 1) { // note that different_count is at least 1, because last_in is initialized to -1
                    int pos = c1 + 1;
                    int nb_src_disp = src_disp + (1+scr_len)*sstr1 + byte; // src_line + (pos+scr_len) * sstr1 + byte
                    int nb_dst_disp = dst_disp + (1)        *dstr1 + byte; // dst_line + (pos) * sstr1 + byte
                    while(pos + scr_len < n1 && src[nb_src_disp] == last_in) {
                        dst[nb_dst_disp] = last_in;
                        pos++;
                        nb_src_disp += sstr1;
                        nb_dst_disp += dstr1;
                    }
                    skipbuf[byte] = pos;
                }
            }
        }
    }
}

static void
gaussian_pass_IIR(Geom::Dim2 d, double deviation, cairo_surface_t *src, cairo_surface_t *dest,
    IIRValue **tmpdata, int num_threads)
{
    // Filter variables
    IIRValue b[N+1];  // scaling coefficient + filter coefficients (can be 10.21 fixed point)
    double bf[N];  // computed filter coefficients
    double M[N*N]; // matrix used for initialization procedure (has to be double)

    // Compute filter
    calcFilter(deviation, bf);
    for(double & i : bf) i = -i;
    b[0] = 1; // b[0] == alpha (scaling coefficient)
    for(size_t i=0; i<N; i++) {
        b[i+1] = bf[i];
        b[0] -= b[i+1];
    }

    // Compute initialization matrix
    calcTriggsSdikaM(bf, M);

    int stride = cairo_image_surface_get_stride(src);
    int w = cairo_image_surface_get_width(src);
    int h = cairo_image_surface_get_height(src);
    if (d != Geom::X) std::swap(w, h);

    // Filter
    switch (cairo_image_surface_get_format(src)) {
    case CAIRO_FORMAT_A8:        ///< Grayscale
        filter2D_IIR<unsigned char,1,false>(
            cairo_image_surface_get_data(dest), d == Geom::X ? 1 : stride, d == Geom::X ? stride : 1,
            cairo_image_surface_get_data(src),  d == Geom::X ? 1 : stride, d == Geom::X ? stride : 1,
            w, h, b, M, tmpdata, num_threads);
        break;
    case CAIRO_FORMAT_ARGB32: ///< Premultiplied 8 bit RGBA
        filter2D_IIR<unsigned char,4,true>(
            cairo_image_surface_get_data(dest), d == Geom::X ? 4 : stride, d == Geom::X ? stride : 4,
            cairo_image_surface_get_data(src),  d == Geom::X ? 4 : stride, d == Geom::X ? stride : 4,
            w, h, b, M, tmpdata, num_threads);
        break;
    default:
        g_warning("gaussian_pass_IIR: unsupported image format");
    };
}

static void
gaussian_pass_FIR(Geom::Dim2 d, double deviation, cairo_surface_t *src, cairo_surface_t *dest,
    int num_threads)
{
    int scr_len = _effect_area_scr(deviation);
    // Filter kernel for x direction
    std::vector<FIRValue> kernel(scr_len + 1);
    _make_kernel(&kernel[0], deviation);

    int stride = cairo_image_surface_get_stride(src);
    int w = cairo_image_surface_get_width(src);
    int h = cairo_image_surface_get_height(src);
    if (d != Geom::X) std::swap(w, h);

    // Filter (x)
    switch (cairo_image_surface_get_format(src)) {
    case CAIRO_FORMAT_A8:        ///< Grayscale
        filter2D_FIR<unsigned char,1>(
            cairo_image_surface_get_data(dest), d == Geom::X ? 1 : stride, d == Geom::X ? stride : 1,
            cairo_image_surface_get_data(src),  d == Geom::X ? 1 : stride, d == Geom::X ? stride : 1,
            w, h, &kernel[0], scr_len, num_threads);
        break;
    case CAIRO_FORMAT_ARGB32: ///< Premultiplied 8 bit RGBA
        filter2D_FIR<unsigned char,4>(
            cairo_image_surface_get_data(dest), d == Geom::X ? 4 : stride, d == Geom::X ? stride : 4,
            cairo_image_surface_get_data(src),  d == Geom::X ? 4 : stride, d == Geom::X ? stride : 4,
            w, h, &kernel[0], scr_len, num_threads);
        break;
    default:
        g_warning("gaussian_pass_FIR: unsupported image format");
    };
}

void FilterGaussian::render_cairo(FilterSlot &slot)
{
    cairo_surface_t *in = slot.getcairo(_input);
    if (!(in && ink_cairo_surface_get_width(in) && ink_cairo_surface_get_height(in))) {
        return;
    }

    // We may need to transform input surface to correct color interpolation space. The input surface
    // might be used as input to another primitive but it is likely that all the primitives in a given
    // filter use the same color interpolation space so we don't copy the input before converting.
    SPColorInterpolation ci_fp = SP_CSS_COLOR_INTERPOLATION_AUTO;
    if( _style ) {
        ci_fp = (SPColorInterpolation)_style->color_interpolation_filters.computed;
    }
    set_cairo_surface_ci( in, ci_fp );

    // zero deviation = no change in output
    if (_deviation_x <= 0 && _deviation_y <= 0) {
        cairo_surface_t *cp = ink_cairo_surface_copy(in);
        slot.set(_output, cp);
        cairo_surface_destroy(cp);
        return;
    }

    // Handle bounding box case.
    double dx = _deviation_x;
    double dy = _deviation_y;
    if( slot.get_units().get_primitive_units() == SP_FILTER_UNITS_OBJECTBOUNDINGBOX ) {
        Geom::OptRect const bbox = slot.get_units().get_item_bbox();
        if( bbox ) {
            dx *= (*bbox).width();
            dy *= (*bbox).height();
        }
    }

    Geom::Affine trans = slot.get_units().get_matrix_user2pb();

    double deviation_x_orig = dx * trans.expansionX();
    double deviation_y_orig = dy * trans.expansionY();

    int device_scale = slot.get_device_scale();

    deviation_x_orig *= device_scale;
    deviation_y_orig *= device_scale;

    cairo_format_t fmt = cairo_image_surface_get_format(in);
    int bytes_per_pixel = 0;
    switch (fmt) {
        case CAIRO_FORMAT_A8:
            bytes_per_pixel = 1; break;
        case CAIRO_FORMAT_ARGB32:
        default:
            bytes_per_pixel = 4; break;
    }

#if HAVE_OPENMP
    Inkscape::Preferences *prefs = Inkscape::Preferences::get();
    int threads = prefs->getIntLimited("/options/threading/numthreads", omp_get_num_procs(), 1, 256);
#else
    int threads = 1;
#endif

    int quality = slot.get_blurquality();
    int x_step = 1 << _effect_subsample_step_log2(deviation_x_orig, quality);
    int y_step = 1 << _effect_subsample_step_log2(deviation_y_orig, quality);
    bool resampling = x_step > 1 || y_step > 1;
    int w_orig = ink_cairo_surface_get_width(in);  // Pixels
    int h_orig = ink_cairo_surface_get_height(in);
    int w_downsampled = resampling ? static_cast<int>(ceil(static_cast<double>(w_orig)/x_step))+1 : w_orig;
    int h_downsampled = resampling ? static_cast<int>(ceil(static_cast<double>(h_orig)/y_step))+1 : h_orig;
    double deviation_x = deviation_x_orig / x_step;
    double deviation_y = deviation_y_orig / y_step;
    int scr_len_x = _effect_area_scr(deviation_x);
    int scr_len_y = _effect_area_scr(deviation_y);

    // Decide which filter to use for X and Y
    // This threshold was determined by trial-and-error for one specific machine,
    // so there's a good chance that it's not optimal.
    // Whatever you do, don't go below 1 (and preferably not even below 2), as
    // the IIR filter gets unstable there.
    bool use_IIR_x = deviation_x > 3;
    bool use_IIR_y = deviation_y > 3;

    // Temporary storage for IIR filter
    // NOTE: This can be eliminated, but it reduces the precision a bit
    IIRValue * tmpdata[threads];
    std::fill_n(tmpdata, threads, (IIRValue*)0);
    if ( use_IIR_x || use_IIR_y ) {
        for(int i = 0; i < threads; ++i) {
            tmpdata[i] = new IIRValue[std::max(w_downsampled,h_downsampled)*bytes_per_pixel];
        }
    }

    cairo_surface_t *downsampled = nullptr;
    if (resampling) {
        // Divide by device scale as w_downsampled is in pixels while
        // cairo_surface_create_similar() uses device units.
        downsampled = cairo_surface_create_similar(in, cairo_surface_get_content(in),
            w_downsampled/device_scale, h_downsampled/device_scale);
        cairo_t *ct = cairo_create(downsampled);
        cairo_scale(ct, static_cast<double>(w_downsampled)/w_orig, static_cast<double>(h_downsampled)/h_orig);
        cairo_set_source_surface(ct, in, 0, 0);
        cairo_paint(ct);
        cairo_destroy(ct);
    } else {
        downsampled = ink_cairo_surface_copy(in);
    }
    cairo_surface_flush(downsampled);

    if (scr_len_x > 0) {
        if (use_IIR_x) {
            gaussian_pass_IIR(Geom::X, deviation_x, downsampled, downsampled, tmpdata, threads);
        } else {
            gaussian_pass_FIR(Geom::X, deviation_x, downsampled, downsampled, threads);
        }
    }

    if (scr_len_y > 0) {
        if (use_IIR_y) {
            gaussian_pass_IIR(Geom::Y, deviation_y, downsampled, downsampled, tmpdata, threads);
        } else {
            gaussian_pass_FIR(Geom::Y, deviation_y, downsampled, downsampled, threads);
        }
    }

    // free the temporary data
    if ( use_IIR_x || use_IIR_y ) {
        for(int i = 0; i < threads; ++i) {
            delete[] tmpdata[i];
        }
    }

    cairo_surface_mark_dirty(downsampled);
    if (resampling) {
        cairo_surface_t *upsampled = cairo_surface_create_similar(downsampled, cairo_surface_get_content(downsampled),
            w_orig/device_scale, h_orig/device_scale);
        cairo_t *ct = cairo_create(upsampled);
        cairo_scale(ct, static_cast<double>(w_orig)/w_downsampled, static_cast<double>(h_orig)/h_downsampled);
        cairo_set_source_surface(ct, downsampled, 0, 0);
        cairo_paint(ct);
        cairo_destroy(ct);

        set_cairo_surface_ci( upsampled, ci_fp );

        slot.set(_output, upsampled);
        cairo_surface_destroy(upsampled);
        cairo_surface_destroy(downsampled);
    } else {
        set_cairo_surface_ci( downsampled, ci_fp );

        slot.set(_output, downsampled);
        cairo_surface_destroy(downsampled);
    }
}

void FilterGaussian::area_enlarge(Geom::IntRect &area, Geom::Affine const &trans)
{
    int area_x = _effect_area_scr(_deviation_x * trans.expansionX());
    int area_y = _effect_area_scr(_deviation_y * trans.expansionY());
    // 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.expandBy(area_max);
}

bool FilterGaussian::can_handle_affine(Geom::Affine const &)
{
    // Previously we tried to be smart and return true for rotations.
    // However, the transform passed here is NOT the total transform
    // from filter user space to screen.
    // TODO: fix this, or replace can_handle_affine() with isotropic().
    return false;
}

double FilterGaussian::complexity(Geom::Affine const &trans)
{
    int area_x = _effect_area_scr(_deviation_x * trans.expansionX());
    int area_y = _effect_area_scr(_deviation_y * trans.expansionY());
    return 2.0 * area_x * area_y;
}

void FilterGaussian::set_deviation(double deviation)
{
    if(std::isfinite(deviation) && deviation >= 0) {
        _deviation_x = _deviation_y = deviation;
    }
}

void FilterGaussian::set_deviation(double x, double y)
{
    if(std::isfinite(x) && x >= 0 && std::isfinite(y) && y >= 0) {
        _deviation_x = x;
        _deviation_y = y;
    }
}

} /* namespace Filters */
} /* namespace Inkscape */

/*
  Local Variables:
  mode:c++
  c-file-style:"stroustrup"
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
  indent-tabs-mode:nil
  fill-column:99
  End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :