Linear indexing in 3D matrix

Hi all,

I am trying to optimize some code to avoid using for loops.

The goal is to extract the values of a 3d matrix M for a list of coordinates (x,y) : so that values = M(coordinates 1:n)

Using linear indexing, I can extract all values for the coordinates (x,y) for 1 slice using sub2ind(y,x,m).

This looks something like this:

linearnanIndices = sub2ind(size(M),y,x,m*ones(length(x)));

I would like to avoid doing it slice by slice and instead extract all values for a given point, up to slice n.

My attempt:

linearnanIndices = sub2ind(size(M),y,x,linspace(1,n,n));

This does not work because the vector do not match.

Thanks in advance!

[James Tursa]

In what order do you want them extracted? All of the x-y values for slice 1, followed by all the x-y values for slice 2, etc.? Or the x(1)-y(1) values for slices 1:n followed by the x(2)-y(2) values for slices 1:n, etc.? I.e., do you want the result in memory order, or in x-y pair order?

[Vincent den Ronden]

I think the second. So that when indexing M(linearindices) returns an matrix with a row for every point and the columns should be the slice 1:n.

Answer 1

[Jan]

Look into the code of sub2ind. It is easy to create an inline version of what you want:

s = [2, 3, 4];
x = reshape(1:prod(s), s);

index = (1:s(1)).' + ...        % Or: reshape(1:s(1), s(1), 1)
        s(1) * (0:s(2)-1) + ... % Or: reshape(0:s(2)-1, 1, s(2))
        s(1)*s(2) * reshape(0:s(3)-1, 1, 1, s(3));

Now x is equal to x(index).

You can change the order of the indices by changing the order of terms in the formula considering the shapes: The i.th term must be a vector in the i.th dimension.

[Vincent den Ronden]

Thank you for the response.

I have quickly tested it and see that index is also a s(1) * s(2) * s(3) matrix with 1 through 200. Meaning index(end,end,end) = 200

I dont understand however how i can use this to extract the slice M(row,col,1:n) into a list. Or ideally,

extract multiple points so that X(# of points,1:n) = M(rows,cols,1:n)

[Jan]

s = [3, 4, 5];
x = reshape(1:prod(s), s) + 0.5;  % Arbitrary test data

rows  = [1,2];
cols  = [2,3];
pages = [3,4];

index = (rows).' + ...
        s(1) * (cols - 1) + ...
        s(1)*s(2) * reshape(pages - 1, 1, 1, numel(pages));
x(index)