I'd like to make a faster version of the Gaussian blur for Inkscape. Here's what I found various open source projects use:

  1. librsvg: 3 box blurs when stdev > 2. Else, separable FIR filter.
  2. Krita:  FFT convolution or separable FIR. KisGaussianKernel::applyGaussian
  3. Gimp (GEGL): IIR (aka. recursive filter)
  4. Inkscape: FIR or IIR, depending on size
  5. OpenCV, CPU: uses IPP, separable FIR filter? GaussianBlurFixedPoint

The problem with the current IIR implementation (van Vliet et al. (1998)) is

  1. The forward pass output needs to be a higher precision than the input or else the output won't be smooth - will get mach bands. The output of the backward pass is the same precision as the input, but that could cause accuracy problems later when filtering in the 2nd dimension.

  2. You need a high precision accumulator or else the output won't be smooth. This happens when sigma is very large. I found for sigma > 30, float32 is no longer enough and had to switch to float64 (int32 should work, but harder to do), or else you'll see mach bands. When sigma = 30, the filter is

    y[t] = 0.0001395 x[t] + 2.880 y[t-1] -2.767 y[t-2] + 0.8866 y[t-3]

0.0001395 x[t] is much smaller than 2.880 y[t-1], meaning small changes in x[t] are lost. Effectively, it's as if the input image had been quantized. Changing the order of the accumulation by adding the x[t] term last helps some.

Here're some other formulations I read about,

  1. Deriche (1992). I'm guessing adding the output of 2 filters (forward & backward) is less sensitive to rounding errors than cascading them (van Vliet et al. (1998)). But it says the accuracy is good only for sigma < 32.

  2. D. Demigny, L. Kessal & J. Pons (2002). This one claims, no backward pass is needed to keep the phase constant! 1 minor disadvantage is that sigma can only change in steps of 0.32:

sigma = 0.3217 w + 0.481, % where w is an integer, 
                          % approximately 1/2 the 
                          % width of the gaussian

So, I'm asking is there a good method that can be accurate for large sigma (< 200 pixels) and avoids needing to store high precision intermediate results, which increases cache foot print and reduces SIMD width.

  • $\begingroup$ (A) I doubt there's one best answer, because with image processing not just the instruction set, but also the memory organization matters. Also, the best use of a modern CPU in this case may be to manage the GPU that does the actual computation. (B) your code snippet starting with y[t] = is committing a cardinal sin -- IIR filters should always be split into cascaded second-order sections. If you do it with this filter, you should find that you can use a higher sigma (although still not infinite) before you run into numerical issues. $\endgroup$
    – TimWescott
    Commented Dec 29, 2023 at 18:34
  • $\begingroup$ Great suggestion. I didn't know about factoring the transfer function into lower order polynomials. For the example above, the x coefficient can now be sqrt(0.0001395) = 0.0118, which will require much less precision. Even in the worst case of sigma = 200, where the filter is 4.960e-07 x[t] +2.982 y[t-1] -2.964 y[t-2] + 0.9820 y[t-3], the square root will make float32 enough. But I'll need to understand that Triggs, Sdika method for initializing the filter state and adjust it, now that there's 2 recurrences. $\endgroup$
    – Yale Zhang
    Commented Jan 1 at 11:46


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.