ps_weed.m

来自「StaMps最新测试版」· M 代码 · 共 477 行

function []=ps_weed(all_da_flag,no_weed_adjacent,no_weed_noisy)
%PS_WEED weeds out neighboring PS and save those kept to new version
%   PS_weed(all_da_flag,no_weed_adjacent,no_weed_noisy)
%
%   Andy Hooper, June 2006
%
%   ==========================================================
%   09/2006 AH: create all workspace files directly
%   09/2006 AH: drop noisy pixels
%   09/2006 AH: add small baselines 
%   01/2007 AH: drop pixels with duplicate lon/lat 
%   05/2007 AH: optional weeding of pixels with zero elevation added 
%   ==========================================================

fprintf('Weeding selected pixels...\n')

if nargin<1
    all_da_flag=0;
end

weed_alpha=getparm('weed_alpha',1);
weed_standard_dev=getparm('weed_standard_dev',1);
weed_zero_elevation=getparm('weed_zero_elevation',1);
small_baseline_flag=getparm('small_baseline_flag',1);

if nargin<2
    no_weed_adjacent=0;
end
if nargin<3
    if weed_standard_dev>=10
        no_weed_noisy=1;
    else
        no_weed_noisy=0;
    end
end

load psver

psname=['ps',num2str(psver)];
pmname=['pm',num2str(psver)];
phname=['ph',num2str(psver)];
selectname=['select',num2str(psver)];
hgtname=['hgt',num2str(psver),'.mat'];
laname=['la',num2str(psver),'.mat'];
bpname=['bp',num2str(psver),'.mat'];

psothername=['ps_other'];
psothername=['pm_other'];
selectothername=['select_other'];
hgtothername=['hgt_other'];
laothername=['la_other'];
bpothername=['bp_other'];

ps=load(psname);
sl=load(selectname);

if exist([phname,'.mat'],'file')
    phin=load(phname);
    ph=phin.ph;
    clear phin
else
    ph=ps.ph;
end

day=ps.day;
bperp=ps.bperp;
master_day=ps.master_day;

if isfield(sl,'keep_ix')
    ix2=sl.ix(sl.keep_ix);
    K_ps2=sl.K_ps2(sl.keep_ix);
    C_ps2=sl.C_ps2(sl.keep_ix);
    coh_ps2=sl.coh_ps2(sl.keep_ix);
else
    ix2=sl.ix2;
    K_ps2=sl.K_ps2;
    C_ps2=sl.C_ps2;
    coh_ps2=sl.coh_ps2;
end

ij2=ps.ij(ix2,:);
xy2=ps.xy(ix2,:);
ph2=ph(ix2,:);
lonlat2=ps.lonlat(ix2,:);

pm=load(pmname);
ph_patch2=pm.ph_patch(ix2,:); % use original patch phase, with PS left in
if isfield(sl,'ph_res2')
    ph_res2=sl.ph_res2(sl.keep_ix,:);
else
    ph_res2=[];
end
%snr2=pm.snr(sl.ix2,:);
clear pm

clear sl

clear ph
if isfield(ps,'ph')
    ps=rmfield(ps,'ph');
end
ps=rmfield(ps,{'xy','ij','lonlat','sort_ix'});

if all_da_flag~=0
    pso=load(psothername);
    slo=load(selectothername);
    ix_other=slo.ix_other;
    n_ps_other=sum(ix_other);
    K_ps_other2=pso.K_ps_other(ix_other);
    C_ps_other2=pso.C_ps_other(ix_other);
    coh_ps_other2=pso.coh_ps_other(ix_other);
    %ph_patch_other2=pso.ph_patch_other(ix_other,:);
    ph_res_other2=pso.ph_res_other(ix_other,:);
    ij2=[ij2;pso.ij_other(ix_other,:)];
    xy2=[xy2;pso.xy_other(ix_other,:)];
    ph2=[ph2;pso.ph_other(ix_other,:)];
    lonlat2=[lonlat2;pso.lonlat_other(ix_other,:)];
    clear pso slo

    pmo=load(pmothername);
    ph_patch_other2=pmo.ph_patch_other(ix_other,:);
    %snr_other2=pmo.snr_other(ix_other);
    clear pm

    K_ps2=[K_ps2;K_ps_other2];
    C_ps2=[C_ps2;C_ps_other2];
    coh_ps2=[coh_ps2;coh_ps_other2];
    ph_patch2=[ph_patch2;ph_patch_other2];
    %snr2=[snr2;snr_other2];
    ph_res2=[ph_res2;ph_res_other2];
else
    n_ps_other=0;
end

if exist(hgtname,'file')
    ht=load(hgtname);
    hgt=ht.hgt(ix2);
    clear ht
    if all_da_flag~=0
        hto=load(hgtothername);
        hgt=[hgt;hto.hgt_other(ix_other)];
        clear hto
    end
end


n_ps_low_D_A=length(ix2);
n_ps=n_ps_low_D_A + n_ps_other;
ix_weed=logical(ones(n_ps,1));
disp([num2str(n_ps_low_D_A),' low D_A PS, ',num2str(n_ps_other),' high D_A PS']);




if no_weed_adjacent==0
	step_name='INITIALISE NEIGHBOUR MATRIX'
    
    ij_shift=ij2(:,2:3)+repmat([2,2]-min(ij2(:,2:3)),n_ps,1);
	%neigh_ix=sparse(max(ij_shift(:,1))+1,max(ij_shift(:,2))+1);
	neigh_ix=zeros(max(ij_shift(:,1))+1,max(ij_shift(:,2))+1);
    miss_middle=logical(ones(3));
    miss_middle(2,2)=0;
    
	for i=1:n_ps
        neigh_this=neigh_ix(ij_shift(i,1)-1:ij_shift(i,1)+1,ij_shift(i,2)-1:ij_shift(i,2)+1);
        neigh_this(neigh_this==0&miss_middle)=i;
        neigh_ix(ij_shift(i,1)-1:ij_shift(i,1)+1,ij_shift(i,2)-1:ij_shift(i,2)+1)=neigh_this;
        
        if i/100000==floor(i/100000)
            disp([num2str(i),' PS processed'])
            save log step_name i
            !sync
        end
    end
    
	step_name='FIND NEIGHBOURS'
    
    neigh_ps=cell(n_ps,1);
	for i=1:n_ps
        my_neigh_ix=neigh_ix(ij_shift(i,1),ij_shift(i,2));
        if my_neigh_ix~=0
            neigh_ps{my_neigh_ix}=[neigh_ps{my_neigh_ix},i];
        end    
        if i/100000==floor(i/100000)
            disp([num2str(i),' PS processed'])
            save log step_name i
            !sync
        end

    end	
    
    clear neigh_ix
    
	step_name='SELECT BEST'
	
	for i=1:n_ps
        if ~isempty(neigh_ps{i})
            same_ps=i;
            i2=1;
            while i2 <= length(same_ps)
                ps_i=same_ps(i2);
                same_ps=[same_ps,neigh_ps{ps_i}];
                neigh_ps{ps_i}=[];
                i2=i2+1;
            end
            same_ps=unique(same_ps);
            [dummy,high_coh]=max(coh_ps2(same_ps));
            low_coh_ix=logical(ones(size(same_ps)));
            low_coh_ix(high_coh)=0;
            %ix(same_ps(low_coh_ix))=nan;
            ix_weed(same_ps(low_coh_ix))=0;
        end
        if i/100000==floor(i/100000)
            disp([num2str(i),' PS processed'])
            save log step_name i
            !sync
        end
	
	end
	


    disp([num2str(sum(ix_weed)),' PS kept after dropping adjacent pixels']);

end

    if strcmpi(weed_zero_elevation,'y') & exist('hgt','var')
        sea_ix=hgt==0;
        ix_weed(sea_ix)=false;
    end


    xy_weed=xy2(ix_weed,:);

    %%%%%%%%% Some non-adjacent pixels are allocated the same lon/lat by DORIS. %%%%%%%%%
    %%%%%%%%% If duplicates occur, the pixel with the highest coherence is kept.%%%%%%%%% 
    ix_weed_num=find(ix_weed); 
    [dummy,I]=unique(xy_weed(:,2:3),'rows');
    dups=setxor(I,[1:sum(ix_weed)]'); % pixels with duplicate lon/lat

    for i=1:length(dups)
        dups_ix_weed=find(xy_weed(:,2)==xy_weed(dups(i),2)&xy_weed(:,3)==xy_weed(dups(i),3));
        dups_ix=ix_weed_num(dups_ix_weed);
        [dummy,I]=max(coh_ps2(dups_ix));
        ix_weed(dups_ix([1:end]~=I))=0; % drop dups with lowest coh
    end

    if ~isempty(dups)
        xy_weed=xy2(ix_weed,:);
        fprintf('%d PS with duplicate lon/lat dropped\n\n',length(dups)')
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

n_ps=sum(ix_weed);
ps_std=zeros(n_ps,1);

if no_weed_noisy==0
    
    pack
    step_name='DROP NOISY'

    nodename=['psweed.1.node'];
    fid=fopen(nodename,'w');
    fprintf(fid,'%d 2 0 0\n',n_ps);

    for i=1:n_ps
        fprintf(fid,'%d %f %f\n',i,xy_weed(i,2),xy_weed(i,3));
    end

    fclose(fid);

    !triangle -e psweed.1.node

    fid=fopen('psweed.2.edge','r');
    header=str2num(fgetl(fid));
    N=header(1);
    edges=zeros(N,4);
    for i=1:N
        edges(i,:)=str2num(fgetl(fid));
    end
    fclose(fid);
    n_edge=size(edges,1);
    ph_weed=ph2(ix_weed,:).*exp(-j*(K_ps2(ix_weed)*bperp'));  % subtract range error 
    edge_std=zeros(n_edge,1);
    dph_space=(ph_weed(edges(:,3),:).*conj(ph_weed(edges(:,2),:)));

    if ~strcmpi(small_baseline_flag,'y')

        n_pad=2^ceil(log(ps.n_ifg*2)/log(2));
        dph_padded=zeros(1,n_pad);
        smooth_win=gausswin(n_pad,weed_alpha)';

        for i=1:n_edge

            dph_padded(1:ps.n_ifg)=dph_space(i,:);
            dph_fft=fftshift(fft(dph_padded));
            dph_smooth=ifft(ifftshift(dph_fft.*smooth_win));
            edge_std(i)=std(angle(dph_space(i,:).*conj(dph_smooth(1:ps.n_ifg))));

            if i/10000==floor(i/10000)
                fprintf('   %d edges of %d processed\n',i,n_edge)
            end
        end

    else
         edge_std=std(angle(dph_space),0,2);
    %     G=zeros(ps.n_ifg,ps.n_image-1);
    %     for i=1:ps.n_ifg
    %         G(i,ps.ifgday_ix(i,1):ps.ifgday_ix(i,2)-1)=1;
    %     end
    %     
    %     n_pad=2^ceil(log(ps.n_image*2)/log(2));
    %     dph_padded=zeros(1,n_pad);
    %     smooth_win=gausswin(n_pad,weed_alpha)';
    % 
    %     dph_space=zeros(n_edge,ps.n_ifg);
    % 
    %     for i=1:n_edge
    %         dph_space(i,:)=(ph_weed(edges(i,3),:).*conj(ph_weed(edges(i,2),:)));
    %         dph_space(i,:)=dph_space(i,:)./abs(dph_space(i,:));
    %         dph_uw=angle(dph_space(i,:));
    %         if dph_uw==0
    %             loop_sw=0;
    %             dph_space_step_time=zeros(1,ps.n_image-1);
    %         else
    %             loop_sw=1;
    %         end
    %         ii=1;
    %         max_diff=0.7*pi;
    %         r=zeros(1000,1);
    %         dph_uw_ok=[];
    %         
    %         while loop_sw==1
    % 
    %             dph_space_step_time=(G\dph_uw')';
    %             dph_hat=(G*dph_space_step_time')';
    %             dph_r=dph_hat-dph_uw;
    %             if max(abs(dph_r))>max_diff & ii<100
    %                 %if max(abs(dph_r))<pi
    %                 %    dph_uw_ok=dph_uw;
    %                 %end
    %                 [Y,I]=max(abs(dph_r));
    %                 dph_uw(I)=dph_uw(I)+2*pi*sign(dph_r(I));
    %                 r(ii)=sum(dph_uw).^2;
    %                 if r(ii)<min(r(1:ii-1))
    %                     dph_uw_ok=dph_uw;
    %                 end
    % 
    %                 if sum(r(ii)-r(1:ii-1)==0)>0 % been here before
    %                     dph_uw(I)=dph_uw(I)-2*pi*sign(dph_r(I));
    %                     dph_r(I)=0;
    %                     [Y,I]=max(abs(dph_r));
    %                     dph_uw(I)=dph_uw(I)+2*pi*sign(dph_r(I));
    %                     r(ii)=sum(dph_uw).^2;
    %                 end
    %                 if sum(r(ii)-r(1:ii-1)==0)>0 % been here too many times
    %                     if ~isempty(dph_uw_ok)
    %                         dph_uw=dph_uw_ok; % previously found OK solution
    %                     end
    %                     loop_sw=0;
    %                 end
    %                 ii=ii+1;
    %             else
    %                 loop_sw=0;
    %             end
    %         end
    %         dph_padded(1:ps.n_image)=exp(j*cumsum([0,dph_space_step_time]));
    %         dph_fft=fftshift(fft(dph_padded));
    %         dph_smooth_image=ifft(ifftshift(dph_fft.*smooth_win));
    %         dph_smooth=exp(j*(G*angle(dph_smooth_image(2:ps.n_image).*conj(dph_smooth_image(1:ps.n_image-1)))')');
    %         edge_std(i)=std(angle(dph_space(i,:).*conj(dph_smooth)));
    %      
    %         if i/10000==floor(i/10000)
    %             fprintf('   %d edges of %d processed\n',i,n_edge)
    %         end
    % 
    %     end  
    end



    %nodeedge_ix=sortrows([edges(:,[2,1]);edges(:,[3,1]]);

    for i=1:n_ps
      edge_ix=[find(edges(:,2)==i);find(edges(:,3)==i)];
      ps_std(i)=min(edge_std(edge_ix)); % least noisy 
      if i/10000==floor(i/10000)
          fprintf('   %d PS of %d processed\n',i,n_ps)
      end
    end

    ix_weed(ix_weed)=ps_std<weed_standard_dev;
    n_ps=sum(ix_weed);

    disp([num2str(n_ps),' PS kept after dropping noisy pixels']);

end

weedname=['weed',num2str(psver)];
save(weedname,'ix_weed','ps_std')


coh_ps=coh_ps2(ix_weed);
K_ps=K_ps2(ix_weed);
C_ps=C_ps2(ix_weed);
ph_patch=ph_patch2(ix_weed,:);
if ~isempty(ph_res2)
    ph_res=ph_res2(ix_weed,:);
else
    ph_res=ph_res2;
end

pmname=['pm',num2str(psver+1)];
save(pmname,'ph_patch','ph_res','coh_ps','K_ps','C_ps')
clear ph_patch ph_res coh_ps K_ps C_ps ph_patch2 ph_res2 coh_ps2 K_ps2 C_ps2

ph2=ph2(ix_weed,:);
ph=ph2;
phname=['ph',num2str(psver+1)];
save(phname,'ph')
clear ph

xy2=xy2(ix_weed,:);
ij2=ij2(ix_weed,:);
lonlat2=lonlat2(ix_weed,:);
ps.xy=xy2;
ps.ij=ij2;
ps.lonlat=lonlat2;
ps.n_ps=size(ph2,1);
psname=['ps',num2str(psver+1)];
save(psname,'-struct','ps');
clear ps xy2 ij2 lonlat2

if exist(hgtname,'file')
    hgt=hgt(ix_weed);
    save(['hgt',num2str(psver+1),'.mat'],'hgt');
    clear hgt
end

if exist(laname,'file')
    la=load(laname);
    la=[la.la(ix2)];
    if all_da_flag~=0
        lao=load(laothername);
        la=[la;lao.la_other(ix_other)];
        clear lao
    end
    la=la(ix_weed);
    save(['la',num2str(psver+1),'.mat'],'la');
    clear la
end

if exist(bpname,'file')
    bp=load(bpname);
    bperp_mat=[bp.bperp_mat(ix2,:)];
    clear bp
    if all_da_flag~=0
        bpo=load(bpothername);
        bperp_mat=[bperp_mat;bpo.bperp_other(ix_other,:)];
        clear bpo
    end
    bperp_mat=bperp_mat(ix_weed,:);
    save(['bp',num2str(psver+1),'.mat'],'bperp_mat');
end

if exist(['scla',num2str(psver+1),'.mat'],'file')
    delete(['scla',num2str(psver+1),'.mat'])
end
if exist(['aps',num2str(psver+1),'.mat'],'file')
    delete(['aps',num2str(psver+1),'.mat'])
end
if exist(['scn',num2str(psver+1),'.mat'],'file')
    delete(['scn',num2str(psver+1),'.mat'])
end

setpsver(psver+1)