Logo ROOT  
Reference Guide
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 <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 developped 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~~~~
181fAxix[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 :
195`fNNodes = (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) | ...`
251`fBoundaries` : 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 separatly 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(0x0)
280 ,fValue(0x0)
281 ,fRange(0x0)
282 ,fData(0x0)
283 ,fBoundaries(0x0)
284 ,fIndPoints(0x0)
285 ,fRowT0(0)
286 ,fCrossNode(0)
287 ,fOffset(0)
288{
289}
290
291template <typename Index, typename Value>
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(0x0)
302 ,fValue(0x0)
303 ,fRange(0x0)
304 ,fData(0x0)
305 ,fBoundaries(0x0)
306 ,fIndPoints(0x0)
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/// Parameteres:
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>
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(0x0)
354 ,fValue(0x0)
355 ,fRange(0x0)
356 ,fData(data) //Columnwise!!!!!
357 ,fBoundaries(0x0)
358 ,fIndPoints(0x0)
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 alogrithm
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 Int_t iter =0;
459 rowStack[0] = 0;
460 nodeStack[0] = 0;
461 npointStack[0] = fNPoints;
462 posStack[0] = 0;
463 //
464 Int_t nbucketsall =0;
465 while (currentIndex>=0){
466 iter++;
467 //
468 Int_t npoints = npointStack[currentIndex];
469 if (npoints<=fBucketSize) {
470 //printf("terminal node : index %d iter %d\n", currentIndex, iter);
471 currentIndex--;
472 nbucketsall++;
473 continue; // terminal node
474 }
475 Int_t crow = rowStack[currentIndex];
476 Int_t cpos = posStack[currentIndex];
477 Int_t cnode = nodeStack[currentIndex];
478 //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
479 //
480 // divide points
481 Int_t nbuckets0 = npoints/fBucketSize; //current number of buckets
482 if (npoints%fBucketSize) nbuckets0++; //
483 Int_t restRows = fRowT0-rowStack[currentIndex]; // rest of fully occupied node row
484 if (restRows<0) restRows =0;
485 for (;nbuckets0>(2<<restRows); restRows++) {}
486 Int_t nfull = 1<<restRows;
487 Int_t nrest = nbuckets0-nfull;
488 Int_t nleft =0, nright =0;
489 //
490 if (nrest>(nfull/2)){
491 nleft = nfull*fBucketSize;
492 nright = npoints-nleft;
493 }else{
494 nright = nfull*fBucketSize/2;
495 nleft = npoints-nright;
496 }
497
498 //
499 //find the axis with biggest spread
500 Value maxspread=0;
501 Value tempspread, min, max;
502 Index axspread=0;
503 Value *array;
504 for (Int_t idim=0; idim<fNDim; idim++){
505 array = fData[idim];
506 Spread(npoints, array, fIndPoints+cpos, min, max);
507 tempspread = max - min;
508 if (maxspread < tempspread) {
509 maxspread=tempspread;
510 axspread = idim;
511 }
512 if(cnode) continue;
513 //printf("set %d %6.3f %6.3f\n", idim, min, max);
514 fRange[2*idim] = min; fRange[2*idim+1] = max;
515 }
516 array = fData[axspread];
517 KOrdStat(npoints, array, nleft, fIndPoints+cpos);
518 fAxis[cnode] = axspread;
519 fValue[cnode] = array[fIndPoints[cpos+nleft]];
520 //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
521 //
522 //
523 npointStack[currentIndex] = nleft;
524 rowStack[currentIndex] = crow+1;
525 posStack[currentIndex] = cpos;
526 nodeStack[currentIndex] = cnode*2+1;
527 currentIndex++;
528 npointStack[currentIndex] = nright;
529 rowStack[currentIndex] = crow+1;
530 posStack[currentIndex] = cpos+nleft;
531 nodeStack[currentIndex] = (cnode*2)+2;
532 //
533 if (0){
534 // consistency check
535 Info("Build()", "%s", Form("points %d left %d right %d", npoints, nleft, nright));
536 if (nleft<nright) Warning("Build", "Problem Left-Right");
537 if (nleft<0 || nright<0) Warning("Build()", "Problem Negative number");
538 }
539 }
540}
541
542////////////////////////////////////////////////////////////////////////////////
543///Find kNN nearest neighbors to the point in the first argument
544///Returns 1 on success, 0 on failure
545///Arrays ind and dist are provided by the user and are assumed to be at least kNN elements long
546
547template <typename Index, typename Value>
549{
550
551 if (!ind || !dist) {
552 Error("FindNearestNeighbors", "Working arrays must be allocated by the user!");
553 return;
554 }
555 for (Int_t i=0; i<kNN; i++){
556 dist[i]=std::numeric_limits<Value>::max();
557 ind[i]=-1;
558 }
559 MakeBoundariesExact();
560 UpdateNearestNeighbors(0, point, kNN, ind, dist);
561
562}
563
564////////////////////////////////////////////////////////////////////////////////
565///Update the nearest neighbors values by examining the node inode
566
567template <typename Index, typename Value>
569{
570 Value min=0;
571 Value max=0;
572 DistanceToNode(point, inode, min, max);
573 if (min > dist[kNN-1]){
574 //there are no closer points in this node
575 return;
576 }
577 if (IsTerminal(inode)) {
578 //examine points one by one
579 Index f1, l1, f2, l2;
580 GetNodePointsIndexes(inode, f1, l1, f2, l2);
581 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
582 Double_t d = Distance(point, fIndPoints[ipoint]);
583 if (d<dist[kNN-1]){
584 //found a closer point
585 Int_t ishift=0;
586 while(ishift<kNN && d>dist[ishift])
587 ishift++;
588 //replace the neighbor #ishift with the found point
589 //and shift the rest 1 index value to the right
590 for (Int_t i=kNN-1; i>ishift; i--){
591 dist[i]=dist[i-1];
592 ind[i]=ind[i-1];
593 }
594 dist[ishift]=d;
595 ind[ishift]=fIndPoints[ipoint];
596 }
597 }
598 return;
599 }
600 if (point[fAxis[inode]]<fValue[inode]){
601 //first examine the node that contains the point
602 UpdateNearestNeighbors(GetLeft(inode), point, kNN, ind, dist);
603 UpdateNearestNeighbors(GetRight(inode), point, kNN, ind, dist);
604 } else {
605 UpdateNearestNeighbors(GetRight(inode), point, kNN, ind, dist);
606 UpdateNearestNeighbors(GetLeft(inode), point, kNN, ind, dist);
607 }
608}
609
610////////////////////////////////////////////////////////////////////////////////
611///Find the distance between point of the first argument and the point at index value ind
612///Type argument specifies the metric: type=2 - L2 metric, type=1 - L1 metric
613
614template <typename Index, typename Value>
616{
617 Double_t dist = 0;
618 if (type==2){
619 for (Int_t idim=0; idim<fNDim; idim++){
620 dist+=(point[idim]-fData[idim][ind])*(point[idim]-fData[idim][ind]);
621 }
622 return TMath::Sqrt(dist);
623 } else {
624 for (Int_t idim=0; idim<fNDim; idim++){
625 dist+=TMath::Abs(point[idim]-fData[idim][ind]);
626 }
627
628 return dist;
629 }
630 return -1;
631
632}
633
634////////////////////////////////////////////////////////////////////////////////
635///Find the minimal and maximal distance from a given point to a given node.
636///Type argument specifies the metric: type=2 - L2 metric, type=1 - L1 metric
637///If the point is inside the node, both min and max are set to 0.
638
639template <typename Index, typename Value>
641{
642 Value *bound = GetBoundaryExact(inode);
643 min = 0;
644 max = 0;
645 Double_t dist1, dist2;
646
647 if (type==2){
648 for (Int_t idim=0; idim<fNDimm; idim+=2){
649 dist1 = (point[idim/2]-bound[idim])*(point[idim/2]-bound[idim]);
650 dist2 = (point[idim/2]-bound[idim+1])*(point[idim/2]-bound[idim+1]);
651 //min+=TMath::Min(dist1, dist2);
652 if (point[idim/2]<bound[idim] || point[idim/2]>bound[idim+1])
653 min+= (dist1>dist2)? dist2 : dist1;
654 // max+=TMath::Max(dist1, dist2);
655 max+= (dist1>dist2)? dist1 : dist2;
656 }
657 min = TMath::Sqrt(min);
658 max = TMath::Sqrt(max);
659 } else {
660 for (Int_t idim=0; idim<fNDimm; idim+=2){
661 dist1 = TMath::Abs(point[idim/2]-bound[idim]);
662 dist2 = TMath::Abs(point[idim/2]-bound[idim+1]);
663 //min+=TMath::Min(dist1, dist2);
664 min+= (dist1>dist2)? dist2 : dist1;
665 // max+=TMath::Max(dist1, dist2);
666 max+= (dist1>dist2)? dist1 : dist2;
667 }
668 }
669}
670
671////////////////////////////////////////////////////////////////////////////////
672/// returns the index of the terminal node to which point belongs
673/// (index in the fAxis, fValue, etc arrays)
674/// returns -1 in case of failure
675
676template <typename Index, typename Value>
678{
679 Index stackNode[128], inode;
680 Int_t currentIndex =0;
681 stackNode[0] = 0;
682 while (currentIndex>=0){
683 inode = stackNode[currentIndex];
684 if (IsTerminal(inode)) return inode;
685
686 currentIndex--;
687 if (point[fAxis[inode]]<=fValue[inode]){
688 currentIndex++;
689 stackNode[currentIndex]=(inode<<1)+1; //GetLeft()
690 }
691 if (point[fAxis[inode]]>=fValue[inode]){
692 currentIndex++;
693 stackNode[currentIndex]=(inode+1)<<1; //GetRight()
694 }
695 }
696
697 return -1;
698}
699
700
701
702////////////////////////////////////////////////////////////////////////////////
703///
704/// find the index of point
705/// works only if we keep fData pointers
706
707template <typename Index, typename Value>
709 Int_t stackNode[128];
710 Int_t currentIndex =0;
711 stackNode[0] = 0;
712 iter =0;
713 //
714 while (currentIndex>=0){
715 iter++;
716 Int_t inode = stackNode[currentIndex];
717 currentIndex--;
718 if (IsTerminal(inode)){
719 // investigate terminal node
720 Int_t indexIP = (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNNodes)*fBucketSize+fOffset;
721 printf("terminal %d indexP %d\n", inode, indexIP);
722 for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
723 Bool_t isOK = kTRUE;
724 indexIP+=ibucket;
725 printf("ibucket %d index %d\n", ibucket, indexIP);
726 if (indexIP>=fNPoints) continue;
727 Int_t index0 = fIndPoints[indexIP];
728 for (Int_t idim=0;idim<fNDim;idim++) if (fData[idim][index0]!=point[idim]) isOK = kFALSE;
729 if (isOK) index = index0;
730 }
731 continue;
732 }
733
734 if (point[fAxis[inode]]<=fValue[inode]){
735 currentIndex++;
736 stackNode[currentIndex]=(inode*2)+1;
737 }
738 if (point[fAxis[inode]]>=fValue[inode]){
739 currentIndex++;
740 stackNode[currentIndex]=(inode*2)+2;
741 }
742 }
743 //
744 // printf("Iter\t%d\n",iter);
745}
746
747////////////////////////////////////////////////////////////////////////////////
748///Find all points in the sphere of a given radius "range" around the given point
749///1st argument - the point
750///2nd argument - radius of the shere
751///3rd argument - a vector, in which the results will be returned
752
753template <typename Index, typename Value>
754void TKDTree<Index, Value>::FindInRange(Value * point, Value range, std::vector<Index> &res)
755{
756 MakeBoundariesExact();
757 UpdateRange(0, point, range, res);
758}
759
760////////////////////////////////////////////////////////////////////////////////
761///Internal recursive function with the implementation of range searches
762
763template <typename Index, typename Value>
764void TKDTree<Index, Value>::UpdateRange(Index inode, Value* point, Value range, std::vector<Index> &res)
765{
766 Value min, max;
767 DistanceToNode(point, inode, min, max);
768 if (min>range) {
769 //all points of this node are outside the range
770 return;
771 }
772 if (max<range && max>0) {
773 //all points of this node are inside the range
774
775 Index f1, l1, f2, l2;
776 GetNodePointsIndexes(inode, f1, l1, f2, l2);
777
778 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
779 res.push_back(fIndPoints[ipoint]);
780 }
781 for (Int_t ipoint=f2; ipoint<=l2; ipoint++){
782 res.push_back(fIndPoints[ipoint]);
783 }
784 return;
785 }
786
787 //this node intersects with the range
788 if (IsTerminal(inode)){
789 //examine the points one by one
790 Index f1, l1, f2, l2;
791 Double_t d;
792 GetNodePointsIndexes(inode, f1, l1, f2, l2);
793 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
794 d = Distance(point, fIndPoints[ipoint]);
795 if (d <= range){
796 res.push_back(fIndPoints[ipoint]);
797 }
798 }
799 return;
800 }
801 if (point[fAxis[inode]]<=fValue[inode]){
802 //first examine the node that contains the point
803 UpdateRange(GetLeft(inode),point, range, res);
804 UpdateRange(GetRight(inode),point, range, res);
805 } else {
806 UpdateRange(GetLeft(inode),point, range, res);
807 UpdateRange(GetRight(inode),point, range, res);
808 }
809}
810
811////////////////////////////////////////////////////////////////////////////////
812///return the indices of the points in that terminal node
813///for all the nodes except last, the size is fBucketSize
814///for the last node it's fOffset%fBucketSize
815
816template <typename Index, typename Value>
818{
819 if (!IsTerminal(node)){
820 printf("GetPointsIndexes() only for terminal nodes, use GetNodePointsIndexes() instead\n");
821 return 0;
822 }
823 Int_t offset = (node >= fCrossNode) ? (node-fCrossNode)*fBucketSize : fOffset+(node-fNNodes)*fBucketSize;
824 return &fIndPoints[offset];
825}
826
827////////////////////////////////////////////////////////////////////////////////
828///Return the indices of points in that node
829///Indices are returned as the first and last value of the part of indices array, that belong to this node
830///Sometimes points are in 2 intervals, then the first and last value for the second one are returned in
831///third and fourth parameter, otherwise first2 is set to 0 and last2 is set to -1
832///To iterate over all the points of the node `#inode`, one can do, for example:
833///Index *indices = kdtree->GetPointsIndexes();
834///Int_t first1, last1, first2, last2;
835///kdtree->GetPointsIndexes(inode, first1, last1, first2, last2);
836///for (Int_t ipoint=first1; ipoint<=last1; ipoint++){
837/// point = indices[ipoint];
838/// //do something with point;
839///}
840///for (Int_t ipoint=first2; ipoint<=last2; ipoint++){
841/// point = indices[ipoint];
842/// //do something with point;
843///}
844
845template <typename Index, typename Value>
846void TKDTree<Index, Value>::GetNodePointsIndexes(Int_t node, Int_t &first1, Int_t &last1, Int_t &first2, Int_t &last2) const
847{
848
849 if (IsTerminal(node)){
850 //the first point in the node is computed by the following formula:
851 Index offset = (node >= fCrossNode) ? (node-fCrossNode)*fBucketSize : fOffset+(node-fNNodes)*fBucketSize;
852 first1 = offset;
853 last1 = offset + GetNPointsNode(node)-1;
854 first2 = 0;
855 last2 = -1;
856 return;
857 }
858
859 Index firsttermnode = fNNodes;
860 Index ileft = node;
861 Index iright = node;
862 Index f1, l1, f2, l2;
863//this is the left-most node
864 while (ileft<firsttermnode)
865 ileft = GetLeft(ileft);
866//this is the right-most node
867 while (iright<firsttermnode)
868 iright = GetRight(iright);
869
870 if (ileft>iright){
871// first1 = firsttermnode;
872// last1 = iright;
873// first2 = ileft;
874// last2 = fTotalNodes-1;
875 GetNodePointsIndexes(firsttermnode, f1, l1, f2, l2);
876 first1 = f1;
877 GetNodePointsIndexes(iright, f1, l1, f2, l2);
878 last1 = l1;
879 GetNodePointsIndexes(ileft, f1, l1, f2, l2);
880 first2 = f1;
881 GetNodePointsIndexes(fTotalNodes-1, f1, l1, f2, l2);
882 last2 = l1;
883
884 } else {
885 GetNodePointsIndexes(ileft, f1, l1, f2, l2);
886 first1 = f1;
887 GetNodePointsIndexes(iright, f1, l1, f2, l2);
888 last1 = l1;
889 first2 = 0;
890 last2 = -1;
891 }
892}
893
894////////////////////////////////////////////////////////////////////////////////
895///Get number of points in this node
896///for all the terminal nodes except last, the size is fBucketSize
897///for the last node it's fOffset%fBucketSize, or if fOffset%fBucketSize==0, it's also fBucketSize
898
899template <typename Index, typename Value>
901{
902 if (IsTerminal(inode)){
903
904 if (inode!=fTotalNodes-1) return fBucketSize;
905 else {
906 if (fOffset%fBucketSize==0) return fBucketSize;
907 else return fOffset%fBucketSize;
908 }
909 }
910
911 Int_t f1, l1, f2, l2;
912 GetNodePointsIndexes(inode, f1, l1, f2, l2);
913 Int_t sum = l1-f1+1;
914 sum += l2-f2+1;
915 return sum;
916}
917
918
919////////////////////////////////////////////////////////////////////////////////
920/// Set the data array. See the constructor function comments for details
921
922template <typename Index, typename Value>
924{
925// TO DO
926//
927// Check reconstruction/reallocation of memory of data. Maybe it is not
928// necessary to delete and realocate space but only to use the same space
929
930 Clear();
931
932 //Columnwise!!!!
933 fData = data;
934 fNPoints = npoints;
935 fNDim = ndim;
936 fBucketSize = bsize;
937
938 Build();
939}
940
941////////////////////////////////////////////////////////////////////////////////
942///Set the coordinate `#ndim` of all points (the column `#ndim` of the data matrix)
943///After setting all the data columns, proceed by calling Build() function
944///Note, that calling this function after Build() is not possible
945///Note also, that no checks on the array sizes is performed anywhere
946
947template <typename Index, typename Value>
949{
950 if (fAxis || fValue) {
951 Error("SetData", "The tree has already been built, no updates possible");
952 return 0;
953 }
954
955 if (!fData) {
956 fData = new Value*[fNDim];
957 }
958 fData[idim]=data;
959 fDataOwner = 2;
960 return 1;
961}
962
963
964////////////////////////////////////////////////////////////////////////////////
965///Calculate spread of the array a
966
967template <typename Index, typename Value>
969{
970 Index i;
971 min = a[index[0]];
972 max = a[index[0]];
973 for (i=0; i<ntotal; i++){
974 if (a[index[i]]<min) min = a[index[i]];
975 if (a[index[i]]>max) max = a[index[i]];
976 }
977}
978
979
980////////////////////////////////////////////////////////////////////////////////
981///
982///copy of the TMath::KOrdStat because I need an Index work array
983
984template <typename Index, typename Value>
986{
987 Index i, ir, j, l, mid;
988 Index arr;
989 Index temp;
990
991 Index rk = k;
992 l=0;
993 ir = ntotal-1;
994 for(;;) {
995 if (ir<=l+1) { //active partition contains 1 or 2 elements
996 if (ir == l+1 && a[index[ir]]<a[index[l]])
997 {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
998 Value tmp = a[index[rk]];
999 return tmp;
1000 } else {
1001 mid = (l+ir) >> 1; //choose median of left, center and right
1002 {temp = index[mid]; index[mid]=index[l+1]; index[l+1]=temp;}//elements as partitioning element arr.
1003 if (a[index[l]]>a[index[ir]]) //also rearrange so that a[l]<=a[l+1]
1004 {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
1005
1006 if (a[index[l+1]]>a[index[ir]])
1007 {temp=index[l+1]; index[l+1]=index[ir]; index[ir]=temp;}
1008
1009 if (a[index[l]]>a[index[l+1]])
1010 {temp = index[l]; index[l]=index[l+1]; index[l+1]=temp;}
1011
1012 i=l+1; //initialize pointers for partitioning
1013 j=ir;
1014 arr = index[l+1];
1015 for (;;) {
1016 do i++; while (a[index[i]]<a[arr]);
1017 do j--; while (a[index[j]]>a[arr]);
1018 if (j<i) break; //pointers crossed, partitioning complete
1019 {temp=index[i]; index[i]=index[j]; index[j]=temp;}
1020 }
1021 index[l+1]=index[j];
1022 index[j]=arr;
1023 if (j>=rk) ir = j-1; //keep active the partition that
1024 if (j<=rk) l=i; //contains the k_th element
1025 }
1026 }
1027}
1028
1029////////////////////////////////////////////////////////////////////////////////
1030/// Build boundaries for each node. Note, that the boundaries here are built
1031/// based on the splitting planes of the kd-tree, and don't necessarily pass
1032/// through the points of the original dataset. For the latter functionality
1033/// see function MakeBoundariesExact()
1034/// Boundaries can be retrieved by calling GetBoundary(inode) function that would
1035/// return an array of boundaries for the specified node, or GetBoundaries() function
1036/// that would return the complete array.
1037
1038template <typename Index, typename Value>
1040{
1041
1042 if(range) memcpy(fRange, range, fNDimm*sizeof(Value));
1043 // total number of nodes including terminal nodes
1044 Int_t totNodes = fNNodes + fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);
1045 fBoundaries = new Value[totNodes*fNDimm];
1046 //Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes));
1047
1048
1049 // loop
1050 Value *tbounds = 0x0, *cbounds = 0x0;
1051 Int_t cn;
1052 for(int inode=fNNodes-1; inode>=0; inode--){
1053 tbounds = &fBoundaries[inode*fNDimm];
1054 memcpy(tbounds, fRange, fNDimm*sizeof(Value));
1055
1056 // check left child node
1057 cn = (inode<<1)+1;
1058 if(IsTerminal(cn)) CookBoundaries(inode, kTRUE);
1059 cbounds = &fBoundaries[fNDimm*cn];
1060 for(int idim=0; idim<fNDim; idim++) tbounds[idim<<1] = cbounds[idim<<1];
1061
1062 // check right child node
1063 cn = (inode+1)<<1;
1064 if(IsTerminal(cn)) CookBoundaries(inode, kFALSE);
1065 cbounds = &fBoundaries[fNDimm*cn];
1066 for(int idim=0; idim<fNDim; idim++) tbounds[(idim<<1)+1] = cbounds[(idim<<1)+1];
1067 }
1068}
1069
1070////////////////////////////////////////////////////////////////////////////////
1071/// define index of this terminal node
1072
1073template <typename Index, typename Value>
1075{
1076 Int_t index = (node<<1) + (LEFT ? 1 : 2);
1077 //Info("CookBoundaries()", Form("Node %d", index));
1078
1079 // build and initialize boundaries for this node
1080 Value *tbounds = &fBoundaries[fNDimm*index];
1081 memcpy(tbounds, fRange, fNDimm*sizeof(Value));
1082 Bool_t flag[256]; // cope with up to 128 dimensions
1083 memset(flag, kFALSE, fNDimm);
1084 Int_t nvals = 0;
1085
1086 // recurse parent nodes
1087 Int_t pn = node;
1088 while(pn >= 0 && nvals < fNDimm){
1089 if(LEFT){
1090 index = (fAxis[pn]<<1)+1;
1091 if(!flag[index]) {
1092 tbounds[index] = fValue[pn];
1093 flag[index] = kTRUE;
1094 nvals++;
1095 }
1096 } else {
1097 index = fAxis[pn]<<1;
1098 if(!flag[index]) {
1099 tbounds[index] = fValue[pn];
1100 flag[index] = kTRUE;
1101 nvals++;
1102 }
1103 }
1104 LEFT = pn&1;
1105 pn = (pn - 1)>>1;
1106 }
1107}
1108
1109////////////////////////////////////////////////////////////////////////////////
1110/// Build boundaries for each node. Unlike MakeBoundaries() function
1111/// the boundaries built here always pass through a point of the original dataset
1112/// So, for example, for a terminal node with just one point minimum and maximum for each
1113/// dimension are the same.
1114/// Boundaries can be retrieved by calling GetBoundaryExact(inode) function that would
1115/// return an array of boundaries for the specified node, or GetBoundaries() function
1116/// that would return the complete array.
1117
1118template <typename Index, typename Value>
1120{
1121
1122 // total number of nodes including terminal nodes
1123 //Int_t totNodes = fNNodes + fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);
1124 if (fBoundaries){
1125 //boundaries were already computed for this tree
1126 return;
1127 }
1128 fBoundaries = new Value[fTotalNodes*fNDimm];
1129 Value *min = new Value[fNDim];
1130 Value *max = new Value[fNDim];
1131 for (Index inode=fNNodes; inode<fTotalNodes; inode++){
1132 //go through the terminal nodes
1133 for (Index idim=0; idim<fNDim; idim++){
1134 min[idim]= std::numeric_limits<Value>::max();
1135 max[idim]=-std::numeric_limits<Value>::max();
1136 }
1137 Index *points = GetPointsIndexes(inode);
1138 Index npoints = GetNPointsNode(inode);
1139 //find max and min in each dimension
1140 for (Index ipoint=0; ipoint<npoints; ipoint++){
1141 for (Index idim=0; idim<fNDim; idim++){
1142 if (fData[idim][points[ipoint]]<min[idim])
1143 min[idim]=fData[idim][points[ipoint]];
1144 if (fData[idim][points[ipoint]]>max[idim])
1145 max[idim]=fData[idim][points[ipoint]];
1146 }
1147 }
1148 for (Index idim=0; idim<fNDimm; idim+=2){
1149 fBoundaries[inode*fNDimm + idim]=min[idim/2];
1150 fBoundaries[inode*fNDimm + idim+1]=max[idim/2];
1151 }
1152 }
1153
1154 delete [] min;
1155 delete [] max;
1156
1157 Index left, right;
1158 for (Index inode=fNNodes-1; inode>=0; inode--){
1159 //take the min and max of left and right
1160 left = GetLeft(inode)*fNDimm;
1161 right = GetRight(inode)*fNDimm;
1162 for (Index idim=0; idim<fNDimm; idim+=2){
1163 //take the minimum on each dimension
1164 fBoundaries[inode*fNDimm+idim]=TMath::Min(fBoundaries[left+idim], fBoundaries[right+idim]);
1165 //take the maximum on each dimension
1166 fBoundaries[inode*fNDimm+idim+1]=TMath::Max(fBoundaries[left+idim+1], fBoundaries[right+idim+1]);
1167
1168 }
1169 }
1170}
1171
1172////////////////////////////////////////////////////////////////////////////////
1173///
1174/// find the smallest node covering the full range - start
1175///
1176
1177template <typename Index, typename Value>
1179 inode =0;
1180 for (;inode<fNNodes;){
1181 if (TMath::Abs(point[fAxis[inode]] - fValue[inode])<delta[fAxis[inode]]) break;
1182 inode = (point[fAxis[inode]] < fValue[inode]) ? (inode*2)+1: (inode*2)+2;
1183 }
1184}
1185
1186////////////////////////////////////////////////////////////////////////////////
1187/// Get the boundaries
1188
1189template <typename Index, typename Value>
1191{
1192 if(!fBoundaries) MakeBoundaries();
1193 return fBoundaries;
1194}
1195
1196
1197////////////////////////////////////////////////////////////////////////////////
1198/// Get the boundaries
1199
1200template <typename Index, typename Value>
1202{
1203 if(!fBoundaries) MakeBoundariesExact();
1204 return fBoundaries;
1205}
1206
1207////////////////////////////////////////////////////////////////////////////////
1208/// Get a boundary
1209
1210template <typename Index, typename Value>
1212{
1213 if(!fBoundaries) MakeBoundaries();
1214 return &fBoundaries[node*2*fNDim];
1215}
1216
1217////////////////////////////////////////////////////////////////////////////////
1218/// Get a boundary
1219
1220template <typename Index, typename Value>
1222{
1223 if(!fBoundaries) MakeBoundariesExact();
1224 return &fBoundaries[node*2*fNDim];
1225}
1226
1227
1228////////////////////////////////////////////////////////////////////////////////
1229///
1230/// Example function to
1231///
1232
1233TKDTreeIF *TKDTreeTestBuild(const Int_t npoints, const Int_t bsize){
1234 Float_t *data0 = new Float_t[npoints*2];
1235 Float_t *data[2];
1236 data[0] = &data0[0];
1237 data[1] = &data0[npoints];
1238 for (Int_t i=0;i<npoints;i++) {
1239 data[1][i]= gRandom->Rndm();
1240 data[0][i]= gRandom->Rndm();
1241 }
1242 TKDTree<Int_t, Float_t> *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
1243 return kdtree;
1244}
1245
1246
1247
1248template class TKDTree<Int_t, Float_t>;
1249template class TKDTree<Int_t, Double_t>;
#define d(i)
Definition: RSha256.hxx:102
int Int_t
Definition: RtypesCore.h:45
unsigned char UChar_t
Definition: RtypesCore.h:38
const Bool_t kFALSE
Definition: RtypesCore.h:101
float Float_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define templateClassImp(name)
Definition: Rtypes.h:420
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition: TError.cxx:220
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition: TError.cxx:231
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(const Int_t npoints, const Int_t bsize)
Example function to.
Definition: TKDTree.cxx:1233
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:2447
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:708
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:923
~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:846
Value * GetBoundaryExact(const Int_t node)
Get a boundary.
Definition: TKDTree.cxx:1221
void MakeBoundaries(Value *range=0x0)
Build boundaries for each node.
Definition: TKDTree.cxx:1039
void FindBNodeA(Value *point, Value *delta, Int_t &inode)
find the smallest node covering the full range - start
Definition: TKDTree.cxx:1178
void UpdateRange(Index inode, Value *point, Value range, std::vector< Index > &res)
Internal recursive function with the implementation of range searches.
Definition: TKDTree.cxx:764
Index FindNode(const Value *point) const
returns the index of the terminal node to which point belongs (index in the fAxis,...
Definition: TKDTree.cxx:677
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:985
Value * GetBoundariesExact()
Get the boundaries.
Definition: TKDTree.cxx:1201
void MakeBoundariesExact()
Build boundaries for each node.
Definition: TKDTree.cxx:1119
void CookBoundaries(const Int_t node, Bool_t left)
define index of this terminal node
Definition: TKDTree.cxx:1074
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:568
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:615
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:817
void Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const
Calculate spread of the array a.
Definition: TKDTree.cxx:968
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:640
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:754
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:900
Value * GetBoundaries()
Get the boundaries.
Definition: TKDTree.cxx:1190
Value * GetBoundary(const Int_t node)
Get a boundary.
Definition: TKDTree.cxx:1211
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:548
kNN::Event describes point in input variable vector-space, with additional functionality like distanc...
Mother of all ROOT objects.
Definition: TObject.h:37
Double_t Rndm() override
Machine independent random number generator.
Definition: TRandom.cxx:552
RooCmdArg Index(RooCategory &icat)
TF1 * f1
Definition: legend1.C:11
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
AxisAngle::Scalar Distance(const AxisAngle &r1, const R &r2)
Distance between two rotations.
Definition: AxisAngle.h:321
Element KOrdStat(Size n, const Element *a, Size k, Size *work=0)
Returns k_th order statistic of the array a of size n (k_th smallest element out of n elements).
Definition: TMath.h:1356
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:659
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
TArc a
Definition: textangle.C:12
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345