3

The code calculates the minimum value at each row and picks the next minimum by scanning the nearby element on the same and the next row. Instead, I want the code to start with minimum value of the first row and then progress by scanning the nearby elements. I don't want it to calculate the minimum value for each row. The outputs are attached.

import numpy as np
from scipy.ndimage import minimum_filter as mf

Pe = np.random.rand(5,5) 



b = np.zeros((Pe.shape[0], 2))

#the footprint of the window, i.e., we do not consider the value itself or value in the row above
ft = np.asarray([[0, 0, 0],
                 [1, 0, 1],
                 [1, 1, 1]])
#applying scipy's minimum filter
#mode defines what should be considered as values at the edges
#setting the edges to INF
Pe_min = mf(Pe, footprint=ft, mode="constant", cval=np.inf)

#finding rowwise index of minimum value
idx = Pe.argmin(axis=1)

#retrieving minimum values and filtered values
b[:, 0] = np.take_along_axis(Pe, idx[None].T, 1).T[0]
b[:, 1] = np.take_along_axis(Pe_min, idx[None].T, 1).T[0]

print(b)

Present Output:

Desired Output:

  • 1
    You need to create a directed graph (with NetworkX for example) where neighboring cells are connected with each other. A graph should be created in a way that you can move along rows and down to the next row. Then you can use the BFS algorithm or something like that. See this question https://stackoverflow.com/questions/59524874/python-iterate-through-connected-components-in-grayscale-image/59561214#59561214 – Mykola Zotko Feb 07 '22 at 07:25
  • @MykolaZotko Thanks. I was wondering if there's a way to modify the above code so that I can get the desired output. –  Feb 07 '22 at 08:16
  • It looks like a graph problem and in my opinion, you can solve it only with graphs. – Mykola Zotko Feb 07 '22 at 08:20
  • When does the algorithm is supposed to end precisely? Is it when the neighbour value are greater than the current one? Is it also when there is no neighbour for the current cell? – Jérôme Richard Feb 07 '22 at 17:50
  • 1
    @JérômeRichard The algorithm should end when it has reached the last row of the matrix and scanned the left and/or the right elements for the minimum. –  Feb 07 '22 at 17:58

1 Answers1

2

You can solve this using a simple while loop: for a given current location, each step of the loop iterates over the neighborhood to find the smallest value amongst all the valid next locations and then update/write it.

Since this can be pretty inefficient in pure Numpy, you can use Numba so the code can be executed efficiently. Here is the implementation:

import numpy as np
import numba as nb

Pe = np.random.rand(5,5)
# array([[0.58268917, 0.99645225, 0.06229945, 0.5741654 , 0.41407074],
#        [0.4933553 , 0.93253261, 0.1485588 , 0.00133828, 0.09301049],
#        [0.49055436, 0.53794993, 0.81358814, 0.25031136, 0.76174586],
#        [0.69885908, 0.90878292, 0.25387689, 0.25735301, 0.63913838],
#        [0.33781117, 0.99406778, 0.49133067, 0.95026241, 0.14237322]])

@nb.njit('int_[:,:](float64[:,::1])', boundscheck=True)
def minValues(arr):
    n, m = arr.shape
    assert n >= 1 and m >= 2
    res = []
    i, j = 0, np.argmin(arr[0,:])
    res.append((i, j))
    iPrev = jPrev = -1
    while iPrev < n-1:
        cases = [(i, j-1), (i, j+1), (i+1, j-1), (i+1, j), (i+1, j+1)]
        minVal = np.inf
        iMin = jMin = -1
        # Find the best candidate (smalest value)
        for (i2, j2) in cases:
            if i2 == iPrev and j2 == jPrev: # No cycles
                continue
            if i2 < 0 or i2 >= n or j2 < 0 or j2 >= m: # No out-of-bounds
                continue
            if arr[i2, j2] < minVal:
                iMin, jMin = i2, j2
                minVal = arr[i2, j2]
        assert not np.isinf(minVal)
        # Store it and update the values
        res.append((iMin, jMin))
        iPrev, jPrev = i, j
        i, j = iMin, jMin
    return np.array(res)

minValues(Pe)
# array([[0, 2],
#        [1, 3],
#        [1, 4],
#        [2, 3],
#        [3, 2],
#        [3, 3],
#        [4, 4],
#        [4, 3]], dtype=int32)

The algorithm is relatively fast as it succeeds to find the result of a path of length 141_855 in a Pe array of shape (100_000, 1_000) in only 15 ms on my machine (although it can be optimized). The same code using only CPython (ie. without the Numba JIT) takes 591 ms.

Jérôme Richard
  • 41,678
  • 6
  • 29
  • 59
  • How do I get the matrix values? For example, in this format, ```[[0.06229945 0.00133828] [0.00133828 0.09301049]...]``` –  Feb 08 '22 at 05:18
  • 1
    Then you can just replace the first line `res.append((i, j))` with `res.append(arr[i, j])` and the `res.append((iMin, jMin))` line with `res.append(minVal)` as well as the output+`res` type+shape (from `int_[:,:]` to `float[:]`). Then you can transform the output array provided by Numba with Numpy so to get the pairs easily with something like `np.hstack([output[:-1].reshape(-1, 1), output[1:].reshape(-1, 1)])` – Jérôme Richard Feb 08 '22 at 11:56