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

    function varargout = pRBFSurfaceGen(axisbp, pts, x, action, SurfaceGenData)
%PRBFSURFACEGEN RBF Surface Generation function
%
%   PRBFSURFACEGEN uses RBF extrapolation to generate a surface. This
%   function is called by CGTABGRADCONSRAINT using the syntaxes below.
%
%   SURFACEGENDATA = PRBFSURFACEGEN(AXISBP, PTS, X, 'FIXEDDATA') calculates
%   the interpolation matrices and the duplicate transform. These matrices
%   only depend on the fixed variables and will be stored in the working
%   store of the table gradient constraint. 
%
%   [Y, YGRAD] = PRBFSURFACEGEN(AXISBP, PTS, X, 'EVALUATE', SURFACEGENDATA)
%   returns an extrapolation of X from PTS to the axis breakpoints in
%   AXISBP. The gradient of this function at X is returned in YGRAD. Note
%   that the axis variables are scaled onto [0, 1].
%
%   [Y, YGRAD] = PRBFSURFACEGEN(AXISBP, PTS, X, 'EVALUATE', SURFACEGENDATA)
%   returns a extrapolation of multiple values of X in the case where X is
%   a NPTS-by-NRUN matrix. In this case, the axis values, PTS, are the same
%   for each column of X. Y is a NTAB-by-NRUN matrix of extrapolations of X
%   and YGRAD is a NTAB-by-NPTS-by-NRUN matrix of gradient evaluations.

%   This function uses an RBF interpolant. A thin plate spline interpolant
%   is used for 2-d tables and a cubic spline interpolant for 1-d tables.
%
%   See also PGETDG

%   Copyright 2006-2007 The MathWorks, Inc.

switch lower(action)
    case 'fixeddata'
        varargout{1} = i_FixedData(axisbp, pts);
    case 'evaluate'
        [varargout{1:nargout}] = i_Evaluate(axisbp, pts, x, SurfaceGenData);
end

%--------------------------------------------------------------------------
function SurfaceGenData = i_FixedData(axisbp, pts)
%--------------------------------------------------------------------------

SurfaceGenData = cell(1, 3);

persistent fRemoveDups1 fRemoveDups2 

if isempty(fRemoveDups1)
    % store function handles for cgmathobjects function in persistent data
    % to avoid unnecessary construction of these handles
    fRemoveDups1 = gethandle(cgmathsobject,'removeDuplicatePts1');
    fRemoveDups2 = gethandle(cgmathsobject,'removeDuplicatePts2');
end

% Determine RBF Kernel
NAXES = length(axisbp);
if NAXES == 2
    rbfkernel = 'thinplate';
elseif NAXES == 1
    rbfkernel = 'cubicrbf';
else
    % Do not calculate regression matrix if object has not been configured
    return
end

% Find the unique axis points. Repeats are replaced by the average of the
% output at that point. 
dummyX = ones(size(pts, 1));
if size(pts, 2) > 1
    % Tolerance set at 1e-8, determined by an external
    % study to determine the tolerance for this RBF Kernel.
    [newX, newY, unused, DupTransMat] = fRemoveDups2(pts(:, 1), ...
        pts(:, 2), dummyX, 1e-8);
    pts = [newX, newY];
else
    % Tolerance set at 1e-7, after a trial and error study in CAGE to
    % determine the tolerance for the cubic RBF Kernel.
    [pts, unused, DupTransMat] = fRemoveDups1(pts, dummyX, 1e-7);
end    

% Create the Vandermonde matrix
P = i_calcVandermonde(pts);

% Do the QR decomposition for the linear algebra trick for thin plate
% splines. Note need to return second argument, as there is a special
% calling case of QR with one return argument.
[Q, unused] = qr(P);
Q2 = Q(:, NAXES+2:end);

% Calculate Xp
Phi = xregrbfeval(rbfkernel, [], pts, [], []);
Xp = [Phi*Q2 P];

% Calculate Xt
% Calculate grid points
if NAXES == 2
    % Grid points
    [xx, yy] = ndgrid(axisbp{:});
    tpts = [xx(:), yy(:)];
elseif NAXES == 1
    % Grid points
    tpts = axisbp{1}(:);   
else
    tpts = [];
end
Phi = xregrbfeval(rbfkernel, tpts, pts, [], []);
P = i_calcVandermonde(tpts);
Xt = [Phi*Q2 P];

% Generate data to be cached in the working store
sizeXp = size(Xp);
if sizeXp(1) == sizeXp(2)
    SurfaceGenData{1} = Xt/Xp;
else
    [Qp, Rp] = qr(Xp, 0);
    SurfaceGenData{1} = Xt*(Rp\Qp');
end
SurfaceGenData{2} = Xt*(Xp\DupTransMat);
SurfaceGenData{3} = DupTransMat;

%--------------------------------------------------------------------------
function varargout = i_Evaluate(axisbp, pts, x, SurfaceGenData)
%--------------------------------------------------------------------------

% Map x onto unique x
x = SurfaceGenData{3}*x;

% Return transformation points
varargout{1} = SurfaceGenData{1}*x;

% Return gradient
if nargout > 1
    varargout{2} = repmat(SurfaceGenData{2}, [1 1 size(x, 2)]);
end

%--------------------------------------------------------------------------
function P = i_calcVandermonde(pts)
%--------------------------------------------------------------------------

NPTS = size(pts, 1);
P = [ones(NPTS, 1), pts];