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