1 内容介绍
本文提出的Step过滤和多项式过滤方法的实现 提供对部分曲线结构的增强。
该方法应用于地球物理图像，以识别线性模式存在于考古遗址 [1] 和故障检测 [2] 中的地下建筑结构。该方法有效地结合了旋转和尺度不变滤波器和像素标记方法，提供对大多数线结构的强大增强和检测分别为二维灰度图像。
2 仿真代码
% I : given image
% grammes: size,
% DTheta: angle vector
% DPaxos: vector of width of main wave
% typeOfFilter : type of step filter in {0,1} : 1 step function, 0 Polynomial filter
% Y: filter output
function [Y]= step_pol_Filtering(I,grammes, DTheta,DPaxos,typeOfFilter)
gabors = length(DTheta);
gabors2 = length(DPaxos);
for j=1:gabors2,
bw = DPaxos(j);
for i=1:gabors,
theta = DTheta(i);
if typeOfFilter == 1,
gb=mask_fn(grammes,bw,theta,0);
else
gb=mask_fn_smooth(grammes,bw,theta,0);
end
gb = gb / sum(sum(gb.^2)); %normalazation in energy
G = conv2(I,gb,'same'); %idio size me I
if i == 1 && j == 1,
Y = abs(G);
else
Y = max(abs(Y),abs(G));
end
end
end
del = 8;
Y(1:del,:) = 0;
Y(:,1:del) = 0;
Y(size(I,1)-[0:del-1],:) = 0;
Y(:,size(I,2)-[0:del-1]) = 0;
i = 1;
if typeOfFilter == 1,
gb=mask_fn(grammes,bw,pi/4,1);
else
gb=mask_fn_smooth(grammes,bw,pi/4,1);
end
function gb=mask_fn(grammes,bw,theta,toPlot)
% bw = bandwidth, (1)
% gamma = aspect ratio, (0.5)
%
% lambda= wave length, (>=2)
% theta = angle in rad, [0 pi)
sz=grammes;
if mod(sz,2)==0, sz=sz+1;end
[x y]=meshgrid(-fix(sz/2):1:fix(sz/2),fix(sz/2):-1:fix(-sz/2));
% Rotation
x_theta=x*cos(theta)+y*sin(theta);
y_theta=-x*sin(theta)+y*cos(theta);
for i=1:length(x),
for j=1:length(y),
d = abs(x_theta(i,j));
u = x_theta(i,j);
v = y_theta(i,j);
gb(i,j) = 0;
if d < bw,
gb(i,j) = 1;
elseif d <= 3*bw
gb(i,j) = -1;
end
end
end
[n1 n2] = find(gb == -1);
[p1 p2] = find(gb == 1);
nv_n1 = -length(p1)/length(n1);
for i=1:length(n1),
gb(n1(i),n2(i))= nv_n1;
end
gb = gb/sum(sum(gb.^2));
%gb=exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+psi);
if toPlot == 1,
figure;
imagesc(gb);
colormap(gray);
title('Filter');
mean(mean(gb))
end
function gb=mask_fn_smooth(grammes,bw,theta,toPlot)
% bw = bandwidth, (1)
% grammes = filter size
%
% toPlot= 1 for plotting
% theta = angle in rad, [0 pi)
sz=grammes;
if mod(sz,2)==0, sz=sz+1;end
[x y]=meshgrid(-fix(sz/2):1:fix(sz/2),fix(sz/2):-1:fix(-sz/2));
% Rotation
x_theta=x*cos(theta)+y*sin(theta);
y_theta=-x*sin(theta)+y*cos(theta);
for i=1:length(x),
for j=1:length(y),
d = abs(x_theta(i,j));
u = x_theta(i,j);
v = y_theta(i,j);
gb(i,j) = 0;
if d < bw,
gb(i,j) = 1-((d/bw)^2);
elseif d <= 3*bw
gb(i,j) = 0.5*(-1+(((d/bw)-2)^2));
end
end
end
[u1 v1] = find(gb > 0);
[u2 v2] = find(gb < 0);
m1 = 0;
N1 = length(u1);
for i=1:length(u1),
m1 = m1+gb(u1(i),v1(i));
end
m1 = m1/N1;
m2 = 0;
N2 = length(u2);
for i=1:length(u2),
m2 = m2+gb(u2(i),v2(i));
end
m2 = m2/N2;
c = -m1*N1/(m2*N2);
for i=1:length(u2),
gb(u2(i),v2(i)) = c*gb(u2(i),v2(i));
end
gb = gb/sum(sum(gb.^2));
%gb=exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+psi);
if toPlot == 1,
figure;
imagesc(gb);
colormap(gray);
title('Filter');
figure;
mesh(gb);
title('Filter in 3D view');
mean(mean(gb))
end
3 运行结果
编辑
编辑
编辑
4 参考文献
[1] C. Panagiotakis, E. Kokinou and A. Sarris, Curvilinear Structure Enhancement and Detection in Geophysical
Images, IEEE Trans. on Geoscience and Remote Sensing, vol. 49, no. 6, pp. 2040-2048, 2011.
[2] C. Panagiotakis and E. Kokinou, Linear Pattern Detection of Sea Faults via a Topology and Shape Optimization Method,
IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 8, no. 1, pp. 3-11, 2015.
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。
