50 template<
class NodeData,
class Real>
int OctNode<NodeData,Real>::UseAlloc=0;
53 template<
class NodeData,
class Real>
59 internalAllocator.set(blockSize);
63 template<
class NodeData,
class Real>
66 template <
class NodeData,
class Real>
69 d=off[0]=off[1]=off[2]=0;
72 template <
class NodeData,
class Real>
74 if(!UseAlloc){
if(children){
delete[] children;}}
77 template <
class NodeData,
class Real>
81 if( !children ) initChildren();
82 for(
int i=0 ; i<8 ; i++ ) children[i].setFullDepth( maxDepth-1 );
86 template <
class NodeData,
class Real>
90 if(UseAlloc){children=internalAllocator.newElements(8);}
92 if(children){
delete[] children;}
97 fprintf(stderr,
"Failed to initialize children in OctNode::initChildren\n");
102 depthAndOffset(d,off);
107 children[idx].parent=
this;
108 children[idx].children=NULL;
110 off2[0]=(off[0]<<1)+i;
111 off2[1]=(off[1]<<1)+j;
112 off2[2]=(off[2]<<1)+k;
113 Index(d+1,off2,children[idx].d,children[idx].off);
119 template <
class NodeData,
class Real>
122 off[0]=short((1<<depth)+offset[0]-1);
123 off[1]=short((1<<depth)+offset[1]-1);
124 off[2]=short((1<<depth)+offset[2]-1);
127 template<
class NodeData,
class Real>
130 offset[0]=(int(off[0])+1)&(~(1<<depth));
131 offset[1]=(int(off[1])+1)&(~(1<<depth));
132 offset[2]=(int(off[2])+1)&(~(1<<depth));
134 template<
class NodeData,
class Real>
136 template<
class NodeData,
class Real>
138 depth=int(index&DepthMask);
139 offset[0]=(int((index>>OffsetShift1)&OffsetMask)+1)&(~(1<<depth));
140 offset[1]=(int((index>>OffsetShift2)&OffsetMask)+1)&(~(1<<depth));
141 offset[2]=(int((index>>OffsetShift3)&OffsetMask)+1)&(~(1<<depth));
143 template<
class NodeData,
class Real>
145 template <
class NodeData,
class Real>
149 offset[0]=(int(off[0])+1)&(~(1<<depth));
150 offset[1]=(int(off[1])+1)&(~(1<<depth));
151 offset[2]=(int(off[2])+1)&(~(1<<depth));
152 width=
Real(1.0/(1<<depth));
153 for(
int dim=0;dim<DIMENSION;dim++){center.
coords[dim]=
Real(0.5+offset[dim])*width;}
155 template<
class NodeData ,
class Real >
160 centerAndWidth( c , w );
162 return (c[0]-w)<p[0] && p[0]<=(c[0]+w) && (c[1]-w)<p[1] && p[1]<=(c[1]+w) && (c[2]-w)<p[2] && p[2]<=(c[2]+w);
164 template <
class NodeData,
class Real>
167 depth=index&DepthMask;
168 offset[0]=(int((index>>OffsetShift1)&OffsetMask)+1)&(~(1<<depth));
169 offset[1]=(int((index>>OffsetShift2)&OffsetMask)+1)&(~(1<<depth));
170 offset[2]=(int((index>>OffsetShift3)&OffsetMask)+1)&(~(1<<depth));
171 width=
Real(1.0/(1<<depth));
172 for(
int dim=0;dim<DIMENSION;dim++){center.
coords[dim]=
Real(0.5+offset[dim])*width;}
175 template <
class NodeData,
class Real>
177 if(!children){
return 0;}
181 d=children[i].maxDepth();
187 template <
class NodeData,
class Real>
189 if(!children){
return 1;}
196 template <
class NodeData,
class Real>
198 if(!children){
return 1;}
205 template<
class NodeData,
class Real>
207 if(depth()>maxDepth){
return 0;}
208 if(!children){
return 1;}
211 for(
int i=0;i<
Cube::CORNERS;i++){c+=children[i].maxDepthLeaves(maxDepth);}
215 template <
class NodeData,
class Real>
223 template <
class NodeData,
class Real>
226 if( !current->
parent || current==
this )
return NULL;
228 else return current+1;
230 template <
class NodeData,
class Real>
232 if(!current->
parent || current==
this){
return NULL;}
234 else{
return current+1;}
236 template<
class NodeData ,
class Real >
239 if( !current->
parent || current==
this )
return NULL;
241 else return current-1;
243 template<
class NodeData ,
class Real >
246 if( !current->
parent || current==
this )
return NULL;
248 else return current-1;
250 template <
class NodeData,
class Real>
258 const OctNode* temp=nextBranch(current);
259 if(!temp){
return NULL;}
262 template <
class NodeData,
class Real>
270 OctNode* temp=nextBranch(current);
271 if(!temp){
return NULL;}
275 template <
class NodeData,
class Real>
278 if( !current )
return this;
280 else return nextBranch(current);
282 template <
class NodeData,
class Real>
285 if( !current )
return this;
287 else return nextBranch( current );
290 template <
class NodeData,
class Real>
294 centerAndWidth(center,width);
295 for(
int dim=0;dim<DIMENSION;dim++){
296 printf(
"%[%f,%f]",center.
coords[dim]-width/2,center.
coords[dim]+width/2);
297 if(dim<DIMENSION-1){printf(
"x");}
302 template <
class NodeData,
class Real>
305 template <
class NodeData,
class Real>
306 template<
class NodeAdjacencyFunction>
308 if(processCurrent){F->Function(
this,node);}
309 if(children){__processNodeNodes(node,F);}
311 template <
class NodeData,
class Real>
312 template<
class NodeAdjacencyFunction>
314 if(processCurrent){F->Function(
this,node);}
318 __processNodeFaces(node,F,c1,c2,c3,c4);
321 template <
class NodeData,
class Real>
322 template<
class NodeAdjacencyFunction>
324 if(processCurrent){F->Function(
this,node);}
328 __processNodeEdges(node,F,c1,c2);
331 template <
class NodeData,
class Real>
332 template<
class NodeAdjacencyFunction>
334 if(processCurrent){F->Function(
this,node);}
338 F->Function(temp,node);
341 template <
class NodeData,
class Real>
342 template<
class NodeAdjacencyFunction>
344 F->Function(&children[0],node);
345 F->Function(&children[1],node);
346 F->Function(&children[2],node);
347 F->Function(&children[3],node);
348 F->Function(&children[4],node);
349 F->Function(&children[5],node);
350 F->Function(&children[6],node);
351 F->Function(&children[7],node);
352 if(children[0].children){children[0].__processNodeNodes(node,F);}
353 if(children[1].children){children[1].__processNodeNodes(node,F);}
354 if(children[2].children){children[2].__processNodeNodes(node,F);}
355 if(children[3].children){children[3].__processNodeNodes(node,F);}
356 if(children[4].children){children[4].__processNodeNodes(node,F);}
357 if(children[5].children){children[5].__processNodeNodes(node,F);}
358 if(children[6].children){children[6].__processNodeNodes(node,F);}
359 if(children[7].children){children[7].__processNodeNodes(node,F);}
361 template <
class NodeData,
class Real>
362 template<
class NodeAdjacencyFunction>
364 F->Function(&children[cIndex1],node);
365 F->Function(&children[cIndex2],node);
366 if(children[cIndex1].children){children[cIndex1].__processNodeEdges(node,F,cIndex1,cIndex2);}
367 if(children[cIndex2].children){children[cIndex2].__processNodeEdges(node,F,cIndex1,cIndex2);}
369 template <
class NodeData,
class Real>
370 template<
class NodeAdjacencyFunction>
372 F->Function(&children[cIndex1],node);
373 F->Function(&children[cIndex2],node);
374 F->Function(&children[cIndex3],node);
375 F->Function(&children[cIndex4],node);
376 if(children[cIndex1].children){children[cIndex1].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
377 if(children[cIndex2].children){children[cIndex2].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
378 if(children[cIndex3].children){children[cIndex3].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
379 if(children[cIndex4].children){children[cIndex4].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
381 template<
class NodeData,
class Real>
382 template<
class NodeAdjacencyFunction>
384 int c1[3],c2[3],w1,w2;
387 w1=node1->
width(maxDepth+1);
388 w2=node2->
width(maxDepth+1);
390 ProcessNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,F,processCurrent);
392 template<
class NodeData,
class Real>
393 template<
class NodeAdjacencyFunction>
397 NodeAdjacencyFunction* F,
int processCurrent){
398 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
399 if(processCurrent){F->Function(node2,node1);}
401 __ProcessNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,F);
403 template<
class NodeData,
class Real>
404 template<
class TerminatingNodeAdjacencyFunction>
406 int c1[3],c2[3],w1,w2;
409 w1=node1->
width(maxDepth+1);
410 w2=node2->
width(maxDepth+1);
412 ProcessTerminatingNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,F,processCurrent);
414 template<
class NodeData,
class Real>
415 template<
class TerminatingNodeAdjacencyFunction>
419 TerminatingNodeAdjacencyFunction* F,
int processCurrent)
421 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
422 if(processCurrent){F->Function(node2,node1);}
424 __ProcessTerminatingNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,F);
426 template<
class NodeData,
class Real>
427 template<
class Po
intAdjacencyFunction>
432 w2 = node2->
width( maxDepth+1 );
433 ProcessPointAdjacentNodes( c1[0]-c2[0] , c1[1]-c2[1] , c1[2]-c2[2] , node2 , (width2*w2)>>1 , w2 , F , processCurrent );
435 template<
class NodeData,
class Real>
436 template<
class Po
intAdjacencyFunction>
439 PointAdjacencyFunction* F,
int processCurrent)
441 if( !Overlap(dx,dy,dz,radius2) )
return;
442 if( processCurrent ) F->Function(node2);
444 __ProcessPointAdjacentNodes( -dx , -dy , -dz , node2 , radius2 , width2>>1 , F );
446 template<
class NodeData,
class Real>
447 template<
class NodeAdjacencyFunction>
451 int depth,NodeAdjacencyFunction* F,
int processCurrent)
453 int c1[3],c2[3],w1,w2;
456 w1=node1->
width(maxDepth+1);
457 w2=node2->
width(maxDepth+1);
459 ProcessFixedDepthNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,depth,F,processCurrent);
461 template<
class NodeData,
class Real>
462 template<
class NodeAdjacencyFunction>
466 int depth,NodeAdjacencyFunction* F,
int processCurrent)
468 int d=node2->
depth();
470 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
471 if(d==depth){
if(processCurrent){F->Function(node2,node1);}}
474 __ProcessFixedDepthNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,depth-1,F);
477 template<
class NodeData,
class Real>
478 template<
class NodeAdjacencyFunction>
482 int depth,NodeAdjacencyFunction* F,
int processCurrent)
484 int c1[3],c2[3],w1,w2;
487 w1=node1->
width(maxDepth+1);
488 w2=node2->
width(maxDepth+1);
489 ProcessMaxDepthNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,depth,F,processCurrent);
491 template<
class NodeData,
class Real>
492 template<
class NodeAdjacencyFunction>
496 int depth,NodeAdjacencyFunction* F,
int processCurrent)
498 int d=node2->
depth();
500 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
501 if(processCurrent){F->Function(node2,node1);}
502 if(d<depth && node2->children){__ProcessMaxDepthNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2>>1,depth-1,F);}
504 template <
class NodeData,
class Real>
505 template<
class NodeAdjacencyFunction>
508 OctNode* node2,
int radius2,
int cWidth2,
509 NodeAdjacencyFunction* F)
511 int cWidth=cWidth2>>1;
512 int radius=radius2>>1;
513 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
521 if(o& 1){F->Function(&node2->
children[0],node1);
if(node2->
children[0].
children){__ProcessNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->
children[0],radius,cWidth,F);}}
522 if(o& 2){F->Function(&node2->
children[1],node1);
if(node2->
children[1].
children){__ProcessNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->
children[1],radius,cWidth,F);}}
523 if(o& 4){F->Function(&node2->
children[2],node1);
if(node2->
children[2].
children){__ProcessNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->
children[2],radius,cWidth,F);}}
524 if(o& 8){F->Function(&node2->
children[3],node1);
if(node2->
children[3].
children){__ProcessNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->
children[3],radius,cWidth,F);}}
525 if(o& 16){F->Function(&node2->
children[4],node1);
if(node2->
children[4].
children){__ProcessNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->
children[4],radius,cWidth,F);}}
526 if(o& 32){F->Function(&node2->
children[5],node1);
if(node2->
children[5].
children){__ProcessNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->
children[5],radius,cWidth,F);}}
527 if(o& 64){F->Function(&node2->
children[6],node1);
if(node2->
children[6].
children){__ProcessNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->
children[6],radius,cWidth,F);}}
528 if(o&128){F->Function(&node2->
children[7],node1);
if(node2->
children[7].
children){__ProcessNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->
children[7],radius,cWidth,F);}}
531 template <
class NodeData,
class Real>
532 template<
class TerminatingNodeAdjacencyFunction>
535 OctNode* node2,
int radius2,
int cWidth2,
536 TerminatingNodeAdjacencyFunction* F)
538 int cWidth=cWidth2>>1;
539 int radius=radius2>>1;
540 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
548 if(o& 1){
if(F->Function(&node2->
children[0],node1) && node2->
children[0].
children){__ProcessTerminatingNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->
children[0],radius,cWidth,F);}}
549 if(o& 2){
if(F->Function(&node2->
children[1],node1) && node2->
children[1].
children){__ProcessTerminatingNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->
children[1],radius,cWidth,F);}}
550 if(o& 4){
if(F->Function(&node2->
children[2],node1) && node2->
children[2].
children){__ProcessTerminatingNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->
children[2],radius,cWidth,F);}}
551 if(o& 8){
if(F->Function(&node2->
children[3],node1) && node2->
children[3].
children){__ProcessTerminatingNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->
children[3],radius,cWidth,F);}}
552 if(o& 16){
if(F->Function(&node2->
children[4],node1) && node2->
children[4].
children){__ProcessTerminatingNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->
children[4],radius,cWidth,F);}}
553 if(o& 32){
if(F->Function(&node2->
children[5],node1) && node2->
children[5].
children){__ProcessTerminatingNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->
children[5],radius,cWidth,F);}}
554 if(o& 64){
if(F->Function(&node2->
children[6],node1) && node2->
children[6].
children){__ProcessTerminatingNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->
children[6],radius,cWidth,F);}}
555 if(o&128){
if(F->Function(&node2->
children[7],node1) && node2->
children[7].
children){__ProcessTerminatingNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->
children[7],radius,cWidth,F);}}
558 template <
class NodeData,
class Real>
559 template<
class Po
intAdjacencyFunction>
561 OctNode* node2,
int radius2,
int cWidth2,
562 PointAdjacencyFunction* F)
564 int cWidth=cWidth2>>1;
565 int radius=radius2>>1;
566 int o=ChildOverlap(dx,dy,dz,radius,cWidth);
585 template <
class NodeData,
class Real>
586 template<
class NodeAdjacencyFunction>
589 OctNode* node2,
int radius2,
int cWidth2,
590 int depth,NodeAdjacencyFunction* F)
592 int cWidth=cWidth2>>1;
593 int radius=radius2>>1;
594 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
602 if(node2->
depth()==depth){
603 if(o& 1){F->Function(&node2->
children[0],node1);}
604 if(o& 2){F->Function(&node2->
children[1],node1);}
605 if(o& 4){F->Function(&node2->
children[2],node1);}
606 if(o& 8){F->Function(&node2->
children[3],node1);}
607 if(o& 16){F->Function(&node2->
children[4],node1);}
608 if(o& 32){F->Function(&node2->
children[5],node1);}
609 if(o& 64){F->Function(&node2->
children[6],node1);}
610 if(o&128){F->Function(&node2->
children[7],node1);}
613 if(o& 1){
if(node2->
children[0].
children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->
children[0],radius,cWidth,depth,F);}}
614 if(o& 2){
if(node2->
children[1].
children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->
children[1],radius,cWidth,depth,F);}}
615 if(o& 4){
if(node2->
children[2].
children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->
children[2],radius,cWidth,depth,F);}}
616 if(o& 8){
if(node2->
children[3].
children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->
children[3],radius,cWidth,depth,F);}}
617 if(o& 16){
if(node2->
children[4].
children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->
children[4],radius,cWidth,depth,F);}}
618 if(o& 32){
if(node2->
children[5].
children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->
children[5],radius,cWidth,depth,F);}}
619 if(o& 64){
if(node2->
children[6].
children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->
children[6],radius,cWidth,depth,F);}}
620 if(o&128){
if(node2->
children[7].
children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->
children[7],radius,cWidth,depth,F);}}
624 template <
class NodeData,
class Real>
625 template<
class NodeAdjacencyFunction>
628 OctNode* node2,
int radius2,
int cWidth2,
629 int depth,NodeAdjacencyFunction* F)
631 int cWidth=cWidth2>>1;
632 int radius=radius2>>1;
633 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
641 if(node2->
depth()<=depth){
642 if(o& 1){F->Function(&node2->
children[0],node1);}
643 if(o& 2){F->Function(&node2->
children[1],node1);}
644 if(o& 4){F->Function(&node2->
children[2],node1);}
645 if(o& 8){F->Function(&node2->
children[3],node1);}
646 if(o& 16){F->Function(&node2->
children[4],node1);}
647 if(o& 32){F->Function(&node2->
children[5],node1);}
648 if(o& 64){F->Function(&node2->
children[6],node1);}
649 if(o&128){F->Function(&node2->
children[7],node1);}
651 if(node2->
depth()<depth){
652 if(o& 1){
if(node2->
children[0].
children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->
children[0],radius,cWidth,depth,F);}}
653 if(o& 2){
if(node2->
children[1].
children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->
children[1],radius,cWidth,depth,F);}}
654 if(o& 4){
if(node2->
children[2].
children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->
children[2],radius,cWidth,depth,F);}}
655 if(o& 8){
if(node2->
children[3].
children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->
children[3],radius,cWidth,depth,F);}}
656 if(o& 16){
if(node2->
children[4].
children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->
children[4],radius,cWidth,depth,F);}}
657 if(o& 32){
if(node2->
children[5].
children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->
children[5],radius,cWidth,depth,F);}}
658 if(o& 64){
if(node2->
children[6].
children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->
children[6],radius,cWidth,depth,F);}}
659 if(o&128){
if(node2->
children[7].
children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->
children[7],radius,cWidth,depth,F);}}
663 template <
class NodeData,
class Real>
671 if(dx<w2 && dx>-w1){test =1;}
672 if(dx<w1 && dx>-w2){test|=2;}
675 if(dz<w2 && dz>-w1){test1 =test;}
676 if(dz<w1 && dz>-w2){test1|=test<<4;}
678 if(!test1){
return 0;}
679 if(dy<w2 && dy>-w1){overlap =test1;}
680 if(dy<w1 && dy>-w2){overlap|=test1<<2;}
684 template <
class NodeData,
class Real>
690 if(!children){
return this;}
691 centerAndWidth(center,width);
694 cIndex=CornerIndex(center,p);
697 if(cIndex&1){center.
coords[0]+=width/2;}
698 else {center.
coords[0]-=width/2;}
699 if(cIndex&2){center.
coords[1]+=width/2;}
700 else {center.
coords[1]-=width/2;}
701 if(cIndex&4){center.
coords[2]+=width/2;}
702 else {center.
coords[2]-=width/2;}
706 template <
class NodeData,
class Real>
710 if(!children){
return this;}
713 if(!i || temp<dist2){
718 return children[nearest].getNearestLeaf(p);
721 template <
class NodeData,
class Real>
723 int o1,o2,i1,i2,j1,j2;
727 if(o1!=o2){
return 0;}
733 case 0: dir[0]=1; dir[1]=2;
break;
734 case 1: dir[0]=0; dir[1]=2;
break;
735 case 2: dir[0]=0; dir[1]=1;
break;
737 int d1,d2,off1[3],off2[3];
740 idx1[0]=off1[dir[0]]+(1<<d1)+i1;
741 idx1[1]=off1[dir[1]]+(1<<d1)+j1;
742 idx2[0]=off2[dir[0]]+(1<<d2)+i2;
743 idx2[1]=off2[dir[1]]+(1<<d2)+j2;
752 if(idx1[0]==idx2[0] && idx1[1]==idx2[1]){
return 1;}
755 template<
class NodeData,
class Real>
763 template <
class NodeData,
class Real>
764 template<
class NodeData2>
767 if(children){
delete[] children;}
771 for(i=0;i<DIMENSION;i++){this->offset[i] = node.offset[i];}
778 template <
class NodeData,
class Real>
782 template<
class NodeData ,
class Real >
787 if( n1->
d!=n2->
d )
return int(n1->
d)-int(n2->
d);
788 else if( n1->
off[0]!=n2->
off[0] )
return int(n1->
off[0]) - int(n2->
off[0]);
789 else if( n1->
off[1]!=n2->
off[1] )
return int(n1->
off[1]) - int(n2->
off[1]);
790 else if( n1->
off[2]!=n2->
off[2] )
return int(n1->
off[2]) - int(n2->
off[2]);
797 for(
int i=0 ; i<32 ; i++ ) key |= ( ( p[0] & (1<<i) )<<(2*i) ) | ( ( p[1] & (1<<i) )<<(2*i+1) ) | ( ( p[2] & (1<<i) )<<(2*i+2) );
800 template <
class NodeData,
class Real>
805 int d1 , off1[3] , d2 , off2[3];
807 if ( d1>d2 )
return 1;
808 else if( d1<d2 )
return -1;
810 if ( k1>k2 )
return 1;
811 else if( k1<k2 )
return -1;
815 template <
class NodeData,
class Real>
820 if(n1->
d!=n2->
d){
return int(n1->
d)-int(n2->
d);}
826 if(n1->
off[0]!=n2->
off[0]){
return int(n1->
off[0])-int(n2->
off[0]);}
827 if(n1->
off[1]!=n2->
off[1]){
return int(n1->
off[1])-int(n2->
off[1]);}
828 return int(n1->
off[2])-int(n2->
off[2]);
831 template <
class NodeData,
class Real>
835 template <
class NodeData,
class Real>
839 template <
class NodeData,
class Real>
842 Real w=multiplier2+multiplier1*(1<<d);
845 fabs(
Real(offSet2[0]-(offSet1[0]<<d))-w2)>=w ||
846 fabs(
Real(offSet2[1]-(offSet1[1]<<d))-w2)>=w ||
847 fabs(
Real(offSet2[2]-(offSet1[2]<<d))-w2)>=w
851 template <
class NodeData,
class Real>
853 if(c1>=dWidth || c1<=-dWidth || c2>=dWidth || c2<=-dWidth || c3>=dWidth || c3<=-dWidth){
return 0;}
856 template <
class NodeData,
class Real>
858 template <
class NodeData,
class Real>
860 template <
class NodeData,
class Real>
862 if(!parent){
return NULL;}
863 int pIndex=int(
this-parent->children);
865 if((pIndex & (1<<dir))==(off<<dir)){
return &parent->children[pIndex];}
867 OctNode* temp=parent->__faceNeighbor(dir,off,forceChildren);
868 if(!temp){
return NULL;}
876 template <
class NodeData,
class Real>
878 if(!parent){
return NULL;}
879 int pIndex=int(
this-parent->children);
881 if((pIndex & (1<<dir))==(off<<dir)){
return &parent->children[pIndex];}
883 const OctNode* temp=parent->__faceNeighbor(dir,off);
884 if(!temp || !temp->
children){
return temp;}
885 else{
return &temp->
children[pIndex];}
889 template <
class NodeData,
class Real>
894 case 0: idx[0]=1; idx[1]=2;
break;
895 case 1: idx[0]=0; idx[1]=2;
break;
896 case 2: idx[0]=0; idx[1]=1;
break;
898 return __edgeNeighbor(o,i,idx,forceChildren);
900 template <
class NodeData,
class Real>
905 case 0: idx[0]=1; idx[1]=2;
break;
906 case 1: idx[0]=0; idx[1]=2;
break;
907 case 2: idx[0]=0; idx[1]=1;
break;
909 return __edgeNeighbor(o,i,idx);
911 template <
class NodeData,
class Real>
913 if(!parent){
return NULL;}
914 int pIndex=int(
this-parent->children);
915 int aIndex,x[DIMENSION];
918 aIndex=(~((i[0] ^ x[idx[0]]) | ((i[1] ^ x[idx[1]])<<1))) & 3;
919 pIndex^=(7 ^ (1<<o));
921 const OctNode* temp=parent->__faceNeighbor(idx[0],i[0]);
922 if(!temp || !temp->
children){
return NULL;}
923 else{
return &temp->
children[pIndex];}
926 const OctNode* temp=parent->__faceNeighbor(idx[1],i[1]);
927 if(!temp || !temp->
children){
return NULL;}
928 else{
return &temp->
children[pIndex];}
934 const OctNode* temp=parent->__edgeNeighbor(o,i,idx);
935 if(!temp || !temp->
children){
return temp;}
936 else{
return &temp->
children[pIndex];}
940 template <
class NodeData,
class Real>
942 if(!parent){
return NULL;}
943 int pIndex=int(
this-parent->children);
944 int aIndex,x[DIMENSION];
947 aIndex=(~((i[0] ^ x[idx[0]]) | ((i[1] ^ x[idx[1]])<<1))) & 3;
948 pIndex^=(7 ^ (1<<o));
950 OctNode* temp=parent->__faceNeighbor(idx[0],i[0],0);
951 if(!temp || !temp->
children){
return NULL;}
952 else{
return &temp->
children[pIndex];}
955 OctNode* temp=parent->__faceNeighbor(idx[1],i[1],0);
956 if(!temp || !temp->
children){
return NULL;}
957 else{
return &temp->
children[pIndex];}
963 OctNode* temp=parent->__edgeNeighbor(o,i,idx,forceChildren);
964 if(!temp){
return NULL;}
974 template <
class NodeData,
class Real>
977 if(!parent){
return NULL;}
979 pIndex=int(
this-parent->children);
980 aIndex=(cornerIndex ^ pIndex);
983 return &parent->children[pIndex];
987 if(!temp || !temp->
children){
return temp;}
988 else{
return &temp->
children[pIndex];}
991 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(0,cornerIndex & 1);
992 if(!temp || !temp->
children){
return NULL;}
993 else{
return & temp->
children[pIndex];}
996 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(1,(cornerIndex & 2)>>1);
997 if(!temp || !temp->
children){
return NULL;}
998 else{
return & temp->
children[pIndex];}
1001 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(2,(cornerIndex & 4)>>2);
1002 if(!temp || !temp->
children){
return NULL;}
1003 else{
return & temp->
children[pIndex];}
1007 if(!temp || !temp->
children){
return NULL;}
1008 else{
return & temp->
children[pIndex];}
1012 if(!temp || !temp->
children){
return NULL;}
1013 else{
return & temp->
children[pIndex];}
1017 if(!temp || !temp->
children){
return NULL;}
1018 else{
return & temp->
children[pIndex];}
1022 template <
class NodeData,
class Real>
1024 int pIndex,aIndex=0;
1025 if(!parent){
return NULL;}
1027 pIndex=int(
this-parent->children);
1028 aIndex=(cornerIndex ^ pIndex);
1031 return &parent->children[pIndex];
1035 if(!temp){
return NULL;}
1043 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(0,cornerIndex & 1,0);
1044 if(!temp || !temp->
children){
return NULL;}
1045 else{
return & temp->
children[pIndex];}
1048 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(1,(cornerIndex & 2)>>1,0);
1049 if(!temp || !temp->
children){
return NULL;}
1050 else{
return & temp->
children[pIndex];}
1053 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(2,(cornerIndex & 4)>>2,0);
1054 if(!temp || !temp->
children){
return NULL;}
1055 else{
return & temp->
children[pIndex];}
1059 if(!temp || !temp->
children){
return NULL;}
1060 else{
return & temp->
children[pIndex];}
1064 if(!temp || !temp->
children){
return NULL;}
1065 else{
return & temp->
children[pIndex];}
1069 if(!temp || !temp->
children){
return NULL;}
1070 else{
return & temp->
children[pIndex];}
1077 template<
class NodeData,
class Real>
1079 template<
class NodeData,
class Real>
1081 for(
int i=0;i<3;i++){
for(
int j=0;j<3;j++){
for(
int k=0;k<3;k++){
neighbors[i][j][k]=NULL;}}}
1083 template<
class NodeData,
class Real>
1085 template<
class NodeData,
class Real>
1092 template<
class NodeData,
class Real>
1100 template<
class NodeData ,
class Real >
1110 Neighbors3& temp = setNeighbors( root , p , d-1 );
1112 int i , j , k , x1 , y1 , z1 , x2 , y2 , z2;
1121 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1166 i=x1<<1 , j=y1<<1 , k=z1<<1;
1176 template<
class NodeData ,
class Real >
1186 Neighbors3& temp = getNeighbors( root , p , d-1 );
1188 int i , j , k , x1 , y1 , z1 , x2 , y2 , z2;
1198 fprintf( stderr ,
"[ERROR] Couldn't find node at appropriate depth\n" );
1201 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1228 i=x1<<1 , j=y1<<1 , k=z1<<1;
1236 template<
class NodeData ,
class Real >
1239 int d = node->
depth();
1243 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
if( !neighbors[d].neighbors[i][j][k] ) reset =
true;
1244 if( reset ) neighbors[
d].neighbors[1][1][1] = NULL;
1246 if( node!=neighbors[d].neighbors[1][1][1] )
1248 neighbors[
d].clear();
1250 if( !node->
parent ) neighbors[
d].neighbors[1][1][1] = node;
1253 int i,j,k,x1,y1,z1,x2,y2,z2;
1254 int idx=int(node-node->
parent->children);
1301 i=x1<<1; j=y1<<1; k=z1<<1;
1308 return neighbors[
d];
1312 template<
class NodeData ,
class Real >
1315 int d = node->
depth();
1319 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
if( flags[i][j][k] && !neighbors[d].neighbors[i][j][k] ) reset =
true;
1320 if( reset ) neighbors[
d].neighbors[1][1][1] = NULL;
1322 if( node!=neighbors[d].neighbors[1][1][1] )
1324 neighbors[
d].clear();
1326 if( !node->
parent ) neighbors[
d].neighbors[1][1][1] = node;
1329 int x1,y1,z1,x2,y2,z2;
1330 int idx=int(node-node->
parent->children);
1333 for(
int i=0 ; i<2 ; i++ )
1334 for(
int j=0 ; j<2 ; j++ )
1335 for(
int k=0 ; k<2 ; k++ )
1368 int i=x1<<1 , j=y1<<1;
1376 int i=x1<<1 , k=z1<<1;
1384 int j=y1<<1 , k=z1<<1;
1394 int i=x1<<1 , j=y1<<1 , k=z1<<1;
1403 return neighbors[
d];
1406 template<
class NodeData,
class Real>
1408 int d=node->
depth();
1410 neighbors[
d].clear();
1412 if(!node->
parent){neighbors[
d].neighbors[1][1][1]=node;}
1414 int i,j,k,x1,y1,z1,x2,y2,z2;
1415 int idx=int(node-node->
parent->children);
1456 i=x1<<1; j=y1<<1; k=z1<<1;
1462 return neighbors[node->
depth()];
1468 template<
class NodeData,
class Real>
1470 template<
class NodeData,
class Real>
1472 for(
int i=0;i<3;i++){
for(
int j=0;j<3;j++){
for(
int k=0;k<3;k++){
neighbors[i][j][k]=NULL;}}}
1474 template<
class NodeData,
class Real>
1476 template<
class NodeData,
class Real>
1482 template<
class NodeData,
class Real>
1489 template<
class NodeData,
class Real>
1491 int d=node->
depth();
1493 neighbors[
d].clear();
1495 if(!node->
parent){neighbors[
d].neighbors[1][1][1]=node;}
1497 int i,j,k,x1,y1,z1,x2,y2,z2;
1498 int idx=int(node-node->
parent->children);
1539 i=x1<<1; j=y1<<1; k=z1<<1;
1545 return neighbors[node->
depth()];
1547 template<
class NodeData,
class Real>
1550 int d=node->
depth();
1551 if( d<minDepth ) fprintf( stderr ,
"[ERROR] Node depth lower than min-depth: %d < %d\n" , d , minDepth ) , exit( 0 );
1554 neighbors[
d].clear();
1556 if( d==minDepth ) neighbors[
d].neighbors[1][1][1]=node;
1559 int i,j,k,x1,y1,z1,x2,y2,z2;
1560 int idx = int(node-node->
parent->children);
1567 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1597 i=x1<<1 , j=y1<<1 , k=z1<<1;
1602 return neighbors[node->
depth()];
1607 template<
class NodeData ,
class Real >
1610 for(
int i=0 ; i<5 ; i++ )
for(
int j=0 ; j<5 ; j++ )
for(
int k=0 ; k<5 ; k++ )
neighbors[i][j][k] = NULL;
1612 template<
class NodeData ,
class Real >
1615 for(
int i=0 ; i<5 ; i++ )
for(
int j=0 ; j<5 ; j++ )
for(
int k=0 ; k<5 ; k++ )
neighbors[i][j][k] = NULL;
1617 template<
class NodeData ,
class Real >
1623 template<
class NodeData ,
class Real >
1629 template<
class NodeData ,
class Real >
1635 template<
class NodeData ,
class Real >
1641 template<
class NodeData ,
class Real >
1650 template<
class NodeData ,
class Real >
1659 template<
class NodeData ,
class Real >
1662 int d=node->
depth();
1665 neighbors[
d].clear();
1667 if( !node->
parent ) neighbors[
d].neighbors[2][2][2]=node;
1670 getNeighbors( node->
parent );
1672 int x1 , y1 , z1 , x2 , y2 , z2;
1679 int fx0 = x2+1 , fy0 = y2+1 , fz0 = z2+1;
1680 int cx1 = x1*2+1 , cy1 = y1*2+1 , cz1 = z1*2+1;
1681 int cx2 = x2*2+1 , cy2 = y2*2+1 , cz2 = z2*2+1;
1682 int fx1 = x1*3 , fy1 = y1*3 , fz1 = z1*3;
1683 int fx2 = x2*4 , fy2 = y2*4 , fz2 = z2*4;
1686 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1691 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1694 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1697 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1700 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1703 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1706 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1711 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1714 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1717 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1720 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1723 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1726 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1729 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1732 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1735 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1738 for( k=0 ; k<2 ; k++ )
1741 for( j=0 ; j<2 ; j++ )
1744 for( i=0 ; i<2 ; i++ )
1749 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1752 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1755 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1758 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1761 for( i=0 ; i<2 ; i++ )
1764 for( j=0 ; j<2 ; j++ )
1767 for( k=0 ; k<2 ; k++ )
1773 return neighbors[
d];
1775 template<
class NodeData ,
class Real >
1778 int d=node->
depth();
1781 neighbors[
d].clear();
1783 if( !node->
parent ) neighbors[
d].neighbors[2][2][2]=node;
1786 setNeighbors( node->
parent , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1788 int x1 , y1 , z1 , x2 , y2 , z2 , ii , jj , kk;
1792 for(
int i=xStart ; i<xEnd ; i++ )
1797 for(
int j=yStart ; j<yEnd ; j++ )
1802 for(
int k=zStart ; k<zEnd ; k++ )
1817 return neighbors[
d];
1819 template<
class NodeData ,
class Real >
1822 int d=node->
depth();
1825 neighbors[
d].clear();
1827 if(!node->
parent) neighbors[
d].neighbors[2][2][2]=node;
1830 getNeighbors( node->
parent );
1832 int x1,y1,z1,x2,y2,z2,ii,jj,kk;
1836 for(
int i=0;i<5;i++)
1841 for(
int j=0;j<5;j++)
1846 for(
int k=0;k<5;k++)
1858 return neighbors[
d];
1862 template <
class NodeData,
class Real>
1864 FILE* fp=fopen(fileName,
"wb");
1870 template <
class NodeData,
class Real>
1876 template <
class NodeData,
class Real>
1878 FILE* fp=fopen(fileName,
"rb");
1884 template <
class NodeData,
class Real>
1898 template<
class NodeData,
class Real>
1901 return 1<<(maxDepth-
d);
1903 template<
class NodeData,
class Real>
const OctNode * neighbors[5][5][5]
OctNode * neighbors[3][3][3]
static void FactorCornerIndex(int idx, int &x, int &y, int &z)
void processNodeNodes(OctNode *node, NodeAdjacencyFunction *F, int processCurrent=1)
int width(int maxDepth) const
void printRange(void) const
const OctNode * root(void) const
static const int OffsetMask
OctNode * edgeNeighbor(int edgeIndex, int forceChildren=0)
long long _InterleaveBits(int p[3])
static int CompareBackwardDepths(const void *v1, const void *v2)
static void ProcessMaxDepthNodeAdjacentNodes(int maxDepth, OctNode *node1, int width1, OctNode *node2, int width2, int depth, NodeAdjacencyFunction *F, int processCurrent=1)
static int CompareBackwardPointerDepths(const void *v1, const void *v2)
const OctNode * nextBranch(const OctNode *current) const
static Allocator< OctNode > internalAllocator
static int CommonEdge(const OctNode *node1, int eIndex1, const OctNode *node2, int eIndex2)
static void ProcessNodeAdjacentNodes(int maxDepth, OctNode *node1, int width1, OctNode *node2, int width2, NodeAdjacencyFunction *F, int processCurrent=1)
const OctNode * neighbors[3][3][3]
bool isInside(Point3D< Real > p) const
ConstNeighbors3 & getNeighbors(const OctNode *node)
const OctNode * prevBranch(const OctNode *current) const
OctNode * faceNeighbor(int faceIndex, int forceChildren=0)
static void FaceCorners(int idx, int &c1, int &c2, int &c3, int &c4)
static int CornerIndex(const Point3D< Real > ¢er, const Point3D< Real > &p)
int maxDepthLeaves(int maxDepth) const
static int UseAllocator(void)
void processNodeFaces(OctNode *node, NodeAdjacencyFunction *F, int fIndex, int processCurrent=1)
OctNode * cornerNeighbor(int cornerIndex, int forceChildren=0)
int write(const char *fileName) const
static void FactorEdgeIndex(int idx, int &orientation, int &i, int &j)
void centerAndWidth(Point3D< Real > ¢er, Real &width) const
int read(const char *fileName)
void setFullDepth(int maxDepth)
static int CompareForwardDepths(const void *v1, const void *v2)
static const int DepthShift
static const int OffsetShift2
static void Index(int depth, const int offset[3], short &d, short off[DIMENSION])
static void EdgeCorners(int idx, int &c1, int &c2)
const OctNode * nextLeaf(const OctNode *currentLeaf=NULL) const
static void ProcessFixedDepthNodeAdjacentNodes(int maxDepth, OctNode *node1, int width1, OctNode *node2, int width2, int depth, NodeAdjacencyFunction *F, int processCurrent=1)
OctNode * getNearestLeaf(const Point3D< Real > &p)
void processNodeEdges(OctNode *node, NodeAdjacencyFunction *F, int eIndex, int processCurrent=1)
static int CompareByDepthAndXYZ(const void *v1, const void *v2)
static int CompareForwardPointerDepths(const void *v1, const void *v2)
static int CornerIndex(int depth, int offSet)
ConstNeighbors5 & getNeighbors(const OctNode *node)
OctNode & operator=(const OctNode< NodeData2, Real > &node)
static int Depth(const long long &index)
static void DepthAndOffset(const long long &index, int &depth, int offset[DIMENSION])
void processNodeCorners(OctNode *node, NodeAdjacencyFunction *F, int cIndex, int processCurrent=1)
static void CenterAndWidth(const long long &index, Point3D< Real > ¢er, Real &width)
static const int OffsetShift
static void ProcessTerminatingNodeAdjacentNodes(int maxDepth, OctNode *node1, int width1, OctNode *node2, int width2, TerminatingNodeAdjacencyFunction *F, int processCurrent=1)
Neighbors3 & setNeighbors(OctNode *root, Point3D< Real > p, int d)
static int Overlap2(const int &depth1, const int offSet1[DIMENSION], const Real &multiplier1, const int &depth2, const int offSet2[DIMENSION], const Real &multiplier2)
static void ProcessPointAdjacentNodes(int maxDepth, const int center1[3], OctNode *node2, int width2, PointAdjacencyFunction *F, int processCurrent=1)
static void SetAllocator(int blockSize)
Neighbors3 & getNeighbors(OctNode *root, Point3D< Real > p, int d)
static const int OffsetShift1
Neighbors5 & getNeighbors(OctNode *node)
static int CornerIndex(int x, int y, int z)
void depthAndOffset(int &depth, int offset[DIMENSION]) const
void centerIndex(int maxDepth, int index[DIMENSION]) const
static const int OffsetShift3
const OctNode * nextNode(const OctNode *currentNode=NULL) const
static int CompareByDepthAndZIndex(const void *v1, const void *v2)
double SquareDistance(const Point3D< Real > &p1, const Point3D< Real > &p2)
OctNode * neighbors[5][5][5]
static const int DepthMask
Neighbors5 & setNeighbors(OctNode *node, int xStart=0, int xEnd=5, int yStart=0, int yEnd=5, int zStart=0, int zEnd=5)