如何提高sobel边缘检测器的效率

Question

我正在用Python从头开始编写计算机视觉库，以与rpi相机配合使用。目前，我已经实现了转换为greyscale以及其他一些基本img操作，它们在我的上都运行得相对较快model B rpi3。

但是，我的使用sobel运算符的边缘检测功能（维基百科描述）比其他功能慢很多，尽管它确实起作用。这里是：

def sobel(img):
    xKernel = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
    yKernel = np.array([[-1,-2,-1],[0,0,0],[1,2,1]])
    sobelled = np.zeros((img.shape[0]-2, img.shape[1]-2, 3), dtype="uint8")
    for y in range(1, img.shape[0]-1):
        for x in range(1, img.shape[1]-1):
            gx = np.sum(np.multiply(img[y-1:y+2, x-1:x+2], xKernel))
            gy = np.sum(np.multiply(img[y-1:y+2, x-1:x+2], yKernel))
            g = abs(gx) + abs(gy) #math.sqrt(gx ** 2 + gy ** 2) (Slower)
            g = g if g > 0 and g < 255 else (0 if g < 0 else 255)
            sobelled[y-1][x-2] = g
    return sobelled

并使用greyscale猫的图片运行它：

我得到此答复，这似乎是正确的：

该库的应用程序，尤其是此功能，是在下棋机器人上进行的，其中边缘检测将有助于识别棋子的位置。问题在于>15运行需要几秒钟，这是一个重大问题，因为它将增加机器人大量移动的时间。

我的问题是：如何加快速度？

到目前为止，我已经尝试了一些方法：

1. 而不是squaringthen adding，然后square rooting使用gx和gy值获取总梯度，而我只是使用sum这些absolute值。这将速度提高了很多。

2. 使用相机的较低resolution图像rpi。显然，这是使这些操作更快运行的简单方法，但是它并不是那么可行，因为在最低可用分辨率下480x360，它仍然非常慢，而从相机的max最大值开始，该分辨率已大大降低3280x2464。

3. 编写嵌套的for循环来matrix convolutions代替np.sum(np.multiply(...))。最终速度稍慢，我很惊讶，因为np.multiply返回了一个新数组，所以我认为使用它应该会更快loops。我认为尽管这可能是由于numpy大多数情况下写入的事实，C或者实际上并未存储新数组，所以并不需要很长时间，但是我不太确定。

任何帮助将不胜感激-我认为主要的改进是要点3，即matrix乘法和求和。

Answer 1

即使您正在构建自己的库，您也绝对应该使用库进行卷积，它们将在后端用C或Fortran进行结果运算，这将大大加快速度。

但是，如果您愿意，可以自己做，使用线性可分离滤波器。这是想法：

图片：

1 2 3 4 5
2 3 4 5 1
3 4 5 1 2

Sobel x内核：

-1 0 1
-2 0 2
-1 0 1

结果：

8, 3, -7

在卷积的第一个位置，您将计算9个值。首先，为什么？您永远不会添加中间列，不必费心将其相乘。但这不是线性可分离滤波器的重点。这个想法很简单。将内核放在第一个位置时，将第三列乘以[1, 2, 1]。但是随后的两步，您将第三列乘以[-1, -2, -1]。真是浪费！您已经计算过了，现在只需要取反即可。这就是线性可分离滤波器的想法。请注意，您可以将过滤器分解为两个向量的矩阵外积：

[1]
[2]  *  [-1, 0, 1]
[1]

将外部乘积在这里得到相同的矩阵。因此，这里的想法是将操作分为两部分。首先将整个图像乘以行向量，然后乘以列向量。以行向量

-1 0 1

在整个图像上，我们最终得到

2  2  2
2  2 -3
2 -3 -3

然后将列向量传递给乘法和求和，我们再次得到

8, 3, -7

另一种可能有用或无效的妙招（取决于您在内存和效率之间的权衡）：

请注意，在单行乘法中，您忽略中间值，而只从左边的值中减去右边的值。这意味着您正在有效地执行以下操作减去这两个图像：

3 4 5     1 2 3
4 5 1  -  2 3 4
5 1 2     3 4 5

如果将图像的前两列剪掉，则会得到左矩阵，如果将后两列剪掉，则会得到右矩阵。因此，您只需将卷积的第一部分计算为

result_h = img[:,2:] - img[:,:-2]

然后，您可以遍历sobel运算符的其余列。或者，您甚至可以继续进行，然后做我们刚才做的相同的事情。这次是垂直情况，您只需要添加第一行和第三行，第二行加倍即可；或者，使用numpy加法：

result_v = result_h[:-2] + result_h[2:] + 2*result_h[1:-1]

大功告成！我可能会在不久的将来在此处添加一些时间。对于包络的某些计算（例如，在1000x1000图像上草草的Jupyter笔记本计时）：

新方法（图像总和）：每个循环8.18 ms±399 µs（平均±标准偏差，共运行7次，每个循环100个）

旧方法（双循环）：每个循环7.32 s±207毫秒（平均±标准偏差，共运行7次，每个循环1次）

是的，您没看错：1000倍加速。

这是比较两者的一些代码：

import numpy as np

def sobel_x_orig(img):
    xKernel = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
    sobelled = np.zeros((img.shape[0]-2, img.shape[1]-2))
    for y in range(1, img.shape[0]-1):
        for x in range(1, img.shape[1]-1):
            sobelled[y-1, x-1] = np.sum(np.multiply(img[y-1:y+2, x-1:x+2], xKernel))
    return sobelled

def sobel_x_new(img):
    result_h = img[:,2:] - img[:,:-2]
    result_v = result_h[:-2] + result_h[2:] + 2*result_h[1:-1]
    return result_v

img = np.random.rand(1000, 1000)
sobel_new = sobel_x_new(img)
sobel_orig = sobel_x_orig(img)

assert (np.abs(sobel_new-sobel_orig) < 1e-12).all()

当然，1e-12有一些严重的公差，但这是每个元素的，所以应该可以。但是我也有float图像，您在图像上当然会有更大的差异uint8。

请注意，您可以对任何线性可分离滤波器执行此操作！包括高斯滤波器。还要注意，通常，这需要很多操作。在C或Fortran或任何其他语言中，通常仅将其实现为单行/列向量的两个卷积，因为最后，它实际上仍然需要遍历每个矩阵的每个元素；不管您是将它们相加还是相乘，因此在C中用这种方法添加图像值并不比卷积要快。但是遍历numpy数组非常慢，因此在Python中这种方法要快得多。