Reputation: 778
I wish to obtain a matrix consisting of every 3x3 square within a matrix. The matrix is a #channels by m by n. The end result is a matrix of size #channels x (size1-2)*(size2-2) x 9
I can do this in a non-vector fashion via the following code:
clear
size1 = 10;
size2 = 10;
matrix = 1:size1*size2;
matrix =reshape(matrix,[size1 size2]);
matrix_withdraw = 1:(88*size1*size2);
matrix_withdraw = reshape(matrix_withdraw,[88 size1 size2]);
tic
iter = 1;
for ii = 1:size1-2
for jj = 1:size2-2
locs(ii,jj,:,:) = matrix(ii:ii+2,jj:jj+2);
method1(:,iter,:) = reshape(matrix_withdraw(:,ii:ii+2,jj:jj+2),[88,9]);
iter = iter+1;
end
end
locs = permute(locs,[3 4 1 2]);
This is obviously fairly slow for large sizes of matrices. I'm looking to vectorize this. I managed to obtain a solution which works but it's not very clean
locs2 = ones(size1*(size2-2),1)+(0:size1*(size2-2)-1)';
temp = 1:size1*(size2-2);
temp = mod(temp,size1);
temp = (temp>(size1-2)) | (temp==0);
locs2(temp,:) = [];
locs3 = reshape(locs2,[1 1 (size1-2) (size2-2)]);
locs3(2,1,:,:) = reshape(locs2+1,[1 1 (size1-2) (size2-2)]);
locs3(3,1,:,:) = reshape(locs2+2,[1 1 (size1-2) (size2-2)]);
locs3(:,2,:,:) = locs3(:,1,:,:)+size1;
locs3(:,3,:,:) = locs3(:,2,:,:)+size1;
locs3 = permute(locs3,[1 2 4 3]);
locs3 = vec(locs3);
method2 = matrix_withdraw(:,locs3);
method2 = reshape(method2,[88,9,64]);
method2 = permute(method2,[1 3 2]);
Method1 and method2 are equivalent and render the exact same results. Method2 is also 10x faster than method1. My question is, is there a cleaner way to do this?
Upvotes: 1
Views: 35
Reputation: 112749
If you have the Image Processing Toolbox, you can use im2col
to get each block as a column, and then reshape
as a 4-D array:
blockSize = [3 3];
locs = reshape(im2col(matrix, blockSize), [blockSize size(matrix)-blockSize+1]);
Or you can do it directly, building a 4-D indexing array using implicit singleton expansion. This doesn't require any toolboxes:
m = 3; n = 3; % block size
[M, N] = size(matrix);
ind = (1:m).'+(0:n-1)*M + reshape(0:M-m, 1, 1, []) + reshape((0:N-n)*M, 1, 1, 1, []);
locs = matrix(ind);
Upvotes: 3