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