compute_homogeneous_length_scale(I)

function [ResStats] = compute_homogeneous_length_scale(I)

% Particles are white, i.e. value = 1
% Matrix is black, i.e. value = 0
% Modifying

I = double(I)/max(max(double(I)));
if mean(mean(I(:,:))) >= 0.5; I = 1 - I; end

% Rest of code

ImageX = double(I);
ImageY = double(I);
Image_Box = double(I);

ImagePower = log2(length(I(1,:)));
ResStats = zeros(ImagePower,9);

% Subsequent passes
for n = 1:ImagePower

% Accumulate statistics
    ResStats(n,1) = length(Image_Box(:,1));
    ResStats(n,2) = std(double(Image_Box(:)));
    ResStats(n,3) = mean(double(Image_Box(:)));
    ResStats(n,4) = 2^ImagePower/length(Image_Box(:,1));
    ResStats(n,5) = ResStats(n,2)/ResStats(n,3);
    ResStats(n,6) = std(double(ImageX(:)));
    ResStats(n,7) = ResStats(n,6)/ResStats(n,3);
    ResStats(n,8) = std(double(ImageY(:)));
    ResStats(n,9) = ResStats(n,8)/ResStats(n,3);

% Calculate new images
    Image_Box = 0.25*(Image_Box(1:2:end,1:2:end)+Image_Box(1:2:end,2:2:end)...
        +Image_Box(2:2:end,1:2:end)+Image_Box(2:2:end,2:2:end));
    ImageX = 0.5*(ImageX(1:2:end,:)+ImageX(2:2:end,:));
    ImageY = 0.5*(ImageY(:,1:2:end)+ImageY(:,2:2:end));

end