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
- 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).
- 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.
- 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;
boxGrid = [-boxRadius:boxRadius].';
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);
vKernelGrid = [-kernelRadius:kernelRadius].';
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} $$
