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