uw_stat_costs.m

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

function []=uw_stat_costs(unwrap_method,subset_ifg_index);
%UW_STAT_COSTS Find unwrapped solutions using MAP cost functions
%
%   Andy Hooper May 2007
%
%   ======================================================================
%   03/2008 AH: 2D unwrapping option added
%   ======================================================================


tic
if nargin<1
    unwrap_method='3D';
end

costscale=100;
nshortcycle=200;
maxshort=32000;

fprintf('Unwrapping in space...\n')

uw=load('uw_grid');
ui=load('uw_interp');
ut=load('uw_space_time');

if nargin<2
    subset_ifg_index=[1:size(uw.ph,2)];
end

[nrow,ncol]=size(uw.nzix);

[y,x]=find(uw.nzix);
nzix=find(uw.nzix);
z=[1:uw.n_ps];
 

colix=ui.colix;
rowix=ui.rowix;
Z=ui.Z;


grid_edges=[colix(abs(colix)>0);rowix(abs(rowix)>0)];
n_edges=hist(abs(grid_edges),[1:ui.n_edge])';


if strcmpi(unwrap_method,'2D')
      sigsq_defo=0;
      dph_smooth=zeros(size(ut.dph_space));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      n_win=32;
      n_fft=ceil(n_win*1.25);
      clap_alpha=1;
      clap_beta=0.3;
      clap_low_pass_wavelength=800;
      freq0=1/clap_low_pass_wavelength;
      freq_i=-(n_fft)/uw.pix_size/n_fft/2:1/uw.pix_size/n_fft:(n_fft-2)/uw.pix_size/n_fft/2;
      butter_i=1./(1+(freq_i/freq0).^(2*5));
      low_pass=butter_i'*butter_i;
      low_pass=fftshift(low_pass);
      i1=zeros(uw.n_ps,1);
      i2=i1;
      j1=i1;
      j2=i1;
      ps_bit_i=i1;
      ps_bit_j=i1;
      [ps_i,ps_j]=find(uw.nzix); % coords of non-zero grid phase values

      for i=1:uw.n_ps
        i_min=max(ps_i(i)-n_win/2+1,1);
        i_max=min(ps_i(i)+n_win/2,uw.n_i);
        j_min=max(ps_j(i)-n_win/2+1,1);
        j_max=min(ps_j(i)+n_win/2,uw.n_j);
        i1(i)=i_min;
        i2(i)=i_max;
        j1(i)=j_min;
        j2(i)=j_max;
        ps_bit_i(i)=ps_i(i)-i_min+1;
        ps_bit_j(i)=ps_j(i)-j_min+1;
      end

      ifgw_filt=zeros([nrow,ncol,uw.n_ifg],'single');
      ph_var=zeros(uw.n_ps,uw.n_ifg); 
      ifgw=zeros(nrow,ncol,'single');
      x=[0:n_win/2-1];
      [X,Y]=meshgrid(x,x);
      X=X+Y;
      weight_func=[X,fliplr(X)];
      weight_func=[weight_func;flipud(weight_func)]; % pyramidal weighting func
      hw=floor(n_win/2)+1; % half window plus 1
      sa=1; % discontinuity search area in pixels either side
      for i_ifg=1:uw.n_ifg
        ifgw(uw.nzix)=uw.ph(:,i_ifg);
        for i=1:uw.n_ps
          ph_fft=zeros(n_fft,'single');
          ph_bit=ifgw(i1(i):i2(i),j1(i):j2(i));
          ph_fft(1:i2(i)-i1(i)+1,1:j2(i)-j1(i)+1)=ph_bit;
          nz_ix=ph_fft~=0;
          nz_bit_ix=ph_bit~=0;
          n_nz=sum(nz_ix(:));
          if n_nz<=7 
            ph_var(i,i_ifg)=nan;
          else
            ph_filt=clap_filt_patch(ph_fft,clap_alpha,clap_beta,low_pass);
            ph_filt_bit=ph_filt(1:i2(i)-i1(i)+1,1:j2(i)-j1(i)+1);
            wf=weight_func(hw-ps_bit_i(i):i2(i)-i1(i)+hw-ps_bit_i(i),hw-ps_bit_j(i):j2(i)-j1(i)+hw-ps_bit_j(i)); % in case window not full size
            ph_var(i,i_ifg)=sum(angle(ph_filt(nz_ix).*conj(ph_fft(nz_ix))).^2.*wf(nz_bit_ix))/sum(wf(nz_bit_ix))*n_nz/(n_nz-7); % assume 21 degrees of freedom
            ifgw_filt(i1(i):i2(i),j1(i):j2(i),i_ifg)=ifgw_filt(i1(i):i2(i),j1(i):j2(i),i_ifg)+ph_filt_bit.*wf;
          end
        end
        ifgw_hp=angle(ifgw.*conj(ifgw_filt(:,:,i_ifg)));
        ifgw_hp_nz=ifgw_hp(uw.nzix);
        ifgw_hp=reshape(ifgw_hp_nz(Z,i_ifg),nrow,ncol);
        ifgw_hp(ifgw_hp>0)=1;
        ifgw_hp(ifgw_hp<0)=-1;
        col_disc=edge(angle(ifgw_hp),'sobel',0.8,'vertical');
        col_disc(:,1:end-1)=col_disc(:,2:end);
        row_disc=edge(angle(ifgw_hp),'sobel',0.8,'horizontal');
        row_disc(1:end-1,:)=row_disc(2:end,:);
        poss_disc=edge(angle(ifgw_hp),'canny',[0.2],0.4);
       
        fprintf('Stats for IFG %d of %d estimated\n',i_ifg,uw.n_ifg)
      end
      sigsq_noise=zeros(ui.n_edge,uw.n_ifg);
      ps_ij=[ps_i,ps_j];
      for i=1:ui.n_edge
          ij1=ps_ij(ui.edges(i,2),:);
          ij2=ps_ij(ui.edges(i,3),:);
          i1=min([ij1(1),ij2(1)]);
          i2=max([ij1(1),ij2(1)]);
          j1=min([ij1(2),ij2(2)]);
          j2=max([ij1(2),ij2(2)]);
          ph_int_i=squeeze(sum(angle(ifgw_filt(i1+1:i2,j1,:).*conj(ifgw_filt(i1:i2-1,j1,:))),1));
          ph_int_j=squeeze(sum(angle(ifgw_filt(i2,j1+1:j2,:).*conj(ifgw_filt(i2,j1:j2-1,:))),2));
          if ij1(1)<ij2(1)
              ph_int=ph_int_i;
          else
              ph_int=-ph_int_i;
          end
          if ij1(2)<ij2(2)
              dph_smooth(i,:)=(ph_int+ph_int_j).';
          else
              dph_smooth(i,:)=(ph_int-ph_int_j).';
          end
          sigsq_noise(i,:)=(ph_var(ui.edges(i,2),:)+ph_var(ui.edges(i,3),:)/4/pi/pi);
          if round(i/100000)==i/100000 | i==ui.n_edge
             fprintf('Stats for %d of %d edges estimated\n',i,ui.n_edge)
          end
      end
      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else % 3D
  sigsq_noise=(std(ut.dph_noise,0,2)/2/pi).^2;
  sigsq_defo=(std(ut.dph_space_uw-ut.dph_noise,0,2)/2/pi).^2;
  dph_smooth=ut.dph_space_uw-ut.dph_noise;

  nostats_ix=find(isnan(sigsq_noise))'; %noise set to nans in uw_sb_unwrap_time
  for i=nostats_ix
    rowix(abs(rowix)==i)=nan;
    colix(abs(colix)==i)=nan;
  end

  sigsq=int16(round(((sigsq_noise+sigsq_defo)*nshortcycle^2)/costscale.*n_edges)); % weight by number of occurences
  sigsq(sigsq<1)=1; % zero causes snaphu to crash

  rowcost=zeros((nrow-1),ncol*4,'int16');
  colcost=zeros((nrow),(ncol-1)*4,'int16');

  nzrowix=abs(rowix)>0; % 0 = same node, NaN = no stats
  stdgrid=ones(size(rowix),'int16');
  stdgrid(nzrowix)= sigsq(abs(rowix(nzrowix)));
  rowcost(:,2:4:end)= stdgrid; % sigsq

  nzcolix=abs(colix)>0;
  stdgrid=ones(size(colix),'int16');
  stdgrid(nzcolix)= sigsq(abs(colix(nzcolix)));
  colcost(:,2:4:end)= stdgrid; % sigsq

  rowcost(:,3:4:end)= maxshort;% dzmax - value at which cost soars proportional to square of value exceeding dzmax
  colcost(:,3:4:end)= maxshort; % dzmax

  stats_ix=~isnan(rowix); 
  rowcost(:,4:4:end)= int16(stats_ix)*(-1-maxshort)+1; % laycost, -32000 signifies no cost shelf 
  stats_ix=~isnan(colix); 
  colcost(:,4:4:end)= int16(stats_ix)*(-1-maxshort)+1;
end
ph_uw=zeros(uw.n_ps,uw.n_ifg,'single');
ifguw=zeros(nrow,ncol);
for i1=subset_ifg_index
    fprintf('\nProcessing IFG %d of %d\n',i1,length(subset_ifg_index));
    if strcmpi(unwrap_method,'2D')
      nostats_ix=find(isnan(sigsq_noise(:,i1)))'; %noise set to nans in uw_sb_unwrap_time
      for i=nostats_ix
        rowix(abs(rowix)==i)=nan;
        colix(abs(colix)==i)=nan;
      end
      sigsq=((sigsq_noise(:,i1))*nshortcycle^2)/costscale; % don't weight by number of occurences
      sigsq=int16(round(sigsq)); 
      rowcost=zeros((nrow-1),ncol*4,'int16');
      colcost=zeros((nrow),(ncol-1)*4,'int16');
      nzrowix=abs(rowix)>0; % 0 = same node, NaN = no stats
      w_disc=0.5;
      stdgrid=ones(size(rowix),'int16');
      stdgrid(nzrowix)= sigsq(abs(rowix(nzrowix)));
      stdgrid=stdgrid+int16(row_disc(1:end-1,:)*nshortcycle^2*w_disc/costscale);
      rowcost(:,2:4:end)= stdgrid; % sigsq
      nzcolix=abs(colix)>0;
      stdgrid=ones(size(colix),'int16');
      stdgrid(nzcolix)= sigsq(abs(colix(nzcolix)));
      stdgrid=stdgrid+int16(col_disc(:,1:end-1)*nshortcycle^2*w_disc/costscale);
      colcost(:,2:4:end)= stdgrid; % sigsq
      rowcost(:,3:4:end)= maxshort;% dzmax - value at which cost soars proportional to square of value exceeding dzmax
      colcost(:,3:4:end)= maxshort; % dzmax
      stats_ix=~isnan(rowix); 
      rowcost(:,4:4:end)= int16(stats_ix)*(-1-maxshort)+1; % laycost, -32000 signifies no cost shelf 
      stats_ix=~isnan(colix); 
      colcost(:,4:4:end)= int16(stats_ix)*(-1-maxshort)+1;
    end

    offset_cycle=(angle(ut.dph_space(:,i1))-dph_smooth(:,i1))/2/pi;
    offgrid=zeros(size(rowix),'int16');
    offgrid(nzrowix)=round(offset_cycle(abs(rowix(nzrowix))).*sign(rowix(nzrowix))*nshortcycle);
    rowcost(:,1:4:end)= -offgrid; % offset
    offgrid=zeros(size(colix),'int16');
    offgrid(nzcolix)=round(offset_cycle(abs(colix(nzcolix))).*sign(colix(nzcolix))*nshortcycle);
    colcost(:,1:4:end)= offgrid; % offset
    fid=fopen('snaphu.costinfile','w');
    fwrite(fid,rowcost','int16');
    fwrite(fid,colcost','int16');
    fclose(fid);
    ifgw=reshape(uw.ph(Z,i1),nrow,ncol);
    writecpx('snaphu.in',ifgw)
    cmdstr=['!snaphu -d -f $STAMPS/snaphu.conf ',num2str(ncol),' >> snaphu.log'];
    eval(cmdstr)
    fid=fopen('snaphu.out');
    ifguw=fread(fid,[ncol,inf],'float');
    fclose(fid);
    ifguw=ifguw';
    ph_uw(:,i1)=ifguw(uw.nzix);
end


save('uw_phaseuw','ph_uw')