获得两个图像之间的高斯模糊的西格玛

Question

假设我有一个图像A，我用对其进行了高斯模糊处理，Sigam=3所以我得到了另一个图像B。如果给定A，B，是否有办法知道应用的西格玛？

进一步澄清：

图片A：

图片B：

我想编写一个接受A,B并返回的函数Sigma：

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

有什么建议么？

Answer 1

编辑1： 建议的方法实际上无法以其原始形式（即，仅对3 x 3内核使用9个方程式）工作，后来我意识到了这一点。有关说明，请参见下面的EDIT1，有关有效方法，请参见EDIT2。

编辑2： 如Humam所建议，我使用最小二乘估计（LSE）来找到系数。

我认为您可以在这种情况下通过求解线性方程组来估计滤波器内核。线性滤波器通过其系数加权窗口中的像素，然后取它们的总和并将该值分配给结果图像中窗口的中心像素。因此，对于3 x 3像

滤波后图像中的最终像素值

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)

a's原始图像窗口中的像素值在哪里？在这里，对于3 x 3滤波器，您有9个未知数，因此需要9个方程式。您可以在结果图像中使用9个像素获得这9个方程。然后，您可以形成一个Ax = b系统并求解x以获得滤波器系数。有了可用的系数，我想您可以找到sigma。

在下面的示例中，我使用所示的非重叠窗口来获取方程式。

您不必知道过滤器的大小。如果使用较大的大小，则不相关的系数将接近零。

您的结果图像大小与输入图像不同，因此我没有将该图像用于后续计算。我使用您的输入图像并应用自己的过滤器。

我在Octave中测试过。如果您有Octave / Matlab，则可以快速运行它。对于八度，您需要加载图像包。

我正在使用以下内核来模糊图像：

h =

   0.10963   0.11184   0.10963
   0.11184   0.11410   0.11184
   0.10963   0.11184   0.10963

当我使用窗口大小5估算时，得到以下结果。正如我所说，不相关的系数接近于零。

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： 首先，我道歉。

• 这种方法在实践中实际上并不起作用。我filt = conv2(a, h, 'same');在代码中使用了。在这种情况下，生成的图像数据类型为double，而在实际图像中，数据类型通常为uint8，因此会丢失信息，我们可以将其视为噪声。我用较小的修改对此进行了模拟filt = floor(conv2(a, h, 'same'));，然后没有得到预期的结果。
• 采样方法并不理想，因为它可能导致系统退化。更好的方法是使用随机采样，避免边界，并确保b向量中的条目唯一。在理想情况下，如我的代码所示，我们将确保系统Ax = b以这种方式具有唯一的解决方案。
• 一种方法是将其重新构造为Mv = 0系统，并尝试最小化Mv约束下的平方范数squared-norm v = 1，我们可以使用SVD解决此问题。我在这里可能是错的，而且我还没有尝试过。
• 另一种方法是使用高斯核的对称性。这样，一个3x3内核将只有3个未知数，而不是9个。我想，这样，我们对上v一段施加了额外的约束。
• 即使我没有得到预期的结果，我也会尝试这些并发布结果。

编辑2：

• 使用LSE，我们可以找到的滤波器系数pinv(A'A)A'b。为了完成，我添加了一个简单（缓慢）的LSE代码。

初始八度音阶代码：

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方法：

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)