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