% 相位拉伸matlab实现
% 作者:美国加州大学洛杉矶分校电子与计算机工程系 Jalali实验室M. Asghari
function [out, PST_Kernel]= PST(I,handles,Morph_flag)
%% 相位拉伸变换是在图像中发现特征的一种算子,它的输入是图像I的强度值,输出是一个和I同样大小的二值图像,
%二值图像种的1表示图像I中强度显著变化的点,0表示其它的地方。PST函数在没有阈值化的情况下,在灰度图像中也
%可以输出检测特征点。
% 相位拉伸变换操作步骤是:首先使用高斯低通滤波器对输入图像进行平滑滤波,去除图像噪音,
%接下来使用PST相位核描述与输入图像频率相关的非线性相位,最后变换输出的是空间域相位谱。
% 这其中最主要的操作是建立2D相位函数(即PST相位核),它在图像的频率域使用.
% 使用相位角的数量与频率相关,图像中频率特征越高,使用的相位数量就越多.由于显著变化,
% 比如边缘和角变换显著,它们包含了较高的频率,PST强调的时边缘信息.通过阈值化和形态学操作可以进一步增强边缘特征信息
%
%% 定义一个2维的笛卡尔向量, X and Y
Image_orig_size=size(I);
L=0.5;
x=linspace(-L,L,Image_orig_size(1));%创建一个线性空间向量,大小等于原始图像的宽度
y=linspace(-L,L,Image_orig_size(2));%创建一个线性空间向量,大小等于原始图像的高度
[X,Y]=ndgrid(x,y);%创建一个2维向量的矩形网格
% 把笛卡尔坐标下X,Y向量转换为极坐标下 THETA 和 RHO
[THETA,RHO] = cart2pol(X,Y);
% 定义2维的笛卡尔频率向量, FX and FY
X_step=x(2)-x(1);%步长
fx=linspace(-0.5/X_step,0.5/X_step,length(x));%创建一个线性空间向量,大小等于原始图像的宽度
Y_step=y(2)-y(1);
fy=linspace(-0.5/Y_step,0.5/Y_step,length(y));
[FX,FY]=ndgrid(fx,fy);
% 把FX,FY转换成极坐标, FTHETA 和 FRHO
[FTHETA,FRHO] = cart2pol(FX,FY);
% 原始图像缩减噪音的低通滤波器
%Image_orig_f=((fft2(I)));
%sigma=(handles.LPF)^2/log(2);
%Image_orig_f=Image_orig_f.*fftshift(exp(-(RHO/sqrt(sigma)).^2));
%Image_orig_filtered=real(ifft2((Image_orig_f)));
% Image_orig_filtered = L0Smoothing(I,0.01);%L0滤波替代高斯低通滤波
% 构造PST核
PST_Kernel=(RHO*handles.Warp_strength.*atan(RHO*handles.Warp_strength)-0.5*log(1+(RHO*handles.Warp_strength).^2));
%PST_Kernel=log(1+exp(RHO*handles.Warp_strength))-0.5*RHO*handles.Warp_strength;
PST_Kernel=PST_Kernel/max(max(PST_Kernel))*handles.Phase_strength;
% 使用PST核
%
temp=(frft2d(I,[1,1])).*exp(-1j*PST_Kernel);
Image_orig_filtered_PST=frft2d(temp,[-1,-1]);
% 计算变换图像的相位
PHI_features=angle(Image_orig_filtered_PST);
if Morph_flag ==0
out=PHI_features;
else
%通过阈值化相位,找到图像中显著变化的特征
features=zeros(size(PHI_features));
features(find(PHI_features>handles.Thresh_max))=1;
features(find(PHI_features<handles.Thresh_min))=1; % 因为输出相位具有正负值
features(find(I<max(max(I))/20))=0; % 忽略图效忠非常暗区域的特征
% 应用二值形态学特征,清理转换图像中单独点
out=features;
out = bwmorph(out,'thin',1);
out = bwperim(out,4);
out = bwmorph(out,'thin',1);
out = bwmorph(out,'clean',30); %清除孤立点
out = bwmorph(out,'remove',5);%移除区域内部像素
end