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