Firstly, let's say that in order to smooth an image, I convolve it with a Gaussian function having standard deviation $\sigma_x$ and $\sigma_y$. I am now interested in knowing if there exist methods for estimating $\sigma_x$ and $\sigma_y$ from the smoothed image?

Secondly, to potentially make things more complicated, let's say the image pixels instead are drawn from a zero mean unit variance normal distribution (top image), which is then convolved as above (bottom image), could the standard deviations be estimated from the smoothed image?

Here is the MATLAB code used to generate the images. What I would like to achieve is to estimate the standard deviations, which in this case were $\sigma_x=2$ and $\sigma_y=2$.

img = randn(128,128);
figure(1); imagesc(img); colormap(gray); axis off; title('Noise image')
nimg = imgaussfilt(img, [2  2]);
figure(2); imagesc(nimg); colormap(gray); axis off; title('Smoothed noise image')

enter image description here enter image description here


You can exploit the convolution theorem: Convolving with a Gaussian corresponds to Multiplication with a Gaussian of inverse $\sigma$ in frequency. Let x be the input signal and y be the output signal. Further, Let g be the Gaussian filter.

$$ y = g*x $$

Or, in frequency domain (using captial letters) $$ Y=GX $$

Hence, if both Y and X is known, you can do $$ G=Y/X $$ to get your Gaussian Kernel. However, due to boundary conditions etc this will not yield the perfect Kernel. You have to fit a Gaussian function to it to get its standard deviation.

Additionally, if you have no knowledge of the input signal, but know its white noise, you can also fit a Gaussian Curve on the FFT of the output (since the input is white noise and hence flat):

N = 16*128;
img = randn(N,N);
s = 2;

G = fspecial('gaussian',[25 25],s);
Ig = imfilter(img,G,'replicate');
imagesc(img);colormap(gray); axis off; title('Noise image')
imagesc(Ig);colormap(gray); axis off; title('filtered image')

% Perform FFT
Img = fft2(img);
GImg = fft2(Ig);


fy = -N/2:(N/2-1);
hold off;
plot(fy, cumsum(fftshift(abs(GImg(:,1)))));
hold on;
plot(fy, normcdf(fy, 0, N/(2*pi*s))*sum(abs(GImg(:,1))), 'r');
hold off;

div = abs(GImg ./ Img);
title('Division of known input with output');

program output

  • $\begingroup$ Thank you! Your method worked well. To extract the standard deviation I used lsqcurvefit in Matlab to fit the normcdf to the curve seen in your subplot(326), which also worked for different sigma values in the Gaussian filter. $\endgroup$ – Smajjk Nov 23 '16 at 21:19

Let's consider your input image as x ,output image as y and g as the gaussian kernel. The output image is formed by convolution of gaussian kernel and the input image. In frequency domain it will change to multiplication.

In frequency domain,

Y(u,v) = H(u,v)X(u,v)

H(u,v) = Y(u,v)/X(u,v)

Gaussian in spatial domain is same as in frequency domain.

$H(u,v) = e^{-D^2(u,v)/2Sigma^2}$ from which you can estimate the value of sigma.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.