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

77 阅读5分钟

1.算法运行效果图预览

 

1.jpeg

2.jpeg

3.jpeg

4.jpeg

5.jpeg

6.png  

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:2n1n2;

       

        if i<=n1*n2

           k      = i;

           folder = 'save_images\1';

           lists  = dir('save_images\1*.jpg');       

        end

       

        if i<=2n1n2 & 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('q和阿尔法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('q和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('q和阿尔法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('q和f(a) ');

    xlabel('权重因子q');

    ylabel('多重分行谱f(a)');

    grid on;  

   

   

   

   

    for i = 1:length(Feature)

        P(i) =  Feature{i}(3);

    end

 

    T =  [1ones(1,length(Feature)/2),2ones(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`