Here's one way, but you'll have to change it to have pixels as length dimensions:
%% *fractalDimension*
%% *definition*
function [fractalDim,Hurst,outImage] = fractalDimension(inImage,epsilon,window)
%% *purpose*
% to compute the fractal dimension of an image
%% *example*
%{
MaxLevel = 12; % size of image is 2^MaxLevel+1
seed = 8675309; % seed enables repeatability
H = 0.5; % Hurst parameters a values between 0 and 1
myImage = midpoint(MaxLevel,H,seed);
N = 2.0^MaxLevel;
figure('Color','white');
imagesc(-N/2:N/2,-N/2:N/2,h01,[-3 3]);
title({'Fractional Brownian Motion';['Hurst =' num2str(H) ...
' Fractal Dimension =' num2str(3-H)]},'fontsize',14);
axis equal
axis tight
colormap(bone);
colorbar
set(gca,'fontweight','bold');
epsilon = 11;
window = 21;
[fractalDim,Hurst,outImage] = fractalDimension(myImage,epsilon,window);
%}
%% *inputs*
% inImage - input image
% epsilon - scaling parameter for search (typically between 3 and 11)
% window - dimension of fractal homogeneity (local fractal dimension)
%% *outputs*
% outImage - local fractal dimension
%% *history*
% when who why
% 20200115 mnoah original code
%%
% compute window region and log epsilon
halfmask = floor(window/2);
window = 2*halfmask + 1;
nmask = double(window)^2;
mask = ones(window,window)./nmask;
log_epsilon = log(double(epsilon));
%% allocate space for temporary arrays
[nrow,ncol] = size(inImage);
idata2e = zeros(nrow,ncol);
idata2r = zeros(nrow,ncol);
outImage = zeros(nrow,ncol);
%% create difference arrays
idata2r = ...
abs(inImage - circshift(inImage, 1,1)) + ...
abs(inImage - circshift(inImage,-1,1)) + ...
abs(inImage - circshift(inImage, 1,2)) + ...
abs(inImage - circshift(inImage,-1,2));
idata2e = ...
abs(inImage - circshift(inImage, epsilon,1)) + ...
abs(inImage - circshift(inImage,-epsilon,1)) + ...
abs(inImage - circshift(inImage, epsilon,2)) + ...
abs(inImage - circshift(inImage,-epsilon,2));
datavalr = conv2(idata2r,mask,'same');
datavale = conv2(idata2e,mask,'same');
idx = (datavalr > 0.0 & datavale > 0.0);
outImage(idx) = log(datavalr(idx)) - log(datavale(idx));
outImage(idx) = 3.0 + outImage(idx)/log_epsilon;
fractalDim = mean(mean(outImage(halfmask:end-halfmask,halfmask:end-halfmask)));
Hurst = 3 - fractalDim;
disp(['Image average fractal dimension = ' num2str(fractalDim)]);
disp(['Image average Hurst Parameter = ' num2str(Hurst)]);
end
The example uses a function I shared on matlab central
