function [S,Savg,x,y,xavg,yavg] =  blockAvg_vectorized(Nx,Ny,Nxavg,Nyavg)
% For example, Nx = Ny = 1000, Nxavg = Nyavg = 25;
  
Lx = 10;
Ly = 10;

%% Compute the function to downsample in meshgrid
x = linspace(1,Lx,Nx);
y = linspace(1,Ly,Ny);
[X,Y] = meshgrid(x,y);
S = (5*cos((X+Y)*2*pi) + 2*sin(2*X.*pi) + 2*cos(2*X*pi))';


%% The averaging
% how many steps you're taking in each direction
xsteps = floor(Nx/Nxavg);
ysteps = floor(Ny/Nyavg);

% create the new x and y axis by selecting every N_avg:th point, and
% shifting them left half the length of the downsampling interval -- i.e. 
% the new samplepoint is the middle of the rectangle. 
% N.B. requires linear intervals.
xavg = x(Nxavg:Nxavg:end) - 0.5*(x(Nxavg)-x(1));
yavg = y(Nxavg:Nxavg:end) - 0.5*(y(Nxavg)-y(1));

% xsteps is floor(Nx/Nxavg), so xsteps*Nxavg = Nx - remainder
newX = (1:xsteps*Nxavg);  % removing excess points 
newY = (1:ysteps*Nyavg);  % not used in averaging

% S(newX,newY) selects from S only the elements that would affect the
% computation; reshape(S,nxa,[]) reshapes the S into nxa by N matrix, where 
% N is determined by size of S. This is guaranteed to work since we
% selected from S the correct amount of entries. 
Savg = sum(reshape(S(newX,newY),Nxavg,[]));
% Savg is now the column sums of the previous reshape - we compute the
% average of each of those, and reshape the matrix into xsteps by Ny matrix 
Savg = reshape((1/Nxavg)* Savg,xsteps,[]);

% repeat the previous, but in y-direction
Savg = sum(reshape(Savg',Nyavg,[]));
Savg = reshape((1/Nyavg)* Savg,ysteps,[])';
