QTNode

0001 classdef QTNode
0002     %Quadtree Partition of Unity RBF Node
0003     
0004     properties
0005         x; % Center of the node in X
0006         y; % Center of the node in Y
0007         w; % Width of the node
0008         h; % Height of the node
0009         pts = []; % Points in this node (of size numPts x 3)
0010         childs = []; % Children of the node (up to 4)
0011         rbfInterp = []; % The interpolator function (just computed at leaves)
0012         weightingRBF = []; % The weighting CSRBF
0013         distFun; % Distance function between points in the XY plane
0014         overlap; % Amount of overlap between circles (scalar between 0 and 1). The base radius of a cell (overlap = 0) in the QuadTree corresponds to half the length of its diagonal. This parameter is a multiplicative factor applied to this base radius.
0015     end
0016     
0017     methods
0018         function obj = QTNode(x, y, w, h, pts, distFun, overlap)
0019             % Constructor of the class
0020             obj.x = x; 
0021             obj.y = y;
0022             obj.w = w;
0023             obj.h = h;            
0024             obj.childs = [];      
0025             obj.distFun = distFun;
0026             obj.overlap = overlap;
0027             
0028             % From the input points, just retain those within the limits
0029             ind = obj.ptsInNodeCircle(pts);
0030             obj.pts = pts(ind, :);            
0031         end
0032                 
0033         function center = getCenter(obj)
0034             %getCenter Gets the center point of the node
0035             center = [obj.x+(obj.w/2), obj.y+(obj.h/2)];
0036         end
0037         
0038         function radius = getRadius(obj)
0039             %getRadius Gets the radius of the node using the internal
0040             %distance function (required while evaluating the RBF)
0041             %   it also applies the "overlap" factor
0042             diag = getDiagonalLength(obj);
0043             radius = diag*0.5;
0044             radius = radius+radius*obj.overlap;
0045         end
0046         
0047         function diagLength = getDiagonalLength(obj)
0048             %getDiagonalLength Gets the length of the diagonal of the square represented by
0049             % the node
0050             diagLength = obj.distFun([obj.x+obj.w, obj.y+obj.h], [obj.x, obj.y]);
0051         end
0052         
0053         function diagLength = getEuclideanDiagonalLength(obj)
0054             %getDiagonalLength Gets the length of the diagonal of the square represented by
0055             % the node
0056             diagLength = norm([obj.x+obj.w, obj.y+obj.h]-[obj.x, obj.y]);
0057         end
0058         
0059         function radius = getEuclideanRadius(obj)
0060             %getRadius Gets the radius of the node using Euclidean distance
0061             %on the XY plane
0062             %   it also applies the "overlap" factor
0063             diag = getEuclideanDiagonalLength(obj);
0064             radius = diag*0.5;
0065             radius = radius+radius*obj.overlap;            
0066         end
0067         
0068         function b = isLeaf(obj)
0069             b = isempty(obj.childs);
0070         end
0071         
0072         function b = ptsInNodeCircle(obj, pts)
0073 %             center = obj.getCenter();
0074 %             radius = obj.getRadius(overlap);
0075 %             rads = obj.distFun([x, y], repmat(center, numel(x), 1), 2, 2);
0076 %             b = rads <= radius;
0077             
0078                         % Check if the points fall within the circle using Euclidean
0079             % distance.
0080             
0081             % Keep just those points falling within a circle
0082             center = obj.getCenter();
0083             radius = obj.getEuclideanRadius();
0084             
0085             rads = vecnorm(pts(:, 1:2)-center, 2, 2);
0086             b = rads <= radius;
0087         end
0088         
0089         function nodes = getLeaves(node)
0090             % If the children nodes are empty, this is a leaf
0091             if node.isLeaf()
0092                 nodes = node;
0093             else
0094                 nodes = [];
0095                 for i = 1:4
0096                     nodes = [nodes, getLeaves(node.childs(i))];
0097                 end
0098             end            
0099         end
0100         
0101         function obj = subdivide(obj, minPts, minCellSideLength)
0102             %subdivide Recursive subdivision of the Quadtree
0103             
0104             % End of recursion if we have less than minimum number of points in the cell
0105             if size(obj.pts, 1) < minPts
0106                 return;
0107             end
0108             
0109             % The width/height of children nodes is halved
0110             w = obj.w/2;
0111             h = obj.h/2;
0112             
0113             % End of recursion if children will have a length smaller than the minimum cell length
0114             if w < minCellSideLength || h < minCellSideLength
0115                 return;
0116             end
0117             
0118             % Create the 4 childrens and span subdivision
0119             % South-West node
0120             swNode = QTNode(obj.x, obj.y, w, h, obj.pts, obj.distFun, obj.overlap);             
0121             if size(swNode.pts, 1) < minPts
0122                 % We also end recursion if one of the childrens to span has
0123                 % less than minPts, because this would mean that we should
0124                 % not be dividing the current node
0125                 return;
0126             end
0127             % North-West node
0128             nwNode = QTNode(obj.x+w, obj.y, w, h, obj.pts, obj.distFun, obj.overlap);             
0129             if size(nwNode.pts, 1) < minPts
0130                 return;
0131             end
0132             % South-East node
0133             seNode = QTNode(obj.x, obj.y+h, w, h, obj.pts, obj.distFun, obj.overlap);             
0134             if size(seNode.pts, 1) < minPts
0135                 return;
0136             end
0137             % North-East node
0138             neNode = QTNode(obj.x+w, obj.y+h, w, h, obj.pts, obj.distFun, obj.overlap); 
0139             if size(neNode.pts, 1) < minPts
0140                 return;
0141             end
0142             
0143             % Span subdivision
0144             swNode = swNode.subdivide(minPts, minCellSideLength);
0145             nwNode = nwNode.subdivide(minPts, minCellSideLength);
0146             seNode = seNode.subdivide(minPts, minCellSideLength);
0147             neNode = neNode.subdivide(minPts, minCellSideLength);
0148             
0149             obj.childs = [swNode, seNode, nwNode, neNode];
0150 %             obj.pts = []; % Free space, points are stored only at leaves
0151         end
0152         
0153         function obj = computeRBFAtLeafs(obj, varargin)
0154             if obj.isLeaf()
0155                 % Compute the local RBF interpolant corresponding to this node
0156                 obj.rbfInterp = RBFInterpolant(obj.pts(:, 1), obj.pts(:, 2), obj.pts(:, 3), varargin{:});
0157                 % Create the weighting RBF
0158                 r = obj.getRadius();
0159                 obj.weightingRBF = rbfTypeToFunctor('wendland', r);
0160             else
0161                 % Recurse down the tree
0162                 for i = 1:4
0163                     obj.childs(i) = obj.childs(i).computeRBFAtLeafs(varargin{:});
0164                 end
0165             end
0166         end
0167         
0168         function [f, w] = evalRBF(obj, x, y)
0169             % Evaluates the RBF value of a point at each intersecting leaf
0170             % It also returns the weight at the evaluation point
0171             
0172             if obj.ptsInNodeCircle([x, y])            
0173                 if obj.isLeaf()
0174                     % Evaluate the RBF
0175                     f = obj.rbfInterp.interpolate(x, y);
0176                     w = obj.weightingRBF(obj.distFun([x, y], obj.getCenter()));
0177                 else
0178                     % Recurse down the tree
0179                     f = [];
0180                     w = [];
0181                     for i = 1:4
0182                         [fc, wc] = evalRBF(obj.childs(i), x, y);
0183                         f = [f; fc];
0184                         w = [w; wc];
0185                     end                
0186                 end
0187             else
0188                 f = [];
0189                 w = [];
0190             end
0191         end
0192         
0193         function obj = freeMemory(obj)
0194             % Remove stored points in non-leaf nodes
0195             %  Use this function after creating the QuadTree (i.e., after
0196             %  subdivide and correct methods
0197             
0198             if isempty(obj.childs)
0199                 % Leaf node, end recursion
0200                 return;
0201             else
0202                 % Non-leaf node, eliminate points
0203                 obj.pts = [];
0204                 % And follow traversal
0205                 obj.childs(1) = obj.childs(1).freeMemory();
0206                 obj.childs(2) = obj.childs(2).freeMemory();
0207                 obj.childs(3) = obj.childs(3).freeMemory();
0208                 obj.childs(4) = obj.childs(4).freeMemory();
0209             end
0210             
0211         end
0212     
0213     end
0214     
0215     methods (Static)
0216         function pts = ptsInNode(x, y, w, h, pts)
0217             % Keep just those points falling within the node (quad-shape)
0218             ptsInNodeInds = pts(:, 1) >= x & pts(:, 1) <= x+w & pts(:, 2) >= y & pts(:, 2) <= y+h;
0219             pts = pts(ptsInNodeInds, :);
0220         end
0221         
0222 %         function pts = ptsInNodeCircle(x, y, w, h, pts, overlap)
0223 %             % Check if the points fall within the circle using Euclidean
0224 %             % distance.
0225 %
0226 %             % Keep just those points falling within a circle
0227 %             center = [x+(w/2), y+(h/2)];
0228 %             diag = norm([x+w, y+h]-[x, y]);
0229 %             radius = diag*0.5;
0230 %             radius = radius+radius*overlap;
0231 %
0232 %             rads = vecnorm(pts(:, 1:2)-center, 2, 2);
0233 %             pts = pts(rads <= radius, :);
0234 %         end
0235     end
0236 end
0237