// - after are the terminal nodes
// fNodes : array of primary nodes
///////////////////////////////////////////////////////////////////////
+
+
//_________________________________________________________________
template <typename Index, typename Value>
TKDTree<Index, Value>::TKDTree() :
,kDataOwner(kFALSE)
,fNnodes(0)
,fNDim(0)
+ ,fNDimm(0)
,fNpoints(0)
,fBucketSize(0)
- ,fData(0x0)
+ ,fAxis(0x0)
+ ,fValue(0x0)
,fRange(0x0)
+ ,fData(0x0)
,fBoundaries(0x0)
- ,fNodes(0x0)
+ ,fkNNdim(0)
,fkNN(0x0)
+ ,fkNNdist(0x0)
+ ,fDistBuffer(0x0)
+ ,fIndBuffer(0x0)
,fIndPoints(0x0)
,fRowT0(0)
,fCrossNode(0)
,kDataOwner(kFALSE)
,fNnodes(0)
,fNDim(ndim)
+ ,fNDimm(2*ndim)
,fNpoints(npoints)
,fBucketSize(bsize)
- ,fData(data) //Columnwise!!!!!
+ ,fAxis(0x0)
+ ,fValue(0x0)
,fRange(0x0)
+ ,fData(data) //Columnwise!!!!!
,fBoundaries(0x0)
- ,fNodes(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.
-
+// Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
+
Build();
}
template <typename Index, typename Value>
TKDTree<Index, Value>::~TKDTree()
{
+ if (fIndBuffer) delete [] fIndBuffer;
+ if (fDistBuffer) delete [] fDistBuffer;
+ if (fkNNdist) delete [] fkNNdist;
if (fkNN) delete [] fkNN;
- if (fNodes) delete [] fNodes;
+ if (fAxis) delete [] fAxis;
+ if (fValue) delete [] fValue;
if (fIndPoints) delete [] fIndPoints;
if (fRange) delete [] fRange;
if (fBoundaries) delete [] fBoundaries;
fRange = new Value[2*fNDim];
fIndPoints= new Index[fNpoints];
for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i;
- fNodes = new TKDNode[fNnodes];
+ fAxis = new UChar_t[fNnodes];
+ fValue = new Value[fNnodes];
//
fCrossNode = (1<<(fRowT0+1))-1;
if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1;
Int_t crow = rowStack[currentIndex];
Int_t cpos = posStack[currentIndex];
Int_t cnode = nodeStack[currentIndex];
- TKDNode * node = &(fNodes[cnode]);
//printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
//
// divide points
maxspread=tempspread;
axspread = idim;
}
- if(cnode==0) {fRange[2*idim] = min; fRange[2*idim+1] = max;}
+ 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);
- node->fAxis = axspread;
- node->fValue = array[fIndPoints[cpos+nleft]];
+ fAxis[cnode] = axspread;
+ fValue[cnode] = array[fIndPoints[cpos+nleft]];
//printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
//
//
//_________________________________________________________________
template <typename Index, typename Value>
-Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value &d)
+Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value *&d)
{
// Find "k" nearest neighbors to "point".
//
// increasing order of the distance to "point" and the maxim distance
// in "d".
//
-// The array "in" is managed internally by the TKDTree.
+// The arrays "in" and "d" are managed internally by the TKDTree.
Index inode = FindNode(point);
if(inode < fNnodes){
return kFALSE;
}
+ UInt_t debug = 0;
+ Int_t nCalculations = 0;
+
// allocate working memory
- if(!fkNN) fkNN = new Index[k];
- if(sizeof(fkNN) < k*sizeof(Index)) fkNN = new(fkNN) Index[k];
- for(int i=0; i<k; i++) fkNN[i] = -1;
- Index *fkNN_tmp = new Index[k];
- Value *fkNN_dist = new Value[k];
- for(int i=0; i<k; i++) fkNN_dist[i] = 9999.;
- Value *fkNN_dist_tmp = new Value[k];
- Value *dist = new Value[fBucketSize];
- Index *index = new Index[fBucketSize];
-
+ 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.
- Index nBounds = 0;
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]) nBounds++;
- if(bounds[2*idim+1] == fRange[2*idim+1]) nBounds++;
+ 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
- TKDNode *node;
+ 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;
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
- node = &fNodes[tnode/2 + (tnode%2) - 1];
- if((TMath::Abs(point[node->fAxis] - node->fValue) > fkNN_dist[k-1])
- &&
- ((point[node->fAxis] > node->fValue && tnode%2) ||
- (point[node->fAxis] < node->fValue && !tnode%2))) {
- //printf("\tREMOVE NODE\n");
+ 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++;
}
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 nPoints = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
- for(int idx=0; idx<nPoints; idx++){
+ 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
- dist[idx] = 0.;
- for(int idim=0; idim<fNDim; idim++) dist[idx] += TMath::Abs(point[idim] - fData[idim][indexPoints[idx]]);
+ 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++;
+ }
}
- // arrange increasingly distances array and update kNN indexes
- TMath::Sort(nPoints, dist, index, kFALSE);
- for(int ibrowse=0; ibrowse<nPoints; ibrowse++){
- // exit if the current distance calculated in this node is
- // larger than the largest distance stored in the kNN array
- if(dist[index[ibrowse]] >= fkNN_dist[k-1]) break;
- for(int iNN=0; iNN<k; iNN++){
- if(dist[index[ibrowse]] < fkNN_dist[iNN]){
- //insert neighbor. In principle this is only one call to STL vector. Maybe we can change to it ?!
- //printf("\t\tinsert data %d @ %d distance %f\n", indexPoints[index[ibrowse]], iNN, dist[index[ibrowse]]);
-
- memcpy(fkNN_tmp, &fkNN[iNN], (k-iNN-1)*sizeof(Index));
- fkNN[iNN] = indexPoints[index[ibrowse]];
- memcpy(&fkNN[iNN+1], fkNN_tmp, (k-iNN-1)*sizeof(Index));
-
- memcpy(fkNN_dist_tmp, &fkNN_dist[iNN], (k-iNN-1)*sizeof(Value));
- fkNN_dist[iNN] = dist[index[ibrowse]];
- memcpy(&fkNN_dist[iNN+1], fkNN_dist_tmp, (k-iNN-1)*sizeof(Value));
- break;
- }
+ 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 parent_node = tnode/2 + (tnode%2) - 1;
- if(parent_node >= 0 && entry != parent_node){
+ Int_t pn = (tnode-1)>>1;//tnode/2 + (tnode%2) - 1;
+ if(pn >= 0 && entry != pn){
// check if parent node is eligible at all
- node = &fNodes[parent_node];
- if(TMath::Abs(point[node->fAxis] - node->fValue) > fkNN_dist[k-1]){
+ 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]=tnode/2 + (tnode%2) - 1;
+ nodeStack[currentIndex]=pn;
nodeIn[currentIndex]=tnode;
- //printf("\tregister %d\n", nodeStack[currentIndex]);
+ if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
}
-
+ if(IsTerminal(tnode)) continue;
+
// register children nodes
- Int_t child_node;
+ Int_t cn;
Bool_t kAllow[] = {kTRUE, kTRUE};
- node = &fNodes[tnode];
- if(TMath::Abs(point[node->fAxis] - node->fValue) > fkNN_dist[k-1]){
- if(point[node->fAxis] > node->fValue) kAllow[0] = kFALSE;
+ 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;
- child_node = (tnode*2)+ic;
- if(child_node < fNnodes + GetNTerminalNodes() && entry != child_node){
+ cn = (tnode<<1)+ic;
+ if(cn < nAllNodes && entry != cn){
currentIndex++;
- nodeStack[currentIndex] = child_node;
+ nodeStack[currentIndex] = cn;
nodeIn[currentIndex]=tnode;
- //printf("\tregister %d\n", nodeStack[currentIndex]);
+ if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
}
- }
+ }
}
// save results
in = fkNN;
- d = fkNN_dist[k-1];
+ d = fkNNdist;
- delete [] index;
- delete [] dist;
- delete [] fkNN_tmp;
- delete [] fkNN_dist;
- delete [] fkNN_dist_tmp;
-
return kTRUE;
}
Index TKDTree<Index, Value>::FindNode(const Value * point){
//
// 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];
- currentIndex--;
if (IsTerminal(inode)) return inode;
- TKDNode & node = fNodes[inode];
- if (point[node.fAxis]<=node.fValue){
+ currentIndex--;
+ if (point[fAxis[inode]]<=fValue[inode]){
currentIndex++;
- stackNode[currentIndex]=(inode*2)+1;
+ stackNode[currentIndex]=(inode<<1)+1;
}
- if (point[node.fAxis]>=node.fValue){
+ if (point[fAxis[inode]]>=fValue[inode]){
currentIndex++;
- stackNode[currentIndex]=(inode*2)+2;
+ stackNode[currentIndex]=(inode+1)<<1;
}
}
continue;
}
- TKDNode & node = fNodes[inode];
- if (point[node.fAxis]<=node.fValue){
+ if (point[fAxis[inode]]<=fValue[inode]){
currentIndex++;
stackNode[currentIndex]=(inode*2)+1;
}
- if (point[node.fAxis]>=node.fValue){
+ if (point[fAxis[inode]]>=fValue[inode]){
currentIndex++;
stackNode[currentIndex]=(inode*2)+2;
}
currentIndex--;
if (!IsTerminal(inode)){
// not terminal
- TKDNode * node = &(fNodes[inode]);
- if (point[node->fAxis] - delta[node->fAxis] < node->fValue) {
+ //TKDNode<Index, Value> * node = &(fNodes[inode]);
+ if (point[fAxis[inode]] - delta[fAxis[inode]] < fValue[inode]) {
currentIndex++;
stackNode[currentIndex]= (inode*2)+1;
}
- if (point[node->fAxis] + delta[node->fAxis] >= node->fValue){
+ if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
currentIndex++;
stackNode[currentIndex]= (inode*2)+2;
}
}
}
//
- TKDNode * node = &(fNodes[inode]);
- if (point[node->fAxis] - delta[node->fAxis] < node->fValue) {
+ // TKDNode<Index, Value> * node = &(fNodes[inode]);
+ if (point[fAxis[inode]] - delta[fAxis[inode]] < fValue[inode]) {
currentIndex++;
stackNode[currentIndex]= (inode<<1)+1;
- if (point[node->fAxis] + delta[node->fAxis] > node->fValue)
- stackStatus[currentIndex]= status | (1<<(2*node->fAxis));
+ if (point[fAxis[inode]] + delta[fAxis[inode]] > fValue[inode])
+ stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]));
}
- if (point[node->fAxis] + delta[node->fAxis] >= node->fValue){
+ if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
currentIndex++;
stackNode[currentIndex]= (inode<<1)+2;
- if (point[node->fAxis] - delta[node->fAxis]<node->fValue)
- stackStatus[currentIndex]= status | (1<<(2*node->fAxis+1));
+ if (point[fAxis[inode]] - delta[fAxis[inode]]<fValue[inode])
+ stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]+1));
}
}
if (fData){
// Build boundaries for each node.
- if(range) memcpy(fRange, range, 2*fNDim*sizeof(Value));
+ if(range) memcpy(fRange, range, fNDimm*sizeof(Value));
// total number of nodes including terminal nodes
- Int_t totNodes = fNnodes + GetNTerminalNodes();
- fBoundaries = new Value[totNodes*2*fNDim];
- printf("Allocate %d nodes\n", totNodes);
+ Int_t totNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
+ fBoundaries = new Value[totNodes*fNDimm];
+ //Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes));
+
+
// loop
- Int_t child_index;
+ Value *tbounds = 0x0, *cbounds = 0x0;
+ Int_t cn;
for(int inode=fNnodes-1; inode>=0; inode--){
- //printf("\tAllocate node %d\n", inode);
- memcpy(&fBoundaries[inode*2*fNDim], fRange, 2*fNDim*sizeof(Value));
+ tbounds = &fBoundaries[inode*fNDimm];
+ memcpy(tbounds, fRange, fNDimm*sizeof(Value));
// check left child node
- child_index = 2*inode+1;
- if(IsTerminal(child_index)) CookBoundariesTerminal(inode, kTRUE);
- for(int idim=0; idim<fNDim; idim++) fBoundaries[2*fNDim*inode+2*idim] = fBoundaries[2*fNDim*child_index+2*idim];
+ cn = (inode<<1)+1;
+ if(IsTerminal(cn)) CookBoundaries(inode, kTRUE);
+ cbounds = &fBoundaries[fNDimm*cn];
+ for(int idim=0; idim<fNDim; idim++) tbounds[idim<<1] = cbounds[idim<<1];
// check right child node
- child_index = 2*inode+2;
- if(IsTerminal(child_index)) CookBoundariesTerminal(inode, kFALSE);
- for(int idim=0; idim<fNDim; idim++) fBoundaries[2*fNDim*inode+2*idim+1] = fBoundaries[2*fNDim*child_index+2*idim+1];
+ cn = (inode+1)<<1;
+ if(IsTerminal(cn)) CookBoundaries(inode, kFALSE);
+ cbounds = &fBoundaries[fNDimm*cn];
+ for(int idim=0; idim<fNDim; idim++) tbounds[(idim<<1)+1] = cbounds[(idim<<1)+1];
}
}
//_________________________________________________________________
template <typename Index, typename Value>
-void TKDTree<Index, Value>::CookBoundariesTerminal(Int_t parent_node, Bool_t LEFT)
+void TKDTree<Index, Value>::CookBoundaries(const Int_t node, Bool_t LEFT)
{
// define index of this terminal node
- Int_t index = 2 * parent_node + (LEFT ? 1 : 2);
+ Int_t index = (node<<1) + (LEFT ? 1 : 2);
+ //Info("CookBoundaries()", Form("Node %d", index));
// build and initialize boundaries for this node
- //printf("\tAllocate terminal node %d\n", index);
- memcpy(&fBoundaries[2*fNDim*index], fRange, 2*fNDim*sizeof(Value));
+ Value *tbounds = &fBoundaries[fNDimm*index];
+ memcpy(tbounds, fRange, fNDimm*sizeof(Value));
Bool_t flag[256]; // cope with up to 128 dimensions
- memset(flag, kFALSE, 2*fNDim);
+ memset(flag, kFALSE, fNDimm);
Int_t nvals = 0;
// recurse parent nodes
- TKDNode *node = 0x0;
- while(parent_node >= 0 && nvals < 2 * fNDim){
- node = &(fNodes[parent_node]);
+ Int_t pn = node;
+ while(pn >= 0 && nvals < fNDimm){
if(LEFT){
- if(!flag[2*node->fAxis+1]) {
- fBoundaries[2*fNDim*index+2*node->fAxis+1] = node->fValue;
- flag[2*node->fAxis+1] = kTRUE;
+ index = (fAxis[pn]<<1)+1;
+ if(!flag[index]) {
+ tbounds[index] = fValue[pn];
+ flag[index] = kTRUE;
nvals++;
}
} else {
- if(!flag[2*node->fAxis]) {
- fBoundaries[2*fNDim*index+2*node->fAxis] = node->fValue;
- flag[2*node->fAxis] = kTRUE;
+ index = fAxis[pn]<<1;
+ if(!flag[index]) {
+ tbounds[index] = fValue[pn];
+ flag[index] = kTRUE;
nvals++;
}
}
- LEFT = parent_node%2;
- parent_node = parent_node/2 + (parent_node%2) - 1;
+ LEFT = pn&1;
+ pn = (pn - 1)>>1;
}
}