#include "TKDTree.h"
+#include "TRandom.h"
#include "TString.h"
#include <string.h>
#ifndef R__ALPHA
templateClassImp(TKDTree)
#endif
+
+/*
+
+
+
+Content:
+1. What is kd-tree
+2. How to cosntruct kdtree - Pseudo code
+3. TKDTree implementation
+4. User interface
+
+
+
+1. What is kdtree ? ( http://en.wikipedia.org/wiki/Kd-tree )
+
+In computer science, a kd-tree (short for k-dimensional tree) is a space-partitioning data structure for organizing points in a k-dimensional space. kd-trees are a useful data structure for several applications, such as searches involving a multidimensional search key (e.g. range searches and nearest neighbour searches). kd-trees are a special case of BSP trees.
+
+A kd-tree uses only splitting planes that are perpendicular to one of the coordinate system axes. This differs from BSP trees, in which arbitrary splitting planes can be used. In addition, in the typical definition every node of a kd-tree, from the root to the leaves, stores a point.[1] This differs from BSP trees, in which leaves are typically the only nodes that contain points (or other geometric primitives). As a consequence, each splitting plane must go through one of the points in the kd-tree. kd-tries are a variant that store data only in leaf nodes. It is worth noting that in an alternative definition of kd-tree the points are stored in its leaf nodes only, although each splitting plane still goes through one of the points.[2]
+<img src=".....">
+
+2. Constructing a kd-tree ( Pseudo code)
+
+Since there are many possible ways to choose axis-aligned splitting planes, there are many different ways to construct kd-trees. The canonical method of kd-tree construction has the following constraints:
+
+ * As one moves down the tree, one cycles through the axes used to select the splitting planes. (For example, the root would have an x-aligned plane, the root's children would both have y-aligned planes, the root's grandchildren would all have z-aligned planes, and so on.)
+ * At each step, the point selected to create the splitting plane is the median of the points being put into the kd-tree, with respect to their coordinates in the axis being used. (Note the assumption that we feed the entire set of points into the algorithm up-front.)
+
+This method leads to a balanced kd-tree, in which each leaf node is about the same distance from the root. However, balanced trees are not necessarily optimal for all applications.
+
+Note also that it is not required to select the median point. In that case, the result is simply that there is no guarantee that the tree will be balanced. A simple heuristic to avoid coding a complex linear-time median-finding algorithm nor using an O(n log n) sort is to use sort to find the median of a fixed number of randomly selected points to serve as the cut line. Practically this technique results in nicely balanced trees.
+
+Given a list of n points, the following algorithm will construct a balanced kd-tree containing those points.
+<img src=".....">
+
+function kdtree (list of points pointList, int depth)
+{
+ if pointList is empty
+ return nil;
+ else
+ {
+ // Select axis based on depth so that axis cycles through all valid values
+ var int axis := depth mod k;
+
+ // Sort point list and choose median as pivot element
+ select median from pointList;
+
+ // Create node and construct subtrees
+ var tree_node node;
+ node.location := median;
+ node.leftChild := kdtree(points in pointList before median, depth+1);
+ node.rightChild := kdtree(points in pointList after median, depth+1);
+ return node;
+ }
+}
+
+
+3. TKDtree implementation.
+
+// TKDtree is optimized to minimize memory consumption.
+// TKDTree is by default balanced. Nodes are mapped to the 1 dimensional arrays of size
+// fNnodes.
+//
+// fAxix[fNnodes] - Division axis (0,1,2,3 ...)
+// fValue[fNnodes] - Division value
+//
+// Navigation to the Left and right doughter node is expressed by analytical formula
+// Let's parent node id is inode
+// Then:
+// Left = inode*2+1
+// Right = (inode+1)*2
+// Let's daughter node id is inode
+// Then:
+// Parent = inode/2
+//
+//
+// The mapping of the nodes to the 1D array enables us to use the indexing mechanism.
+// This is important if some additional kind of information
+// (not directly connected to kd-tree, e.g function fit parameters in the node)
+// follow binary tree structure. The mapping can be reused.
+//
+// Drawback: Insertion to the TKDtree is not supported.
+// Advantage: Random access is supported - simple formulas
+//
+// Number of division nodes and number of terminals :
+// fNnodes = (fNPoints/fBucketSize)
+//
+// TKDTree is mapped to one dimensional array.
+// Mapping: (see for example the TKDTree::FindNode)
+//
+// The nodes are filled always from left side to the right side:
+//
+//
+// TERMINOLOGY:
+// - inode - index of node
+// - irow - index of row
+
+// Ideal case:
+// Number of terminal nodes = 2^N - N=3
+//
+// INode
+// irow 0 0 - 1 inode
+// irow 1 1 2 - 2 inodes
+// irow 2 3 4 5 6 - 4 inodes
+// irow 3 7 8 9 10 11 12 13 14 - 8 inodes
+//
+//
+// Non ideal case:
+// Number of terminal nodes = 2^N+k - N=3 k=1
+
+// INode
+// irow 0 0 - 1 inode
+// irow 1 1 2 - 2 inodes
+// irow 2 3 4 5 6 - 3 inodes
+// irow 3 7 8 9 10 11 12 13 14 - 8 inodes
+// irow 4 15 16 - 2 inodes
+//
+//
+// The division algorithm:
+// The n points are divided to tho groups.
+//
+// n=2^k+rest
+//
+// Left = 2^k-1 +(rest>2^k-2) ? 2^k-2 : rest
+// Right = 2^k-1 +(rest>2^k-2) ? rest-2^k-2 : 0
+//
+//
+
+
+
+
+*/
+
+
+
//////////////////////////////////////////////////////////////////////
// Memory setup of protected data members:
//
//_________________________________________________________________
template <typename Index, typename Value>
TKDTree<Index, Value>::TKDTree() :
- TObject()
- ,fDataOwner(kFALSE)
- ,fNnodes(0)
- ,fNDim(0)
- ,fNDimm(0)
- ,fNpoints(0)
- ,fBucketSize(0)
- ,fAxis(0x0)
- ,fValue(0x0)
- ,fRange(0x0)
- ,fData(0x0)
- ,fBoundaries(0x0)
- ,fkNNdim(0)
- ,fkNN(0x0)
- ,fkNNdist(0x0)
- ,fDistBuffer(0x0)
- ,fIndBuffer(0x0)
- ,fIndPoints(0x0)
- ,fRowT0(0)
- ,fCrossNode(0)
- ,fOffset(0)
+ TObject()
+ ,fDataOwner(kFALSE)
+ ,fNnodes(0)
+ ,fNDim(0)
+ ,fNDimm(0)
+ ,fNpoints(0)
+ ,fBucketSize(0)
+ ,fAxis(0x0)
+ ,fValue(0x0)
+ ,fRange(0x0)
+ ,fData(0x0)
+ ,fBoundaries(0x0)
+ ,fkNNdim(0)
+ ,fkNN(0x0)
+ ,fkNNdist(0x0)
+ ,fDistBuffer(0x0)
+ ,fIndBuffer(0x0)
+ ,fIndPoints(0x0)
+ ,fRowT0(0)
+ ,fCrossNode(0)
+ ,fOffset(0)
{
// Default constructor. Nothing is built
}
//_________________________________________________________________
template <typename Index, typename Value>
TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data) :
- TObject()
- ,fDataOwner(kFALSE)
- ,fNnodes(0)
- ,fNDim(ndim)
- ,fNDimm(2*ndim)
- ,fNpoints(npoints)
- ,fBucketSize(bsize)
- ,fAxis(0x0)
- ,fValue(0x0)
- ,fRange(0x0)
- ,fData(data) //Columnwise!!!!!
- ,fBoundaries(0x0)
- ,fkNNdim(0)
- ,fkNN(0x0)
- ,fkNNdist(0x0)
- ,fDistBuffer(0x0)
- ,fIndBuffer(0x0)
- ,fIndPoints(0x0)
- ,fRowT0(0)
- ,fCrossNode(0)
- ,fOffset(0)
+ TObject()
+ ,fDataOwner(kFALSE)
+ ,fNnodes(0)
+ ,fNDim(ndim)
+ ,fNDimm(2*ndim)
+ ,fNpoints(npoints)
+ ,fBucketSize(bsize)
+ ,fAxis(0x0)
+ ,fValue(0x0)
+ ,fRange(0x0)
+ ,fData(data) //Columnwise!!!!!
+ ,fBoundaries(0x0)
+ ,fkNNdim(0)
+ ,fkNN(0x0)
+ ,fkNNdist(0x0)
+ ,fDistBuffer(0x0)
+ ,fIndBuffer(0x0)
+ ,fIndPoints(0x0)
+ ,fRowT0(0)
+ ,fCrossNode(0)
+ ,fOffset(0)
{
-// Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
-
- Build();
+ // Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
+
+ Build();
}
//_________________________________________________________________
TKDTree<Index, Value>::~TKDTree()
{
// Destructor
- if (fIndBuffer) delete [] fIndBuffer;
- if (fDistBuffer) delete [] fDistBuffer;
- if (fkNNdist) delete [] fkNNdist;
- if (fkNN) delete [] fkNN;
- if (fAxis) delete [] fAxis;
- if (fValue) delete [] fValue;
- if (fIndPoints) delete [] fIndPoints;
- if (fRange) delete [] fRange;
- if (fBoundaries) delete [] fBoundaries;
- if (fDataOwner && fData){
- for(int idim=0; idim<fNDim; idim++) delete [] fData[idim];
- delete [] fData;
- }
+ if (fIndBuffer) delete [] fIndBuffer;
+ if (fDistBuffer) delete [] fDistBuffer;
+ if (fkNNdist) delete [] fkNNdist;
+ if (fkNN) delete [] fkNN;
+ if (fAxis) delete [] fAxis;
+ if (fValue) delete [] fValue;
+ if (fIndPoints) delete [] fIndPoints;
+ if (fRange) delete [] fRange;
+ if (fBoundaries) delete [] fBoundaries;
+ if (fDataOwner && fData){
+ for(int idim=0; idim<fNDim; idim++) delete [] fData[idim];
+ delete [] fData;
+ }
}
//_________________________________________________________________
template <typename Index, typename Value>
void TKDTree<Index, Value>::Build(){
- //
- // Build binning
- //
- // 1. calculate number of nodes
- // 2. calculate first terminal row
- // 3. initialize index array
- // 4. non recursive building of the binary tree
- //
- //1.
- fNnodes = fNpoints/fBucketSize-1;
- if (fNpoints%fBucketSize) fNnodes++;
-
- //2.
- fRowT0=0;
- for (;(fNnodes+1)>(1<<fRowT0);fRowT0++);
- fRowT0-=1;
- // 2 = 2**0 + 1
- // 3 = 2**1 + 1
- // 4 = 2**1 + 2
- // 5 = 2**2 + 1
- // 6 = 2**2 + 2
- // 7 = 2**2 + 3
- // 8 = 2**2 + 4
-
- //3.
- // allocate space for boundaries
- fRange = new Value[2*fNDim];
- fIndPoints= new Index[fNpoints];
- for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i;
- fAxis = new UChar_t[fNnodes];
- fValue = new Value[fNnodes];
- //
- fCrossNode = (1<<(fRowT0+1))-1;
- if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1;
- //
- // fOffset = (((fNnodes+1)-(1<<fRowT0)))*2;
- Int_t over = (fNnodes+1)-(1<<fRowT0);
- Int_t filled = ((1<<fRowT0)-over)*fBucketSize;
- fOffset = fNpoints-filled;
- //
-// printf("Row0 %d\n", fRowT0);
-// printf("CrossNode %d\n", fCrossNode);
-// printf("Offset %d\n", fOffset);
- //
- //
- //4.
- // stack for non recursive build - size 128 bytes enough
- Int_t rowStack[128];
- Int_t nodeStack[128];
- Int_t npointStack[128];
- Int_t posStack[128];
- Int_t currentIndex = 0;
- Int_t iter =0;
- rowStack[0] = 0;
- nodeStack[0] = 0;
- npointStack[0] = fNpoints;
- posStack[0] = 0;
- //
- Int_t nbucketsall =0;
- while (currentIndex>=0){
- iter++;
- //
- Int_t npoints = npointStack[currentIndex];
- if (npoints<=fBucketSize) {
- //printf("terminal node : index %d iter %d\n", currentIndex, iter);
- currentIndex--;
- nbucketsall++;
- continue; // terminal node
- }
- Int_t crow = rowStack[currentIndex];
- Int_t cpos = posStack[currentIndex];
- Int_t cnode = nodeStack[currentIndex];
- //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
- //
- // divide points
- Int_t nbuckets0 = npoints/fBucketSize; //current number of buckets
- if (npoints%fBucketSize) nbuckets0++; //
- Int_t restRows = fRowT0-rowStack[currentIndex]; // rest of fully occupied node row
- if (restRows<0) restRows =0;
- for (;nbuckets0>(2<<restRows); restRows++);
- Int_t nfull = 1<<restRows;
- Int_t nrest = nbuckets0-nfull;
- Int_t nleft =0, nright =0;
- //
- if (nrest>(nfull/2)){
- nleft = nfull*fBucketSize;
- nright = npoints-nleft;
- }else{
- nright = nfull*fBucketSize/2;
- nleft = npoints-nright;
- }
-
- //
- //find the axis with biggest spread
- Value maxspread=0;
- Value tempspread, min, max;
- Index axspread=0;
- Value *array;
- for (Int_t idim=0; idim<fNDim; idim++){
- array = fData[idim];
- Spread(npoints, array, fIndPoints+cpos, min, max);
- tempspread = max - min;
- if (maxspread < tempspread) {
- maxspread=tempspread;
- axspread = idim;
- }
- if(cnode) continue;
- //printf("set %d %6.3f %6.3f\n", idim, min, max);
- fRange[2*idim] = min; fRange[2*idim+1] = max;
- }
- array = fData[axspread];
- KOrdStat(npoints, array, nleft, fIndPoints+cpos);
- fAxis[cnode] = axspread;
- fValue[cnode] = array[fIndPoints[cpos+nleft]];
- //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
- //
- //
- npointStack[currentIndex] = nleft;
- rowStack[currentIndex] = crow+1;
- posStack[currentIndex] = cpos;
- nodeStack[currentIndex] = cnode*2+1;
- currentIndex++;
- npointStack[currentIndex] = nright;
- rowStack[currentIndex] = crow+1;
- posStack[currentIndex] = cpos+nleft;
- nodeStack[currentIndex] = (cnode*2)+2;
- //
- if (0){
- // consistency check
- Info("Build()", Form("points %d left %d right %d", npoints, nleft, nright));
- if (nleft<nright) Warning("Build", "Problem Left-Right");
- if (nleft<0 || nright<0) Warning("Build()", "Problem Negative number");
- }
- }
-
- //printf("NBuckets\t%d\n", nbucketsall);
- //fData = 0;
+ //
+ // Build binning
+ //
+ // 1. calculate number of nodes
+ // 2. calculate first terminal row
+ // 3. initialize index array
+ // 4. non recursive building of the binary tree
+
+ //
+ // Teh tree is recursivaly divided. The division point is choosen in such a way, that the left side
+ // alway has 2^k points, while at the same time trying
+ //
+ //1.
+ fNnodes = fNpoints/fBucketSize-1;
+ if (fNpoints%fBucketSize) fNnodes++;
+
+ //2.
+ fRowT0=0;
+ for (;(fNnodes+1)>(1<<fRowT0);fRowT0++);
+ fRowT0-=1;
+ // 2 = 2**0 + 1
+ // 3 = 2**1 + 1
+ // 4 = 2**1 + 2
+ // 5 = 2**2 + 1
+ // 6 = 2**2 + 2
+ // 7 = 2**2 + 3
+ // 8 = 2**2 + 4
+
+ //3.
+ // allocate space for boundaries
+ fRange = new Value[2*fNDim];
+ fIndPoints= new Index[fNpoints];
+ for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i;
+ fAxis = new UChar_t[fNnodes];
+ fValue = new Value[fNnodes];
+ //
+ fCrossNode = (1<<(fRowT0+1))-1;
+ if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1;
+ //
+ // fOffset = (((fNnodes+1)-(1<<fRowT0)))*2;
+ Int_t over = (fNnodes+1)-(1<<fRowT0);
+ Int_t filled = ((1<<fRowT0)-over)*fBucketSize;
+ fOffset = fNpoints-filled;
+ //
+ // printf("Row0 %d\n", fRowT0);
+ // printf("CrossNode %d\n", fCrossNode);
+ // printf("Offset %d\n", fOffset);
+ //
+ //
+ //4.
+ // stack for non recursive build - size 128 bytes enough
+ Int_t rowStack[128];
+ Int_t nodeStack[128];
+ Int_t npointStack[128];
+ Int_t posStack[128];
+ Int_t currentIndex = 0;
+ Int_t iter =0;
+ rowStack[0] = 0;
+ nodeStack[0] = 0;
+ npointStack[0] = fNpoints;
+ posStack[0] = 0;
+ //
+ Int_t nbucketsall =0;
+ while (currentIndex>=0){
+ iter++;
+ //
+ Int_t npoints = npointStack[currentIndex];
+ if (npoints<=fBucketSize) {
+ //printf("terminal node : index %d iter %d\n", currentIndex, iter);
+ currentIndex--;
+ nbucketsall++;
+ continue; // terminal node
+ }
+ Int_t crow = rowStack[currentIndex];
+ Int_t cpos = posStack[currentIndex];
+ Int_t cnode = nodeStack[currentIndex];
+ //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
+ //
+ // divide points
+ Int_t nbuckets0 = npoints/fBucketSize; //current number of buckets
+ if (npoints%fBucketSize) nbuckets0++; //
+ Int_t restRows = fRowT0-rowStack[currentIndex]; // rest of fully occupied node row
+ if (restRows<0) restRows =0;
+ for (;nbuckets0>(2<<restRows); restRows++);
+ Int_t nfull = 1<<restRows;
+ Int_t nrest = nbuckets0-nfull;
+ Int_t nleft =0, nright =0;
+ //
+ if (nrest>(nfull/2)){
+ nleft = nfull*fBucketSize;
+ nright = npoints-nleft;
+ }else{
+ nright = nfull*fBucketSize/2;
+ nleft = npoints-nright;
+ }
+
+ //
+ //find the axis with biggest spread
+ Value maxspread=0;
+ Value tempspread, min, max;
+ Index axspread=0;
+ Value *array;
+ for (Int_t idim=0; idim<fNDim; idim++){
+ array = fData[idim];
+ Spread(npoints, array, fIndPoints+cpos, min, max);
+ tempspread = max - min;
+ if (maxspread < tempspread) {
+ maxspread=tempspread;
+ axspread = idim;
+ }
+ if(cnode) continue;
+ //printf("set %d %6.3f %6.3f\n", idim, min, max);
+ fRange[2*idim] = min; fRange[2*idim+1] = max;
+ }
+ array = fData[axspread];
+ KOrdStat(npoints, array, nleft, fIndPoints+cpos);
+ fAxis[cnode] = axspread;
+ fValue[cnode] = array[fIndPoints[cpos+nleft]];
+ //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
+ //
+ //
+ npointStack[currentIndex] = nleft;
+ rowStack[currentIndex] = crow+1;
+ posStack[currentIndex] = cpos;
+ nodeStack[currentIndex] = cnode*2+1;
+ currentIndex++;
+ npointStack[currentIndex] = nright;
+ rowStack[currentIndex] = crow+1;
+ posStack[currentIndex] = cpos+nleft;
+ nodeStack[currentIndex] = (cnode*2)+2;
+ //
+ if (0){
+ // consistency check
+ Info("Build()", Form("points %d left %d right %d", npoints, nleft, nright));
+ if (nleft<nright) Warning("Build", "Problem Left-Right");
+ if (nleft<0 || nright<0) Warning("Build()", "Problem Negative number");
+ }
+ }
}
template <typename Index, typename Value>
Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value *&d)
{
-// Find "k" nearest neighbors to "point".
-//
-// Return true in case of success and false in case of failure.
-// The indexes of the nearest k points are stored in the array "in" in
-// increasing order of the distance to "point" and the maxim distance
-// in "d".
-//
-// The arrays "in" and "d" are managed internally by the TKDTree.
-
- Index inode = FindNode(point);
- if(inode < fNnodes){
- Error("FindNearestNeighbors()", "Point outside data range.");
- return kFALSE;
- }
-
- UInt_t debug = 0;
- Int_t nCalculations = 0;
-
- // allocate working memory
- if(!fDistBuffer){
- fDistBuffer = new Value[fBucketSize];
- fIndBuffer = new Index[fBucketSize];
- }
- if(fkNNdim < k){
- if(debug>=1){
- if(fkNN) printf("Reallocate memory %d -> %d \n", fkNNdim, 2*k);
- else printf("Allocate %d memory\n", 2*k);
- }
- fkNNdim = 2*k;
- if(fkNN){
- delete [] fkNN;
- delete [] fkNNdist;
- }
- fkNN = new Index[fkNNdim];
- fkNNdist = new Value[fkNNdim];
+ //
+ // NOT MI - to be checked
+ //
+ //
+ // Find "k" nearest neighbors to "point".
+ //
+ // Return true in case of success and false in case of failure.
+ // The indexes of the nearest k points are stored in the array "in" in
+ // increasing order of the distance to "point" and the maxim distance
+ // in "d".
+ //
+ // The arrays "in" and "d" are managed internally by the TKDTree.
+
+ Index inode = FindNode(point);
+ if(inode < fNnodes){
+ Error("FindNearestNeighbors()", "Point outside data range.");
+ return kFALSE;
+ }
+
+ UInt_t debug = 0;
+ Int_t nCalculations = 0;
+
+ // allocate working memory
+ if(!fDistBuffer){
+ fDistBuffer = new Value[fBucketSize];
+ fIndBuffer = new Index[fBucketSize];
+ }
+ if(fkNNdim < k){
+ if(debug>=1){
+ if(fkNN) printf("Reallocate memory %d -> %d \n", fkNNdim, 2*k);
+ else printf("Allocate %d memory\n", 2*k);
+ }
+ fkNNdim = 2*k;
+ if(fkNN){
+ delete [] fkNN;
+ delete [] fkNNdist;
+ }
+ fkNN = new Index[fkNNdim];
+ fkNNdist = new Value[fkNNdim];
+ }
+ memset(fkNN, -1, k*sizeof(Index));
+ for(int i=0; i<k; i++) fkNNdist[i] = 9999.;
+ Index itmp, jtmp; Value ftmp, gtmp;
+
+ // calculate number of boundaries with the data domain.
+ if(!fBoundaries) MakeBoundaries();
+ Value *bounds = &fBoundaries[inode*2*fNDim];
+ Index nBounds = 0;
+ for(int idim=0; idim<fNDim; idim++){
+ if((bounds[2*idim] - fRange[2*idim]) < 1.E-10) nBounds++;
+ if((bounds[2*idim+1] - fRange[2*idim+1]) < 1.E-10) nBounds++;
+ }
+ if(debug>=1) printf("Calculate boundaries [nBounds = %d]\n", nBounds);
+
+
+ // traverse tree
+ UChar_t ax;
+ Value val, dif;
+ Int_t nAllNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
+ Int_t nodeStack[128], nodeIn[128];
+ Index currentIndex = 0;
+ nodeStack[0] = inode;
+ nodeIn[0] = inode;
+ while(currentIndex>=0){
+ Int_t tnode = nodeStack[currentIndex];
+ Int_t entry = nodeIn[currentIndex];
+ currentIndex--;
+ if(debug>=1) printf("Analyse tnode %d entry %d\n", tnode, entry);
+
+ // check if node is still eligible
+ Int_t inode = (tnode-1)>>1; //tnode/2 + (tnode%2) - 1;
+ ax = fAxis[inode];
+ val = fValue[inode];
+ dif = point[ax] - val;
+ if((TMath::Abs(dif) > fkNNdist[k-1]) &&
+ ((dif > 0. && (tnode&1)) || (dif < 0. && !(tnode&1)))) {
+ if(debug>=1) printf("\tremove %d\n", (tnode-1)>>1);
+
+ // mark bound
+ nBounds++;
+ // end all recursions
+ if(nBounds==2 * fNDim) break;
+ continue;
+ }
+
+ if(IsTerminal(tnode)){
+ if(debug>=2) printf("\tProcess terminal node %d\n", tnode);
+ // Link data to terminal node
+ Int_t offset = (tnode >= fCrossNode) ? (tnode-fCrossNode)*fBucketSize : fOffset+(tnode-fNnodes)*fBucketSize;
+ Index *indexPoints = &fIndPoints[offset];
+ Int_t nbs = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
+ nbs = nbs ? nbs : fBucketSize;
+ nCalculations += nbs;
+
+ Int_t npoints = 0;
+ for(int idx=0; idx<nbs; idx++){
+ // calculate distance in the L1 metric
+ fDistBuffer[npoints] = 0.;
+ for(int idim=0; idim<fNDim; idim++) fDistBuffer[npoints] += TMath::Abs(point[idim] - fData[idim][indexPoints[idx]]);
+ // register candidate neighbor
+ if(fDistBuffer[npoints] < fkNNdist[k-1]){
+ fIndBuffer[npoints] = indexPoints[idx];
+ npoints++;
}
- memset(fkNN, -1, k*sizeof(Index));
- for(int i=0; i<k; i++) fkNNdist[i] = 9999.;
- Index itmp, jtmp; Value ftmp, gtmp;
+ }
+ for(int ibrowse=0; ibrowse<npoints; ibrowse++){
+ if(fDistBuffer[ibrowse] >= fkNNdist[k-1]) continue;
+ //insert neighbor
+ int iNN=0;
+ while(fDistBuffer[ibrowse] >= fkNNdist[iNN]) iNN++;
+ if(debug>=2) printf("\t\tinsert data %d @ %d distance %f\n", fIndBuffer[ibrowse], iNN, fDistBuffer[ibrowse]);
- // calculate number of boundaries with the data domain.
- if(!fBoundaries) MakeBoundaries();
- Value *bounds = &fBoundaries[inode*2*fNDim];
- Index nBounds = 0;
- for(int idim=0; idim<fNDim; idim++){
- if((bounds[2*idim] - fRange[2*idim]) < 1.E-10) nBounds++;
- if((bounds[2*idim+1] - fRange[2*idim+1]) < 1.E-10) nBounds++;
+ itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
+ fkNN[iNN] = fIndBuffer[ibrowse];
+ fkNNdist[iNN] = fDistBuffer[ibrowse];
+ for(int ireplace=iNN+1; ireplace<k; ireplace++){
+ jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
+ fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp;
+ itmp = jtmp; ftmp = gtmp;
+ if(ftmp == 9999.) break;
}
- if(debug>=1) printf("Calculate boundaries [nBounds = %d]\n", nBounds);
-
-
- // traverse tree
- UChar_t ax;
- Value val, dif;
- Int_t nAllNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
- Int_t nodeStack[128], nodeIn[128];
- Index currentIndex = 0;
- nodeStack[0] = inode;
- nodeIn[0] = inode;
- while(currentIndex>=0){
- Int_t tnode = nodeStack[currentIndex];
- Int_t entry = nodeIn[currentIndex];
- currentIndex--;
- if(debug>=1) printf("Analyse tnode %d entry %d\n", tnode, entry);
-
- // check if node is still eligible
- Int_t inode = (tnode-1)>>1; //tnode/2 + (tnode%2) - 1;
- ax = fAxis[inode];
- val = fValue[inode];
- dif = point[ax] - val;
- if((TMath::Abs(dif) > fkNNdist[k-1]) &&
- ((dif > 0. && (tnode&1)) || (dif < 0. && !(tnode&1)))) {
- if(debug>=1) printf("\tremove %d\n", (tnode-1)>>1);
-
- // mark bound
- nBounds++;
- // end all recursions
- if(nBounds==2 * fNDim) break;
- continue;
- }
-
- if(IsTerminal(tnode)){
- if(debug>=2) printf("\tProcess terminal node %d\n", tnode);
- // Link data to terminal node
- Int_t offset = (tnode >= fCrossNode) ? (tnode-fCrossNode)*fBucketSize : fOffset+(tnode-fNnodes)*fBucketSize;
- Index *indexPoints = &fIndPoints[offset];
- Int_t nbs = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
- nbs = nbs ? nbs : fBucketSize;
- nCalculations += nbs;
-
- Int_t npoints = 0;
- for(int idx=0; idx<nbs; idx++){
- // calculate distance in the L1 metric
- fDistBuffer[npoints] = 0.;
- for(int idim=0; idim<fNDim; idim++) fDistBuffer[npoints] += TMath::Abs(point[idim] - fData[idim][indexPoints[idx]]);
- // register candidate neighbor
- if(fDistBuffer[npoints] < fkNNdist[k-1]){
- fIndBuffer[npoints] = indexPoints[idx];
- npoints++;
- }
- }
- for(int ibrowse=0; ibrowse<npoints; ibrowse++){
- if(fDistBuffer[ibrowse] >= fkNNdist[k-1]) continue;
- //insert neighbor
- int iNN=0;
- while(fDistBuffer[ibrowse] >= fkNNdist[iNN]) iNN++;
- if(debug>=2) printf("\t\tinsert data %d @ %d distance %f\n", fIndBuffer[ibrowse], iNN, fDistBuffer[ibrowse]);
-
- itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
- fkNN[iNN] = fIndBuffer[ibrowse];
- fkNNdist[iNN] = fDistBuffer[ibrowse];
- for(int ireplace=iNN+1; ireplace<k; ireplace++){
- jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
- fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp;
- itmp = jtmp; ftmp = gtmp;
- if(ftmp == 9999.) break;
- }
- if(debug>=3){
- for(int i=0; i<k; i++){
- if(fkNNdist[i] == 9999.) break;
- printf("\t\tfkNNdist[%d] = %f\n", i, fkNNdist[i]);
- }
- }
- }
- }
-
- // register parent
- Int_t pn = (tnode-1)>>1;//tnode/2 + (tnode%2) - 1;
- if(pn >= 0 && entry != pn){
- // check if parent node is eligible at all
- if(TMath::Abs(point[fAxis[pn]] - fValue[pn]) > fkNNdist[k-1]){
- // mark bound
- nBounds++;
- // end all recursions
- if(nBounds==2 * fNDim) break;
- }
- currentIndex++;
- nodeStack[currentIndex]=pn;
- nodeIn[currentIndex]=tnode;
- if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
- }
- if(IsTerminal(tnode)) continue;
-
- // register children nodes
- Int_t cn;
- Bool_t kAllow[] = {kTRUE, kTRUE};
- ax = fAxis[tnode];
- val = fValue[tnode];
- if(TMath::Abs(point[ax] - val) > fkNNdist[k-1]){
- if(point[ax] > val) kAllow[0] = kFALSE;
- else kAllow[1] = kFALSE;
- }
- for(int ic=1; ic<=2; ic++){
- if(!kAllow[ic-1]) continue;
- cn = (tnode<<1)+ic;
- if(cn < nAllNodes && entry != cn){
- currentIndex++;
- nodeStack[currentIndex] = cn;
- nodeIn[currentIndex]=tnode;
- if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
- }
- }
- }
- // save results
- in = fkNN;
- d = fkNNdist;
-
- return kTRUE;
+ if(debug>=3){
+ for(int i=0; i<k; i++){
+ if(fkNNdist[i] == 9999.) break;
+ printf("\t\tfkNNdist[%d] = %f\n", i, fkNNdist[i]);
+ }
+ }
+ }
+ }
+
+ // register parent
+ Int_t pn = (tnode-1)>>1;//tnode/2 + (tnode%2) - 1;
+ if(pn >= 0 && entry != pn){
+ // check if parent node is eligible at all
+ if(TMath::Abs(point[fAxis[pn]] - fValue[pn]) > fkNNdist[k-1]){
+ // mark bound
+ nBounds++;
+ // end all recursions
+ if(nBounds==2 * fNDim) break;
+ }
+ currentIndex++;
+ nodeStack[currentIndex]=pn;
+ nodeIn[currentIndex]=tnode;
+ if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
+ }
+ if(IsTerminal(tnode)) continue;
+
+ // register children nodes
+ Int_t cn;
+ Bool_t kAllow[] = {kTRUE, kTRUE};
+ ax = fAxis[tnode];
+ val = fValue[tnode];
+ if(TMath::Abs(point[ax] - val) > fkNNdist[k-1]){
+ if(point[ax] > val) kAllow[0] = kFALSE;
+ else kAllow[1] = kFALSE;
+ }
+ for(int ic=1; ic<=2; ic++){
+ if(!kAllow[ic-1]) continue;
+ cn = (tnode<<1)+ic;
+ if(cn < nAllNodes && entry != cn){
+ currentIndex++;
+ nodeStack[currentIndex] = cn;
+ nodeIn[currentIndex]=tnode;
+ if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
+ }
+ }
+ }
+ // save results
+ in = fkNN;
+ d = fkNNdist;
+ return kTRUE;
}
//
// find the terminal node to which point belongs
- Index stackNode[128], inode;
- Int_t currentIndex =0;
- stackNode[0] = 0;
- while (currentIndex>=0){
- inode = stackNode[currentIndex];
- if (IsTerminal(inode)) return inode;
-
- currentIndex--;
- if (point[fAxis[inode]]<=fValue[inode]){
- currentIndex++;
- stackNode[currentIndex]=(inode<<1)+1;
- }
- if (point[fAxis[inode]]>=fValue[inode]){
- currentIndex++;
- stackNode[currentIndex]=(inode+1)<<1;
- }
- }
-
- return -1;
+ Index stackNode[128], inode;
+ Int_t currentIndex =0;
+ stackNode[0] = 0;
+ while (currentIndex>=0){
+ inode = stackNode[currentIndex];
+ if (IsTerminal(inode)) return inode;
+
+ currentIndex--;
+ if (point[fAxis[inode]]<=fValue[inode]){
+ currentIndex++;
+ stackNode[currentIndex]=(inode<<1)+1;
+ }
+ if (point[fAxis[inode]]>=fValue[inode]){
+ currentIndex++;
+ stackNode[currentIndex]=(inode+1)<<1;
+ }
+ }
+
+ return -1;
}
// find the index of point
// works only if we keep fData pointers
- Int_t stackNode[128];
- Int_t currentIndex =0;
- stackNode[0] = 0;
- iter =0;
- //
- while (currentIndex>=0){
- iter++;
- Int_t inode = stackNode[currentIndex];
- currentIndex--;
- if (IsTerminal(inode)){
- // investigate terminal node
- Int_t indexIP = (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset;
- printf("terminal %d indexP %d\n", inode, indexIP);
- for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
- Bool_t isOK = kTRUE;
- indexIP+=ibucket;
- printf("ibucket %d index %d\n", ibucket, indexIP);
- if (indexIP>=fNpoints) continue;
- Int_t index0 = fIndPoints[indexIP];
- for (Int_t idim=0;idim<fNDim;idim++) if (fData[idim][index0]!=point[idim]) isOK = kFALSE;
- if (isOK) index = index0;
- }
- continue;
- }
-
- if (point[fAxis[inode]]<=fValue[inode]){
- currentIndex++;
- stackNode[currentIndex]=(inode*2)+1;
- }
- if (point[fAxis[inode]]>=fValue[inode]){
- currentIndex++;
- stackNode[currentIndex]=(inode*2)+2;
- }
- }
+ Int_t stackNode[128];
+ Int_t currentIndex =0;
+ stackNode[0] = 0;
+ iter =0;
+ //
+ while (currentIndex>=0){
+ iter++;
+ Int_t inode = stackNode[currentIndex];
+ currentIndex--;
+ if (IsTerminal(inode)){
+ // investigate terminal node
+ Int_t indexIP = (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset;
+ printf("terminal %d indexP %d\n", inode, indexIP);
+ for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
+ Bool_t isOK = kTRUE;
+ indexIP+=ibucket;
+ printf("ibucket %d index %d\n", ibucket, indexIP);
+ if (indexIP>=fNpoints) continue;
+ Int_t index0 = fIndPoints[indexIP];
+ for (Int_t idim=0;idim<fNDim;idim++) if (fData[idim][index0]!=point[idim]) isOK = kFALSE;
+ if (isOK) index = index0;
+ }
+ continue;
+ }
+
+ if (point[fAxis[inode]]<=fValue[inode]){
+ currentIndex++;
+ stackNode[currentIndex]=(inode*2)+1;
+ }
+ if (point[fAxis[inode]]>=fValue[inode]){
+ currentIndex++;
+ stackNode[currentIndex]=(inode*2)+2;
+ }
+ }
//
// printf("Iter\t%d\n",iter);
}
}
}
+
+
+//______________________________________________________________________
+TKDTreeIF *TKDTreeTestBuild(const Int_t npoints, const Int_t bsize){
+ //
+ // Example function to
+ //
+ Float_t *data0 = new Float_t[npoints*2];
+ Float_t *data[2];
+ data[0] = &data0[0];
+ data[1] = &data0[npoints];
+ for (Int_t i=0;i<npoints;i++) {
+ data[1][i]= gRandom->Rndm();
+ data[0][i]= gRandom->Rndm();
+ }
+ TKDTree<Int_t, Float_t> *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
+ return kdtree;
+}
+
+
+
template class TKDTree<Int_t, Float_t>;
template class TKDTree<Int_t, Double_t>;
template <typename Index, typename Value> class TKDTree : public TObject
{
public:
- enum{
- kDimMax = 6
- };
-
- TKDTree();
- TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data);
- ~TKDTree();
-
- // getters
- Index* GetPointsIndexes(Int_t node) const {
- if(node < fNnodes) return 0x0;
- Int_t offset = (node >= fCrossNode) ? (node-fCrossNode)*fBucketSize : fOffset+(node-fNnodes)*fBucketSize;
- return &fIndPoints[offset];
- }
- UChar_t GetNodeAxis(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fAxis[id];}
- Value GetNodeValue(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fValue[id];}
- Int_t GetNNodes() const {return fNnodes;}
- Value* GetBoundaries();
- Value* GetBoundary(const Int_t node);
- static Int_t GetIndex(Int_t row, Int_t collumn){return collumn+(1<<row);}
- static void GetCoord(Int_t index, Int_t &row, Int_t &collumn){for (row=0; index>=(16<<row);row+=4); for (; index>=(2<<row);row++);collumn= index-(1<<row);};
- Bool_t FindNearestNeighbors(const Value *point, const Int_t kNN, Index *&i, Value *&d);
- Index FindNode(const Value * point);
- void FindPoint(Value * point, Index &index, Int_t &iter);
- void FindInRangeA(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode);
- void FindInRangeB(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode);
- void FindBNodeA(Value * point, Value * delta, Int_t &inode);
- Bool_t IsTerminal(Index inode) const {return (inode>=fNnodes);}
- Value KOrdStat(Index ntotal, Value *a, Index k, Index *index) const;
- void MakeBoundaries(Value *range = 0x0);
- void SetData(Index npoints, Index ndim, UInt_t bsize, Value **data);
- void Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const;
-
-protected:
- void Build(); // build tree
-
-private:
- TKDTree(const TKDTree &); // not implemented
- TKDTree<Index, Value>& operator=(const TKDTree<Index, Value>&); // not implemented
- void CookBoundaries(const Int_t node, Bool_t left);
-
-
-protected:
- Bool_t fDataOwner; //! Toggle ownership of the data
- Int_t fNnodes; // size of node array
- Index fNDim; // number of dimensions
- Index fNDimm; // dummy 2*fNDim
- Index fNpoints; // number of multidimensional points
- Index fBucketSize; // limit statistic for nodes
- UChar_t *fAxis; //[fNnodes] nodes cutting axis
- Value *fValue; //[fNnodes] nodes cutting value
- Value *fRange; //[fNDimm] range of data for each dimension
- Value **fData; //! data points
- Value *fBoundaries;//! nodes boundaries - check class doc
-
-// kNN related data
-protected:
- Int_t fkNNdim; //! current kNN arrays allocated dimension
- Index *fkNN; //! k nearest neighbors indexes
- Value *fkNNdist; //! k nearest neighbors distances
- Value *fDistBuffer;//! working space for kNN
- Index *fIndBuffer; //! working space for kNN
-private:
- Index *fIndPoints; //! array of points indexes
- Int_t fRowT0; //! smallest terminal row
- Int_t fCrossNode; //! cross node
- Int_t fOffset; //! offset in fIndPoints
-
- ClassDef(TKDTree, 1) // KD tree
+ TKDTree();
+ TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data);
+ ~TKDTree();
+ // Get indexes of left and right daughter nodes
+ Int_t GetLeft(Int_t inode) const {return inode*2+1;}
+ Int_t GetRight(Int_t inode) const {return (inode+1)*2;}
+ Int_t GetParent(Int_t inode) const {return inode/2;}
+ //
+ // getters
+ Index* GetPointsIndexes(Int_t node) const {
+ if(node < fNnodes) return 0x0;
+ Int_t offset = (node >= fCrossNode) ? (node-fCrossNode)*fBucketSize : fOffset+(node-fNnodes)*fBucketSize;
+ return &fIndPoints[offset];
+ }
+ UChar_t GetNodeAxis(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fAxis[id];}
+ Value GetNodeValue(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fValue[id];}
+ Int_t GetNNodes() const {return fNnodes;}
+ Value* GetBoundaries();
+ Value* GetBoundary(const Int_t node);
+ static Int_t GetIndex(Int_t row, Int_t collumn){return collumn+(1<<row);}
+ static void GetCoord(Int_t index, Int_t &row, Int_t &collumn){for (row=0; index>=(16<<row);row+=4); for (; index>=(2<<row);row++);collumn= index-(1<<row);};
+ Bool_t FindNearestNeighbors(const Value *point, const Int_t kNN, Index *&i, Value *&d);
+ Index FindNode(const Value * point);
+ void FindPoint(Value * point, Index &index, Int_t &iter);
+ void FindInRangeA(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode);
+ void FindInRangeB(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode);
+ void FindBNodeA(Value * point, Value * delta, Int_t &inode);
+ Bool_t IsTerminal(Index inode) const {return (inode>=fNnodes);}
+ Value KOrdStat(Index ntotal, Value *a, Index k, Index *index) const;
+ void MakeBoundaries(Value *range = 0x0);
+ void SetData(Index npoints, Index ndim, UInt_t bsize, Value **data);
+ void Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const;
+
+
+ protected:
+ void Build(); // build tree
+
+ private:
+ TKDTree(const TKDTree &); // not implemented
+ TKDTree<Index, Value>& operator=(const TKDTree<Index, Value>&); // not implemented
+ void CookBoundaries(const Int_t node, Bool_t left);
+
+
+ protected:
+ Bool_t fDataOwner; //! Toggle ownership of the data
+ Int_t fNnodes; // size of node array
+ Index fNDim; // number of dimensions
+ Index fNDimm; // dummy 2*fNDim
+ Index fNpoints; // number of multidimensional points
+ Index fBucketSize; // limit statistic for nodes
+ UChar_t *fAxis; //[fNnodes] nodes cutting axis
+ Value *fValue; //[fNnodes] nodes cutting value
+ //
+ Value *fRange; //[fNDimm] range of data for each dimension
+ Value **fData; //! data points
+ Value *fBoundaries;//! nodes boundaries - check class doc
+
+ // kNN related data
+ protected:
+ Int_t fkNNdim; //! current kNN arrays allocated dimension
+ Index *fkNN; //! k nearest neighbors indexes
+ Value *fkNNdist; //! k nearest neighbors distances
+ Value *fDistBuffer;//! working space for kNN
+ Index *fIndBuffer; //! working space for kNN
+
+ private:
+ Index *fIndPoints; //! array of points indexes
+ Int_t fRowT0; //! smallest terminal row
+ Int_t fCrossNode; //! cross node
+ Int_t fOffset; //! offset in fIndPoints
+
+ ClassDef(TKDTree, 1) // KD tree
};
typedef TKDTree<Int_t, Double_t> TKDTreeID;
typedef TKDTree<Int_t, Float_t> TKDTreeIF;
+// Test functions:
+TKDTreeIF * TKDTreeTestBuild();
+
+
+
//_________________________________________________________________
template <typename Index, typename Value>
void TKDTree<Index, Value>::FindBNodeA(Value *point, Value *delta, Int_t &inode){
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff