1

Given a matrix A, and a list of row indices, and a list of column indices, how to efficiently extract the squared submatrices with size k centered by the row and column indices?

For example:

A = array([[12,  6, 14,  8,  4,  1],
       [18, 13,  8, 10,  9, 19],
       [ 8, 15,  6,  5,  6, 18],
       [ 3,  0,  2, 14, 13, 12],
       [ 4,  4,  5, 19,  0, 14],
       [16,  8,  7,  7, 11,  0],
       [ 3, 11,  2, 19, 11,  5],
       [ 4,  2,  1,  9, 12, 12]])
r = np.array([2, 5])
c = np.array([3, 2])
k = 3

The output should be A[1:4, 2:5] and A[4:7, 1:4]. So basically, the outputs are squared submatrices in size kxk and centered on the [r,c] elements (A[2,3] and A[5,2] in this case)

How to do this efficiently and elegantly? Thanks

tczj
  • 438
  • 4
  • 17
  • I don't think there's any special trick. For each pair of values in `r` and `c` determine the relevant slice (just some basic math), and do the slices. – hpaulj Nov 06 '19 at 03:45
  • Yes, sure, that will do it. But what if the length of r and c is very large, then loop through each case one by one might be slow. – tczj Nov 06 '19 at 04:21
  • The is a way of `viewing` the array as a bunch (possibly overlapping) windows. From those you could select a subset. It uses a `as_strided` function. It's efficient, but a bit hard to comprehend and do right. – hpaulj Nov 06 '19 at 04:49
  • So, would all the submatrices be of the same shape? – Divakar Nov 06 '19 at 05:03

2 Answers2

2

You mean something like this?

for x,y in zip(r,c):
    s = k // 2
    print("position:",[x - s,x + s + 1], [y - s,y + s + 1])
    print(A[x - s:x + s + 1,y - s:y + s + 1])
    print()

Output: 

position: [1, 4] [2, 5]
[[ 8 10  9]
 [ 6  5  6]
 [ 2 14 13]]

position: [4, 7] [1, 4]
[[ 4  5 19]
 [ 8  7  7]
 [11  2 19]]

Note that k should be odd here

ExplodingGayFish
  • 2,807
  • 1
  • 5
  • 14
  • This will do the job, but it's not efficient as it requires a for loop to loop through each element in r and c. It will become slow when the length of r and c is very big. I am wondering if there is a nice slicing way to efficiently get the results even when the number of elements in r and c is big. – tczj Nov 06 '19 at 04:19
2

For the case when the submatrices be of the same shape, we can get sliding windows and then index those with the start indices along the rows and cols for our desired output. To get those windows, we can leverage np.lib.stride_tricks.as_strided based scikit-image's view_as_windows. More info on use of as_strided based view_as_windows -

from skimage.util.shape import view_as_windows

# Get all sliding windows
w = view_as_windows(A,(k,k))

# Select relevant ones for final o/p
out = w[r-k//2,c-k//2]
Divakar
  • 218,885
  • 19
  • 262
  • 358
  • Thanks. This is a very nice solution. Just one more question, since the returned value of `view_as_windows` returns all the slicing windows in the 2D array. This might consume extra memory. Is there a way to get slicing windows only on the given r and c? – tczj Nov 06 '19 at 13:10
  • @Tao Those are sliding windows as views into the input array and hence, no extra memory is used there. So, we are good there. You can read the last third link in the post for more info. – Divakar Nov 06 '19 at 13:18
  • Thanks. But when I print the shape of w. It says (6, 4, 3, 3). The shape of A matrix is (8, 6). So I think all the possible sliding windows are indeed generated. – tczj Nov 08 '19 at 04:54
  • @Tao Use `np.shares_memory(A,w)`. It will print `True`, meaning `w` is a view into `A`. For more info on `np.shares_memory`, you can google for its docs. – Divakar Nov 08 '19 at 04:56
  • I see. Thank you. – tczj Nov 09 '19 at 01:24