g=wft2f(type,f,sigmax,wxl,wxi,wxh,sigmay,wyl,wyi,wyh,thr)

function g=wft2f(type,f,sigmax,wxl,wxi,wxh,sigmay,wyl,wyi,wyh,thr)
%This function and all the input and output are basicly the same as 
%in the function "wft2" shown in the OLEN paper. 
%Only the realization of convolution is slighly different.
%Fourier transfomrs are used to realize conv2.
%wft2f comsumes about 75% of wft2s time, and much less than wft2.
%However, wft2f uses more memory than wft2s.
%further, sigma is changed to sigmax and sigmay so that the parameter
%setting can be more general
%--------------------------------------------------------------------------
%type:  'wff' or 'wfr'
%f:     2D input signal, (1) exponential field from phase-shifting technique 
%       and digital holography; (2)carrier fringe from carrier technique, 
%       (3) it can also be a closed fringe pattern
%sigmax: window size along x, recomended value: 10
%wxl:   low bound of freuqency in x
%wxi:   increasement of wx;
%wxh:   high bound of frequency in x
%sigmay: window size along y, recomended value: 10
%wyl:   low bound of freuqency in y
%wyi:   increasement of wx
%wyh:   high bound of frequency in y
%thr:   threshold for 'wff', not needed for 'wfr', recomended
%       value:3*standard deviation of noise
%g:     For 'wff', g.filtered is 2D filtered signal, for phase-shifting and
%       carrier technique, use "phase=angle(g)" to get the phase map; 
%       for closed fringe, use "h=real(g.filtered)" to get smoothed 
%       fringe pattern. wxi<=1/sigmax and wyi<=1/sigmay are required.
%       For 'wfr', g is a structure, ridge value, wx, wy and phase are
%       instored as g.r,g.wx,g.wy and g.phase. 
%       wxi and wyi are set according to required frequency resolution
%interesting: fft2 0f 297*297 is much faster than 296*296 and
%298*298
%--------------------------------------------------------------------------
%Example:   g=wft2f('wff',f,10,-0.5,0.1,0.5,10,-0.5,0.1,0.5,6);
%           g=wft2f('wfr',f,10,-0.5,0.1,0.5,10,-0.5,0.1,0.5);
%--------------------------------------------------------------------------
%References:
%[1]Windowed Fourier transform for fringe pattern analysis,
%   Applied Optics, 43(13):2695-2702, 2004
%[2]Two-dimensional windowed Fourier transform for fringe pattern analysis: 
%   principles, applications and implementations, 
%   Optics and Lasers in Engineering, 45(2): 304-317, 2007.
%[3]Windowed Fourier transform for fringe pattern analysis: 
%   theoretical analyses, Applied Optics, 47(29): 5408-5419, 2008
%[4]A windowed Fourier filtered and quality guided phase unwrapping 
%   algorithm, Applied Optics, 47(29):5420-5428, 2008
%--------------------------------------------------------------------------
%Last update: 27 Nov 2008
%Contact: [e-mail] (Dr Qian Kemao, Nanyang Technological Univ.)
%All Copyrights reserved.
%--------------------------------------------------------------------------

%the purpose is to make the wff faster by choosing smaller window size
%which does not affect its accuracy.
if strcmp(type,'wff')
    %half window size along x, by default 2*sigmax; window size: 2*sx+1
    sx=round(2*sigmax); 
    %half window size along y, by default 2*sigmay; window size: 2*sy+1
    sy=round(2*sigmay); 
elseif strcmp(type,'wfr')
    %half window size along x, by default 3*sigmax; window size: 2*sx+1
    sx=round(3*sigmax); 
    %half window size along y, by default 3*sigmay; window size: 2*sy+1
    sy=round(3*sigmay); 
end

%image size
[m n]=size(f); 
%expanded size: size(A)+size(B)-1
mm=m+2*sx;nn=n+2*sy;
%expand f to size [mm nn]
f=fexpand(f,mm,nn); 
%pre-compute the spectrum of f; 
Ff=fft2(f); 
%meshgrid (2D index) for window
[y x]=meshgrid(-sy:sy,-sx:sx); 
%generate a window w0
w0=exp(-x.*x/2/sigmax/sigmax-y.*y/2/sigmay/sigmay); 
%norm2 normalization
w0=w0/sqrt(sum(sum(w0.*w0))); 

if strcmp(type,'wff')
    %to store result
    g.filtered=zeros(m,n);
    for wyt=wyl:wyi:wyh
        for wxt=wxl:wxi:wxh
            %WFT basis
            w=w0.*exp(j*wxt*x+j*wyt*y);
            %expand w to size [mm nn]
            w=fexpand(w,mm,nn); 
            %spectrum of w
            Fw=fft2(w);  
            %implement of WFT: conv2(f,w)=ifft2(Ff*Fw);
            sf=ifft2(Ff.*Fw); 
            %cut to get desired data size
            sf=sf(1+sx:m+sx,1+sy:n+sy); 
            %threshold the spectrum
            sf=sf.*(abs(sf)>=thr); 
            %exapand sf to size [mm nn]
            sf=fexpand(sf,mm,nn); 
            %implement of IWFT: conv2(sf,w);
            gtemp=ifft2(fft2(sf).*Fw); 
            %update
            g.filtered=g.filtered+gtemp(1+sx:m+sx,1+sy:n+sy); 
        end
    end
    %scale the data
    g.filtered=g.filtered/4/pi/pi*wxi*wyi; 
elseif strcmp(type,'wfr')
    %to store wx, wy, phase and ridge values
    g.wx=zeros(m,n); g.wy=zeros(m,n); g.phase=zeros(m,n); g.r=zeros(m,n);
    for wyt=wyl:wyi:wyh
        for wxt=wxl:wxi:wxh
            %WFT basis
            w=w0.*exp(j*wxt*x+j*wyt*y);
            %expand w to size [mm nn]
            w=fexpand(w,mm,nn); 
            %spectrum of w 
            Fw=fft2(w); 
            %implement of WFT: conv2(f,w)=ifft2(Ff*Fw);
            sf=ifft2(Ff.*Fw); 
            %cut to get desired data size
            sf=sf(1+sx:m+sx,1+sy:n+sy);
            %indicate where to update
            t=(abs(sf)>g.r); 
            %update r
            g.r=g.r.*(1-t)+abs(sf).*t; 
            %update wx
            g.wx=g.wx.*(1-t)+wxt*t; 
            %update wy
            g.wy=g.wy.*(1-t)+wyt*t; 
            %update phase
            g.phase=g.phase.*(1-t)+angle(sf).*t; 
        end
    end
end

function f=fexpand(f,mm,nn)
%expand f to [m n]
%this function can be realized by padarray, but is slower

%size f
[m n]=size(f); 
%store f
f0=f;
%generate a larger matrix with size [mm nn]
f=zeros(mm,nn);
%copy original data
f(1:m,1:n)=f0;