盒子计数法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)))]);