大尺寸 Sobel 滤波器核
Posted
技术标签:
【中文标题】大尺寸 Sobel 滤波器核【英文标题】：Sobel filter kernel of large size 【发布时间】：2012-03-22 23:56:46 【问题描述】：我正在使用大小为 3x3 的 sobel 滤波器来计算图像导数。看网上的一些文章，尺寸 5x5 和 7x7 的 sobel 过滤器的内核似乎也很常见，但我找不到它们的内核值。
有人可以告诉我尺寸为 5x5 和 7x7 的 sobel 滤波器的内核值吗？此外，如果有人可以分享一种生成内核值的方法，那将非常有用。
提前致谢。
【问题讨论】：
【参考方案1】：2018 年 4 月 23 日更新：似乎下面链接中定义的内核不是真正的 Sobel 内核（适用于 5x5 及更高版本）——它们可能在边缘检测方面做得不错，但它们不应该是称为 Sobel 核。请参阅Daniel’s answer 以获得更准确和全面的摘要。 （我将把这个答案留在这里，因为 (a) 它是从各个地方链接到的，并且 (b) 接受的答案不容易被删除。）
Google 似乎提供了很多结果，例如
http://rsbweb.nih.gov/nih-image/download/user-macros/slowsobel.macro 建议 3x3、5x5、7x7 和 9x9 使用以下内核：
3x3：
1 0 -1
2 0 -2
1 0 -1
5x5：
2 1 0 -1 -2
3 2 0 -2 -3
4 3 0 -3 -4
3 2 0 -2 -3
2 1 0 -1 -2
7x7：
3 2 1 0 -1 -2 -3
4 3 2 0 -2 -3 -4
5 4 3 0 -3 -4 -5
6 5 4 0 -4 -5 -6
5 4 3 0 -3 -4 -5
4 3 2 0 -2 -3 -4
3 2 1 0 -1 -2 -3
9x9：
4 3 2 1 0 -1 -2 -3 -4
5 4 3 2 0 -2 -3 -4 -5
6 5 4 3 0 -3 -4 -5 -6
7 6 5 4 0 -4 -5 -6 -7
8 7 6 5 0 -5 -6 -7 -8
7 6 5 4 0 -4 -5 -6 -7
6 5 4 3 0 -3 -4 -5 -6
5 4 3 2 0 -2 -3 -4 -5
4 3 2 1 0 -1 -2 -3 -4
【讨论】：
谢谢，有生成器函数的链接吗？ 我不认为这是一个。 Sobel 实际上只为 3x3 定义，更大的内核似乎是在 ad hoc 的基础上确定的。它们只是微分器，因此为任何大小的内核生成系数都应该相当容易。 抱歉，这些不是 Sobel 内核。 @Daniel：你愿意扩展一下吗？ 检查我的答案，写了一段时间:-)【参考方案2】：其他来源似乎对较大的内核给出了不同的定义。例如，Intel IPP library 将 5x5 内核设为
1 2 0 -2 -1
4 8 0 -8 -4
6 12 0 -12 -6
4 8 0 -8 -4
1 2 0 -2 -1
直观地说，这对我来说更有意义，因为您更加关注靠近中心的元素。它还具有 3x3 内核的自然定义，易于扩展以生成更大的内核。也就是说，在我的简短搜索中，我发现了 5x5 内核的 3 种不同定义 - 所以我怀疑（正如 Paul 所说）较大的内核是临时的，所以这绝不是明确的答案。
3x3 内核是平滑内核和梯度内核的外积，在 Matlab 中类似于
sob3x3 = [ 1 2 1 ]' * [1 0 -1]
可以通过将 3x3 内核与另一个平滑内核进行卷积来定义更大的内核
sob5x5 = conv2( [ 1 2 1 ]' * [1 2 1], sob3x3 )
您可以重复该过程以获得逐渐变大的内核
sob7x7 = conv2( [ 1 2 1 ]' * [1 2 1], sob5x5 )
sob9x9 = conv2( [ 1 2 1 ]' * [1 2 1], sob7x7 )
...
还有很多其他的写法，但我认为这最能准确地解释正在发生的事情。基本上，您从一个方向的平滑内核和另一个方向的导数的有限差分估计开始，然后只应用平滑，直到获得所需的内核大小。
因为它只是一系列卷积，所以所有好的属性（交换性、关联性等）都对您的实现有用。例如，您可以简单地将 5x5 内核分离为其平滑和导数组件：
sob5x5 = conv([1 2 1],[1 2 1])' * conv([1 2 1],[-1 0 1])
请注意，为了成为“正确”的导数估计器，3x3 Sobel 应按 1/8 倍缩放：
sob3x3 = 1/8 * [ 1 2 1 ]' * [1 0 -1]
每个较大的内核都需要额外缩放 1/16 倍（因为平滑内核未归一化）：
sob5x5 = 1/16 * conv2( [ 1 2 1 ]' * [1 2 1], sob3x3 )
sob7x7 = 1/16 * conv2( [ 1 2 1 ]' * [1 2 1], sob5x5 )
...
【讨论】：
但是为什么 3x3 内核的比例因子是 1/8，而 5x5 的比例因子是 1/16？ 一维有限差分内核为1/2 * [1 0 -1]，一维平滑内核为1/4 * [ 1 2 1 ]。所以 3x3 需要1/2 * 1/4 = 1/8 的比例，而较大的内核在从 1 维移动到 2 维时需要1/4 * 1/4 = 1/16。
如果您仍然对该主题感兴趣，如果您对我的任意大小和旋转的 Sobel 内核的解决方案提供反馈，我将不胜感激。【参考方案3】：
根据@Paul R 给出的示例，我快速破解了一个算法来生成任何大于 1 的奇数大小的 Sobel 内核：
public static void CreateSobelKernel(int n, ref float[][] Kx, ref float[][] Ky)
int side = n * 2 + 3;
int halfSide = side / 2;
for (int i = 0; i < side; i++)
int k = (i <= halfSide) ? (halfSide + i) : (side + halfSide - i - 1);
for (int j = 0; j < side; j++)
if (j < halfSide)
Kx[i][j] = Ky[j][i] = j - k;
else if (j > halfSide)
Kx[i][j] = Ky[j][i] = k - (side - j - 1);
else
Kx[i][j] = Ky[j][i] = 0;
希望对你有帮助。
【讨论】：
【参考方案4】：谢谢大家，我将尝试@Adam Bowen 的第二个变体，为 Sobel5x5、7x7、9x9 获取 C# 代码......这个变体的矩阵生成（可能有错误，如果你发现错误或可以优化代码 - 写它那里）：
static void Main(string[] args)
float[,] Sobel3x3 = new float[,]
-1, 0, 1,
-2, 0, 2,
-1, 0, 1;
float[,] Sobel5x5 = Conv2DforSobelOperator(Sobel3x3);
float[,] Sobel7x7 = Conv2DforSobelOperator(Sobel5x5);
Console.ReadKey();
public static float[,] Conv2DforSobelOperator(float[,] Kernel)
if (Kernel == null)
throw new Exception("Kernel = null");
if (Kernel.GetLength(0) != Kernel.GetLength(1))
throw new Exception("Kernel matrix must be Square matrix!");
float[,] BaseMatrix = new float[,]
1, 2, 1,
2, 4, 2,
1, 2, 1;
int KernelSize = Kernel.GetLength(0);
int HalfKernelSize = KernelSize / 2;
int OutSize = KernelSize + 2;
if ((KernelSize & 1) == 0) // Kernel_Size must be: 3, 5, 7, 9 ...
throw new Exception("Kernel size must be odd (3x3, 5x5, 7x7...)");
float[,] Out = new float[OutSize, OutSize];
float[,] InMatrix = new float[OutSize, OutSize];
for (int x = 0; x < BaseMatrix.GetLength(0); x++)
for (int y = 0; y < BaseMatrix.GetLength(1); y++)
InMatrix[HalfKernelSize + x, HalfKernelSize + y] = BaseMatrix[x, y];
for (int x = 0; x < OutSize; x++)
for (int y = 0; y < OutSize; y++)
for (int Kx = 0; Kx < KernelSize; Kx++)
for (int Ky = 0; Ky < KernelSize; Ky++)
int X = x + Kx - HalfKernelSize;
int Y = y + Ky - HalfKernelSize;
if (X >= 0 && Y >= 0 && X < OutSize && Y < OutSize)
Out[x, y] += InMatrix[X, Y] * Kernel[KernelSize - 1 - Kx, KernelSize - 1 - Ky];
return Out;
Results (NormalMap) 或 it copy there，这里方法 - №2，@Paul R 方法 - №1。现在我使用最后一个，因为它可以提供更流畅的结果，并且使用this 代码很容易生成内核。
【讨论】：
【参考方案5】：任意 Sobel 核大小和角度的完整解决方案
tl;dr：跳到“示例”部分
添加另一个解决方案，扩展 this document（它的质量不是特别高，但从第 2 页底部开始显示一些可用的图形和矩阵）。
目标
我们要做的是估计图像在位置 (x,y) 处的局部梯度。梯度是由x和y方向的分量gx和gy组成的向量。
现在，假设我们想根据我们的像素 (x,y) 及其邻居作为核操作（3x3、5x5 或任何大小）来近似梯度。
解决思路
我们可以通过对所有相邻中心对在梯度方向上的投影求和来近似梯度。 （Sobel 的核只是一种特殊的加权不同贡献的方法，Prewitt 基本上也是如此）。
3x3 的显式中间步骤
这是本地图像，中心像素 (x,y) 标记为 'o' (center)
a b c
d o f
g h i
假设我们想要正 x 方向的梯度。 x 正方向的单位向量是 (1,0) [我稍后会使用正 y 方向向下的约定，即 (0,1)，并且 (0,0) 是图像的左上角) .]
从o到f（简称'of'）的向量是(1,0)。 “of”方向的梯度为 (f - o) / 1（此处像素处的图像值表示 f 减去中心 o 处的值，除以这些像素之间的距离）。如果我们通过点积将特定邻居梯度的单位向量投影到我们想要的梯度方向 (1,0) 上，我们得到 1。这是一个包含所有邻居贡献的小表格，从更简单的情况开始。注意对于对角线，它们的距离是sqrt2，对角线方向的单位向量是1/sqrt2 * (+/-1, +/-1)
f: (f-o)/1 * 1
d: (d-o)/1 * -1 because (-1, 0) dot (1, 0) = -1
b: (b-o)/1 * 0 because (0, -1) dot (1, 0) = 0
h: (h-o)/1 * 0 (as per b)
a: (a-o)/sqrt2 * -1/sqrt2 distance is sqrt2, and 1/sqrt2*(-1,-1) dot (1,0) = -1/sqrt2
c: (c-o)/sqrt2 * +1/sqrt2 ...
g: (g-o)/sqrt2 * -1/sqrt2 ...
i: (i-o)/sqrt2 * +1/sqrt2 ...
编辑澄清： 1/sqrt(2) 有两个因数，原因如下：
我们感兴趣的是在特定方向（这里是x）对梯度的贡献，所以我们需要将从中心像素到相邻像素的方向梯度投影到该方向上我们感兴趣。这是通过在各个方向上取单位向量的标量积来实现的，它引入了第一个因子 1/L（这里是对角线的 1/sqrt(2)）。
梯度测量一点的无穷小变化，我们用有限差分来近似。就线性方程而言，m = (y2-y1)/(x2-x1)。出于这个原因，从中心像素到相邻像素 (y2-y1) 的值差必须分布在它们的距离上（对应于 x2-x1），以便获得每个距离单位的上升单位。这产生了 1/L 的第二个因子（这里对角线为 1/sqrt(2)）
好的，现在我们知道贡献了。让我们通过组合相反的像素贡献对来简化这个表达式。我将从 d 和 f 开始：
(f-o)/1 * 1 + (d-o)/1 * -1
= f - o - (d - o)
= f - d
现在是第一个对角线：
(c-o)/sqrt2 * 1/sqrt2 + (g-o)/sqrt2 * -1/sqrt2
= (c - o)/2 - (g - o)/2
= (c - g)/2
第二条对角线的贡献是 (i - a)/2。垂直方向贡献为零。请注意，中心像素“o”的所有贡献都消失了。
我们现在已经计算了像素 (x,y) 处所有最近邻对正 x 方向梯度的贡献，因此我们对 x 方向梯度的总近似就是它们的总和：
gx(x,y) = f - d + (c - g)/2 + (i - a)/2
我们可以通过使用卷积核来获得相同的结果，其中系数被写入相应的相邻像素的位置：
-1/2 0 1/2
-1 0 1
-1/2 0 1/2
如果您不想处理分数，请将其乘以 2 并获得著名的 Sobel 3x3 内核。
-1 0 1
G_x = -2 0 2
-1 0 1
乘以二只是为了得到方便的整数。输出图像的缩放基本上是任意的，大多数时候您将其标准化为您的图像范围，无论如何（以获得清晰可见的结果）。
通过与上述相同的推理，您可以通过将邻居贡献投影到正 y 方向 (0,1) 上的单位向量上来获得垂直梯度 gy 的内核
-1 -2 -1
G_y = 0 0 0
1 2 1
任意大小内核的公式
如果你想要 5x5 或更大的内核，你只需要注意距离，例如
A B 2 B A
B C 1 C B
2 1 - 1 2
B C 1 C B
A B 2 B A
在哪里
A = 2 * sqrt2
B = sqrt5
C = sqrt2.
如果连接任意两个像素的向量的长度为 L，则该方向上的单位向量的前置因子为 1/L。出于这个原因，任何像素'k'对（比如说）x梯度（1,0）的贡献可以简化为“（平方距离上的值差）乘以（非归一化方向向量'ok'与梯度向量的DotProduct , 例如 (1,0) )"
gx_k = (k - o)/(pixel distance^2) ['ok' dot (1,0)].
因为连接向量与x单位向量的点积选择了对应的向量入口，所以位置k对应的G_x核入口正好
i / (i*i + j*j)
其中 i 和 j 是从中心像素到像素 k 在 x 和 y 方向上的步数。在上述 3x3 计算中，像素 'a' 将具有 i = -1（左侧 1），j = -1（顶部 1），因此 'a' 内核条目为 -1 / (1 + 1 ) = -1/2。
G_y 内核的条目是
j/(i*i + j*j).
如果我想要内核的整数值，请按照以下步骤操作：
检查输出图像的可用范围 通过应用浮点内核计算可能的最高结果（即假设所有正内核条目下的最大输入值，因此输出值为（所有正内核值的总和）*（最大可能的输入图像值）。如果您有签名输入,你也需要考虑负值。最坏的情况是所有正值的总和+负条目的所有abs值的总和（如果在正数下的最大输入，在负数下的最大输入）。编辑：所有的总和abs 值也被恰当地称为内核的权重 计算内核允许的最大放大倍数（不溢出输出图像范围） 对于浮点内核的所有整数倍数（从 2 到最大值以上）：检查哪个具有最小的绝对舍入误差总和并使用此内核总之：
Gx_ij = i / (i*i + j*j)
Gy_ij = j / (i*i + j*j)
其中 i,j 是从中心算起的内核中的位置。根据需要缩放内核条目以获得整数（或至少接近近似值）。
这些公式适用于所有内核大小。
示例
-2/8 -1/5 0 1/5 2/8 -5 -4 0 4 5
-2/5 -1/2 0 1/2 2/5 -8 -10 0 10 8
G_x (5x5) -2/4 -1/1 0 1/1 2/4 (*20) = -10 -20 0 20 10
-2/5 -1/2 0 1/2 2/5 -8 -10 0 10 8
-2/8 -1/5 0 1/5 2/8 -5 -4 0 4 5
请注意，浮点表示法中 5x5 内核的中心 3x3 像素只是 3x3 内核，即较大的内核表示具有附加但权重较低的数据的持续近似值。这继续到更大的内核大小：
-3/18 -2/13 -1/10 0 1/10 2/13 3/18
-3/13 -2/8 -1/5 0 1/5 2/8 3/13
-3/10 -2/5 -1/2 0 1/2 2/5 3/10
G_x (7x7) -3/9 -2/4 -1/1 0 1/1 2/4 3/9
-3/10 -2/5 -1/2 0 1/2 2/5 3/10
-3/13 -2/8 -1/5 0 1/5 2/8 3/13
-3/18 -2/13 -1/10 0 1/10 2/13 3/18
此时，精确的整数表示变得不切实际。
据我所知（无法访问原始论文），“Sobel”部分对贡献进行了适当的加权。 Prewitt解可以省略距离权重，只在内核中适当输入i和j。
奖励：任意方向的 Sobel 内核
所以我们可以近似图像梯度的 x 和 y 分量（它实际上是一个向量，如开头所述）。通过将梯度向量投影到 alpha-梯度单位向量上，可以获得任意方向 alpha 上的梯度（在数学上测量为正，在这种情况下为顺时针方向，因为正 y 向下）。
alpha 单位向量是 (cos alpha, sin alpha)。对于 alpha = 0°，您可以获得 gx 的结果，对于 alpha = 90°，您可以获得 gy。
g_alpha = (alpha-unit vector) dot (gx, gy)
= (cos a, sin a) dot (gx, gy)
= cos a * gx + sin a * gy
如果您费心将 gx 和 gy 写为邻居贡献的总和，您会意识到您可以通过适用于同一邻居像素的术语对生成的长表达式进行分组，然后将其重写为带有条目的单个卷积核
G_alpha_ij = (i * cos a + j * sin a)/(i*i + j*j)
如果您想要最接近的整数近似值，请按照上述步骤操作。
【讨论】：
感谢您的全面报道——这显然比我的回答所依据的 NIH 文件中的近似值要严格得多。 @Daniel 我对此有不同的答案。根据 Sobel 的描述，我不确定所有根是如何脱落的，以及为什么 G(1,2)=-4 和 G(2,1)=-8 的权重不同。它们都与测量导数的中心等距，并且应该具有相同的权重，不是吗？这是我尝试导出 NxN 2D Sobel 算子时得到的结果，也许我在某个地方出错了：reddit.com/r/computervision/comments/8e234p 嗨！ -4 和 -8 不是由于不同的距离权重，而是由于我所说的 i 和 j，即距离中心元素梯度方向的距离。如果条目相等，您将测量对角梯度（在这种情况下，水平和垂直边界会从内核中引发相同的响应强度）。第二个根因子应该来自于以反距离加权像素贡献，第一个来自对“梯度贡献向量”的长度进行归一化。 哦哦。 “梯度方向的距离......”我明白你在说什么。因此为什么中心行/列为零（我希望投影到梯度方向）。我能够通过手动完成所有数学运算部分成功地重新导出 3x3 案例，但我有 sqrt(2) 和 2 而不是 1 和 2。这只是为了近似还是我错过了一步？当我试图将它们分开时，一切都取消了，我想出了 Prewitt。完全感谢。我试图对它进行一些小改动以换取其他东西，这真的很有帮助。 有一点需要注意：如果图像中的渐变是均匀的，则距离因子越小，就可以平衡这样一个事实，即在渐变方向上越远，值差异就越大。在这种特殊情况下，梯度方向上的所有像素都贡献了相同的数量。在更常见的嘈杂或不完全均匀梯度的情况下，更大的内核实际上可以提高平滑度。就我而言，到目前为止，5x5 通常就足够了。【参考方案6】：Sobel 梯度滤波器生成器
（这个答案指的是上面@Daniel给出的analysis。）
Gx[i,j] = i / (i*i + j*j)
Gy[i,j] = j / (i*i + j*j)
这是一个重要的结果，并且比在original paper 中找到的解释更好。它应该写在Wikipedia 或其他地方，因为它似乎也优于互联网上任何其他关于该问题的讨论。
然而，正如声称的那样，对于尺寸大于 5*5 的过滤器，整数值表示实际上是不切实际的。使用 64 位整数，可以精确表示最大为 15*15 的 Sobel 滤波器。
这是前四个；结果应该除以“权重”，以便如下图的图像区域的梯度被归一化为值 1。
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
Gx(3)：
-1/2 0/1 1/2 -1 0 1
-1/1 0 1/1 * 2 = -2 0 2
-1/2 0/1 1/2 -1 0 1
weight = 4 weight = 8
Gx(5)：
-2/8 -1/5 0/4 1/5 2/8 -5 -4 0 4 5
-2/5 -1/2 0/1 1/2 2/5 -8 -10 0 10 8
-2/4 -1/1 0 1/1 2/4 * 20 = -10 -20 0 20 10
-2/5 -1/2 0/1 1/2 2/5 -8 -10 0 10 8
-2/8 -1/5 0/4 1/5 2/8 -5 -4 0 4 5
weight = 12 weight = 240
Gx(7)：
-3/18 -2/13 -1/10 0/9 1/10 2/13 3/18 -130 -120 -78 0 78 120 130
-3/13 -2/8 -1/5 0/4 1/5 2/8 3/13 -180 -195 -156 0 156 195 180
-3/10 -2/5 -1/2 0/1 1/2 2/5 3/10 -234 -312 -390 0 390 312 234
-3/9 -2/4 -1/1 0 1/1 2/4 3/9 * 780 = -260 -390 -780 0 780 390 260
-3/10 -2/5 -1/2 0/1 1/2 2/5 3/10 -234 -312 -390 0 390 312 234
-3/13 -2/8 -1/5 0/4 1/5 2/8 3/13 -180 -195 -156 0 156 195 180
-3/18 -2/13 -1/10 0/9 1/10 2/13 3/18 -130 -120 -78 0 78 120 130
weight = 24 weight = 18720
Gx(9)：
-4/32 -3/25 -2/20 -1/17 0/16 1/17 2/20 3/25 4/32 -16575 -15912 -13260 -7800 0 7800 13260 15912 16575
-4/25 -3/18 -2/13 -1/10 0/9 1/10 2/13 3/18 4/25 -21216 -22100 -20400 -13260 0 13260 20400 22100 21216
-4/20 -3/13 -2/8 -1/5 0/4 1/5 2/8 3/13 4/20 -26520 -30600 -33150 -26520 0 26520 33150 30600 26520
-4/17 -3/10 -2/5 -1/2 0/1 1/2 2/5 3/10 4/17 -31200 -39780 -53040 -66300 0 66300 53040 39780 31200
-4/16 -3/9 -2/4 -1/1 0 1/1 2/4 3/9 4/16 * 132600 = -33150 -44200 -66300 -132600 0 132600 66300 44200 33150
-4/17 -3/10 -2/5 -1/2 0/1 1/2 2/5 3/10 4/17 -31200 -39780 -53040 -66300 0 66300 53040 39780 31200
-4/20 -3/13 -2/8 -1/5 0/4 1/5 2/8 3/13 4/20 -26520 -30600 -33150 -26520 0 26520 33150 30600 26520
-4/25 -3/18 -2/13 -1/10 0/9 1/10 2/13 3/18 4/25 -21216 -22100 -20400 -13260 0 13260 20400 22100 21216
-4/32 -3/25 -2/20 -1/17 0/16 1/17 2/20 3/25 4/32 -16575 -15912 -13260 -7800 0 7800 13260 15912 16575
weight = 40 weight = 5304000
下面附加的 Ruby 程序将计算任意大小的 Sobel 过滤器和相应的权重，尽管整数值过滤器不太可能对大于 15*15 的尺寸有用。
#!/usr/bin/ruby
# Sobel image gradient filter generator
# by <ian_bruce@mail.ru> -- Sept 2017
# reference:
# https://***.com/questions/9567882/sobel-filter-kernel-of-large-size
if (s = ARGV[0].to_i) < 3 || (s % 2) == 0
$stderr.puts "invalid size"
exit false
end
s /= 2
n = 1
# find least-common-multiple of all fractional denominators
(0..s).each |j|
(1..s).each |i|
d = i*i + j*j
n = n.lcm(d / d.gcd(i))
fw1 = format("%d/%d", s, 2*s*s).size + 2
fw2 = format("%d", n).size + 2
weight = 0
s1 = ""
s2 = ""
(-s..s).each |y|
(-s..s).each |x|
i, j = x, y # "i, j = y, x" for transpose
d = i*i + j*j
if (i != 0)
if (n * i % d) != 0 # this should never happen
$stderr.puts "inexact division: #n * #i / ((#i)^2 + (#j)^2)"
exit false
end
w = n * i / d
weight += i * w
else
w = 0
end
s1 += "%*s" % [fw1, d > 0 ? "%d/%d" % [i, d] : "0"]
s2 += "%*d" % [fw2, w]
s1 += "\n" ; s2 += "\n"
f = n.gcd(weight)
puts s1
puts "\nweight = %d%s" % [weight/f, f < n ? "/%d" % (n/f) : ""]
puts "\n* #n =\n\n"
puts s2
puts "\nweight = #weight"
【讨论】：
感谢您的启发，这里是c++中的例程（整数和浮点版本）：gist.github.com/eudoxos/b88f36e27eb782c6c60f25481c642aaf【参考方案7】：这是一个使用 numpy 和 @Daniel 答案使用 python 3 制作的简单解决方案。
def custom_sobel(shape, axis):
"""
shape must be odd: eg. (5,5)
axis is the direction, with 0 to positive x and 1 to positive y
"""
k = np.zeros(shape)
p = [(j,i) for j in range(shape[0])
for i in range(shape[1])
if not (i == (shape[1] -1)/2. and j == (shape[0] -1)/2.)]
for j, i in p:
j_ = int(j - (shape[0] -1)/2.)
i_ = int(i - (shape[1] -1)/2.)
k[j,i] = (i_ if axis==0 else j_)/float(i_*i_ + j_*j_)
return k
它像这样返回内核 (5,5)：
Sobel x:
[[-0.25 -0.2 0. 0.2 0.25]
[-0.4 -0.5 0. 0.5 0.4 ]
[-0.5 -1. 0. 1. 0.5 ]
[-0.4 -0.5 0. 0.5 0.4 ]
[-0.25 -0.2 0. 0.2 0.25]]
Sobel y:
[[-0.25 -0.4 -0.5 -0.4 -0.25]
[-0.2 -0.5 -1. -0.5 -0.2 ]
[ 0. 0. 0. 0. 0. ]
[ 0.2 0.5 1. 0.5 0.2 ]
[ 0.25 0.4 0.5 0.4 0.25]]
如果有人知道在 python 中执行此操作的更好方法，请告诉我。我还是个新手;)
【讨论】：
【参考方案8】：TL;DR：改用高斯导数运算符。
作为Adam Bowen explained in his answer，Sobel 核是沿一个轴的平滑和沿另一轴的中心差分导数的组合：
sob3x3 = [1 2 1]' * [1 0 -1]
平滑增加了正则化（降低对噪声的敏感性）。
（我在这篇文章中省略了所有因素1/8，as did Sobel himself，这意味着运算符确定了缩放的导数。另外，* 在这篇文章中总是意味着卷积。）
让我们概括一下：
deriv_kernel = smoothing_kernel * d/dx
卷积的特性之一是
d/dx f = d/dx * f
也就是说，用元素导数算子对图像进行卷积会产生图像的导数。还要注意卷积是可交换的，
deriv_kernel = d/dx * smoothing_kernel = d/dx smoothing_kernel
也就是说，导数核是平滑核的导数。
请注意，通过卷积将这样的内核应用于图像：
image * deriv_kernel = image * smoothing_kernel * d/dx = d/dx (image * smoothing_kernel)
也就是说，通过这个广义的、理想化的导数核，我们可以计算平滑图像的真导数。索贝尔核当然不是这种情况，因为它使用了对导数的中心差分近似。
但是选择一个更好的smoothing_kernel，这个是可以实现的。高斯核是这里的理想选择，因为它在空间域的紧凑性（小内核占用空间）与频域的紧凑性（良好的平滑性）之间提供了最佳折衷。此外，高斯是完全各向同性和可分离的。使用高斯导数核会产生最佳的正则化导数算子。
因此，如果您正在寻找更大的 Sobel 算子，因为您需要更多的正则化，请改用高斯导数算子。
让我们进一步分析一下 Sobel 内核。
平滑内核是三角形的，样本[1 2 1]。这是一个三角函数，通过采样得到这三个值：
2 + x , if -2 < x < 0
h = 2 , if x = 0
2 - x , if 0 < x < 2
它的导数是：
1 , if -2 < x < 0
d/dx h = 0 , if x = 0 (not really, but it's the sensible solution)
-1 , if 0 < x < 2
因此，我们可以看到中心差分导数近似可以看作是用于平滑的同一三角函数的解析导数的采样。因此我们有：
sob3x3 = [1 2 1]' * d/dx [1 2 1] = d/dx ( [1 2 1]' * [1 2 1] )
所以，如果你想让这个内核更大，只需放大平滑内核即可：
sob5x5 = d/dx ( [1 2 3 2 1]' * [1 2 3 2 1] ) = [1 2 3 2 1]' * [1 1 0 -1 -1]
sob7x7 = d/dx ( [1 2 3 4 3 2 1]' * [1 2 3 4 3 2 1] ) = [1 2 3 4 3 2 1]' * [1 1 1 0 -1 -1 -1]
这与advice given by Adam Bowen 完全不同，advice given by Adam Bowen 建议将内核与沿每个维度的 3-tab 三角形内核进行卷积：[1 2 1] * [1 2 1] = [1 4 6 4 1] 和 [1 2 1] * [1 0 -1] = [1 2 0 -2 -1]。请注意，由于中心极限定理，将这个三角形核与其自身进行卷积会导致滤波器更接近高斯。我们通过与自身重复卷积创建的内核越大，我们就越接近这个高斯。所以，不使用这种方法，还不如直接对高斯函数进行采样。
Daniel has a long post 在其中他建议以另一种方式扩展 Sobel 内核。这里的平滑核的形状与高斯近似不同，我没有尝试研究它的性质。
请注意，这三个可能的 Sobel 内核扩展都不是真正的 Sobel 内核，因为 Sobel 内核明确地是一个 3x3 内核（see an historical note by Sobel about his operator，他从未真正发布过）。
还要注意，我并不是在提倡从这里派生的扩展 Sobel 内核。使用高斯导数！
【讨论】：
【参考方案9】：Daniel's answer的Matlab实现：
kernel_width = 9;
halfway = floor(kernel_width/2);
step = -halfway : halfway;
i_matrix = repmat(step,[kernel_width 1]);
j_matrix = i_matrix';
gx = i_matrix ./ ( i_matrix.*i_matrix + j_matrix.*j_matrix );
gx(halfway+1,halfway+1) = 0; % deals with NaN in middle
gy = gx';
【讨论】：
【参考方案10】：我做了一个Daniel's answer 的 Python NumPy 实现。它似乎比 Joao Ponte's implementation 快 3 倍左右。
def calc_sobel_kernel(target_shape: tuple[int, int]):
assert target_shape[0] % 2 != 0
assert target_shape[1] % 2 != 0
gx = np.zeros(target_shape, dtype=np.float32)
gy = np.zeros(target_shape, dtype=np.float32)
indices = np.indices(target_shape, dtype=np.float32)
cols = indices[0] - target_shape[0] // 2
rows = indices[1] - target_shape[1] // 2
squared = cols ** 2 + rows ** 2
np.divide(cols, squared, out=gy, where=squared!=0)
np.divide(rows, squared, out=gx, where=squared!=0)
return gx, gy
【讨论】：
以上是关于大尺寸 Sobel 滤波器核的主要内容，如果未能解决你的问题，请参考以下文章