Obtain sigma of gaussian blur between two images

Suppose I have an image A, I applied Gaussian Blur on it with Sigam=3 So I got another Image B. Is there a way to know the applied sigma if A,B is given?

Further clarification:

Image A:

Image B:

I want to write a function that take A,B and return Sigma:

double get_sigma(cv::Mat const& A,cv::Mat const& B);

Any suggestions?

Solution

EDIT1: The suggested approach doesn't work in practice in its original form(i.e. using only 9 equations for a 3 x 3 kernel), and I realized this later. See EDIT1 below for an explanation and EDIT2 for a method that works.

EDIT2: As suggested by Humam, I used the Least Squares Estimate (LSE) to find the coefficients.

I think you can estimate the filter kernel by solving a linear system of equations in this case. A linear filter weighs the pixels in a window by its coefficients, then take their sum and assign this value to the center pixel of the window in the result image. So, for a 3 x 3 filter like

the resulting pixel value in the filtered image

result_pix_value = h11 * a(y, x) + h12 * a(y, x+1) + h13 * a(y, x+2) +
                   h21 * a(y+1, x) + h22 * a(y+1, x+1) + h23 * a(y+1, x+2) +
                   h31 * a(y+2, x) + h32 * a(y+2, x+1) + h33 * a(y+2, x+2)

where a's are the pixel values within the window in the original image. Here, for the 3 x 3 filter you have 9 unknowns, so you need 9 equations. You can obtain those 9 equations using 9 pixels in the resulting image. Then you can form an Ax = b system and solve for x to obtain the filter coefficients. With the coefficients available, I think you can find the sigma.

In the following example I'm using non-overlapping windows as shown to obtain the equations.

You don't have to know the size of the filter. If you use a larger size, the coefficients that are not relevant will be close to zero.

Your result image size is different than the input image, so i didn't use that image for following calculation. I use your input image and apply my own filter.

I tested this in Octave. You can quickly run it if you have Octave/Matlab. For Octave, you need to load the image package.

I'm using the following kernel to blur the image:

h =

   0.10963   0.11184   0.10963
   0.11184   0.11410   0.11184
   0.10963   0.11184   0.10963

When I estimate it using a window size 5, I get the following. As I said, the coefficients that are not relevant are close to zero.

g =

  9.5787e-015  -3.1508e-014  1.2974e-015  -3.4897e-015  1.2739e-014
  -3.7248e-014  1.0963e-001  1.1184e-001  1.0963e-001  1.8418e-015
  4.1825e-014  1.1184e-001  1.1410e-001  1.1184e-001  -7.3554e-014
  -2.4861e-014  1.0963e-001  1.1184e-001  1.0963e-001  9.7664e-014
  1.3692e-014  4.6182e-016  -2.9215e-014  3.1305e-014  -4.4875e-014

EDIT1: First of all, my apologies.

This approach doesn't really work in the practice. I've used the filt = conv2(a, h, 'same'); in the code. The resulting image data type in this case is double, whereas in the actual image the data type is usually uint8, so there's loss of information, which we can think of as noise. I simulated this with the minor modification filt = floor(conv2(a, h, 'same'));, and then I don't get the expected results.
The sampling approach is not ideal, because it's possible that it results in a degenerated system. Better approach is to use random sampling, avoiding the borders and making sure the entries in the b vector are unique. In the ideal case, as in my code, we are making sure the system Ax = b has a unique solution this way.
One approach would be to reformulate this as Mv = 0 system and try to minimize the squared norm of Mv under the constraint squared-norm v = 1, which we can solve using SVD. I could be wrong here, and I haven't tried this.
Another approach is to use the symmetry of the Gaussian kernel. Then a 3x3 kernel will have only 3 unknowns instead of 9. I think, this way we impose additional constraints on v of the above paragraph.
I'll try these out and post the results, even if I don't get the expected results.

EDIT2:

Using the LSE, we can find the filter coefficients as pinv(A'A)A'b. For completion, I'm adding a simple (and slow) LSE code.

Initial Octave Code:

clear all

im = double(imread('I2vxD.png'));

k = 5;
r = floor(k/2);

a = im(:, :, 1);  % take the red channel
h = fspecial('gaussian', [3 3], 5); % filter with a 3x3 gaussian
filt = conv2(a, h, 'same');

% use non-overlapping windows to for the Ax = b syatem
% NOTE: boundry error checking isn't performed in the code below
s = floor(size(a)/2);
y = s(1);
x = s(2);

w = k*k;
y1 = s(1)-floor(w/2) + r;
y2 = s(1)+floor(w/2);
x1 = s(2)-floor(w/2) + r;
x2 = s(2)+floor(w/2);

b = [];
A = [];
for y = y1:k:y2
  for x = x1:k:x2
    b = [b; filt(y, x)];
    f = a(y-r:y+r, x-r:x+r);
    A = [A; f(:)'];
  end
end

% estimated filter kernel
g = reshape(A\b, k, k)

LSE method:

clear all

im = double(imread('I2vxD.png'));

k = 5;
r = floor(k/2);

a = im(:, :, 1);  % take the red channel
h = fspecial('gaussian', [3 3], 5); % filter with a 3x3 gaussian
filt = floor(conv2(a, h, 'same'));

s = size(a);
y1 = r+2; y2 = s(1)-r-2;
x1 = r+2; x2 = s(2)-r-2;

b = [];
A = [];

for y = y1:2:y2
  for x = x1:2:x2
    b = [b; filt(y, x)];
    f = a(y-r:y+r, x-r:x+r);
    f = f(:)';
    A = [A; f];
  end
end

g = reshape(A\b, k, k) % A\b returns the least squares solution
%g = reshape(pinv(A'*A)*A'*b, k, k)