function CtSent_grid_wrapped(savepath,phsbas,xy,pix_size,prefilt_win,goldfilt_flag,lowfilt_flag,gold_alpha)
%CtSent_GRID_WRAPPED resample unwrapped phase to a grid and filter
%
%   Andy Hooper, June 2006
%   Zhangfeng Ma, April 2021
% ============================================================
% 10/2008 AH: Amended to deal with zero phase values
% 08/2009 AH: Goldstein alpha value added to calling parms
% 02/2012 AH: save ij
% 03/2012 AH: Allow for non-complex wrapped phase
% 01/2013 AH: No gridding option (pix_size=0)
% 04/2013 AH: Correct bug calling wrap_filt
% 02/2014 AH: Add predefined ph_uw option
% 01/2017 DB: add stamps_save for large variables
% 04/2021 ZM: Mofied for CtSent NRT and Added SPAFU for windows
% ============================================================
% ll0 = (max(lonlat)+min(lonlat))/2;
% xy = llh2local(lonlat',ll0)*1000;
% xy = xy';
fprintf('   Resampling phase to grid...\n')

if nargin<3
    error('not enough arguments')
end
if nargin<4
    pix_size = 50;
end
if nargin<5
    prefilt_win = 20;
end
if nargin<6
    goldfilt_flag = 'y';
end
if nargin<7
    lowfilt_flag = 'n';
end
if nargin<8
    gold_alpha = 1;
end

[n_ps,n_ifg]=size(phsbas);

fprintf('   Number of interferograms  : %d\n',n_ifg)
fprintf('   Number of points per ifg  : %d\n',n_ps)

if ~isreal(phsbas) && sum(phsbas(:)==0)>0
    error('Containing Zero Phase Values!')
end
% Define the resampled grid size
if pix_size==0
    grid_x_min = 1;
    grid_y_min = 1;
    n_i = max(xy(:,2));
    n_j = max(xy(:,1));
    grid_ij = [xy(:,2),xy(:,1)];
else
    grid_x_min = min(xy(:,1));
    grid_y_min = min(xy(:,2));

    grid_ij(:,1)=ceil((xy(:,2)-grid_y_min+1e-3)/pix_size);
    grid_ij(grid_ij(:,1)==max(grid_ij(:,1)),1)=max(grid_ij(:,1))-1;
    grid_ij(:,2)=ceil((xy(:,1)-grid_x_min+1e-3)/pix_size);
    grid_ij(grid_ij(:,2)==max(grid_ij(:,2)),2)=max(grid_ij(:,2))-1;

    n_i = max(grid_ij(:,1));
    n_j = max(grid_ij(:,2));
end

ph_grid = zeros(n_i,n_j,'single');
 
if min(size(ph_grid))<prefilt_win
    error(['Minimum dimension of the resampled grid (',num2str(min(size(ph_grid))),...
        ' pixels) is less than prefilter window size (',num2str(prefilt_win),')'])
end

for i1=1:n_ifg
    if isreal(phsbas)
        ph_this = exp(1i*phsbas(:,i1));
    else
        ph_this = phsbas(:,i1);
    end   
    ph_grid(:) = 0;
    if pix_size==0
        ph_grid((xy(:,1)-1)*n_i+xy(:,2)) = ph_this;
        if predef_flag=='y'
           ph_grid_uw((xy(:,1)-1)*n_i+xy(:,2)) = ph_this_uw;
        end
    else
        for i=1:n_ps     
            ph_grid(grid_ij(i,1),grid_ij(i,2)) = ph_grid(grid_ij(i,1),grid_ij(i,2))+ph_this(i);
        end
    end
  
    if i1==1
        nzix=ph_grid~=0;
        n_ps_grid=sum(nzix(:));
        ph=zeros(n_ps_grid,n_ifg,'single');
        if strcmpi(lowfilt_flag,'y')
            ph_lowpass=ph;
        else
            ph_lowpass=[];
        end
    end
    if strcmpi(goldfilt_flag,'y') || strcmpi(lowfilt_flag,'y')
        [ph_this_gold,ph_this_low] = wrap_filt(ph_grid,prefilt_win,gold_alpha,[],lowfilt_flag);
        if strcmpi(lowfilt_flag,'y')
            ph_lowpass(:,i1)=ph_this_low(nzix);
        end
    end
    if strcmpi(goldfilt_flag,'y')
        ph(:,i1)=ph_this_gold(nzix);
    else
        ph(:,i1)=ph_grid(nzix);
    end
end

n_ps = n_ps_grid;

fprintf('   Number of resampled points: %d\n',n_ps)

[nz_i,nz_j]=find(ph_grid~=0);
if pix_size==0
    xy = xy;
else
    xy=[[1:n_ps]',(nz_j-0.5)*pix_size,(nz_i-0.5)*pix_size];
end
ij=[nz_i,nz_j];
oldFolder = cd(savepath);
save('phwgrid.mat','ph','ph_lowpass','xy','ij','nzix','grid_x_min',...
    'grid_y_min','n_i','n_j','n_ifg','n_ps','grid_ij','pix_size','-v7.3');
newfolder = cd(oldFolder);
end
function [ph_out,ph_out_low] = wrap_filt(ph,n_win,alpha,n_pad,low_flag)
%WRAP_FILT Goldstein adaptive and lowpass filtering
%   [ph_out]=wrap_filt(ph,n_win,alpha,n_pad)
%
%   Andy Hooper, June 2006
%   
%   Change Log:
%   07/2006 AH: Added zero padding 
%   02/2012 AH: Set magnitude to original values

if nargin<4 || isempty(n_pad)
    n_pad=round(n_win*0.25);
end

if nargin<5
    low_flag='n';
end
    
[n_i,n_j]=size(ph);
n_inc=floor(n_win/2);
n_win_i=ceil(n_i/n_inc)-1;
n_win_j=ceil(n_j/n_inc)-1;

ph_out=zeros(size(ph));
if strcmpi(low_flag,'y')
    ph_out_low=ph_out;
else
    ph_out_low=[];
end
x=[1:n_win/2];
[X,Y]=meshgrid(x,x);
X=X+Y;
wind_func=[X,fliplr(X)];
wind_func=[wind_func;flipud(wind_func)];


ph(isnan(ph))=0;
B=gausswin(7)*gausswin(7)';
ph_bit=zeros(n_win+n_pad);

L=ifftshift(gausswin(n_win+n_pad,16)*gausswin(n_win+n_pad,16)');

for ix1=1:n_win_i
    wf=wind_func;
    i1=(ix1-1)*n_inc+1;
    i2=i1+n_win-1;
    if i2>n_i
        i_shift=i2-n_i;
        i2=n_i;
        i1=n_i-n_win+1;
        wf=[zeros(i_shift,n_win);wf(1:n_win-i_shift,:)];
    end
    for ix2=1:n_win_j
        wf2=wf;
        j1=(ix2-1)*n_inc+1;
        j2=j1+n_win-1;
        if j2>n_j
           j_shift=j2-n_j;
           j2=n_j;
           j1=n_j-n_win+1;
           wf2=[zeros(n_win,j_shift),wf2(:,1:n_win-j_shift)];
        end
        ph_bit(1:n_win,1:n_win)=ph(i1:i2,j1:j2);
        ph_fft=fft2(ph_bit);
        H=abs(ph_fft);
        H=ifftshift(filter2(B,fftshift(H))); % smooth response
        meanH=median(H(:));
        if meanH~=0
            H=H/meanH;
        end
        H=H.^alpha;
        ph_filt=ifft2(ph_fft.*H);
        ph_filt=ph_filt(1:n_win,1:n_win).*wf2;
        if strcmpi(low_flag,'y')
            ph_filt_low=ifft2(ph_fft.*L);
            ph_filt_low=ph_filt_low(1:n_win,1:n_win).*wf2;
        end
        if isnan(ph_filt(1,1))
            disp('filtered phase contains NaNs in goldstein_filt')
            keyboard
        end
        ph_out(i1:i2,j1:j2)=ph_out(i1:i2,j1:j2)+ph_filt;
        if strcmpi(low_flag,'y')
            ph_out_low(i1:i2,j1:j2)=ph_out_low(i1:i2,j1:j2)+ph_filt_low;
        end
    end
end

ph_out=abs(ph).*exp(1j*angle(ph_out)); % reset magnitude
if strcmpi(low_flag,'y')
    ph_out_low=abs(ph).*exp(1j*angle(ph_out_low)); % reset magnitude
end
end
