CODEWITHSUNDEEP
pythonnumpyfiltering
Eshel
Eshel

Reputation: 11

What is the efficient way to perform 2D matrix filtering in python?

I am trying to implement a matrix filtering in Python, and so far the implementation appears to be very slow and inefficient. I wonder if there is an efficient way of performing such filtering.

Provided a large matrix A and a filtering matrix M, the function should return a "remixed" matrix R which is obtained by multiplying each element (i,j) of A by M, then the result is superposed/inserted into R at position (i,j). Please find below the code that is expected to do this.

The example below takes about 68 sec (!) at my computer, which seems very inefficient.

I would be very grateful if you could recommend the way to speed-up this function. Many thanks in advance!


import numpy as np
import time

nx = ny = 1500
n_mix = 50

# matrix to be filtered
A = np.random.random_sample( (nx, ny) )

# filter to be applied to each point:
M = np.random.random_sample( (2*n_mix+1, 2*n_mix+1) )

# the result is stored in "remix":
remix = np.zeros_like(A)

start = time.time()

for i in range(n_mix, nx-n_mix):
    for j in range(n_mix, ny-n_mix):
        remix[i - n_mix:i + n_mix + 1, j - n_mix:j + n_mix + 1 ] += M * A[i,j]

print remix

duration = time.time() - start
print(round(duration))

UPDATE

In fact the ndimage package in scipy has the general convolution function that does the job. I post below the 3 variants of doing the filtering, with respected times. The fastest is with ndimage.convolution (24 seconds vs. 56 and 68 by other methods). However, it still seems rather slow...

import numpy as np
from scipy import ndimage
import time
import sys


def remix_function(A, M):
    n = (np.shape(M)[0]-1)/2
    R = np.zeros_like(A)

    for k in range(-n, n+1):
        for l in range(-n, n+1):
            # Ak  = np.roll(A, -k, axis = 0)
            # Akl = np.roll(Ak, -l, axis = 1)
            R += np.roll(A, (-k,-l), axis = (0,1) ) * M[n-k, n-l]
    return R

if __name__ == '__main__':
    np.set_printoptions(precision=2)

    nx = ny = 1500
    n_mix = 50
    nb = 2*n_mix+1

    # matrix to be filtered
    A = np.random.random_sample( (nx, ny) )
    # filter to be applied to each point:
    M = np.random.random_sample( (nb, nb) )


    # the result is stored in "remix":
    remix1 = np.zeros_like(A)
    remix2 = np.zeros_like(A)
    remix3 = np.zeros_like(A)


#------------------------------------------------------------------------------
# var 1
#------------------------------------------------------------------------------
    start = time.time()

    remix1 = remix_function(A, M)

    duration = time.time() - start
    print('time for var1 =', round(duration))

#------------------------------------------------------------------------------
# var 2
#------------------------------------------------------------------------------
    start = time.time()

    for i in range(n_mix, nx-n_mix):
        for j in range(n_mix, ny-n_mix):
            remix2[i - n_mix:i + n_mix + 1, j - n_mix:j + n_mix + 1 ] += M * A[i,j]


    duration = time.time() - start
    print('time for var2  =', round(duration))

#------------------------------------------------------------------------------
# var 3
#------------------------------------------------------------------------------
    start = time.time()

    remix3 = ndimage.convolve(A, M)

    duration = time.time() - start
    print('time for var3 (convolution) =', round(duration))

Upvotes: 1

Views: 93

Answers (1)

Bugbeeb
Bugbeeb

Reputation: 2161

I can't comment on posts yet, but your double for loop is the problem. Have you tried defining a function and then using np.vectorize?

Upvotes: 1

Related Questions