Accessing multidimensional array with pairwise indices

Hi everybody

Probably a stupid question but I can't get my head around it. Suppose A is a matrix with dimensions (20,20,1000). I have two vectors x and y with pairwise indices (lets say 5x1). Is there an easy & fast (vectorized) way to extract A(x(1),y(1),:), A(x(2),y(2),:) ... A(x(5),y(5),:)? The only solution I came up with would be sub2ind for the first two dimensions and repmat to expand a binary mask of the indices to the third dimension but that's probably not ideal performance- and memory-wise?

Any ideas?

Thanks in advance

[James Tursa]

How do you want the extracted data (a bunch of 1000 element "vectors") arranged in the result variable?

[Stefan]

The output variable should be formatted as a matrix with numel(x) x size(A,3), so here a matrix with dimensions 5x1000

Answer 1

[Jon]

I think your for loop is only getting data off the top page of the array. When I run your code, B is [ni 1], but all the others results are [ni 2000]. You need to change the last index in A from '1' to ':'. You can also speed up your loop by initializing B before the loop.

B=zeros(numel(ix),size(A,3));
for i=1:numel(ix)
    B(i,:)=A(ix(i),iy(i),:);
end

Also, I found the error in my snippet. The third line should be

ind2=bsxfun(@plus,ind,0:numel(A1):numel(A)-numel(A1));

With these changes, I get the following results for ni=1000.

forloop =0.074519
arrayfun=0.086363
bsxfun1 =0.040455
bsxfun2 =0.077431

Answer 2

[dpb]

One solution; not sure if it's the fastest or not...returns an array of columns length as number of elements i the index arrays by the depth of the planes in the 3D array. Miniature example...

>> x=rand(3,3,2)  % small data array
x(:,:,1) =
  0.1524    0.9961    0.1067
  0.8258    0.0782    0.9619
  0.5383    0.4427    0.0046
x(:,:,2) =
  0.7749    0.0844    0.8001
  0.8173    0.3998    0.4314
  0.8687    0.2599    0.9106
>> ix=[3 1];iy=[2 2];  % x,y coordinates to select
>> squeeze(cell2mat(arrayfun(@(i,j) x(i,j,:),ix,iy,'uniformoutput',0)))
ans =
  0.4427    0.2599
  0.9961    0.0844
>>

Answer 3

[Jon]

I think sub2ind is the recommended approach in most cases. If you want to avoid repmat, you can do something like this.

>> A=rand(3,3,2)
A(:,:,1) =
      0.35166      0.54972       0.7572
      0.83083      0.91719      0.75373
      0.58526      0.28584      0.38045
A(:,:,2) =
      0.56782       0.5308      0.12991
     0.075854      0.77917      0.56882
      0.05395      0.93401      0.46939
  >> ix=[3 1];iy=[2 2];
% Pull out the top page of the array for sizing
  >> A1=A(:,:,1)
% get the indices of the rows/columns you want
  >> ind=sub2ind(size(A1),ix,iy));
% extend the indices to cover all pages of the 3D array
  >> ind2=bsxfun(@sum,ind',0:numel(A1):numel(A)-numel(A1));
% pull out the elements of the array you want 
  >> B=A(ind2)
B =

      0.28584      0.93401
      0.54972       0.5308

You can even put everything in one line if you want.

Answer 4

[Roger Stafford]

If, as you have indicated, x and y are column vectors, do this:

B = reshape(A(bsxfun(@plus,(x+20*(y-1)).',(20*20)*(0:999).')),1,[]);

Answer 5

[Stefan]

Thanks for all the help and replies. I tested all the suggestions with a bigger matrix and compared the computation time. "Unfortunately" the for-loop beats all the other solutions, for larger indexing vectors (larger ni in the code) this might change. For Jon's example I get an error in the bsxfun statement, but it seems to be similar to the last approach. Here are the times I get:

    0.0019704 % for loop
      0.22081 % arrayfun
     0.018394 % bsxfun

and the code:

ni=120;

A=rand(200,200,2000);
ix=floor(rand(ni,1)*100+1);
iy=floor(rand(ni,1)*100+1);

tic;

%for loop
for i=1:numel(ix)
    B(i,:)=A(ix(i),iy(i),1);
end
t(1)=toc;
tic;

%first suggestion
B1=squeeze(cell2mat(arrayfun(@(i,j) A(i,j,:),ix,iy,'uniformoutput',0)));
t(2)=toc;
tic;

%second suggestion, error in bsxfun statement
% A1=A(:,:,1);
% ind=sub2ind(size(A1),ix,iy);
% ind2=bsxfun(@sum,ind',0:numel(A1):numel(A)-numel(A1));
% B2=A(ind2);
% t(3)=toc;

%third suggestion
B3 = permute(A(bsxfun(@plus,(ix+size(A,1)*(iy-1)).',(size(A,1)*size(A,2))*(0:size(A,3)-1).')),[2 1]);
t(3)=toc;

t'

[Stefan]

just to add, for ni=1000, the third code snippet is a clear winner:

       3.6153 % for loop
      0.81556 % arrayfun
      0.18155 % bsxfun