sampleFunctionAvg (a,sampleRate,sizeFull,distDim)

function [ avgA ] = sampleFunctionAvg (a,sampleRate,sizeFull,distDim)
%SAMPLEFUNCTION 

% Figuring out the indices to start the and end the averaging 

starts = [1, 1];
ends = sizeFull - sampleRate + 1;

counter = codcolon(1,max(sizeFull));
counter1 = localPart(counter);

% Moding by sample-rate, only interested in these points
counter2 = counter1(mod((counter1 - 1),sampleRate)==0);

% Don't want to sample beyond the last point, even if it is multiple of
% sampleRate

if labindex == numlabs && counter2(end) + sampleRate - 1 > sizeFull(distDim)
    counter2 = counter2(1:end-1);
end

% Adjusting our counting by the offset of the overlap nodes
starts(distDim) = counter2(1) - counter1(1) + 1 + sampleRate;
ends(distDim) = counter2(end) - counter1(1) + 1 + sampleRate;

% Starting at starts, Ending at ends, Counting by sampleRate
itx = (starts(1):sampleRate:ends(1));
ity = (starts(2):sampleRate:ends(2));

avgA = zeros(length(itx),length(ity));

% Simple nested for loop to iterate over positions and average
% Note: this might run faster if more vectorized, for this example, the speed
% is not an issue

for nx=1:length(itx)
    for ny = 1:length(ity)
        
        idown = itx(nx); 
        jdown = ity(ny);
        
        iup = idown + (sampleRate - 1); % forward averaging 
        jup = jdown + (sampleRate - 1); % forward averaging
        
        avgA(nx,ny) =...
            sum(sum(a(idown:iup,jdown:jup)))/(sampleRate*sampleRate);
        
    end
end