RBFInterpolant

0001 classdef RBFInterpolant < Interpolant
0002     %RBFINTERPOLANT Interpolant using Radial Basis Functions
0003     
0004     properties
0005         weights; % The weights of the RBF equation
0006         poly; % The constant polynomial part
0007         distFun; % Distance functor
0008         rbfFun; % RBF functor
0009     end
0010     
0011     methods
0012         function obj = RBFInterpolant(x, y, z, varargin)
0013             %RBFINTERPOLANT Construct an instance of this class
0014             %   Detailed explanation goes here
0015             
0016             % Validation functions for the input parameters
0017             validDistanceType = @(x) ischar(x) && strcmpi(x, 'euclidean') || strcmpi(x, 'haversine');
0018             validPolynomialDegree = @(x) isscalar(x) && x <= 3;
0019             validRBFType = @(x) isa(x, 'function_handle') || ...
0020                                 (ischar(x) && strcmpi(x, 'linear') || ...                                              
0021                                               strcmpi(x, 'cubic') || ...
0022                                               strcmpi(x, 'quintic') || ...
0023                                               strcmpi(x, 'multiquadric') || ...
0024                                               strcmpi(x, 'inversemultiquadric') || ...
0025                                               strcmpi(x, 'thinplate') || ...
0026                                               strcmpi(x, 'green') || ...
0027                                               strcmpi(x, 'greenwithtension') || ...
0028                                               strcmpi(x, 'greenregularizedwithtension') || ...
0029                                               strcmpi(x, 'tensionspline') || ...
0030                                               strcmpi(x, 'regularizedspline') || ...
0031                                               strcmpi(x, 'gaussian') || ...
0032                                               strcmpi(x, 'wendland'));
0033             validScalarLogical = @(x) isscalar(x) && islogical(x);
0034             
0035             % Check the input parameters using inputParser
0036             p = inputParser;
0037             addParameter(p, 'DistanceType', 'euclidean', validDistanceType);
0038             addParameter(p, 'PolynomialDegree', 0, validPolynomialDegree);
0039             addParameter(p, 'RBF', 'multiquadric', validRBFType);
0040             addParameter(p, 'RBFEpsilon', 0, @isscalar); % The additional parameter of some RBF
0041             addParameter(p, 'RBFEpsilonIsNormalized', false, validScalarLogical); % If set to true, the RBF
0042             addParameter(p, 'Regularization', 0, @isscalar); % Regularization coefficient to avoid matrix being close to singular (specially needed when using gaussianRBF). Can also be seen as a 'smoothing' parameter (set it to something >0 in order to APPROXIMATE instead of INTERPOLATE)
0043             parse(p, varargin{:});
0044             
0045             % If required to, scale the epsilon
0046             epsilon = p.Results.RBFEpsilon;
0047             if p.Results.RBFEpsilonIsNormalized
0048                 if epsilon < 0 || epsilon > 1
0049                     error('RBFEpsilonIsNormalized parameter is set to true, but RBFEpsilon is not between 0 and 1!');
0050                 end
0051                 epsilon = scaleEpsilonAccordingToRbfType(p.Results.RBF, epsilon, x, y, z);
0052             end
0053             
0054             % Superclass constructor
0055             obj@Interpolant(x, y, z);
0056             
0057             % Get the distanceFunctor
0058             distanceType = p.Results.DistanceType;
0059             if strcmpi(distanceType, 'euclidean')
0060                 obj.distFun = @(x, y) vecnorm(x-y, 2, 2);
0061             elseif strcmpi(distanceType, 'haversine')
0062                 obj.distFun = @haversine;
0063             else 
0064                 error('Unknown DistanceType');
0065             end
0066            
0067             % Get the RBF functor
0068             obj.rbfFun = rbfTypeToFunctor(p.Results.RBF, epsilon);                        
0069                         
0070             % Functor
0071             [obj.weights, obj.poly] = RBFInterpolant.solveWeightsAndPoly(x, y, z, obj.distFun, obj.rbfFun, p.Results.PolynomialDegree, p.Results.Regularization);
0072         end
0073         
0074         function z = interpolate(obj, x, y)
0075             %INTERPOLATE Summary of this method goes here
0076             %   Detailed explanation goes here
0077             
0078             % Check input sizes
0079             Interpolant.checkSizes(x, y);
0080             
0081             % Evaluate the RBF at the input points
0082             % Slower version, left here for reference (more readable)
0083 %             numPts = size(obj.data, 1);
0084 %             numQueryPts = numel(x);
0085 %             numPolyCoeffs = numel(obj.poly);
0086 %             z = zeros(size(x));
0087 %             for q = 1:numQueryPts
0088 %                 queryPt = [x(q), y(q)];
0089 %                 % Polynomial part
0090 %                 if numPolyCoeffs > 0
0091 %                     if numPolyCoeffs == 1
0092 %                         z(q) = obj.poly(1);
0093 %                     elseif numPolyCoeffs == 3
0094 %                         z(q) = queryPt*obj.poly(1) + queryPt*obj.poly(2) + obj.poly(3);
0095 %                     end
0096 %                 end
0097 %                 for i = 1:numPts
0098 %                     z(q) = z(q) + obj.weights(i)*obj.rbfFun(obj.distFun(obj.data(i, 1:2), queryPt));
0099 %                 end
0100 %             end
0101 
0102             % Faster version
0103             % Compute all pair-wise distances
0104             dists = pdist2([x(:) y(:)], obj.data(:, 1:2), obj.distFun);
0105             
0106             % Evaluate the RBF for all the distances
0107             A = obj.rbfFun(dists);
0108 %
0109 %             numPolyCoeffs = numel(obj.poly);
0110 %             if numPolyCoeffs > 0
0111 %                 if numPolyCoeffs == 1
0112 %                    A(:, end+1) = 1;
0113 %                    A(end+1, :) = 1;
0114 %                    A(end, end) = 0;
0115 %                 else % numPolyUnknowns == 3
0116 %                    A(:, end+1:end+3) = [x(:), y(:), ones(numel(x), 1)];
0117 %                    A(end+1:end+3, :) = [[x(:), y(:), ones(numel(x), 1)]' zeros(3,3)];
0118 %                 end
0119 %             end
0120 
0121             polyEval = bivariatePolynomialEval(obj.poly, x(:), y(:));
0122             
0123             z = A*obj.weights + polyEval;
0124             z = reshape(z, size(x));
0125         end
0126     end
0127     
0128     methods (Static)
0129         function [weights, poly] = solveWeightsAndPoly(x, y, z, distFun, rbfFun, polynomialDegree, regularizationCoeff)
0130             % Redundant check, since this would have returned a warning in
0131             % the constructor of the superclass, which is called before
0132             % this function. Still, since we made the function static, we
0133             % leave here this check just in case...
0134             if numel(x) ~= numel(y) || numel(x) ~= numel(z)
0135                 error('The number of elements in x, y and z must be the same');
0136             end
0137             numSamples = size(x, 1);
0138             
0139             % Compose the system of equations (slower, but more readable version, we leave it here commented for reference)
0140 %             A = zeros(numSamples, numSamples);
0141 %             for i = 1:numSamples
0142 %                 for j = 1:numSamples
0143 %                     A(i, j) = rbfFun(distFun([x(i), y(i)], [x(j), y(j)]));
0144 %                 end
0145 %             end
0146             
0147             % Compute all pair-wise distances
0148             dists = pdist([x y], distFun);
0149             
0150             % Evaluate the RBF for all the distances
0151             rbfEvals = rbfFun(dists);
0152             
0153             % Compose the system of equations
0154             % - RBF part
0155             A = zeros(numSamples, numSamples);
0156             A(tril(true(numSamples, numSamples), -1)) = rbfEvals; % Fill the lower triangular part of the matrix
0157             A = A + A'; % Mirror over the diagonal (matrix A is symmetric)
0158             % Compute RBF values at the diagonal (rbf(0))
0159             A(logical(eye(numSamples))) = rbfFun(0);
0160 
0161             % Regularizer?
0162             if regularizationCoeff ~= 0
0163                 A = A + eye(size(A, 1))*regularizationCoeff;
0164             end
0165             
0166             b = z(:);
0167             
0168             % - Polynomial part
0169             terms = bivariatePolynomialTerms(polynomialDegree, x(:), y(:));
0170             numTerms = size(terms, 2);
0171             A(:, end+1:end+numTerms) = terms;
0172             A(end+1:end+numTerms, :) = [terms' zeros(numTerms, numTerms)];
0173             b(end+1:end+numTerms) = 0;
0174 
0175             % Solve it
0176             solution = A\b;
0177             
0178             % Recover the results
0179             weights = solution(1:numSamples);
0180             if polynomialDegree >= 0
0181                 poly = solution(numSamples+1:end);
0182             else
0183                 poly = [];
0184             end
0185         end
0186     end
0187 end
0188