【WTMM】基于WTMM的图像多重分形谱matlab仿真

1.软件版本

matlab2013b

2.系统程序

提出一种新的多重分形图象去噪算法,讨论了基于小波极大模的多重分形谱估计算法.在此基础上推导了图像取噪声的谱移位算子.该方法没有对噪声的类型提出任何假设条件,而是通过定义一个变换算子对每一点的 Hausdorf指数进行处理,使处理后的图象的Hausdorf指数接近于2,从而取得最佳效果.实验结果表明,该方法在去除噪声的同时可很好地保留了原始图像的纹理信息.

算法流程图如图所示：

部分程序如下：

clc;
clear;
close all;
warning off;
addpath 'FUNC\'
sel = 0;%1：计算多重分形谱；0：分类
%%
%读取三个图片
I1s = imread('图像\D10.gif');
I1 = imresize(I1s,[512,512]);
I2s = imread('图像\D05.gif');
I2 = imresize(I2s,[512,512]);
figure;
subplot(121);imshow(I1);title('图片1');
subplot(122);imshow(I2);title('图片2');
%%
%分成多个子图片
[M,N] = size(I1);
sub_length = 64;
n1 = floor(M/sub_length);
n2 = floor(N/sub_length);
DF = zeros(1,n1*n2);
for num_image = 1:2
k=1;
for i=1:n1
for j=1:n2
x_start = (i-1)*sub_length+1;
x_end = i*sub_length;
y_start = (j-1)*sub_length+1;
y_end = j*sub_length;
if num_image == 1
sub_image = I1(x_start:x_end,y_start:y_end);
imwrite(sub_image,['save_images\1\I_',num2str(k),'.jpg']);
end
if num_image == 2
sub_image = I2(x_start:x_end,y_start:y_end);
imwrite(sub_image,['save_images\2\I_',num2str(k),'.jpg']);
end
%保存分割的多个子图片
k = k + 1;
end
end
end
%%
%对保存的多张图片读取并调用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},'b-o');
xlabel('奇异指数a');
ylabel('多重分行谱f(a)')
grid on;
figure;
plot(q{K}+2,D{K},'b-o');
xlabel('q');
ylabel('D(q)')
grid on;
figure
plot(q{K},r{K},'b-o');
title('ｑ和阿尔法a');
xlabel('权重因子q');
ylabel('奇异指数a');
grid on;
figure;
plot(q{K},f{K},'b-o');
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},'b-o');
xlabel('奇异指数a');
ylabel('多重分行谱f(a)')
grid on;
figure;
plot(q{K}+2,D{K},'b-o');
xlabel('q');
ylabel('D(q)')
grid on;
figure
plot(q{K},r{K},'b-o');
title('ｑ和阿尔法a');
xlabel('权重因子q');
ylabel('奇异指数a');
grid on;
figure;
plot(q{K},f{K},'b-o');
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,'r-*');hold on
plot(T,'b-','LineWidth',2);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