[SciPy-user] Indexing question

[Angus McMorland]

Hi Alex,

2008/12/17 Alexander Borghgraef <alexander.borghgraef.rma at gmail.com>:
> Hi all,
>
> I'm doing some image processing work in scipy, and I have a question
> regarding indexing 3d arrays. My work is on particle filters (for a
> simple example see the cookbook entry I wrote:
> http://www.scipy.org/Cookbook/ParticleFilter), which involves
> evaluating a list of points (coordinates) in an image (2d array of
> shape (height, width)). In my code I represent the list of n points x
> as an array of shape (2, n) and I get the corresponding array of n
> image intensity values by writing:
>
>  features = im[tuple(x)]
>
> Now I would like to extend the code to deal with vector images of
> shape (2, height, width), but I'm not sure on how to do the evaluation
> of coordinates here, which should in this case return a list of
> vectors in the form of an (2, n) array. I tried:
>
>  features = im[:, tuple(x.T)]
>
> but that was obviously naive and doesn't work. An old-fashioned
> solution would be:
>
>  features = []
>  for point in x.T:
>    features.append(im[:, point[0], point[1]])
>  features = array(features)
>
> But clearly I would prefer to work with the indexing features of
> numpy. The following works, but it's a pretty ugly hack, and it's not
> extentable to length m vectors:
>
>  features = vstack( ( im[ tuple( vstack( (zeros(n, int), x)  ) )  ],
> im[ tuple( vstack( (ones(n, int), x)  ) )  ]  ) )
>
> So, can any of you more knowledgeable regarding numpy indexing point
> me towards a more elegant solution?


>
> --
> Alex Borghgraef

There are two tricks that will help you out here:

(1) indexing with a list is treated differently to indexing with an array
(2) indexing across the first consecutive dimensions and extracting
latter ones is easier than extracting earlier dimensions together
based on indices in latter dimensions

So, this works:

import numpy as np

sz_x = 100
sz_y = 120
n_pts = 20
n_vec_dims = 2

# make a fake vector image
im = np.random.random(size=(n_vec_dims, sz_x, sz_y))

# make some fake indices to extract
id_x = np.random.random_integers(size=n_pts, low=0, high=sz_x - 1)
id_y = np.random.random_integers(size=n_pts, low=0, high=sz_y - 1)
ids = np.concatenate((id_x[:,None], id_y[:,None]), axis=1)

# here's the important bit...
pts = im.transpose(1,2,0)[list(ids.T)]

# and a check that it's working
print np.all(pts[0] == im[:,id_x[0], id_y[0]])

I hope that's what you were after,

Angus.
-- 
AJC McMorland
Post-doctoral research fellow
Neurobiology, University of Pittsburgh