Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TKDTree.cxx
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Authors: Marian Ivanov and Alexandru Bercuci 04/03/2005
3
4/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12#include "TKDTree.h"
13#include "TRandom.h"
14
15#include "TString.h"
16#include <cstring>
17#include <limits>
18
20
21
22/**
23\class TKDTree
24\brief Class implementing a kd-tree
25
26Contents:
271. What is kd-tree
282. How to cosntruct kdtree - Pseudo code
293. Using TKDTree
30 a. Creating the kd-tree and setting the data
31 b. Navigating the kd-tree
324. TKDTree implementation - technical details
33 a. The order of nodes in internal arrays
34 b. Division algorithm
35 c. The order of nodes in boundary related arrays
36
37
38
39### 1. What is kdtree ? ( [http://en.wikipedia.org/wiki/Kd-tree] )
40
41In computer science, a kd-tree (short for k-dimensional tree) is a space-partitioning data structure
42for organizing points in a k-dimensional space. kd-trees are a useful data structure for several
43applications, such as searches involving a multidimensional search key (e.g. range searches and
44nearest neighbour searches). kd-trees are a special case of BSP trees.
45
46A kd-tree uses only splitting planes that are perpendicular to one of the coordinate system axes.
47This differs from BSP trees, in which arbitrary splitting planes can be used.
48In addition, in the typical definition every node of a kd-tree, from the root to the leaves, stores a point.
49This differs from BSP trees, in which leaves are typically the only nodes that contain points
50(or other geometric primitives). As a consequence, each splitting plane must go through one of
51the points in the kd-tree. kd-trees are a variant that store data only in leaf nodes.
52
53### 2. Constructing a classical kd-tree ( Pseudo code)
54
55Since there are many possible ways to choose axis-aligned splitting planes, there are many different ways
56to construct kd-trees. The canonical method of kd-tree construction has the following constraints:
57
58* As one moves down the tree, one cycles through the axes used to select the splitting planes.
59 (For example, the root would have an x-aligned plane, the root's children would both have y-aligned
60 planes, the root's grandchildren would all have z-aligned planes, and so on.)
61* At each step, the point selected to create the splitting plane is the median of the points being
62 put into the kd-tree, with respect to their coordinates in the axis being used. (Note the assumption
63 that we feed the entire set of points into the algorithm up-front.)
64
65This method leads to a balanced kd-tree, in which each leaf node is about the same distance from the root.
66However, balanced trees are not necessarily optimal for all applications.
67The following pseudo-code illustrates this canonical construction procedure (NOTE, that the procedure used
68by the TKDTree class is a bit different, the following pseudo-code is given as a simple illustration of the
69concept):
70
71~~~~
72function kdtree (list of points pointList, int depth)
73{
74 if pointList is empty
75 return nil;
76 else
77 {
78 // Select axis based on depth so that axis cycles through all valid values
79 var int axis := depth mod k;
80
81 // Sort point list and choose median as pivot element
82 select median from pointList;
83
84 // Create node and construct subtrees
85 var tree_node node;
86 node.location := median;
87 node.leftChild := kdtree(points in pointList before median, depth+1);
88 node.rightChild := kdtree(points in pointList after median, depth+1);
89 return node;
90 }
91}
92~~~~
93
94Our construction method is optimized to save memory, and differs a bit from the constraints above.
95In particular, the division axis is chosen as the one with the biggest spread, and the point to create the
96splitting plane is chosen so, that one of the two subtrees contains exactly 2^k terminal nodes and is a
97perfectly balanced binary tree, and, while at the same time, trying to keep the number of terminal nodes
98in the 2 subtrees as close as possible. The following section gives more details about our implementation.
99
100### 3. Using TKDTree
101
102#### 3a. Creating the tree and setting the data
103 The interface of the TKDTree, that allows to set input data, has been developed to simplify using it
104 together with TTree::Draw() functions. That's why the data has to be provided column-wise. For example:
105
106\code{.cpp}
107 {
108 TTree *datatree = ...
109 ...
110 datatree->Draw("x:y:z", "selection", "goff");
111 //now make a kd-tree on the drawn variables
112 TKDTreeID *kdtree = new TKDTreeID(npoints, 3, 1);
113 kdtree->SetData(0, datatree->GetV1());
114 kdtree->SetData(1, datatree->GetV2());
115 kdtree->SetData(2, datatree->GetV3());
116 kdtree->Build();
117 }
118\endcode
119
120 NOTE, that this implementation of kd-tree doesn't support adding new points after the tree has been built
121 Of course, it's not necessary to use TTree::Draw(). What is important, is to have data columnwise.
122 An example with regular arrays:
123
124\code{.cpp}
125 {
126 Int_t npoints = 100000;
127 Int_t ndim = 3;
128 Int_t bsize = 1;
129 Double_t xmin = -0.5;
130 Double_t xmax = 0.5;
131 Double_t *data0 = new Double_t[npoints];
132 Double_t *data1 = new Double_t[npoints];
133 Double_t *data2 = new Double_t[npoints];
134 Double_t *y = new Double_t[npoints];
135 for (Int_t i=0; i<npoints; i++){
136 data0[i]=gRandom->Uniform(xmin, xmax);
137 data1[i]=gRandom->Uniform(xmin, xmax);
138 data2[i]=gRandom->Uniform(xmin, xmax);
139 }
140 TKDTreeID *kdtree = new TKDTreeID(npoints, ndim, bsize);
141 kdtree->SetData(0, data0);
142 kdtree->SetData(1, data1);
143 kdtree->SetData(2, data2);
144 kdtree->Build();
145 }
146\endcode
147
148 By default, the kd-tree doesn't own the data and doesn't delete it with itself. If you want the
149 data to be deleted together with the kd-tree, call TKDTree::SetOwner(kTRUE).
150
151 Most functions of the kd-tree don't require the original data to be present after the tree
152 has been built. Check the functions documentation for more details.
153
154#### 3b. Navigating the kd-tree
155
156 Nodes of the tree are indexed top to bottom, left to right. The root node has index 0. Functions
157 TKDTree::GetLeft(Index inode), TKDTree::GetRight(Index inode) and TKDTree::GetParent(Index inode)
158 allow to find the children and the parent of a given node.
159
160 For a given node, one can find the indexes of the original points, contained in this node,
161 by calling the GetNodePointsIndexes(Index inode) function. Additionally, for terminal nodes,
162 there is a function GetPointsIndexes(Index inode) that returns a pointer to the relevant
163 part of the index array. To find the number of point in the node
164 (not only terminal), call TKDTree::GetNpointsNode(Index inode).
165
166### 4. TKDtree implementation details - internal information, not needed to use the kd-tree.
167
168#### 4a. Order of nodes in the node information arrays:
169
170TKDtree is optimized to minimize memory consumption.
171Nodes of the TKDTree do not store pointers to the left and right children or to the parent node,
172but instead there are several 1-d arrays of size fNNodes with information about the nodes.
173The order of the nodes information in the arrays is described below. It's important to understand
174it, if one's class needs to store some kind of additional information on the per node basis, for
175example, the fit function parameters.
176
177- Drawback: Insertion to the TKDtree is not supported.
178- Advantage: Random access is supported
179
180As noted above, the construction of the kd-tree involves choosing the axis and the point on
181that axis to divide the remaining points approximately in half. The exact algorithm for choosing
182the division point is described in the next section. The sequence of divisions is
183recorded in the following arrays:
184~~~~
185fAxis[fNNodes] - Division axis (0,1,2,3 ...)
186fValue[fNNodes] - Division value
187~~~~
188
189Given the index of a node in those arrays, it's easy to find the indices, corresponding to
190children nodes or the parent node:
191Suppose, the parent node is stored under the index inode. Then:
192- Left child `index = inode*2+1`
193- Right child `index = (inode+1)*2`
194
195Suppose, that the child node is stored under the index inode. Then:
196- Parent `index = inode/2`
197
198Number of division nodes and number of terminals :
199`fNNodes = (fNPoints/fBucketSize)`
200
201The nodes are filled always from left side to the right side:
202Let inode be the index of a node, and irow - the index of a row
203The TKDTree looks the following way:
204Ideal case:
205~~~~
206Number of _terminal_ nodes = 2^N, N=3
207
208 INode
209irow 0 0 - 1 inode
210irow 1 1 2 - 2 inodes
211irow 2 3 4 5 6 - 4 inodes
212irow 3 7 8 9 10 11 12 13 14 - 8 inodes
213~~~~
214
215Non ideal case:
216~~~~
217Number of _terminal_ nodes = 2^N+k, N=3 k=1
218
219 INode
220irow 0 0 - 1 inode
221irow 1 1 2 - 2 inodes
222irow 2 3 4 5 6 - 3 inodes
223irow 3 7 8 9 10 11 12 13 14 - 8 inodes
224irow 4 15 16 - 2 inodes
225~~~~
226
227#### 4b. The division algorithm:
228
229As described above, the kd-tree is built by repeatingly dividing the given set of points into
2302 smaller sets. The cut is made on the axis with the biggest spread, and the value on the axis,
231on which the cut is performed, is chosen based on the following formula:
232Suppose, we want to divide n nodes into 2 groups, left and right. Then the left and right
233will have the following number of nodes:
234
235~~~~
236n=2^k+rest
237
238Left = 2^k-1 + ((rest>2^k-2) ? 2^k-2 : rest)
239Right = 2^k-1 + ((rest>2^k-2) ? rest-2^k-2 : 0)
240~~~~
241
242For example, let `n_nodes=67`. Then, the closest `2^k=64, 2^k-1=32, 2^k-2=16`.
243Left node gets `32+3=35` sub-nodes, and the right node gets 32 sub-nodes
244
245The division process continues until all the nodes contain not more than a predefined number
246of points.
247
248#### 4c. The order of nodes in boundary-related arrays
249
250Some kd-tree based algorithms need to know the boundaries of each node. This information can
251be computed by calling the TKDTree::MakeBoundaries() function. It fills the following arrays:
252
253- `fRange` : array containing the boundaries of the domain:
254 `| 1st dimension (min + max) | 2nd dimension (min + max) | ...`
255`fBoundaries` : nodes boundaries
256 `| 1st node {1st dim * 2 elements | 2nd dim * 2 elements | ...} | 2nd node {...} | ...`
257
258
259The nodes are arranged in the order described in section 3a.
260
261
262- **Note**: the storage of the TKDTree in a file which include also the contained data is not
263 supported. One must store the data separately in a file (e.g. using a TTree) and then
264 re-creating the TKDTree from the data, after having read them from the file
265
266@ingroup Random
267
268
269*/
270////////////////////////////////////////////////////////////////////////////////
271/// Default constructor. Nothing is built
272
273template <typename Index, typename Value>
275 TObject()
276 ,fDataOwner(kFALSE)
277 ,fNNodes(0)
278 ,fTotalNodes(0)
279 ,fNDim(0)
280 ,fNDimm(0)
281 ,fNPoints(0)
282 ,fBucketSize(0)
283 ,fAxis(nullptr)
284 ,fValue(nullptr)
285 ,fRange(nullptr)
286 ,fData(nullptr)
287 ,fBoundaries(nullptr)
288 ,fIndPoints(nullptr)
289 ,fRowT0(0)
290 ,fCrossNode(0)
291 ,fOffset(0)
292{
293}
294
295template <typename Index, typename Value>
297 TObject()
298 ,fDataOwner(0)
299 ,fNNodes(0)
300 ,fTotalNodes(0)
301 ,fNDim(ndim)
302 ,fNDimm(2*ndim)
303 ,fNPoints(npoints)
304 ,fBucketSize(bsize)
305 ,fAxis(nullptr)
306 ,fValue(nullptr)
307 ,fRange(nullptr)
308 ,fData(nullptr)
309 ,fBoundaries(nullptr)
310 ,fIndPoints(nullptr)
311 ,fRowT0(0)
312 ,fCrossNode(0)
313 ,fOffset(0)
314{
315// Create the kd-tree of npoints from ndim-dimensional space. Parameter bsize stands for the
316// maximal number of points in the terminal nodes (buckets).
317// Proceed by calling one of the SetData() functions and then the Build() function
318// Note, that updating the tree with new data after the Build() function has been called is
319// not possible!
320
321}
322
323////////////////////////////////////////////////////////////////////////////////
324/// Create a kd-tree from the provided data array. This function only sets the data,
325/// call Build() to build the tree!!!
326/// Parameters:
327/// - npoints - total number of points. Adding points after the tree is built is not supported
328/// - ndim - number of dimensions
329/// - bsize - maximal number of points in the terminal nodes
330/// - data - the data array
331///
332/// The data should be placed columnwise (like in a TTree).
333/// The columnwise orientation is chosen to simplify the usage together with TTree::GetV1() like functions.
334/// An example of filling such an array for a 2d case:
335/// Double_t **data = new Double_t*[2];
336/// data[0] = new Double_t[npoints];
337/// data[1] = new Double_t[npoints];
338/// for (Int_t i=0; i<npoints; i++){
339/// data[0][i]=gRandom->Uniform(-1, 1); //fill the x-coordinate
340/// data[1][i]=gRandom->Uniform(-1, 1); //fill the y-coordinate
341/// }
342///
343/// By default, the kd-tree doesn't own the data. If you want the kd-tree to delete the data array, call
344/// kdtree->SetOwner(kTRUE).
345///
346
347template <typename Index, typename Value>
349 TObject()
350 ,fDataOwner(0)
351 ,fNNodes(0)
352 ,fTotalNodes(0)
353 ,fNDim(ndim)
354 ,fNDimm(2*ndim)
355 ,fNPoints(npoints)
356 ,fBucketSize(bsize)
357 ,fAxis(nullptr)
358 ,fValue(nullptr)
359 ,fRange(nullptr)
360 ,fData(data) //Columnwise!!!!!
361 ,fBoundaries(nullptr)
362 ,fIndPoints(nullptr)
363 ,fRowT0(0)
364 ,fCrossNode(0)
365 ,fOffset(0)
366{
367
368 //Build();
369}
370
371////////////////////////////////////////////////////////////////////////////////
372/// Destructor
373/// By default, the original data is not owned by kd-tree and is not deleted with it.
374/// If you want to delete the data along with the kd-tree, call SetOwner(kTRUE).
375
376template <typename Index, typename Value>
378{
379 if (fAxis) delete [] fAxis;
380 if (fValue) delete [] fValue;
381 if (fIndPoints) delete [] fIndPoints;
382 if (fRange) delete [] fRange;
383 if (fBoundaries) delete [] fBoundaries;
384 if (fData) {
385 if (fDataOwner==1){
386 //the tree owns all the data
387 for(int idim=0; idim<fNDim; idim++) delete [] fData[idim];
388 }
389 if (fDataOwner>0) {
390 //the tree owns the array of pointers
391 delete [] fData;
392 }
393 }
394// if (fDataOwner && fData){
395// for(int idim=0; idim<fNDim; idim++) delete [] fData[idim];
396// delete [] fData;
397// }
398}
399
400
401////////////////////////////////////////////////////////////////////////////////
402///
403/// Build the kd-tree
404///
405/// 1. calculate number of nodes
406/// 2. calculate first terminal row
407/// 3. initialize index array
408/// 4. non recursive building of the binary tree
409///
410///
411/// The tree is divided recursively. See class description, section 4b for the details
412/// of the division algorithm
413
414template <typename Index, typename Value>
416{
417 //1.
418 fNNodes = fNPoints/fBucketSize-1;
419 if (fNPoints%fBucketSize) fNNodes++;
420 fTotalNodes = fNNodes + fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);
421 //2.
422 fRowT0=0;
423 for ( ;(fNNodes+1)>(1<<fRowT0);fRowT0++) {}
424 fRowT0-=1;
425 // 2 = 2**0 + 1
426 // 3 = 2**1 + 1
427 // 4 = 2**1 + 2
428 // 5 = 2**2 + 1
429 // 6 = 2**2 + 2
430 // 7 = 2**2 + 3
431 // 8 = 2**2 + 4
432
433 //3.
434 // allocate space for boundaries
435 fRange = new Value[2*fNDim];
436 fIndPoints= new Index[fNPoints];
437 for (Index i=0; i<fNPoints; i++) fIndPoints[i] = i;
438 fAxis = new UChar_t[fNNodes];
439 fValue = new Value[fNNodes];
440 //
441 fCrossNode = (1<<(fRowT0+1))-1;
442 if (fCrossNode<fNNodes) fCrossNode = 2*fCrossNode+1;
443 //
444 // fOffset = (((fNNodes+1)-(1<<fRowT0)))*2;
445 Int_t over = (fNNodes+1)-(1<<fRowT0);
446 Int_t filled = ((1<<fRowT0)-over)*fBucketSize;
447 fOffset = fNPoints-filled;
448
449 //
450 // printf("Row0 %d\n", fRowT0);
451 // printf("CrossNode %d\n", fCrossNode);
452 // printf("Offset %d\n", fOffset);
453 //
454 //
455 //4.
456 // stack for non recursive build - size 128 bytes enough
457 Int_t rowStack[128];
458 Int_t nodeStack[128];
459 Int_t npointStack[128];
460 Int_t posStack[128];
462 rowStack[0] = 0;
463 nodeStack[0] = 0;
464 npointStack[0] = fNPoints;
465 posStack[0] = 0;
466 //
467 while (currentIndex>=0){
468 //
470 if (npoints<=fBucketSize) {
471 currentIndex--;
472 continue; // terminal node
473 }
477 //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
478 //
479 // divide points
480 Int_t nbuckets0 = npoints/fBucketSize; //current number of buckets
481 if (npoints%fBucketSize) nbuckets0++; //
482 Int_t restRows = fRowT0-rowStack[currentIndex]; // rest of fully occupied node row
483 if (restRows<0) restRows =0;
484 for (;nbuckets0>(2<<restRows); restRows++) {}
485 Int_t nfull = 1<<restRows;
487 Int_t nleft =0, nright =0;
488 //
489 if (nrest>(nfull/2)){
490 nleft = nfull*fBucketSize;
492 }else{
493 nright = nfull*fBucketSize/2;
495 }
496
497 //
498 //find the axis with biggest spread
500 Value tempspread, min, max;
501 Index axspread=0;
502 Value *array;
503 for (Int_t idim=0; idim<fNDim; idim++){
504 array = fData[idim];
505 Spread(npoints, array, fIndPoints+cpos, min, max);
506 tempspread = max - min;
507 if (maxspread < tempspread) {
509 axspread = idim;
510 }
511 if(cnode) continue;
512 //printf("set %d %6.3f %6.3f\n", idim, min, max);
513 fRange[2*idim] = min; fRange[2*idim+1] = max;
514 }
515 array = fData[axspread];
516 KOrdStat(npoints, array, nleft, fIndPoints+cpos);
517 fAxis[cnode] = axspread;
518 fValue[cnode] = array[fIndPoints[cpos+nleft]];
519 //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
520 //
521 //
526 currentIndex++;
531 //
532 if (false){
533 // consistency check
534 Info("Build()", "%s", Form("points %d left %d right %d", npoints, nleft, nright));
535 if (nleft<nright) Warning("Build", "Problem Left-Right");
536 if (nleft<0 || nright<0) Warning("Build()", "Problem Negative number");
537 }
538 }
539}
540
541////////////////////////////////////////////////////////////////////////////////
542///Find kNN nearest neighbors to the point in the first argument
543///Returns 1 on success, 0 on failure
544///Arrays ind and dist are provided by the user and are assumed to be at least kNN elements long
545
546template <typename Index, typename Value>
547void TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t kNN, Index *ind, Value *dist)
548{
549
550 if (!ind || !dist) {
551 Error("FindNearestNeighbors", "Working arrays must be allocated by the user!");
552 return;
553 }
554 for (Int_t i=0; i<kNN; i++){
555 dist[i]=std::numeric_limits<Value>::max();
556 ind[i]=-1;
557 }
558 MakeBoundariesExact();
559 UpdateNearestNeighbors(0, point, kNN, ind, dist);
560
561}
562
563////////////////////////////////////////////////////////////////////////////////
564///Update the nearest neighbors values by examining the node inode
565
566template <typename Index, typename Value>
568{
569 Value min=0;
570 Value max=0;
571 DistanceToNode(point, inode, min, max);
572 if (min > dist[kNN-1]){
573 //there are no closer points in this node
574 return;
575 }
576 if (IsTerminal(inode)) {
577 //examine points one by one
578 Index f1, l1, f2, l2;
579 GetNodePointsIndexes(inode, f1, l1, f2, l2);
580 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
581 Double_t d = Distance(point, fIndPoints[ipoint]);
582 if (d<dist[kNN-1]){
583 //found a closer point
584 Int_t ishift=0;
585 while(ishift<kNN && d>dist[ishift])
586 ishift++;
587 //replace the neighbor #ishift with the found point
588 //and shift the rest 1 index value to the right
589 for (Int_t i=kNN-1; i>ishift; i--){
590 dist[i]=dist[i-1];
591 ind[i]=ind[i-1];
592 }
593 dist[ishift]=d;
594 ind[ishift]=fIndPoints[ipoint];
595 }
596 }
597 return;
598 }
599 if (point[fAxis[inode]]<fValue[inode]){
600 //first examine the node that contains the point
601 UpdateNearestNeighbors(GetLeft(inode), point, kNN, ind, dist);
602 UpdateNearestNeighbors(GetRight(inode), point, kNN, ind, dist);
603 } else {
604 UpdateNearestNeighbors(GetRight(inode), point, kNN, ind, dist);
605 UpdateNearestNeighbors(GetLeft(inode), point, kNN, ind, dist);
606 }
607}
608
609////////////////////////////////////////////////////////////////////////////////
610///Find the distance between point of the first argument and the point at index value ind
611///Type argument specifies the metric: type=2 - L2 metric, type=1 - L1 metric
612
613template <typename Index, typename Value>
615{
616 Double_t dist = 0;
617 if (type==2){
618 for (Int_t idim=0; idim<fNDim; idim++){
619 dist+=(point[idim]-fData[idim][ind])*(point[idim]-fData[idim][ind]);
620 }
621 return TMath::Sqrt(dist);
622 } else {
623 for (Int_t idim=0; idim<fNDim; idim++){
624 dist+=TMath::Abs(point[idim]-fData[idim][ind]);
625 }
626
627 return dist;
628 }
629 return -1;
630
631}
632
633////////////////////////////////////////////////////////////////////////////////
634///Find the minimal and maximal distance from a given point to a given node.
635///Type argument specifies the metric: type=2 - L2 metric, type=1 - L1 metric
636///If the point is inside the node, both min and max are set to 0.
637
638template <typename Index, typename Value>
640{
641 Value *bound = GetBoundaryExact(inode);
642 min = 0;
643 max = 0;
645
646 if (type==2){
647 for (Int_t idim=0; idim<fNDimm; idim+=2){
648 dist1 = (point[idim/2]-bound[idim])*(point[idim/2]-bound[idim]);
649 dist2 = (point[idim/2]-bound[idim+1])*(point[idim/2]-bound[idim+1]);
650 //min+=TMath::Min(dist1, dist2);
651 if (point[idim/2]<bound[idim] || point[idim/2]>bound[idim+1])
652 min+= (dist1>dist2)? dist2 : dist1;
653 // max+=TMath::Max(dist1, dist2);
654 max+= (dist1>dist2)? dist1 : dist2;
655 }
656 min = TMath::Sqrt(min);
657 max = TMath::Sqrt(max);
658 } else {
659 for (Int_t idim=0; idim<fNDimm; idim+=2){
660 dist1 = TMath::Abs(point[idim/2]-bound[idim]);
661 dist2 = TMath::Abs(point[idim/2]-bound[idim+1]);
662 //min+=TMath::Min(dist1, dist2);
663 min+= (dist1>dist2)? dist2 : dist1;
664 // max+=TMath::Max(dist1, dist2);
665 max+= (dist1>dist2)? dist1 : dist2;
666 }
667 }
668}
669
670////////////////////////////////////////////////////////////////////////////////
671/// returns the index of the terminal node to which point belongs
672/// (index in the fAxis, fValue, etc arrays)
673/// returns -1 in case of failure
674
675template <typename Index, typename Value>
676Index TKDTree<Index, Value>::FindNode(const Value * point) const
677{
678 Index stackNode[128], inode;
680 stackNode[0] = 0;
681 while (currentIndex>=0){
683 if (IsTerminal(inode)) return inode;
684
685 currentIndex--;
686 if (point[fAxis[inode]]<=fValue[inode]){
687 currentIndex++;
688 stackNode[currentIndex]=(inode<<1)+1; //GetLeft()
689 }
690 if (point[fAxis[inode]]>=fValue[inode]){
691 currentIndex++;
692 stackNode[currentIndex]=(inode+1)<<1; //GetRight()
693 }
694 }
695
696 return -1;
697}
698
699
700
701////////////////////////////////////////////////////////////////////////////////
702///
703/// find the index of point
704/// works only if we keep fData pointers
705
706template <typename Index, typename Value>
708 Int_t stackNode[128];
710 stackNode[0] = 0;
711 iter =0;
712 //
713 while (currentIndex>=0){
714 iter++;
716 currentIndex--;
717 if (IsTerminal(inode)){
718 // investigate terminal node
719 Int_t indexIP = (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNNodes)*fBucketSize+fOffset;
720 printf("terminal %d indexP %d\n", inode, indexIP);
721 for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
722 Bool_t isOK = kTRUE;
724 printf("ibucket %d index %d\n", ibucket, indexIP);
725 if (indexIP>=fNPoints) continue;
726 Int_t index0 = fIndPoints[indexIP];
727 for (Int_t idim=0;idim<fNDim;idim++) if (fData[idim][index0]!=point[idim]) isOK = kFALSE;
728 if (isOK) index = index0;
729 }
730 continue;
731 }
732
733 if (point[fAxis[inode]]<=fValue[inode]){
734 currentIndex++;
736 }
737 if (point[fAxis[inode]]>=fValue[inode]){
738 currentIndex++;
740 }
741 }
742 //
743 // printf("Iter\t%d\n",iter);
744}
745
746////////////////////////////////////////////////////////////////////////////////
747///Find all points in the sphere of a given radius "range" around the given point
748///1st argument - the point
749///2nd argument - radius of the shere
750///3rd argument - a vector, in which the results will be returned
751
752template <typename Index, typename Value>
753void TKDTree<Index, Value>::FindInRange(Value * point, Value range, std::vector<Index> &res)
754{
755 MakeBoundariesExact();
756 UpdateRange(0, point, range, res);
757}
758
759////////////////////////////////////////////////////////////////////////////////
760///Internal recursive function with the implementation of range searches
761
762template <typename Index, typename Value>
763void TKDTree<Index, Value>::UpdateRange(Index inode, Value* point, Value range, std::vector<Index> &res)
764{
765 Value min, max;
766 DistanceToNode(point, inode, min, max);
767 if (min>range) {
768 //all points of this node are outside the range
769 return;
770 }
771 if (max<range && max>0) {
772 //all points of this node are inside the range
773
774 Index f1, l1, f2, l2;
775 GetNodePointsIndexes(inode, f1, l1, f2, l2);
776
777 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
778 res.push_back(fIndPoints[ipoint]);
779 }
780 for (Int_t ipoint=f2; ipoint<=l2; ipoint++){
781 res.push_back(fIndPoints[ipoint]);
782 }
783 return;
784 }
785
786 //this node intersects with the range
787 if (IsTerminal(inode)){
788 //examine the points one by one
789 Index f1, l1, f2, l2;
790 Double_t d;
791 GetNodePointsIndexes(inode, f1, l1, f2, l2);
792 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
793 d = Distance(point, fIndPoints[ipoint]);
794 if (d <= range){
795 res.push_back(fIndPoints[ipoint]);
796 }
797 }
798 return;
799 }
800 if (point[fAxis[inode]]<=fValue[inode]){
801 //first examine the node that contains the point
802 UpdateRange(GetLeft(inode),point, range, res);
803 UpdateRange(GetRight(inode),point, range, res);
804 } else {
805 UpdateRange(GetLeft(inode),point, range, res);
806 UpdateRange(GetRight(inode),point, range, res);
807 }
808}
809
810////////////////////////////////////////////////////////////////////////////////
811///return the indices of the points in that terminal node
812///for all the nodes except last, the size is fBucketSize
813///for the last node it's fOffset%fBucketSize
814
815template <typename Index, typename Value>
817{
818 if (!IsTerminal(node)){
819 printf("GetPointsIndexes() only for terminal nodes, use GetNodePointsIndexes() instead\n");
820 return nullptr;
821 }
822 Int_t offset = (node >= fCrossNode) ? (node-fCrossNode)*fBucketSize : fOffset+(node-fNNodes)*fBucketSize;
823 return &fIndPoints[offset];
824}
825
826////////////////////////////////////////////////////////////////////////////////
827///Return the indices of points in that node
828///Indices are returned as the first and last value of the part of indices array, that belong to this node
829///Sometimes points are in 2 intervals, then the first and last value for the second one are returned in
830///third and fourth parameter, otherwise first2 is set to 0 and last2 is set to -1
831///To iterate over all the points of the node `#inode`, one can do, for example:
832///Index *indices = kdtree->GetPointsIndexes();
833///Int_t first1, last1, first2, last2;
834///kdtree->GetPointsIndexes(inode, first1, last1, first2, last2);
835///for (Int_t ipoint=first1; ipoint<=last1; ipoint++){
836/// point = indices[ipoint];
837/// //do something with point;
838///}
839///for (Int_t ipoint=first2; ipoint<=last2; ipoint++){
840/// point = indices[ipoint];
841/// //do something with point;
842///}
843
844template <typename Index, typename Value>
846{
847
848 if (IsTerminal(node)){
849 //the first point in the node is computed by the following formula:
850 Index offset = (node >= fCrossNode) ? (node-fCrossNode)*fBucketSize : fOffset+(node-fNNodes)*fBucketSize;
851 first1 = offset;
852 last1 = offset + GetNPointsNode(node)-1;
853 first2 = 0;
854 last2 = -1;
855 return;
856 }
857
858 Index firsttermnode = fNNodes;
859 Index ileft = node;
860 Index iright = node;
861 Index f1, l1, f2, l2;
862//this is the left-most node
863 while (ileft<firsttermnode)
864 ileft = GetLeft(ileft);
865//this is the right-most node
866 while (iright<firsttermnode)
867 iright = GetRight(iright);
868
869 if (ileft>iright){
870// first1 = firsttermnode;
871// last1 = iright;
872// first2 = ileft;
873// last2 = fTotalNodes-1;
874 GetNodePointsIndexes(firsttermnode, f1, l1, f2, l2);
875 first1 = f1;
876 GetNodePointsIndexes(iright, f1, l1, f2, l2);
877 last1 = l1;
878 GetNodePointsIndexes(ileft, f1, l1, f2, l2);
879 first2 = f1;
880 GetNodePointsIndexes(fTotalNodes-1, f1, l1, f2, l2);
881 last2 = l1;
882
883 } else {
884 GetNodePointsIndexes(ileft, f1, l1, f2, l2);
885 first1 = f1;
886 GetNodePointsIndexes(iright, f1, l1, f2, l2);
887 last1 = l1;
888 first2 = 0;
889 last2 = -1;
890 }
891}
892
893////////////////////////////////////////////////////////////////////////////////
894///Get number of points in this node
895///for all the terminal nodes except last, the size is fBucketSize
896///for the last node it's fOffset%fBucketSize, or if fOffset%fBucketSize==0, it's also fBucketSize
897
898template <typename Index, typename Value>
900{
901 if (IsTerminal(inode)){
902
903 if (inode!=fTotalNodes-1) return fBucketSize;
904 else {
905 if (fOffset%fBucketSize==0) return fBucketSize;
906 else return fOffset%fBucketSize;
907 }
908 }
909
910 Int_t f1, l1, f2, l2;
911 GetNodePointsIndexes(inode, f1, l1, f2, l2);
912 Int_t sum = l1-f1+1;
913 sum += l2-f2+1;
914 return sum;
915}
916
917
918////////////////////////////////////////////////////////////////////////////////
919/// Set the data array. See the constructor function comments for details
920
921template <typename Index, typename Value>
923{
924// TO DO
925//
926// Check reconstruction/reallocation of memory of data. Maybe it is not
927// necessary to delete and realocate space but only to use the same space
928
929 Clear();
930
931 //Columnwise!!!!
932 fData = data;
933 fNPoints = npoints;
934 fNDim = ndim;
935 fBucketSize = bsize;
936
937 Build();
938}
939
940////////////////////////////////////////////////////////////////////////////////
941///Set the coordinate `#ndim` of all points (the column `#ndim` of the data matrix)
942///After setting all the data columns, proceed by calling Build() function
943///Note, that calling this function after Build() is not possible
944///Note also, that no checks on the array sizes is performed anywhere
945
946template <typename Index, typename Value>
948{
949 if (fAxis || fValue) {
950 Error("SetData", "The tree has already been built, no updates possible");
951 return 0;
952 }
953
954 if (!fData) {
955 fData = new Value*[fNDim];
956 }
957 fData[idim]=data;
958 fDataOwner = 2;
959 return 1;
960}
961
962
963////////////////////////////////////////////////////////////////////////////////
964///Calculate spread of the array a
965
966template <typename Index, typename Value>
967void TKDTree<Index, Value>::Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const
968{
969 Index i;
970 min = a[index[0]];
971 max = a[index[0]];
972 for (i=0; i<ntotal; i++){
973 if (a[index[i]]<min) min = a[index[i]];
974 if (a[index[i]]>max) max = a[index[i]];
975 }
976}
977
978
979////////////////////////////////////////////////////////////////////////////////
980///
981///copy of the TMath::KOrdStat because I need an Index work array
982
983template <typename Index, typename Value>
985{
986 Index i, ir, j, l, mid;
987 Index arr;
988 Index temp;
989
990 Index rk = k;
991 l=0;
992 ir = ntotal-1;
993 for(;;) {
994 if (ir<=l+1) { //active partition contains 1 or 2 elements
995 if (ir == l+1 && a[index[ir]]<a[index[l]])
996 {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
997 Value tmp = a[index[rk]];
998 return tmp;
999 } else {
1000 mid = (l+ir) >> 1; //choose median of left, center and right
1001 {temp = index[mid]; index[mid]=index[l+1]; index[l+1]=temp;}//elements as partitioning element arr.
1002 if (a[index[l]]>a[index[ir]]) //also rearrange so that a[l]<=a[l+1]
1003 {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
1004
1005 if (a[index[l+1]]>a[index[ir]])
1006 {temp=index[l+1]; index[l+1]=index[ir]; index[ir]=temp;}
1007
1008 if (a[index[l]]>a[index[l+1]])
1009 {temp = index[l]; index[l]=index[l+1]; index[l+1]=temp;}
1010
1011 i=l+1; //initialize pointers for partitioning
1012 j=ir;
1013 arr = index[l+1];
1014 for (;;) {
1015 do i++; while (a[index[i]]<a[arr]);
1016 do j--; while (a[index[j]]>a[arr]);
1017 if (j<i) break; //pointers crossed, partitioning complete
1018 {temp=index[i]; index[i]=index[j]; index[j]=temp;}
1019 }
1020 index[l+1]=index[j];
1021 index[j]=arr;
1022 if (j>=rk) ir = j-1; //keep active the partition that
1023 if (j<=rk) l=i; //contains the k_th element
1024 }
1025 }
1026}
1027
1028////////////////////////////////////////////////////////////////////////////////
1029/// Build boundaries for each node. Note, that the boundaries here are built
1030/// based on the splitting planes of the kd-tree, and don't necessarily pass
1031/// through the points of the original dataset. For the latter functionality
1032/// see function MakeBoundariesExact()
1033/// Boundaries can be retrieved by calling GetBoundary(inode) function that would
1034/// return an array of boundaries for the specified node, or GetBoundaries() function
1035/// that would return the complete array.
1036
1037template <typename Index, typename Value>
1039{
1040
1041 if(range) memcpy(fRange, range, fNDimm*sizeof(Value));
1042 // total number of nodes including terminal nodes
1043 Int_t totNodes = fNNodes + fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);
1044 fBoundaries = new Value[totNodes*fNDimm];
1045 //Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes));
1046
1047
1048 // loop
1049 Value *tbounds = nullptr, *cbounds = nullptr;
1050 Int_t cn;
1051 for(int inode=fNNodes-1; inode>=0; inode--){
1052 tbounds = &fBoundaries[inode*fNDimm];
1053 memcpy(tbounds, fRange, fNDimm*sizeof(Value));
1054
1055 // check left child node
1056 cn = (inode<<1)+1;
1057 if(IsTerminal(cn)) CookBoundaries(inode, kTRUE);
1058 cbounds = &fBoundaries[fNDimm*cn];
1059 for(int idim=0; idim<fNDim; idim++) tbounds[idim<<1] = cbounds[idim<<1];
1060
1061 // check right child node
1062 cn = (inode+1)<<1;
1063 if(IsTerminal(cn)) CookBoundaries(inode, kFALSE);
1064 cbounds = &fBoundaries[fNDimm*cn];
1065 for(int idim=0; idim<fNDim; idim++) tbounds[(idim<<1)+1] = cbounds[(idim<<1)+1];
1066 }
1067}
1068
1069////////////////////////////////////////////////////////////////////////////////
1070/// define index of this terminal node
1071
1072template <typename Index, typename Value>
1074{
1075 Int_t index = (node<<1) + (LEFT ? 1 : 2);
1076 //Info("CookBoundaries()", Form("Node %d", index));
1077
1078 // build and initialize boundaries for this node
1079 Value *tbounds = &fBoundaries[fNDimm*index];
1080 memcpy(tbounds, fRange, fNDimm*sizeof(Value));
1081 Bool_t flag[256]; // cope with up to 128 dimensions
1082 memset(flag, kFALSE, fNDimm);
1083 Int_t nvals = 0;
1084
1085 // recurse parent nodes
1086 Int_t pn = node;
1087 while(pn >= 0 && nvals < fNDimm){
1088 if(LEFT){
1089 index = (fAxis[pn]<<1)+1;
1090 if(!flag[index]) {
1091 tbounds[index] = fValue[pn];
1092 flag[index] = kTRUE;
1093 nvals++;
1094 }
1095 } else {
1096 index = fAxis[pn]<<1;
1097 if(!flag[index]) {
1098 tbounds[index] = fValue[pn];
1099 flag[index] = kTRUE;
1100 nvals++;
1101 }
1102 }
1103 LEFT = pn&1;
1104 pn = (pn - 1)>>1;
1105 }
1106}
1107
1108////////////////////////////////////////////////////////////////////////////////
1109/// Build boundaries for each node. Unlike MakeBoundaries() function
1110/// the boundaries built here always pass through a point of the original dataset
1111/// So, for example, for a terminal node with just one point minimum and maximum for each
1112/// dimension are the same.
1113/// Boundaries can be retrieved by calling GetBoundaryExact(inode) function that would
1114/// return an array of boundaries for the specified node, or GetBoundaries() function
1115/// that would return the complete array.
1116
1117template <typename Index, typename Value>
1119{
1120
1121 // total number of nodes including terminal nodes
1122 //Int_t totNodes = fNNodes + fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);
1123 if (fBoundaries){
1124 //boundaries were already computed for this tree
1125 return;
1126 }
1127 fBoundaries = new Value[fTotalNodes*fNDimm];
1128 Value *min = new Value[fNDim];
1129 Value *max = new Value[fNDim];
1130 for (Index inode=fNNodes; inode<fTotalNodes; inode++){
1131 //go through the terminal nodes
1132 for (Index idim=0; idim<fNDim; idim++){
1133 min[idim]= std::numeric_limits<Value>::max();
1134 max[idim]=-std::numeric_limits<Value>::max();
1135 }
1136 Index *points = GetPointsIndexes(inode);
1137 Index npoints = GetNPointsNode(inode);
1138 //find max and min in each dimension
1139 for (Index ipoint=0; ipoint<npoints; ipoint++){
1140 for (Index idim=0; idim<fNDim; idim++){
1141 if (fData[idim][points[ipoint]]<min[idim])
1142 min[idim]=fData[idim][points[ipoint]];
1143 if (fData[idim][points[ipoint]]>max[idim])
1144 max[idim]=fData[idim][points[ipoint]];
1145 }
1146 }
1147 for (Index idim=0; idim<fNDimm; idim+=2){
1148 fBoundaries[inode*fNDimm + idim]=min[idim/2];
1149 fBoundaries[inode*fNDimm + idim+1]=max[idim/2];
1150 }
1151 }
1152
1153 delete [] min;
1154 delete [] max;
1155
1156 Index left, right;
1157 for (Index inode=fNNodes-1; inode>=0; inode--){
1158 //take the min and max of left and right
1159 left = GetLeft(inode)*fNDimm;
1160 right = GetRight(inode)*fNDimm;
1161 for (Index idim=0; idim<fNDimm; idim+=2){
1162 //take the minimum on each dimension
1163 fBoundaries[inode*fNDimm+idim]=TMath::Min(fBoundaries[left+idim], fBoundaries[right+idim]);
1164 //take the maximum on each dimension
1165 fBoundaries[inode*fNDimm+idim+1]=TMath::Max(fBoundaries[left+idim+1], fBoundaries[right+idim+1]);
1166
1167 }
1168 }
1169}
1170
1171////////////////////////////////////////////////////////////////////////////////
1172///
1173/// find the smallest node covering the full range - start
1174///
1175
1176template <typename Index, typename Value>
1178 inode =0;
1179 for (;inode<fNNodes;){
1180 if (TMath::Abs(point[fAxis[inode]] - fValue[inode])<delta[fAxis[inode]]) break;
1181 inode = (point[fAxis[inode]] < fValue[inode]) ? (inode*2)+1: (inode*2)+2;
1182 }
1183}
1184
1185////////////////////////////////////////////////////////////////////////////////
1186/// Get the boundaries
1187
1188template <typename Index, typename Value>
1190{
1191 if(!fBoundaries) MakeBoundaries();
1192 return fBoundaries;
1193}
1194
1195
1196////////////////////////////////////////////////////////////////////////////////
1197/// Get the boundaries
1198
1199template <typename Index, typename Value>
1201{
1202 if(!fBoundaries) MakeBoundariesExact();
1203 return fBoundaries;
1204}
1205
1206////////////////////////////////////////////////////////////////////////////////
1207/// Get a boundary
1208
1209template <typename Index, typename Value>
1211{
1212 if(!fBoundaries) MakeBoundaries();
1213 return &fBoundaries[node*2*fNDim];
1214}
1215
1216////////////////////////////////////////////////////////////////////////////////
1217/// Get a boundary
1218
1219template <typename Index, typename Value>
1221{
1222 if(!fBoundaries) MakeBoundariesExact();
1223 return &fBoundaries[node*2*fNDim];
1224}
1225
1226
1227////////////////////////////////////////////////////////////////////////////////
1228///
1229/// Example function to
1230///
1231
1233 Float_t *data0 = new Float_t[npoints*2];
1234 Float_t *data[2];
1235 data[0] = &data0[0];
1236 data[1] = &data0[npoints];
1237 for (Int_t i=0;i<npoints;i++) {
1238 data[1][i]= gRandom->Rndm();
1239 data[0][i]= gRandom->Rndm();
1240 }
1242 return kdtree;
1243}
1244
1245
1246
1247template class TKDTree<Int_t, Float_t>;
1248template class TKDTree<Int_t, Double_t>;
#define d(i)
Definition RSha256.hxx:102
#define a(i)
Definition RSha256.hxx:99
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
unsigned char UChar_t
Unsigned Character 1 byte (unsigned char)
Definition RtypesCore.h:52
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
#define templateClassImp(name)
Definition Rtypes.h:405
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition TError.cxx:241
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:252
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h offset
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t points
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
TKDTreeIF * TKDTreeTestBuild()
TKDTree< Int_t, Float_t > TKDTreeIF
Definition TKDTree.h:105
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2496
Class implementing a kd-tree.
Definition TKDTree.h:10
void FindPoint(Value *point, Index &index, Int_t &iter)
find the index of point works only if we keep fData pointers
Definition TKDTree.cxx:707
void SetData(Index npoints, Index ndim, UInt_t bsize, Value **data)
Set the data array. See the constructor function comments for details.
Definition TKDTree.cxx:922
~TKDTree() override
Destructor By default, the original data is not owned by kd-tree and is not deleted with it.
Definition TKDTree.cxx:377
void GetNodePointsIndexes(Int_t node, Int_t &first1, Int_t &last1, Int_t &first2, Int_t &last2) const
Return the indices of points in that node Indices are returned as the first and last value of the par...
Definition TKDTree.cxx:845
Value * GetBoundaryExact(const Int_t node)
Get a boundary.
Definition TKDTree.cxx:1220
void FindBNodeA(Value *point, Value *delta, Int_t &inode)
find the smallest node covering the full range - start
Definition TKDTree.cxx:1177
void UpdateRange(Index inode, Value *point, Value range, std::vector< Index > &res)
Internal recursive function with the implementation of range searches.
Definition TKDTree.cxx:763
Index FindNode(const Value *point) const
returns the index of the terminal node to which point belongs (index in the fAxis,...
Definition TKDTree.cxx:676
Value KOrdStat(Index ntotal, Value *a, Index k, Index *index) const
copy of the TMath::KOrdStat because I need an Index work array
Definition TKDTree.cxx:984
Value * GetBoundariesExact()
Get the boundaries.
Definition TKDTree.cxx:1200
void MakeBoundariesExact()
Build boundaries for each node.
Definition TKDTree.cxx:1118
void CookBoundaries(const Int_t node, Bool_t left)
define index of this terminal node
Definition TKDTree.cxx:1073
void MakeBoundaries(Value *range=nullptr)
Build boundaries for each node.
Definition TKDTree.cxx:1038
void UpdateNearestNeighbors(Index inode, const Value *point, Int_t kNN, Index *ind, Value *dist)
Update the nearest neighbors values by examining the node inode.
Definition TKDTree.cxx:567
Double_t Distance(const Value *point, Index ind, Int_t type=2) const
Find the distance between point of the first argument and the point at index value ind Type argument ...
Definition TKDTree.cxx:614
Index * GetPointsIndexes(Int_t node) const
return the indices of the points in that terminal node for all the nodes except last,...
Definition TKDTree.cxx:816
void Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const
Calculate spread of the array a.
Definition TKDTree.cxx:967
void Build()
Build the kd-tree.
Definition TKDTree.cxx:415
void DistanceToNode(const Value *point, Index inode, Value &min, Value &max, Int_t type=2)
Find the minimal and maximal distance from a given point to a given node.
Definition TKDTree.cxx:639
TKDTree()
Default constructor. Nothing is built.
Definition TKDTree.cxx:274
void FindInRange(Value *point, Value range, std::vector< Index > &res)
Find all points in the sphere of a given radius "range" around the given point 1st argument - the poi...
Definition TKDTree.cxx:753
Index GetNPointsNode(Int_t node) const
Get number of points in this node for all the terminal nodes except last, the size is fBucketSize for...
Definition TKDTree.cxx:899
Value * GetBoundaries()
Get the boundaries.
Definition TKDTree.cxx:1189
Value * GetBoundary(const Int_t node)
Get a boundary.
Definition TKDTree.cxx:1210
void FindNearestNeighbors(const Value *point, Int_t k, Index *ind, Value *dist)
Find kNN nearest neighbors to the point in the first argument Returns 1 on success,...
Definition TKDTree.cxx:547
kNN::Event describes point in input variable vector-space, with additional functionality like distanc...
Mother of all ROOT objects.
Definition TObject.h:41
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:559
TF1 * f1
Definition legend1.C:11
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:668
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
TLine l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2340