Taking the max of contiguous matrix chunks in MATLAB

Question

Given the matrix:

a =
   1   1   2   2
   1   1   2   2
   3   3   4   4
   3   3   4   4

I would like to get the following four 2x2 matrices:

a1 =
   1   1
   1   1

a2 =
   2   2
   2   2

a3 =
   3   3
   3   3

a4 =
   4   4
   4   4

From there, I would like to take the max of each matrix and then reshape the result into a 2x2 result matrix, like so:

r =
   1   2
   3   4

The location of the result max values relative to their original position in the initial matrix is important.

Currently, I'm using the following code to accomplish this:

w = 2
S = zeros(size(A, 1)/w);
for i = 1:size(S)
  for j = 1:size(S)
    Window = A(i*w-1:i*w, j*w-1:j*w);
    S(i, j) = max(max(Window));
  end
end

This works but it seems like there must be a way that doesn't involve iteration (vectorization).

I tried using reshape like so: reshape(max(max(reshape(A, w, w, []))), w, w, []) however that takes the max of the wrong values and returns:

ans =
   3   4
   3   4

Is there any way to accomplish this without iteration or otherwise improve my iterative method?

Answer 1

UPDATE: I'm not sure how I've ended up with the most votes (as of 2012-10-28). For anyone reading this, please see angainor's or Rody's answers for better solutions that don't require any additional toolboxes.

Here is a horse race of every answer thus far (excluding Nates - sorry, don't have the requisite toolbox):

Z = 1000;

A = [1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4];
w = 2;

%Method 1 (OP method)
tic
for z = 1:Z
S = zeros(size(A, 1)/w);
for i = 1:size(S)
  for j = 1:size(S)
    Window = A(i*w-1:i*w, j*w-1:j*w);
    S(i, j) = max(max(Window));
  end
end
end
toc

%Method 2 (My double loop with improved indexing)
tic
for z = 1:Z
wm = w - 1;
Soln2 = NaN(w, w);
for m = 1:w:size(A, 2)
    for n = 1:w:size(A, 1)
        Soln2((m+1)/2, (n+1)/2) = max(max(A(n:n+wm, m:m+wm)));
    end
end
Soln2 = Soln2';
end
toc


%Method 3 (My one line method)
tic
for z = 1:Z
Soln = cell2mat(cellfun(@max, cellfun(@max, mat2cell(A, [w w], [w w]), 'UniformOutput', false), 'UniformOutput', false));
end
toc

%Method 4 (Rody's method)
tic
for z = 1:Z
b = [A(1:2,:) A(3:4,:)];
reshape(max(reshape(b, 4,[])), 2,2);
end
toc

The results of the speed test (the loop over z) are:

Elapsed time is 0.042246 seconds.
Elapsed time is 0.019071 seconds.
Elapsed time is 0.165239 seconds.
Elapsed time is 0.011743 seconds.

Drat! It appears that Rody (+1) is the winner. :-)

UPDATE: New entrant to the race angainor (+1) takes the lead!

Answer 2

Not very general, but it works for a:

b = [a(1:2,:) a(3:4,:)];
reshape(max(reshape(b, 4,[])), 2,2).'

The general version of this is a bit *ahum* fuglier:

% window size
W = [2 2];

% number of blocks (rows, cols)
nW = size(a)./W;


% indices to first block
ids = bsxfun(@plus, (1:W(1)).', (0:W(2)-1)*size(a,1));

% indices to all blocks in first block-column
ids = bsxfun(@plus, ids(:), (0:nW(1)-1)*W(1));

% indices to all blocks
ids = reshape(bsxfun(@plus, ids(:), 0:nW(1)*prod(W):numel(a)-1), size(ids,1),[]);

% maxima
M = reshape(max(a(ids)), nW)

It can be done a bit more elegantly:

b = kron(reshape(1:prod(nW), nW), ones(W));    
C = arrayfun(@(x) find(b==x), 1:prod(nW), 'uni', false);    
M = reshape(max(a([C{:}])), nW)

but I doubt that's gonna be faster...