CODEWITHSUNDEEP
pythonnumpyscipy
PUJA
PUJA

Reputation: 649

Finding the location of maximum sum in a 3D array

I have a time series of 2D NumPy array, I would like to find out the location (center) of a maximum sum with specific window size (n),

I have tried with conv2 like below:

from scipy.signal import convolve2d as conv2

def largest_sum_pos_app1(a, n):
    idx = conv2(a, np.ones((n,n),dtype=int),'same').argmax()
    return np.unravel_index(idx, a.shape)

However, this would provide with the location of the individual 2D array, now I would like to find out the location of the window based on all-time series. Is there a built-in module within numpy or scipy to handle such 3D array.

Just giving as an example : Input array of

       ([[0 1 4 0] 
        [1 2 5 1]
        [2 3 6 0]],

       [[1 2 9 4]
        [2 4 6 2]
        [1 5 1 3]],

       [[0 2 3 1]
        [0 3 5 0]
        [1 4 6 1]])

Taking window of size 3 x 3, the sum of each window would be:

       [[24 22]
        [31 36]
        [24 25]]

Now when we take the overall sum would be [79 83], so I would pick the second window. This was the simple case, but I have a bigger size of the array and thousand of time steps. Is there a way to handle this without any loop.

Upvotes: 1

Views: 131

Answers (2)

David Hoffman
David Hoffman

Reputation: 2343

Your question could use a little clarification, but I'm going to assume that your "windows" are overlapping. In that case you could do:

import numpy as np
import scipy.ndimage as ndi

def largest_sum_pos_app1(a, n):
    # assumes that your data is arranged as (time, y, x)
    # uniform_filter will essentially calculate the sum of all pixels in a
    # neighborhood around each pixel in your original array
    window_sums = ndi.uniform_filter(a, n)
    # to find the index we use argmax, but that requires a little acrobatics
    max_idx = window_sums.reshape((len(a), -1)).argmax(1)
    # the result of unravel_index is designed to be used for NumPy fancy indexing,
    # so we need to reshape it.
    coords = np.array(np.unravel_index(max_idx, a.shape[1:])).T
    # the result is the y, x coordinates for each time point
    return coords

One thing you'll need to pay attention to is the mode parameter of uniform_filter, it determines how the edges of the image are handled. The default is to just pad the edges with zeros, which may be what you want, or maybe not.

Upvotes: 1

Mad Physicist
Mad Physicist

Reputation: 114300

You probably want oaconvolve, which handles multiple dimensions, and allows you to select which ones you want to operate on. Assuming you have an array a of shape (k, width, height), with k being the number of planes:

from scipy.signal import oaconvolve

c = oaconvolve(a, np.ones((1, n, n)), axes=(-2, -1), mode='same')
idx = c.reshape(a.shape[0], -1).argmax(axis=1)
result = np.unravel_index(idx, a.shape[1:])

This does not allow you to select the method by which the convolution will be done, so it may not be the optimal choice of algorithm.

Upvotes: 1

Related Questions