Logo ROOT   master
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 
26 Contents:
27 1. What is kd-tree
28 2. How to cosntruct kdtree - Pseudo code
29 3. Using TKDTree
30  a. Creating the kd-tree and setting the data
31  b. Navigating the kd-tree
32 4. 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 
41 In computer science, a kd-tree (short for k-dimensional tree) is a space-partitioning data structure
42 for organizing points in a k-dimensional space. kd-trees are a useful data structure for several
43 applications, such as searches involving a multidimensional search key (e.g. range searches and
44 nearest neighbour searches). kd-trees are a special case of BSP trees.
45 
46 A kd-tree uses only splitting planes that are perpendicular to one of the coordinate system axes.
47 This differs from BSP trees, in which arbitrary splitting planes can be used.
48 In addition, in the typical definition every node of a kd-tree, from the root to the leaves, stores a point.
49 This 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
51 the 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 
55 Since there are many possible ways to choose axis-aligned splitting planes, there are many different ways
56 to 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 
65 This method leads to a balanced kd-tree, in which each leaf node is about the same distance from the root.
66 However, balanced trees are not necessarily optimal for all applications.
67 The following pseudo-code illustrates this canonical construction procedure (NOTE, that the procedure used
68 by the TKDTree class is a bit different, the following pseudo-code is given as a simple illustration of the
69 concept):
70 
71 ~~~~
72 function 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 
94 Our construction method is optimized to save memory, and differs a bit from the constraints above.
95 In particular, the division axis is chosen as the one with the biggest spread, and the point to create the
96 splitting plane is chosen so, that one of the two subtrees contains exactly 2^k terminal nodes and is a
97 perfectly balanced binary tree, and, while at the same time, trying to keep the number of terminal nodes
98 in 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 
166 TKDtree is optimized to minimize memory consumption.
167 Nodes of the TKDTree do not store pointers to the left and right children or to the parent node,
168 but instead there are several 1-d arrays of size fNNodes with information about the nodes.
169 The order of the nodes information in the arrays is described below. It's important to understand
170 it, if one's class needs to store some kind of additional information on the per node basis, for
171 example, the fit function parameters.
172 
173 - Drawback: Insertion to the TKDtree is not supported.
174 - Advantage: Random access is supported
175 
176 As noted above, the construction of the kd-tree involves choosing the axis and the point on
177 that axis to divide the remaining points approximately in half. The exact algorithm for choosing
178 the division point is described in the next section. The sequence of divisions is
179 recorded in the following arrays:
180 ~~~~
181 fAxix[fNNodes] - Division axis (0,1,2,3 ...)
182 fValue[fNNodes] - Division value
183 ~~~~
184 
185 Given the index of a node in those arrays, it's easy to find the indices, corresponding to
186 children nodes or the parent node:
187 Suppose, 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 
191 Suppose, that the child node is stored under the index inode. Then:
192 - Parent `index = inode/2`
193 
194 Number of division nodes and number of terminals :
195 `fNNodes = (fNPoints/fBucketSize)`
196 
197 The nodes are filled always from left side to the right side:
198 Let inode be the index of a node, and irow - the index of a row
199 The TKDTree looks the following way:
200 Ideal case:
201 ~~~~
202 Number of _terminal_ nodes = 2^N, N=3
203 
204  INode
205 irow 0 0 - 1 inode
206 irow 1 1 2 - 2 inodes
207 irow 2 3 4 5 6 - 4 inodes
208 irow 3 7 8 9 10 11 12 13 14 - 8 inodes
209 ~~~~
210 
211 Non ideal case:
212 ~~~~
213 Number of _terminal_ nodes = 2^N+k, N=3 k=1
214 
215  INode
216 irow 0 0 - 1 inode
217 irow 1 1 2 - 2 inodes
218 irow 2 3 4 5 6 - 3 inodes
219 irow 3 7 8 9 10 11 12 13 14 - 8 inodes
220 irow 4 15 16 - 2 inodes
221 ~~~~
222 
223 #### 4b. The division algorithm:
224 
225 As described above, the kd-tree is built by repeatingly dividing the given set of points into
226 2 smaller sets. The cut is made on the axis with the biggest spread, and the value on the axis,
227 on which the cut is performed, is chosen based on the following formula:
228 Suppose, we want to divide n nodes into 2 groups, left and right. Then the left and right
229 will have the following number of nodes:
230 
231 ~~~~
232 n=2^k+rest
233 
234 Left = 2^k-1 + ((rest>2^k-2) ? 2^k-2 : rest)
235 Right = 2^k-1 + ((rest>2^k-2) ? rest-2^k-2 : 0)
236 ~~~~
237 
238 For example, let `n_nodes=67`. Then, the closest `2^k=64, 2^k-1=32, 2^k-2=16`.
239 Left node gets `32+3=35` sub-nodes, and the right node gets 32 sub-nodes
240 
241 The division process continues until all the nodes contain not more than a predefined number
242 of points.
243 
244 #### 4c. The order of nodes in boundary-related arrays
245 
246 Some kd-tree based algorithms need to know the boundaries of each node. This information can
247 be 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 
255 The 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 
269 template <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 
291 template <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 
343 template <typename Index, typename Value>
344 TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data) :
345  TObject()
346  ,fDataOwner(0)
347  ,fNNodes(0)
348  ,fTotalNodes(0)
349  ,fNDim(ndim)
350  ,fNDimm(2*ndim)
351  ,fNPoints(npoints)
352  ,fBucketSize(bsize)
353  ,fAxis(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 
372 template <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 
410 template <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 
547 template <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 
567 template <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 
614 template <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 
639 template <typename Index, typename Value>
640 void TKDTree<Index, Value>::DistanceToNode(const Value *point, Index inode, Value &min, Value &max, Int_t type)
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 
676 template <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 
707 template <typename Index, typename Value>
708 void TKDTree<Index, Value>::FindPoint(Value * point, Index &index, Int_t &iter){
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 
753 template <typename Index, typename Value>
754 void 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 
763 template <typename Index, typename Value>
764 void 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 
816 template <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 
845 template <typename Index, typename Value>
846 void 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 
899 template <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 
922 template <typename Index, typename Value>
923 void TKDTree<Index, Value>::SetData(Index npoints, Index ndim, UInt_t bsize, Value **data)
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 
947 template <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 
967 template <typename Index, typename Value>
968 void TKDTree<Index, Value>::Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const
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 
984 template <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 
1038 template <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 
1073 template <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 
1118 template <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 
1177 template <typename Index, typename Value>
1178  void TKDTree<Index, Value>::FindBNodeA(Value *point, Value *delta, Int_t &inode){
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 
1189 template <typename Index, typename Value>
1191 {
1192  if(!fBoundaries) MakeBoundaries();
1193  return fBoundaries;
1194 }
1195 
1196 
1197 ////////////////////////////////////////////////////////////////////////////////
1198 /// Get the boundaries
1199 
1200 template <typename Index, typename Value>
1202 {
1203  if(!fBoundaries) MakeBoundariesExact();
1204  return fBoundaries;
1205 }
1206 
1207 ////////////////////////////////////////////////////////////////////////////////
1208 /// Get a boundary
1209 
1210 template <typename Index, typename Value>
1212 {
1213  if(!fBoundaries) MakeBoundaries();
1214  return &fBoundaries[node*2*fNDim];
1215 }
1216 
1217 ////////////////////////////////////////////////////////////////////////////////
1218 /// Get a boundary
1219 
1220 template <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 
1233 TKDTreeIF *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 
1248 template class TKDTree<Int_t, Float_t>;
1249 template class TKDTree<Int_t, Double_t>;
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
void Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const
Calculate spread of the array a.
Definition: TKDTree.cxx:968
static long int sum(long int i)
Definition: Factory.cxx:2274
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
Index FindNode(const Value *point) const
returns the index of the terminal node to which point belongs (index in the fAxis, fValue, etc arrays) returns -1 in case of failure
Definition: TKDTree.cxx:677
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t mid
Definition: TRolke.cxx:630
Index * GetPointsIndexes(Int_t node) const
return the indices of the points in that terminal node for all the nodes except last, the size is fBucketSize for the last node it&#39;s fOffsetfBucketSize
Definition: TKDTree.cxx:817
TKDTree< Int_t, Float_t > TKDTreeIF
Definition: TKDTree.h:105
float Float_t
Definition: RtypesCore.h:53
void FindBNodeA(Value *point, Value *delta, Int_t &inode)
find the smallest node covering the full range - start
Definition: TKDTree.cxx:1178
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
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
int Int_t
Definition: RtypesCore.h:41
void MakeBoundariesExact()
Build boundaries for each node.
Definition: TKDTree.cxx:1119
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
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
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
void Build()
Build the kd-tree.
Definition: TKDTree.cxx:411
#define templateClassImp(name)
Definition: Rtypes.h:409
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
TKDTree()
Default constructor. Nothing is built.
Definition: TKDTree.cxx:270
void Info(const char *location, const char *msgfmt,...)
TKDTreeIF * TKDTreeTestBuild(const Int_t npoints, const Int_t bsize)
Example function to.
Definition: TKDTree.cxx:1233
void Error(const char *location, const char *msgfmt,...)
Class implementing a kd-tree
Definition: TKDTree.h:9
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
point * points
Definition: X3DBuffer.c:22
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
~TKDTree()
Destructor By default, the original data is not owned by kd-tree and is not deleted with it...
Definition: TKDTree.cxx:373
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
auto * a
Definition: textangle.C:12
PyObject * fValue
unsigned int UInt_t
Definition: RtypesCore.h:42
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
char * Form(const char *fmt,...)
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, 0 on failure Arrays ind and dist are provided by the user and are assumed to be at least kNN elements long.
Definition: TKDTree.cxx:548
void CookBoundaries(const Int_t node, Bool_t left)
define index of this terminal node
Definition: TKDTree.cxx:1074
void Warning(const char *location, const char *msgfmt,...)
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
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
const Bool_t kFALSE
Definition: RtypesCore.h:88
void MakeBoundaries(Value *range=0x0)
Build boundaries for each node.
Definition: TKDTree.cxx:1039
#define d(i)
Definition: RSha256.hxx:102
Value * GetBoundariesExact()
Get the boundaries.
Definition: TKDTree.cxx:1201
int type
Definition: TGX11.cxx:120
Mother of all ROOT objects.
Definition: TObject.h:37
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:1323
auto * l
Definition: textangle.C:4
RooCmdArg Index(RooCategory &icat)
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
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
TF1 * f1
Definition: legend1.C:11
AxisAngle::Scalar Distance(const AxisAngle &r1, const R &r2)
Distance between two rotations.
Definition: AxisAngle.h:326
Value * GetBoundaries()
Get the boundaries.
Definition: TKDTree.cxx:1190
unsigned char UChar_t
Definition: RtypesCore.h:34
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Value * GetBoundaryExact(const Int_t node)
Get a boundary.
Definition: TKDTree.cxx:1221
const Bool_t kTRUE
Definition: RtypesCore.h:87
Value * GetBoundary(const Int_t node)
Get a boundary.
Definition: TKDTree.cxx:1211
const char * Value
Definition: TXMLSetup.cxx:72