src/display/nr-filter-gaussian.cpp

   1 #define __NR_FILTER_GAUSSIAN_CPP__
   2
   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  */
  13
  14 #include <cmath>
  15 #include <glib.h>
  16
  17 using std::isnormal;
  18
  19 #include "display/nr-filter-primitive.h"
  20 #include "display/nr-filter-gaussian.h"
  21 #include "display/nr-filter-types.h"
  22 #include "libnr/nr-pixblock.h"
  23 #include "libnr/nr-matrix.h"
  24 #include "prefs-utils.h"
  25
  26 namespace NR {
  27
  28 FilterGaussian::FilterGaussian()
  29 {
  30     _deviation_x = _deviation_y = prefs_get_double_attribute("options.filtertest", "value", 0.0);
  31 }
  32
  33 FilterPrimitive *FilterGaussian::create()
  34 {
  35     return new FilterGaussian();
  36 }
  37
  38 FilterGaussian::~FilterGaussian()
  39 {
  40     // Nothing to do here
  41 }
  42
  43 int FilterGaussian::_kernel_size(Matrix const &trans)
  44 {
  45     int length_x = _effect_area_scr_x(trans);
  46     int length_y = _effect_area_scr_y(trans);
  47     return _max(length_x, length_y) * 2 + 1;
  48 }
  49
  50 void FilterGaussian::_make_kernel(double *kernel, double deviation, double expansion)
  51 {
  52     double length = deviation * 3.0;
  53     int scr_len = (int)std::floor(length * expansion);
  54     if(scr_len < 1) scr_len = 1;
  55     double d_sq = deviation * deviation * 2;
  56     double step = length / scr_len;
  57
  58     double sum = 0;
  59     for ( int i = 0; i < scr_len * 2 + 1 ; i++ ) {
  60         double i_sq = (step * i - length) * (step * i - length);
  61         sum += (std::exp(-i_sq / d_sq) / std::sqrt(M_PI * d_sq));
  62     }
  63
  64     for ( int i = 0; i < scr_len * 2 + 1 ; i++ ) {
  65         double i_sq = (step * i - length) * (step * i - length);
  66         kernel[i] = (std::exp(-i_sq / d_sq) / std::sqrt(M_PI * d_sq)) / sum;
  67     }
  68 }
  69
  70 int FilterGaussian::_effect_area_scr_x(Matrix const &trans)
  71 {
  72     int ret = (int)std::floor(_deviation_x * 3.0 * trans.expansionX());
  73     if(ret < 1) ret = 1;
  74     return ret;
  75 }
  76
  77 int FilterGaussian::_effect_area_scr_y(Matrix const &trans)
  78 {
  79     int ret = (int)std::floor(_deviation_y * 3.0 * trans.expansionY());
  80     if(ret < 1) ret = 1;
  81     return ret;
  82 }
  83
  84 int FilterGaussian::_effect_subsample_step(int scr_len_x)
  85 {
  86     if (scr_len_x < 16) {
  87         return 1;
  88     } else if (scr_len_x < 80) {
  89         return 4;
  90     } else if (scr_len_x < 160) {
  91         return 8;
  92     } else if (scr_len_x < 320) {
  93         return 32;
  94     } else if (scr_len_x < 640) {
  95         return 64;
  96     } else if (scr_len_x < 1280) {
  97         return 256;
  98     } else if (scr_len_x < 2560) {
  99         return 1024;
 100     } else {
 101         return 65536;
 102     }
 103 }
 104
 105 int FilterGaussian::_effect_subsample_step_log2(int scr_len_x)
 106 {
 107     if (scr_len_x < 16) {
 108         return 0;
 109     } else if (scr_len_x < 80) {
 110         return 2;
 111     } else if (scr_len_x < 160) {
 112         return 3;
 113     } else if (scr_len_x < 320) {
 114         return 5;
 115     } else if (scr_len_x < 640) {
 116         return 6;
 117     } else if (scr_len_x < 1280) {
 118         return 8;
 119     } else if (scr_len_x < 2560) {
 120         return 10;
 121     } else {
 122         return 16;
 123     }
 124 }
 125
 126 /**
 127  * Sanity check function for indexing pixblocks.
 128  * Catches reading and writing outside the pixblock area.
 129  * When enabled, decreases filter rendering speed massively.
 130  */
 131 inline void _check_index(NRPixBlock const * const pb, int const location, int const line)
 132 {
 133     if(true) {
 134         int max_loc = pb->rs * (pb->area.y1 - pb->area.y0);
 135         if (location < 0 || location >= max_loc)
 136             g_warning("Location %d out of bounds (0 ... %d) at line %d", location, max_loc, line);
 137     }
 138 }
 139
 140 int FilterGaussian::render(FilterSlot &slot, Matrix const &trans)
 141 {
 142     /* in holds the input pixblock */
 143     NRPixBlock *in = slot.get(_input);
 144
 145     /* If to either direction, the standard deviation is zero, a transparent
 146      * black image should be returned */
 147     if (_deviation_x <= 0 || _deviation_y <= 0) {
 148         NRPixBlock *out = new NRPixBlock;
 149         nr_pixblock_setup_fast(out, in->mode, in->area.x0, in->area.y0,
 150                                in->area.x1, in->area.y1, true);
 151         out->empty = false;
 152         slot.set(_output, out);
 153         return 0;
 154     }
 155
 156     /* Blur radius in screen units (pixels) */
 157     int scr_len_x = _effect_area_scr_x(trans);
 158     int scr_len_y = _effect_area_scr_y(trans);
 159
 160     // subsampling step; it depends on the radius, but somewhat nonlinearly, to make high zooms
 161     // workable
 162     int stepx = _effect_subsample_step(scr_len_x);
 163     int stepx_l2 = _effect_subsample_step_log2(scr_len_x);
 164     int stepy = _effect_subsample_step(scr_len_y);
 165     int stepy_l2 = _effect_subsample_step_log2(scr_len_y);
 166     int stepx2 = stepx >> 1;
 167     int stepy2 = stepy >> 1;
 168
 169     /* buffer for x-axis blur */
 170     NRPixBlock *bufx = new NRPixBlock;
 171     /* buffer for y-axis blur */
 172     NRPixBlock *bufy = new NRPixBlock;
 173
 174     // boundaries of the subsampled (smaller, unless step==1) buffers
 175     int xd0 = (in->area.x0 >> stepx_l2);
 176     int xd1 = (in->area.x1 >> stepx_l2) + 1;
 177     int yd0 = (in->area.y0 >> stepy_l2);
 178     int yd1 = (in->area.y1 >> stepy_l2) + 1;
 179
 180     // set up subsampled buffers
 181     nr_pixblock_setup_fast(bufx, in->mode, xd0, yd0, xd1, yd1, true);
 182     nr_pixblock_setup_fast(bufy, in->mode, xd0, yd0, xd1, yd1, true);
 183
 184     //mid->visible_area = in->visible_area;
 185     //out->visible_area = in->visible_area;
 186
 187     /* Array for filter kernel, big enough to fit kernels for both X and Y
 188      * direction kernel, one at time */
 189     double kernel[_kernel_size(trans)];
 190
 191     /* 1. Blur in direction of X-axis, from in to bufx (they have different resolution)*/
 192     _make_kernel(kernel, _deviation_x, trans.expansionX());
 193
 194     for ( int y = bufx->area.y0 ; y < bufx->area.y1; y++ ) {
 195
 196         // corresponding line in the source buffer
 197         int in_line;
 198         if ((y << stepy_l2) >= in->area.y1) {
 199             in_line = (in->area.y1 - in->area.y0 - 1) * in->rs;
 200         } else {
 201             in_line = ((y << stepy_l2) - (in->area.y0))  * in->rs;
 202             if (in_line < 0)
 203                 in_line = 0;
 204         }
 205
 206         // current line in bufx
 207         int bufx_line = (y - yd0) * bufx->rs;
 208
 209         int skipbuf[4] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN};
 210
 211         for ( int x = bufx->area.x0 ; x < bufx->area.x1 ; x++ ) {
 212
 213             // for all bytes of the pixel
 214             for ( int byte = 0 ; byte < NR_PIXBLOCK_BPP(in) ; byte++) {
 215
 216                 if(skipbuf[byte] > x) continue;
 217
 218                 double sum = 0;
 219                 int last_in = -1;
 220                 int different_count = 0;
 221
 222                 // go over our point's neighborhood on x axis in the in buffer, with stepx increment
 223                 for ( int i = -scr_len_x ; i <= scr_len_x ; i += stepx ) {
 224
 225                     // the pixel we're looking at
 226                     int x_in = (((x << stepx_l2) + i + stepx2) >> stepx_l2) << stepx_l2;
 227
 228                     // distance from it to the current x,y
 229                     int dist = x_in - (x << stepx_l2);
 230                     if (dist < -scr_len_x)
 231                         dist = -scr_len_x;
 232                     if (dist > scr_len_x)
 233                         dist = scr_len_x;
 234
 235                     if (x_in >= in->area.x1) {
 236                         x_in = (in->area.x1 - in->area.x0 - 1);
 237                     } else {
 238                         x_in = (x_in - in->area.x0);
 239                         if (x_in < 0)
 240                             x_in = 0;
 241                     }
 242
 243                     // value at the pixel
 244                     _check_index(in, in_line + NR_PIXBLOCK_BPP(in) * x_in + byte, __LINE__);
 245                     unsigned char in_byte = NR_PIXBLOCK_PX(in)[in_line + NR_PIXBLOCK_BPP(in) * x_in + byte];
 246
 247                     // is it the same as last one we saw?
 248                     if(in_byte != last_in) different_count++;
 249                     last_in = in_byte;
 250
 251                     // sum pixels weighted by the kernel; multiply by stepx because we're skipping stepx pixels
 252                     sum += stepx * in_byte * kernel[scr_len_x + dist];
 253                 }
 254
 255                 // store the result in bufx
 256                 _check_index(bufx, bufx_line + NR_PIXBLOCK_BPP(bufx) * (x - xd0) + byte, __LINE__);
 257                 NR_PIXBLOCK_PX(bufx)[bufx_line + NR_PIXBLOCK_BPP(bufx) * (x - xd0) + byte] = (unsigned char)sum;
 258
 259                 // optimization: if there was no variation within this point's neighborhood,
 260                 // skip ahead while we keep seeing the same last_in byte:
 261                 // blurring flat color would not change it anyway
 262                 if (different_count <= 1) {
 263                     int pos = x + 1;
 264                     while(((pos << stepx_l2) + scr_len_x) < in->area.x1 &&
 265                           NR_PIXBLOCK_PX(in)[in_line + NR_PIXBLOCK_BPP(in) * ((pos << stepx_l2) + scr_len_x - in->area.x0) + byte] == last_in)
 266                     {
 267                         _check_index(in, in_line + NR_PIXBLOCK_BPP(in) * ((pos << stepx_l2) + scr_len_x - in->area.x0) + byte, __LINE__);
 268                         _check_index(bufx, bufx_line + NR_PIXBLOCK_BPP(bufx) * (pos - xd0) + byte, __LINE__);
 269                         NR_PIXBLOCK_PX(bufx)[bufx_line + NR_PIXBLOCK_BPP(bufx) * (pos - xd0) + byte] = last_in;
 270                         pos++;
 271                     }
 272                     skipbuf[byte] = pos;
 273                 }
 274             }
 275         }
 276     }
 277
 278
 279     /* 2. Blur in direction of Y-axis, from bufx to bufy (they have the same resolution) */
 280     _make_kernel(kernel, _deviation_y, trans.expansionY());
 281
 282     for ( int x = bufy->area.x0 ; x < bufy->area.x1; x++ ) {
 283
 284         int bufy_disp = NR_PIXBLOCK_BPP(bufy) * (x - xd0);
 285         int bufx_disp = NR_PIXBLOCK_BPP(bufx) * (x - xd0);
 286
 287         int skipbuf[4] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN};
 288
 289         for ( int y = bufy->area.y0; y < bufy->area.y1; y++ ) {
 290
 291             int bufy_line = (y - yd0) * bufy->rs;
 292
 293             for ( int byte = 0 ; byte < NR_PIXBLOCK_BPP(bufx) ; byte++) {
 294
 295                 if (skipbuf[byte] > y) continue;
 296
 297                 double sum = 0;
 298                 int last_in = -1;
 299                 int different_count = 0;
 300
 301                 for ( int i = -scr_len_y ; i <= scr_len_y ; i += stepy ) {
 302
 303                     int y_in = ((((y << stepy_l2) + i + stepy2) >> stepy_l2) - yd0);
 304
 305                     int dist = ((y_in + yd0) << stepy_l2) - (y << stepy_l2);
 306                     if (dist < -scr_len_y)
 307                         dist = -scr_len_y;
 308                     if (dist > scr_len_y)
 309                         dist = scr_len_y;
 310
 311                     if (y_in >= (yd1 - yd0)) y_in = (yd1 - yd0) - 1;
 312                     if (y_in < 0) y_in = 0;
 313
 314                     _check_index(bufx, y_in * bufx->rs + NR_PIXBLOCK_BPP(bufx) * (x - xd0) + byte, __LINE__);
 315                     unsigned char in_byte = NR_PIXBLOCK_PX(bufx)[y_in * bufx->rs + NR_PIXBLOCK_BPP(bufx) * (x - xd0) + byte];
 316                     if(in_byte != last_in) different_count++;
 317                     last_in = in_byte;
 318                     sum += stepy * in_byte * kernel[scr_len_y + dist];
 319                 }
 320
 321                 _check_index(bufy, bufy_line + bufy_disp + byte, __LINE__);
 322                 NR_PIXBLOCK_PX(bufy)[bufy_line + bufy_disp + byte] = (unsigned char)sum;
 323
 324                 if (different_count <= 1) {
 325                     int pos = y + 1;
 326                     while((pos + (scr_len_y >> stepy_l2) + 1) < yd1 &&
 327                           NR_PIXBLOCK_PX(bufx)[(pos + (scr_len_y >> stepy_l2) + 1 - yd0) * bufx->rs + bufx_disp + byte] == last_in)
 328                     {
 329                         _check_index(bufx, (pos + (scr_len_y >> stepy_l2) + 1 - yd0) * bufx->rs + bufx_disp + byte, __LINE__);
 330                         _check_index(bufy, (pos - yd0) * bufy->rs + bufy_disp + byte, __LINE__);
 331                         NR_PIXBLOCK_PX(bufy)[(pos - yd0) * bufy->rs + bufy_disp + byte] = last_in;
 332                         pos++;
 333                     }
 334                     skipbuf[byte] = pos;
 335                 }
 336
 337             }
 338         }
 339     }
 340
 341     // we don't need bufx anymore
 342     nr_pixblock_release(bufx);
 343     delete bufx;
 344
 345     // interpolation will need to divide by stepx * stepy
 346     int divisor = stepx_l2 + stepy_l2;
 347
 348     // new buffer for the final output, same resolution as the in buffer
 349     NRPixBlock *out = new NRPixBlock;
 350     nr_pixblock_setup_fast(out, in->mode, in->area.x0, in->area.y0,
 351                            in->area.x1, in->area.y1, true);
 352
 353     for ( int y = yd0 ; y < yd1 - 1; y++ ) {
 354         for ( int x = xd0 ; x < xd1 - 1; x++ ) {
 355             for ( int byte = 0 ; byte < NR_PIXBLOCK_BPP(bufy) ; byte++) {
 356
 357                 // get 4 values at the corners of the pixel from bufy
 358                 _check_index(bufy, ((y - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) + (x - xd0) + byte, __LINE__);
 359                 unsigned char a00 = NR_PIXBLOCK_PX(bufy)[((y - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) * (x - xd0) + byte];
 360                 if (stepx == 1 && stepy == 1) { // if there was no subsampling, just use a00
 361                     _check_index(out, ((y - yd0) * out->rs) + NR_PIXBLOCK_BPP(out) * (x - xd0) + byte, __LINE__);
 362                     NR_PIXBLOCK_PX(out)[((y - yd0) * out->rs) + NR_PIXBLOCK_BPP(out) * (x - xd0) + byte] = a00;
 363                     continue;
 364                 }
 365                 _check_index(bufy, ((y - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) + (x + 1 - xd0) + byte, __LINE__);
 366                 unsigned char a10 = NR_PIXBLOCK_PX(bufy)[((y - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) * (x + 1 - xd0) + byte];
 367                 _check_index(bufy, ((y + 1 - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) + (x - xd0) + byte, __LINE__);
 368                 unsigned char a01 = NR_PIXBLOCK_PX(bufy)[((y + 1 - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) * (x - xd0) + byte];
 369                 _check_index(bufy, ((y + 1 - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) + (x + 1 - xd0) + byte, __LINE__);
 370                 unsigned char a11 = NR_PIXBLOCK_PX(bufy)[((y + 1 - yd0) * bufy->rs) + NR_PIXBLOCK_BPP(bufy) * (x + 1 - xd0) + byte];
 371
 372                 // iterate over the rectangle to be interpolated
 373                 for ( int yi = 0 ; yi < stepy; yi++ ) {
 374                     int iy = stepy - yi;
 375                     int y_out = (y << stepy_l2) + yi;
 376                     if ((y_out < out->area.y0) || (y_out >= out->area.y1))
 377                         continue;
 378                     int out_line  = (y_out - out->area.y0) * out->rs;
 379
 380                     for ( int xi = 0 ; xi < stepx; xi++ ) {
 381                         int ix = stepx - xi;
 382                         int x_out = (x << stepx_l2) + xi;
 383                         if ((x_out < out->area.x0) || (x_out >= out->area.x1))
 384                             continue;
 385
 386                         // simple linear interpolation
 387                         int a = (a00*ix*iy + a10*xi*iy + a01*ix*yi + a11*xi*yi) >> divisor;
 388
 389                         _check_index(out, out_line + NR_PIXBLOCK_BPP(out) * (x_out - out->area.x0) + byte, __LINE__);
 390                         NR_PIXBLOCK_PX(out)[out_line + NR_PIXBLOCK_BPP(out) * (x_out - out->area.x0) + byte] = (unsigned char) a;
 391                     }
 392                 }
 393             }
 394         }
 395     }
 396
 397     nr_pixblock_release(bufy);
 398     delete bufy;
 399
 400     out->empty = FALSE;
 401     slot.set(_output, out);
 402
 403     return 0;
 404 }
 405
 406 int FilterGaussian::get_enlarge(Matrix const &trans)
 407 {
 408     int area_x = _effect_area_scr_x(trans);
 409     int area_y = _effect_area_scr_y(trans);
 410     return _max(area_x, area_y);
 411 }
 412
 413 void FilterGaussian::set_deviation(double deviation)
 414 {
 415     if(isnormal(deviation) && deviation >= 0) {
 416         _deviation_x = _deviation_y = deviation;
 417     }
 418 }
 419
 420 void FilterGaussian::set_deviation(double x, double y)
 421 {
 422     if(isnormal(x) && x >= 0 && isnormal(y) && y >= 0) {
 423         _deviation_x = x;
 424         _deviation_y = y;
 425     }
 426 }
 427
 428 } /* namespace NR */
 429
 430 /*
 431   Local Variables:
 432   mode:c++
 433   c-file-style:"stroustrup"
 434   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 435   indent-tabs-mode:nil
 436   fill-column:99
 437   End:
 438 */
 439 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :