Vectorized Implementation for a Convoluted Numpy Script

Question

I am working on a personal project which involves predicting weather pattern movements from radar. I have three n by m numpy arrays; one with precipitation intensity values, one with the movement (in pixels) in the X direction of that precipitation and one with the movement (in pixels) in the Y direction of that precipitation. I want to use these three arrays to determine the location of the precipitation pixels using the offsets in the other two arrays.

xMax = currentReflectivity.shape[0]
yMax = currentReflectivity.shape[1]                            
for x in xrange(currentReflectivity.shape[0]):
    for y in xrange(currentReflectivity.shape[1]):
        targetPixelX = xOffsetArray[x,y] + x
        targetPixelY = yOffsetArray[x,y] + y
        targetPixelX = int(targetPixelX)
        targetPixelY = int(targetPixelY)
        if targetPixelX < xMax and targetPixelY < yMax:
            interpolatedReflectivity[targetPixelX,targetPixelY] = currentReflectivity[x,y]

I can't think of a way to vectorize this; any ideas?

Thank you!!! This was exactly what I was looking for!! Just marked as an answer. — WXMan
– WXMan, Commented May 17, 2017 at 11:59

Divakar · Accepted Answer · 2017-05-12 20:17:16Z

Here's a vectorized approach making use of broadcasting -

x_arr = np.arange(currentReflectivity.shape[0])[:,None]
y_arr = np.arange(currentReflectivity.shape[1])

targetPixelX_arr = (xOffsetArray[x_arr, y_arr] + x_arr).astype(int)
targetPixelY_arr = (yOffsetArray[x_arr, y_arr] + y_arr).astype(int)

valid_mask = (targetPixelX_arr < xMax) & (targetPixelY_arr < yMax)

R = targetPixelX_arr[valid_mask]
C = targetPixelY_arr[valid_mask]
interpolatedReflectivity[R,C] = currentReflectivity[valid_mask]

Runtime test

Approaches -

def org_app(currentReflectivity, xOffsetArray, yOffsetArray):
    m,n = currentReflectivity.shape
    interpolatedReflectivity = np.zeros((m,n))

    xMax = currentReflectivity.shape[0]
    yMax = currentReflectivity.shape[1]                            
    for x in xrange(currentReflectivity.shape[0]):
        for y in xrange(currentReflectivity.shape[1]):
            targetPixelX = xOffsetArray[x,y] + x
            targetPixelY = yOffsetArray[x,y] + y
            targetPixelX = int(targetPixelX)
            targetPixelY = int(targetPixelY)

            if targetPixelX < xMax and targetPixelY < yMax:
                interpolatedReflectivity[targetPixelX,targetPixelY] = \
                currentReflectivity[x,y]
    return interpolatedReflectivity

def broadcasting_app(currentReflectivity, xOffsetArray, yOffsetArray):
    m,n = currentReflectivity.shape
    interpolatedReflectivity = np.zeros((m,n))

    xMax, yMax = m,n        
    x_arr = np.arange(currentReflectivity.shape[0])[:,None]
    y_arr = np.arange(currentReflectivity.shape[1])

    targetPixelX_arr = (xOffsetArray[x_arr, y_arr] + x_arr).astype(int)
    targetPixelY_arr = (yOffsetArray[x_arr, y_arr] + y_arr).astype(int)

    valid_mask = (targetPixelX_arr < xMax) & (targetPixelY_arr < yMax)     
    R = targetPixelX_arr[valid_mask]
    C = targetPixelY_arr[valid_mask]
    interpolatedReflectivity[R,C] = currentReflectivity[valid_mask]
    return interpolatedReflectivity

Timings and verification -

In [276]: # Setup inputs
     ...: m,n = 100,110  # currentReflectivity.shape
     ...: max_r = 120  # xOffsetArray's extent
     ...: max_c = 130  # yOffsetArray's extent
     ...: 
     ...: currentReflectivity = np.random.rand(m, n)
     ...: xOffsetArray = np.random.randint(0,max_r,(m, n))
     ...: yOffsetArray = np.random.randint(0,max_c,(m, n))
     ...: 

In [277]: out1 = org_app(currentReflectivity, xOffsetArray, yOffsetArray)
     ...: out2 = broadcasting_app(currentReflectivity, xOffsetArray, yOffsetArray)
     ...: print np.allclose(out1, out2)
     ...: 
True

In [278]: %timeit org_app(currentReflectivity, xOffsetArray, yOffsetArray)
100 loops, best of 3: 6.86 ms per loop

In [279]: %timeit broadcasting_app(currentReflectivity, xOffsetArray, yOffsetArray)
1000 loops, best of 3: 212 µs per loop

In [280]: 6860.0/212        # Speedup number
Out[280]: 32.35849056603774

Mad Physicist · Accepted Answer · 2017-05-12 19:49:13Z

1

I am pretty sure that you can vectorize this by just taking everything out of the loop:

targetPixelX = (xOffsetArray + np.arange(xMax).reshape(xMax, 1)).astype(np.int)
targetPixelY = (yOffsetArray + np.arange(yMax)).astype(np.int)
mask = ((targetPixelX < xMax) & (targetPixelY < yMax))
interpolatedReflectivity[mask] = currentReflectivity[mask]

This will be much faster but more memory intensive. Basically, targetPixelX and targetPixelY are now arrays containing the values for each pixel that were before computed on a per-iteration basis.

Only the masked values are set in interpolatedReflectivity, similarly to what the if statement was doing in the loop.

answered May 12, 2017 at 19:49

Mad Physicist

116k29 gold badges202 silver badges292 bronze badges

Collectives™ on Stack Overflow

Vectorized Implementation for a Convoluted Numpy Script

2 Answers 2

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related