Calculate the average of each matrix block

Hi,

I have a square matrix as shown below in my code and I would like to ocalculate the averge of each block of it. My code is follwoing:

close all;
clear all;

A = [1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ];

cell_ndp = 2;
cell_nds = 2;

ObsPoints = cell_ndp*cell_nds; %number of total sp observables in each sub
subs = size(A,1)/ObsPoints; %number of main subdivisions
G_sam = zeros(subs); %the averaged sampled G matrix
% calculate the averages
for icell = 1:subs
    obs_hor_start = (ObsPoints * (icell - 1)) + 1;
    obs_hor_end = ObsPoints * icell;
    for obs_hor = obs_hor_start:obs_hor_end
        obs_ver_start = (ObsPoints * (icell - 1)) + 1;
        obs_ver_end = ObsPoints * icell;
        for obs_ver = obs_ver_start:obs_ver_end
            G_sam(icell) = G_sam(icell) + A(obs_hor,obs_ver);
        end
    end
    G_sam(icell) = G_sam(icell) / (ObsPoints*ObsPoints);
end

figure;
subplot(1,2,1)
spy(A)
subplot(1,2,2)
spy(G_sam)

However, this code calculates the averge of the diagonal blocks and leaving the rest, as shown in the spy figures below.

What I need is something like the following:

where the averages of all blocks are calculated in the same order of appearing.

Any help would be appreicted.

Thanks.

[Stephen23]

This is a simple task, for which MATLAB does not require lots and lots of loops and adding numbers one-at-at-time.

Use arrays to your benefit, not lots of nested loops.

Answer 1

[Vishnu]

Hi Lama Hamadeh,

To calculate the average of each block in the matrix in order, you need to modify the code as follows:

In this modified code, the outer loops iterate over the row and column indices of the blocks, and the inner loops calculate the sum of each block. The average of each block is then calculated and stored in the G_sam matrix. The resulting G_sam matrix will contain the averages of each block in order.

close all;
clear all;

A = [1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ];

cell_ndp = 2;
cell_nds = 2;

ObsPoints = cell_ndp * cell_nds; % number of total sp observables in each sub
subs = size(A, 1) / ObsPoints; % number of main subdivisions
G_sam = zeros(subs); % the averaged sampled G matrix

% calculate the averages
for irow = 1:cell_ndp:subs*cell_ndp
    for icol = 1:cell_nds:subs*cell_nds
        block_sum = 0;
        for obs_hor = irow:irow+cell_ndp-1
            for obs_ver = icol:icol+cell_nds-1
                block_sum = block_sum + A(obs_hor, obs_ver);
            end
        end
        block_avg = block_sum / (cell_ndp * cell_nds);
        block_row = ceil(irow / cell_ndp);
        block_col = ceil(icol / cell_nds);
        G_sam(block_row, block_col) = block_avg;
    end
end

figure;
subplot(1, 2, 1)
spy(A)
subplot(1, 2, 2)
spy(G_sam)

You can also plot only the four slots provided by only taking the subarray of the averages(G_sam).

[LH]

Thanks for your reply.

I am confused though about the G_sam output. Using your modified code, the output average are:

G_sam =

    1.5000    3.5000    1.5000    3.5000
    1.5000    3.5000    1.5000    3.5000
    1.5000    3.5000    1.5000    3.5000
    1.5000    3.5000    1.5000    3.5000

which I think do not represent the averages of each block. I think it should be something like:

    G_sam =

    2.5000    2.5000    2.5000    5.5000
    2.5000    2.5000    2.5000    5.5000
    2.5000    2.5000    2.5000    5.5000
    2.5000    2.5000    2.5000    5.5000

[LH]

I think I solved it. It has to be like the following:

ObsPoints = cell_ndp * cell_nds; % number of total sp observables in each sub
subs = size(A, 1) / ObsPoints; % number of main subdivisions
G_sam = zeros(subs); % the averaged sampled G matrix
% calculate the averages
for irow = 1:subs:subs*ObsPoints
    for icol = 1:subs:subs*ObsPoints
        block_sum = 0;
        for obs_hor = irow:irow+subs-1
            for obs_ver = icol:icol+subs-1
                block_sum = block_sum + A(obs_hor, obs_ver);
            end
        end
        block_avg = block_sum / (ObsPoints * ObsPoints);
        block_row = ceil(irow / ObsPoints);
        block_col = ceil(icol / ObsPoints);
        G_sam(block_row, block_col) = block_avg;
    end
end

Thanks for the help.

Answer 2

[Malay Agarwal]

Hi Lama,

You can use something similar to a convolution operation in Machine Learning. You can start by creating a "filter", which is a 4x4 matrix of ones. Then, you can slide the filter over the 16x16 matrix in "strides" of 4. When the filter is on a 4x4 block in the 16x16 matrix, just multiply (element-wise) the elements in the block with the elements in the filter, take a sum and then divide by 16 to get the average. This opens up the avenue to doing weighted averages. Here's a code snippet which does this:

% Input image
data = randi([0, 255], 16, 16) % Replace with your actual data

data = 16×16

158 99 62 181 204 107 211 41 214 183 102 180 238 235 161 245 73 18 63 255 104 198 249 167 107 193 103 36 218 44 53 135 51 235 100 46 180 242 94 9 99 124 182 157 179 92 59 66 213 125 169 22 85 21 100 125 190 211 134 4 1 68 42 46 71 159 14 254 154 66 252 243 114 193 199 169 70 27 66 232 136 241 32 166 108 16 34 250 8 237 23 85 52 97 149 250 185 243 204 84 40 57 31 192 67 43 235 26 18 6 100 183 178 157 140 153 212 29 77 42 152 191 63 87 93 201 157 82 123 44 223 125 51 27 33 54 59 26 186 101 138 138 91 211 53 171 154 18 149 216 178 177 171 121 170 29 219 30 120 27

% Filter
filter = ones(4, 4);

% Convolution parameters
stride = 4;

% Calculate output size
output_size = floor((size(data) - size(filter)) / stride) + 1

output_size = 1×2

4 4

% Initialize output feature map
output = zeros(output_size);

% Scale
n = size(filter, 1) * size(filter, 2);

% Perform averaging
for i = 1:stride:size(data, 1)-size(filter, 1)+1
    for j = 1:stride:size(data, 2)-size(filter, 2)+1
        % Extract region of interest
        roi = data(i:i+size(filter, 1)-1, j:j+size(filter, 2)-1);
        % Perform element-wise multiplication and sum
        summation = sum(sum(roi .* filter));
        % Then divide
        output((i+stride-1)/stride, (j+stride-1)/stride) = summation / n;
    end
end

output

output = 4×4

116.8750 133.5625 138.6875 117.6250 151.0625 112.6875 118.2500 111.4375 108.7500 113.3750 114.4375 152.6250 113.5625 128.8750 131.1250 120.0000

Answer 3

[Kanishk Singhal]

You can use the mean function for calculating the average of submatrix.

Better way to iterate is using row, col rather than indexing.

It is general code also works for (4,4), (4,2)

close all;
clear all;
A = [1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ];
cell_ndp = 2;
cell_nds = 2;
G_sam = zeros(cell_ndp, cell_nds);
for i=1:cell_ndp:size(A,1)
    for j=1:cell_nds:size(A,2)
        G_sam((i+cell_ndp-1)/cell_ndp, (j+cell_ndp-1)/cell_nds) = mean(A(i:i+cell_ndp-1,j:j+cell_nds-1), "all");
    end
end
figure;
subplot(1,2,1)
spy(A)
subplot(1,2,2)
spy(G_sam)
G_sam

Answer 4

[Stephen23]

Here are some simple MATLAB approaches:

A = repmat([1,2,3,4,1,2,3,4,1,2,3,4,7,6,5,4],16,1)

A = 16×16

1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4

M = mean(permute(reshape(A,4,4,4,4),[2,4,1,3]),3:4)

M = 4×4

2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000

M = blockproc(A,[4,4], @(s)mean(s.data(:))) % requires image toolbox

M = 4×4

2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000

V = 4*ones(1,4);
M = cellfun(@(m)mean(m(:)), mat2cell(A,V,V))

M = 4×4

2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000

Answer 5

[Bruno Luong]

A = [1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ;
    1 2 3 4 1 2 3 4 1 2 3 4 7 6 5 4 ];

M = squeeze(mean(reshape(A,[4 4 4 4]),[1 3]))

M = 4×4

2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000 2.5000 2.5000 2.5000 5.5000

[LH]

Thanks for your answer. May I ask what is the [1 3] inside the reshape function refers to?

[Bruno Luong]

Take the average in first and third dimensions (after reshape).