www.gusucode.com > mbcdata 工具箱 matlab 源码程序 > mbcdata/@cgmathsobject/private/extinterp2.m

    function F = extinterp2(X,Y,Z,U,V)
%EXTINTERP2 Extends interp2 to provide extrapolation outside the bounds
%
%   ZI = EXTINTERP2(X,Y,Z,XI,YI) interpolates to find ZI, the values of the
%   underlying 2-D function Z at the points in matrices XI and YI.
%   Matrices X and Y specify the points at which the data Z is given.
%   
%   XI can be a row vector, in which case it specifies a matrix with
%   constant columns. Similarly, YI can be a column vector and it specifies
%   a matrix with constant rows. 
%
%	The function extrapolates by extend the mesh originally specified by X
%	and Y. Thus if originally X and Y gave an n by m grid, then a new n by
%	m grid will be created where the maximal and minimal values are those
%	largest and smallest of the values in the combined X and Xi entities. Z
%	values are adjusted by extending the previous bilinear interpolations
%	into the new region.
%
%   This function will use bilinear interpolation and extrapolation to
%   obtain estimates of a function at grid points specified by u and v.
%   Unfortunately due to the way this function works, certain error
%   messages may not be as helpful as in interp2. 

%  Copyright 2000-2012 The MathWorks, Inc. and Ford Global Technologies, Inc.



% Preprocess data for extrapolation
[X, Y, Z] = pPrepareForExtrapolate2(X, Y, Z);
if any(size(Z) < 2)
    error(message('mbc:extinterp2:InvalidArgument'));
end

% Find min/max of evaluation locations 
ut=max(U(:));
ub=min(U(:));
vt=max(V(:));
vb=min(V(:));

[m,n]=size(X);
% Get vector of X gridpoints
if m == 1
    Xv = X;
elseif n == 1
    Xv = X';
else 
    Xv=X(1,:);
end
[k,l]=size(Y);
% Do same for Y.
if k == 1
    Yv = Y;
elseif l == 1
    Yv = Y';
else 
    Yv = Y(:,1)';
end

% Find min and max of x/y gridpoints augmented with xi/yi points
xb=min(X(1),ub);
xt=max(max(X(:)),ut);
yb=min(Y(1),vb);
yt=max(max(Y(:)),vt);

% Replace min/max of x/y breakpoints with min/max of augmented points
Xu=[xb Xv(2:length(Xv)-1) xt];
Yu=[yb Yv(2:length(Yv)-1) yt];
   
% Form mesh grid of new breakpoints
[XX,YY]=meshgrid(Xu,Yu);

% Replace elements of Z matrix at min and max X/Y breakpoints
[p,q]=size(Z);
ZZ=zeros(size(Z));
ZZ(2:p-1,2:q-1)=Z(2:p-1,2:q-1);
ZZ(1,2:q-1)=Z(1,2:q-1)+(Z(2,2:q-1)-Z(1,2:q-1))/(Yv(2)-Yv(1))*(Yu(1)-Yv(1));
ZZ(p,2:q-1)=Z(p,2:q-1)+(Z(p,2:q-1)-Z(p-1,2:q-1))/(Yv(p)-Yv(p-1))*(Yu(p)-Yv(p));
ZZ(2:p-1,1)=Z(2:p-1,1)+(Z(2:p-1,2)-Z(2:p-1,1))/(Xv(2)-Xv(1))*(Xu(1)-Xv(1));
ZZ(2:p-1,q)=Z(2:p-1,q)+(Z(2:p-1,q)-Z(2:p-1,q-1))/(Xv(q)-Xv(q-1))*(Xu(q)-Xv(q));
ZZ(1,1)=lin(Z(1,1),Z(1,2),Z(2,1),Z(2,2),Xv(1),Yv(1),Xv(2),Yv(2),Xu(1),Yu(1));
ZZ(1,q)=lin(Z(1,q-1),Z(1,q),Z(2,q-1),Z(2,q),Xv(q-1),Yv(1),Xv(q),Yv(2),Xu(q),Yu(1));
ZZ(p,1)=lin(Z(p-1,1),Z(p-1,2),Z(p,1),Z(p,2),Xv(1),Yv(p-1),Xv(2),Yv(p),Xu(1),Yu(p));
ZZ(p,q)=lin(Z(p-1,q-1),Z(p-1,q),Z(p,q-1),Z(p,q),Xv(q-1),Yv(p-1),Xv(q),Yv(p),Xu(q),Yu(p));

% Now can call interp2
% interp2 does not accept mixed orientation inputs
if ~isequal(size(U), size(V)) && isvector(U) && isvector(V)
    % if U and V are mixed orientation we need to call meshgrid
    [U, V] = meshgrid(U,V);
end    
F=interp2(XX,YY,ZZ,U,V);

function  K=lin(a,b,c,d,u,p,v,q,x,y)

% Compute values of ZZ matrix at corners.
K=-(v*q*a-x*q*a-y*v*a+x*y*a-b*u*q+b*x*q+b*y*u-x*y*b-c*v*p+c*x*p...
    +c*y*v-x*y*c+u*p*d-x*p*d-y*u*d+x*y*d)/(-v*q+v*p+u*q-u*p);