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.2 $  $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