CODEWITHSUNDEEP
pythonnumpysubmatrix
James Jenkinson
James Jenkinson

Reputation: 1651

Numpy submatrix operations

I'm looking for an efficient way to perform submatrix operations over a larger matrix without resorting to for loops.

I'm currently doing the operation (for a 3x3 window):

newMatrix = numpy.zeros([numRows, numCols])
for i in range(1, numRows-1):
    for j in range(1, numCols-1):
        sub = matrix[i-1:i+2, j-1:j+2]
        newMatrix[i][j] = ... #do things with sub matrix

This is considerably slower than normal operations with numpy matrices. Is there anything numpy has to offer to solve this, or is that hoping for too much?

Edit: Specific example

xWeight = numpy.array([[-1./8, 0, 1./8], [-2./8, 0, 2./8], [-1./8, 0, 1./8]])
yWeight = numpy.array([[1./8, 2./8, 1./8], [0, 0, 0], [-1./8, -2./8, -1./8]])

Inside loop:

        dz_dx = numpy.sum(xWeight * sub)
        dz_dy = numpy.sum(yWeight * sub)

Upvotes: 3

Views: 4427

Answers (4)

Viktor Kerkez
Viktor Kerkez

Reputation: 46596

It seems to me that you're trying to do a simple convolution?

def do(m):
    rows, cols = m.shape
    newMatrix = np.zeros_like(m)
    for i in range(1, rows-1):
        for j in range(1, cols-1):
            sub = matrix[i-1:i+2, j-1:j+2]
            newMatrix[i][j] = numpy.sum(xWeight * sub)
    return newMatrix[1:-1, 1:-1]
>>> res1 = do(matrix)
>>> res2 = scipy.signal.convolve2d(matrix, xWeight)[2:-2,2:-2]
>>> np.allclose(np.abs(res1), np.abs(res2))
True

Didn't went into details about the sign, but that should hopefully put you on the right track.

Upvotes: 3

James Jenkinson
James Jenkinson

Reputation: 1651

I found a solution in numpy.lib.stride_tricks

from numpy.lib.stride_tricks import as_strided

In the method:

    expansion = stride.as_strided(matrix, shape = (numRows-2, numCols-2, 3, 3), strides = matrix.strides * 2)
    xWeight = numpy.array([[-1./8, 0, 1./8], [-2./8, 0, 2./8], [-1./8, 0, 1./8]])
    yWeight = numpy.array([[1./8, 2./8, 1./8], [0, 0, 0], [-1./8, -2./8, -1./8]])

    dx = xWeight * expansion
    dy = yWeight * expansion

    dx = numpy.sum(numpy.sum(dx, axis=3), axis=2)
    dy = numpy.sum(numpy.sum(dy, axis=3), axis=2)

There may well be a better solution, but this is sufficiently simple and general purpose for what I was after. This went through a 1600x1200 matrix in 3.41 seconds, vs 188.47 seconds using for loops.

(Feel free to offer said better solution, if you have it)

Upvotes: 3

Saullo G. P. Castro
Saullo G. P. Castro

Reputation: 58895

It seems you can use np.ix_, see this example from the documentation:

a = np.arange(10).reshape(2, 5)
#array([[0, 1, 2, 3, 4],
#       [5, 6, 7, 8, 9]])

ixgrid = np.ix_([0,1], [2,4])

a[ixgrid]
#array([[2, 4],
#       [7, 9]])

Upvotes: 1

Jiminion
Jiminion

Reputation: 5168

Use scipy instead for image processing operations:

http://scipy-lectures.github.io/advanced/image_processing/

Upvotes: 2

Related Questions