ps_merge_patches.m

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

function []=ps_merge_patches(psver)
%PS_MERGE_PATCHES merge patches
%
%   Andy Hooper, September 2006
%  
%   ======================================================================
%   01/2007 AH: scla and scn added
%   01/2007 AH: patch.list added
%   ======================================================================

fprintf('Merging patches...\n')

if nargin < 1
  psver=2;
end

small_baseline_flag=getparm('small_baseline_flag');

psname=['ps',num2str(psver)];
phname=['ph',num2str(psver)];
rcname=['rc',num2str(psver)];
pmname=['pm',num2str(psver)];
phuwname=['phuw',num2str(psver)];
sclaname=['scla',num2str(psver)];
sclasbname=['scla_sb',num2str(psver)];
scnname=['scn',num2str(psver)];
bpname=['bp',num2str(psver)];
laname=['la',num2str(psver)];
hgtname=['hgt',num2str(psver)];

if exist('patch.list','file')
    dirname=struct;
    fid=fopen('patch.list','r');
    i=0;
    while feof(fid)==0
        i=i+1;
        dirname(i).name=fgetl(fid);
    end
    fclose(fid);
else
    dirname=dir('PATCH_*');
    %dirname=flipud(dirname);
end

n_patch=length(dirname);
remove_ix=logical([]);
ij=zeros(0,3);
lonlat=zeros(0,2);
ph=zeros(0,0);
ph_rc=zeros(0,0);
ph_reref=zeros(0,0);
ph_uw=zeros(0,0);
ph_patch=zeros(0,0);
ph_res=zeros(0,0);
ph_scla=zeros(0,0);
ph_scla_sb=zeros(0,0);
ph_scn_master=zeros(0,0);
ph_scn_slave=zeros(0,0);
K_ps=zeros(0,0);
C_ps=zeros(0,0);
coh_ps=zeros(0,0);
K_ps_uw=zeros(0,0);
K_ps_uw_sb=zeros(0,0);
C_ps_uw=zeros(0,0);
C_ps_uw_sb=zeros(0,0);
bperp_mat=zeros(0,0);
la=zeros(0,0);
hgt=zeros(0,0);
amp=zeros(0,0);


for i=1:n_patch
    fprintf('   Processing %s\n',dirname(i).name)
    cd(dirname(i).name);
    ps=load(psname);
    n_ifg=ps.n_ifg;
    if isfield('ps','n_image')
        n_image=ps.n_image;
    else
        n_image=ps.n_ifg;
    end
    rc=load(rcname);
    if exist([phname,'.mat'],'file')
        phin=load(phname);
        ph_w=phin.ph;
        clear phin;
    elseif isfield(ps,'ph')
        ph_w=ps.ph;
    end
    pm=load(pmname);
    if exist([phuwname,'.mat'],'file')
        phuw=load(phuwname);
    end
    if exist([sclaname,'.mat'],'file')
        scla=load(sclaname);
    end
    if exist([sclasbname,'.mat'],'file')
        sclasb=load(sclasbname);
    end
    if exist([scnname,'.mat'],'file')
        scn=load(scnname);
    end
    bp=load(bpname);
    
    patch.ij=load('patch_noover.in');
    ix=ps.ij(:,2)>=patch.ij(3)-1 & ps.ij(:,2)<=patch.ij(4)-1 & ps.ij(:,3)>=patch.ij(1)-1 & ps.ij(:,3)<=patch.ij(2)-1;
    [C,IA,IB]=intersect(ps.ij(ix,2:3),ij(:,2:3),'rows'); 
    remove_ix=[remove_ix;IB]; % now more reliable values for these pixels
    [C,IA,IB]=intersect(ps.ij(:,2:3),ij(:,2:3),'rows');
    ix_ex=true(ps.n_ps,1);
    ix_ex(IA)=0; % exclusive index (non-intersecting)
    ix(ix_ex)=1; % keep those in patch proper and those outside not already kept
    

    ij=[ij;ps.ij(ix,:)];
    lonlat=[lonlat;ps.lonlat(ix,:)];
    if exist('ph_w','var')
        ph=[ph;ph_w(ix,:)];
        clear ph_w
    end
    ph_rc=[ph_rc;rc.ph_rc(ix,:)];
    if ~strcmpi(small_baseline_flag,'y')
        ph_reref=[ph_reref;rc.ph_reref(ix,:)];
    end
    
    ph_patch=[ph_patch;pm.ph_patch(ix,:)];
    if isfield(pm,'ph_res')
        ph_res=[ph_res;pm.ph_res(ix,:)];
    end
    if isfield(pm,'K_ps')
        K_ps=[K_ps;pm.K_ps(ix,:)];
    end
    if isfield(pm,'C_ps')
        C_ps=[C_ps;pm.C_ps(ix,:)];
    end
    if isfield(pm,'coh_ps')
        coh_ps=[coh_ps;pm.coh_ps(ix,:)];
    end
    bperp_mat=[bperp_mat;bp.bperp_mat(ix,:)];
    if exist([laname,'.mat'],'file')
        lain=load(laname);
        la=[la;lain.la(ix,:)];
    end

    if exist([hgtname,'.mat'],'file')
        hgtin=load(hgtname);
        hgt=[hgt;hgtin.hgt(ix,:)];
    end


    %overlap_ij=ps.ij(~ix,:);
    if exist([phuwname,'.mat'],'file')
        if ~isempty(C)
            ph_uw_diff=mean(phuw.ph_uw(IA,:)-ph_uw(IB,:),1);
            if ~strcmpi(small_baseline_flag,'y')
                ph_uw_diff=round(ph_uw_diff/2/pi)*2*pi; % round to nearest 2 pi
            end
        else
            ph_uw_diff=zeros(1,size(phuw.ph_uw,2));
        end
        ph_uw=[ph_uw;phuw.ph_uw(ix,:)-repmat(ph_uw_diff,sum(ix),1)];
    else
        ph_uw=[ph_uw;zeros(sum(ix),n_image,'single')];
    end
    
    if exist([sclaname,'.mat'],'file')
        if ~isempty(C)
            ph_scla_diff=mean(scla.ph_scla(IA,:)-ph_scla(IB,:));
            K_ps_diff=mean(scla.K_ps_uw(IA,:)-K_ps_uw(IB,:));
            C_ps_diff=mean(scla.C_ps_uw(IA,:)-C_ps_uw(IB,:));
        else
            ph_scla_diff=zeros(1,size(scla.ph_scla,2));
            K_ps_diff=0;
            C_ps_diff=0;
        end
        ph_scla=[ph_scla;scla.ph_scla(ix,:)-repmat(ph_scla_diff,sum(ix),1)];
        K_ps_uw=[K_ps_uw;scla.K_ps_uw(ix,:)-repmat(K_ps_diff,sum(ix),1)];
        C_ps_uw=[C_ps_uw;scla.C_ps_uw(ix,:)-repmat(C_ps_diff,sum(ix),1)];
    else
        ph_scla=[ph_scla;zeros(sum(ix),n_image,'single')];
        K_ps_uw=[K_ps_uw;zeros(sum(ix),1,'single')];
        C_ps_uw=[C_ps_uw;zeros(sum(ix),1,'single')];
    end

    if exist([sclasbname,'.mat'],'file')
        ph_scla_diff=mean(sclasb.ph_scla(IA,:)-ph_scla_sb(IB,:));
        K_ps_diff=mean(sclasb.K_ps_uw(IA,:)-K_ps_uw_sb(IB,:));
        C_ps_diff=mean(sclasb.C_ps_uw(IA,:)-C_ps_uw_sb(IB,:));
        ph_scla_sb=[ph_scla_sb;sclasb.ph_scla(ix,:)-repmat(ph_scla_diff,sum(ix),1)];
        K_ps_uw_sb=[K_ps_uw_sb;sclasb.K_ps_uw(ix,:)-repmat(K_ps_diff,sum(ix),1)];
        C_ps_uw_sb=[C_ps_uw_sb;sclasb.C_ps_uw(ix,:)-repmat(C_ps_diff,sum(ix),1)];
    else
        ph_scla_sb=[ph_scla_sb;zeros(sum(ix),n_image,'single')];
        K_ps_uw_sb=[K_ps_uw_sb;zeros(sum(ix),1,'single')];
        C_ps_uw_sb=[C_ps_uw_sb;zeros(sum(ix),1,'single')];
    end

    
    
    if exist([scnname,'.mat'],'file')
%        if ~isempty(C)
%            ph_scn_diff=mean(scn.ph_scn_master(IA,:)-ph_scn_master(IB,:));
%        else
%            ph_scn_diff=0;
%        end
%        ph_scn_master=[ph_scn_master;scn.ph_scn_master(ix)-repmat(ph_scn_diff,sum(ix),1)];

        if ~isempty(C)
            ph_scn_diff=mean(scn.ph_scn_slave(IA,:)-ph_scn_slave(IB,:));
        else
            ph_scn_diff=zeros(1,size(scn.ph_scn_slave,2));
        end
        ph_scn_slave=[ph_scn_slave;scn.ph_scn_slave(ix,:)-repmat(ph_scn_diff,sum(ix),1)];
    else
%        ph_scn_master=[ph_scn_master;zeros(sum(ix),1)];
        ph_scn_slave=[ph_scn_slave;zeros(sum(ix),n_image,'single')];
    end

    if exist('./mean_amp.flt','file')    
        pxy=load('patch.in');
        fid=fopen('mean_amp.flt');
        amp(pxy(1):pxy(2),pxy(3):pxy(4))=fread(fid,[pxy(2)-pxy(1)+1,inf],'float=>single');
        fclose(fid);
    end 
    
    cd ..
end
ps_new=ps;
%bperp=ps.bperp;
%day=ps.day;
%master_day=ps.master_day;
%master_ix=ps.master_ix;
n_ps_orig=size(ij,1); % before duplicates removed
keep_ix=true(n_ps_orig,1);
keep_ix(remove_ix)=0;
lonlat_save=lonlat;
coh_ps_weed=coh_ps(keep_ix);
lonlat=lonlat(keep_ix,:);

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

for i=1:length(dups)
    dups_ix_weed=find(lonlat(:,1)==lonlat(dups(i),1)&lonlat(:,2)==lonlat(dups(i),2));
    dups_ix=keep_ix_num(dups_ix_weed);
    [dummy,I]=max(coh_ps_weed(dups_ix_weed));
    keep_ix(dups_ix([1:end]~=I))=0; % drop dups with lowest coh
end

if ~isempty(dups)
    lonlat=lonlat_save(keep_ix,:);
    fprintf('   %d pixel with duplicate lon/lat dropped\n\n',length(dups)')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

ll0=(max(lonlat)+min(lonlat))/2;
xy=llh2local(lonlat',ll0)*1000;
xy=xy';
sort_x=sortrows(xy,1);
sort_y=sortrows(xy,2);
n_pc=round(size(xy,1)*0.001);
bl=mean(sort_x(1:n_pc,:)); % bottom left corner
tr=mean(sort_x(end-n_pc:end,:)); % top right corner
br=mean(sort_y(1:n_pc,:)); % bottom right  corner
tl=mean(sort_y(end-n_pc:end,:)); % top left corner

heading=getparm('heading');
theta=(180-heading)*pi/180;
if theta>pi
    theta=theta-2*pi;
end

rotm=[cos(theta),sin(theta); -sin(theta),cos(theta)];
xy=xy';
xynew=rotm*xy; % rotate so that scene axes approx align with x=0 and y=0
if max(xynew(1,:))-min(xynew(1,:))<max(xy(1,:))-min(xy(1,:)) &...
   max(xynew(2,:))-min(xynew(2,:))<max(xy(2,:))-min(xy(2,:))
    xy=xynew; % check that rotation is an improvement
    disp(['   Rotating xy by ',num2str(theta*180/pi),' degrees']);
end

xy=single(xy');
[dummy,sort_ix]=sortrows(xy,[2,1]); % sort in ascending y order
xy=xy(sort_ix,:);
xy=[[1:size(xy,1)]',xy];
xy(:,2:3)=round(xy(:,2:3)*1000)/1000; % round to mm
lonlat=lonlat(sort_ix,:);

all_ix=[1:size(ij,1)]';
keep_ix=all_ix(keep_ix);
sort_ix=keep_ix(sort_ix);

n_ps=length(sort_ix);
fprintf('   Merged dataset contains %d pixels\n',n_ps)

ij=ij(sort_ix,:);

ph_rc=ph_rc(sort_ix,:);
if ~strcmpi(small_baseline_flag,'y')
    ph_reref=ph_reref(sort_ix,:);
end
save(rcname,'ph_rc','ph_reref');
clear ph_rc ph_reref

ph_uw=ph_uw(sort_ix,:);
save(phuwname,'ph_uw');
clear ph_uw

ph_patch=ph_patch(sort_ix,:);
if size(ph_res,1)==n_ps_orig
    ph_res=ph_res(sort_ix,:);
else
    ph_res=[];
end
if size(K_ps,1)==n_ps_orig
    K_ps=K_ps(sort_ix,:);
else
    K_ps=[];
end
if size(C_ps,1)==n_ps_orig
    C_ps=C_ps(sort_ix,:);
else
    C_ps=[];
end
if size(coh_ps,1)==n_ps_orig
    coh_ps=coh_ps(sort_ix,:);
else
    coh_ps=[];
end
save(pmname,'ph_patch','ph_res','K_ps','C_ps','coh_ps');
clear ph_patch ph_res K_ps C_ps coh_ps

ph_scla=ph_scla(sort_ix,:);
K_ps_uw=K_ps_uw(sort_ix,:);
C_ps_uw=C_ps_uw(sort_ix,:);
save(sclaname,'ph_scla','K_ps_uw','C_ps_uw');
clear ph_scla K_ps_uw C_ps_uw

ph_scla=ph_scla_sb(sort_ix,:);
K_ps_uw=K_ps_uw_sb(sort_ix,:);
C_ps_uw=C_ps_uw_sb(sort_ix,:);
save(sclasbname,'ph_scla','K_ps_uw','C_ps_uw');
clear ph_scla K_ps_uw C_ps_uw ph_scla_sb K_ps_uw_sb C_ps_uw_sb

%ph_scn_master=ph_scn_master(sort_ix,:);
ph_scn_slave=ph_scn_slave(sort_ix,:);
save(scnname,'ph_scn_slave');
clear ph_scn_slave

if size(ph,1)==n_ps_orig
    ph=ph(sort_ix,:);
else
    ph=[];
end
save(phname,'ph');
clear ph


if size(la,1)==n_ps_orig
    la=la(sort_ix,:);
else
    la=[];
end
save(laname,'la');

if size(hgt,1)==n_ps_orig
    hgt=hgt(sort_ix,:);
else
    hgt=[];
end
save(hgtname,'hgt');

fid=fopen('mean_amp.flt','w');
fwrite(fid,amp,'float');
fclose(fid);

bperp_mat=bperp_mat(sort_ix,:);
save(bpname,'bperp_mat');
clear bperp_mat

ps_new.n_ps=n_ps;
ps_new.ij=ij;
ps_new.xy=xy;
ps_new.lonlat=lonlat;
save(psname,'-struct','ps_new');

save psver psver