Average a subset of a matrix in a loop in matlab

Question

I work with an image that I consider as a matrix.

I want to turn a 800 x 800 matrix (A) into a 400 x 400 matrix (B) where the mean of 4 cells of the A matrix = 1 cell of the B matrix (I know this not a right code line) :

B[1,1] =mean2(A[1,1 + 1,2 + 2,1 + 2,2]) 

and so on for the whole matrix ...

B [1,2]=mean2(A[1,3 + 1,4 + 2,3 + 2,4 ])

I thought to :

1) Reshape the A matrix into a 2 x 320 000 matrix so I get the four cells I need to average next to each other and it is easier to deal with the row number afterwards.

Im4bis=reshape(permute(reshape(Im4,size(Im4,2),2,[]),[2,3,1]),2,[]);

2) Create a cell-array with the 4 cells I need to average (subsetted) and calculate the mean of it. That's where it doesn't work

I{1,160000}=ones,
for k=drange(1:2:319999)
    for n=1:160000
        I{n}=mean2(Im4bis(1:2,k:k+1));
    end
end

I created an empty matrix of 400 x 400 cells (actually a vector of 1 x 160000) and I wanted to fill it with the mean but I get a matrix of 1 x 319 999 cells with one cell out of 2 empty.

Looking for light

My Input Image:

Answer 1

Method 1

Using mat2cell and cellfun

AC = mat2cell(A, repmat(2,size(A,1)/2,1), repmat(2,size(A,2)/2,1));

out = cellfun(@(x) mean(x(:)), AC);

Method 2

using im2col

out = reshape(mean(im2col(A,[2 2],'distinct')),size(A)./2);

Method 3

Using simple for loop

out(size(A,1)/2,size(A,2)/2) = 0;
k = 1;
for i = 1:2:size(A,1)
    l = 1;
    for j = 1:2:size(A,2)
        out(k,l) = mean(mean(A(i:i+1,j:j+1)));
        l = l+1;
    end
    k = k+1;
end

Test on input image:

A = rgb2gray(imread('inputImage.png'));

%// Here, You could use any of the method from any answers 
%// or you could use the best method from the bench-marking tests done by Divakar
out = reshape(mean(im2col(A,[2 2],'distinct')),size(A)./2);  

imshow(uint8(out));  

imwrite(uint8(out),'outputImage.bmp');

Output Image:

Final check by reading the already written image

B = imread('outputImage.bmp');

>> whos B

Name        Size              Bytes  Class    Attributes

B         400x400            160000  uint8              

Answer 2

One approach based on this solution using reshape, sum & squeeze -

sublen = 2; %// subset length
part1 = reshape(sum(reshape(A,sublen,[])),size(A,1)/sublen,sublen,[]);
out = squeeze(sum(part1,2))/sublen^2;

---

Benchmarking

Set #1

Here are the runtime comparisons for the approaches listed so far for a input datasize of 800x 800 -

%// Input
A = rand(800,800);

%// Warm up tic/toc.
for k = 1:50000
    tic(); elapsed = toc();
end

disp('----------------------- With  RESHAPE + SUM + SQUEEZE')
tic
sublen = 2; %// subset length
part1 = reshape(sum(reshape(A,sublen,[])),size(A,1)/sublen,sublen,[]);
out = squeeze(sum(part1,2))/sublen^2;
toc, clear sublen part1 out

disp('----------------------- With  BLOCKPROC')
tic
B = blockproc(A, [2 2], @(x) mean(x.data(:))); %// [m n]
toc, clear B

disp('----------------------- With  PERMUTE + MEAN + RESHAPE')
tic
m = 2;n = 2;
T = permute(reshape(permute(reshape(A, size(A, 1), n, []),...
      [2 1 3]), n, m, []), [2 1 3]);
B = reshape(mean(mean(T,1),2), size(A,1)/m, size(A,2)/m);
toc, clear B T m n

disp('----------------------- With  CONVOLUTION')
tic
m = 2;n = 2;
B = conv2(ones(m,1)/m, ones(1,n)/n ,A,'same');
B = B(m-1:m:end ,n-1:n:end);
toc, clear m n B

disp('----------------------- With  MAT2CELL')
tic
AC = mat2cell(A, repmat(2,size(A,1)/2,1), repmat(2,size(A,2)/2,1));
out = cellfun(@(x) mean(x(:)), AC);
toc

disp('----------------------- With  IM2COL')
tic
out = reshape(mean(im2col(A,[2 2],'distinct')),size(A)./2);
toc

Runtime results -

----------------------- With  RESHAPE + SUM + SQUEEZE
Elapsed time is 0.004702 seconds.
----------------------- With  BLOCKPROC
Elapsed time is 6.039851 seconds.
----------------------- With  PERMUTE + MEAN + RESHAPE
Elapsed time is 0.006015 seconds.
----------------------- With  CONVOLUTION
Elapsed time is 0.002174 seconds.
----------------------- With  MAT2CELL
Elapsed time is 2.362291 seconds.
----------------------- With  IM2COL
Elapsed time is 0.239218 seconds.

To make the runtimes more fair, we can use a number of trials of 1000 on top of the fastest three approaches for the same input datasize of 800 x 800, giving us -

----------------------- With  RESHAPE + SUM + SQUEEZE
Elapsed time is 1.264722 seconds.
----------------------- With  PERMUTE + MEAN + RESHAPE
Elapsed time is 3.986038 seconds.
----------------------- With  CONVOLUTION
Elapsed time is 1.992030 seconds.

Set #2

Here are the runtime comparisons for a larger input datasize of 10000x 10000 for the fastest three approaches -

----------------------- With  RESHAPE + SUM + SQUEEZE
Elapsed time is 0.158483 seconds.
----------------------- With  PERMUTE + MEAN + RESHAPE
Elapsed time is 0.589322 seconds.
----------------------- With  CONVOLUTION
Elapsed time is 0.307836 seconds.

Answer 3

Let A denote your matrix and

m = 2; %// block size: rows
n = 2; %// block size: columns

Method 1

Use blockproc:

B = blockproc(A, [m n], @(x) mean(x.data(:)));

Example:

>> A = magic(6)
A =
    35     1     6    26    19    24
     3    32     7    21    23    25
    31     9     2    22    27    20
     8    28    33    17    10    15
    30     5    34    12    14    16
     4    36    29    13    18    11
>> B = blockproc(A, [m n], @(x) mean(x.data(:)))
B =
   17.7500   15.0000   22.7500
   19.0000   18.5000   18.0000
   18.7500   22.0000   14.7500

Method 2

If you prefer the reshaping way (which is probably faster), use this great answer to organize the matrix into 2x2 blocks tiled along the third dimension, average along the first two dimensions, and reshape the result:

T = permute(reshape(permute(reshape(A, size(A, 1), n, []), [2 1 3]), n, m, []), [2 1 3]);
B = reshape(mean(mean(T,1),2), size(A,1)/m, size(A,2)/n);

Method 3

Apply a 2D convolution (conv2) and then downsample. The convolution computes more entries than are really necessary (hence the downsampling), but on the other hand it can be done separably, which helps speed things up:

B = conv2(ones(m,1)/m, ones(1,n)/n ,A,'same');
B = B(m-1:m:end ,n-1:n:end);