Vectorizing the assignment of values into a 3D array

Question

Here is the simple code:

A=zeros(60,60,30);
a=rand(28,28,30);

for i=1:30
    m=round(rand*32)+1; %because 60-28=32
    n=round(rand*32)+1;
    A(m:m+27,n:n+27,i)=a(:,:,i);
end

What it does is simply take a 28*28 random matrix and "plant" it inside a bigger (zeroes) matrix, and repeat for 30 times. Each time, it randomly selects a different corner position (i.e. 'm' and 'n') to place the small matrix inside the bigger one. I'm sure I can do this without the for loop - just not sure how.

Answer 1

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