How to compare fast the value with the ones of neighbours in a 3D array Question:

I want to compare every item in a 3D array with its first neighborhoods. It's really slow when I have a 500x500x500 (e.g. with only values of 0, 1, 2) ndarray. I post the principle lines here:

import numpy as np

# Create a list to stock all the neighbours' coordinations of the voxel wanted 
def check_neighbor(array, x, y, z):
   #top
   t = array[x, y , z + 1]
   #down
   d = array[x, y , z - 1]
   #left
   l = array[x, y - 1 , z]
   #right
   r = array[x, y + 1 , z]
   #front
   f = array[x - 1, y , z]
   #back
   b = array[x + 1, y , z]
   return [t, d, l, r, f, b]

# Check the voxel with all its neighborhood
def compare_neighbor(array, Value2match, Value2Bmatched):
   for index in np.argwhere(array==Value2match)
      output[index] = 1 if Value2BMatched in checkneighbor(array, index[0],index[1], index[2]) else 0 
   return output

# Main
array = np.random.randint(3, shape=(500, 500, 500))
output = compare_neighbor(array, 1, 2)


This code take me hours only for the n=1 neighbour! Is there an efficient way which can also check the 2, 3... nearest neighbours? Can somebody help me?

Solution 1:

Based on jayprich's answer and n1k31t4's comments, I fused the both function into one and replace the argwhere() by where(). The advantage of this code is that we don't iterate voxel by voxel but do it in a vectorized way:

import numpy as np
# Build a helper function to SHIFT(not roll) a 3Darray
def shift_helper(array, neib_value, shift=0, axis=0):
    #Roll the 3D array along axis with certain unity
    _array = np.roll(array == neib_value, shift=shift, axis=axis)

    # Cancel the last/first slice shifted to the first/last slice
    if axis == 0:
        if shift &gt;= 0:
            _array[:1, :, :] = 0
        else:
            _array[-1:, :, :] = 0
        return _array
    elif axis == 1:
        if shift &gt;= 0:
            _array[:, :1, :] = 0
        else:
            _array[:, -1:, :] = 0
        return _array
    elif axis == 2:
        if shift &gt;= 0:
            _array[:, :, :1] = 0
        else:
            _array[:, :, -1:] = 0
        return _array


def compneib(array, that_value, neib_value):
    _array = np.zeros(array.shape)
    _array[np.where((array == that_value)
                    &amp; (shift_helper(array, neib_value, shift=-1, axis=0)
                    | shift_helper(array, neib_value, shift=1, axis=0)
                    | shift_helper(array, neib_value, shift=-1, axis=1)
                    | shift_helper(array, neib_value, shift=1, axis=1)
                    | shift_helper(array, neib_value, shift=-1, axis=2)
                    | shift_helper(array, neib_value, shift=1, axis=2)
                    ))] = 1
    return _array

# Main
array = np.random.randint(3, shape=(500, 500, 500))
output = compneib(array, 1, 2)


Solution 2:

If there would still be a faster way.
NumPy performs operations in a vectorised manner, you should try and operate on the whole array and avoid explicit loops.

e.g. np.argwhere( (test == 1) &amp; np.roll( test==2 , shift=-1 , axis=0 ) )

can be done for different shift and axis values and results combined