# Gaussian Equivalent of Convolving an Image 50 Times with a Box Filter

I have a image that I'm convolving 50 times with a box filter, and I like to replace this with a single gaussian filter. Because both box filter and gaussian filter are separable I can just study the 1D case.

Therefore...

Let $$l \geq 1$$ be fixed and $$n = 50 \approx +\infty$$ The box filter I'm using is defined as

$$h_l(x) = \frac{1}{2l+1} \begin{cases} 1 & -l\leq x \leq l \\ 0 & \text{otherwise} \end{cases},$$

here $$x \in \mathbb{Z}$$.

After a while I got that the gaussian equivalent for given $$l,n$$ (large $$n$$) is given by

$$f_{n,l}(x) = \frac{1}{\sqrt{2\pi \sigma_{l,n}^2}} e^{-\frac{x^2}{2 \sigma_{l,n}^2}},$$ where $$\sigma_{l,n} = \sqrt{\frac{n(l+1)l}{3}},$$ the proof is in the answer of Fourier Transform of a Sequence and Inverse Fourier Transform.

Therefore I approximated the operator

$$A_{h_l,n}={\underbrace{h_l*h_l*\ldots *h_l}_{\text{n times}} * (\cdot)}$$

with

$$G_{l,n} = f_{l,n}*(\cdot)$$

However for practical purposes I had to consider the $$f_{l,n}$$ defined on a compact support (like $$\left[-k\sigma,k\sigma \right], k \in \left\{1,2,3\right\}$$.

The questions are:

1. Does this approach make sense?
2. If yes, and I've already implemented it for some reason using a single convolution with a gaussian (large kernel though) with 50 times the one with the box filter (you can assume $$l = 5$$ for my experiments) yields worse results, timewise speaking, assuming the approach is correct what can the reason be?

So either my approach is wrong or I'm probably missing something.

Just to clarify I'm trying to replace

for(int i = 0; i < n; ++i)
blur(image,image,Size(2*l+1,2*l+1);


With

GaussianBlur(image,image,Size(2*k*sigma,2*k*sigma),sigma,sigma);


in opencv, and hopefully gain a speed-up.

What you should do is work with the Variance.

The Variance of Discrete Random Variables with support of $$\left[ -l, l \right]$$ is given by (Notice that since the Expected Value is $$0$$, the Variance is given by the 2nd moment)):

$$\sum_{k = -l}^{l} \mathbb{P} \left( k \right) {k}^{2} = \sum_{k = -l}^{l} \frac{1}{2l + 1} {k}^{2} = \frac{1}{2l + 1} \sum_{k = -l}^{l} {k}^{2}$$

Convolution of a function which is the PDF of a random variable is equivalent of addition of the Random Variables (Assuming they are independent).
The Variance is linear with addition of Independent Variables hence the Variance of $$50$$ Convolutions (Which is addition of 51 variables) is given by: $$\frac{50}{2l + 1} \sum_{k = -l}^{l} {k}^{2}$$

Now generate a Gaussian Blur with $$\sigma = \sqrt{ \frac{50}{2l + 1} \sum_{k = -l}^{l} {k}^{2} }$$ apply it on the image and it will be the best approximation for applying 50 times the Box Filter with Radius of $$l$$.

By the Central Limit Theorem we can get better and better approximation as the number of Kernels gets higher.

Remarks

1. Usually Gaussian Blur with large STD is approximated by few (4-5) iterations of Box Blur. Yet in this case since there is a large number of convolutions the proper way is to calculate the STD of equivalent Gaussian Blur and then apply it either directly or by accelerated method (For example, approximation by few box blur iterations).
2. Box Blur can be applied in O(n) manner. Namely complexity is independent of the box kernel radius. It can be done by Integral Image or Running Sum.
3. Pay attention that if we do 50 Convolution of Box Blur with itself it is equivalent of 51 additions of Random Variables. Yet above the first convolution is done between the Box Blur Kernel and the images hence 50 Convolution means addition of 50 Random Variables.

### MATALB Code

boxRadius   = 4;
boxSupport  = (2 * boxRadius) + 1;
vB          = ones([boxSupport, 1]) / boxSupport;

boxStd = sqrt(sum((boxGrid .^ 2) .* vB));

vT = vB;

numConv = 9;

for ii = 1:numConv
vT = conv2(vB, vT, 'full');
end

kernelRadius    = floor(size(vT, 1) / 2);
kernelSupport   = (2 * kernelRadius) + 1;

kernelStd = sqrt(sum((vKernelGrid .^ 2) .* vT));

kernelStd / boxStd %<! Should be 'numConv + 1

vGaussianKernel = (1 / sqrt(2 * pi * kernelStd * kernelStd)) * exp(-(vKernelGrid .^ 2) / (2 * kernelStd * kernelStd));

figure();
plot(vKernelGrid, [vT, vGaussianKernel]);
xlabel('Kernel Support');
ylabel('Kernel Value');
legend({['Convolved Box Blur'], ['Gaussian Blur']});


## Closed Form Formula for the Variance of Centered Uniform Discrete Random Variable

As seen above if $$x$$ is a centered uniform Random Variable on the support $$\left\{ -l, -l + 1, \ldots, 0, 1, \ldots, l - 1, l \right\}$$ then its variance is given by the second moment (As it is centered which means the mean value is $$0$$):

$$\operatorname{Var} \left( x \right) = \sum_{k = -l}^{l} \mathbb{P} \left( x = k \right) {k}^{2} = \frac{2}{2 l + 1} \sum_{k = 1}^{l} {k}^{2}$$

The above sum of squares can be solved either by integral or nice trick form the Binomial Identity:

$${\left( k - 1 \right)}^{3} = {k}^{3} - 3 {k}^{2} + 3 k - 1 \Rightarrow {k}^{3} - {\left( k - 1 \right)}^{3} = 3 {k}^{2} - 3 k + 1$$

Summing each side from $$k = 1$$ until $$n$$:

$$\sum_{k = 1}^{n} \left( {k}^{3} - {\left( k - 1 \right)}^{3} \right) = {n}^{3} = 3 \sum_{k = 1}^{n} {k}^{2} - 3 \sum_{k = 1}^{n} k + \sum_{k = 1}^{n} 1 = \sum_{k = 1}^{n} {k}^{2} - 3 \frac{n \left( n + 1 \right)}{2} - n$$

Rearranging yields:

$$\sum_{k = 1}^{n} = \frac{n \left( n + 1 \right) \left( 2 n + 1 \right)}{6}$$

Hence, for $$x$$ as above:

$$\operatorname{Var} \left( x \right) = \sum_{k = -l}^{l} \mathbb{P} \left( x = k \right) {k}^{2} = \frac{2}{2 l + 1} \sum_{k = 1}^{l} {k}^{2} = \frac{2 l \left( l + 1 \right) \left( 2 l + 1 \right)}{6 \left( 2 l + 1 \right)} = \frac{l \left( l + 1 \right)}{3}$$

• How is your approach different from mine? Work out the math you should end up with my same result. Also, what you computed is $\sigma^2$ not $\sigma$. Jun 1 '18 at 10:30
• I'm not 100% convinced of your $\sigma$ formula yet. For the accelleration bit then, which is more what I'm interested in, trying to speed-up the iterative box blur with a gaussian blur is totally wrong then. How about replacing the 50 times box filter with a single box filter but with a larger kernel, can that speed up things? Jun 1 '18 at 10:46
• In one line you say: The std is $$\frac{50}{2l + 1} \sum_{k = -l}^{l} {k}^{2}$$ In the next you say $$\sigma = 50 \sqrt{ \frac{1}{2l + 1} \sum_{k = -l}^{l} {k}^{2} }$$ Jun 1 '18 at 12:15
• Shouldn't "50" go under the square root? Jun 1 '18 at 12:57
• Just noticed you got rid of the $2500$ under the square root, ok now this matches my original thought. Thank you. Jun 1 '18 at 13:59