1

The code:

A = NxMxLxD matrix
B = NxMxL matrix
index = NxMxL matrix containing values between 1 and D

for i=1:N
 for j=1:M
  for k=1:L
   B(i,j,k) = A(i,j,k,index(i,j,k));
  end
 end
end

How do I write this in vector form, that is using some built it function? I'm basically selecting with the index each point.

user1581390
  • 1,900
  • 5
  • 25
  • 38

2 Answers2

1

It can be done as follows:

ind = bsxfun(@plus, bsxfun(@plus, ...
      (1:N).', N*(0:M-1)), N*M*permute(0:L-1, [3 1 2])) + N*M*L*(index-1);
B = A(ind);

The trick is to build a linear index (ind variable) that corresponds to the indexing (i,j,k,index(i,j,k)), vectorized. I do it with two calls to bsxfun. The key is to remember that linear indices run along columns first, then along rows, then along third-dim slices.

The standard way would be using ndgrid and then sub2ind (see @rayryeng`s answer), but that requires a little more memory.

Community
  • 1
  • 1
Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • +1 - For using `bsxfun`. This will probably be much faster than using `meshgrid` and `sub2ind`! It also looks like I implemented your alternative approach :). – rayryeng Oct 20 '14 at 20:55
  • 1
    Luis, I see you are calculating the linear indices, but you don't actually use them in the final result. I don't even think you need a `reshape`. You can just do `B = A(ind);`, as you have already computed the linear indices! – rayryeng Oct 20 '14 at 21:01
  • @rayryeng Thanks! I always get confused as to when linear indexing preserves shape and when it doesn't – Luis Mendo Oct 20 '14 at 21:05
  • Not a problem :) I get confused too... just takes practice right!? – rayryeng Oct 20 '14 at 21:25
  • Now that I think about it: it _always_ preserves shape, doesn't it? I mean, `size(A(ind))` always equals `size(ind)`. – Luis Mendo Oct 20 '14 at 22:05
  • When you're computing `sub2ind`, whatever size matrices you present to it will give you the same shape yes. Because you were using `bsxfun`, I actually couldn't tell what the shape was at first until I ran the code on my computer. – rayryeng Oct 20 '14 at 22:07
  • I didn't mean `sub2ind`, but right after that, when you use `sub2ind`'s output to index into `A`: that is, the shape of `A(ind)` is the shape of `ind`. – Luis Mendo Oct 20 '14 at 22:09
  • 1
    Oh, in that case yes always! – rayryeng Oct 20 '14 at 22:10
  • 1
    I have just remembered what I think was confusing me about this: `A(ind)` has the same shape as `ind` when `ind` is _integer_; but when `ind` is _logical_ the shape is not kept (it can't be) – Luis Mendo Oct 21 '14 at 13:38
  • 1
    Yes that's exactly it. It will become a single vector that is as long as how many elements are True in the logical array! – rayryeng Oct 21 '14 at 13:45
  • @rayryeng I knew there was some caveat! (I just couldn't remember what it was). Shape is actually [not always preserved](http://blogs.mathworks.com/loren/2006/08/09/essence-of-indexing/): "Indexing with one array C = A(B) produces output the size of B unless both A and B are vectors. _When both A and B are vectors, the number of elements in C is the number of elements in B and with orientation of A_.". Example: try `B = (1:3).'; A = 1:2;`. Then `B(A)` is a column, when it "should" be a row – Luis Mendo Nov 14 '14 at 19:24
1

OK, what this looks like for your output array of B is that for each location in this output, the corresponding location in A has a D-element vector and you want to choose which element to select from this D-element vector based on the value stored in index, which is of the same size as B.

We can achieve this using a combination of meshgrid and sub2ind. First use meshgrid to generate a grid of 3D co-ordinates, which will be the same size as your matrix B. After, use sub2ind to determine a set of linear indices to access the fourth dimension of A. After this, we simply do a straight assignment using the output of sub2ind, which uses a combination of the output of meshgrid and index.

Therefore, try something like this:

[cols,rows,dim] = meshgrid(1:size(A,2), 1:size(A,1), 1:size(A,3));
ind = sub2ind(size(A), rows, cols, dim, index);
B = A(ind);
rayryeng
  • 102,964
  • 22
  • 184
  • 193
  • +1 for readable code. My two `bsxfun`'s are kind of a mess :-) – Luis Mendo Oct 20 '14 at 20:54
  • It would be even more readable with `ndgrid`. Remember [@chappjc's advice](http://stackoverflow.com/a/22461766/2586922)! – Luis Mendo Oct 20 '14 at 20:58
  • @LuisMendo - Good idea :) I actually didn't know that about `ndgrid`. I found it annoying with `meshgrid` to have to specify the second dimension range first. With `ndgrid` I don't have to do this anymore. Nice find Luis! – rayryeng Oct 20 '14 at 20:59
  • 1
    Once you go `ndgrid`, you never go back to `meshgrid` :-D – Luis Mendo Oct 20 '14 at 21:01