There are probably multiple ways of doing this. It could probably be done with a bunch of reshaping and such, but I did it this way. I'm not sure which would be faster or less confusing.
arr_i1 = [111,112,113,114; ...
121,122,123,124];
arr_i2 = [211,212,213,214; ...
221,222,223,224];
arr_i3 = [311,312,313,314; ...
321,322,323,324];
arr_i = cat(3,arr_i1,arr_i2,arr_i3);
% [211,312,113,214;
% 221,322,123,224];
% calculate all 8 subscripts for each axis
xvec = repelem(1:4,2);
yvec = repmat(1:2,[1 4]);
zvec = repelem([2 3 1 2],2);
% convert to linear indices
idx = sub2ind(size(arr_i),yvec,xvec,zvec);
reshape(arr_i(idx),2,[])
Otherwise, trying to do something directly like
arr_i(:,1:4,[2 3 1 2])
Would give the intersection of the subscript ranges. Such a result would be a much larger set than intended, and the desired result would be its block diagonal.
