How to effectively scan 2d numpy array?

Question

I need to scan an image and see if the values in the 3x3 window of each pixel match with a certain pattern. I use the following code

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

But the code above takes forever to finish. I'm new to Python and numpy, how to effectively scan 2d numpy array?

Actually I'm trying to port this thinning code from C++ to Python.

Answer 1

I ended up with scipy.weave to write inline C++ code for iterating the Numpy array. It makes the code runs very fast. The previously naive approach took 134 seconds to finish processing 300x150 image. While this approach takes only 75ms.

Here is the complete thinning code in Python if you interested:

# 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()

Sample source image and the thinning result:

A useful tutorial: Python Numpy Performance

Answer 2

You could do this with three convolutions. Create the three template/mask arrays

1/3 0 0
1/3 0 0
1/3 0 0

0 0 0
0 0 1
0 0 0

0 0 0
0 0 0
0 0 1

perform a convolution with each array. Then your result will be given by:

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

Answer 3

EDIT Given your actual pattern search, I would go with something like:

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:]

Now cond[i, j] has the boolean value for the window centered at im[i+1, j+1], and is two items shorter in each direction than your original image. You could get a boolean array for your whole image as:

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

---

Take a windowed view of your array:

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)

Now win_img[i, j] is a (3, 3) array, with the contents of a 3x3 window of your image with top left corner at i, j.

If the pattern you are after is an array pattern of shape (3, 3), you could simply do:

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

to get a tuple of two arrays, with the rows and columns of the top left corners of windows where your pattern is matched.

Your only problem here is that, when you do win_img == pattern, an array of 9 times the size of your image is created, which may be troublesome if your image is very large. If you have memory problems, divide the pattern checking into several bands and loop through them. A for loop over 10 bands will still be much faster than your current two nested loops over the whole width and height of the image.

Answer 4

You can try the following approach, which:

1. Reads test values directly from the image instead of creating temporary variables;
2. Performs the less expensive tests first, since boolean tests are short-circuited.

.

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