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