Variable radius Gaussian blur, approximating the kernel

I'm writing a Gaussian blur with variable radius (standard deviation), i.e. each pixel of the image is convolved using a different kernel. The standard techniques to compute Gaussian blur don't work here: FFT, axis-separation, repeated box-blur—they all assume that the kernel is the same for the whole image.

Now, I'm trying to approximate it using the following scheme:

Approximate the Gaussian kernel K(x,y) with a piecewise constant function f(x,y) defined by a set N of axis-aligned rectangles R_k and coefficients α_k as:

    f(x,y) = ∑_k=1^N α_k·χ_Rk(x,y)

Let g(x,y) be our image, then

    ∬_ℝ² K(x,y)·g(x,y) dxdy ≈ ∬_ℝ² f(x,y)·g(x,y) dxdy = ∑_k=1^N α_k·∬_Rkg(x,y) dxdy

The integral on the RHS is a simple integral over a rectangle, and as such can be computed in constant time by precomputing the partial sums for the whole image.

The resulting algorithm runs in O(W·H·N) where W and H are the dimensions of the image and N is (AFAIK) inverse proportional to the error of the the approximation.

The remaining part is to find a good approximation function f(x,y). How to find the optimal approximation to the Gaussian when given either the number of rectangles N (m开发者_Python百科inimizing the error) or given the error (minimizing the number of rectangles)?

Given the locations and size of the rectangles, it should be fairly easy to work out your coefficients, so the real problem is working out where to put the rectangles.

Since you are approximating a Gaussian, it seems at least reasonable to restrict our attention to rectangles whose centre coincides with the centre of the Gaussian, so we actually have only a 1-dimensional problem - working out the sizes of a nested set of rectangles which I presume are either squares or are similar to the Gaussian if you have an aspect ratio other than unity.

This can be solved via dynamic programming. Suppose that you work from the outside into the middle. At stage N you have worked out an n x k table that gives you the best possible approximation error coming from 1,2...N rings of outer pixels for up 1,2,..k different rectangles, and the size of the innermost rectangle responsible for that best error. To work out stage N+1 you consider every possible size for what will be an innermost rectangle so far, contributing x rings of pixels to the outer area. You work out the alpha for that rectangle that gives the best fit for the pixels in the new ring and the rings outside it not left to outer rectangles. Using the values in the table already calculated you know the best possible error you will get when you leave up to k outer rectangles to cover those areas, so you can work out the best total error contributed from what is now N+1 rings of pixels. This allows you to fill in the table entries for N+1 outer pixels. When you have worked your way into the middle of the area you will be able to work out the optimal solution for the whole area.