-
Notifications
You must be signed in to change notification settings - Fork 1
/
PST.m
131 lines (119 loc) · 5.91 KB
/
PST.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
function [out PST_Kernel]= PST(I,handles,Morph_flag)
%PST or Phase Stretch Transform is an operator that finds features in an
% image. PST takes an intensity image I as its input, and returns a
% binary image out of the same size as I, with 1's where the function finds sharp transitions in
% I and 0's elsewhere. PST function is also able to return the detected features in gray scale
% level (i.e. without thresholding).
%
% In PST, the image is first filtered by passing through a smoothing
% filter followed by application of a nonlinear frequency-dependent phase
% described by the PST phase kernel. The output of the transform is the
% phase in the spatial domain. The main step is the 2-D phase function
% (PST phase kernel) which is typically applied in the frequency domain.
% The amount of phase applied to the image is frequency dependent with
% higher amount of phase applied to higher frequency features of the
% image. Since sharp transitions, such as edges and corners, contain
% higher frequencies, PST emphasizes the edge information. Features can
% be further enhanced by applying thresholding and morphological operations.
% For more information please visit: https://en.wikipedia.org/wiki/Phase_stretch_transform
%
% [out PST_Kernel]= PST(I,handles, Morph_flag) takes the original image I and applies
% PST to it. PST kernel paramters are given using a handle variable:
%
% handles.LPF : Gaussian low-pass filter Full Width at Half Maximum (FWHM) (min : 0, max : 1)
% handles.Phase_strength : PST kernel Phase Strength (min : 0, max : 1)
% handles.Warp_strength : PST Kernel Warp Strength (min : 0, max : 1)
% handles.Thresh_min : minimum Threshold (min : -1, max : 0)
% handles.Thresh_max : maximum Threshold (min : 0, max : 1)
%
% Morph_flag allows user to compute the analog edge (if Morph_flag=0) or the digital edge (analog edge followed
% by thresholding and morphological operations)(if Morph_flag=1).
%
% Class Support
% -------------
% I is a double precision 2D array. out is of class logical.
%
% Remarks
% -------
% image processing toolbox is needed to run this function, function has
% been tested on MATLAB R2013a on a computer with Windows 7, 64 bits operating sytsem.
% The code uses IMOVERLAY function deveopled by Steven L. Eddins for visulaization of detected features.
%
% Example
% -------
% Example 1: Find the features of the circuit.tif image using PST method:
%
% I = imread('circuit.tif');
% I=double(I);
% handles.LPF=0.21;
% handles.Phase_strength=0.48;
% handles.Warp_strength=12.14;
% handles.Thresh_min=-1;
% handles.Thresh_max=0.004;
% [out PST_Kernel]= PST(I,handles);
% figure, imshow(out)
%
% Example 2: Find the features in the Lena image. See the attached test script.
%
% Copyright
% ---------
% PST function is developed in Jalali Lab at University of California,
% Los Angeles (UCLA). PST is a spin-off from research on the photonic time stretch technique in Jalali lab at UCLA.
% More information about the technique can be found in our group
% website: http://www.photonics.ucla.edu
% This function is provided for research purposes only. A license must be
% obtained from the University of California, Los Angeles for any commercial
% applications. The software is protected under a US patent.
% Citations:
% 1. M. H. Asghari, and B. Jalali, "Edge detection in digital images using dispersive phase stretch," International Journal of Biomedical Imaging, Vol. 2015, Article ID 687819, pp. 1-6 (2015).
% 2. M. H. Asghari, and B. Jalali, "Physics-inspired image edge detection," IEEE Global Signal and Information Processing Symposium (GlobalSIP 2014), paper: WdBD-L.1, Atlanta, December 2014.
%
% Copyright 1992-2016 The MathWorks, Inc.
% $Revision: 0.0.0.1 $ $Date: 2016/02/09 13:20:56 $
% define two dimentional cartesian (rectangular) vectors, 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);
% Convert cartesian X and Y vectors to polar vectors, THETA and RHO
[THETA,RHO] = cart2pol(X,Y);
% define two dimentional cartesian frequency vectors, FX and FY
X_step=x(2)-x(1);
fx=linspace(-0.5/X_step,0.5/X_step,length(x));
fx_step=fx(2)-fx(1);
Y_step=y(2)-y(1);
fy=linspace(-0.5/Y_step,0.5/Y_step,length(y));
fy_step=fy(2)-fy(1);
[FX,FY]=ndgrid(fx,fy);
% Convert cartesian vectors (FX and FY) to polar vectors, FTHETA and FRHO
[FTHETA,FRHO] = cart2pol(FX,FY);
% low pass filter the original image to reduce noise
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)));
% Constrcut the PST Kernel
PST_Kernel=(RHO*handles.Warp_strength.*atan(RHO*handles.Warp_strength)-0.5*log(1+(RHO*handles.Warp_strength).^2));
PST_Kernel=PST_Kernel/max(max(PST_Kernel))*handles.Phase_strength;
% Apply the PST Kernel
temp=(fft2(Image_orig_filtered)).*fftshift(exp(-1j*PST_Kernel));
Image_orig_filtered_PST=ifft2(temp);
% Calculate phase of the transformed image
PHI_features=angle(Image_orig_filtered_PST);
if Morph_flag ==0
out=PHI_features;
else
% find image sharp transitions by thresholding the phase
features=zeros(size(PHI_features));
features(find(PHI_features>handles.Thresh_max))=1;
features(find(PHI_features<handles.Thresh_min))=1; % because output phase has both positive and negative values
features(find(I<max(max(I))/20))=0; % ignore the features in the very dark areas of the image
% apply binary morphological operations to clean the transformed image
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