269 template <
typename Index,
typename Value>
291 template <
typename Index,
typename Value>
343 template <
typename Index,
typename Value>
372 template <
typename Index,
typename Value>
383 for(
int idim=0; idim<
fNDim; idim++)
delete []
fData[idim];
410 template <
typename Index,
typename Value>
454 Int_t nodeStack[128];
455 Int_t npointStack[128];
457 Int_t currentIndex = 0;
464 Int_t nbucketsall =0;
465 while (currentIndex>=0){
468 Int_t npoints = npointStack[currentIndex];
469 if (npoints<=fBucketSize) {
475 Int_t crow = rowStack[currentIndex];
476 Int_t cpos = posStack[currentIndex];
477 Int_t cnode = nodeStack[currentIndex];
482 if (npoints%fBucketSize) nbuckets0++;
483 Int_t restRows =
fRowT0-rowStack[currentIndex];
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;
490 if (nrest>(nfull/2)){
492 nright = npoints-nleft;
494 nright = nfull*fBucketSize/2;
495 nleft = npoints-nright;
501 Value tempspread, min, max;
504 for (Int_t idim=0; idim<
fNDim; idim++){
507 tempspread = max - min;
508 if (maxspread < tempspread) {
509 maxspread=tempspread;
516 array =
fData[axspread];
518 fAxis[cnode] = axspread;
523 npointStack[currentIndex] = nleft;
524 rowStack[currentIndex] = crow+1;
525 posStack[currentIndex] = cpos;
526 nodeStack[currentIndex] = cnode*2+1;
528 npointStack[currentIndex] = nright;
529 rowStack[currentIndex] = crow+1;
530 posStack[currentIndex] = cpos+nleft;
531 nodeStack[currentIndex] = (cnode*2)+2;
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");
547 template <
typename Index,
typename Value>
552 Error(
"FindNearestNeighbors",
"Working arrays must be allocated by the user!");
555 for (Int_t i=0; i<kNN; i++){
556 dist[i]=std::numeric_limits<Value>::max();
567 template <
typename Index,
typename Value>
573 if (min > dist[kNN-1]){
579 Index f1, l1, f2, l2;
581 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
586 while(ishift<kNN && d>dist[ishift])
590 for (Int_t i=kNN-1; i>ishift; i--){
614 template <
typename Index,
typename Value>
619 for (Int_t idim=0; idim<
fNDim; idim++){
620 dist+=(point[idim]-
fData[idim][ind])*(point[idim]-
fData[idim][ind]);
624 for (Int_t idim=0; idim<
fNDim; idim++){
625 dist+=TMath::Abs(point[idim]-
fData[idim][ind]);
639 template <
typename Index,
typename Value>
645 Double_t dist1, dist2;
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]);
652 if (point[idim/2]<bound[idim] || point[idim/2]>bound[idim+1])
653 min+= (dist1>dist2)? dist2 : dist1;
655 max+= (dist1>dist2)? dist1 : dist2;
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]);
664 min+= (dist1>dist2)? dist2 : dist1;
666 max+= (dist1>dist2)? dist1 : dist2;
676 template <
typename Index,
typename Value>
679 Index stackNode[128], inode;
680 Int_t currentIndex =0;
682 while (currentIndex>=0){
683 inode = stackNode[currentIndex];
689 stackNode[currentIndex]=(inode<<1)+1;
693 stackNode[currentIndex]=(inode+1)<<1;
707 template <
typename Index,
typename Value>
709 Int_t stackNode[128];
710 Int_t currentIndex =0;
714 while (currentIndex>=0){
716 Int_t inode = stackNode[currentIndex];
721 printf(
"terminal %d indexP %d\n", inode, indexIP);
722 for (Int_t ibucket=0;ibucket<
fBucketSize;ibucket++){
725 printf(
"ibucket %d index %d\n", ibucket, indexIP);
728 for (Int_t idim=0;idim<
fNDim;idim++)
if (
fData[idim][index0]!=point[idim]) isOK = kFALSE;
729 if (isOK) index = index0;
736 stackNode[currentIndex]=(inode*2)+1;
740 stackNode[currentIndex]=(inode*2)+2;
753 template <
typename Index,
typename Value>
763 template <
typename Index,
typename Value>
772 if (max<range && max>0) {
775 Index f1, l1, f2, l2;
778 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
781 for (Int_t ipoint=f2; ipoint<=l2; ipoint++){
790 Index f1, l1, f2, l2;
793 for (Int_t ipoint=f1; ipoint<=l1; ipoint++){
816 template <
typename Index,
typename Value>
820 printf(
"GetPointsIndexes() only for terminal nodes, use GetNodePointsIndexes() instead\n");
845 template <
typename Index,
typename Value>
862 Index f1, l1, f2, l2;
864 while (ileft<firsttermnode)
867 while (iright<firsttermnode)
899 template <
typename Index,
typename Value>
911 Int_t f1, l1, f2, l2;
922 template <
typename Index,
typename Value>
947 template <
typename Index,
typename Value>
951 Error(
"SetData",
"The tree has already been built, no updates possible");
967 template <
typename Index,
typename Value>
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]];
984 template <
typename Index,
typename Value>
987 Index i, ir, j, l,
mid;
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]];
1002 {temp = index[
mid]; index[
mid]=index[l+1]; index[l+1]=temp;}
1003 if (a[index[l]]>a[index[ir]])
1004 {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
1006 if (a[index[l+1]]>a[index[ir]])
1007 {temp=index[l+1]; index[l+1]=index[ir]; index[ir]=temp;}
1009 if (a[index[l]]>a[index[l+1]])
1010 {temp = index[l]; index[l]=index[l+1]; index[l+1]=temp;}
1016 do i++;
while (a[index[i]]<a[arr]);
1017 do j--;
while (a[index[j]]>a[arr]);
1019 {temp=index[i]; index[i]=index[j]; index[j]=temp;}
1021 index[l+1]=index[j];
1023 if (j>=rk) ir = j-1;
1038 template <
typename Index,
typename Value>
1050 Value *tbounds = 0x0, *cbounds = 0x0;
1052 for(
int inode=
fNNodes-1; inode>=0; inode--){
1060 for(
int idim=0; idim<
fNDim; idim++) tbounds[idim<<1] = cbounds[idim<<1];
1066 for(
int idim=0; idim<
fNDim; idim++) tbounds[(idim<<1)+1] = cbounds[(idim<<1)+1];
1073 template <
typename Index,
typename Value>
1076 Int_t index = (node<<1) + (LEFT ? 1 : 2);
1083 memset(flag, kFALSE,
fNDimm);
1088 while(pn >= 0 && nvals <
fNDimm){
1090 index = (
fAxis[pn]<<1)+1;
1092 tbounds[index] =
fValue[pn];
1093 flag[index] = kTRUE;
1097 index =
fAxis[pn]<<1;
1099 tbounds[index] =
fValue[pn];
1100 flag[index] = kTRUE;
1118 template <
typename Index,
typename Value>
1129 Value *min =
new Value[
fNDim];
1130 Value *max =
new Value[
fNDim];
1133 for (Index idim=0; idim<
fNDim; idim++){
1134 min[idim]= std::numeric_limits<Value>::max();
1135 max[idim]=-std::numeric_limits<Value>::max();
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]];
1148 for (Index idim=0; idim<
fNDimm; idim+=2){
1158 for (Index inode=
fNNodes-1; inode>=0; inode--){
1162 for (Index idim=0; idim<
fNDimm; idim+=2){
1177 template <
typename Index,
typename Value>
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;
1189 template <
typename Index,
typename Value>
1200 template <
typename Index,
typename Value>
1210 template <
typename Index,
typename Value>
1220 template <
typename Index,
typename Value>
1234 Float_t *data0 =
new Float_t[npoints*2];
1236 data[0] = &data0[0];
1237 data[1] = &data0[npoints];
1238 for (Int_t i=0;i<npoints;i++) {
double dist(Rotation3D const &r1, Rotation3D const &r2)
void Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const
Calculate spread of the array a.
Int_t fOffset
cross node - node that begins the last row (with terminal nodes only)
void SetData(Index npoints, Index ndim, UInt_t bsize, Value **data)
Set the data array. See the constructor function comments for details.
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
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t mid
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
TKDTree< Int_t, Float_t > TKDTreeIF
Bool_t IsTerminal(Index inode) const
void FindBNodeA(Value *point, Value *delta, Int_t &inode)
find the smallest node covering the full range - start
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...
void MakeBoundariesExact()
Build boundaries for each node.
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...
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 ...
void Build()
Build the kd-tree.
Int_t GetRight(Int_t inode) const
void FindPoint(Value *point, Index &index, Int_t &iter)
find the index of point works only if we keep fData pointers
TKDTree()
Default constructor. Nothing is built.
TKDTreeIF * TKDTreeTestBuild(const Int_t npoints, const Int_t bsize)
Example function to.
Index * fIndPoints
nodes boundaries
Class implementing a kd-tree.
virtual Double_t Rndm()
Machine independent random number generator.
void UpdateNearestNeighbors(Index inode, const Value *point, Int_t kNN, Index *ind, Value *dist)
Update the nearest neighbors values by examining the node inode.
Int_t fRowT0
array of points indexes
~TKDTree()
Destructor By default, the original data is not owned by kd-tree and is not deleted with it...
Value KOrdStat(Index ntotal, Value *a, Index k, Index *index) const
copy of the TMath::KOrdStat because I need an Index work array
Int_t GetLeft(Int_t inode) const
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...
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.
void CookBoundaries(const Int_t node, Bool_t left)
define index of this terminal node
R__EXTERN TRandom * gRandom
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.
void MakeBoundaries(Value *range=0x0)
Build boundaries for each node.
Value * GetBoundariesExact()
Get the boundaries.
Int_t fCrossNode
smallest terminal row - first row that contains terminal nodes
void UpdateRange(Index inode, Value *point, Value range, std::vector< Index > &res)
Internal recursive function with the implementation of range searches.
Int_t fNNodes
0 - not owner, 2 - owner of the pointer array, 1 - owner of the whole 2-d array
Value * GetBoundaries()
Get the boundaries.
Double_t Sqrt(Double_t x)
Value * GetBoundaryExact(const Int_t node)
Get a boundary.
Value * fBoundaries
data points
Value * GetBoundary(const Int_t node)
Get a boundary.