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