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
  16 using std::isnormal;
  17
  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"
  23
  24 namespace NR {
  25
  26 FilterGaussian::FilterGaussian()
  27 {
  28     _deviation_x = _deviation_y = prefs_get_double_attribute("options.filtertest", "value", 0.0);
  29 }
  30
  31 FilterPrimitive *FilterGaussian::create()
  32 {
  33     return new FilterGaussian();
  34 }
  35
  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 }
  42
  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;
  50
  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     }
  56
  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 }
  62
  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 }
  69
  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 }
  76
  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 }
  97
  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 }
 118
 119
 120 int FilterGaussian::render(NRPixBlock **pb, Matrix const &trans)
 121 {
 122     /* in holds the input pixblock */
 123     NRPixBlock *in = pb[0];
 124
 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);
 128
 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;
 137
 138     /* buffer for x-axis blur */
 139     NRPixBlock *bufx = new NRPixBlock;
 140     /* buffer for y-axis blur */
 141     NRPixBlock *bufy = new NRPixBlock;
 142
 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;
 148
 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);
 152
 153     //mid->visible_area = in->visible_area;
 154     //out->visible_area = in->visible_area;
 155
 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)];
 159
 160     /* 1. Blur in direction of X-axis, from in to bufx (they have different resolution)*/
 161     _make_kernel(kernel, _deviation_x, trans.expansionX());
 162
 163     for ( int y = bufx->area.y0 ; y < bufx->area.y1; y++ ) {
 164
 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         }
 174
 175         // current line in bufx
 176         int bufx_line = (y - yd0) * bufx->rs;
 177
 178         int skipbuf[4] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN};
 179
 180         for ( int x = bufx->area.x0 ; x < bufx->area.x1 ; x++ ) {
 181
 182             // for all bytes of the pixel
 183             for ( int byte = 0 ; byte < NR_PIXBLOCK_BPP(in) ; byte++) {
 184
 185                 if(skipbuf[byte] > x) continue;
 186
 187                 double sum = 0;
 188                 int last_in = -1;
 189                 int different_count = 0;
 190
 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 ) {
 193
 194                     // the pixel we're looking at
 195                     int x_in = (((x << stepx_l2) + i + stepx2) >> stepx_l2) << stepx_l2;
 196
 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;
 203
 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                     }
 211
 212                     // value at the pixel
 213                     unsigned char in_byte = NR_PIXBLOCK_PX(in)[in_line + NR_PIXBLOCK_BPP(in) * x_in + byte];
 214
 215                     // is it the same as last one we saw?
 216                     if(in_byte != last_in) different_count++;
 217                     last_in = in_byte;
 218
 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                 }
 222
 223                 // store the result in bufx
 224                 NR_PIXBLOCK_PX(bufx)[bufx_line + NR_PIXBLOCK_BPP(bufx) * (x - xd0) + byte] = (unsigned char)sum;
 225
 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     }
 242
 243
 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());
 246
 247     for ( int x = bufy->area.x0 ; x < bufy->area.x1; x++ ) {
 248
 249         int bufy_disp = NR_PIXBLOCK_BPP(bufy) * (x - xd0);
 250         int bufx_disp = NR_PIXBLOCK_BPP(bufx) * (x - xd0);
 251
 252         int skipbuf[4] = {INT_MIN, INT_MIN, INT_MIN, INT_MIN};
 253
 254         for ( int y = bufy->area.y0; y < bufy->area.y1; y++ ) {
 255
 256             int bufy_line = (y - yd0) * bufy->rs;
 257
 258             for ( int byte = 0 ; byte < NR_PIXBLOCK_BPP(bufx) ; byte++) {
 259
 260                 if (skipbuf[byte] > y) continue;
 261
 262                 double sum = 0;
 263                 int last_in = -1;
 264                 int different_count = 0;
 265
 266                 for ( int i = -scr_len_y ; i <= scr_len_y ; i += stepy ) {
 267
 268                     int y_in = ((((y << stepy_l2) + i + stepy2) >> stepy_l2) - yd0);
 269
 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;
 275
 276                     if (y_in > (yd1 - yd0)) y_in = (yd1 - yd0);
 277                     if (y_in < 0) y_in = 0;
 278
 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                 }
 284
 285                 NR_PIXBLOCK_PX(bufy)[bufy_line + bufy_disp + byte] = (unsigned char)sum;
 286
 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                 }
 297
 298             }
 299         }
 300     }
 301
 302     // we don't need bufx anymore
 303     nr_pixblock_release(bufx);
 304     delete bufx;
 305
 306     // interpolation will need to divide by stepx * stepy
 307     int divisor = stepx_l2 + stepy_l2;
 308
 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);
 313
 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++) {
 317
 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];
 327
 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;
 335
 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;
 341
 342                         // simple linear interpolation
 343                         int a = (a00*ix*iy + a10*xi*iy + a01*ix*yi + a11*xi*yi) >> divisor;
 344
 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     }
 351
 352     nr_pixblock_release(bufy);
 353     delete bufy;
 354
 355     out->empty = FALSE;
 356     pb[1] = out;
 357
 358     return 0;
 359 }
 360
 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 }
 367
 368 void FilterGaussian::set_deviation(double deviation)
 369 {
 370     if(isnormal(deviation) && deviation >= 0) {
 371         _deviation_x = _deviation_y = deviation;
 372     }
 373 }
 374
 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 }
 382
 383 } /* namespace NR */
 384
 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 :