Approximating Gaussian Blur Using Extended Box Blur

Question

The problem is a s following, how to approximate a Gaussian Blur Filter with a given STD using Box Blur / Extended Box Blur.

More specifically, I know this is the way Photoshop applies its Gaussian Blur.

First, an article about "Extended Box Blur can be seen here - Theoretical Foundations of Gaussian Convolution by Extended Box Filtering.

The problem I'm having is with Figure 2 in the article.
The best way to explain this would be using an example.

Let's say we need to approximate a Gaussian Blur with STD of 15.4 -> Var = 237.16.
In order to have a good approximation we'll do that with 6 iterations of a Box Blur.

Now, How do I choose the length of the Box Blur (We'll do it in a separable manner, namely, working in 1D)?
Should I chose different lengths (It seems I have to)?
The target is matching the GB Level of Blur (Which is its STD / VAR).

Thank You.

P.S.
I'm working on MATLAB, so code is easy :-).

Answer 1

This is my MATLAB implementation of the article:

```

function [ vBoxBlurKernel ] = GenerateBoxBlurKernel( boxBlurVar, numIterations )
% ----------------------------------------------------------------------------------------------- %
% [ boxBlurKernel ] = GenerateBoxBlurKernel( boxBlurVar, numIterations )
%   Approximates 1D Gaussian Kernel by iterative convolutions of "Extended Box Filter".
% Input:
%   - boxBlurVar        -   BoxFilter Varaiance.
%                           The variance of the output Box Filter.
%                           Scalar, Floating Point (0, inf).
%   - numIterations     -   Number of Iterations.
%                           The number of convolution iterations in order
%                           to produce the output Box Filter.
%                           Scalar, Floating Point [1, inf), Integer.
% Output:
%   - vBoxBlurKernel    -   Output Box Filter.
%                           The Box Filter with 'boxBlurVar' Variance.
%                           Vector, Floating Point, (0, 1).
% Remarks:
%   1.  The output Box Filter has a variance of '' as if it is treated as
%       Discrete Probability Function.
%   2.  References: "Theoretical Foundations of Gaussian Convolution by Extended Box Filtering"
%   3.  Prefixes:
%       -   'm' - Matrix.
%       -   'v' - Vector.
% TODO:
%   1.  F
%   Release Notes:
%   -   1.0.001     07/05/2014  xxxx xxxxxx
%       *   Accurate calculation of the "Extended Box Filter" length as in
%           the reference.
%   -   1.0.000     06/05/2014  xxxx xxxxxx
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

boxBlurLength = sqrt(((12 * boxBlurVar) / numIterations) + 1);
boxBlurRadius = (boxBlurLength - 1) / 2;

% 'boxBlurRadiusInt' -> 'l' in the reference
boxBlurRadiusInt    = floor(boxBlurRadius);
% boxBlurRadiusFrac   = boxBlurRadius - boxBlurRadiusInt;

% The length of the "Integer" part of the filter.
% 'boxBlurLengthInt' -> 'L' in the reference
boxBlurLengthInt = 2 * boxBlurRadiusInt + 1;

a1 = ((2 * boxBlurRadiusInt) + 1);
a2 = (boxBlurRadiusInt * (boxBlurRadiusInt + 1)) - ((3 * boxBlurVar) / numIterations);
a3 = (6 * ((boxBlurVar / numIterations) - ((boxBlurRadiusInt + 1) ^ 2)));

alpha = a1 * (a2 / a3);
ww = alpha / ((2 * boxBlurRadiusInt) + 1 + (2 * alpha));

% The length of the "Extended Box Filter".
% 'boxBlurLength' -> '\Gamma' in the reference.
boxBlurLength = (2 * (alpha + boxBlurRadiusInt)) + 1;

% The "Single Box Filter" with Varaince - boxBlurVar / numIterations
% It is normalized by definition.
vSingleBoxBlurKernel = [ww, (ones(1, boxBlurLengthInt) / boxBlurLength), ww];
% vBoxBlurKernel = vBoxBlurKernel / sum(vBoxBlurKernel);

vBoxBlurKernel = vSingleBoxBlurKernel;

% singleBoxKernelVar = sum(([-(boxBlurRadiusInt + 1):(boxBlurRadiusInt + 1)] .^ 2) .* boxBlurKernel)
% boxKernelVar = numIterations * singleBoxKernelVar


for iIter = 2:numIterations
    vBoxBlurKernel = conv2(vBoxBlurKernel, vSingleBoxBlurKernel, 'full');
end


end

Here's a demo to try it:

% Box Blur Demo

gaussianKernelStd = 9.6;
gaussianKernelVar = gaussianKernelStd * gaussianKernelStd;
gaussianKernelRadius = ceil(6 * gaussianKernelStd);
gaussianKernel = exp(-([-gaussianKernelRadius:gaussianKernelRadius] .^ 2) / (2 * gaussianKernelVar));
gaussianKernel = gaussianKernel / sum(gaussianKernel);

boxBlurKernel = GenerateBoxBlurKernel(gaussianKernelVar, 6);
boxBlurKernelRadius = (length(boxBlurKernel) - 1) / 2;

figure();
plot([-gaussianKernelRadius:gaussianKernelRadius], gaussianKernel, [-boxBlurKernelRadius:boxBlurKernelRadius], boxBlurKernel);

sum(([-boxBlurKernelRadius:boxBlurKernelRadius] .^ 2) .* boxBlurKernel)
sum(([-gaussianKernelRadius:gaussianKernelRadius] .^ 2) .* gaussianKernel)

The tricky part is the calculation of the effective length of the "Extended Box Filter".
It's not the length by the calculation of the length using the variance of a regular "Box Filter".

The article is great and this method is amazing.