Fast Gaussian blur IIR filter. Building on - “Recursive implementation of the Gaussian filter” by Ian T. Young, Lucas J. van Vliet - “Boundary Conditions for Young - van Vliet Recursive Filtering” by Bill Triggs, Michael Sdika - “Boundary Treatment for Young–van Vliet Recursive Zero-Mean Gabor Filtering” by Vladimir Ulman

fast_gaussian_blur.cpp

#include <algorithm>
#include <array>
#include <memory>
#include <stdexcept>
#include <vector>

typedef std::array<double, 5> b_coefficient_type;
b_coefficient_type b_coefficients(double const sigma)
{
    /* See the Young-van Vliet paper:
    * I.T. Young, L.J. van Vliet, M. van Ginkel, Recursive Gabor filtering.
    * IEEE Trans. Sig. Proc., vol. 50, pp. 2799-2805, 2002.
    *
    * (this is an improvement over Young-Van Vliet, Sig. Proc. 44, 1995)
    *
    * Calculation of q goes according to: http://staff.science.uva.nl/~mark/downloads/anigauss.c
    */

    if (sigma < 0.5)
    {
        throw invalid_argument("Fast Gaussian cannot handle sigma < 0.5.");
    }

    double const q = sigma < 3.556
        ? -0.2568 + 0.5784 * sigma + 0.0561 * SQR(sigma)
        : 2.5091 + 0.9804 * (sigma - 3.556);
    double const q2 = SQR(q);

    double const m0 = 1.16680;
    double const m1 = 1.10783;
    double const m2 = 1.40586;
    double const m1sq = SQR(m1);
    double const m2sq = SQR(m2);

    b_coefficient_type b;
    // [0] is scale...
    b[0] = 1.0 / ((m0 + q) * (m1sq + m2sq + 2.0 * m1 * q + q2));

    b[1] = q * (2.0 * m0 * m1 + m1sq + m2sq + (2.0 * m0 + 4.0 * m1) * q + 3.0 * q2)* b[0];
    b[2] = -q2 * (m0 + 2.0 * m1 + 3.0 * q)* b[0];
    b[3] = q2 * q* b[0];

    // ... and [4] is B.
    b[4] = m0 * (m1sq + m2sq)* b[0];
    b[4] *= b[4];

    b[0] = 1.0 / (1.0 - b[1] - b[2] - b[3]);

    return b;
}

typedef std::array<double, 9> M_matrix_type;
M_matrix_type M_matrix(b_coefficient_type const b)
{
    double const scale = b[0] / ((1.0 + b[1] - b[2] + b[3]) * (1.0 + b[2] + (b[1] - b[3]) * b[3]));

    double const b2sq = SQR(b[2]);
    double const b3sq = SQR(b[3]);
    double const b31 = b[3] * b[1];
    double const b32 = b[3] * b[2];
    double const sb3 = scale * b[3];
    // Helper matrix (3x3) for the forward-backward transition.
    M_matrix_type M;
    M[0] = scale * (-b31 + 1.0 - b3sq - b[2]);
    M[1] = scale * (b[3] + b[1]) * (b[2] + b31);
    M[3] = scale * (b[1] + b32);
    M[2] = M[3] * b[3];
    M[4] = -scale * (b[2] - 1.0) * (b[2] + b31);
    M[5] = -sb3 * (b31 + b3sq + b[2] - 1.0);
    M[6] = scale * (b31 + b[2] + b[1] * b[1] - b2sq);
    M[7] = scale * (b[1] * b[2] + b[3] * b2sq - b[1] * b3sq - b3sq * b[3] - b32 + b[3]);
    M[8] = sb3 * (b[1] + b32);

    return M;
}

#define FAST_GAUSS_PADDING 3
template<typename Iter>
void forward_backward_iir(Iter const begin, Iter const end, b_coefficient_type const & b, M_matrix_type const & M)
{
    //////////////0/////////1/////////2/////////3
    //////////////0123456789012345678901234567890123456
    char msg[] = "The IIR Gauss filter needs at least D pixels to work.";
    msg[36] = '0' + FAST_GAUSS_PADDING + 1;

    if (begin == end)
    {
        throw std::length_error(msg);
    }

    double buffer[FAST_GAUSS_PADDING];

    // Establish left "safe access" padding.
    double const u0 = *begin * b[0];

    Iter cursor = begin;
    double last = *cursor;
    buffer[0] = *cursor += b[1] * u0 + b[2] * u0 + b[3] * u0;
    if (++cursor == end)
    {
        throw std::length_error(msg);
    }
    last = *cursor;
    buffer[1] = *cursor += b[1] * buffer[0] + b[2] * u0 + b[3] * u0;
    if (++cursor == end)
    {
        throw std::length_error(msg);
    }
    last = *cursor;
    buffer[2] = *cursor += b[1] * buffer[1] + b[2] * buffer[0] + b[3] * u0;
    if (++cursor == end)
    {
        throw std::length_error(msg);
    }

    // Forward pass.
    int ring = 3;
    for (; cursor != end; ++cursor)
    {
        last = *cursor;
        *cursor += b[1] * buffer[(ring - 1) % FAST_GAUSS_PADDING] + b[2] * buffer[(ring - 2) % FAST_GAUSS_PADDING] + b[3] * buffer[ring /* - 3 */ % FAST_GAUSS_PADDING];
        buffer[ring++ % FAST_GAUSS_PADDING] = *cursor;
    }

    // Prepare right "safe access" padding.
    double const uplus = last * b[0];
    double const vplus = uplus * b[0];

    double rightBuffer[3];
    for (int i = FAST_GAUSS_PADDING - 1; i >= 0; i--)
    {
        rightBuffer[i] = vplus +
            M[i * 3 + 0] * (buffer[(ring - 1) % FAST_GAUSS_PADDING] - uplus) +
            M[i * 3 + 1] * (buffer[(ring - 2) % FAST_GAUSS_PADDING] - uplus) +
            M[i * 3 + 2] * (buffer[(ring - 3) % FAST_GAUSS_PADDING] - uplus);
    }

    // Establish right "safe access" padding.
    --cursor;
    --ring;
    buffer[ring % FAST_GAUSS_PADDING] = *cursor = rightBuffer[0];
    --cursor;
    --ring;
    buffer[ring % FAST_GAUSS_PADDING] = *cursor += b[1] * buffer[(ring + 1) % FAST_GAUSS_PADDING] + b[2] * rightBuffer[1] + b[3] * rightBuffer[2];
    --cursor;
    --ring;
    buffer[ring % FAST_GAUSS_PADDING] = *cursor += b[1] * buffer[(ring + 1) % FAST_GAUSS_PADDING] + b[2] *  buffer[(ring + 2) % FAST_GAUSS_PADDING] + b[3] * rightBuffer[1];

    // Backward pass.
    do
    {
        --cursor;
        --ring;
        buffer[ring % FAST_GAUSS_PADDING] = *cursor += b[1] * buffer[(ring + 1) % FAST_GAUSS_PADDING] + b[2] * buffer[(ring + 2) % FAST_GAUSS_PADDING] + b[3] * buffer[(ring + 3) % FAST_GAUSS_PADDING];
    } while (cursor != begin);

    // Normalize
    std::transform(begin, end, begin, [&b] (double x)
    {
        return b[4] * x;
    });
}

void blur(std::vector<double> & numbers, double sigma)
{
    b_coefficient_type const b = b_coefficients(sigma);
    forward_backward_iir(numbers.begin(), numbers.end(), b, M_matrix(b));
}