making coarse matrix from fine resolution matrix

Hi, I am trying to make a coarse resolution matrix of 3600x1800 from a 8640x4320 matrix by summing up the elements of the matrix.

Hi am trying the following: However this doesnot work with fractions (here, 2.4)

%%%%Z1 is the  original 8640x4320 matrix
abc = blockproc(Z1,[2.4,2.4],@(x)sum(x.data));
abc1 = blockproc(abc,[1,2.4],@(x)sum(x.data));

Answer 1

[Matt J]

A 3rd approach, more memory conserving and faster,.

Z1=randi(100,8640,4320);

u = 5;       %upsampling factor
d = 12;      %downsampling factor
t = d/u;     %reduction factor

[m,n]=size(Z1);
L1=speye(m); L2=speye(round(m/t))/u;
R1=speye(n); R2=speye(round(n/t))/u;

L=repelem(L2,1,d) * repelem(L1,u,1);
R=(repelem(R2,1,d) * repelem(R1,u,1)).';

tic
abc4= L*(Z1*R);
toc

Note as well that the matrices L and R are reusable on other Z1 of the same size that one might wish to downsample later on.

Answer 2

[Matt J]

If you have the Image Processing Toolbox,

abc1=imresize(Z1,[3600,1800])

[Image Analyst]

Or if you want it to look coarse/chunky rather than smooth, use the 'nearest' option.

abc = imresize(Z1, [3600,1800], 'nearest');

Though I'm not sure it would be that noticeable - depends on your image.

[SChow]

The imresize function interpolates between the elements of the matrix, I want them to be summed up,

[Matt J]

No, imresize will do antialiasing (pre-smoothing) of the image before it interpolates, unless you turn it off or except as noted here.

[Image Analyst]

To sum up an image, you can use conv2()

kernel = ones(3,3); % Sum up in a local 3x3 neighborhood.
sumImage = conv2(double(grayImage), kernel);

[SChow]

@Matt J, @Image Analyst

The code that I wrote (embeded in the question) effciently sums up the elements of Z1 to abc1 if the factors in the square bracket were not in fractions. Like the code works in this form (below) to give 360x180 matrix.

However I need a 3600x1800 matrix, that means I need to sum every 2.4x2.4 cells

abc = blockproc (Z1, [24,24], @ (x) sum (x.data));
abc1 = blockproc (abc, [1,24], @ (x) sum (x.data));

[Matt J]

However I need a 3600x1800 matrix, that means I need to sum every 2.4x2.4 cells.

I suspect what you are pursuing might not be different enough from imresize to be worth the trouble. The only difference between what you are doing, and what imresize does is that you are using a rectangular anti-aliasing filter window, while imresize uses some other low pass filter (Gaussian?). Does the difference really matter to you? imresize was written to do resolution reduction in a pretty smart way.

The code that I wrote (embeded in the question) effciently sums up the elements of Z1 to abc1 if the factors in the square bracket were not in fractions.

Your blockproc method is quite inefficient, I'm afraid. Compare to sepblockfun (Download)

Z1=rand(8640,4320);

tic;
abc = blockproc(Z1,[4,1],@(x)sum(x.data));
abc1 = blockproc(abc,[1,4],@(x)sum(x.data));
toc; %Elapsed time is 218.637263 seconds.

tic
 abc2=sepblockfun(Z1,[4,4],'sum');
toc; %Elapsed time is 0.086412 seconds.

[SChow]

Yes I agree, Thanks for pointing to sepblockfun

Now I can get the required matrix in matter of seconds,

scaling = 5;  %%%I scale Z1 by factor of 5 so that the sepblockfun works,
S=repelem(Z1/scaling^2, scaling, scaling);%% where S is 43200x21600
abc2 = sepblockfun (S, [12,12], 'sum' ); %%% abc2 is 3600x1800

[SChow]

I suspect what you are pursuing might not be different enough from imresize to be worth the trouble. The only difference between what you are doing, and what imresize does is that you are using a rectangular anti-aliasing filter window, while imresize uses some other low pass filter (Gaussian?). Does the difference really matter to you? imresize was written to do resolution reduction in a pretty smart way

I tried it in both ways and compared with Z1: The sum of matrix Z1 was 5.8x10^8 . the target matrix intended to have the same sum.

[Z1 contains gridded population]

Using repelem and sepblockfun :

scaling = 5;  %%% I scale Z1 by factor of 5 so that the sepblockfun works,
S = repelem (Z1 / scaling ^ 2, scaling, scaling); %% where S is 43200x21600
abc2 = sepblockfun (S, [12,12], 'sum' ); %%% abc2 is 3600x1800

%%%% SUM of abc2 - 5.8x10^8

Using imresize :

abc3=imresize(Z1,0.4166667)
%%%% SUM of abc3 - 1.0269x10^8

[Matt J]

This should give proper agreement in the sums

abc3=imresize(Z1,1/2.4)*2.4^2;

or

abc3=imresize(Z1,1/2.4);
 abc3=abc3/sum(abc3(:))*sum(Z1(:));

[SChow]

Thanks, they work well!!

Answer 3

[SChow]

scaling = 5;  %%% I scale Z1 by factor of 5 so that sepblockfun works,
%%%%% sepblockfun is developed by MattJ and is available for
%downlaod at https://de.mathworks.com/matlabcentral/fileexchange/48089-separable-block-wise-operations
S = repelem (Z1 / scaling ^ 2, scaling, scaling); %% where S is 43200x21600
abc2 = sepblockfun (S, [12,12], 'sum' ); %%% abc2 is 3600x1800

[Matt J]

I would recommend doing this in horizontal and vertical passes, so that the intermediate matrix S won't be so large,

S = sepblockfun(  repelem (Z1 , scaling, 1)      ,   [12,1], 'sum');

abc2 = sepblockfun(  repelem (S, 1, scaling)      ,   [1,12], 'sum')/scaling^2;

[SChow]

@MattJ

I completely agree!!

Thank you so much for your help

[SChow]

Hi @Matt J,

Is there a way to use sepblockfun if the matrix contains NaN values?

in other words, it does not seem to support functions - 'nansum' / 'nanmean'

[Matt J]

You can do,

S = sepblockfun(  repelem (Z1 , scaling, 1)      ,   [12,1], @nansum);
abc2 = sepblockfun(  repelem (S, 1, scaling)      ,   [1,12], @nansum)/scaling^2;

You cannot apply sepblockfun to nanmean directly, because nanmean is not separable, e.g.

>> A=rand(5); A(1:3)=nan

A =

       NaN    0.9730    0.8253    0.8314    0.4168
       NaN    0.6490    0.0835    0.8034    0.6569
       NaN    0.8003    0.1332    0.0605    0.6280
    0.6596    0.4538    0.1734    0.3993    0.2920
    0.5186    0.4324    0.3909    0.5269    0.4317


>> nanmean(nanmean(A,1),2)

ans =

    0.5163
    
>> nanmean(nanmean(A,2),1)

ans =

    0.5142

However, you can implement block-wise nanmean indirectly by using the separability of nansum,

>> sepblockfun(A,[5,5],@nansum)./sepblockfun(~isnan(A),[5,5],@nansum)

ans =

    0.5063
    
>> nanmean(A(:))

ans =

    0.5063