MLSInterpolant

0001 classdef MLSInterpolant < Interpolant
0002     %MLSINTERPOLANT Local interpolant using Moving Least Squares
0003     
0004     properties        
0005         polyDeg; % The polynomial degree to be fitted
0006         rbfFun; % RBF functor
0007         searchObj; % The search object
0008         minSamples; % The minimum number of samples required in the neighborhood of a point so that the fit is possible
0009     end
0010     
0011     methods
0012         function obj = MLSInterpolant(x, y, z, varargin)
0013             %MLSINTERPOLANT Construct an instance of this class
0014             %   Detailed explanation goes here
0015             
0016             % Superclass constructor
0017             obj@Interpolant(x, y, z);
0018             
0019             % Validation functions for the input parameters
0020             validDistanceType = @(x) ischar(x) && strcmpi(x, 'euclidean') || strcmpi(x, 'haversine');
0021             validPolynomialDegree = @(x) isscalar(x) && x >= 0 && x < 3; % Interpolate up to 3rd degree polynomial
0022             validRBFType = @(x) isa(x, 'function_handle') || ...
0023                                 (ischar(x) && strcmpi(x, 'wendland'));
0024             validMinNumSamples = @(x) isscalar(x) && x > 0;
0025             
0026             % Check the input parameters using inputParser
0027             p = inputParser;
0028             addParameter(p, 'DistanceType', 'euclidean', validDistanceType);
0029             addParameter(p, 'PolynomialDegree', 0, validPolynomialDegree);            
0030             addParameter(p, 'RBF', 'wendland', validRBFType); % Note: this should be a compactly-supported RBF
0031             addParameter(p, 'RBFEpsilon', 1, @isscalar); % In this case, this parameter represents the support of the RBF and the range search query!
0032             addParameter(p, 'MinSamples', 10, validMinNumSamples);
0033             parse(p, varargin{:});
0034             
0035             % Get the RBF functor
0036             obj.rbfFun = rbfTypeToFunctor(p.Results.RBF, p.Results.RBFEpsilon);
0037             
0038             % Get the polynomial degree to fit
0039             obj.polyDeg = p.Results.PolynomialDegree;
0040             
0041             % Get the minimum number of points
0042             obj.minSamples = p.Results.MinSamples;
0043             
0044             % Check that the minimum number of samples desired is larger
0045             % than the minimum number of samples that correspond to the
0046             % selected polynomial degree
0047             coeffs = bivariatePolynomialTerms(obj.polyDeg, 1, 1);
0048             if obj.minSamples < numel(coeffs)
0049                 warning('The minimum number of samples set (%d) is less than the minimum number of samples that correspond to the selected polynomial degree. Settin the minimum number of samples to %d!', obj.minSamples, numel(coeffs));
0050                 obj.minSamples = numel(coeffs);
0051             end
0052             
0053             % Create the search object
0054             obj.searchObj = QueryNeighborhood(obj.data, 'Radius', p.Results.RBFEpsilon, ...
0055                                                         'SearchType', 'radial', ...
0056                                                         'DistanceType', p.Results.DistanceType);
0057         end
0058         
0059         function z = interpolate(obj, x, y)
0060             %INTERPOLATE MLS interpolates
0061             
0062             % Check input sizes
0063             Interpolant.checkSizes(x, y);
0064             
0065             % Get the neighbor points within the support for each input query point
0066             [ind, distances] = obj.searchObj.getNeighbors(x, y);
0067             
0068             z = zeros(size(x));
0069             for i = 1:numel(ind)
0070                 % The set of neighboring points
0071                 neighData = obj.data(ind{i}, :);
0072                 
0073                 % Check if there is a minimum number of points in the
0074                 % neighborhood
0075                 if size(neighData, 1) < obj.minSamples
0076                     z(i) = NaN; % Undefined value
0077                     continue;
0078                 end
0079                 
0080                 % Solve the local MLS fitting
0081                 polyCoeffs = obj.solveSingleMLS(neighData(:, 1), neighData(:, 2), neighData(:, 3), distances{i});
0082                 
0083                 % Evaluate the local polynomial at the query point
0084                 polyTerms = bivariatePolynomialTerms(obj.polyDeg, x(i), y(i));
0085                 z(i) = polyTerms*polyCoeffs;
0086             end        
0087         end
0088         
0089         function polyCoeffs = solveSingleMLS(obj, x, y, z, dist)
0090             
0091             numSamples = size(x, 1);
0092             
0093             % Construct matrix containing the terms
0094             P = bivariatePolynomialTerms(obj.polyDeg, x(:), y(:));
0095             
0096             % and the matrix constructed with weights
0097             W = zeros(numSamples);
0098             W(logical(eye(numSamples))) = obj.rbfFun(dist);
0099             
0100             % Solve the local fit
0101             A = P'*W*P;
0102             B = P'*W;
0103             polyCoeffs = inv(A)*B*z;
0104             
0105         end
0106     end
0107 end
0108