mph.m

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

        scale        = scale_in;
      else
        wrong = 1;
      end;
      phase = angle(mag);
      if (exp~=1) mag = abs(mag).^exp;
      else        mag = abs(mag); end;
    end
  %%% mag/phase/exp/scale/cmap OR m/p/exp/cmap/sc OR mag/phase/cmap/exp/sc
  case 5
    if (isreal(mag) & isreal(phase))
      if (size(exp_in,2)==3)            phasecolors = exp_in;
      elseif (size(scale_in,2)==3)      phasecolors = scale_in;
      elseif (size(cmap_phase_in,2)==3) phasecolors = cmap_phase_in;
      else wrong=1; end;
      if (prod(size(exp_in))==1 & prod(size(scale_in))==1)
        exp          = exp_in;
        scale        = scale_in;
      elseif (prod(size(exp_in))==1 & prod(size(cmap_phase_in))==1)
        exp          = exp_in;
        scale        = cmap_phase_in;
      elseif (prod(size(scale_in))==1 & prod(size(cmap_phase_in))==1)
        exp          = scale_in;
        scale        = cmap_phase_in;
      else
        wrong = 1;
      end;
    else
      wrong = 1;
    end
    if (exp~=1) mag = mag.^exp; end;
  %%% More input than I can handle.
  otherwise
    wrong = 1;
end;
if (wrong~=0) helphelp; break; end;


%%% Internal variables.
numcolorsphase = size(phasecolors,1);
numcolorsmag   = 16;%		or input variable?
if (size(CMAP_MAG_MPH,2)==3)
  disp('Using global variable CMAP_MAG_MPH for magnitude');
  numcolorsmag=size(CMAP_MAG_MPH,1);
end;
disp(['number of colors for phase/mag: ',int2str(numcolorsphase), ...
      ' ',int2str(numcolorsmag), ' (', ...
      int2str(numcolorsphase*numcolorsmag),')']);


%%% After input handling magnitude data is: mag.^exp
%%% Rescale data as: scale*150./mean(mag.^exp)
meanmag = mean(mag(:));% or approximate value for this

% should threshold top 3% or so.
%%% This seems to work fine, fixed threshold, scale influences number
%%% of pixels that are thresholded.
%%mag     = mag./meanmag;%	normalize around one.
%%mag     = 150.*mag;%		around 0:150:1000
%%mag     = scale.*mag;%	around scale*[0:150:1000]
%%mag(find(mag<16))  = 16;%	lower threshold
%%mag(find(mag>255)) = 255;%	and upper threshold
mag     = (150.*scale./meanmag).*mag;% i.e. rescale
q1      = find(mag<16);%	lower threshold
q2      = find(mag>255);%	and upper threshold
mag(q1) = 16;
mag(q2) = 255;
disp(['scale: ',num2str(scale),' --> thresholding bottom: ', ...
       num2str(100.*length(q1)./prod(size(mag))),'%']);
disp(['scale: ',num2str(scale),' --> thresholding top:    ', ...
       num2str(100.*length(q2)./prod(size(mag))),'%']);
%%% Now bin mag in integers [0:numcolorsmag-1].
%mag = floor(mag./(255./numcolorsmag));% i.e. 1:16?
mag = floor(mag./(256./numcolorsmag));% i.e. 1:15

%
if (1==2)% --COMMENTED OUT--
%maxmag  = max(mag(:));%	or an approximate value for max
%minmag  = min(mag(:));%	or an approximate value for min
threshi = meanmag*thresh
%threslo = meanmag/thresh
threslo = 0
q1  =find(mag>threshi);
%q2 =find(mag<threslo);
%length(q1)
disp(['threshold: ',num2str(threshi), ...
       ' --> thresholding top:    ', ...
       num2str(100.*length(q1)./prod(size(mag))),'%']);
%disp(['thresholding bottom: ',num2str(100.*length(q2)./prod(size(mag))),'%']);
mag(q1)=threshi;
%mag(q2)=threslo;
%%% Rescale mag [0:numcolorsmag-1], threslo=black, threshi=white
mag = floor((mag-threslo)./((threshi-threslo)./numcolorsmag));
%max(mag(:))
%min(mag(:))
%mean(mag(:))
%figure; hist(mag(:))
end


%%% Rescale phase to integers [0:ncolorspha-1].
%%% Make sure data is in this interval.
%%% First convert phase to integers 0:15 (floats 0:15.999, floor them)
%%% Be sure minphase and maxphase are the correct values,
%%% otherwise check values later on with find<0, find>max.
maxphase = max(phase(:));%	or an approximate value for max
minphase = min(phase(:));%	or an approximate value for min
phasenan = isnan(phase);
phase(phasenan) = 0;%	restore later, but else NaN -> computation wrong
phase    = floor((phase-minphase)./((maxphase-minphase+0.001)./(numcolorsphase)));
%phase   = int8((phase-minphase) ./((maxphase-minphase)./(numcolorsphase)));
%phase(find(phase>numcolorsphase-1))=numcolorsphase-1;  
%phase(find(phase<0))=0;  

%%% Rescale magnitude to integers [0:ncolmag-1]. (floats 0:ncol, floor this)
%%% Make sure data is in this interval.
%mag = int8((mag-minmag) ./ ((maxmag-minmag)./(numcolorsmag)));
%mag = floor((mag-minmag) ./ ((maxmag-minmag)./(numcolorsmag)));
%mag(find(mag>numcolorsmag-1))=numcolorsmag-1;  
%mag(find(mag<0))=0;  

%%% Now create new data array with correct entries in 
%%% data = 16*mag + phase == e[0:255] (uchar)
%%% Make sure data is in this interval.
%%% colormap [0:255]eN, colorindex = 16*I + H
%%% for sure colorbar lookup exact, otherwise not garantueed
%%% or use: caxis([0 255]);
data      = numcolorsphase.*mag + phase;%	transformation formula
data(1,1) = 0;
data(1,2) = numcolorsphase*numcolorsmag-1;%

%%% Create colormap, entries r,g,b, all e[0,1]
cmap_mixed = [];
cmap_phase = [];
cmap_mag   = [];
whitemap   = ones(size(phasecolors));
%
for magnitudelevel = 0:numcolorsmag-1
  % Level scales colormap RGB for phase (shade)
  % Adapt level, leave all black out.
  %q1    = numcolorsmag/2;% control whiteness min. ampl.
  %q2    = numcolorsmag/4;% control whiteness max. ampl.
  q1         = 0;
  q2         = 0;
  level      = (magnitudelevel+q1)./(numcolorsmag-1+q1+q2);
  cmap_mag   = [cmap_mag;   level.*whitemap];
  cmap_phase = [cmap_phase; phasecolors];%	simply repmat...
  cmap_mixed = [cmap_mixed; level.*phasecolors];
end

%%% Display magnitude where phase is NaN.
%%% Comment this out if undesired.
if (1==1)% do put magnitude where phase has NaN.
if (~isempty(phasenan))
  disp('% Setting phase NaN area to magnitude.');
  %%% change cmap, use lowest level of magnitude for this (excluded: mag>=16)
  if (numcolorsmag>numcolorsphase)
    warning('more mag colors than phasecolors, might give trouble.');
  end
  map_x       = gray(numcolorsmag);
  cmap_mixed(1:numcolorsmag,:) = map_x;%	extended with gray for mag.
  mag_scaled  = floor(data(phasenan)./numcolorsphase);% get rescaled amplitude [0:15];
  data(phasenan) = mag_scaled;%                        [-1,-16]
end
else
  phase(phasenan) = NaN;%restore later, but else NaN -> computation wrong
end

%%% Output.
if (nargout==0)
  %%% Plot for now
  %figure(1);
  %  imagesc(phase); colormap(phasecolors); colorbar
  %figure(2);
  %  imagesc(mag); colormap(gray(numcolorsmag)); colorbar
  figure(10011);
    imagesc(data,[0 numcolorsmag*numcolorsphase-1]);
    colormap(cmap_mixed);
    colorbar;
    title('rescaled data: mix magnitude,phase');
    if (exist('pixval'))
      set(10011,'units','pixels')
      pixval;
    end; %images toolbox
    %caxis([0 size(cmap_mixed,1)-1]);
    %caxis([0 numcolorsmag*numcolorsphase-1]);
  figure(10012);
    imagesc(data); colormap(cmap_phase); colorbar
    title('rescaled data: phase');
    if (exist('pixval'))
      set(10012,'units','pixels')
      pixval;
    end; %images toolbox
  figure(10013);
    imagesc(data); colormap(cmap_mag); colorbar
    title('rescaled data: magnitude');
    if (exist('pixval'))
      set(10013,'units','pixels')
      pixval;
    end; %images toolbox
  clear data cmap_mixed cmap_mag cmap_phase
end;
if (nargout==1)
  clear cmap_mixed cmap_mag cmap_phase
end
if (nargout==2)
  clear cmap_mag cmap_phase
end
if (nargout==3)
  clear cmap_phase
end

%%% EOF.