RbfInterpolant.m

classdef RbfInterpolant < handle
    % implements rbf interpolation by thin plate splines or gaussian
    
    properties (SetAccess=private,Hidden=true)
        data_points = [];
        data_values = [];
        weights;
    end
    
    properties (Access=private)
        dilation_ = 0.1;
    end
    
    properties (SetAccess=private)
        num_points = 0;
        dimension;
        kernel_type = 'tps';
    end
    
    properties (Dependent)
        dilation;
    end
    
    methods
        
        % constructor: points is a (d x n)-matrix, values is a (n x 1)-matrix
        function this = RbfInterpolant(points, values, kernel_type, dilation)
            this.dimension = size(points, 1);
            if nargin >= 3
                if any(strcmp(kernel_type, {'tps', 'gaussian', 'gaussian_loc'}))
                    this.kernel_type = kernel_type;
                else
                    warning('RbfInterpolant:RbfInterpolant', 'kernel type not supported');
                end
            end
            if nargin == 4
                this.dilation_ = dilation;
            end
            if ~isempty(points)
                add_data(this, points, values, []);
            end
        end
        
        function set.dilation(this, dilation)
            this.dilation_ = dilation;
            update(this);
        end
        
        function dilation = get.dilation(this)
            dilation = this.dilation_;
        end
        
        % points is a (d x n)-matrix, res is a (1 x n)-matrix
        function res = evaluate(this, points)
            if this.num_points > 0
                res = this.weights * this.evaluate_basis(points);
            else
                res = zeros(0, size(points, 2));
            end
        end
        
        % points is a (d x n)-matrix, values is a (n x 1)-matrix
        function this = add_data(this, varargin)
            this.data_points = [this.data_points, varargin{1}];
            this.data_values = [this.data_values; varargin{2}];
            this.num_points = size(this.data_points, 2);
            if strcmp(this.kernel_type, 'gaussian_loc')
                this.dilation_ = [this.dilation_; varargin{3}(:)];
            end
            update(this);
        end
        
        % if indices is not specified, the last added data is removed
        function this = remove_data(this, indices)
            if nargin < 2
                indices = this.num_points;
            elseif any(indices > this.num_points)
                error('rbmatlab:RbfInterpolant:remove_data', ...
                    'Index exceeds number of data points.');
            end
            remain = setdiff(1:this.num_points, indices);
            this.data_points = this.data_points(:, remain);
            this.data_values = this.data_values(remain);
            this.num_points = length(remain);
            if strcmp(this.kernel_type, 'gaussian_loc')
                this.dilation = this.dilation_(remain);
            end
            update(this);
        end
        
        % points is a (d x n)-matrix, res is a (d x n)-matrix
        % can be vectorized
        function res = evaluate_gradient(this, points)
            r = this.compute_radii(points);
            res = zeros(this.dimension, size(points, 2));
            for i = 1:this.dimension
                inner = repmat(points(i, :), this.num_points, 1) ...
                    - repmat(this.data_points(i, :), size(points, 2), 1)';
                outer = (r>eps) .* (2*log(r + (r<=eps)) + 1);
                res(i, :) = this.weights * (outer .* inner);
            end
        end
        
        % points is a (d x n)-matrix, res is a (1 x n)-matrix
        % returns -inf if r=0
        % can be vectorized
        function res = evaluate_laplacian(this, points)
            r = this.compute_radii(points);
            res = zeros(1, size(points, 2));
            for i = 1:this.dimension
                inner = repmat(points(i, :), this.num_points, 1) ...
                    - repmat(this.data_points(i, :), size(points, 2), 1)';
                res = res + this.weights * (2*inner.^2./r + 2*log(r) + 1);
            end
        end
            
    end
    
    methods (Hidden=true)
        
        % vectorized thin plate spline basis
        function res = evaluate_basis(this, x)
            r = this.compute_radii(x);
            switch this.kernel_type
                case 'tps'
                    res = (r>eps) .* r.^2 .* log(r + (r<=eps));
                case 'gaussian'
                    res = exp( -(r/this.dilation_).^2 );
                case 'gaussian_loc'
                    res = exp( -(r./repmat(this.dilation_, 1, size(r, 2))).^2 );
            end
        end
        
        % overload subsref for evaluation
        function res = subsref(this, S)
            if strcmp(S(1).type, '()')
                res = evaluate(this, S(1).subs{1});
            else
                res = builtin('subsref', this, S);
            end
        end
        
        % compute weights from data
        function this = update(this)
            if this.num_points > 0
                this.weights = ((1e-9*eye(this.num_points) ...
                                 + this.evaluate_basis(this.data_points)') ...
                                \ this.data_values)';
            end
        end
        
        % compute radii for a point-matrix
        function r = compute_radii(this, x)
            x_shift = permute(repmat(x, 1, 1, this.num_points), [1 3 2]) ...
                - repmat(this.data_points, 1, 1, size(x, 2));
            r = sqrt(sum(x_shift.^2, 1));
            r = shiftdim(r);
        end
            
    end
    
end