function [fractalDim,Hurst,outImage] = fractalDimension(inImage,epsilon,window)
halfmask = floor(window/2);
window = 2*halfmask + 1;
nmask = double(window)^2;
mask = ones(window,window)./nmask;
log_epsilon = log(double(epsilon));
[nrow,ncol] = size(inImage);
idata2e = zeros(nrow,ncol);
idata2r = zeros(nrow,ncol);
outImage = zeros(nrow,ncol);
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
Best Answer