Given a 2D array of size (W,H), for example (3,4). You can think of each point in that 2D array as a corner of a rectangular patch. I want to divide that array up into its individual patches as shown in the exmaple below. The order by which this division happen is row by row and for each patch we start from the top left pixel and move counter clockwise. Meaning patch 1 is [A[0,0],A[1,0],A[1,1],A[0,1]]. Then patch two is [A[0,1],A[1,1],A[1,2],A[0,2]]. The number of patches is (w-1)(h-1) and the final ouput is another 2D array of size ((w-1)(h-1),4) (e.g. (6,4) here), where 6 is the number of patches and 4 are array values for the patch corners in counterclockwise fashion as described above.
Now my question is about the most efficient way to produce these patches in python only using simply python functions and numpy.
Thank you so much.
I can write a double for loop that does the job as below, but I'm wondering if there is a better algorithm that can provide faster performance:
import numpy as np
A=np.array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
H=A.shape[0]
W=A.shape[1]
patches=[]
for y in range(H-1):
for x in range(W-1):
patches.append([A[0+y,0+x],A[1+y,0+x],A[1+y,1+x],A[0+y,1+x]])
patches =np.array(patches)
print(patches)
Assuming you made a typo and meant
for y in range(H-1):
for x in range(W-1):
Then you can get the same result with:
# indeces offsets for a counterclockwise 2x2 block
i = np.array([0,1,1,0])
j = np.array([0,0,1,1])
# top corner indeces for the (H-1)x(W-1) 2D array of blocks
h, w = np.meshgrid(np.arange(H-1), np.arange(W-1), indexing='ij')
# turn to a flat 1D array of blocks
h = h.flatten()
w = w.flatten()
# add the ofsets for every block (create extra axis and broadcast)
h = h[:,None]+i
w = w[:,None]+j
# use advanced indexing to "sample at given indeces
A[h,w]
Comparison (mine is ~3x slower on that toy example, but ~10x faster on bigger arrays):
def yours(A):
H=A.shape[0]
W=A.shape[1]
grids=[]
for y in range(H-1):
for x in range(W-1):
grids.append([A[0+y,0+x],A[1+y,0+x],A[1+y,1+x],A[0+y,1+x]])
grids=np.array(grids)
return grids
def mine(A):
i = np.array([0,1,1,0])
j = np.array([0,0,1,1])
H = A.shape[0]
W = A.shape[1]
h, w = np.meshgrid(np.arange(H-1), np.arange(W-1), indexing='ij')
h = h.flatten()[:,None]+i
w = w.flatten()[:,None]+j
return A[h,w]
A = np.arange(30*20).reshape((30,20))
assert np.all(yours(A) == mine(A))
%timeit yours(A)
# 261 µs ± 229 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
%timeit mine(A)
# 27.4 µs ± 46.4 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)