Octave: Multiple submatrices from a matrix

Question

I have a large matrix from which I would like to gather a collection of submatrices. If my matrix is NxN and the submatrix size is MxM, I want to collect I=(N - M + 1)^2 submatrices. In other words I want one MxM submatrix for each element in the original matrix that can be in the top-left corner of such a matrix.

Here's the code I have:

for y = 1:I
    for x = 1:I
        index = (y - 1) * I + x;
        block_set(index) = big_mat(x:x+M-1, y:y+M-1)
    endfor
 endfor

The output if a) wrong, and b) implying there is something in the big_mat(x:x+M-1, y:y+M-1) expression that can get me what I want without needing the two for loops. Any help would be much appreciated

It appears that you are doing this in Octave, but maybe this MATLAB question will help to give you some ideas: http://stackoverflow.com/questions/2678857/general-method-for-making-sub-arrays-around-a-particular-element. — gnovice, Apr 26 '10 at 20:33

Jonas · Accepted Answer · 2010-04-27T12:53:50.933

There seem to be a few things wrong in your code. Here's how I'd do it if I were to use the double loop:

M = someNumber;
N = size(big_mat,1); %# I assume big_mat is square here

%# you need different variables for maxCornerCoord and nSubMatrices (your I)
%# otherwise, you are going to index outside the image in the loops!
maxCornerCoord = N-M+1;
nSubMatrices = maxCornerCoord^2;

%# if you want a vector of submatrices, you have to use a cell array...
block_set = cell(nSubMatrices,1); 
%# ...or a M-by-M-by-nSubMatrices array...
block_set = zeros(M,M,nSubMatrices);
%# ...or a nSubMatrices-by-M^2 array
block_set = zeros(nSubMatrices,M^2);

for y = 1:maxCornerCoord
    for x = 1:maxCornerCoord
        index = (y - 1) * maxCornerCoord + x; 
        %# use this line if block_set is a cell array
        block_set{index} = big_mat(x:x+M-1, y:y+M-1);
        %# use this line if block_set is a M-by-M-by-nSubMatrices array
        block_set(:,:,index) = big_mat(x:x+M-1, y:y+M-1);
        %# use this line if block_set is a nSubMatrices-by-M^2 array
        block_set(index,:) = reshape(big_mat(x:x+M-1, y:y+M-1),1,M^2);
    endfor
 endfor

EDIT

I just saw that there is an implementation of im2col for Octave. Thus, you can rewrite the double-loop as

%# block_set is a M^2-by-nSubMatrices array
block_set = im2col(big_mat,[M,M],'sliding');

%# if you want, you can reshape the result to a M-by-M-by-nSubMatrices array
block_set = reshape(block_set,M,M,[]);

This is probably faster, and saves lots of digital trees.

Thanks; is there a way to do it without the double-loop? I'm assuming if there is the performance of such would be much better than the looped version... — fbrereto, Apr 26 '10 at 23:42
@fbereto: It turns out you can write this as a 1-liner. See my edit. — Jonas, Apr 27 '10 at 00:30
@Jonas: oh it seems you already thought of IM2COL, I guess I didnt see your edit... Consider adding permute/reshape to turn it into the desired shape: `block_set = reshape(permute(im2col(big_mat,[M M],'sliding'),[1 3 2]),2,2,[]);` — Amro, Apr 27 '10 at 01:09
BTW the proper link to the function documentation is: http://octave.sourceforge.net/image/function/im2col.html — Amro, Apr 27 '10 at 01:12
@Amro: Thanks for the suggestions! I think the call to permute isn't necessary, though - or is this something special about Octave? — Jonas, Apr 27 '10 at 01:34

Davorak · Answer 2 · 2010-04-27T21:29:50.527

With Mathematica: This code makes a matrix where each element is a matrix of MxM with each element in the original matrix in the top-left corner of such a matrix.

The matrix elements along the right and bottom are padded with x.

Partition[big_mat, {M, M}, {1, 1}, {1, 1}, x]

Example: alt text http://img130.imageshack.us/img130/6203/partitionf.png

If you leave off x argument, then then it automatically samples periodically.

score 0 · Answer 3 · answered Apr 26 '10 at 20:39

0

About your output being wrong, it might be because of the assignment. You're trying to assign a matrix to a vector position. Try using

block_set(:,:,index) = big_mat(x:x+M-1, y:y+M-1)

instead.

answered Apr 26 '10 at 20:39

Ofri Raviv

24,375
3
55
55

@Ofri: thanks for the tip. What if I _do_ want a vector of matrices? – fbrereto Apr 26 '10 at 20:48
Technically, this is what you'd get with this code. its a vector of matrices, but the index to get to different matrices is the 3rd index. the 1st and 2nd indexes are the indexes inside the matrix. Another possibility is to use cell arrays, but I think in this case a 3D matrix would be better. – Ofri Raviv Apr 27 '10 at 04:07

Octave: Multiple submatrices from a matrix

3 Answers3