]>
Commit | Line | Data |
---|---|---|
f2040a8f | 1 | #include "TKDTree.h" |
2 | ||
3 | #include "TString.h" | |
4 | #include <string.h> | |
5 | ||
6 | #ifndef R__ALPHA | |
7 | templateClassImp(TKDTree) | |
8 | #endif | |
9 | ////////////////////////////////////////////////////////////////////// | |
10 | // Memory setup of protected data members: | |
11 | // | |
12 | // kDataOwner; // Toggle ownership of the data | |
13 | // fNnodes: | |
14 | // size of branch node array, and index of first terminal node | |
15 | // fNDim; // number of variables | |
16 | // fNpoints; // number of multidimensional points | |
17 | // fBucketSize; // number of data points per terminal node | |
18 | // | |
19 | // fIndPoints : array containing rearanged indexes according to nodes: | |
20 | // | point indexes from 1st terminal node (fBucketSize elements) | ... | point indexes from the last terminal node (<= fBucketSize elements) | |
21 | // | |
22 | // Value **fData; | |
23 | // fRange : array containing the boundaries of the domain: | |
24 | // | 1st dimension (min + max) | 2nd dimension (min + max) | ... | |
25 | // fBoundaries : nodes boundaries | |
26 | // | 1st node {1st dim * 2 elements | 2nd dim * 2 elements | ...} | 2nd node {...} | ... | |
27 | // the nodes are arranged in the following order : | |
28 | // - first fNnodes are primary nodes | |
29 | // - after are the terminal nodes | |
30 | // fNodes : array of primary nodes | |
31 | /////////////////////////////////////////////////////////////////////// | |
117c62ab | 32 | |
117c62ab | 33 | |
f2040a8f | 34 | //_________________________________________________________________ |
35 | template <typename Index, typename Value> | |
36 | TKDTree<Index, Value>::TKDTree() : | |
37 | TObject() | |
38 | ,kDataOwner(kFALSE) | |
39 | ,fNnodes(0) | |
40 | ,fNDim(0) | |
316a7f5a | 41 | ,fNDimm(0) |
f2040a8f | 42 | ,fNpoints(0) |
43 | ,fBucketSize(0) | |
117c62ab | 44 | ,fAxis(0x0) |
45 | ,fValue(0x0) | |
f2040a8f | 46 | ,fRange(0x0) |
316a7f5a | 47 | ,fData(0x0) |
f2040a8f | 48 | ,fBoundaries(0x0) |
316a7f5a | 49 | ,fkNNdim(0) |
f2040a8f | 50 | ,fkNN(0x0) |
316a7f5a | 51 | ,fkNNdist(0x0) |
52 | ,fDistBuffer(0x0) | |
53 | ,fIndBuffer(0x0) | |
a9c20b1f | 54 | ,fIndPoints(0x0) |
55 | ,fRowT0(0) | |
56 | ,fCrossNode(0) | |
57 | ,fOffset(0) | |
f2040a8f | 58 | { |
59 | // Default constructor. Nothing is built | |
60 | } | |
61 | ||
62 | ||
63 | //_________________________________________________________________ | |
64 | template <typename Index, typename Value> | |
65 | TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data) : | |
66 | TObject() | |
67 | ,kDataOwner(kFALSE) | |
68 | ,fNnodes(0) | |
69 | ,fNDim(ndim) | |
316a7f5a | 70 | ,fNDimm(2*ndim) |
f2040a8f | 71 | ,fNpoints(npoints) |
72 | ,fBucketSize(bsize) | |
117c62ab | 73 | ,fAxis(0x0) |
74 | ,fValue(0x0) | |
f2040a8f | 75 | ,fRange(0x0) |
316a7f5a | 76 | ,fData(data) //Columnwise!!!!! |
f2040a8f | 77 | ,fBoundaries(0x0) |
316a7f5a | 78 | ,fkNNdim(0) |
f2040a8f | 79 | ,fkNN(0x0) |
316a7f5a | 80 | ,fkNNdist(0x0) |
81 | ,fDistBuffer(0x0) | |
82 | ,fIndBuffer(0x0) | |
a9c20b1f | 83 | ,fIndPoints(0x0) |
84 | ,fRowT0(0) | |
85 | ,fCrossNode(0) | |
86 | ,fOffset(0) | |
f2040a8f | 87 | { |
316a7f5a | 88 | // Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method. |
89 | ||
f2040a8f | 90 | Build(); |
91 | } | |
92 | ||
f2040a8f | 93 | //_________________________________________________________________ |
94 | template <typename Index, typename Value> | |
95 | TKDTree<Index, Value>::~TKDTree() | |
96 | { | |
316a7f5a | 97 | if (fIndBuffer) delete [] fIndBuffer; |
98 | if (fDistBuffer) delete [] fDistBuffer; | |
99 | if (fkNNdist) delete [] fkNNdist; | |
f2040a8f | 100 | if (fkNN) delete [] fkNN; |
117c62ab | 101 | if (fAxis) delete [] fAxis; |
102 | if (fValue) delete [] fValue; | |
f2040a8f | 103 | if (fIndPoints) delete [] fIndPoints; |
104 | if (fRange) delete [] fRange; | |
105 | if (fBoundaries) delete [] fBoundaries; | |
106 | if (kDataOwner && fData){ | |
107 | for(int idim=0; idim<fNDim; idim++) delete [] fData[idim]; | |
108 | delete [] fData; | |
109 | } | |
110 | } | |
111 | ||
112 | ||
113 | //_________________________________________________________________ | |
114 | template <typename Index, typename Value> | |
115 | void TKDTree<Index, Value>::Build(){ | |
116 | // | |
117 | // Build binning | |
118 | // | |
119 | // 1. calculate number of nodes | |
120 | // 2. calculate first terminal row | |
121 | // 3. initialize index array | |
122 | // 4. non recursive building of the binary tree | |
123 | // | |
124 | //1. | |
125 | fNnodes = fNpoints/fBucketSize-1; | |
126 | if (fNpoints%fBucketSize) fNnodes++; | |
127 | ||
128 | //2. | |
129 | fRowT0=0; | |
130 | for (;(fNnodes+1)>(1<<fRowT0);fRowT0++); | |
131 | fRowT0-=1; | |
132 | // 2 = 2**0 + 1 | |
133 | // 3 = 2**1 + 1 | |
134 | // 4 = 2**1 + 2 | |
135 | // 5 = 2**2 + 1 | |
136 | // 6 = 2**2 + 2 | |
137 | // 7 = 2**2 + 3 | |
138 | // 8 = 2**2 + 4 | |
139 | ||
140 | //3. | |
141 | // allocate space for boundaries | |
142 | fRange = new Value[2*fNDim]; | |
143 | fIndPoints= new Index[fNpoints]; | |
144 | for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i; | |
117c62ab | 145 | fAxis = new UChar_t[fNnodes]; |
146 | fValue = new Value[fNnodes]; | |
f2040a8f | 147 | // |
148 | fCrossNode = (1<<(fRowT0+1))-1; | |
149 | if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1; | |
150 | // | |
151 | // fOffset = (((fNnodes+1)-(1<<fRowT0)))*2; | |
152 | Int_t over = (fNnodes+1)-(1<<fRowT0); | |
153 | Int_t filled = ((1<<fRowT0)-over)*fBucketSize; | |
154 | fOffset = fNpoints-filled; | |
155 | // | |
156 | // printf("Row0 %d\n", fRowT0); | |
157 | // printf("CrossNode %d\n", fCrossNode); | |
158 | // printf("Offset %d\n", fOffset); | |
159 | // | |
160 | // | |
161 | //4. | |
162 | // stack for non recursive build - size 128 bytes enough | |
163 | Int_t rowStack[128]; | |
164 | Int_t nodeStack[128]; | |
165 | Int_t npointStack[128]; | |
166 | Int_t posStack[128]; | |
167 | Int_t currentIndex = 0; | |
168 | Int_t iter =0; | |
169 | rowStack[0] = 0; | |
170 | nodeStack[0] = 0; | |
171 | npointStack[0] = fNpoints; | |
172 | posStack[0] = 0; | |
173 | // | |
174 | Int_t nbucketsall =0; | |
175 | while (currentIndex>=0){ | |
176 | iter++; | |
177 | // | |
178 | Int_t npoints = npointStack[currentIndex]; | |
179 | if (npoints<=fBucketSize) { | |
180 | //printf("terminal node : index %d iter %d\n", currentIndex, iter); | |
181 | currentIndex--; | |
182 | nbucketsall++; | |
183 | continue; // terminal node | |
184 | } | |
185 | Int_t crow = rowStack[currentIndex]; | |
186 | Int_t cpos = posStack[currentIndex]; | |
187 | Int_t cnode = nodeStack[currentIndex]; | |
f2040a8f | 188 | //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode); |
189 | // | |
190 | // divide points | |
191 | Int_t nbuckets0 = npoints/fBucketSize; //current number of buckets | |
192 | if (npoints%fBucketSize) nbuckets0++; // | |
193 | Int_t restRows = fRowT0-rowStack[currentIndex]; // rest of fully occupied node row | |
194 | if (restRows<0) restRows =0; | |
195 | for (;nbuckets0>(2<<restRows); restRows++); | |
196 | Int_t nfull = 1<<restRows; | |
197 | Int_t nrest = nbuckets0-nfull; | |
198 | Int_t nleft =0, nright =0; | |
199 | // | |
200 | if (nrest>(nfull/2)){ | |
201 | nleft = nfull*fBucketSize; | |
202 | nright = npoints-nleft; | |
203 | }else{ | |
204 | nright = nfull*fBucketSize/2; | |
205 | nleft = npoints-nright; | |
206 | } | |
207 | ||
208 | // | |
209 | //find the axis with biggest spread | |
210 | Value maxspread=0; | |
211 | Value tempspread, min, max; | |
212 | Index axspread=0; | |
213 | Value *array; | |
214 | for (Int_t idim=0; idim<fNDim; idim++){ | |
215 | array = fData[idim]; | |
216 | Spread(npoints, array, fIndPoints+cpos, min, max); | |
217 | tempspread = max - min; | |
218 | if (maxspread < tempspread) { | |
219 | maxspread=tempspread; | |
220 | axspread = idim; | |
221 | } | |
117c62ab | 222 | if(cnode) continue; |
223 | //printf("set %d %6.3f %6.3f\n", idim, min, max); | |
224 | fRange[2*idim] = min; fRange[2*idim+1] = max; | |
f2040a8f | 225 | } |
226 | array = fData[axspread]; | |
227 | KOrdStat(npoints, array, nleft, fIndPoints+cpos); | |
117c62ab | 228 | fAxis[cnode] = axspread; |
229 | fValue[cnode] = array[fIndPoints[cpos+nleft]]; | |
f2040a8f | 230 | //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue); |
231 | // | |
232 | // | |
233 | npointStack[currentIndex] = nleft; | |
234 | rowStack[currentIndex] = crow+1; | |
235 | posStack[currentIndex] = cpos; | |
236 | nodeStack[currentIndex] = cnode*2+1; | |
237 | currentIndex++; | |
238 | npointStack[currentIndex] = nright; | |
239 | rowStack[currentIndex] = crow+1; | |
240 | posStack[currentIndex] = cpos+nleft; | |
241 | nodeStack[currentIndex] = (cnode*2)+2; | |
242 | // | |
243 | if (0){ | |
244 | // consistency check | |
245 | Info("Build()", Form("points %d left %d right %d", npoints, nleft, nright)); | |
246 | if (nleft<nright) Warning("Build", "Problem Left-Right"); | |
247 | if (nleft<0 || nright<0) Warning("Build()", "Problem Negative number"); | |
248 | } | |
249 | } | |
f2040a8f | 250 | |
251 | //printf("NBuckets\t%d\n", nbucketsall); | |
252 | //fData = 0; | |
253 | } | |
254 | ||
255 | ||
256 | //_________________________________________________________________ | |
257 | template <typename Index, typename Value> | |
117c62ab | 258 | Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value *&d) |
f2040a8f | 259 | { |
260 | // Find "k" nearest neighbors to "point". | |
261 | // | |
262 | // Return true in case of success and false in case of failure. | |
263 | // The indexes of the nearest k points are stored in the array "in" in | |
264 | // increasing order of the distance to "point" and the maxim distance | |
265 | // in "d". | |
266 | // | |
117c62ab | 267 | // The arrays "in" and "d" are managed internally by the TKDTree. |
f2040a8f | 268 | |
269 | Index inode = FindNode(point); | |
270 | if(inode < fNnodes){ | |
271 | Error("FindNearestNeighbors()", "Point outside data range."); | |
272 | return kFALSE; | |
273 | } | |
274 | ||
316a7f5a | 275 | UInt_t debug = 0; |
276 | Int_t nCalculations = 0; | |
277 | ||
f2040a8f | 278 | // allocate working memory |
316a7f5a | 279 | if(!fDistBuffer){ |
280 | fDistBuffer = new Value[fBucketSize]; | |
281 | fIndBuffer = new Index[fBucketSize]; | |
282 | } | |
283 | if(fkNNdim < k){ | |
117c62ab | 284 | if(debug>=1){ |
316a7f5a | 285 | if(fkNN) printf("Reallocate memory %d -> %d \n", fkNNdim, 2*k); |
286 | else printf("Allocate %d memory\n", 2*k); | |
117c62ab | 287 | } |
316a7f5a | 288 | fkNNdim = 2*k; |
289 | if(fkNN){ | |
290 | delete [] fkNN; | |
291 | delete [] fkNNdist; | |
292 | } | |
293 | fkNN = new Index[fkNNdim]; | |
294 | fkNNdist = new Value[fkNNdim]; | |
295 | } | |
296 | memset(fkNN, -1, k*sizeof(Index)); | |
297 | for(int i=0; i<k; i++) fkNNdist[i] = 9999.; | |
298 | Index itmp, jtmp; Value ftmp, gtmp; | |
299 | ||
f2040a8f | 300 | // calculate number of boundaries with the data domain. |
f2040a8f | 301 | if(!fBoundaries) MakeBoundaries(); |
302 | Value *bounds = &fBoundaries[inode*2*fNDim]; | |
117c62ab | 303 | Index nBounds = 0; |
f2040a8f | 304 | for(int idim=0; idim<fNDim; idim++){ |
117c62ab | 305 | if((bounds[2*idim] - fRange[2*idim]) < 1.E-10) nBounds++; |
306 | if((bounds[2*idim+1] - fRange[2*idim+1]) < 1.E-10) nBounds++; | |
f2040a8f | 307 | } |
316a7f5a | 308 | if(debug>=1) printf("Calculate boundaries [nBounds = %d]\n", nBounds); |
117c62ab | 309 | |
310 | ||
f2040a8f | 311 | // traverse tree |
117c62ab | 312 | UChar_t ax; |
313 | Value val, dif; | |
a3408ed3 | 314 | Int_t nAllNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0); |
f2040a8f | 315 | Int_t nodeStack[128], nodeIn[128]; |
316 | Index currentIndex = 0; | |
317 | nodeStack[0] = inode; | |
318 | nodeIn[0] = inode; | |
319 | while(currentIndex>=0){ | |
320 | Int_t tnode = nodeStack[currentIndex]; | |
321 | Int_t entry = nodeIn[currentIndex]; | |
322 | currentIndex--; | |
316a7f5a | 323 | if(debug>=1) printf("Analyse tnode %d entry %d\n", tnode, entry); |
324 | ||
f2040a8f | 325 | // check if node is still eligible |
117c62ab | 326 | Int_t inode = (tnode-1)>>1; //tnode/2 + (tnode%2) - 1; |
327 | ax = fAxis[inode]; | |
328 | val = fValue[inode]; | |
329 | dif = point[ax] - val; | |
330 | if((TMath::Abs(dif) > fkNNdist[k-1]) && | |
331 | ((dif > 0. && (tnode&1)) || (dif < 0. && !(tnode&1)))) { | |
332 | if(debug>=1) printf("\tremove %d\n", (tnode-1)>>1); | |
f2040a8f | 333 | |
334 | // mark bound | |
335 | nBounds++; | |
336 | // end all recursions | |
337 | if(nBounds==2 * fNDim) break; | |
338 | continue; | |
339 | } | |
340 | ||
341 | if(IsTerminal(tnode)){ | |
316a7f5a | 342 | if(debug>=2) printf("\tProcess terminal node %d\n", tnode); |
f2040a8f | 343 | // Link data to terminal node |
344 | Int_t offset = (tnode >= fCrossNode) ? (tnode-fCrossNode)*fBucketSize : fOffset+(tnode-fNnodes)*fBucketSize; | |
345 | Index *indexPoints = &fIndPoints[offset]; | |
316a7f5a | 346 | Int_t nbs = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize; |
347 | nbs = nbs ? nbs : fBucketSize; | |
348 | nCalculations += nbs; | |
349 | ||
350 | Int_t npoints = 0; | |
351 | for(int idx=0; idx<nbs; idx++){ | |
f2040a8f | 352 | // calculate distance in the L1 metric |
316a7f5a | 353 | fDistBuffer[npoints] = 0.; |
354 | for(int idim=0; idim<fNDim; idim++) fDistBuffer[npoints] += TMath::Abs(point[idim] - fData[idim][indexPoints[idx]]); | |
355 | // register candidate neighbor | |
356 | if(fDistBuffer[npoints] < fkNNdist[k-1]){ | |
357 | fIndBuffer[npoints] = indexPoints[idx]; | |
358 | npoints++; | |
359 | } | |
f2040a8f | 360 | } |
316a7f5a | 361 | for(int ibrowse=0; ibrowse<npoints; ibrowse++){ |
362 | if(fDistBuffer[ibrowse] >= fkNNdist[k-1]) continue; | |
363 | //insert neighbor | |
364 | int iNN=0; | |
365 | while(fDistBuffer[ibrowse] >= fkNNdist[iNN]) iNN++; | |
366 | if(debug>=2) printf("\t\tinsert data %d @ %d distance %f\n", fIndBuffer[ibrowse], iNN, fDistBuffer[ibrowse]); | |
367 | ||
368 | itmp = fkNN[iNN]; ftmp = fkNNdist[iNN]; | |
369 | fkNN[iNN] = fIndBuffer[ibrowse]; | |
370 | fkNNdist[iNN] = fDistBuffer[ibrowse]; | |
371 | for(int ireplace=iNN+1; ireplace<k; ireplace++){ | |
372 | jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace]; | |
373 | fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp; | |
374 | itmp = jtmp; ftmp = gtmp; | |
375 | if(ftmp == 9999.) break; | |
f2040a8f | 376 | } |
316a7f5a | 377 | if(debug>=3){ |
378 | for(int i=0; i<k; i++){ | |
379 | if(fkNNdist[i] == 9999.) break; | |
380 | printf("\t\tfkNNdist[%d] = %f\n", i, fkNNdist[i]); | |
381 | } | |
382 | } | |
f2040a8f | 383 | } |
384 | } | |
385 | ||
386 | // register parent | |
117c62ab | 387 | Int_t pn = (tnode-1)>>1;//tnode/2 + (tnode%2) - 1; |
388 | if(pn >= 0 && entry != pn){ | |
f2040a8f | 389 | // check if parent node is eligible at all |
117c62ab | 390 | if(TMath::Abs(point[fAxis[pn]] - fValue[pn]) > fkNNdist[k-1]){ |
f2040a8f | 391 | // mark bound |
392 | nBounds++; | |
393 | // end all recursions | |
394 | if(nBounds==2 * fNDim) break; | |
395 | } | |
396 | currentIndex++; | |
117c62ab | 397 | nodeStack[currentIndex]=pn; |
f2040a8f | 398 | nodeIn[currentIndex]=tnode; |
316a7f5a | 399 | if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]); |
f2040a8f | 400 | } |
117c62ab | 401 | if(IsTerminal(tnode)) continue; |
316a7f5a | 402 | |
f2040a8f | 403 | // register children nodes |
117c62ab | 404 | Int_t cn; |
f2040a8f | 405 | Bool_t kAllow[] = {kTRUE, kTRUE}; |
117c62ab | 406 | ax = fAxis[tnode]; |
407 | val = fValue[tnode]; | |
408 | if(TMath::Abs(point[ax] - val) > fkNNdist[k-1]){ | |
409 | if(point[ax] > val) kAllow[0] = kFALSE; | |
f2040a8f | 410 | else kAllow[1] = kFALSE; |
411 | } | |
412 | for(int ic=1; ic<=2; ic++){ | |
413 | if(!kAllow[ic-1]) continue; | |
117c62ab | 414 | cn = (tnode<<1)+ic; |
415 | if(cn < nAllNodes && entry != cn){ | |
f2040a8f | 416 | currentIndex++; |
117c62ab | 417 | nodeStack[currentIndex] = cn; |
f2040a8f | 418 | nodeIn[currentIndex]=tnode; |
316a7f5a | 419 | if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]); |
f2040a8f | 420 | } |
316a7f5a | 421 | } |
f2040a8f | 422 | } |
423 | // save results | |
424 | in = fkNN; | |
316a7f5a | 425 | d = fkNNdist; |
f2040a8f | 426 | |
117c62ab | 427 | return kTRUE; |
f2040a8f | 428 | } |
429 | ||
430 | ||
431 | ||
432 | //_________________________________________________________________ | |
433 | template <typename Index, typename Value> | |
434 | Index TKDTree<Index, Value>::FindNode(const Value * point){ | |
435 | // | |
436 | // find the terminal node to which point belongs | |
117c62ab | 437 | |
f2040a8f | 438 | Index stackNode[128], inode; |
439 | Int_t currentIndex =0; | |
440 | stackNode[0] = 0; | |
f2040a8f | 441 | while (currentIndex>=0){ |
442 | inode = stackNode[currentIndex]; | |
f2040a8f | 443 | if (IsTerminal(inode)) return inode; |
444 | ||
117c62ab | 445 | currentIndex--; |
446 | if (point[fAxis[inode]]<=fValue[inode]){ | |
f2040a8f | 447 | currentIndex++; |
117c62ab | 448 | stackNode[currentIndex]=(inode<<1)+1; |
f2040a8f | 449 | } |
117c62ab | 450 | if (point[fAxis[inode]]>=fValue[inode]){ |
f2040a8f | 451 | currentIndex++; |
117c62ab | 452 | stackNode[currentIndex]=(inode+1)<<1; |
f2040a8f | 453 | } |
454 | } | |
a9c20b1f | 455 | |
456 | return -1; | |
f2040a8f | 457 | } |
458 | ||
459 | ||
a9c20b1f | 460 | |
f2040a8f | 461 | //_________________________________________________________________ |
462 | template <typename Index, typename Value> | |
463 | void TKDTree<Index, Value>::FindPoint(Value * point, Index &index, Int_t &iter){ | |
464 | // | |
465 | // find the index of point | |
466 | // works only if we keep fData pointers | |
467 | ||
468 | Int_t stackNode[128]; | |
469 | Int_t currentIndex =0; | |
470 | stackNode[0] = 0; | |
471 | iter =0; | |
472 | // | |
473 | while (currentIndex>=0){ | |
474 | iter++; | |
475 | Int_t inode = stackNode[currentIndex]; | |
476 | currentIndex--; | |
477 | if (IsTerminal(inode)){ | |
478 | // investigate terminal node | |
479 | Int_t indexIP = (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset; | |
480 | printf("terminal %d indexP %d\n", inode, indexIP); | |
481 | for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){ | |
482 | Bool_t isOK = kTRUE; | |
483 | indexIP+=ibucket; | |
484 | printf("ibucket %d index %d\n", ibucket, indexIP); | |
485 | if (indexIP>=fNpoints) continue; | |
486 | Int_t index0 = fIndPoints[indexIP]; | |
487 | for (Int_t idim=0;idim<fNDim;idim++) if (fData[idim][index0]!=point[idim]) isOK = kFALSE; | |
488 | if (isOK) index = index0; | |
489 | } | |
490 | continue; | |
491 | } | |
492 | ||
117c62ab | 493 | if (point[fAxis[inode]]<=fValue[inode]){ |
f2040a8f | 494 | currentIndex++; |
495 | stackNode[currentIndex]=(inode*2)+1; | |
496 | } | |
117c62ab | 497 | if (point[fAxis[inode]]>=fValue[inode]){ |
f2040a8f | 498 | currentIndex++; |
499 | stackNode[currentIndex]=(inode*2)+2; | |
500 | } | |
501 | } | |
502 | // | |
503 | // printf("Iter\t%d\n",iter); | |
504 | } | |
505 | ||
506 | //_________________________________________________________________ | |
507 | template <typename Index, typename Value> | |
508 | void TKDTree<Index, Value>::FindInRangeA(Value * point, Value * delta, Index *res , Index &npoints, Index & iter, Int_t bnode) | |
509 | { | |
510 | // | |
511 | // Find all points in the range specified by (point +- range) | |
512 | // res - Resulting indexes are stored in res array | |
513 | // npoints - Number of selected indexes in range | |
514 | // NOTICE: | |
515 | // For some cases it is better to don't keep data - because of memory consumption | |
516 | // If the data are not kept - only boundary conditions are investigated | |
517 | // some of the data can be outside of the range | |
518 | // What is guranteed in this mode: All of the points in the range are selected + some fraction of others (but close) | |
519 | ||
520 | Index stackNode[128]; | |
521 | iter=0; | |
522 | Index currentIndex = 0; | |
523 | stackNode[currentIndex] = bnode; | |
524 | while (currentIndex>=0){ | |
525 | iter++; | |
526 | Int_t inode = stackNode[currentIndex]; | |
527 | // | |
528 | currentIndex--; | |
529 | if (!IsTerminal(inode)){ | |
530 | // not terminal | |
117c62ab | 531 | //TKDNode<Index, Value> * node = &(fNodes[inode]); |
532 | if (point[fAxis[inode]] - delta[fAxis[inode]] < fValue[inode]) { | |
f2040a8f | 533 | currentIndex++; |
534 | stackNode[currentIndex]= (inode*2)+1; | |
535 | } | |
117c62ab | 536 | if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){ |
f2040a8f | 537 | currentIndex++; |
538 | stackNode[currentIndex]= (inode*2)+2; | |
539 | } | |
540 | }else{ | |
541 | Int_t indexIP = | |
542 | (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset; | |
543 | for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){ | |
544 | if (indexIP+ibucket>=fNpoints) break; | |
545 | res[npoints] = fIndPoints[indexIP+ibucket]; | |
546 | npoints++; | |
547 | } | |
548 | } | |
549 | } | |
550 | if (fData){ | |
551 | // | |
552 | // compress rest if data still accesible | |
553 | // | |
554 | Index npoints2 = npoints; | |
555 | npoints=0; | |
556 | for (Index i=0; i<npoints2;i++){ | |
557 | Bool_t isOK = kTRUE; | |
558 | for (Index idim = 0; idim< fNDim; idim++){ | |
559 | if (TMath::Abs(fData[idim][res[i]]- point[idim])>delta[idim]) | |
560 | isOK = kFALSE; | |
561 | } | |
562 | if (isOK){ | |
563 | res[npoints] = res[i]; | |
564 | npoints++; | |
565 | } | |
566 | } | |
567 | } | |
568 | } | |
569 | ||
570 | ||
571 | //_________________________________________________________________ | |
572 | template <typename Index, typename Value> | |
573 | void TKDTree<Index, Value>::FindInRangeB(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode) | |
574 | { | |
575 | // | |
576 | Long64_t goldStatus = (1<<(2*fNDim))-1; // gold status | |
577 | Index stackNode[128]; | |
578 | Long64_t stackStatus[128]; | |
579 | iter=0; | |
580 | Index currentIndex = 0; | |
581 | stackNode[currentIndex] = bnode; | |
582 | stackStatus[currentIndex] = 0; | |
583 | while (currentIndex>=0){ | |
584 | Int_t inode = stackNode[currentIndex]; | |
585 | Long64_t status = stackStatus[currentIndex]; | |
586 | currentIndex--; | |
587 | iter++; | |
588 | if (IsTerminal(inode)){ | |
589 | Int_t indexIP = | |
590 | (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset; | |
591 | for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){ | |
592 | if (indexIP+ibucket>=fNpoints) break; | |
593 | res[npoints] = fIndPoints[indexIP+ibucket]; | |
594 | npoints++; | |
595 | } | |
596 | continue; | |
597 | } | |
598 | // not terminal | |
599 | if (status == goldStatus){ | |
600 | Int_t ileft = inode; | |
601 | Int_t iright = inode; | |
602 | for (;ileft<fNnodes; ileft = (ileft<<1)+1); | |
603 | for (;iright<fNnodes; iright = (iright<<1)+2); | |
604 | Int_t indexL = | |
605 | (ileft >= fCrossNode) ? (ileft-fCrossNode)*fBucketSize : (ileft-fNnodes)*fBucketSize+fOffset; | |
606 | Int_t indexR = | |
607 | (iright >= fCrossNode) ? (iright-fCrossNode)*fBucketSize : (iright-fNnodes)*fBucketSize+fOffset; | |
608 | if (indexL<=indexR){ | |
609 | Int_t endpoint = indexR+fBucketSize; | |
610 | if (endpoint>fNpoints) endpoint=fNpoints; | |
611 | for (Int_t ipoint=indexL;ipoint<endpoint;ipoint++){ | |
612 | res[npoints] = fIndPoints[ipoint]; | |
613 | npoints++; | |
614 | } | |
615 | continue; | |
616 | } | |
617 | } | |
618 | // | |
117c62ab | 619 | // TKDNode<Index, Value> * node = &(fNodes[inode]); |
620 | if (point[fAxis[inode]] - delta[fAxis[inode]] < fValue[inode]) { | |
f2040a8f | 621 | currentIndex++; |
622 | stackNode[currentIndex]= (inode<<1)+1; | |
117c62ab | 623 | if (point[fAxis[inode]] + delta[fAxis[inode]] > fValue[inode]) |
624 | stackStatus[currentIndex]= status | (1<<(2*fAxis[inode])); | |
f2040a8f | 625 | } |
117c62ab | 626 | if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){ |
f2040a8f | 627 | currentIndex++; |
628 | stackNode[currentIndex]= (inode<<1)+2; | |
117c62ab | 629 | if (point[fAxis[inode]] - delta[fAxis[inode]]<fValue[inode]) |
630 | stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]+1)); | |
f2040a8f | 631 | } |
632 | } | |
633 | if (fData){ | |
634 | // compress rest | |
635 | Index npoints2 = npoints; | |
636 | npoints=0; | |
637 | for (Index i=0; i<npoints2;i++){ | |
638 | Bool_t isOK = kTRUE; | |
639 | for (Index idim = 0; idim< fNDim; idim++){ | |
640 | if (TMath::Abs(fData[idim][res[i]]- point[idim])>delta[idim]) | |
641 | isOK = kFALSE; | |
642 | } | |
643 | if (isOK){ | |
644 | res[npoints] = res[i]; | |
645 | npoints++; | |
646 | } | |
647 | } | |
648 | } | |
649 | } | |
650 | ||
651 | ||
652 | ||
653 | //_________________________________________________________________ | |
654 | template <typename Index, typename Value> | |
655 | void TKDTree<Index, Value>::SetData(Index npoints, Index ndim, UInt_t bsize, Value **data) | |
656 | { | |
657 | // TO DO | |
658 | // | |
659 | // Check reconstruction/reallocation of memory of data. Maybe it is not | |
660 | // necessary to delete and realocate space but only to use the same space | |
661 | ||
662 | Clear(); | |
663 | ||
664 | //Columnwise!!!! | |
665 | fData = data; | |
666 | fNpoints = npoints; | |
667 | fNDim = ndim; | |
668 | fBucketSize = bsize; | |
669 | ||
670 | Build(); | |
671 | } | |
672 | ||
673 | ||
674 | ||
675 | //_________________________________________________________________ | |
676 | template <typename Index, typename Value> | |
677 | void TKDTree<Index, Value>::Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) | |
678 | { | |
679 | //Value min, max; | |
680 | Index i; | |
681 | min = a[index[0]]; | |
682 | max = a[index[0]]; | |
683 | for (i=0; i<ntotal; i++){ | |
684 | if (a[index[i]]<min) min = a[index[i]]; | |
685 | if (a[index[i]]>max) max = a[index[i]]; | |
686 | } | |
687 | } | |
688 | ||
689 | ||
690 | //_________________________________________________________________ | |
691 | template <typename Index, typename Value> | |
692 | Value TKDTree<Index, Value>::KOrdStat(Index ntotal, Value *a, Index k, Index *index) | |
693 | { | |
694 | // | |
695 | //copy of the TMath::KOrdStat because I need an Index work array | |
696 | ||
697 | Index i, ir, j, l, mid; | |
698 | Index arr; | |
699 | Index temp; | |
700 | ||
701 | Index rk = k; | |
702 | l=0; | |
703 | ir = ntotal-1; | |
704 | for(;;) { | |
705 | if (ir<=l+1) { //active partition contains 1 or 2 elements | |
706 | if (ir == l+1 && a[index[ir]]<a[index[l]]) | |
707 | {temp = index[l]; index[l]=index[ir]; index[ir]=temp;} | |
708 | Value tmp = a[index[rk]]; | |
709 | return tmp; | |
710 | } else { | |
711 | mid = (l+ir) >> 1; //choose median of left, center and right | |
712 | {temp = index[mid]; index[mid]=index[l+1]; index[l+1]=temp;}//elements as partitioning element arr. | |
713 | if (a[index[l]]>a[index[ir]]) //also rearrange so that a[l]<=a[l+1] | |
714 | {temp = index[l]; index[l]=index[ir]; index[ir]=temp;} | |
715 | ||
716 | if (a[index[l+1]]>a[index[ir]]) | |
717 | {temp=index[l+1]; index[l+1]=index[ir]; index[ir]=temp;} | |
718 | ||
719 | if (a[index[l]]>a[index[l+1]]) | |
720 | {temp = index[l]; index[l]=index[l+1]; index[l+1]=temp;} | |
721 | ||
722 | i=l+1; //initialize pointers for partitioning | |
723 | j=ir; | |
724 | arr = index[l+1]; | |
725 | for (;;) { | |
726 | do i++; while (a[index[i]]<a[arr]); | |
727 | do j--; while (a[index[j]]>a[arr]); | |
728 | if (j<i) break; //pointers crossed, partitioning complete | |
729 | {temp=index[i]; index[i]=index[j]; index[j]=temp;} | |
730 | } | |
731 | index[l+1]=index[j]; | |
732 | index[j]=arr; | |
733 | if (j>=rk) ir = j-1; //keep active the partition that | |
734 | if (j<=rk) l=i; //contains the k_th element | |
735 | } | |
736 | } | |
737 | } | |
738 | ||
739 | //_________________________________________________________________ | |
740 | template <typename Index, typename Value> | |
741 | void TKDTree<Index, Value>::MakeBoundaries(Value *range) | |
742 | { | |
743 | // Build boundaries for each node. | |
744 | ||
745 | ||
117c62ab | 746 | if(range) memcpy(fRange, range, fNDimm*sizeof(Value)); |
f2040a8f | 747 | // total number of nodes including terminal nodes |
a3408ed3 | 748 | Int_t totNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0); |
117c62ab | 749 | fBoundaries = new Value[totNodes*fNDimm]; |
750 | //Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes)); | |
751 | ||
316a7f5a | 752 | |
f2040a8f | 753 | // loop |
117c62ab | 754 | Value *tbounds = 0x0, *cbounds = 0x0; |
755 | Int_t cn; | |
f2040a8f | 756 | for(int inode=fNnodes-1; inode>=0; inode--){ |
117c62ab | 757 | tbounds = &fBoundaries[inode*fNDimm]; |
758 | memcpy(tbounds, fRange, fNDimm*sizeof(Value)); | |
f2040a8f | 759 | |
760 | // check left child node | |
117c62ab | 761 | cn = (inode<<1)+1; |
762 | if(IsTerminal(cn)) CookBoundaries(inode, kTRUE); | |
763 | cbounds = &fBoundaries[fNDimm*cn]; | |
764 | for(int idim=0; idim<fNDim; idim++) tbounds[idim<<1] = cbounds[idim<<1]; | |
f2040a8f | 765 | |
766 | // check right child node | |
117c62ab | 767 | cn = (inode+1)<<1; |
768 | if(IsTerminal(cn)) CookBoundaries(inode, kFALSE); | |
769 | cbounds = &fBoundaries[fNDimm*cn]; | |
770 | for(int idim=0; idim<fNDim; idim++) tbounds[(idim<<1)+1] = cbounds[(idim<<1)+1]; | |
f2040a8f | 771 | } |
772 | } | |
773 | ||
774 | //_________________________________________________________________ | |
775 | template <typename Index, typename Value> | |
117c62ab | 776 | void TKDTree<Index, Value>::CookBoundaries(const Int_t node, Bool_t LEFT) |
f2040a8f | 777 | { |
778 | // define index of this terminal node | |
117c62ab | 779 | Int_t index = (node<<1) + (LEFT ? 1 : 2); |
780 | //Info("CookBoundaries()", Form("Node %d", index)); | |
f2040a8f | 781 | |
782 | // build and initialize boundaries for this node | |
117c62ab | 783 | Value *tbounds = &fBoundaries[fNDimm*index]; |
784 | memcpy(tbounds, fRange, fNDimm*sizeof(Value)); | |
f2040a8f | 785 | Bool_t flag[256]; // cope with up to 128 dimensions |
117c62ab | 786 | memset(flag, kFALSE, fNDimm); |
f2040a8f | 787 | Int_t nvals = 0; |
788 | ||
789 | // recurse parent nodes | |
117c62ab | 790 | Int_t pn = node; |
791 | while(pn >= 0 && nvals < fNDimm){ | |
f2040a8f | 792 | if(LEFT){ |
117c62ab | 793 | index = (fAxis[pn]<<1)+1; |
794 | if(!flag[index]) { | |
795 | tbounds[index] = fValue[pn]; | |
796 | flag[index] = kTRUE; | |
f2040a8f | 797 | nvals++; |
798 | } | |
799 | } else { | |
117c62ab | 800 | index = fAxis[pn]<<1; |
801 | if(!flag[index]) { | |
802 | tbounds[index] = fValue[pn]; | |
803 | flag[index] = kTRUE; | |
f2040a8f | 804 | nvals++; |
805 | } | |
806 | } | |
117c62ab | 807 | LEFT = pn&1; |
808 | pn = (pn - 1)>>1; | |
f2040a8f | 809 | } |
810 | } | |
811 | ||
f2040a8f | 812 | template class TKDTree<Int_t, Float_t>; |
813 | template class TKDTree<Int_t, Double_t>; | |
814 |