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

Question

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

但是，我使用 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. 不是平方然后相加，然后是平方根 gx 和gy 值以获得总梯度，我只是求和 绝对 值。这大大提高了速度。

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

3. 编写嵌套的 for 循环来代替 np.sum(np.multiply(...)) 进行矩阵卷积。这最终变得稍微慢，这让我感到惊讶，因为自从 np.multiply 返回一个新数组以来，我认为使用 应该更快循环。我认为这可能是由于 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 notebook 时间仓促):

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 中要快得多。