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