如何提高索贝尔边缘检测器的效率

我正在写一个计算机视觉库 https://github.com/RoadKillCat/PiCamVision从头开始使用 Python 来使用rpi相机。目前，我已经实现了转换为greyscale以及其他一些基本的img在我的设备上运行速度相对较快的操作model B rpi3.

然而，我的边缘检测功能使用sobel操作员 （维基百科描述 https://en.wikipedia.org/wiki/Sobel_operator）尽管它确实有效，但比其他函数慢得多。这里是：

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. 代替squaring then adding, then square rooting the gx and gy值来获得总梯度，我只是sum the absolute价值观。这大大提高了速度。

2. 使用较低的resolution图像来自rpi相机。这显然是一种使这些操作运行得更快的简单方法，但它并不是那么可行，因为在最低可用分辨率下它仍然相当慢480x360这比相机的最大3280x2464.

3. 编写嵌套的 for 循环来执行以下操作matrix convolutions代替np.sum(np.multiply(...))。这最终有点slower我对此感到惊讶，因为从那时起np.multiply返回一个新数组，我认为这样做应该更快loops。我认为这可能是因为numpy大部分写在C或者新数组实际上并未存储，因此不需要很长时间，但我不太确定。

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

即使您正在构建自己的库，您确实绝对应该使用库进行卷积，它们将在后端用 C 或 Fortran 执行结果操作，这会快得多。

但如果您愿意的话，可以使用线性可分离滤波器自行完成。想法是这样的：

Image:

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

Sobel x kernel:

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

Result:

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]

然后您可以循环遍历索贝尔运算符的剩余列。或者，您甚至可以进一步进行，做我们刚刚做的同样的事情。这次对于垂直情况，您只需添加第一行和第三行以及第二行的两倍即可；或者，使用 numpy 加法：

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

你就完成了！我可能会在不久的将来在这里添加一些时间安排。对于一些粗略计算（即 1000x1000 图像上的仓促 Jupyter 笔记本计时）：

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

旧方法（双 for 循环）：每个循环 7.32 s ± 207 ms（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 images.

请注意，您可以这样做任何线性可分离滤波器！其中包括高斯滤波器。另请注意，一般来说，这需要大量操作。在 C 或 Fortran 或其他语言中，它通常只是实现为单行/列向量的两个卷积，因为最终它需要实际循环每个矩阵的每个元素；无论您只是将它们相加还是相乘，因此在 C 中将图像值相加并不比仅进行卷积更快。但循环遍历numpy数组非常慢，所以这个方法在 Python 中要快得多。