4

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 Rk and coefficients αk as:

    f(x,y) = ∑k=1N αk·χRk(x,y)

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

    ∬ℝ2 K(x,y)·g(x,y) dxdy ≈ ∬ℝ2 f(x,y)·g(x,y) dxdy = ∑k=1N α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 (minimizing the error) or given the error (minimizing the number of rectangles)?

Yakov Galka
  • 70,775
  • 16
  • 139
  • 220
  • How does the variable radius enter into the approximation to the gaussian kernel? Also, your equation for the convolution (blur) seems odd: do you mean to write `G(x0, y0) = ∬ K(x-x0,y-y0)·g(x,y) dxdy`? In any case, another approach you can take is to perform several convolutions with different kernels and then linearly interpolate between them in a position dependent way. – Kipton Barros Sep 25 '11 at 01:06
  • 1
    @Kipton: I assume that the kernel is centered about (0,0) and has a fixed radius. Translating and scaling the set of the rectangles that approximate it can be done in constant time per pixel without changing the relative error rate. – Yakov Galka Sep 25 '11 at 07:07
  • @Kipton: Your approach has linearithmic time and space complexity, which becomes quadratic when tending the error to zero. – Yakov Galka Sep 25 '11 at 07:13

1 Answers1

0

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.

mcdowella
  • 19,301
  • 2
  • 19
  • 25
  • It's quite obvious that the configuration is symmetric with respect to the origin, and without loss of generality the rectangles are centered (if not, there are 4 symmetric rectangles that can be replaced by a combination of 4 centered rectangles: one big, two in a `+` figure and a small one). However I don't understand why you say it's a 1D problem. Who said that all the rectangles have the same aspect ratio? Why can't they be overlapping but not nested? In fact I'm almost sure the optimal solution will contain non-nested rectangles. – Yakov Galka Sep 25 '11 at 10:41
  • What I meant by a 1D problem is that - under the constraints I assumed - there is only one degree of freedom in a rectangle: its size. In turn, this means that a dynamic programming approach can work its way along a 1-dimensional line - increasing N in my answer - which is generally the difference between a dynamic programming problem being tractable and not. If the global optimum solution cannot be found by considering nested rectangles - or at least a set of non-intersecting figures of predetermined shape - then this approach is probably not practical - sorry. – mcdowella Sep 25 '11 at 11:25