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