基于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ε→0​logεlog⁡|Wf​(x,y)|​

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

4.部分核心程序

%%
%对保存的多张图片读取并调用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