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