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 <vector>
#include <stdexcept>
#include <iterator>

using namespace std;

namespace
{
	inline double sqr(double const x)
	{
		return x * x;
	}

	enum
	{
		FAST_GAUSS_PADDING = 3
	};

	typedef 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 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;
	}

	template<typename Iter>
	void forward_backward_iir(Iter const begin, Iter const end, b_coefficient_type const & b, M_matrix_type const & M)
	{
		if (distance(begin, end) <= FAST_GAUSS_PADDING)
		{
			//////////////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;

			throw 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;
		++cursor == end;
		last = *cursor;
		buffer[1] = *cursor += b[1] * buffer[0] + b[2] * u0 + b[3] * u0;
		++cursor == end;
		last = *cursor;
		buffer[2] = *cursor += b[1] * buffer[1] + b[2] * buffer[0] + b[3] * u0;
		++cursor == end;

		// 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[] = { vplus, vplus, vplus };
		for (int i = FAST_GAUSS_PADDING - 1; i >= 0; --i)
		{
			rightBuffer[i] +=
				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
		transform(begin, end, begin, [&b](double const x)
		{
			return b[4] * x;
		});
	}
}

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