盒子计数法matlab程序

%基于盒计数法，编制了二维图像分形维数计算程序
clc
clear
clf
c= imread('text.jpg');
axis square %将当前坐标系图形设置为方形。

横轴及纵轴比例是1：1 figure(1) % figure 1
imshow(c) % show image
figure(2)
%将图像转换为二进制图像，基于阈值，0为白色，1为黑色
bw = im2bw(c, graythresh(c));
[siz_x,siz_y]=size(bw);
imshow(bw)
figure(3) %在灰度图像中寻找边缘
e = edge(double(bw));
imshow(e)
c=bw;
%figure(4)
%c=~bw;
%imshow(c)
%求几维数组
dim = ndims(c);
width = max(size(c));
p = log(width)/log(2);
if p~=round(p) || any(size(c)~=width) % p是否为整数
%将A的元素舍入到大于或等于A的最近整数
p = ceil(p);
width = 2^p; % Let width is the power of 2
%设置标准小盒子
mz = zeros(width, width);
mz(1:size(c,1), 1:size(c,2)) = c;
c = mz;
end
%预先分配大小为r的框数
n=zeros(1,p+1);
r=zeros(1,p+1);
n(p+1)=nnz(c); %非零矩阵元素的数目
r(p+1)=2^(-p);
Nx = size(c,1);
Ny = size(c,2);
for p=(p-1):-1:0 %从p-1开始，增量是-1，直到0
numBlocks = 2^p;
sizeBlocks_x = floor(Nx./numBlocks);
sizeBlocks_y = floor(Ny./numBlocks);
flag = zeros(numBlocks,numBlocks);
for i = 1:numBlocks
for j = 1:numBlocks
xStart = (i-1)*sizeBlocks_x + 1;
xEnd = i*sizeBlocks_x;
yStart = (j-1)*sizeBlocks_y + 1;
yEnd = j*sizeBlocks_y;
block = c(xStart:xEnd, yStart:yEnd);
%如果块的任何部分为真，则标记此
flag(i,j) = any(block(:));
end
end
r(p+1)=2^(-p);
n(p+1)=nnz(flag);
figure(4)
hold on
imshow(flag)
end
figure(5)
x=log(r);
y=log(n);
%最小二乘法曲线拟合
p=polyfit(x,y,1);
% 返回在x处求出的n次多项式的值y。

BestFit = polyval(p,x);
plot(x,y,'o','LineWidth',1)
hold on
plot(x,BestFit,'-','LineWidth',2)
k=-std(abs(BestFit-y))
xlabel('In \it\epsilon','FontSize',18)
ylabel('In\itN(\epsilon)','FontSize',18)
FractalDimension =p(:,1) %分形维数
figure(6)
loglog(r, n,'bo-')
grid on
df = -diff(log(n))./diff(log(r))
disp(['Fractal dimension, Df = ' num2str(mean(df(4:7))) ' +/- ' num2str(std(df(4:7)))]);。