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