如何提高sobel边d的效率

Question

我正在用Python从头开始编写一个computer vision library来使用rpi相机。目前，我已经实现了对greyscale的转换和其他一些基本的img操作，它们在我的model Brpi3上运行得比较快。在

然而，我使用sobel运算符（wikipedia description）的边缘检测函数比其他函数慢得多，尽管它确实可以工作。这里是：

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

用这张猫的图片运行它：

我得到的回答似乎是正确的：

这个库的应用，尤其是这个功能，是在一个下棋机器人上，边缘检测将帮助识别棋子的位置。问题是它需要>15秒来运行，这是一个重要的问题，因为它会增加机器人移动所需的时间。在

我的问题是：我怎样才能加快速度？

到目前为止，我尝试了几件事：

1. 不是squaring然后adding，然后square rooting得到总梯度，我只需要sum和{}值。这大大提高了速度。

2. 使用来自rpi相机的较低resolution图像。这显然是一个简单的方法，使这些操作运行得更快，但它并不是真正可行的，因为它仍然是相当缓慢的最低可用分辨率480x360，这是从相机的最大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

Sobelx内核：

^{pr2}$

结果：

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笔记本计时）：

new method (sums of the images): 8.18 ms ± 399 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

old method (double for-loop):7.32 s ± 207 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

是的，你读对了：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中要快得多。在