MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | List of all members
KDtree< DIM > Struct Template Reference

#include <kdtree.h>

Inheritance diagram for KDtree< DIM >:
[legend]
Collaboration diagram for KDtree< DIM >:
[legend]

Public Member Functions

 KDtree (std::vector< Point< DIM > > &pts)
 
 ~KDtree ()
 
MFloat disti (MInt jpt, MInt kpt)
 
MInt locate (Point< DIM > pt)
 
MInt locate (MInt jpt)
 
MInt nearest (Point< DIM > pt, MFloat &dist)
 
void nnearest (MInt jpt, MInt *nn, MFloat *dn, MInt n)
 
MInt locatenear (Point< DIM > pt, MFloat r, MInt *list, MInt nmax, MBool returnCellId=true)
 
MInt locatenearest (Point< DIM > pt, MFloat r, MInt *list, MFloat *dn, MInt nmax, MBool &overflow)
 

Static Public Member Functions

static void sift_down (MFloat *heap, MInt *ndx, MInt nn)
 

Public Attributes

MLong nboxes
 
std::vector< Point< DIM > > ptss
 
Boxnode< DIM > * boxes
 
MLong npts
 
std::vector< MInt > ptindx
 
std::vector< MInt > rptindx
 

Static Public Attributes

static const MFloat BIG
 

Detailed Description

template<MInt DIM>
struct KDtree< DIM >

Definition at line 73 of file kdtree.h.

Constructor & Destructor Documentation

◆ KDtree()

template<MInt DIM>
KDtree< DIM >::KDtree ( std::vector< Point< DIM > > &  pts)

Definition at line 101 of file kdtree.h.

102 : ptss(pts), boxes(nullptr), npts(pts.size()), ptindx(npts), rptindx(npts) {
103 MInt ntmp, m, k, kk, j, nowtask, jbox, np, tmom, tdim, ptlo, pthi;
104 MInt* hp;
105 MFloat* cp;
106 MInt taskmom[50], taskdim[50];
107 if(npts == 0) return;
108 for(k = 0; k < npts; k++)
109 ptindx[k] = k;
110 m = 1;
111 for(ntmp = npts; ntmp; ntmp >>= 1) {
112 m <<= 1;
113 }
114 nboxes = 2 * npts - (m >> 1);
115 if(m < nboxes) nboxes = m;
116 nboxes--;
117 if(nboxes > 0) {
118 boxes = new Boxnode<DIM>[nboxes];
119 MFloat* coords = new MFloat[DIM * npts];
120 for(j = 0, kk = 0; j < DIM; j++, kk += npts) {
121 for(k = 0; k < npts; k++)
122 coords[kk + k] = pts[k].x[j];
123 }
124 Point<DIM> lo(-BIG, -BIG, -BIG), hi(BIG, BIG, BIG);
125 boxes[0] = Boxnode<DIM>(lo, hi, 0, 0, 0, 0, npts - 1);
126 jbox = 0;
127 taskmom[1] = 0;
128 taskdim[1] = 0;
129 nowtask = 1;
130 while(nowtask) {
131 tmom = taskmom[nowtask];
132 tdim = taskdim[nowtask--];
133 ptlo = boxes[tmom].ptlo;
134 pthi = boxes[tmom].pthi;
135 hp = &ptindx[ptlo];
136 cp = &coords[tdim * npts];
137 np = pthi - ptlo + 1;
138 kk = (np - 1) / 2;
139 (void)selecti<DIM>(kk, hp, np, cp);
140 hi = boxes[tmom].hi;
141 lo = boxes[tmom].lo;
142 hi.x[tdim] = lo.x[tdim] = coords[tdim * npts + hp[kk]];
143 if(jbox < nboxes - 1) boxes[++jbox] = Boxnode<DIM>(boxes[tmom].lo, hi, tmom, 0, 0, ptlo, ptlo + kk);
144 if(jbox < nboxes - 1) boxes[++jbox] = Boxnode<DIM>(lo, boxes[tmom].hi, tmom, 0, 0, ptlo + kk + 1, pthi);
145 if(tmom < nboxes) boxes[tmom].dau1 = jbox - 1;
146 if(tmom < nboxes) boxes[tmom].dau2 = jbox;
147 if(kk > 1) {
148 taskmom[++nowtask] = jbox - 1;
149 taskdim[nowtask] = (tdim + 1) % DIM;
150 }
151 if(np - kk > 3) {
152 taskmom[++nowtask] = jbox;
153 taskdim[nowtask] = (tdim + 1) % DIM;
154 }
155 }
156 for(j = 0; j < npts; j++)
157 rptindx[ptindx[j]] = j;
158 delete[] coords;
159 }
160}
MInt
int32_t MInt
Definition: maiatypes.h:62
MFloat
double MFloat
Definition: maiatypes.h:52
Boxnode
Definition: kdtree.h:28
KDtree::boxes
Boxnode< DIM > * boxes
Definition: kdtree.h:78
KDtree::ptss
std::vector< Point< DIM > > ptss
Definition: kdtree.h:77
KDtree::npts
MLong npts
Definition: kdtree.h:79
KDtree::nboxes
MLong nboxes
Definition: kdtree.h:75
KDtree::ptindx
std::vector< MInt > ptindx
Definition: kdtree.h:81
KDtree::BIG
static const MFloat BIG
Definition: kdtree.h:74
KDtree::rptindx
std::vector< MInt > rptindx
Definition: kdtree.h:81
Point
Definition: pointbox.h:20

◆ ~KDtree()

template<MInt DIM>
KDtree< DIM >::~KDtree ( )
inline

Definition at line 84 of file kdtree.h.

84 {
85 if(boxes != nullptr) delete[] boxes;
86 ptss.clear();
87 }

Member Function Documentation

◆ disti()

template<MInt DIM>
MFloat KDtree< DIM >::disti ( MInt  jpt,
MInt  kpt 
)

Definition at line 162 of file kdtree.h.

162 {
163 if(jpt == kpt)
164 return BIG;
165 else
166 return dist(ptss[jpt], ptss[kpt]);
167}
dist
MFloat dist(const Point< DIM > &p, const Point< DIM > &q)
Definition: pointbox.h:54

◆ locate() [1/2]

template<MInt DIM>
MInt KDtree< DIM >::locate ( MInt  jpt)

Definition at line 185 of file kdtree.h.

185 {
186 MInt nb, d1, jh;
187 jh = rptindx[jpt];
188 nb = 0;
189 while(boxes[nb].dau1) {
190 d1 = boxes[nb].dau1;
191 if(jh <= boxes[d1].pthi)
192 nb = d1;
193 else
194 nb = boxes[nb].dau2;
195 }
196 return nb;
197}

◆ locate() [2/2]

template<MInt DIM>
MInt KDtree< DIM >::locate ( Point< DIM >  pt)

Definition at line 169 of file kdtree.h.

169 {
170 MInt nb, d1, jdim;
171 nb = jdim = 0;
172 while(boxes[nb].dau1) {
173 d1 = boxes[nb].dau1;
174 if(pt.x[jdim] <= boxes[d1].hi.x[jdim])
175 nb = d1;
176 else
177 nb = boxes[nb].dau2;
178 // jdim = ++jdim % DIM;
179 ++jdim;
180 jdim %= DIM;
181 }
182 return nb;
183}
Point::x
MFloat x[DIM]
Definition: pointbox.h:21

◆ locatenear()

template<MInt DIM>
MInt KDtree< DIM >::locatenear ( Point< DIM >  pt,
MFloat  r,
MInt *  list,
MInt  nmax,
MBool  returnCellId = true 
)

Definition at line 320 of file kdtree.h.

320 {
321 MInt k, i, nb, nbold, nret, ntask, jdim, d1, d2;
322 MInt task[50];
323 nb = jdim = nret = 0;
324 ASSERT(!(r < 0.0), "radius must be nonnegative");
325 ASSERT(nmax > 0, "");
326 if(npts == 0) {
327 return 0;
328 } else if(npts == 1) {
329 if(dist(ptss[ptindx[0]], pt) < r && nmax > 0) list[nret++] = ptss[ptindx[0]].cellId;
330 } else if(npts == 2) {
331 if(dist(ptss[ptindx[0]], pt) < r && nmax > 0) list[nret++] = ptss[ptindx[0]].cellId;
332 if(dist(ptss[ptindx[1]], pt) < r && nmax > 1) list[nret++] = ptss[ptindx[1]].cellId;
333 } else {
334 while(boxes[nb].dau1) {
335 nbold = nb;
336 d1 = boxes[nb].dau1;
337 d2 = boxes[nb].dau2;
338 if(pt.x[jdim] + r <= boxes[d1].hi.x[jdim])
339 nb = d1;
340 else if(pt.x[jdim] - r >= boxes[d2].lo.x[jdim])
341 nb = d2;
342 // jdim = ++jdim % DIM;//undefined behavior
343 jdim = (jdim + 1) % DIM;
344 if(nb == nbold) break;
345 }
346 task[1] = nb;
347 ntask = 1;
348 while(ntask) {
349 k = task[ntask--];
350 if(dist(boxes[k], pt) > r) continue;
351 if(boxes[k].dau1) {
352 task[++ntask] = boxes[k].dau1;
353 task[++ntask] = boxes[k].dau2;
354 } else {
355 for(i = boxes[k].ptlo; i <= boxes[k].pthi; i++) {
356 if(dist(ptss[ptindx[i]], pt) <= r && nret < nmax) {
357 if(returnCellId) {
358 list[nret++] = ptss[ptindx[i]].cellId;
359 } else {
360 list[nret++] = ptindx[i];
361 }
362 }
363
364 if(nret == nmax) return nmax;
365 }
366 }
367 }
368 }
369 return nret;
370}

◆ locatenearest()

template<MInt DIM>
MInt KDtree< DIM >::locatenearest ( Point< DIM >  pt,
MFloat  r,
MInt *  list,
MFloat *  dn,
MInt  nmax,
MBool &  overflow 
)

Definition at line 372 of file kdtree.h.

372 {
373 MInt k, i, nb, nbold, nret, ntask, jdim, d1, d2;
374 MInt task[50];
375 nb = jdim = nret = 0;
376 overflow = false;
377 ASSERT(!(r < 0.0), "radius must be nonnegative");
378 ASSERT(nmax > 0, "");
379 if(npts == 0) {
380 return 0;
381 } else if(npts == 1) {
382 if(dist(ptss[ptindx[0]], pt) < r && nmax > 0) list[nret++] = ptss[ptindx[0]].cellId;
383 } else if(npts == 2) {
384 MInt smaller = (dist(ptss[ptindx[0]], pt) < dist(ptss[ptindx[1]], pt)) ? 0 : 1;
385 MInt other = (smaller + 1) % 2;
386 if(dist(ptss[ptindx[smaller]], pt) < r && nmax > 0) list[nret++] = ptss[ptindx[smaller]].cellId;
387 if(dist(ptss[ptindx[other]], pt) < r && nmax > 1) list[nret++] = ptss[ptindx[other]].cellId;
388 } else {
389 while(boxes[nb].dau1) {
390 nbold = nb;
391 d1 = boxes[nb].dau1;
392 d2 = boxes[nb].dau2;
393 if(pt.x[jdim] + r <= boxes[d1].hi.x[jdim])
394 nb = d1;
395 else if(pt.x[jdim] - r >= boxes[d2].lo.x[jdim])
396 nb = d2;
397 // jdim = ++jdim % DIM;//undefined behavior
398 jdim = (jdim + 1) % DIM;
399 if(nb == nbold) break;
400 }
401 task[1] = nb;
402 ntask = 1;
403 while(ntask) {
404 k = task[ntask--];
405 if(dist(boxes[k], pt) > r) continue;
406 if(boxes[k].dau1) {
407 task[++ntask] = boxes[k].dau1;
408 task[++ntask] = boxes[k].dau2;
409 } else {
410 for(i = boxes[k].ptlo; i <= boxes[k].pthi; i++) {
411 MFloat dst = dist(ptss[ptindx[i]], pt);
412 if(dst < r) {
413 if(nret == nmax) overflow = true;
414 MInt s = 0;
415 for(s = 0; s < nret; s++) {
416 if(dst < dn[s]) {
417 for(MInt t = std::min(nret, nmax - 1); t > s; t--) {
418 list[t] = list[t - 1];
419 dn[t] = dn[t - 1];
420 }
421 break;
422 }
423 }
424 if(s < nmax) {
425 list[s] = ptss[ptindx[i]].cellId;
426 dn[s] = dst;
427 }
428 if(nret < nmax) nret++;
429 }
430 }
431 }
432 }
433 }
434 return nret;
435}

◆ nearest()

template<MInt DIM>
MInt KDtree< DIM >::nearest ( Point< DIM >  pt,
MFloat &  dist 
)

Definition at line 199 of file kdtree.h.

199 {
200 MInt i, k, ntask;
201 MInt nrst = 0;
202 MInt task[50];
203 MFloat dnrst = BIG, d;
204 k = locate(pt);
205 for(i = boxes[k].ptlo; i <= boxes[k].pthi; i++) {
206 d = dist(ptss[ptindx[i]], pt);
207 if(d < dnrst) {
208 nrst = ptindx[i];
209 dnrst = d;
210 }
211 }
212 task[1] = 0;
213 ntask = 1;
214 while(ntask) {
215 k = task[ntask--];
216 if(dist(boxes[k], pt) < dnrst) {
217 if(boxes[k].dau1) {
218 task[++ntask] = boxes[k].dau1;
219 task[++ntask] = boxes[k].dau2;
220 } else {
221 for(i = boxes[k].ptlo; i <= boxes[k].pthi; i++) {
222 d = dist(ptss[ptindx[i]], pt);
223 if(d < dnrst) {
224 nrst = ptindx[i];
225 dnrst = d;
226 }
227 }
228 }
229 }
230 }
231 dist1 = dist(ptss[nrst], pt); // gives the distance!
232 return ptss[nrst].cellId; // nrst;
233}
KDtree::locate
MInt locate(Point< DIM > pt)
Definition: kdtree.h:169

◆ nnearest()

template<MInt DIM>
void KDtree< DIM >::nnearest ( MInt  jpt,
MInt *  nn,
MFloat *  dn,
MInt  n 
)

Definition at line 235 of file kdtree.h.

235 {
236 MInt i, k, ntask, kp;
237 MInt task[50];
238 MFloat d;
239 ASSERT(!(n > npts - 1), "too many neighbors requested");
240 for(i = 0; i < n; i++)
241 nn[i] = -1;
242 for(i = 0; i < n; i++)
243 dn[i] = BIG;
244 if(npts == 0) {
245 return;
246 }
247 if(npts == 1) {
248 nn[0] = ptss[ptindx[0]].cellId;
249 dn[0] = disti(ptindx[0], jpt);
250 return;
251 } else if(npts == 2) {
252 MInt smaller = (disti(ptindx[0], jpt) < disti(ptindx[1], jpt)) ? 0 : 1;
253 MInt other = (smaller + 1) % 2;
254 nn[0] = ptss[ptindx[smaller]].cellId;
255 dn[0] = disti(ptindx[smaller], jpt);
256 if(n > 1) {
257 nn[1] = ptss[ptindx[other]].cellId;
258 dn[1] = disti(ptindx[other], jpt);
259 }
260 return;
261 }
262 kp = boxes[locate(jpt)].mom;
263 while(boxes[kp].pthi - boxes[kp].ptlo < n)
264 kp = boxes[kp].mom;
265 for(i = boxes[kp].ptlo; i <= boxes[kp].pthi; i++) {
266 if(jpt == ptindx[i]) continue;
267 d = disti(ptindx[i], jpt);
268 if(d < dn[0]) {
269 dn[0] = d;
270 nn[0] = ptindx[i];
271 if(n > 1) sift_down(dn, nn, n);
272 }
273 }
274 task[1] = 0;
275 ntask = 1;
276 while(ntask) {
277 k = task[ntask--];
278 if(k == kp) continue;
279 if(dist(boxes[k], ptss[jpt]) < dn[0]) {
280 if(boxes[k].dau1) {
281 task[++ntask] = boxes[k].dau1;
282 task[++ntask] = boxes[k].dau2;
283 } else {
284 for(i = boxes[k].ptlo; i <= boxes[k].pthi; i++) {
285 d = disti(ptindx[i], jpt);
286 if(d < dn[0]) {
287 dn[0] = d;
288 nn[0] = ptindx[i];
289 if(n > 1) sift_down(dn, nn, n);
290 }
291 }
292 }
293 }
294 }
295 for(i = 0; i < n; i++)
296 nn[i] = (nn[i] > -1) ? ptss[nn[i]].cellId : -1;
297 return;
298}
KDtree::sift_down
static void sift_down(MFloat *heap, MInt *ndx, MInt nn)
Definition: kdtree.h:300
KDtree::disti
MFloat disti(MInt jpt, MInt kpt)
Definition: kdtree.h:162

◆ sift_down()

template<MInt DIM>
void KDtree< DIM >::sift_down ( MFloat *  heap,
MInt *  ndx,
MInt  nn 
)
static

Definition at line 300 of file kdtree.h.

300 {
301 MInt n = nn - 1;
302 MInt j, jold, ia;
303 MFloat a;
304 a = heap[0];
305 ia = ndx[0];
306 jold = 0;
307 j = 1;
308 while(j <= n) {
309 if(j < n && heap[j] < heap[j + 1]) j++;
310 if(a >= heap[j]) break;
311 heap[jold] = heap[j];
312 ndx[jold] = ndx[j];
313 jold = j;
314 j = 2 * j + 1;
315 }
316 heap[jold] = a;
317 ndx[jold] = ia;
318}
a
Definition: contexttypes.h:19

Member Data Documentation

◆ BIG

template<MInt DIM>
const MFloat KDtree< DIM >::BIG
static

Definition at line 74 of file kdtree.h.

◆ boxes

template<MInt DIM>
Boxnode<DIM>* KDtree< DIM >::boxes

Definition at line 78 of file kdtree.h.

◆ nboxes

template<MInt DIM>
MLong KDtree< DIM >::nboxes

Definition at line 75 of file kdtree.h.

◆ npts

template<MInt DIM>
MLong KDtree< DIM >::npts

Definition at line 79 of file kdtree.h.

◆ ptindx

template<MInt DIM>
std::vector<MInt> KDtree< DIM >::ptindx

Definition at line 81 of file kdtree.h.

◆ ptss

template<MInt DIM>
std::vector<Point<DIM> > KDtree< DIM >::ptss

Definition at line 77 of file kdtree.h.

◆ rptindx

template<MInt DIM>
std::vector<MInt> KDtree< DIM >::rptindx

Definition at line 81 of file kdtree.h.

The documentation for this struct was generated from the following file: