Estimating the standard deviation of Gaussian filter from smoothed image

Question

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')

Answer 1

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');
subplot(3,2,1);
imagesc(img);colormap(gray); axis off; title('Noise image')
subplot(3,2,2);
imagesc(Ig);colormap(gray); axis off; title('filtered image')

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

subplot(3,2,3);
imagesc(abs(Img));
subplot(3,2,4)
imagesc(abs(GImg));

subplot(3,2,5);
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;


subplot(3,2,6);
div = abs(GImg ./ Img);
plot(fftshift(abs(div(1,:))));
title('Division of known input with output');

Answer 2

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.