小波模极大值用于边缘提取.txt

97 law to enhance the classic procedure Ridge wavelet extraction Modulus maximum for the wavelet

%  小波模极大值用于边缘提取

 
clc;clear
%  下载图像
load woman
% X=double(imread('1.bmp'));
SIZE=length(X);  %  图像尺寸
    
%  多尺度
m=1.0;
delta=2^m;
%  构造高斯函数的偏导
N=20;  %  滤波器长度（需要调整，必须是偶数）
A=-1/sqrt(2*pi);  %  幅度
for index_x=1:N;
    for index_y=1:N;
        x=index_x-(N+1)/2;
        y=index_y-(N+1)/2;
        phi_x(index_x,index_y)=A*(x/delta^2).*exp(-(x.*x+y.*y)/(2*delta^2));
        phi_y(index_x,index_y)=A*(y/delta^2).*exp(-(x.*x+y.*y)/(2*delta^2));
    end
end;
phi_x=phi_x/norm(phi_x);  %  能量归一化
phi_y=phi_y/norm(phi_y);  %  能量归一化
%  对图象做行列卷积
Gx=conv2(X,phi_x,'same');
Gy=conv2(X,phi_y,'same');
%  求梯度
Grads=sqrt((Gx.*Gx)+(Gy.*Gy));
%  求幅角（梯度方向）
angle_array=zeros(SIZE,SIZE);  %  角度
%  遍历
for i=1:SIZE;
    for j=1:SIZE
        if (abs(Gx(i,j))>eps*100)  %  x的绝对值足够大
            p=atan(Gy(i,j)/Gx(i,j))*180/pi;  %  反正切求角度值(1,4象限)
            if (p<0)        %  负的幅角（4象限）
                p=p+360;
            end;
            if (Gx(i,j)<0 & p>180)     %  2象限的特殊处理
                p=p-180;
            elseif (Gx(i,j)<0 & p<180) %  3象限的特殊处理
                p=p+180;
            end
        else  %  90或270度
            p=90;
        end
        angle_array(i,j)=p;  %  幅角赋值
    end
end;
%  找边缘
edge_array=zeros(SIZE,SIZE);
%  遍历
for i=2:SIZE-1
    for j=2:SIZE-1
        if ((angle_array(i,j)>=(-22.5) & angle_array(i,j)<=22.5) | ...
            (angle_array(i,j)>=(180-22.5) & angle_array(i,j)<=(180+22.5)))     %  0/180
            if (Grads(i,j)>Grads(i+1,j) & Grads(i,j)>Grads(i-1,j))
                edge_array(i,j)=Grads(i,j);
            end
        elseif ((angle_array(i,j)>=(90-22.5) & angle_array(i,j)<=(90+22.5)) | ...
                (angle_array(i,j)>=(270-22.5) & angle_array(i,j)<=(270+22.5))) %  90/270
            if (Grads(i,j)>Grads(i,j+1) & Grads(i,j)>Grads(i,j-1))
                edge_array(i,j)=Grads(i,j);
            end
        elseif ((angle_array(i,j)>=(45-22.5) & angle_array(i,j)<=(45+22.5)) | ...
                (angle_array(i,j)>=(225-22.5) & angle_array(i,j)<=(225+22.5))) %  45/225
            if (Grads(i,j)>Grads(i+1,j+1) & Grads(i,j)>Grads(i-1,j-1))
                edge_array(i,j)=Grads(i,j);
            end
        else  %  135/215
            if (Grads(i,j)>Grads(i+1,j-1) & Grads(i,j)>Grads(i-1,j+1))
                edge_array(i,j)=Grads(i,j);
            end
        end
    end
end
%  去除伪边缘
MAX_E=max(max(edge_array).');     %  最大幅度值
edge_array=edge_array/MAX_E;      %  最大幅度值
threshold=0.2;                    %  阈值（需要调整）
%  遍历
for m=1:SIZE
    for n=1:SIZE
        if (edge_array(m,n)>threshold)
            edge_array(m,n)=1;
        else
            edge_array(m,n)=0;
        end
    end
end
%  显示图像和边缘
figure(1)
subplot(1,2,1)
imshow(X,map)
title('图像')
subplot(1,2,2)
imshow(edge_array)
title('边缘')
% 进一步的工作，连线...
对于信号的奇异性检测，也是一个道理。

%  小波变换用于奇异检测
%  编程人: 沙威(Wei Sha) 安徽大学(Anhui University) [e-mail] 

clc;clear; 

%  下载信号
load freqbrk;
s=freqbrk;  %  s 原始信号 

%  原始信号长度
n=length(s); 

%  构造小波函数
delta=-2;           %  delta 小波变换的尺度
A=-2/sqrt(pi);      %  幅度
delta_t=2^(delta);  %  尺度
N=20;               %  N 高斯函数的长度 

for index_x=1:N;
    x=index_x-(N+1)/2;
    phi_x(index_x)=A*(x/delta_t.^(3/2)).*exp(-(x.*x)/(delta_t^2));
end;
phi_x=phi_x/norm(phi_x);  %  能量归一化 

%  对信号做卷积
g=conv(s,phi_x);    %  卷积
g=wkeep(g,n);       %  保持信号长度 

%  模极大值
maxima_t=abs(g);    %  取模
maxima=zeros(n,1);  %  模极大值数组 

for i=2:n-1;
    if ((maxima_t(i)>=maxima_t(i-1)) & (maxima_t(i)>=maxima_t(i+1)))
        maxima(i)=maxima_t(i);
    end
end; 

% 去除小的模极大值
maxima=maxima/max(maxima);
threshold=0.6;  %  阈值 

for i=2:n-1;
    if (maxima(i)<threshold)
        maxima(i)=0;
    end
end; 

% 显示结果
figure(1); 

subplot(2,1,1);
plot(0:n-1,s); % 原始信号
xlabel('离散时间');
ylabel('信号幅值');
title('原始信号') 

subplot(2,1,2);
hold on;
plot(0:n-1,maxima); % 信号经小波变换后的模极大值 

xlabel('离散时间');
ylabel('模极大值');