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);