function [imgf_temporal,Filter] = DirectionFilter_Gaussian(img,TotalDirectionNumber,DirectionIndex,cutin,cutoff,dbBandwith)
% Spatio-Temporal directional filtering
% lgdf = DirectionFilter(img,nOrt,Orz,cutin,cutoff,nOrd)
% img: input 3D displacement snapshots [x,y,t]
% nOrt: number of directions of signal decomposition
% Orz: interested direction
% ! it is counter clockwise
% cutin: optional low-frequency cutting ratio (0-0.5)
% cutoff: optional high-frequency cutting ratio (0-0.5)
% nOrd: optional order of butterworth bandpass filter
% LI Bing Nan @ NUS, March 2010
% Reference: Manduca et al (2003) MIA 465-73
% Hu Liangliang@ HFUT ,April 2017
% Add Time direction
if nargin<6, dbBandwith=3; end
if nargin<5, cutoff=0.35; end
if nargin<4, cutin=0.01; end
if nargin<3, error('Too few input parameters'); end
if ~isa(img,'double'), img=double(img); end
nOrt=TotalDirectionNumber;
Orz=DirectionIndex;
[rn,cn,sn]=size(img);
CenterFreq=(cutin+cutoff)/2;
SingalSideBandwith=(cutoff-cutin)/2;
Sigma2=SingalSideBandwith.^2/(dbBandwith/10*log(10));
x = (ones(rn,1) * (1:cn) - (fix(cn/2)+1))/cn;
y = ((1:rn)' * ones(1,cn) - (fix(rn/2)+1))/rn;
radius = sqrt(x.^2 + y.^2);
% %Butterworth filter
% fo = 1 ./ (1.0 + (radius ./ cutoff).^(2*nOrd));
% fi = 1 ./ (1.0 + (radius ./ cutin).^(2*nOrd));
% lgdf=fo-fi;
% Gaussian Filter
lgdf=exp(-(radius-CenterFreq).^2/(2*Sigma2));
rdtheta=atan2(-y,x);
sintheta=sin(rdtheta);
costheta=cos(rdtheta);
clear x; clear y; clear rdtheta; % save a little memory
dThetaOnSigma = sqrt(nOrt); % angular ratio modulator
thetaSigma = 2*pi/nOrt/dThetaOnSigma;
theta=(2*pi/nOrt)*Orz;
ds=sintheta.*cos(theta)-costheta.*sin(theta);
dc=costheta.*cos(theta)+sintheta.*sin(theta);
dtheta=abs(atan2(ds,dc));
spread = exp((-dtheta.^2) / (2*thetaSigma^2)); % Calculate the angular filter component.
lgdf=lgdf.*spread;
imgf=zeros(size(img));
% Add by HLL
if (sn>3) %At least 4 frame per-peroid
img_temporalF=fftshift(fft(img,sn,3),3);% Add by HLL
else
img_temporalF=img;
end
for k=1:sn
% imgs(:,:)=img(:,:,k);% org
imgs(:,:)=img_temporalF(:,:,k);% Changed
imgfft=fft2(imgs);
imgfft=fftshift(imgfft);
imgft=imgfft.*lgdf;
imgft=ifftshift(imgft);
imgs=ifft2(imgft);
imgf(:,:,k)=imgs;
end
Filter=lgdf;
if (sn>3) %At least 4 frame per-peroid
imgf_temporal=ifftshift(imgf,3);%ADD by HLL
else
imgf_temporal=imgf;
end
评论3