nr-filter-gaussian.cpp revision 6a5721c2a53716efbebc3f6c2666a5d6ffbdeeb2
#define __NR_FILTER_GAUSSIAN_CPP__
/*
* Gaussian blur renderer
*
* Authors:
* Niko Kiirala <niko@kiirala.com>
* bulia byak
* Jasper van de Gronde <th.v.d.gronde@hccnet.nl>
*
* Copyright (C) 2006-2008 authors
*
* Released under GNU GPL, read the file 'COPYING' for more information
*/
#include "config.h" // Needed for HAVE_OPENMP
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdlib>
#include <glib.h>
#include <limits>
#if HAVE_OPENMP
#include <omp.h>
#endif //HAVE_OPENMP
#include "display/cairo-utils.h"
#include "display/nr-filter-primitive.h"
#include "display/nr-filter-gaussian.h"
#include "display/nr-filter-types.h"
#include "display/nr-filter-units.h"
#include "display/nr-filter-slot.h"
#include "util/fixed_point.h"
#include "preferences.h"
#ifndef INK_UNUSED
#define INK_UNUSED(x) ((void)(x))
#endif
// IIR filtering method based on:
// L.J. van Vliet, I.T. Young, and P.W. Verbeek, Recursive Gaussian Derivative Filters,
// in: A.K. Jain, S. Venkatesh, B.C. Lovell (eds.),
// ICPR'98, Proc. 14th Int. Conference on Pattern Recognition (Brisbane, Aug. 16-20),
// IEEE Computer Society Press, Los Alamitos, 1998, 509-514.
//
// Using the backwards-pass initialization procedure from:
// Boundary Conditions for Young - van Vliet Recursive Filtering
// Bill Triggs, Michael Sdika
// IEEE Transactions on Signal Processing, Volume 54, Number 5 - may 2006
// Number of IIR filter coefficients used. Currently only 3 is supported.
// "Recursive Gaussian Derivative Filters" says this is enough though (and
// some testing indeed shows that the quality doesn't improve much if larger
// filters are used).
static size_t const N = 3;
}
// Type used for IIR filter coefficients (can be 10.21 signed fixed point, see Anisotropic Gaussian Filtering Using Fixed Point Arithmetic, Christoph H. Lampert & Oliver Wirjadi, 2006)
typedef double IIRValue;
// Type used for FIR filter coefficients (can be 16.16 unsigned fixed point, should have 8 or more bits in the fractional part, the integer part should be capable of storing approximately 20*255)
template<typename T> static inline T sqr(T const& v) { return v*v; }
template<typename T> static inline T clip(T const& v, T const& a, T const& b) {
if ( v < a ) return a;
if ( v > b ) return b;
return v;
}
}
/*
template<>
inline unsigned char round_cast(double v) {
// This (fast) rounding method is based on:
#if G_BYTE_ORDER==G_LITTLE_ENDIAN
double const dmr = 6755399441055744.0;
v = v + dmr;
return ((unsigned char*)&v)[0];
#elif G_BYTE_ORDER==G_BIG_ENDIAN
double const dmr = 6755399441055744.0;
v = v + dmr;
return ((unsigned char*)&v)[7];
#else
static double const rndoffset(.5);
return static_cast<unsigned char>(v+rndoffset);
#endif
}*/
if ( v < minval ) return minval_rounded;
if ( v > maxval ) return maxval_rounded;
return round_cast<Tt>(v);
}
if ( v < minval ) return minval_rounded;
if ( v > maxval ) return maxval_rounded;
return round_cast<Tt>(v);
}
namespace Inkscape {
namespace Filters {
{
}
{
return new FilterGaussian();
}
{
// Nothing to do here
}
static int
_effect_area_scr(double const deviation)
{
}
static void
{
double k[scr_len+1]; // This is only called for small kernel sizes (above approximately 10 coefficients the IIR filter is used)
// Compute kernel and sum of coefficients
// Note that actually only half the kernel is computed, as it is symmetric
double sum = 0;
for ( int i = scr_len; i >= 0 ; i-- ) {
if ( i > 0 ) sum += k[i];
}
// the sum of the complete kernel is twice as large (plus the center element which we skipped above to prevent counting it twice)
// Normalize kernel (making sure the sum is exactly 1)
double ksum = 0;
for ( int i = scr_len; i >= 1 ; i-- ) {
}
}
// Return value (v) should satisfy:
// 2^(2*v)*255<2^32
// 255<2^(32-2*v)
// 2^8<=2^(32-2*v)
// 8<=32-2*v
// 2*v<=24
// v<=12
static int
{
// To make sure FIR will always be used (unless the kernel is VERY big):
// deviation/3 <= step
// log(deviation/3) <= log(step)
// So when x below is >= 1/3 FIR will almost always be used.
// This means IIR is almost only used with the modes BETTER or BEST.
int stepsize_l2;
switch (quality) {
case BLUR_QUALITY_WORST:
// 2 == log(x*8/3))
// 2^2 == x*2^3/3
// x == 3/2
break;
case BLUR_QUALITY_WORSE:
// 2 == log(x*16/3))
// 2^2 == x*2^4/3
// x == 3/2^2
break;
case BLUR_QUALITY_BETTER:
// 2 == log(x*32/3))
// 2 == x*2^5/3
// x == 3/2^4
break;
case BLUR_QUALITY_BEST:
stepsize_l2 = 0; // no subsampling at all
break;
case BLUR_QUALITY_NORMAL:
default:
// 2 == log(x*16/3))
// 2 == x*2^4/3
// x == 3/2^3
break;
}
return stepsize_l2;
}
static void calcFilter(double const sigma, double b[N]) {
assert(N==3);
double const d3_org = 1.85132;
double qbeg = 1; // Don't go lower than sigma==2 (we'd probably want a normal convolution in that case anyway)
double s;
do { // Binary search for right q (a linear interpolation scheme is suggested, but this should work fine as well)
// Compute scaled filter coefficients
// Compute actual sigma^2
qbeg = q;
} else {
qend = q;
}
// Compute filter coefficients
b[2] = -bscale;
}
static void calcTriggsSdikaM(double const b[N], double M[N*N]) {
assert(N==3);
for(unsigned int i=0; i<9; i++) M[i] *= Mscale;
}
template<unsigned int SIZE>
static void calcTriggsSdikaInitialization(double const M[N*N], IIRValue const uold[N][SIZE], IIRValue const uplus[SIZE], IIRValue const vplus[SIZE], IIRValue const alpha, IIRValue vold[N][SIZE]) {
for(unsigned int c=0; c<SIZE; c++) {
double uminp[N];
for(unsigned int i=0; i<N; i++) {
double voldf = 0;
for(unsigned int j=0; j<N; j++) {
}
// Properly takes care of the scaling coefficient alpha and vplus (which is already appropriately scaled)
// This was arrived at by starting from a version of the blur filter that ignored the scaling coefficient
// (and scaled the final output by alpha^2) and then gradually reintroducing the scaling coefficient.
}
}
}
// Filters over 1st dimension
static void
{
#if G_BYTE_ORDER == G_LITTLE_ENDIAN
#else
static unsigned int const alpha_PC = 0;
#endif
#if HAVE_OPENMP
#else
#endif // HAVE_OPENMP
#if HAVE_OPENMP
unsigned int tid = omp_get_thread_num();
#else
unsigned int tid = 0;
#endif // HAVE_OPENMP
// corresponding line in the source and output buffer
// Border constants
// Forward pass
for(unsigned int c=0; c<PC; c++) u[0][c] *= b[0];
for(unsigned int i=1; i<N+1; i++) {
for(unsigned int c=0; c<PC; c++) u[0][c] += u[i][c]*b[i];
}
}
// Backward pass
if ( PREMULTIPLIED_ALPHA ) {
PREMUL_ALPHA_LOOP dstimg[c] = clip_round_cast_varmax<PT>(v[0][c], v[0][alpha_PC], dstimg[alpha_PC]);
} else {
}
while(c1-->0) {
for(unsigned int c=0; c<PC; c++) v[0][c] *= b[0];
for(unsigned int i=1; i<N+1; i++) {
for(unsigned int c=0; c<PC; c++) v[0][c] += v[i][c]*b[i];
}
if ( PREMULTIPLIED_ALPHA ) {
PREMUL_ALPHA_LOOP dstimg[c] = clip_round_cast_varmax<PT>(v[0][c], v[0][alpha_PC], dstimg[alpha_PC]);
} else {
}
}
}
}
// Filters over 1st dimension
// Assumes kernel is symmetric
// Kernel should have scr_len+1 elements
static void
{
// Past pixels seen (to enable in-place operation)
#if HAVE_OPENMP
#else
#endif // HAVE_OPENMP
// corresponding line in the source buffer
// current line in the output buffer
// history initialization
// update history
// for all bytes of the pixel
int last_in = -1;
int different_count = 0;
// go over our point's neighbours in the history
for ( int i = 0 ; i <= scr_len ; i++ ) {
// value at the pixel
// is it the same as last one we saw?
// sum pixels weighted by the kernel
}
// go over our point's neighborhood on x axis in the in buffer
for ( int i = 1 ; i <= scr_len ; i++ ) {
// the pixel we're looking at
} else {
nb_src_disp += sstr1;
}
// value at the pixel
// is it the same as last one we saw?
// sum pixels weighted by the kernel
}
// store the result in bufx
// optimization: if there was no variation within this point's neighborhood,
// skip ahead while we keep seeing the same last_in byte:
// blurring flat color would not change it anyway
if (different_count <= 1) { // note that different_count is at least 1, because last_in is initialized to -1
pos++;
nb_src_disp += sstr1;
nb_dst_disp += dstr1;
}
}
}
}
}
}
static void
{
// Filter variables
double bf[N]; // computed filter coefficients
double M[N*N]; // matrix used for initialization procedure (has to be double)
// Compute filter
b[0] = 1; // b[0] == alpha (scaling coefficient)
for(size_t i=0; i<N; i++) {
b[i+1] = bf[i];
b[0] -= b[i+1];
}
// Compute initialization matrix
calcTriggsSdikaM(bf, M);
int w = cairo_image_surface_get_width(src);
int h = cairo_image_surface_get_height(src);
// Filter
switch (cairo_image_surface_get_format(src)) {
case CAIRO_FORMAT_A8: ///< Grayscale
filter2D_IIR<unsigned char,1,false>(
w, h, b, M, tmpdata, num_threads);
break;
case CAIRO_FORMAT_ARGB32: ///< Premultiplied 8 bit RGBA
filter2D_IIR<unsigned char,4,true>(
w, h, b, M, tmpdata, num_threads);
break;
default:
assert(false);
};
}
static void
int num_threads)
{
// Filter kernel for x direction
int w = cairo_image_surface_get_width(src);
int h = cairo_image_surface_get_height(src);
// Filter (x)
switch (cairo_image_surface_get_format(src)) {
case CAIRO_FORMAT_A8: ///< Grayscale
filter2D_FIR<unsigned char,1>(
break;
case CAIRO_FORMAT_ARGB32: ///< Premultiplied 8 bit RGBA
filter2D_FIR<unsigned char,4>(
break;
default:
assert(false);
};
}
{
if (!in) return;
// zero deviation = transparent black as output
if (_deviation_x <= 0 || _deviation_y <= 0) {
return;
}
int bytes_per_pixel = 0;
switch (fmt) {
case CAIRO_FORMAT_A8:
bytes_per_pixel = 1; break;
case CAIRO_FORMAT_ARGB32:
default:
bytes_per_pixel = 4; break;
}
#if HAVE_OPENMP
#else
int threads = 1;
#endif
int w_downsampled = resampling ? static_cast<int>(ceil(static_cast<double>(w_orig)/x_step))+1 : w_orig;
int h_downsampled = resampling ? static_cast<int>(ceil(static_cast<double>(h_orig)/y_step))+1 : h_orig;
// Decide which filter to use for X and Y
// This threshold was determined by trial-and-error for one specific machine,
// so there's a good chance that it's not optimal.
// Whatever you do, don't go below 1 (and preferrably not even below 2), as
// the IIR filter gets unstable there.
// Temporary storage for IIR filter
// NOTE: This can be eliminated, but it reduces the precision a bit
for(int i = 0; i < threads; ++i) {
}
}
if (resampling) {
cairo_scale(ct, static_cast<double>(w_downsampled)/w_orig, static_cast<double>(h_downsampled)/h_orig);
} else {
}
if (scr_len_x > 0) {
if (use_IIR_x) {
} else {
}
}
if (scr_len_y > 0) {
if (use_IIR_y) {
} else {
}
}
if (resampling) {
cairo_surface_t *upsampled = cairo_surface_create_similar(downsampled, cairo_surface_get_content(downsampled),
cairo_scale(ct, static_cast<double>(w_orig)/w_downsampled, static_cast<double>(h_orig)/h_downsampled);
} else {
}
}
{
// maximum is used because rotations can mix up these directions
// TODO: calculate a more tight-fitting rendering area
}
{
// Previously we tried to be smart and return true for rotations.
// However, the transform passed here is NOT the total transform
// from filter user space to screen.
// TODO: fix this, or replace can_handle_affine() with isotropic().
return false;
}
{
}
}
void FilterGaussian::set_deviation(double x, double y)
{
_deviation_x = x;
_deviation_y = y;
}
}
} /* namespace Filters */
} /* namespace Inkscape */
/*
Local Variables:
mode:c++
c-file-style:"stroustrup"
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
indent-tabs-mode:nil
fill-column:99
End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :