如何有效地扫描二维 numpy 数组？答案

【问题标题】：How to effectively scan 2d numpy array?如何有效地扫描二维 numpy 数组？
【发布时间】：2013-03-22 14:55:37
【问题描述】：

我需要扫描图像并查看每个像素的 3x3 窗口中的值是否与特定模式匹配。我使用以下代码

import numpy as np
import cv2

im = cv2.imread("image.png")
h, w = im.shape[:2]

for i in range(1, h-1):
    for j in range(1, w-1):
        p2 = im[i-1, j]
        p3 = im[i-1, j+1]
        p4 = im[i, j+1]
        p5 = im[i+1, j+1]
        p6 = im[i+1, j]
        p7 = im[i+1, j-1]
        p8 = im[i, j-1]
        p9 = im[i-1, j-1]
        # code for checking the pattern looks something like this:
        if (p2 + p3 + p9) == 1 and p4 == 0 and p5 == 1:
            val = True

但是上面的代码需要很长时间才能完成。我是 Python 和 numpy 的新手，如何有效地扫描 2d numpy 数组？

实际上，我正在尝试将此 thinning code 从 C++ 移植到 Python。

【问题讨论】：

@MrE 我没有对图像进行卷积。每个 3x3 窗口扫描的结果都是一个布尔值。
你如何测试模式？简单的相等，还是更复杂的？
@WarrenWeckesser 查看我编辑的代码。
这称为模板匹配，大多数科学图书馆都有一些现成的例程。
看看这里：stackoverflow.com/a/13240281/401828。如果您需要一个布尔值，则必须对图像应用一些阈值，但这样做应该很容易，因为完美匹配比近似匹配提供更高的峰值。

标签： python opencv numpy

【解决方案1】：

我最终用scipy.weave 编写了用于迭代 Numpy 数组的内联 C++ 代码。它使代码运行得非常快。以前幼稚的方法需要 134 秒才能完成 300x150 图像的处理。虽然这种方法只需要 75 毫秒。

如果您有兴趣，这里是完整的 Python 细化代码：

# Code for thinning a binary image using Zhang-Suen algorithm
from scipy import weave
import numpy as np
import cv2
import sys

def _thinningIteration(im, iter):
    I, M = im, np.zeros(im.shape, np.uint8)
    expr = """
    for (int i = 1; i < NI[0]-1; i++) {
        for (int j = 1; j < NI[1]-1; j++) {
            int p2 = I2(i-1, j);
            int p3 = I2(i-1, j+1);
            int p4 = I2(i, j+1);
            int p5 = I2(i+1, j+1);
            int p6 = I2(i+1, j);
            int p7 = I2(i+1, j-1);
            int p8 = I2(i, j-1);
            int p9 = I2(i-1, j-1);

            int A  = (p2 == 0 && p3 == 1) + (p3 == 0 && p4 == 1) +
                     (p4 == 0 && p5 == 1) + (p5 == 0 && p6 == 1) +
                     (p6 == 0 && p7 == 1) + (p7 == 0 && p8 == 1) +
                     (p8 == 0 && p9 == 1) + (p9 == 0 && p2 == 1);
            int B  = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;
            int m1 = iter == 0 ? (p2 * p4 * p6) : (p2 * p4 * p8);
            int m2 = iter == 0 ? (p4 * p6 * p8) : (p2 * p6 * p8);

            if (A == 1 && B >= 2 && B <= 6 && m1 == 0 && m2 == 0) {
                M2(i,j) = 1;
            }
        }
    } 
    """
    weave.inline(expr, ["I", "iter", "M"])
    return (I & ~M)


def thinning(src):
    dst = src.copy() / 255
    prev = np.zeros(src.shape[:2], np.uint8)
    diff = None

    while True:
        dst = _thinningIteration(dst, 0)
        dst = _thinningIteration(dst, 1)
        diff = np.absolute(dst - prev)
        prev = dst.copy()
        if np.sum(diff) == 0:
            break

    return dst * 255

if __name__ == "__main__":
    src = cv2.imread("image.png")
    if src == None:
        sys.exit()
    bw = cv2.cvtColor(src, cv2.cv.CV_BGR2GRAY)
    _, bw2 = cv2.threshold(bw, 10, 255, cv2.THRESH_BINARY)
    bw2 = thinning(bw2)
    cv2.imshow("src", bw)
    cv2.imshow("thinning", bw2)
    cv2.waitKey()

示例源图像和细化结果：

实用教程：Python Numpy Performance

【讨论】：

【解决方案2】：

你可以用三个卷积来做到这一点。创建三个模板/掩码数组

对每个数组执行卷积。那么您的结果将由以下方式给出：

output = (convolved_with_first == 1) & (convolved_with_second == 0) & ...

【讨论】：

【解决方案3】：

编辑鉴于您的实际模式搜索，我会选择类似：

from numpy.lib.stride_tricks import as strided
win_img = as_strided(im, shape=(h, w - 3 + 1, 3),
                     strides=im.strides + im.strides[-1:])
cond_1 = np.sum(win_img, axis=-1) == 1
cond_2 = im == 0
cond_3 = im == 1

cond = cond_1[:-2, :] & cond_2[1:-1, 2:] & cond_3[2:, 2:]

现在cond[i, j] 具有以im[i+1, j+1] 为中心的窗口的布尔值，并且在每个方向上比原始图像短两个项目。您可以为整个图像获取一个布尔数组：

cond_im = np.zeros_like(im, dtype=bool)
cond_im[1:-1, 1:-1] = cond

查看阵列的窗口视图：

from numpy.lib.stride_tricks import as strided
win_img = as_strided(im, shape=(h - 3 + 1, w - 3+ 1 , 3, 3),
                     strides=im.strides * 2)

现在win_img[i, j] 是一个(3, 3) 数组，其中包含图像的3x3 窗口的内容，左上角位于i, j。

如果您所追求的模式是形状为(3, 3) 的数组pattern，您可以这样做：

np.where(np.all(np.all(win_img == pattern, axis=-1), axis=-1))

获取两个数组的元组，其中包含与您的模式匹配的窗口左上角的行和列。

这里唯一的问题是，当你执行win_img == pattern 时，会创建一个 9 倍于你的图像大小的数组，如果你的图像非常大，这可能会很麻烦。如果您有记忆问题，请将模式检查分成几个波段并循环遍历它们。超过 10 个波段的 for 循环仍将比您当前在图像的整个宽度和高度上的两个嵌套循环快得多。

【讨论】：

【解决方案4】：

您可以尝试以下方法，其中：

直接从图像中读取测试值，而不是创建临时变量；
首先执行成本较低的测试，因为布尔测试是短路的。

result_array = numpy.zeros((h-2, w-2)).astype(bool)

for i in xrange(1, h-1):
    for j in xrange(1, w-1):
        if (im[i, j+1] == 0 and
            im[i+1, j+1] == 1 and
            im[i-1,j] + im[i-1,j+1] + im[i-1, j-1]):
                result_array[i,j] = True

【讨论】：