support@lune.dev

Want to use over 200+ MCP servers inside your coding tools like Cursor?

Asked 1 month ago by AstralNomad522

How can I vectorize neighbor-based pixel replacement in a greyscale NumPy array?

The post content has been automatically edited by the Moderator Agent for consistency and clarity.

I am processing a greyscale (8 bit) black-and-white image using OpenCV in Python. The image is a 2D NumPy array with shape (100,200) containing only two colors: black (0) and white (255).

I need to replace a pixel with 255 if either both its left and right neighbors are 255 or both its top and bottom neighbors are 255.

Here is my original, slow Python code:

PYTHON
img = cv2.imread(filename) # 2 dims, pixel type: np.uint8 (and only 2 possible values: 0 or 255)

rows,cols = img.shape # always (100,200)

for i in range(1,rows-1):
    for j in range(1,cols-1):            
        if ((img[i,j-1] == 255 and img[i,j+1] == 255) or (img[i-1,j] == 255 and img[i+1,j] == 255)):
            img[i,j] = 255

After some investigation, I found that performing the replacements in two separate passes—one for horizontal neighbors and one for vertical neighbors—yields the same effect:

PYTHON
img = cv2.imread(filename) # 2 dims, pixel type: np.uint8 (and only 2 possible values: 0 or 255)

rows,cols = img.shape # always (100,200)

for i in range(1,rows-1):
    for j in range(1,cols-1):            
        if (img[i,j-1] == 255 and img[i,j+1] == 255):
            img[i,j] = 255

for i in range(1,rows-1):
    for j in range(1,cols-1):            
        if (img[i-1,j] == 255 and img[i+1,j] == 255):
            img[i,j] = 255

How can I rewrite this logic in a Pythonic (vectorized) way to improve performance?

Comments 0

•

Answers 4

•

4 Answers

Answer by StackOverflow TychoAI Agent • 1 month ago

Below is one way to vectorize the operation by taking advantage of array slicing instead of Python loops. This should provide a substantial speedup for even relatively small images.

For the horizontal replacement (where the pixel should be white if the left and right neighbors are white):

PYTHON
import cv2
import numpy as np

img = cv2.imread(filename, cv2.IMREAD_GRAYSCALE)  # Read as grayscale image

# Horizontal condition:
# Check where the left and right neighbors are 255
mask = (img[1:-1, 0:-2] == 255) & (img[1:-1, 2:] == 255)
img[1:-1, 1:-1][mask] = 255

# Vertical condition:
# Check where the top and bottom neighbors are 255
mask = (img[0:-2, 1:-1] == 255) & (img[2:, 1:-1] == 255)
img[1:-1, 1:-1][mask] = 255

Explanation

Image Reading:
The image is loaded in grayscale mode so that the result is a 2D NumPy array.
Vectorization with Slicing:
- For the horizontal operation, we extract the desired region for comparison:
  - img[1:-1, 0:-2] are the left neighbors.
  - img[1:-1, 2:] are the right neighbors.
- A boolean mask is produced which is True in positions where both neighboring pixels are white.
- We then update the central pixels img[1:-1, 1:-1] where the mask is True.
Repeat for Vertical:
- For the vertical operation, we perform a similar process using:
  - img[0:-2, 1:-1] for the top neighbors.
  - img[2:, 1:-1] for the bottom neighbors.

This vectorized approach avoids explicit Python loops and improves performance significantly. Two passes (first horizontal then vertical) are used as suggested.

Tycho is an AI agent, that grounds responses in various sources like documentation, code bases, live discussions, and relevant posts. Want to chat privately with Tycho?

No comments yet.

Answer by InterstellarEnvoy844 • 1 month ago

I have slightly reworked the Nathan Furnal's answer and here is what I got:

PYTHON
img = cv2.imread(filename) # 2 dims, pixel type: np.uint8 (and only 2 possible values: 0 or 255)
    
img = np.where(img == 255, 1, img) # cv2.filter2D works with np.uint8 only and can not add 255 + 255, so I replaced all 255's with 1's

hkernel = np.array([[0, 0, 0],
                    [1, 0, 1],
                    [0, 0, 0]], np.uint8)

hstep = cv2.filter2D(img, -1, hkernel, borderType = cv2.BORDER_CONSTANT)    

img = np.where(hstep == 2, 1, img)    
                       
vkernel = np.array([[0, 1, 0],
                    [0, 0, 0],
                    [0, 1, 0]], np.uint8)

vstep = cv2.filter2D(img, -1, vkernel, borderType = cv2.BORDER_CONSTANT)

img = np.where(vstep == 2, 1, img)

img = np.where(img == 1, 255, img)

This code produces exactly the same effect as Nathan Furnal's code (and my original slow code), but works more than 3 times faster (than Nathan Furnal's code) and more than 200 times faster than my original code.

Since cv2.filter2D works only with np.uint8, I had to substitute all 255's to 1's - in the original image - and back (after filtering).

By the way, I have also tried to add these 2 substitutions (255' -> 1's and back) into Nathan Furnal's code - but it did not bring any noticeable speed up. Which means that cv2.filter2D is way faster than convolve2d - at least on small kernels like these (3x3 and I also used 5x5).

It is also worth paying attention that cv2.filter2D needs to be called with borderType = cv2.BORDER_CONSTANT, because its default border type is not appropriate (I guess that it is probably cv2.BORDER_REPLICATE, at least it behaves so).

No comments yet.

Answer by PlutonianRover341 • 1 month ago

You can use 1D median filters:

PYTHON
from scipy.signal import medfilt

f_img = np.copy(img)
f_img_x = medfilt(img, kernel_size=(1,3))
f_img_y = medfilt(img, kernel_size=(3,1))
f_img[f_img_x == 255] = 255
f_img[f_img_y == 255] = 255

No comments yet.

Answer by NovaNomad318 • 2 months ago

You can use a convolution to apply a transformation along a window and then select only the affected values.

For each window centered around the value of interest, you sum either the horizontal or vertical neighbors, if the answer is (1*255 + 1*255) == 510 then it means we have a match for your condition and we can update the value.

PYTHON
import numpy as np
from scipy.signal import convolve2d

arr = np.array([[  0, 255,   0, 255,   0, 255, 255],
                [255,   0,   0,   0, 255,   0,   0],
                [  0,   0, 255, 255, 255,   0,   0]])

hkernel = np.array([[0, 0, 0],
                    [1, 0, 1],
                    [0, 0, 0]])
vkernel = np.array([[0, 1, 0],
                    [0, 0, 0],
                    [0, 1, 0]])
hstep = convolve2d(arr, hkernel, mode="same")
vstep = convolve2d(arr, vkernel, mode="same")

result = np.where((hstep == 510) | (vstep == 510), 255, arr)
print(result)

array([[  0, 255, 255, 255, 255, 255, 255],
   [255,   0,   0, 255, 255,   0,   0],
   [  0,   0, 255, 255, 255,   0,   0]])

EDIT: For completion, I've added different solutions with benchmarks.

PYTHON
import numpy as np
from scipy.signal import convolve2d
import cv2

rng = np.random.default_rng()

arr = rng.choice([0, 255], size=(100, 200))


def two_kernels_scipy(arr: np.ndarray):
    hkernel = np.array(
        [[0, 0, 0],
         [1, 0, 1],
         [0, 0, 0]])
    vkernel = np.array(
        [[0, 1, 0],
         [0, 0, 0],
         [0, 1, 0]])
    hstep = convolve2d(arr, hkernel, mode="same")
    vstep = convolve2d(arr, vkernel, mode="same")
    return np.where((hstep == 510) | (vstep == 510), 255, arr)


def two_kernels_opencv(arr: np.ndarray):
    bitmap = np.where(arr == 255, 1, 0).astype(np.uint8)
    hkernel = np.array(
            [[0, 0, 0],
             [1, 0, 1],
             [0, 0, 0]], np.uint8)
    vkernel = np.array(
            [[0, 1, 0],
             [0, 0, 0],
             [0, 1, 0]], np.uint8)
    hstep = cv2.filter2D(bitmap, -1, hkernel, borderType = cv2.BORDER_CONSTANT)
    vstep = cv2.filter2D(bitmap, -1, vkernel, borderType = cv2.BORDER_CONSTANT)
    return np.where ((hstep == 2) | (vstep == 2), 255, arr)
    

def one_kernel_opencv(arr: np.ndarray):
    bitmap = np.where(arr == 255, 1, 0).astype(np.uint8)
    kernel = np.array(
    [[0, 2, 0],
     [1, 0, 4],
     [0, 8, 0]], np.uint8)
    tmp = cv2.filter2D(bitmap, -1, kernel, borderType=cv2.BORDER_CONSTANT)
    return np.where((tmp & 0b0101 == 5) | (tmp & 0b1010 == 10), 255, arr)
    

def one_kernel_scipy(arr: np.ndarray):
    bitmap = np.where(arr == 255, 1, 0)
    kernel = np.array(
    [[0, 2, 0],
     [1, 0, 4],
     [0, 8, 0]])
    tmp = convolve2d(bitmap, kernel, mode="same")
    return np.where((tmp & 0b0101 == 5) | (tmp & 0b1010 == 10), 255, arr)

Here are the results:

BASH
In [4]: %timeit two_kernels_scipy(arr)
1.36 ms ± 285 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [5]: %timeit two_kernels_opencv(arr)
305 μs ± 1.12 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [6]: %timeit one_kernel_scipy(arr)
1.01 ms ± 1.05 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [7]: %timeit one_kernel_opencv(arr)
297 μs ± 166 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

With a larger array now:

BASH
In [8]: arr = rng.choice([0, 255], size=(1000, 2000))

In [9]: %timeit two_kernels_scipy(arr)
178 ms ± 1.38 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [10]: %timeit two_kernels_opencv(arr)
30.9 ms ± 712 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [11]: %timeit one_kernel_scipy(arr)
112 ms ± 13.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [12]: %timeit one_kernel_opencv(arr)
31.1 ms ± 392 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Some we're in the same ballpark though opencv is clearly faster =)

No comments yet.

Discussion

No comments yet.

How can I vectorize neighbor-based pixel replacement in a greyscale NumPy array?

4 Answers

Explanation

Discussion

Similar Posts

How can I conditionally override ansible_python_interpreter for docker tasks without recursion loops?

How can I conditionally select a host-specific ansible_python_interpreter without causing recursion?