基于WTMM算法的图像多重分形谱计算matlab仿真

1.算法运行效果图预览

2.算法运行软件版本

matlab2022a

3.算法理论概述

基于WTMM算法的图像多重分形谱计算是一种利用小波变换模极大值（WTMM）方法，对图像进行多重分形分析的方法。下面将详细介绍这种方法的原理和数学公式。

3.1、WTMM算法概述

分形理论是一种研究自然界中不规则、复杂现象的数学工具，而多重分形则是分形理论的一个重要分支，用于描述具有不同奇异程度的分形结构。在图像处理中，多重分形分析可以帮助我们更好地理解图像的纹理、边缘等特征，以及它们在不同尺度下的表现。

WTMM算法是一种基于小波变换模极大值的方法，用于计算图像的多重分形谱。该方法主要利用小波变换对图像进行多尺度分解，提取出图像在不同尺度下的边缘信息。然后，通过对这些边缘信息进行统计分析，计算出图像的多重分形谱。

具体来说，WTMM算法的计算步骤如下：

对图像进行二维小波变换，得到一系列小波系数。

对每个尺度下的小波系数进行模极大值检测，提取出图像的边缘信息。

对提取出的边缘信息进行统计分析，计算出图像的多重分形谱。

3.2、WTMM算法原理

WTMM算法的数学公式主要包括以下几个部分：

3.2.1 二维小波变换

对图像f(x,y)进行二维小波变换，可以得到一系列小波系数Wf(x,y)，其中下标f表示小波变换的类型，如Haar小波、Daubechies小波等。二维小波变换的数学公式可以表示为：

Wf(x,y)=∫∫f(u,v)ψf(x−u,y−v)dudvWf(x,y) = \int \int f(u,v) \psi_f(x-u,y-v) du dvWf(x,y)=∫∫f(u,v)ψf(x−u,y−v)dudv

其中，ψf(x,y)是小波基函数。

3.2.2 模极大值检测

对每个尺度下的小波系数进行模极大值检测，可以提取出图像的边缘信息。具体地，对于每个像素位置(x,y)，如果满足以下两个条件：

|Wf(x,y)|≥|Wf(x+1,y)|,|Wf(x,y)|≥|Wf(x−1,y)|,|Wf(x,y)|≥|Wf(x,y+1)|,|Wf(x,y)|≥|Wf(x,y−1)||W_f(x,y)| \geq |W_f(x+1,y)|, |W_f(x,y)| \geq |W_f(x-1,y)|,|W_f(x,y)| \geq |W_f(x,y+1)|, |W_f(x,y)| \geq |W_f(x,y-1)||Wf(x,y)|≥|Wf(x+1,y)|,|Wf(x,y)|≥|Wf(x−1,y)|,|Wf(x,y)|≥|Wf(x,y+1)|,|Wf(x,y)|≥|Wf(x,y−1)|

则称该像素位置为模极大值点。

3.2.3 多重分形谱计算

通过对提取出的边缘信息进行统计分析，可以计算出图像的多重分形谱。具体地，可以用以下公式计算多重分形谱：

α=lim⁡ε→0log⁡|Wf(x,y)|log⁡ε\alpha = \lim_{\varepsilon \to 0} \frac{\log |W_f(x,y)|}{\log \varepsilon}α=limε→0logεlog⁡|Wf(x,y)|

其中，ε是小波变换的尺度参数，α是奇异指数，用于描述图像在不同尺度下的奇异程度。通过对所有模极大值点的奇异指数进行统计分析，可以得到图像的多重分形谱。

4.部分核心程序

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

%%
%对保存的多张图片读取并调用WTMM方法求图像的多重分形谱，对得到的结果保存其特征值
if sel == 1
   k = 1;
  
    for i=1:2*n1*n2;
         
        if i<=n1*n2
           k      = i;
           folder = 'save_images\1\';
           lists  = dir('save_images\1\*.jpg');        
        end
         
        if i<=2*n1*n2 & i>n1*n2
           k      = i - n1*n2;
           folder = 'save_images\2\';
           lists  = dir('save_images\2\*.jpg');        
        end       
               
         
        i
        %read an image
        I                               = imread(fullfile(folder,lists(k).name));
        %调用分形函数
        [qt,rt,ft,fft,Dt,feature_data]  = func_Wavelet_multifractal(I);
        q{i}                            = qt;
        r{i}                            = rt;    
        f{i}                            = ft;   
        ff{i}                           = fft;   
        D{i}                            = Dt;
        Feature{i}                      = feature_data;
    end
  
    save result.mat q r f ff D Feature
  
    K = 120;
    figure;
    plot(r{K},f{K},'-r>',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.9,0.0]);
    xlabel('奇异指数a');
    ylabel('多重分行谱f(a)') 
    grid on;
     
    figure;
    plot(q{K}+2,D{K},'-r>',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.9,0.0]);
    xlabel('q');
    ylabel('D(q)') 
    grid on;
     
    figure
    plot(q{K},r{K},'-r>',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.9,0.0]);
    title('ｑ和阿尔法a'); 
    xlabel('权重因子q');
    ylabel('奇异指数a');
    grid on;
     
    figure;
    plot(q{K},f{K},'-r>',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.9,0.0]);
    title('ｑ和f(a) '); 
    xlabel('权重因子q');
    ylabel('多重分行谱f(a)'); 
    grid on;
end
  
  
  
  
%%
%调用分类器对特征参数进行分类
if sel == 0
    load result.mat %q r f ff Feature
  
    K = 120;
    figure;
    plot(r{K},f{K},'-r>',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.9,0.0]);
    xlabel('奇异指数a');
    ylabel('多重分行谱f(a)') 
    grid on;
     
    figure;
    plot(q{K}+2,D{K},'-r>',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.9,0.0]);
    xlabel('q');
    ylabel('D(q)') 
    grid on;
     
    figure
    plot(q{K},r{K},'-r>',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.9,0.0]);
    title('ｑ和阿尔法a'); 
    xlabel('权重因子q');
    ylabel('奇异指数a');
    grid on;
     
    figure;
    plot(q{K},f{K},'-r>',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.9,0.0]);
    title('ｑ和f(a) '); 
    xlabel('权重因子q');
    ylabel('多重分行谱f(a)'); 
    grid on;   
     
     
     
     
    for i = 1:length(Feature)
        P(i) =  Feature{i}(3);
    end
  
    T =  [1*ones(1,length(Feature)/2),2*ones(1,length(Feature)/2)];
  
    t1                      = clock;                              %计时开始
    net                     = fitnet(65);
    net.trainParam.epochs   = 1000;                               %设置训练次数
    net.trainParam.goal     = 0.0001;                             %设置性能函数
    net.trainParam.show     = 1;                                  %每10显示
    net.trainParam.Ir       = 0.005;                              %设置学习速率
    net                     = train(net,P,T);                     %训练BP网络
    datat                   = etime(clock,t1);
    Nets                    = net;
    view(Nets);
    figure;
    plot(P,'b-*');
  
    y = sim(net,P);  
  
    figure;
    stem(y,'-bs',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9,0.0,0.0]);
    hold on
    plot(T,'-mo',...
    'LineWidth',1,...
    'MarkerSize',6,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.5,0.9,0.0]);
    hold on
    legend('预测数据','实际数据');
    title('输出1为第一类，输出2为第二类（即可对比实际的健康部分和肿瘤部分）');
     
     
    disp('预测正确率');
    error = 0;
    for i = 1:length(y)
        if i <= length(y)/2 
           if y(i) > 1.5
              error = error + 1;
           end
        else
           if y(i) < 1.5
              error = error + 1;
           end         
        end
    end
    1-error/length(y)
end