resize_refl.m

大气模型计算

  error('M must be either a scalar multiplier or a 1-by-2 size vector.');
end

if nargin<4,
  nn = 11; h = []; % Default filter size
else
  if length(arg4)==1, nn = h; h = []; else nn = 0; h = arg4; end
end

[m,n] = size(arg1);

if nn>0 & case0(1)=='b',  % Design anti-aliasing filter if necessary
  if bsize(1)<m, h1 = fir1(nn-1,bsize(1)/m); else h1 = 1; end
  if bsize(2)<n, h2 = fir1(nn-1,bsize(2)/n); else h2 = 1; end
  if length(h1)>1 | length(h2)>1, h = h1'*h2; else h = []; end
end

if ~isempty(h), % Anti-alias filter A before interpolation
  a = filter2(h,arg1);
else
  a = arg1;
end

if case0(1)=='n', % Nearest neighbor interpolation
  dx = n/bsize(2); dy = m/bsize(1);
  uu = (dx/2+.5):dx:n+.5; vv = (dy/2+.5):dy:m+.5;
elseif all(case0 == 'bil') | all(case0 == 'bic'),
  uu = 1:(n-1)/(bsize(2)-1):n; vv = 1:(m-1)/(bsize(1)-1):m;
else
  error(['Unknown interpolation method: ',method]);
end

%
% Interpolate in blocks
%
nu = length(uu); nv = length(vv);
blk = bestblk([nv nu]);
nblks = floor([nv nu]./blk); nrem = [nv nu] - nblks.*blk;
mblocks = nblks(1); nblocks = nblks(2);
mb = blk(1); nb = blk(2);

rows = 1:blk(1); b = zeros(nv,nu);
for i=0:mblocks,
  if i==mblocks, rows = (1:nrem(1)); end
  for j=0:nblocks,
    if j==0, cols = 1:blk(2); elseif j==nblocks, cols=(1:nrem(2)); end
    if ~isempty(rows) & ~isempty(cols)
      [u,v] = meshgrid(uu(j*nb+cols),vv(i*mb+rows));
      % Interpolate points
      if case0(1) == 'n', % Nearest neighbor interpolation
        b(i*mb+rows,j*nb+cols) = interp2(a,u,v,'nearest');
      elseif all(case0 == 'bil'), % Bilinear interpolation
         b(i*mb+rows,j*nb+cols) = interp2(a,u,v,'linear');
     elseif all(case0 == 'bic'), % Bicubic interpolation
        b(i*mb+rows,j*nb+cols) = interp2(a,u,v,'cubic');
      end
    end
  end
end

if nargout==0,
  if isgray(b), imshow(b,size(colormap,1)), else imshow(b), end
  return
end

if isgray(arg1), rout = max(0,min(b,1)); else rout = b; end


%*************************************************************************
%*************************************************************************
%*************************************************************************
%
%  Following functions ISGRAY, BESTBLK & FIR1 are called by
%  IMRESIZE.  They have been made local to RESIZEM
%  since they have been taken from other toolboxes for the
%  sole purpose of making RESIZEM work in the Mapping Toolbox.
%
%*************************************************************************
%*************************************************************************
%*************************************************************************


function y = isgray(x)
%ISGRAY True for intensity images.
%	ISGRAY(A) returns 1 if A is an intensity image and 0 otherwise.
%	An intensity image contains values between 0.0 and 1.0.
%
%	See also ISIND, ISBW.

%  Written by: Clay M. Thompson 2-25-93

y = min(min(x))>=0 & max(max(x))<=1;


%*************************************************************************
%*************************************************************************
%*************************************************************************


function [mb,nb] = bestblk(siz,k)
%BESTBLK Best block size for block processing.
%	BLK = BESTBLK([M N],K) returns the 1-by-2 block size BLK
%	closest to but smaller than K-by-K for block processing.
%
%	[MB,NB] = BESTBLK([M N],K) returns the best block size
%	as the two scalars MB and NB.
%
%	[...] = BESTBLK([M N]) returns the best block size smaller
%	than 100-by-100.
%
%	BESTBLK returns the M or N when they are already smaller
%	than K.
%
%	See also BLKPROC, SIZE.

%  Written by: Clay M. Thompson

if nargin==1, k = 100; end % Default block size

%
% Find possible factors of siz that make good blocks
%

% Define acceptable block sizes
m = floor(k):-1:floor(min(ceil(siz(1)/10),k/2));
n = floor(k):-1:floor(min(ceil(siz(2)/10),k/2));

% Choose that largest acceptable block that has the minimum padding.
[dum,ndx] = min(ceil(siz(1)./m).*m-siz(1)); blk(1) = m(ndx);
[dum,ndx] = min(ceil(siz(2)./n).*n-siz(2)); blk(2) = n(ndx);

if nargout==2,
  mb = blk(1); nb = blk(2);
else
  mb = blk;
end


%*************************************************************************
%*************************************************************************
%*************************************************************************


function [b,a] = fir1(N,Wn,Ftype,Wind)
%FIR1	FIR filter design using the window method.
%	B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter
%	and returns the filter coefficients in length N+1 vector B.
%	The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0
%	corresponding to half the sample rate.
%
%	If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an
%	order N bandpass filter with passband  W1 < W < W2.
%	B = FIR1(N,Wn,'high') designs a highpass filter.
%	B = FIR1(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].
%	For highpass and bandstop filters, N must be even.
%
%	By default FIR1 uses a Hamming window.  Other available windows,
%	including Boxcar, Hanning, Bartlett, Blackman, Kaiser and Chebwin
%	can be specified with an optional trailing argument.  For example,
%	B = FIR1(N,Wn,bartlett(N+1)) uses a Bartlett window.
%	B = FIR1(N,Wn,'high',chebwin(N+1,R)) uses a Chebyshev window.
%
%	FIR1 is an M-file implementation of program 5.2 in the IEEE
%	Programs for Digital Signal Processing tape.  See also FIR2,
%	FIRLS, REMEZ, BUTTER, CHEBY1, CHEBY2, YULEWALK, FREQZ and FILTER.

%  Written by: L. Shure

%	Reference(s):
%	  [1] "Programs for Digital Signal Processing", IEEE Press
%	      John Wiley & Sons, 1979, pg. 5.2-1.

nw = 0;
a = 1;
if nargin == 3
	if ~isstr(Ftype);
		nw = max(size(Ftype));
      		Wind = Ftype;
		Ftype = [];
	end
elseif nargin == 4
   nw = max(size(Wind));
else
   Ftype = [];
end

Btype = 1;
if nargin > 2 & max(size(Ftype)) > 0
	Btype = 3;
end
if max(size(Wn)) == 2
	Btype = Btype + 1;
end

N = N + 1;
odd = rem(N, 2);
if (Btype == 3 | Btype == 4)
   if (~odd)
      disp('For highpass and bandstop filters, order must be even.')
      disp('Order is being increased by 1.')
      N = N + 1;
      odd = 1;
   end
end
if nw ~= 0 & nw ~= N
   error('The window length must be the same as the filter length.')
end
if nw > 0
   wind = Wind;
else        	 % replace the following with the default window of your choice.
   wind = hamming(N);
end
%
% to use Chebyshev window, you must specify ripple
% ripple=60;
% wind=chebwin(N, ripple);
%
% to use Kaiser window, beta must be supplied
% att = 60; % dB of attenuation desired in sidelobe
% beta = .1102*(att-8.7);
% wind = kaiser(N,beta);

fl = Wn(1)/2;
if (Btype == 2 | Btype == 4)
   fh = Wn(2)/2;
   if (fl >= .5 | fl <= 0 | fh >= .5 | fh <= 0.)
      error('Frequencies must fall in range between 0 and 1.')
   end
   c1=fh-fl;
   if (c1 <= 0)
      error('Wn(1) must be less than Wn(2).')
   end
else
   c1=fl;
   if (fl >= .5 | fl <= 0)
      error('Frequency must lie between 0 and 1')
   end
end

nhlf = fix((N + 1)/2);
i1=1 + odd;

if Btype == 1			% lowpass
if odd
   b(1) = 2*c1;
end
xn=(odd:nhlf-1) + .5*(1-odd);
c=pi*xn;
c3=2*c1*c;
b(i1:nhlf)=(sin(c3)./c);
b = real([b(nhlf:-1:i1) b(1:nhlf)].*wind(:)');
gain = abs(polyval(b,1));
b = b/gain;
return;

elseif Btype ==2		% bandpass
b(1) = 2*c1;
xn=(odd:nhlf-1)+.5*(1-odd);
c=pi*xn;
c3=c*c1;
b(i1:nhlf)=2*sin(c3).*cos(c*(fl+fh))./c;
b=real([b(nhlf:-1:i1) b(1:nhlf)].*wind(:)');
gain = abs(polyval(b,exp(sqrt(-1)*pi*(fl+fh))));
b = b/gain;
return;

elseif Btype == 3		% highpass
b(1)=2*c1;
xn=(odd:nhlf-1);
c=pi*xn;
c3=2*c1*c;
b(i1:nhlf)=sin(c3)./c;
att=60.;
b=real([b(nhlf:-1:i1) b(1:nhlf)].*wind(:)');
b(nhlf)=1-b(nhlf);
b(1:nhlf-1)=-b(1:nhlf-1);
b(nhlf+1:N)=-b(nhlf+1:N);
gain = abs(polyval(b,-1));
b = b/gain;
return;

elseif Btype == 4		% bandstop
b(1) = 2*c1;
xn=(odd:nhlf-1)+.5*(1-odd);
c=pi*xn;
c3=c*c1;
b(i1:nhlf)=2*sin(c3).*cos(c*(fl+fh))./c;
b=real([b(nhlf:-1:i1) b(1:nhlf)].*wind(:)');
b(nhlf)=1-b(nhlf);
b(1:nhlf-1)=-b(1:nhlf-1);
b(nhlf+1:N)=-b(nhlf+1:N);
gain = abs(polyval(b,1));
b = b/gain;
return;
end

return