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