QuadTreePURBFInterpolant

0001 classdef QuadTreePURBFInterpolant < Interpolant
0002     %QUADTREE QuadTreePURBFInterpolator class
0003     % This simplified implementation assumes all the points are known at construction
0004     % time, and no new insertions are expected after the initial construction
0005     
0006     properties
0007         root; % Root node of the tree (of type QTNode)
0008         weightingRBF; % The RBF used to compute the Sheppard's weights of the Partition of Unity
0009         distFun; % Distance function between points in the XY plane
0010     end
0011     
0012     methods
0013         function obj = QuadTreePURBFInterpolant(x, y, z, varargin)
0014             %QUADTREE Constructor of the class
0015             
0016             % Superclass constructor
0017             obj@Interpolant(x, y, z);
0018             
0019             % Check the input parameters using inputParser
0020             paramsToDelete = {'PolynomialDegree', 'RBF', 'RBFEpsilon', 'Smooth', 'Regularization'};
0021             vararginA = deleteParamsFromVarargin(paramsToDelete, varargin);
0022                         
0023             p = inputParser;
0024             validOverlap = @(x) isscalar(x) && x >= 0;
0025             validDistanceType = @(x) ischar(x) && strcmpi(x, 'euclidean') || strcmpi(x, 'haversine');
0026             validMinCellSizePercent = @(x) isscalar(x) && x <= 1 && x >= 0;
0027             validDomain = @(x) numel(x) == 4;
0028             addParameter(p, 'MinPointsInCell', 25, @isscalar); % The minimum number of points in a QuadTree cell
0029             addParameter(p, 'MinCellSizePercent', 0.1, validMinCellSizePercent); % The minimum cell size, expressed as a percentage of the maximum side length of the bounding box of the query domain
0030             addParameter(p, 'Overlap', 0.25, validOverlap); % Overlap between circles
0031             addParameter(p, 'DistanceType', 'euclidean', validDistanceType);
0032             addParameter(p, 'Domain', [], validDomain); % 4-elements vector specifying the bounding box of the query domain [xmin, ymin, width, height], so that the QuadTree is forced to cover the defined area and no query point gets out of the domain
0033             parse(p, vararginA{:});
0034             
0035             minPts = p.Results.MinPointsInCell;
0036             overlap = p.Results.Overlap;
0037             domain = p.Results.Domain;
0038             minCellSizePercent = p.Results.MinCellSizePercent;
0039             
0040             distanceType = p.Results.DistanceType;
0041             if strcmpi(distanceType, 'euclidean')
0042                 obj.distFun = @(x, y) vecnorm(x-y, 2, 2);
0043             elseif strcmpi(distanceType, 'haversine')
0044                 obj.distFun = @haversine;
0045             else 
0046                 error('Unknown DistanceType');
0047             end
0048             
0049             % Remove this parameters from varargin, as that variable will
0050             % be used later on to create the local RBF interpolator objects
0051             paramsToDelete = {'MinPointsInCell', 'MinCellSizePercent', 'Overlap', 'Domain'};
0052             vararginB = deleteParamsFromVarargin(paramsToDelete, varargin);
0053             
0054             % Compute the extents of the root node
0055             if isempty(domain)
0056                 % Compute the extend of the domain to be that of the input
0057                 % points
0058                 minX = min(obj.data(:, 1)); maxX = max(obj.data(:, 1));
0059                 minY = min(obj.data(:, 2)); maxY = max(obj.data(:, 2));
0060                 w = maxX-minX;
0061                 h = maxY-minY;
0062             else
0063                 minX = domain(1);
0064                 minY = domain(2);
0065                 w = domain(3);
0066                 h = domain(4);
0067             end
0068             
0069             % We want the space to be divided in squares, so we force the
0070             % width and height of the root node to be square!
0071             wh = max([w, h]);
0072             
0073             obj.root = QTNode(minX, minY, wh, wh, obj.data, obj.distFun, overlap);
0074             
0075             % Subdivide the root node (and effectively create the tree top to bottom)
0076             obj.root = obj.root.subdivide(minPts, wh*minCellSizePercent);
0077             
0078             % Remove stored points in non-leaf nodes, from now on we will
0079             % just use leaves
0080             obj.root.freeMemory();
0081             
0082             % Compute a local RBF interpolator for each leaf
0083             obj.root = obj.root.computeRBFAtLeafs(vararginB{:});            
0084         end
0085         
0086         function plot(obj, newFigure)
0087             % Plot the quadtree
0088             if nargin < 2
0089                 newFigure = false;
0090             end
0091             if newFigure
0092                 figure();
0093             end
0094             hold on;
0095             nodes = getLeaves(obj.root);
0096             for i = 1:numel(nodes)
0097                 hold on;
0098                 rectangle('Position', [nodes(i).x, nodes(i).y, nodes(i).w, nodes(i).h]);
0099                 hold off;
0100                 cx = nodes(i).x + (nodes(i).w/2);
0101                 cy = nodes(i).y + (nodes(i).h/2);                
0102                 r = nodes(i).getEuclideanRadius();
0103                 
0104                 plotCircle(cx, cy, r);
0105 %                 hold on;
0106 %                 plot(nodes(i).pts(:, 1), nodes(i).pts(:, 2), '.');
0107 %                 hold off;
0108 %                 pause;
0109             end
0110             hold on;
0111             plot(obj.data(:, 1), obj.data(:, 2), '.k');
0112             hold off;      
0113             axis equal;
0114         end        
0115         
0116         function z = interpolate(obj, x, y)
0117             %interpolate Interpolates the z value at the x, y points.
0118             % Note: A NaN value will be generated for those points not
0119             % covered in the domain of the QuadTree.
0120             
0121             % Check input sizes
0122             Interpolant.checkSizes(x, y);
0123             
0124             % Reshape input, for convenience
0125             orSize = size(x);
0126             x = x(:);
0127             y = y(:);
0128             
0129 %             % Evaluate the RBF of each point (version traversing the tree for each point, very slow, left here for reference)
0130 %             numPts = numel(x);
0131 %             z = ones(numel(x), 1);
0132 %             for i = 1:numPts
0133 %                 % Evaluate the point in the RBF of each point
0134 %                 [f, w] = obj.root.evalRBF(x(i), y(i), obj.overlap);
0135 %
0136 %                 if isempty(f)
0137 %                     % The query point is out of range, set it to an invalid
0138 %                     % value
0139 %                     z(i) = NaN;
0140 %                     continue;
0141 %                 end
0142 %
0143 %                 % Compute the Weighting function for each
0144 %                 w = w ./ sum(w);
0145 %
0146 %                 % Weight the contribution of each RBF evaluation
0147 %                 z(i) = sum(w.*f);
0148 %             end
0149 
0150             % Evaluate the RBF of each point (version accessing the leaves at once, way faster)
0151 
0152             % Get the leaf nodes
0153             leaves = getLeaves(obj.root);
0154             
0155             % Compute the centers and radius of each leaf
0156             numLeaves = numel(leaves);
0157             centers = zeros(numLeaves, 2);
0158             radius = zeros(numLeaves, 1);
0159             for i = 1:numLeaves
0160                 centers(i, :) = leaves(i).getCenter();
0161                 radius(i) = leaves(i).getRadius();                
0162             end
0163             
0164             dists = pdist2([x, y], centers, obj.distFun); 
0165 %             dists = pdist2([x, y], centers);
0166             
0167             f = zeros(numel(x), 1); % We will store here the accumulation of each RBF contribution (weighted locally)
0168             w = zeros(numel(x), 1); % And here the accomulation of the weighting function
0169             for i = 1:numLeaves
0170                 % Find those input points falling in the current leaf
0171                 ind = dists(:, i) <= radius(i);
0172                 
0173                 % Compute the RBF value at those points
0174                 rbfEval = leaves(i).rbfInterp.interpolate(x(ind), y(ind));
0175                 
0176                 % Compute the weighting function at those points
0177                 weights = leaves(i).weightingRBF(dists(ind, i));
0178                 
0179                 % Apply the weights to the corresponding function and
0180                 % accumulate
0181                 f(ind) = f(ind) + rbfEval.*weights;
0182                 
0183                 % Accumulate the weights for the final division
0184                 w(ind) = w(ind)+weights;
0185             end
0186             z = f./w; % Here a division by w = 0 will occur for those points outside the domain covered by the quadtree. Since this will result in a NaN, we use this value as an indicator that the z is undefined at that point.
0187             
0188             % Get z back to the original shape of the input
0189             z = reshape(z, orSize);
0190         end
0191     end
0192 end
0193