Possible to vectorize matrix split and stack?

Let's say that I have a 2D matrix that I want to divide into sections that have the same numbers of rows and columns, then stack the slices and store the result in a 3D matrix. I have included below a bit of test code using a for-loop that does this successfully, but is there a way to fully vectorize this splitting operation?

This is a task that is repeated many times in my code and is the only task in the workflow that isn't vectorized. I'm using GPUs and the Parallelization Toolbox for speed, so I was hoping to speed up this operation as well if possible. Thanks in advance.

My test code:

% my test matrix
a= rand(6,6);

% my "coordinates", stored in "b"
% lets assume I know exactly how large each "frame" taken out of my test
% matrix should be, which in this case will be a 4x4 pixel block
%
% first value is beginning row index
% second value is ending row index
% third value is beginning column index
% fourth value is ending row index
%
% each row will contain the coordinates for a "frame",
% so here I have specified 2 the coordinates for 2 frames.
b = [[1,4,1,4];
     [2,5,2,5]];

% number of arrays to extract - same as # rows in b
num = size(b,1);

% initialize the final, stacked array
% size is known, as stated before
c = zeros( 4 ,4,num );

% extraction loop
for frame_num=1:num

      c(:,:,frame_num) = a( b(frame_num,1):b(frame_num,2) , b(frame_num,3):b(frame_num,4) );

end

Answer 1

[Matt J]

 m=4; n=4; %desired block size

 Q=reshape(1:numel(a), size(a));
 t=Q(1:m,1:n)-1;

 idx=sub2ind(size(a), b(:,1), b(:,3)).'; 

 c=reshape( a( bsxfun(@plus,idx,t(:))) , m ,n ,[]);

[Benjamin Pound]

Thank you for your quick replay, Matt.

I wanted to look at timing results more realistic data set:

m=128;
n=128;
a = rand(200,200);

begr=ceil(rand(900,1)*70);
endr=begr+m-1;
begc=ceil(rand(900,1)*70);
endc=begc+n-1;

b=[begr,endr,begc,endc];

I passed m,n,a, and b into the for-loop and vectorized functions (where I actually made them into functions, not the script that I posted before), and then ran both codes through the profiler.

The vectorized code seemed to run two or three times longer than the for-loop code, and it seemed to be because of the "a(idx)" part in the reshape command. Does that seem plausible to you? Thanks for your help.

[Matt J]

Yep, it seems plausible to me for large m,n. Although, I would have expected the bottleneck to come from,

bsxfun(@plus,idx,t(:))

Since you say you have GPUs at your disposal, you should make all these into gpuArrays.

Also, i hope you are not doing

   Q=reshape(1:numel(a), size(a));
  t=Q(1:m,1:n)-1;

within the loop if m,n are fixed. That was meant to be a pre-computation.

[Benjamin Pound]

I did all of that analysis on your previous set of code. When I split up the reshape command into its constituents the the bsxfun statement and the indexing of the matrix 'a' takes about the same amount of time. All this stays true when made into gpuArrays.

But anyway, thank you for your help! I would have never come up with that solution on my own.