src/display/nr-filter-gaussian.cpp

   1 #define __NR_FILTER_GAUSSIAN_CPP__
   2
   3 /*
   4  * Gaussian blur renderer
   5  *
   6  * Authors:
   7  *   Niko Kiirala <niko@kiirala.com>
   8  *   bulia byak
   9  *   Jasper van de Gronde <th.v.d.gronde@hccnet.nl>
  10  *
  11  * Copyright (C) 2006 authors
  12  *
  13  * Released under GNU GPL, read the file 'COPYING' for more information
  14  */
  15
  16 #include <algorithm>
  17 #include <cmath>
  18 #include <complex>
  19 #include <glib.h>
  20 #include <cstdlib>
  21 #include <limits>
  22
  23 #include "2geom/isnan.h"
  24
  25 #include "display/nr-filter-primitive.h"
  26 #include "display/nr-filter-gaussian.h"
  27 #include "display/nr-filter-types.h"
  28 #include "display/nr-filter-units.h"
  29 #include "libnr/nr-pixblock.h"
  30 #include "libnr/nr-matrix.h"
  31 #include "libnr/nr-matrix-fns.h"
  32 #include "util/fixed_point.h"
  33 #include "preferences.h"
  34
  35 // IIR filtering method based on:
  36 // L.J. van Vliet, I.T. Young, and P.W. Verbeek, Recursive Gaussian Derivative Filters,
  37 // in: A.K. Jain, S. Venkatesh, B.C. Lovell (eds.),
  38 // ICPR'98, Proc. 14th Int. Conference on Pattern Recognition (Brisbane, Aug. 16-20),
  39 // IEEE Computer Society Press, Los Alamitos, 1998, 509-514.
  40 //
  41 // Using the backwards-pass initialization procedure from:
  42 // Boundary Conditions for Young - van Vliet Recursive Filtering
  43 // Bill Triggs, Michael Sdika
  44 // IEEE Transactions on Signal Processing, Volume 54, Number 5 - may 2006
  45
  46 // Number of IIR filter coefficients used. Currently only 3 is supported.
  47 // "Recursive Gaussian Derivative Filters" says this is enough though (and
  48 // some testing indeed shows that the quality doesn't improve much if larger
  49 // filters are used).
  50 static size_t const N = 3;
  51
  52 template<typename InIt, typename OutIt, typename Size>
  53 void copy_n(InIt beg_in, Size N, OutIt beg_out) {
  54     std::copy(beg_in, beg_in+N, beg_out);
  55 }
  56
  57 // 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)
  58 typedef double IIRValue;
  59
  60 // 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)
  61 typedef Inkscape::Util::FixedPoint<unsigned int,16> FIRValue;
  62
  63 template<typename T> static inline T sqr(T const& v) { return v*v; }
  64
  65 template<typename T> static inline T clip(T const& v, T const& a, T const& b) {
  66     if ( v < a ) return a;
  67     if ( v > b ) return b;
  68     return v;
  69 }
  70
  71 template<typename Tt, typename Ts>
  72 static inline Tt round_cast(Ts const& v) {
  73     static Ts const rndoffset(.5);
  74     return static_cast<Tt>(v+rndoffset);
  75 }
  76
  77 template<typename Tt, typename Ts>
  78 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()) {
  79     if ( v < minval ) return minval;
  80     if ( v > maxval ) return maxval;
  81     return round_cast<Tt>(v);
  82 }
  83
  84 namespace NR {
  85
  86 FilterGaussian::FilterGaussian()
  87 {
  88     _deviation_x = _deviation_y = 0.0;
  89 }
  90
  91 FilterPrimitive *FilterGaussian::create()
  92 {
  93     return new FilterGaussian();
  94 }
  95
  96 FilterGaussian::~FilterGaussian()
  97 {
  98     // Nothing to do here
  99 }
 100
 101 static int
 102 _effect_area_scr(double const deviation)
 103 {
 104     return (int)std::ceil(deviation * 3.0);
 105 }
 106
 107 static void
 108 _make_kernel(FIRValue *const kernel, double const deviation)
 109 {
 110     int const scr_len = _effect_area_scr(deviation);
 111     double const d_sq = sqr(deviation) * 2;
 112     double k[scr_len+1]; // This is only called for small kernel sizes (above approximately 10 coefficients the IIR filter is used)
 113
 114     // Compute kernel and sum of coefficients
 115     // Note that actually only half the kernel is computed, as it is symmetric
 116     double sum = 0;
 117     for ( int i = scr_len; i >= 0 ; i-- ) {
 118         k[i] = std::exp(-sqr(i) / d_sq);
 119         if ( i > 0 ) sum += k[i];
 120     }
 121     // the sum of the complete kernel is twice as large (plus the center element which we skipped above to prevent counting it twice)
 122     sum = 2*sum + k[0];
 123
 124     // Normalize kernel (making sure the sum is exactly 1)
 125     double ksum = 0;
 126     FIRValue kernelsum = 0;
 127     for ( int i = scr_len; i >= 1 ; i-- ) {
 128         ksum += k[i]/sum;
 129         kernel[i] = ksum-static_cast<double>(kernelsum);
 130         kernelsum += kernel[i];
 131     }
 132     kernel[0] = FIRValue(1)-2*kernelsum;
 133 }
 134
 135 // Return value (v) should satisfy:
 136 //  2^(2*v)*255<2^32
 137 //  255<2^(32-2*v)
 138 //  2^8<=2^(32-2*v)
 139 //  8<=32-2*v
 140 //  2*v<=24
 141 //  v<=12
 142 static int
 143 _effect_subsample_step_log2(double const deviation, int const quality)
 144 {
 145     // To make sure FIR will always be used (unless the kernel is VERY big):
 146     //  deviation/step <= 3
 147     //  deviation/3 <= step
 148     //  log(deviation/3) <= log(step)
 149     // So when x below is >= 1/3 FIR will almost always be used.
 150     // This means IIR is almost only used with the modes BETTER or BEST.
 151     int stepsize_l2;
 152     switch (quality) {
 153         case BLUR_QUALITY_WORST:
 154             // 2 == log(x*8/3))
 155             // 2^2 == x*2^3/3
 156             // x == 3/2
 157             stepsize_l2 = clip(static_cast<int>(log(deviation*(3./2.))/log(2.)), 0, 12);
 158             break;
 159         case BLUR_QUALITY_WORSE:
 160             // 2 == log(x*16/3))
 161             // 2^2 == x*2^4/3
 162             // x == 3/2^2
 163             stepsize_l2 = clip(static_cast<int>(log(deviation*(3./4.))/log(2.)), 0, 12);
 164             break;
 165         case BLUR_QUALITY_BETTER:
 166             // 2 == log(x*32/3))
 167             // 2 == x*2^5/3
 168             // x == 3/2^4
 169             stepsize_l2 = clip(static_cast<int>(log(deviation*(3./16.))/log(2.)), 0, 12);
 170             break;
 171         case BLUR_QUALITY_BEST:
 172             stepsize_l2 = 0; // no subsampling at all
 173             break;
 174         case BLUR_QUALITY_NORMAL:
 175         default:
 176             // 2 == log(x*16/3))
 177             // 2 == x*2^4/3
 178             // x == 3/2^3
 179             stepsize_l2 = clip(static_cast<int>(log(deviation*(3./8.))/log(2.)), 0, 12);
 180             break;
 181     }
 182     return stepsize_l2;
 183 }
 184
 185 /**
 186  * Sanity check function for indexing pixblocks.
 187  * Catches reading and writing outside the pixblock area.
 188  * When enabled, decreases filter rendering speed massively.
 189  */
 190 static inline void
 191 _check_index(NRPixBlock const * const pb, int const location, int const line)
 192 {
 193     if (false) {
 194         int max_loc = pb->rs * (pb->area.y1 - pb->area.y0);
 195         if (location < 0 || location >= max_loc)
 196             g_warning("Location %d out of bounds (0 ... %d) at line %d", location, max_loc, line);
 197     }
 198 }
 199
 200 static void calcFilter(double const sigma, double b[N]) {
 201     assert(N==3);
 202     std::complex<double> const d1_org(1.40098,  1.00236);
 203     double const d3_org = 1.85132;
 204     double qbeg = 1; // Don't go lower than sigma==2 (we'd probably want a normal convolution in that case anyway)
 205     double qend = 2*sigma;
 206     double const sigmasqr = sqr(sigma);
 207     double s;
 208     do { // Binary search for right q (a linear interpolation scheme is suggested, but this should work fine as well)
 209         double const q = (qbeg+qend)/2;
 210         // Compute scaled filter coefficients
 211         std::complex<double> const d1 = pow(d1_org, 1.0/q);
 212         double const d3 = pow(d3_org, 1.0/q);
 213         double const absd1sqr = std::norm(d1);
 214         double const re2d1 = 2*d1.real();
 215         double const bscale = 1.0/(absd1sqr*d3);
 216         b[2] = -bscale;
 217         b[1] =  bscale*(d3+re2d1);
 218         b[0] = -bscale*(absd1sqr+d3*re2d1);
 219         // Compute actual sigma^2
 220         double const ssqr = 2*(2*(d1/sqr(d1-1.)).real()+d3/sqr(d3-1.));
 221         if ( ssqr < sigmasqr ) {
 222             qbeg = q;
 223         } else {
 224             qend = q;
 225         }
 226         s = sqrt(ssqr);
 227     } while(qend-qbeg>(sigma/(1<<30)));
 228 }
 229
 230 static void calcTriggsSdikaM(double const b[N], double M[N*N]) {
 231     assert(N==3);
 232     double a1=b[0], a2=b[1], a3=b[2];
 233     double const Mscale = 1.0/((1+a1-a2+a3)*(1-a1-a2-a3)*(1+a2+(a1-a3)*a3));
 234     M[0] = 1-a2-a1*a3-sqr(a3);
 235     M[1] = (a1+a3)*(a2+a1*a3);
 236     M[2] = a3*(a1+a2*a3);
 237     M[3] = a1+a2*a3;
 238     M[4] = (1-a2)*(a2+a1*a3);
 239     M[5] = a3*(1-a2-a1*a3-sqr(a3));
 240     M[6] = a1*(a1+a3)+a2*(1-a2);
 241     M[7] = a1*(a2-sqr(a3))+a3*(1+a2*(a2-1)-sqr(a3));
 242     M[8] = a3*(a1+a2*a3);
 243     for(unsigned int i=0; i<9; i++) M[i] *= Mscale;
 244 }
 245
 246 template<unsigned int SIZE>
 247 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]) {
 248     for(unsigned int c=0; c<SIZE; c++) {
 249         double uminp[N];
 250         for(unsigned int i=0; i<N; i++) uminp[i] = uold[i][c] - uplus[c];
 251         for(unsigned int i=0; i<N; i++) {
 252             double voldf = 0;
 253             for(unsigned int j=0; j<N; j++) {
 254                 voldf += uminp[j]*M[i*N+j];
 255             }
 256             // Properly takes care of the scaling coefficient alpha and vplus (which is already appropriately scaled)
 257             // This was arrived at by starting from a version of the blur filter that ignored the scaling coefficient
 258             // (and scaled the final output by alpha^2) and then gradually reintroducing the scaling coefficient.
 259             vold[i][c] = voldf*alpha;
 260             vold[i][c] += vplus[c];
 261         }
 262     }
 263 }
 264
 265 // Filters over 1st dimension
 266 template<typename PT, unsigned int PC, bool PREMULTIPLIED_ALPHA>
 267 static void
 268 filter2D_IIR(PT *const dest, int const dstr1, int const dstr2,
 269              PT const *const src, int const sstr1, int const sstr2,
 270              int const n1, int const n2, IIRValue const b[N+1], double const M[N*N],
 271              IIRValue *const tmpdata)
 272 {
 273     for ( int c2 = 0 ; c2 < n2 ; c2++ ) {
 274         // corresponding line in the source and output buffer
 275         PT const * srcimg = src  + c2*sstr2;
 276         PT       * dstimg = dest + c2*dstr2 + n1*dstr1;
 277         // Border constants
 278         IIRValue imin[PC];  copy_n(srcimg + (0)*sstr1, PC, imin);
 279         IIRValue iplus[PC]; copy_n(srcimg + (n1-1)*sstr1, PC, iplus);
 280         // Forward pass
 281         IIRValue u[N+1][PC];
 282         for(unsigned int i=0; i<N; i++) copy_n(imin, PC, u[i]);
 283         for ( int c1 = 0 ; c1 < n1 ; c1++ ) {
 284             for(unsigned int i=N; i>0; i--) copy_n(u[i-1], PC, u[i]);
 285             copy_n(srcimg, PC, u[0]);
 286             srcimg += sstr1;
 287             for(unsigned int c=0; c<PC; c++) u[0][c] *= b[0];
 288             for(unsigned int i=1; i<N+1; i++) {
 289                 for(unsigned int c=0; c<PC; c++) u[0][c] += u[i][c]*b[i];
 290             }
 291             copy_n(u[0], PC, tmpdata+c1*PC);
 292         }
 293         // Backward pass
 294         IIRValue v[N+1][PC];
 295         calcTriggsSdikaInitialization<PC>(M, u, iplus, iplus, b[0], v);
 296         dstimg -= dstr1;
 297         if ( PREMULTIPLIED_ALPHA ) {
 298             dstimg[PC-1] = clip_round_cast<PT>(v[0][PC-1]);
 299             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]);
 300         } else {
 301             for(unsigned int c=0; c<PC; c++) dstimg[c] = clip_round_cast<PT>(v[0][c]);
 302         }
 303         int c1=n1-1;
 304         while(c1-->0) {
 305             for(unsigned int i=N; i>0; i--) copy_n(v[i-1], PC, v[i]);
 306             copy_n(tmpdata+c1*PC, PC, v[0]);
 307             for(unsigned int c=0; c<PC; c++) v[0][c] *= b[0];
 308             for(unsigned int i=1; i<N+1; i++) {
 309                 for(unsigned int c=0; c<PC; c++) v[0][c] += v[i][c]*b[i];
 310             }
 311             dstimg -= dstr1;
 312             if ( PREMULTIPLIED_ALPHA ) {
 313                 dstimg[PC-1] = clip_round_cast<PT>(v[0][PC-1]);
 314                 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]);
 315             } else {
 316                 for(unsigned int c=0; c<PC; c++) dstimg[c] = clip_round_cast<PT>(v[0][c]);
 317             }
 318         }
 319     }
 320 }
 321
 322 // Filters over 1st dimension
 323 // Assumes kernel is symmetric
 324 // scr_len should be size of kernel - 1
 325 template<typename PT, unsigned int PC>
 326 static void
 327 filter2D_FIR(PT *const dst, int const dstr1, int const dstr2,
 328              PT const *const src, int const sstr1, int const sstr2,
 329              int const n1, int const n2, FIRValue const *const kernel, int const scr_len)
 330 {
 331     // Past pixels seen (to enable in-place operation)
 332     PT history[scr_len+1][PC];
 333
 334     for ( int c2 = 0 ; c2 < n2 ; c2++ ) {
 335
 336         // corresponding line in the source buffer
 337         int const src_line = c2 * sstr2;
 338
 339         // current line in the output buffer
 340         int const dst_line = c2 * dstr2;
 341
 342         int skipbuf[4] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN};
 343
 344         // history initialization
 345         PT imin[PC]; copy_n(src + src_line, PC, imin);
 346         for(int i=0; i<scr_len; i++) copy_n(imin, PC, history[i]);
 347
 348         for ( int c1 = 0 ; c1 < n1 ; c1++ ) {
 349
 350             int const src_disp = src_line + c1 * sstr1;
 351             int const dst_disp = dst_line + c1 * sstr1;
 352
 353             // update history
 354             for(int i=scr_len; i>0; i--) copy_n(history[i-1], PC, history[i]);
 355             copy_n(src + src_disp, PC, history[0]);
 356
 357             // for all bytes of the pixel
 358             for ( unsigned int byte = 0 ; byte < PC ; byte++) {
 359
 360                 if(skipbuf[byte] > c1) continue;
 361
 362                 FIRValue sum = 0;
 363                 int last_in = -1;
 364                 int different_count = 0;
 365
 366                 // go over our point's neighbours in the history
 367                 for ( int i = 0 ; i <= scr_len ; i++ ) {
 368                     // value at the pixel
 369                     PT in_byte = history[i][byte];
 370
 371                     // is it the same as last one we saw?
 372                     if(in_byte != last_in) different_count++;
 373                     last_in = in_byte;
 374
 375                     // sum pixels weighted by the kernel
 376                     sum += in_byte * kernel[i];
 377                 }
 378
 379                 // go over our point's neighborhood on x axis in the in buffer
 380                 int nb_src_disp = src_disp + byte;
 381                 for ( int i = 1 ; i <= scr_len ; i++ ) {
 382                     // the pixel we're looking at
 383                     int c1_in = c1 + i;
 384                     if (c1_in >= n1) {
 385                         c1_in = n1 - 1;
 386                     } else {
 387                         nb_src_disp += sstr1;
 388                     }
 389
 390                     // value at the pixel
 391                     PT in_byte = src[nb_src_disp];
 392
 393                     // is it the same as last one we saw?
 394                     if(in_byte != last_in) different_count++;
 395                     last_in = in_byte;
 396
 397                     // sum pixels weighted by the kernel
 398                     sum += in_byte * kernel[i];
 399                 }
 400
 401                 // store the result in bufx
 402                 dst[dst_disp + byte] = round_cast<PT>(sum);
 403
 404                 // optimization: if there was no variation within this point's neighborhood,
 405                 // skip ahead while we keep seeing the same last_in byte:
 406                 // blurring flat color would not change it anyway
 407                 if (different_count <= 1) {
 408                     int pos = c1 + 1;
 409                     int nb_src_disp = src_disp + (1+scr_len)*sstr1 + byte; // src_line + (pos+scr_len) * sstr1 + byte
 410                     int nb_dst_disp = dst_disp + (1)        *dstr1 + byte; // dst_line + (pos) * sstr1 + byte
 411                     while(pos + scr_len < n1 && src[nb_src_disp] == last_in) {
 412                         dst[nb_dst_disp] = last_in;
 413                         pos++;
 414                         nb_src_disp += sstr1;
 415                         nb_dst_disp += sstr1;
 416                     }
 417                     skipbuf[byte] = pos;
 418                 }
 419             }
 420         }
 421     }
 422 }
 423
 424 template<typename PT, unsigned int PC>
 425 static void
 426 downsample(PT *const dst, int const dstr1, int const dstr2, int const dn1, int const dn2,
 427            PT const *const src, int const sstr1, int const sstr2, int const sn1, int const sn2,
 428            int const step1_l2, int const step2_l2)
 429 {
 430     unsigned int const divisor_l2 = step1_l2+step2_l2; // step1*step2=2^(step1_l2+step2_l2)
 431     unsigned int const round_offset = (1<<divisor_l2)/2;
 432     int const step1 = 1<<step1_l2;
 433     int const step2 = 1<<step2_l2;
 434     int const step1_2 = step1/2;
 435     int const step2_2 = step2/2;
 436     for(int dc2 = 0 ; dc2 < dn2 ; dc2++) {
 437         int const sc2_begin = (dc2<<step2_l2)-step2_2;
 438         int const sc2_end = sc2_begin+step2;
 439         for(int dc1 = 0 ; dc1 < dn1 ; dc1++) {
 440             int const sc1_begin = (dc1<<step1_l2)-step1_2;
 441             int const sc1_end = sc1_begin+step1;
 442             unsigned int sum[PC];
 443             std::fill_n(sum, PC, 0);
 444             for(int sc2 = sc2_begin ; sc2 < sc2_end ; sc2++) {
 445                 for(int sc1 = sc1_begin ; sc1 < sc1_end ; sc1++) {
 446                     for(unsigned int ch = 0 ; ch < PC ; ch++) {
 447                         sum[ch] += src[clip(sc2,0,sn2-1)*sstr2+clip(sc1,0,sn1-1)*sstr1+ch];
 448                     }
 449                 }
 450             }
 451             for(unsigned int ch = 0 ; ch < PC ; ch++) {
 452                 dst[dc2*dstr2+dc1*dstr1+ch] = static_cast<PT>((sum[ch]+round_offset)>>divisor_l2);
 453             }
 454         }
 455     }
 456 }
 457
 458 template<typename PT, unsigned int PC>
 459 static void
 460 upsample(PT *const dst, int const dstr1, int const dstr2, unsigned int const dn1, unsigned int const dn2,
 461          PT const *const src, int const sstr1, int const sstr2, unsigned int const sn1, unsigned int const sn2,
 462          unsigned int const step1_l2, unsigned int const step2_l2)
 463 {
 464     assert(((sn1-1)<<step1_l2)>=dn1 && ((sn2-1)<<step2_l2)>=dn2); // The last pixel of the source image should fall outside the destination image
 465     unsigned int const divisor_l2 = step1_l2+step2_l2; // step1*step2=2^(step1_l2+step2_l2)
 466     unsigned int const round_offset = (1<<divisor_l2)/2;
 467     unsigned int const step1 = 1<<step1_l2;
 468     unsigned int const step2 = 1<<step2_l2;
 469     for ( unsigned int sc2 = 0 ; sc2 < sn2-1 ; sc2++ ) {
 470         unsigned int const dc2_begin = (sc2 << step2_l2);
 471         unsigned int const dc2_end = std::min(dn2, dc2_begin+step2);
 472         for ( unsigned int sc1 = 0 ; sc1 < sn1-1 ; sc1++ ) {
 473             unsigned int const dc1_begin = (sc1 << step1_l2);
 474             unsigned int const dc1_end = std::min(dn1, dc1_begin+step1);
 475             for ( unsigned int byte = 0 ; byte < PC ; byte++) {
 476
 477                 // get 4 values at the corners of the pixel from src
 478                 PT a00 = src[sstr2* sc2    + sstr1* sc1    + byte];
 479                 PT a10 = src[sstr2* sc2    + sstr1*(sc1+1) + byte];
 480                 PT a01 = src[sstr2*(sc2+1) + sstr1* sc1    + byte];
 481                 PT a11 = src[sstr2*(sc2+1) + sstr1*(sc1+1) + byte];
 482
 483                 // initialize values for linear interpolation
 484                 unsigned int a0 = a00*step2/*+a01*0*/;
 485                 unsigned int a1 = a10*step2/*+a11*0*/;
 486
 487                 // iterate over the rectangle to be interpolated
 488                 for ( unsigned int dc2 = dc2_begin ; dc2 < dc2_end ; dc2++ ) {
 489
 490                     // prepare linear interpolation for this row
 491                     unsigned int a = a0*step1/*+a1*0*/+round_offset;
 492
 493                     for ( unsigned int dc1 = dc1_begin ; dc1 < dc1_end ; dc1++ ) {
 494
 495                         // simple linear interpolation
 496                         dst[dstr2*dc2 + dstr1*dc1 + byte] = static_cast<PT>(a>>divisor_l2);
 497
 498                         // compute a = a0*(ix-1)+a1*(xi+1)+round_offset
 499                         a = a - a0 + a1;
 500                     }
 501
 502                     // compute a0 = a00*(iy-1)+a01*(yi+1) and similar for a1
 503                     a0 = a0 - a00 + a01;
 504                     a1 = a1 - a10 + a11;
 505                 }
 506             }
 507         }
 508     }
 509 }
 510
 511 int FilterGaussian::render(FilterSlot &slot, FilterUnits const &units)
 512 {
 513     /* in holds the input pixblock */
 514     NRPixBlock *in = slot.get(_input);
 515     if (!in) {
 516         g_warning("Missing source image for feGaussianBlur (in=%d)", _input);
 517         return 1;
 518     }
 519
 520     Matrix trans = units.get_matrix_primitiveunits2pb();
 521
 522     /* If to either direction, the standard deviation is zero or
 523      * input image is not defined,
 524      * a transparent black image should be returned. */
 525     if (_deviation_x <= 0 || _deviation_y <= 0 || in == NULL) {
 526         NRPixBlock *out = new NRPixBlock;
 527         if (in == NULL) {
 528             // A bit guessing here, but source graphic is likely to be of
 529             // right size
 530             in = slot.get(NR_FILTER_SOURCEGRAPHIC);
 531         }
 532         nr_pixblock_setup_fast(out, in->mode, in->area.x0, in->area.y0,
 533                                in->area.x1, in->area.y1, true);
 534         if (out->data.px != NULL) {
 535             out->empty = false;
 536             slot.set(_output, out);
 537         }
 538         return 0;
 539     }
 540
 541     // Some common constants
 542     int const width_org = in->area.x1-in->area.x0, height_org = in->area.y1-in->area.y0;
 543     double const deviation_x_org = _deviation_x * NR::expansionX(trans);
 544     double const deviation_y_org = _deviation_y * NR::expansionY(trans);
 545     int const PC = NR_PIXBLOCK_BPP(in);
 546
 547     // Subsampling constants
 548     Inkscape::Preferences *prefs = Inkscape::Preferences::get();
 549     int const quality = prefs->getInt("/options/blurquality/value");
 550     int const x_step_l2 = _effect_subsample_step_log2(deviation_x_org, quality);
 551     int const y_step_l2 = _effect_subsample_step_log2(deviation_y_org, quality);
 552     int const x_step = 1<<x_step_l2;
 553     int const y_step = 1<<y_step_l2;
 554     bool const resampling = x_step > 1 || y_step > 1;
 555     int const width = resampling ? static_cast<int>(ceil(static_cast<double>(width_org)/x_step))+1 : width_org;
 556     int const height = resampling ? static_cast<int>(ceil(static_cast<double>(height_org)/y_step))+1 : height_org;
 557     double const deviation_x = deviation_x_org / x_step;
 558     double const deviation_y = deviation_y_org / y_step;
 559     int const scr_len_x = _effect_area_scr(deviation_x);
 560     int const scr_len_y = _effect_area_scr(deviation_y);
 561
 562     // Decide which filter to use for X and Y
 563     // This threshold was determined by trial-and-error for one specific machine,
 564     // so there's a good chance that it's not optimal.
 565     // Whatever you do, don't go below 1 (and preferrably not even below 2), as
 566     // the IIR filter gets unstable there.
 567     bool const use_IIR_x = deviation_x > 3;
 568     bool const use_IIR_y = deviation_y > 3;
 569
 570     // new buffer for the subsampled output
 571     NRPixBlock *out = new NRPixBlock;
 572     nr_pixblock_setup_fast(out, in->mode, in->area.x0/x_step,       in->area.y0/y_step,
 573                                           in->area.x0/x_step+width, in->area.y0/y_step+height, true);
 574     if (out->size != NR_PIXBLOCK_SIZE_TINY && out->data.px == NULL) {
 575         // alas, we've accomplished a lot, but ran out of memory - so abort
 576         return 0;
 577     }
 578     // Temporary storage for IIR filter
 579     // NOTE: This can be eliminated, but it reduces the precision a bit
 580     IIRValue * tmpdata = 0;
 581     if ( use_IIR_x || use_IIR_y ) {
 582         tmpdata = new IIRValue[std::max(width,height)*PC];
 583         if (tmpdata == NULL) {
 584             nr_pixblock_release(out);
 585             delete out;
 586             return 0;
 587         }
 588     }
 589     NRPixBlock *ssin = in;
 590     if ( resampling ) {
 591         ssin = out;
 592         // Downsample
 593         switch(in->mode) {
 594         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
 595             downsample<unsigned char,1>(NR_PIXBLOCK_PX(out), 1, out->rs, width, height, NR_PIXBLOCK_PX(in), 1, in->rs, width_org, height_org, x_step_l2, y_step_l2);
 596             break;
 597         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
 598             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);
 599             break;
 600         case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
 601             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);
 602             break;
 603         case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
 604             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);
 605             break;
 606         default:
 607             assert(false);
 608         };
 609     }
 610
 611     if (use_IIR_x) {
 612         // Filter variables
 613         IIRValue b[N+1];  // scaling coefficient + filter coefficients (can be 10.21 fixed point)
 614         double bf[N];  // computed filter coefficients
 615         double M[N*N]; // matrix used for initialization procedure (has to be double)
 616
 617         // Compute filter (x)
 618         calcFilter(deviation_x, bf);
 619         for(size_t i=0; i<N; i++) bf[i] = -bf[i];
 620         b[0] = 1; // b[0] == alpha (scaling coefficient)
 621         for(size_t i=0; i<N; i++) {
 622             b[i+1] = bf[i];
 623             b[0] -= b[i+1];
 624         }
 625
 626         // Compute initialization matrix (x)
 627         calcTriggsSdikaM(bf, M);
 628
 629         // Filter (x)
 630         switch(in->mode) {
 631         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
 632             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);
 633             break;
 634         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
 635             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);
 636             break;
 637         case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
 638             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);
 639             break;
 640         case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
 641             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);
 642             break;
 643         default:
 644             assert(false);
 645         };
 646     } else { // !use_IIR_x
 647         // Filter kernel for x direction
 648         FIRValue kernel[scr_len_x];
 649         _make_kernel(kernel, deviation_x);
 650
 651         // Filter (x)
 652         switch(in->mode) {
 653         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
 654             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);
 655             break;
 656         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
 657             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);
 658             break;
 659         case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
 660             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);
 661             break;
 662         case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
 663             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);
 664             break;
 665         default:
 666             assert(false);
 667         };
 668     }
 669
 670     if (use_IIR_y) {
 671         // Filter variables
 672         IIRValue b[N+1];  // scaling coefficient + filter coefficients (can be 10.21 fixed point)
 673         double bf[N];  // computed filter coefficients
 674         double M[N*N]; // matrix used for initialization procedure (has to be double)
 675
 676         // Compute filter (y)
 677         calcFilter(deviation_y, bf);
 678         for(size_t i=0; i<N; i++) bf[i] = -bf[i];
 679         b[0] = 1; // b[0] == alpha (scaling coefficient)
 680         for(size_t i=0; i<N; i++) {
 681             b[i+1] = bf[i];
 682             b[0] -= b[i+1];
 683         }
 684
 685         // Compute initialization matrix (y)
 686         calcTriggsSdikaM(bf, M);
 687
 688         // Filter (y)
 689         switch(in->mode) {
 690         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
 691             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);
 692             break;
 693         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
 694             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);
 695             break;
 696         case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
 697             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);
 698             break;
 699         case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
 700             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);
 701             break;
 702         default:
 703             assert(false);
 704         };
 705     } else { // !use_IIR_y
 706         // Filter kernel for y direction
 707         FIRValue kernel[scr_len_y];
 708         _make_kernel(kernel, deviation_y);
 709
 710         // Filter (y)
 711         switch(in->mode) {
 712         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
 713             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);
 714             break;
 715         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
 716             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);
 717             break;
 718         case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
 719             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);
 720             break;
 721         case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
 722             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);
 723             break;
 724         default:
 725             assert(false);
 726         };
 727     }
 728
 729     delete[] tmpdata; // deleting a nullptr has no effect, so this is save
 730
 731     if ( !resampling ) {
 732         // No upsampling needed
 733         out->empty = FALSE;
 734         slot.set(_output, out);
 735     } else {
 736         // New buffer for the final output, same resolution as the in buffer
 737         NRPixBlock *finalout = new NRPixBlock;
 738         nr_pixblock_setup_fast(finalout, in->mode, in->area.x0, in->area.y0,
 739                                                    in->area.x1, in->area.y1, true);
 740         if (finalout->size != NR_PIXBLOCK_SIZE_TINY && finalout->data.px == NULL) {
 741             // alas, we've accomplished a lot, but ran out of memory - so abort
 742             nr_pixblock_release(out);
 743             delete out;
 744             return 0;
 745         }
 746
 747         // Upsample
 748         switch(in->mode) {
 749         case NR_PIXBLOCK_MODE_A8:        ///< Grayscale
 750             upsample<unsigned char,1>(NR_PIXBLOCK_PX(finalout), 1, finalout->rs, width_org, height_org, NR_PIXBLOCK_PX(out), 1, out->rs, width, height, x_step_l2, y_step_l2);
 751             break;
 752         case NR_PIXBLOCK_MODE_R8G8B8:    ///< 8 bit RGB
 753             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);
 754             break;
 755         case NR_PIXBLOCK_MODE_R8G8B8A8N: ///< Normal 8 bit RGBA
 756             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);
 757             break;
 758         case NR_PIXBLOCK_MODE_R8G8B8A8P:  ///< Premultiplied 8 bit RGBA
 759             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);
 760             break;
 761         default:
 762             assert(false);
 763         };
 764
 765         // We don't need the out buffer anymore
 766         nr_pixblock_release(out);
 767         delete out;
 768
 769         // The final out buffer gets returned
 770         finalout->empty = FALSE;
 771         slot.set(_output, finalout);
 772     }
 773
 774     return 0;
 775 }
 776
 777 void FilterGaussian::area_enlarge(NRRectL &area, Matrix const &trans)
 778 {
 779     int area_x = _effect_area_scr(_deviation_x * NR::expansionX(trans));
 780     int area_y = _effect_area_scr(_deviation_y * NR::expansionY(trans));
 781     // maximum is used because rotations can mix up these directions
 782     // TODO: calculate a more tight-fitting rendering area
 783     int area_max = std::max(area_x, area_y);
 784     area.x0 -= area_max;
 785     area.x1 += area_max;
 786     area.y0 -= area_max;
 787     area.y1 += area_max;
 788 }
 789
 790 FilterTraits FilterGaussian::get_input_traits() {
 791     return TRAIT_PARALLER;
 792 }
 793
 794 void FilterGaussian::set_deviation(double deviation)
 795 {
 796     if(IS_FINITE(deviation) && deviation >= 0) {
 797         _deviation_x = _deviation_y = deviation;
 798     }
 799 }
 800
 801 void FilterGaussian::set_deviation(double x, double y)
 802 {
 803     if(IS_FINITE(x) && x >= 0 && IS_FINITE(y) && y >= 0) {
 804         _deviation_x = x;
 805         _deviation_y = y;
 806     }
 807 }
 808
 809 } /* namespace NR */
 810
 811 /*
 812   Local Variables:
 813   mode:c++
 814   c-file-style:"stroustrup"
 815   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 816   indent-tabs-mode:nil
 817   fill-column:99
 818   End:
 819 */
 820 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :