All of the trick here is about linear indexing as used in MATLAB. We create offset linear indices keeping in mind the final output's size at the back of mind. To start off, offset indices are to be created for indexing into one 3D slice of A
and then onto the entire 3D array. These are named as offset2D
and offset3D
respectively in the code. Finally, we add the linear indices as obtained from m
and n
.
Assuming you have the m
and n
indices saved into two separate 1D arrays named, say m_arr
and n_arr
, you would have the final vectorized implementation with bsxfun
like so -
%// Say you have the m,n arrays are created like this -
m_arr = round(rand(30,1)*32)+1;
n_arr = round(rand(30,1)*32)+1;
%// Get linear indices with m,n
mn_arr = (n_arr-1)*size(A,1) + m_arr;
%// Calculate offset indices for 2D and then 3D versions
offset2D = bsxfun(@plus,[0:27]',[0:27]*size(A,1)); %//'
offset3D = bsxfun(@plus,offset2D(:),[0:30-1]*numel(A(:,:,1)));
%// Incorporate m,n indices into offset to get final linear indices to index into A
lidx = bsxfun(@plus,mn_arr(:).',offset3D); %//'
%// Initialize output array and then index into A to values from a
A = zeros(60,60,30);
A(lidx) = a;
For future readers, here's parameterized versions of the for-loop and vectorized codes -
%// Parameters
M = 4;
D = 3;
mx = 3;
%// For-loop code
A = zeros(M+mx,M+mx,D);
a = rand(M,M,D);
m_arr = round(rand(D,1)*mx)+1;
n_arr = round(rand(D,1)*mx)+1;
for i=1:D
m = m_arr(i);
n = n_arr(i);
A(m:m+M-1,n:n+M-1,i) = a(:,:,i);
end
%// Vectorized code
mn_arr = (n_arr-1)*size(A,1) + m_arr;
offset2D = bsxfun(@plus,[0:M-1]',[0:M-1]*size(A,1)); %//'
offset3D = bsxfun(@plus,offset2D(:),[0:D-1]*numel(A(:,:,1)));
lidx = bsxfun(@plus,mn_arr(:).',offset3D); %//'
A_vectorized = zeros(M+mx,M+mx,D);
A_vectorized(lidx) = a;
Finally, let's compare the outputs from loopy and vectorized codes -
>> max(abs(A(:)-A_vectorized(:)))
ans =
0