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

    function V = extrapolate_RBF_points(X, Y, Z, Xi, Yi)
%EXTRAPOLATE_RBF_POINTS Extrapolate RBF at a set of points
%
%   V = EXTRAPOLATE_RBF_POINTS(X, Y, Z, BPx, BPy) extrapolates from a
%   defined set of values (X, Y, Z) to generate the values of Z at the
%   points [Xi, Yi].
%
%   For the different extrapolation cases, see the help for
%   CGMATHSOBJECT/EXTRAPOLATE_RBF.
% 
%   See also CGMATHSOBJECT/EXTRAPOLATE_RBF

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


if numel(X)~=numel(Z) || numel(Y)~=numel(Z)
    error(message('mbc:extrapolate_RBF:InvalidArgument'));
end

X = X(:);
Y = Y(:);
Z = Z(:);
Xi = Xi(:);
Yi = Yi(:);

% Preprocess data for extrapolation
[X, Y, Z] = pPrepareForExtrapolate(X, Y, Z);

% Perform special case extrapolations
[V, Done] = pExtrapolateSpecialCases(X, Y, Z, Xi, Yi);
if Done
    return
end

% Use an RBF if none of the special cases were picked up

% Scale all X and Y values
[Smn,Smx] = i_getScale([X; Xi]);
X = i_scale(X, Smn,Smx);
Xi = i_scale(Xi, Smn,Smx);

[Smn,Smx] = i_getScale([Y; Yi(:)]);
Y = i_scale(Y, Smn,Smx);
Yi = i_scale(Yi, Smn,Smx);

% Remove duplicate points
[X, Y, Z] = removeDuplicatePts2(X, Y, Z, 1e-8);

% RBF part
Centers = [X, Y];
Phi = xregrbfeval('thinplate', [], Centers, [], []);

% Polynomial part
P = [ones(length(X), 1), X, Y];

% Solve system
A = [Phi, P; P', zeros(3)];
[beta,rc] = linsolve(A,[Z; zeros(3,1)]);

if rc>eps(max(A(:)))*size(A,1)*2
    % Evaluate RBF system at the required points
    rbfPart = xregrbfeval('thinplate', [Xi, Yi], Centers, [], beta(1:length(X)));
    polyPart = [ones(numel(Xi),1), Xi, Yi] * beta((length(X)+1):end);
    V = rbfPart+polyPart;
else
    % return NaNs if the system is singular
    V = NaN(size(Xi,1),1);
end

% Helper functions that do scaling of the X/Y values
function Out = i_scale(In, mn,mx)
if nargin<2
    [mn,mx] = i_getScale(In);
end
Out = (In-mn)/(mx-mn);

function [mn,mx] = i_getScale(In)
mx = max(In);
mn = min(In);
if mx == mn
    mn = mn-1;
    mx = mx+1;
end