www.gusucode.com > mbcdata 工具箱 matlab 源码程序 > mbcdata/cgThinplateSpline.m
classdef cgThinplateSpline %cgThinplateSpline simple thinplate spline for use with application %point sets in optimizations % % obj = cgThinplateSpline.fit(X); % This class is intended to provide interpolants for a series of % responses and so works with matrices rather than vectors as the % intended application is to use the matrices as linear operators. % % cgThinplateSpline methods: % fit - fit thinplate spline for inputs (Static) % evaluate - evaluate thinplate spline matrix % fitData - fit thinplate spline to data (Static) % fitMode - fit thinplate spline to data for each mode (Static) % Copyright 2009 The MathWorks, Inc. and Ford Global Technologies, Inc. properties (Access=protected) Ainv Centers Scale = 1; Offset = 0; end methods function Y = evaluate(obj,X,Y) %evaluate thinplate spline if isempty(obj.Ainv) error(message('mbc:cgThinplateSpline:InvalidState')) end Atps = TPSmatrix(obj,X)*obj.Ainv; if nargin>2 Y = Atps*Y; else Y = Atps; end end end methods(Access=protected) function obj = cgThinplateSpline end function A = TPSmatrix(obj,X) %TPSmatrix thinplate spline matrix % RBF part X = scaleData(obj,X); Phi = xregrbfeval('thinplate', X, obj.Centers, [], []); % add polynomial part A = [Phi, ones(size(X,1), 1), X ]; end function X= scaleData(obj,X) %scaleData scale data in preparation for fitting and evaluation X = bsxfun(@minus,X,obj.Offset); X = bsxfun(@rdivide,X,obj.Scale); end function obj = calcScale(obj,X) obj.Offset=min(X,[],1); obj.Scale = max(X,[],1)-obj.Offset; obj.Scale(obj.Scale==0)=1; end end methods (Static) function [obj,OK] = fit(X) %FIT fit a thinplate spline to data % [obj,OK] = fit(X) obj = cgThinplateSpline; obj = calcScale(obj,X); obj.Centers = scaleData(obj,X); Atps = TPSmatrix(obj,X); n = size(X,2); A = [Atps; Atps(:,end-n:end)' zeros(n+1)]; rc = rcond(A); OK = rc > length(A)*eps(max(abs(A(:)))); if OK obj.Ainv = inv(A); %discard last n columns obj.Ainv = obj.Ainv(:,1:end-n-1); else obj.Ainv = []; end end function [OK,InterpMatrix]=interpData(X,Xinterp) %INTERPDATA fit thinplate spline to data % [OK,InterpMatrix]=interpData(X,Xinterp) % X input data to base interpolation on % Xinterp input data to interpolate % InterpMatrix interpolation matrix such that % Yinterp = InterpMatrix*Y [TPS,OK] = cgThinplateSpline.fit(X); InterpMatrix = []; if OK && nargout>1 InterpMatrix = evaluate(TPS, Xinterp ); end end function [OK,InterpMatrix]=interpByMode(X,Mode,Xinterp,XinterpMode) %INTERPBYMODE fit thinplate spline to data for each mode % [OK,InterpMatrix]=interpByMode(X,Mode,Y,YMode) % X input data to base interpolation on % Mode mode or group data for inputs % Xinterp input data to interpolate % XinterpMode mode data for interpolation data % InterpMatrix interpolation matrix such that % Yinterp = InterpMatrix*Y OK = true; InterpMatrix = zeros(length(Xinterp),size(X,1)); Done = false(length(Xinterp),1); for i = 1:max(Mode) SelX = i==Mode; SelInterp = XinterpMode==i; Done(SelInterp)=true; if OK && any(SelX) && any(SelInterp) % form interpolation for mode i if required Xi = X(SelX,:); [TPS,OK] = cgThinplateSpline.fit(Xi); if OK % form interpolation matrix InterpMatrix(SelInterp,SelX) = evaluate(TPS, Xinterp(SelInterp,:) ); end end end % All outputs must be estimated OK = OK && all( Done ); end end end