siminterf.m

荷兰Delft大学开发的insar（干涉合成孔径雷达）图像处理部分源代码

disp(['First range bin at:     ', num2str(round(rangegrid(1))./1e3), ' km']);
disp(['Last range bin at       ', num2str(round(rangegrid(numpixels))./1e3), ' km']);


%%% Compute range diff to slave satellite
%%% Assume range to first bin of DEM = R1
%%% use local coord. sytem to compute range to master
%satH = cos(theta) .* R1 + DEM(1,1);
B      = Bperp ./ cos(theta-alpha);
sat2_x = B .* cos(alpha);
sat2_y = B .* sin(alpha) + sat1_y;
x      = x0:dx:x0+dx*(numpixelsdem-1);%	x coord. w.r.t. sat1
x2sqr  = (x - sat2_x).^2;
xsqr   = x.^2;


%%% Compute range for whole line in 1 go.
% and compute interferometric phase, flat earth corrected.
% refphase may be wrong, but trend should be ok this way.
%%%range2master = sqrt(sat1_y.^2 + xsqr);
%%%range2slave  = sqrt(sat2_y.^2 + x2sqr);
%%%refphase     = pi4divlam .* (range2slave-range2master);
if ( dorefpha == 0 )
%%%  %refphase = zeros(size(refphase));
%%%  refphase = 0.;
  disp('Reference phase:        not subtracted')
else
  disp('Reference phase:        subtracted')
end


%%% Compute range for all lines, same refphase
oldperc = -1;
for linenum=1:size(DEM,1)
  newperc = floor((100*linenum)./numlines);
  if (newperc ~= oldperc)
    oldperc = newperc;
    %comment out on x86 PCs where '\r' does not work...
    fprintf(1,'\rRadarcoding DEM:        %3.0f%%', newperc);
  end
  y            = sat1_y-DEM(linenum,:);
  range2master = sqrt(y.^2+xsqr);
  y2           = sat2_y-DEM(linenum,:);
  range2slave  = sqrt(y2.^2+x2sqr);
  phase = pi4divlam .* (range2slave-range2master);

  if ( dorefpha ~= 0 )% 			subtract it
    tantheta = x./y2;
    deltax   = DEM(linenum,:) ./ tantheta;%	far field approx.
    x2_0     = x - deltax;
    %refpharangemaster = range2master; (approx...)
    refpharangemaster = sqrt(sat1_y.^2 +  x2_0.^2);
    refpharangeslave  = sqrt(sat2_y.^2 + (x2_0-sat2_x).^2);
    refphase = pi4divlam .* (refpharangeslave-refpharangemaster);

    phase = phase - refphase;
%not ok.
%    phase = pi4divlam .* (refpharange2slave-range2slave);
%  else
%    phase = pi4divlam .* (range2slave-range2master);
  end

  %%% Interpolate p to grid rangebins
  %%% range is not always increasing due to foreshortning
  [range2master,sortindex] = sort(range2master);
  phase                    = phase(sortindex);         
  interf_nonoise(linenum,:) = interp1(range2master,phase,rangegrid,method);

  if ( donoise ~= 0 )
    %%% Method must be nearest! (why?)
    slopeDEM = atan2((diff(DEM(linenum,:),1,2)),dx);% 1 kleiner...
    slopeDEM = [slopeDEM,0];
    slopeDEM = slopeDEM(sortindex);
    slope(linenum,:) = interp1(range2master,slopeDEM,rangegrid,'nearest');
  end
end
disp(' ');
disp(' ');

%%% Get watermask for radarcoded matrices.
watermask = [];
if ( fracdim ~= 999 )		% input matrix
  watermask = find(interf_nonoise<=10*eps);
end

%%% Adding noise
if ( donoise ~= 0 )
  disp('Adding to phase:        geometric decorrelation noise')
  [noise,coherence]  = simnoise(slope,Bperp);% rad
else
  noise     = zeros(size(interf_nonoise));
  coherence = zeros(size(interf_nonoise));
end
interferogram = interf_nonoise + noise;
%%% Tidy up a little.
clear slope;


%%% Adjusting the coherence of watermask
%%% Coherence of watermask is calculated wrong using the formula :
%%% Bcritical = lambda*(Bw/c)*R1*tan(theta-slope)
%%% since watermask doen't have any slope, Bcritical is constant
%%% watermask is given random coherence between 0.0 and 0.2
disp('Water area:             Coherence [0.0,0.2]');
coherence(watermask)=rand(size(watermask))*0.2;


if ( donoise ~= 0 )
  disp('Water area:             Phase uniform noise.');
  interferogram(watermask)=rand(size(watermask))*2*pi-pi;
else
  disp('Water area:             Phase = 0.');
  interferogram(watermask)=0;
end


%%% Plot some.
if ( doplot ~= 0 )
disp(' ');
disp('Plotting:               DEM (3D)');
%j       = jet(64);
%topomap = [(j(1,:));j(64:-1:15,:)];
%phasemap = ph(256);
j       = flipud(jet(256));
topomap = [0,0,0.5625;j(1:200,:)];
phasemap = deos(256);
figure(1);
  clf;
  %h=surf(DEM);
  %set (h,'edgecolor','none')
  mesh(DEM);
  colormap(topomap);
  l=line(round([-0.1*numpixels -0.1*numpixels]), ...
         round([-0.3*numlines 1.2*numlines]) ,...
         round([1.2*maxheight 1.2*maxheight]));
  set(l,'LineWidth',3);
  %set(l,'Linestyle','--');
  set(l,'color',[1 0 0]);
  set(l,'markersize',6);
  set(l,'marker','<');
  text(round(-0.25*numpixels),round(numlines),round(1.2*maxheight),'orbit')
  axis ('tight');
  set(gca,'Zlim',[0 1.6*maxheight]);
  %
  %view(20,65);
  %view(-15,40);
  %view(10,50);
  view(25,60);
  t=title (['DEM (height [0:',num2str(maxheight),'], ',num2str(water),'% water)']);
  set(t,'fontweight','bold');
  xlabel ('range pixels');
  ylabel ('azimuth lines');
  zlabel ('height');
  c=colorbar;
  set(get(c,'title'),'string','[meter]');


disp('Plotting:               DEM (2D)');
figure(2);
  clf;
  imagesc(DEM);
  t=title ('Digital Elevation Model');
  set(t,'fontweight','bold');
  xlabel ('range pixels');
  ylabel ('azimuth lines');
  zlabel ('height');
  axis image;
  colormap(topomap);
  c=colorbar;
  set(get(c,'title'),'string','[meter]');

disp('Plotting:               Unwrapped radarcoded DEM w/o noise');
figure(3);
  clf;
  mesh(interf_nonoise);
  axis ('tight');
  set(gca,'Zlim',[0 1.6*max(interf_nonoise(:))]);
  colormap(phasemap);
  view(25,60);
  %g = surf(interf_nonoise);
  %set (g,'edgecolor','none')
  %t=title ('Radarcoded DEM (B\perp=',num2str(Bperp),')');
  t=title (['Radarcoded DEM (Bperp=',num2str(Bperp),')']);
  set(t,'fontweight','bold');
  xlabel ('range pixels');
  ylabel ('azimuth lines');
  zlabel ('height');
  %axis image;
  colormap(topomap);
  c=colorbar;
  set(get(c,'title'),'string','[rad]');

disp('Plotting:               Unwrapped radarcoded DEM with noise');
figure(4);
  clf;
  imagesc(interferogram);
  colormap(phasemap);
  axis image;
  c=colorbar;
  set(get(c,'title'),'string','[rad]');
  t=title (['Phase of radarcoded DEM (Bperp=',num2str(Bperp),')']);
  set(t,'fontweight','bold');
  xlabel ('range pixels');
  ylabel ('azimuth lines');

disp('Plotting:               Wrapped radarcoded DEM with noise');
figure(5);
  clf;
  imagesc(wrap(interferogram));
  colormap(phasemap);
  axis image;
  c=colorbar;
  set(get(c,'title'),'string','[rad]');
  %t=title (['Wrapped phase (B\perp =',num2str(Bperp),')']);
  t=title (['Wrapped phase (Bperp =',num2str(Bperp),')']);
  set(t,'fontweight','bold');
  xlabel ('range pixels');
  ylabel ('azimuth lines');

disp('Plotting:               Coherence map');
figure(6);
  clf;
  imagesc(coherence,[0 1]);
  colormap(gray);
  t=title ('Coherence');
  set(t,'fontweight','bold');
  xlabel ('range pixels');
  ylabel ('azimuth lines');
  axis image;
  c=colorbar;
  set(get(c,'title'),'string','[-]');

% make figures active
if (exist('tip','file'))
  tip;
else
  for ii=6:-1:1
    figure(ii);
  end;
end;
end

%%% EOF
if (nargout ~= 0) 
  interf=[];
  [interf, interferogram] = deal(interferogram,interf);
end
  
more on;