src/display/nr-filter-gaussian.cpp

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