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

小波模极大值原理在图像边缘提取和信号奇异点检测中的应用《转》
2007-05-26 10:18
不做小波很久，陆续接到网友的很多询问，不少信件关于这个话题。

本不想花功夫写程序，因为毕竟研究方向是计算电磁学，然对小波的好奇仍是一种抗拒不了的力量。

再加上网友的一遍遍不厌其烦的请教，我也就利用半天时间，将这一话题做了一个程序，拿出来分享。

1。

什么是模极大值？一般信号的主要信息，由拐点（二阶导数为零的点）确定，而由于噪声的影响，直接求拐点显然困难。

于是，我们求一阶导数的模的极大值。

2。

什么是小波模极大值？就是先将小波函数和原信号卷积（连续小波变换），然后对结果取模，最后找到极大值。

上述步骤，也就等价于：先把某一光滑函数求导（求导后满足积分为零的条件成为小波函数），然后卷积源信号，接着取模，最后发现极大值。

3。

图像处理的操作。

a、给定某一尺度，求出二维高斯函数沿x和沿y方向的导数Phi_x,Phi_y。

这两个函数就等价于小波函数。

b、用Phi_x,Phi_y分别与图像卷积得到Gx，Gy。

c、求出每一个像素点的梯度大小G=(Gx*Gx+Gy*Gy).^(1/2)，用反正切求梯度方向或者称幅角atan(Gy/Gx)。

这里，注意的是反正切只能求出一、四象限的角度，其它象限要分别处理。

且Gx为一个很小的数值时，也要处理。

d、把求得幅角，分成四种方向。

第一种0或180方向（水平），第二种90或270方向（垂直），第三种45或225方向（正对角线），第四种135或315方向（负对角线）。

也就是说，看看你求出幅角的大小与上面的哪个方向最接近。

e、依次检测每一个像素点，看看在它对应“幅角最接近的方向上”是否是极大值。

如果是，纪录该梯度值。

若不是，把梯度值置零。

f、找到记录梯度值中的最大值，然后以该值做归一化。

比较每一个像素归一化的梯度值，当该梯度值大于某个阈值的时候，就是真正边缘，否则认为是伪边缘。

4。

实际上这个算法和canny算子本质上等价的。

让我们再来回顾canny本人经典的原话，来体会边缘提取的目标到底是什么。

a、好的检测性能。

不漏检真实边缘，也不把非边缘点作为边缘点检出，使输出的信噪比最大。

b、好的定位性能。

检测到的边缘点与实际边缘点位置最近。

c、唯一性。

对于单个边缘点仅有一个响应。

沙威(gjsdgjsd)
安徽大学
2007年4月22日
% 小波模极大值用于边缘提取
% 沙威(gjsdgjsd)安徽大学
% 2007年4月22日
clc;clear
% 下载图像
load woman
% X=double(imread('1.bmp'));
SIZE=length(X); % 图像尺寸
% 多尺度
m=1.0;
delta=2^m;
% 构造高斯函数的偏导
N=20; % 滤波器长度（需要调整，必须是偶数）
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)=(x/delta^2).*exp(-(x.*x+y.*y) /(2*delta^2));
phi_y(index_x,index_y)=(y/delta^2).*exp(-(x.*x+y.*y) /(2*delta^2));
end
end;
% 对图象做行列卷积
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('边缘')
% 进一步的工作，连线...
对于信号的奇异性检测，也是一个道理。

% 小波变换用于奇异检测
% 安徽大学沙威
% 2007年4月27日
clc;clear
% 下载信号
load freqbrk;
s=freqbrk;
n=length(s); % 长度
% 多尺度
m=-2.0;
delta=2^m;
% 构造高斯函数的偏导
N=20; % 滤波器长度（需要调整，必须是偶数）
A=-1/sqrt(2*pi); % 幅度
% 构造高斯函数
for index_x=1:N;
x=index_x-(N+1)/2;
phi_x(index_x)=A*(x/delta^2).*exp(-(x.*x)/(2*delta^2)); end;
phi_x=phi_x/norm(phi_x); % 能量归一化
% 对信号做卷积
g=conv(s,phi_x); % 卷积
g=wkeep(g,n); % 保持信号长度
% 模极大值
m=abs(g); % 取模
[V,P]=max(m); % 最大值位置
% 显示结果
figure(1);
subplot(2,1,1);
plot(0:n-1,s); % 原始信号
xlabel('离散时间');
ylabel('信号幅值');
title('原始信号')
subplot(2,1,2);
hold on;
plot(0:n-1,m); % 信号经小波变换后的模值
plot(P-1,V,'r*') % 小波模最大值
title(['小波模最大值（信号突变点）对应的时刻是',num2str(P-1)]) xlabel('离散时间');
ylabel('小波变换后的模值');。