1 #define __NR_FILTER_GAUSSIAN_CPP__
3 /*
4 * Gaussian blur renderer
5 *
6 * Author:
7 * Niko Kiirala <niko@kiirala.com>
8 *
9 * Copyright (C) 2006 Niko Kiirala
10 *
11 * Released under GNU GPL, read the file 'COPYING' for more information
12 */
14 #include <cmath>
16 using std::isnormal;
18 #include "display/nr-filter-primitive.h"
19 #include "display/nr-filter-gaussian.h"
20 #include "libnr/nr-pixblock.h"
21 #include "libnr/nr-matrix.h"
22 #include "prefs-utils.h"
24 namespace NR {
26 FilterGaussian::FilterGaussian()
27 {
28 _deviation_x = _deviation_y = prefs_get_double_attribute("options.filtertest", "value", 0.0);
29 }
31 FilterPrimitive *FilterGaussian::create()
32 {
33 return new FilterGaussian();
34 }
36 int FilterGaussian::_kernel_size(Matrix const &trans)
37 {
38 int length_x = _effect_area_scr_x(trans);
39 int length_y = _effect_area_scr_y(trans);
40 return _max(length_x, length_y) * 2 + 1;
41 }
43 void FilterGaussian::_make_kernel(double *kernel, double deviation, double expansion)
44 {
45 double length = deviation * 3.0;
46 int scr_len = (int)std::floor(length * expansion);
47 if(scr_len < 1) scr_len = 1;
48 double d_sq = deviation * deviation * 2;
49 double step = length / scr_len;
51 double sum = 0;
52 for ( int i = 0; i < scr_len * 2 + 1 ; i++ ) {
53 double i_sq = (step * i - length) * (step * i - length);
54 sum += (std::exp(-i_sq / d_sq) / std::sqrt(M_PI * d_sq));
55 }
57 for ( int i = 0; i < scr_len * 2 + 1 ; i++ ) {
58 double i_sq = (step * i - length) * (step * i - length);
59 kernel[i] = (std::exp(-i_sq / d_sq) / std::sqrt(M_PI * d_sq)) / sum;
60 }
61 }
63 int FilterGaussian::_effect_area_scr_x(Matrix const &trans)
64 {
65 int ret = (int)std::floor(_deviation_x * 3.0 * trans.expansionX());
66 if(ret < 1) ret = 1;
67 return ret;
68 }
70 int FilterGaussian::_effect_area_scr_y(Matrix const &trans)
71 {
72 int ret = (int)std::floor(_deviation_y * 3.0 * trans.expansionY());
73 if(ret < 1) ret = 1;
74 return ret;
75 }
77 int FilterGaussian::_effect_subsample_step(int scr_len_x)
78 {
79 if (scr_len_x < 16) {
80 return 1;
81 } else if (scr_len_x < 80) {
82 return 4;
83 } else if (scr_len_x < 160) {
84 return 8;
85 } else if (scr_len_x < 320) {
86 return 32;
87 } else if (scr_len_x < 640) {
88 return 64;
89 } else if (scr_len_x < 1280) {
90 return 256;
91 } else if (scr_len_x < 2560) {
92 return 1024;
93 } else {
94 return 65536;
95 }
96 }
98 int FilterGaussian::_effect_subsample_step_log2(int scr_len_x)
99 {
100 if (scr_len_x < 16) {
101 return 0;
102 } else if (scr_len_x < 80) {
103 return 2;
104 } else if (scr_len_x < 160) {
105 return 3;
106 } else if (scr_len_x < 320) {
107 return 5;
108 } else if (scr_len_x < 640) {
109 return 6;
110 } else if (scr_len_x < 1280) {
111 return 8;
112 } else if (scr_len_x < 2560) {
113 return 10;
114 } else {
115 return 16;
116 }
117 }
120 int FilterGaussian::render(NRPixBlock **pb, Matrix const &trans)
121 {
122 /* in holds the input pixblock */
123 NRPixBlock *in = pb[0];
125 /* Blur radius in screen units (pixels) */
126 int scr_len_x = _effect_area_scr_x(trans);
127 int scr_len_y = _effect_area_scr_y(trans);
129 // subsampling step; it depends on the radius, but somewhat nonlinearly, to make high zooms
130 // workable
131 int stepx = _effect_subsample_step(scr_len_x);
132 int stepx_l2 = _effect_subsample_step_log2(scr_len_x);
133 int stepy = _effect_subsample_step(scr_len_y);
134 int stepy_l2 = _effect_subsample_step_log2(scr_len_y);
135 int stepx2 = stepx >> 1;
136 int stepy2 = stepy >> 1;
138 /* buffer for x-axis blur */
139 NRPixBlock *bufx = new NRPixBlock;
140 /* buffer for y-axis blur */
141 NRPixBlock *bufy = new NRPixBlock;
143 // boundaries of the subsampled (smaller, unless step==1) buffers
144 int xd0 = (in->area.x0 >> stepx_l2);
145 int xd1 = (in->area.x1 >> stepx_l2) + 1;
146 int yd0 = (in->area.y0 >> stepy_l2);
147 int yd1 = (in->area.y1 >> stepy_l2) + 1;
149 // set up subsampled buffers
150 nr_pixblock_setup_fast(bufx, in->mode, xd0, yd0, xd1, yd1, true);
151 nr_pixblock_setup_fast(bufy, in->mode, xd0, yd0, xd1, yd1, true);
153 //mid->visible_area = in->visible_area;
154 //out->visible_area = in->visible_area;
156 /* Array for filter kernel, big enough to fit kernels for both X and Y
157 * direction kernel, one at time */
158 double kernel[_kernel_size(trans)];
160 /* 1. Blur in direction of X-axis, from in to bufx (they have different resolution)*/
161 _make_kernel(kernel, _deviation_x, trans.expansionX());
163 for ( int y = bufx->area.y0 ; y < bufx->area.y1; y++ ) {
165 // corresponding line in the source buffer
166 int in_line;
167 if ((y << stepy_l2) >= in->area.y1) {
168 in_line = (in->area.y1 - in->area.y0 - 1) * in->rs;
169 } else {
170 in_line = ((y << stepy_l2) - (in->area.y0)) * in->rs;
171 if (in_line < 0)
172 in_line = 0;
173 }
175 // current line in bufx
176 int bufx_line = (y - yd0) * bufx->rs;
178 int skipbuf[4] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN};
180 for ( int x = bufx->area.x0 ; x < bufx->area.x1 ; x++ ) {
182 // for all bytes of the pixel
183 for ( int byte = 0 ; byte < NR_PIXBLOCK_BPP(in) ; byte++) {
185 if(skipbuf[byte] > x) continue;
187 double sum = 0;
188 int last_in = -1;
189 int different_count = 0;
191 // go over our point's neighborhood on x axis in the in buffer, with stepx increment
192 for ( int i = -scr_len_x ; i <= scr_len_x ; i += stepx ) {
194 // the pixel we're looking at
195 int x_in = (((x << stepx_l2) + i + stepx2) >> stepx_l2) << stepx_l2;
197 // distance from it to the current x,y
198 int dist = x_in - (x << stepx_l2);
199 if (dist < -scr_len_x)
200 dist = -scr_len_x;
201 if (dist > scr_len_x)
202 dist = scr_len_x;
204 if (x_in >= in->area.x1) {
205 x_in = (in->area.x1 - in->area.x0 - 1);
206 } else {
207 x_in = (x_in - in->area.x0);
208 if (x_in < 0)
209 x_in = 0;
210 }
212 // value at the pixel
213 unsigned char in_byte = NR_PIXBLOCK_PX(in)[in_line + NR_PIXBLOCK_BPP(in) * x_in + byte];
215 // is it the same as last one we saw?
216 if(in_byte != last_in) different_count++;
217 last_in = in_byte;
219 // sum pixels weighted by the kernel; multiply by stepx because we're skipping stepx pixels
220 sum += stepx * in_byte * kernel[scr_len_x + dist];
221 }
223 // store the result in bufx
224 NR_PIXBLOCK_PX(bufx)[bufx_line + NR_PIXBLOCK_BPP(bufx) * (x - xd0) + byte] = (unsigned char)sum;
226 // optimization: if there was no variation within this point's neighborhood,
227 // skip ahead while we keep seeing the same last_in byte:
228 // blurring flat color would not change it anyway
229 if (different_count <= 1) {
230 int pos = x + 1;
231 while(((pos << stepx_l2) + scr_len_x) < in->area.x1 &&
232 NR_PIXBLOCK_PX(in)[in_line + NR_PIXBLOCK_BPP(in) * ((pos << stepx_l2) + scr_len_x - in->area.x0) + byte] == last_in)
233 {
234 NR_PIXBLOCK_PX(bufx)[bufx_line + NR_PIXBLOCK_BPP(bufx) * (pos - xd0) + byte] = last_in;
235 pos++;
236 }
237 skipbuf[byte] = pos;
238 }
239 }
240 }
241 }
244 /* 2. Blur in direction of Y-axis, from bufx to bufy (they have the same resolution) */
245 _make_kernel(kernel, _deviation_y, trans.expansionY());
247 for ( int x = bufy->area.x0 ; x < bufy->area.x1; x++ ) {
249 int bufy_disp = NR_PIXBLOCK_BPP(bufy) * (x - xd0);
250 int bufx_disp = NR_PIXBLOCK_BPP(bufx) * (x - xd0);
252 int skipbuf[4] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN};
254 for ( int y = bufy->area.y0; y < bufy->area.y1; y++ ) {
256 int bufy_line = (y - yd0) * bufy->rs;
258 for ( int byte = 0 ; byte < NR_PIXBLOCK_BPP(bufx) ; byte++) {
260 if (skipbuf[byte] > y) continue;
262 double sum = 0;
263 int last_in = -1;
264 int different_count = 0;
266 for ( int i = -scr_len_y ; i <= scr_len_y ; i += stepy ) {
268 int y_in = ((((y << stepy_l2) + i + stepy2) >> stepy_l2) - yd0);
270 int dist = ((y_in + yd0) << stepy_l2) - (y << stepy_l2);
271 if (dist < -scr_len_y)
272 dist = -scr_len_y;
273 if (dist > scr_len_y)
274 dist = scr_len_y;
276 if (y_in > (yd1 - yd0)) y_in = (yd1 - yd0);
277 if (y_in < 0) y_in = 0;
279 unsigned char in_byte = NR_PIXBLOCK_PX(bufx)[y_in * bufx->rs + NR_PIXBLOCK_BPP(bufx) * (x - xd0) + byte];
280 if(in_byte != last_in) different_count++;
281 last_in = in_byte;
282 sum += stepy * in_byte * kernel[scr_len_y + dist];
283 }
285 NR_PIXBLOCK_PX(bufy)[bufy_line + bufy_disp + byte] = (unsigned char)sum;
287 if (different_count <= 1) {
288 int pos = y + 1;
289 while((pos + (scr_len_y >> stepy_l2) + 1) < yd1 &&
290 NR_PIXBLOCK_PX(bufx)[(pos + (scr_len_y >> stepy_l2) + 1 - yd0) * bufx->rs + bufx_disp + byte] == last_in)
291 {
292 NR_PIXBLOCK_PX(bufy)[(pos - yd0) * bufy->rs + bufy_disp + byte] = last_in;
293 pos++;
294 }
295 skipbuf[byte] = pos;
296 }
298 }
299 }
300 }
302 // we don't need bufx anymore
303 nr_pixblock_release(bufx);
304 delete bufx;
306 // interpolation will need to divide by stepx * stepy
307 int divisor = stepx_l2 + stepy_l2;
309 // new buffer for the final output, same resolution as the in buffer
310 NRPixBlock *out = new NRPixBlock;
311 nr_pixblock_setup_fast(out, in->mode, in->area.x0, in->area.y0,
312 in->area.x1, in->area.y1, true);
314 for ( int y = yd0 ; y < yd1 - 1; y++ ) {
315 for ( int x = xd0 ; x < xd1 - 1; x++ ) {
316 for ( int byte = 0 ; byte < NR_PIXBLOCK_BPP(bufy) ; byte++) {
318 // get 4 values at the corners of the pixel from bufy
319 unsigned char a00 = NR_PIXBLOCK_PX(bufy)[((y - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) * (x - xd0) + byte];
320 if (stepx == 1 && stepy == 1) { // if there was no subsampling, just use a00
321 NR_PIXBLOCK_PX(out)[((y - yd0) * out->rs) + NR_PIXBLOCK_BPP(out) * (x - xd0) + byte] = a00;
322 continue;
323 }
324 unsigned char a10 = NR_PIXBLOCK_PX(bufy)[((y - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) * (x + 1 - xd0) + byte];
325 unsigned char a01 = NR_PIXBLOCK_PX(bufy)[((y + 1 - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) * (x - xd0) + byte];
326 unsigned char a11 = NR_PIXBLOCK_PX(bufy)[((y + 1 - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) * (x + 1 - xd0) + byte];
328 // iterate over the rectangle to be interpolated
329 for ( int yi = 0 ; yi < stepy; yi++ ) {
330 int iy = stepy - yi;
331 int y_out = (y << stepy_l2) + yi;
332 if ((y_out < out->area.y0) || (y_out >= out->area.y1))
333 continue;
334 int out_line = (y_out - out->area.y0) * out->rs;
336 for ( int xi = 0 ; xi < stepx; xi++ ) {
337 int ix = stepx - xi;
338 int x_out = (x << stepx_l2) + xi;
339 if ((x_out < out->area.x0) || (x_out >= out->area.x1))
340 continue;
342 // simple linear interpolation
343 int a = (a00*ix*iy + a10*xi*iy + a01*ix*yi + a11*xi*yi) >> divisor;
345 NR_PIXBLOCK_PX(out)[out_line + NR_PIXBLOCK_BPP(out) * (x_out - out->area.x0) + byte] = (unsigned char) a;
346 }
347 }
348 }
349 }
350 }
352 nr_pixblock_release(bufy);
353 delete bufy;
355 out->empty = FALSE;
356 pb[1] = out;
358 return 0;
359 }
361 int FilterGaussian::get_enlarge(Matrix const &trans)
362 {
363 int area_x = _effect_area_scr_x(trans);
364 int area_y = _effect_area_scr_y(trans);
365 return _max(area_x, area_y);
366 }
368 void FilterGaussian::set_deviation(double deviation)
369 {
370 if(isnormal(deviation) && deviation >= 0) {
371 _deviation_x = _deviation_y = deviation;
372 }
373 }
375 void FilterGaussian::set_deviation(double x, double y)
376 {
377 if(isnormal(x) && x >= 0 && isnormal(y) && y >= 0) {
378 _deviation_x = x;
379 _deviation_y = y;
380 }
381 }
383 } /* namespace NR */
385 /*
386 Local Variables:
387 mode:c++
388 c-file-style:"stroustrup"
389 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
390 indent-tabs-mode:nil
391 fill-column:99
392 End:
393 */
394 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :