ROOT   6.10/09 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 /*************************************************************************
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>
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.
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;
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
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;
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  }
517  KOrdStat(npoints, array, nleft, fIndPoints+cpos);
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  }
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 {
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>
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) {
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  }
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>
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 {
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 {
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>;
virtual void Clear(Option_t *="")
Definition: TObject.h:91
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:2162
Int_t fOffset
cross node - node that begins the last row (with terminal nodes only)
Definition: TKDTree.h:95
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:847
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
Int_t fDataOwner
Definition: TKDTree.h:77
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
Bool_t IsTerminal(Index inode) const
Definition: TKDTree.h:56
Int_t fTotalNodes
Definition: TKDTree.h:79
void FindBNodeA(Value *point, Value *delta, Int_t &inode)
find the smallest node covering the full range - start
Definition: TKDTree.cxx:1178
Index fNPoints
Definition: TKDTree.h:82
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:168
int Int_t
Definition: RtypesCore.h:41
Index fBucketSize
Definition: TKDTree.h:83
bool Bool_t
Definition: RtypesCore.h:59
void MakeBoundariesExact()
Build boundaries for each node.
Definition: TKDTree.cxx:1119
TArc * a
Definition: textangle.C:12
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
Value * fRange
Definition: TKDTree.h:87
Float_t delta
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
void Build()
Build the kd-tree.
Definition: TKDTree.cxx:411
Index fNDim
Definition: TKDTree.h:80
#define templateClassImp(name)
Definition: Rtypes.h:382
Int_t GetRight(Int_t inode) const
Definition: TKDTree.h:25
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
Int_t bsize[]
Definition: SparseFit4.cxx:31
TKDTree()
Default constructor. Nothing is built.
Definition: TKDTree.cxx:270
TKDTreeIF * TKDTreeTestBuild(const Int_t npoints, const Int_t bsize)
Example function to.
Definition: TKDTree.cxx:1233
Index * fIndPoints
nodes boundaries
Definition: TKDTree.h:92
Class implementing a kd-tree.
Definition: TKDTree.h:9
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:512
point * points
Definition: X3DBuffer.c:20
UChar_t * fAxis
Definition: TKDTree.h:84
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
Int_t fRowT0
array of points indexes
Definition: TKDTree.h:93
~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
PyObject * fValue
Int_t GetLeft(Int_t inode) const
Definition: TKDTree.h:24
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
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:873
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
TLine * l
Definition: textangle.C:4
void CookBoundaries(const Int_t node, Bool_t left)
define index of this terminal node
Definition: TKDTree.cxx:1074
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:92
Value * fValue
Definition: TKDTree.h:85
void MakeBoundaries(Value *range=0x0)
Build boundaries for each node.
Definition: TKDTree.cxx:1039
RooCmdArg Index(RooCategory &icat)
Value * GetBoundariesExact()
Get the boundaries.
Definition: TKDTree.cxx:1201
double Double_t
Definition: RtypesCore.h:55
Int_t fCrossNode
smallest terminal row - first row that contains terminal nodes
Definition: TKDTree.h:94
int type
Definition: TGX11.cxx:120
Index fNDimm
Definition: TKDTree.h:81
Value ** fData
Definition: TKDTree.h:88
Mother of all ROOT objects.
Definition: TObject.h:37
double f2(const double *x)
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
Int_t fNNodes
0 - not owner, 2 - owner of the pointer array, 1 - owner of the whole 2-d array
Definition: TKDTree.h:78
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
TF1 * f1
Definition: legend1.C:11
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:591
Value * GetBoundaryExact(const Int_t node)
Get a boundary.
Definition: TKDTree.cxx:1221
Value * fBoundaries
data points
Definition: TKDTree.h:89
const Bool_t kTRUE
Definition: RtypesCore.h:91
Value * GetBoundary(const Int_t node)
Get a boundary.
Definition: TKDTree.cxx:1211
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:859
const char * Value
Definition: TXMLSetup.cxx:73