%function multifractal(A)
A=imread('2.jpg');%读入图像，并转化为矩阵；
format long g
L=length(A);
i=1;
modify=1;
tmin=2; % 边框间距，“※”
tmax=10;
ttmin=-10;
ttmax=10; % 自定义 q 的范围
for r=tmin:1:tmax
c(i,1)=mod(L,r);
i=i+1;
end
c'; % 计算不能被边长r整除的余数
a=L-c'; % 计算并剔除掉不能被边长r整除的原始数据
n=length(a); % 求解格网化边长的个数，即为 n
TT=[];
j=1;
r=tmin; % 自定义项，“※-2”
for i=1:1:n % 即n=25-10+1，自定义的结果
B=A(1:a(i),1);
U=reshape(B,r,length(B)/r);
T=mean(U);
T=T'.*r^3;
TT(1: length(T),i)=[T];
modifying(modify,1)=length(T);
r=r+1;
if r>= tmin+n % 自定义项，“※-3”
break; % 边长超过10+n，超过初始限制，则程序自动终止
end
modify=modify+1;
end
modifying;
TT= nthroot(TT,1);
% 或者不缩小
% TT;
TT=nonzeros(TT);
for cugb=1:1:n
modifying_modifying(cugb,1)=sum(modifying(1:cugb)); % 有影响的新加卷 % & * % ￥ # @ ！） ……
end
modifying_modifying;
j=1;
% q 为任意数，这里取1到n，为n，与 k取值保持一致，q过大，计算机无法
%识别，默认为无穷大，q过小，结果接近0，则意义不明确
for q=ttmin:ttmax %这里取 q=-10:1:10
for k=1:1:n
X=TT(1:modifying_modifying(k,1),1).^q;
if k>1
X=TT(modifying_modifying(k-1,1)+1:modifying_modifying(k,1),1).^q;
end
t=sum(X);
XX(k,j)=[t]; % 这里用到两个循环，即考虑到了幂函数，又需考虑求和
end
j=j+1;
end
XX; % 得到质量分配函数，Xq(ξ)，
X=log(tmin:1:tmax);
% X=log(tmax:-1:tmin); % 此系以前的自定义输入结果，“※—4”
Y=log(XX);
% figure(1)
% plot(X',Y,'o-k') % 至此，计算多重分形谱的第一步，分配函数构建完毕
side_length= tmin:1:tmax; % 自定义网格边长，“※—5”
side_length=side_length';
q=ttmin:ttmax; % q=-5:1:n-10，q=-10:1:10
m=1;
[ha,hb]=size(XX);
for i=1:1:hb
% XX=XX’; % or not
s=XX(:,i); % XX
b=polyfit(log(side_length),log(s),1); % 在对数尺度下计算斜率
slope(m,1)=b(1,1);
m=m+1;
end
slope; % 这里的Slope即为质量指数，τ(q)
N=polyfit(q', slope,1);
plot(q', slope) % 此步是考察τ(q)-q 之间的关系，
a=diff(slope)./diff(q'); % 第三步计算，diff函数求偏导确实少一列
q=q';
f_a=a.*q(1:end-1)-slope(1:end-1);
% figure(2)
% a=sort(a,'ascend');
plot(a,f_a,'o-k');
xlabel('α','FontSize',12);
ylabel('f(α)','FontSize',12);
% polyfit(a,f_a,3)
a=sort(a,'ascend');
% a+1
% f_a+1
%a % 奇异性指数
%f_a % 不同层次的分维数
