title: Variable radius Gaussian blur, approximating the kernel
body: <p>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.</p>
  
  <p>Now, I'm trying to approximate it using the following scheme:</p>
  
  <p>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<sub>k</sub> and coefficients α<sub>k</sub> as:</p>
  
  <p>    f(x,y) = ∑<sub>k=1</sub><sup>N</sup> α<sub>k</sub>·χ<sub>R<sub>k</sub></sub>(x,y)</p>
  
  <p>Let g(x,y) be our image, then </p>
  
  <p>    ∬<sub>ℝ<sup>2</sup></sub> K(x,y)·g(x,y) dxdy ≈ 
  ∬<sub>ℝ<sup>2</sup></sub> f(x,y)·g(x,y) dxdy = ∑<sub>k=1</sub><sup>N</sup> α<sub>k</sub>·∬<sub>R<sub>k</sub></sub>g(x,y) dxdy</p>
  
  <p>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.</p>
  
  <p>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.</p>
  
  <p>The remaining part is to find a good approximation function f(x,y). <strong>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)?</strong></p>