Numpy: 3D index array with the last axis index stored in the 2D array

Question

I have an ndarray of shape(z,y,x) containing values. I am trying to index this array with another ndarray of shape(y,x) that contains the z-index of the value I'm interested in.

 import numpy as np val_arr = np.arange(27).reshape(3,3,3) z_indices = np.array([[1,0,2], [0,0,1], [2,0,1]]) 

Since my arrays are quite large, I tried using np.take to avoid unnecessary copies of the array, but just can't wrap it around indexing three-dimensional arrays.

How do I need to index val_arr with z_indices to get the values ​​at the desired position along the z axis? Expected Result:

 result_arr = np.array([[9,1,20], [3,4,14], [24,7,17]]) 

Answer 1

You can use choose to make a choice:

 >>> z_indices.choose(val_arr) array([[ 9, 1, 20], [ 3, 4, 14], [24, 7, 17]]) 

The choose function is incredibly useful, but can be a little difficult to understand. In fact, given the array ( val_arr ), we can make a series of choices ( z_indices ) from each n-dimensional slice along the first axis.

In addition: any fancy indexing operation will create a new array, not a representation of the original data. It is not possible to index val_arr with z_indices without creating a new array.

Answer 2

With readability, np.choose definitely looks great.

If performance is essentially, you can calculate linear indexes and then use np.take or use the flattened version with .ravel() and extract these specific elements from val_arr . The implementation will look something like this:

 def linidx_take(val_arr,z_indices): # Get number of columns and rows in values array _,nC,nR = val_arr.shape # Get linear indices and thus extract elements with np.take idx = nC*nR*z_indices + nR*np.arange(nR)[:,None] + np.arange(nC) return np.take(val_arr,idx) # Or val_arr.ravel()[idx] 

Run-time checks and verification of results -

The Ogrid solution from here turns into a generic version for these tests, for example:

 In [182]: def ogrid_based(val_arr,z_indices): ...: v_shp = val_arr.shape ...: y,x = np.ogrid[0:v_shp[1], 0:v_shp[2]] ...: return val_arr[z_indices, y, x] ...: 

Case # 1: Less data

 In [183]: val_arr = np.random.rand(30,30,30) ...: z_indices = np.random.randint(0,30,(30,30)) ...: In [184]: np.allclose(z_indices.choose(val_arr),ogrid_based(val_arr,z_indices)) Out[184]: True In [185]: np.allclose(z_indices.choose(val_arr),linidx_take(val_arr,z_indices)) Out[185]: True In [187]: %timeit z_indices.choose(val_arr) 1000 loops, best of 3: 230 µs per loop In [188]: %timeit ogrid_based(val_arr,z_indices) 10000 loops, best of 3: 54.1 µs per loop In [189]: %timeit linidx_take(val_arr,z_indices) 10000 loops, best of 3: 30.3 µs per loop 

Case # 2: More Data

 In [191]: val_arr = np.random.rand(300,300,300) ...: z_indices = np.random.randint(0,300,(300,300)) ...: In [192]: z_indices.choose(val_arr) # Seems like there is some limitation here with bigger arrays. Traceback (most recent call last): File "<ipython-input-192-10c3bb600361>", line 1, in <module> z_indices.choose(val_arr) ValueError: Need between 2 and (32) array objects (inclusive). In [194]: np.allclose(linidx_take(val_arr,z_indices),ogrid_based(val_arr,z_indices)) Out[194]: True In [195]: %timeit ogrid_based(val_arr,z_indices) 100 loops, best of 3: 3.67 ms per loop In [196]: %timeit linidx_take(val_arr,z_indices) 100 loops, best of 3: 2.04 ms per loop 

Answer 3

Inspired by this thread using np.ogrid :

 y,x = np.ogrid[0:3, 0:3] print [z_indices, y, x] [array([[1, 0, 2], [0, 0, 1], [2, 0, 1]]), array([[0], [1], [2]]), array([[0, 1, 2]])] print val_arr[z_indices, y, x] [[ 9 1 20] [ 3 4 14] [24 7 17]] 

I have to admit that multidimensional fantasy indexing can be messy and confusing :)

Answer 4

If you have numpy> = 1.15.0, you can use numpy.take_along_axis . In your case:

 result_array = numpy.take_along_axis(val_arr, z_indices.reshape((3,3,1)), axis=2) 

This should give you the result you want in one neat line of code. Pay attention to the size of the index array. It should have the same number of dimensions as your val_arr (and the same size in the first two dimensions).