33 #include "poisson_exceptions.h" 51 template<
class NodeData,
class Real>
int OctNode<NodeData,Real>::UseAlloc=0;
54 template<
class NodeData,
class Real>
60 internalAllocator.set(blockSize);
65 template<
class NodeData,
class Real>
68 template <
class NodeData,
class Real>
71 d=off[0]=off[1]=off[2]=0;
74 template <
class NodeData,
class Real>
76 if(!UseAlloc){
delete[] children;}
80 template <
class NodeData,
class Real>
84 if( !children ) initChildren();
85 for(
int i=0 ; i<8 ; i++ ) children[i].setFullDepth( maxDepth-1 );
89 template <
class NodeData,
class Real>
93 if(UseAlloc){children=internalAllocator.newElements(8);}
103 depthAndOffset(d,off);
108 children[idx].parent=
this;
109 children[idx].children=NULL;
111 off2[0]=(off[0]<<1)+i;
112 off2[1]=(off[1]<<1)+j;
113 off2[2]=(off[2]<<1)+k;
114 Index(d+1,off2,children[idx].d,children[idx].off);
121 template <
class NodeData,
class Real>
124 off[0]=short((1<<depth)+offset[0]-1);
125 off[1]=short((1<<depth)+offset[1]-1);
126 off[2]=short((1<<depth)+offset[2]-1);
129 template<
class NodeData,
class Real>
132 offset[0]=(int(off[0])+1)&(~(1<<depth));
133 offset[1]=(int(off[1])+1)&(~(1<<depth));
134 offset[2]=(int(off[2])+1)&(~(1<<depth));
137 template<
class NodeData,
class Real>
139 template<
class NodeData,
class Real>
141 depth=int(index&DepthMask);
142 offset[0]=(int((index>>OffsetShift1)&OffsetMask)+1)&(~(1<<depth));
143 offset[1]=(int((index>>OffsetShift2)&OffsetMask)+1)&(~(1<<depth));
144 offset[2]=(int((index>>OffsetShift3)&OffsetMask)+1)&(~(1<<depth));
147 template<
class NodeData,
class Real>
149 template <
class NodeData,
class Real>
153 offset[0]=(int(off[0])+1)&(~(1<<depth));
154 offset[1]=(int(off[1])+1)&(~(1<<depth));
155 offset[2]=(int(off[2])+1)&(~(1<<depth));
156 width=
Real(1.0/(1<<depth));
157 for(
int dim=0;dim<DIMENSION;dim++){center.
coords[dim]=
Real(0.5+offset[dim])*width;}
160 template<
class NodeData ,
class Real >
165 centerAndWidth( c , w );
167 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);
170 template <
class NodeData,
class Real>
173 depth=index&DepthMask;
174 offset[0]=(int((index>>OffsetShift1)&OffsetMask)+1)&(~(1<<depth));
175 offset[1]=(int((index>>OffsetShift2)&OffsetMask)+1)&(~(1<<depth));
176 offset[2]=(int((index>>OffsetShift3)&OffsetMask)+1)&(~(1<<depth));
177 width=
Real(1.0/(1<<depth));
178 for(
int dim=0;dim<DIMENSION;dim++){center.
coords[dim]=
Real(0.5+offset[dim])*width;}
181 template <
class NodeData,
class Real>
183 if(!children){
return 0;}
187 d=children[i].maxDepth();
194 template <
class NodeData,
class Real>
196 if(!children){
return 1;}
204 template <
class NodeData,
class Real>
206 if(!children){
return 1;}
214 template<
class NodeData,
class Real>
216 if(depth()>maxDepth){
return 0;}
217 if(!children){
return 1;}
220 for(
int i=0;i<
Cube::CORNERS;i++){c+=children[i].maxDepthLeaves(maxDepth);}
225 template <
class NodeData,
class Real>
232 template <
class NodeData,
class Real>
235 if( !current->
parent || current==
this )
return NULL;
237 else return current+1;
240 template <
class NodeData,
class Real>
242 if(!current->
parent || current==
this){
return NULL;}
244 else{
return current+1;}
247 template<
class NodeData ,
class Real >
250 if( !current->
parent || current==
this )
return NULL;
252 else return current-1;
255 template<
class NodeData ,
class Real >
258 if( !current->
parent || current==
this )
return NULL;
260 else return current-1;
263 template <
class NodeData,
class Real>
271 const OctNode* temp=nextBranch(current);
272 if(!temp){
return NULL;}
276 template <
class NodeData,
class Real>
284 OctNode* temp=nextBranch(current);
285 if(!temp){
return NULL;}
289 template <
class NodeData,
class Real>
292 if( !current )
return this;
294 else return nextBranch(current);
297 template <
class NodeData,
class Real>
300 if( !current )
return this;
302 else return nextBranch( current );
305 template <
class NodeData,
class Real>
309 centerAndWidth(center,width);
310 for(
int dim=0;dim<DIMENSION;dim++){
311 printf(
"%[%f,%f]",center.
coords[dim]-width/2,center.
coords[dim]+width/2);
312 if(dim<DIMENSION-1){printf(
"x");}
317 template <
class NodeData,
class Real>
320 template <
class NodeData,
class Real>
321 template<
class NodeAdjacencyFunction>
323 if(processCurrent){F->Function(
this,node);}
324 if(children){__processNodeNodes(node,F);}
327 template <
class NodeData,
class Real>
328 template<
class NodeAdjacencyFunction>
330 if(processCurrent){F->Function(
this,node);}
334 __processNodeFaces(node,F,c1,c2,c3,c4);
338 template <
class NodeData,
class Real>
339 template<
class NodeAdjacencyFunction>
341 if(processCurrent){F->Function(
this,node);}
345 __processNodeEdges(node,F,c1,c2);
349 template <
class NodeData,
class Real>
350 template<
class NodeAdjacencyFunction>
352 if(processCurrent){F->Function(
this,node);}
356 F->Function(temp,node);
360 template <
class NodeData,
class Real>
361 template<
class NodeAdjacencyFunction>
363 F->Function(&children[0],node);
364 F->Function(&children[1],node);
365 F->Function(&children[2],node);
366 F->Function(&children[3],node);
367 F->Function(&children[4],node);
368 F->Function(&children[5],node);
369 F->Function(&children[6],node);
370 F->Function(&children[7],node);
371 if(children[0].children){children[0].__processNodeNodes(node,F);}
372 if(children[1].children){children[1].__processNodeNodes(node,F);}
373 if(children[2].children){children[2].__processNodeNodes(node,F);}
374 if(children[3].children){children[3].__processNodeNodes(node,F);}
375 if(children[4].children){children[4].__processNodeNodes(node,F);}
376 if(children[5].children){children[5].__processNodeNodes(node,F);}
377 if(children[6].children){children[6].__processNodeNodes(node,F);}
378 if(children[7].children){children[7].__processNodeNodes(node,F);}
381 template <
class NodeData,
class Real>
382 template<
class NodeAdjacencyFunction>
383 void OctNode<NodeData,Real>::__processNodeEdges(OctNode* node,NodeAdjacencyFunction* F,
int cIndex1,
int cIndex2){
384 F->Function(&children[cIndex1],node);
385 F->Function(&children[cIndex2],node);
386 if(children[cIndex1].children){children[cIndex1].__processNodeEdges(node,F,cIndex1,cIndex2);}
387 if(children[cIndex2].children){children[cIndex2].__processNodeEdges(node,F,cIndex1,cIndex2);}
390 template <
class NodeData,
class Real>
391 template<
class NodeAdjacencyFunction>
392 void OctNode<NodeData,Real>::__processNodeFaces(OctNode* node,NodeAdjacencyFunction* F,
int cIndex1,
int cIndex2,
int cIndex3,
int cIndex4){
393 F->Function(&children[cIndex1],node);
394 F->Function(&children[cIndex2],node);
395 F->Function(&children[cIndex3],node);
396 F->Function(&children[cIndex4],node);
397 if(children[cIndex1].children){children[cIndex1].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
398 if(children[cIndex2].children){children[cIndex2].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
399 if(children[cIndex3].children){children[cIndex3].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
400 if(children[cIndex4].children){children[cIndex4].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
403 template<
class NodeData,
class Real>
404 template<
class NodeAdjacencyFunction>
406 int c1[3],c2[3],w1,w2;
409 w1=node1->
width(maxDepth+1);
410 w2=node2->
width(maxDepth+1);
412 ProcessNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,F,processCurrent);
415 template<
class NodeData,
class Real>
416 template<
class NodeAdjacencyFunction>
420 NodeAdjacencyFunction* F,
int processCurrent){
421 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
422 if(processCurrent){F->Function(node2,node1);}
424 __ProcessNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,F);
427 template<
class NodeData,
class Real>
428 template<
class TerminatingNodeAdjacencyFunction>
430 int c1[3],c2[3],w1,w2;
433 w1=node1->
width(maxDepth+1);
434 w2=node2->
width(maxDepth+1);
436 ProcessTerminatingNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,F,processCurrent);
439 template<
class NodeData,
class Real>
440 template<
class TerminatingNodeAdjacencyFunction>
444 TerminatingNodeAdjacencyFunction* F,
int processCurrent)
446 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
447 if(processCurrent){F->Function(node2,node1);}
449 __ProcessTerminatingNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,F);
452 template<
class NodeData,
class Real>
453 template<
class Po
intAdjacencyFunction>
458 w2 = node2->
width( maxDepth+1 );
459 ProcessPointAdjacentNodes( c1[0]-c2[0] , c1[1]-c2[1] , c1[2]-c2[2] , node2 , (width2*w2)>>1 , w2 , F , processCurrent );
462 template<
class NodeData,
class Real>
463 template<
class Po
intAdjacencyFunction>
466 PointAdjacencyFunction* F,
int processCurrent)
468 if( !Overlap(dx,dy,dz,radius2) )
return;
469 if( processCurrent ) F->Function(node2);
471 __ProcessPointAdjacentNodes( -dx , -dy , -dz , node2 , radius2 , width2>>1 , F );
474 template<
class NodeData,
class Real>
475 template<
class NodeAdjacencyFunction>
479 int depth,NodeAdjacencyFunction* F,
int processCurrent)
481 int c1[3],c2[3],w1,w2;
484 w1=node1->
width(maxDepth+1);
485 w2=node2->
width(maxDepth+1);
487 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);
490 template<
class NodeData,
class Real>
491 template<
class NodeAdjacencyFunction>
495 int depth,NodeAdjacencyFunction* F,
int processCurrent)
497 int d=node2->
depth();
499 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
500 if(d==depth){
if(processCurrent){F->Function(node2,node1);}}
503 __ProcessFixedDepthNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,depth-1,F);
507 template<
class NodeData,
class Real>
508 template<
class NodeAdjacencyFunction>
512 int depth,NodeAdjacencyFunction* F,
int processCurrent)
514 int c1[3],c2[3],w1,w2;
517 w1=node1->
width(maxDepth+1);
518 w2=node2->
width(maxDepth+1);
519 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);
522 template<
class NodeData,
class Real>
523 template<
class NodeAdjacencyFunction>
527 int depth,NodeAdjacencyFunction* F,
int processCurrent)
529 int d=node2->
depth();
531 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
532 if(processCurrent){F->Function(node2,node1);}
533 if(d<depth && node2->children){__ProcessMaxDepthNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2>>1,depth-1,F);}
536 template <
class NodeData,
class Real>
537 template<
class NodeAdjacencyFunction>
540 OctNode* node2,
int radius2,
int cWidth2,
541 NodeAdjacencyFunction* F)
543 int cWidth=cWidth2>>1;
544 int radius=radius2>>1;
545 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
553 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);}}
554 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);}}
555 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);}}
556 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);}}
557 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);}}
558 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);}}
559 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);}}
560 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);}}
564 template <
class NodeData,
class Real>
565 template<
class TerminatingNodeAdjacencyFunction>
566 void OctNode<NodeData,Real>::__ProcessTerminatingNodeAdjacentNodes(
int dx,
int dy,
int dz,
567 OctNode* node1,
int radius1,
568 OctNode* node2,
int radius2,
int cWidth2,
569 TerminatingNodeAdjacencyFunction* F)
571 int cWidth=cWidth2>>1;
572 int radius=radius2>>1;
573 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
581 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);}}
582 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);}}
583 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);}}
584 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);}}
585 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);}}
586 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);}}
587 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);}}
588 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);}}
592 template <
class NodeData,
class Real>
593 template<
class Po
intAdjacencyFunction>
594 void OctNode<NodeData,Real>::__ProcessPointAdjacentNodes(
int dx,
int dy,
int dz,
595 OctNode* node2,
int radius2,
int cWidth2,
596 PointAdjacencyFunction* F)
598 int cWidth=cWidth2>>1;
599 int radius=radius2>>1;
600 int o=ChildOverlap(dx,dy,dz,radius,cWidth);
609 if(o& 1){F->Function(&node2->children[0]);
if(node2->children[0].children){__ProcessPointAdjacentNodes(dx1,dy1,dz1,&node2->children[0],radius,cWidth,F);}}
610 if(o& 2){F->Function(&node2->children[1]);
if(node2->children[1].children){__ProcessPointAdjacentNodes(dx2,dy1,dz1,&node2->children[1],radius,cWidth,F);}}
611 if(o& 4){F->Function(&node2->children[2]);
if(node2->children[2].children){__ProcessPointAdjacentNodes(dx1,dy2,dz1,&node2->children[2],radius,cWidth,F);}}
612 if(o& 8){F->Function(&node2->children[3]);
if(node2->children[3].children){__ProcessPointAdjacentNodes(dx2,dy2,dz1,&node2->children[3],radius,cWidth,F);}}
613 if(o& 16){F->Function(&node2->children[4]);
if(node2->children[4].children){__ProcessPointAdjacentNodes(dx1,dy1,dz2,&node2->children[4],radius,cWidth,F);}}
614 if(o& 32){F->Function(&node2->children[5]);
if(node2->children[5].children){__ProcessPointAdjacentNodes(dx2,dy1,dz2,&node2->children[5],radius,cWidth,F);}}
615 if(o& 64){F->Function(&node2->children[6]);
if(node2->children[6].children){__ProcessPointAdjacentNodes(dx1,dy2,dz2,&node2->children[6],radius,cWidth,F);}}
616 if(o&128){F->Function(&node2->children[7]);
if(node2->children[7].children){__ProcessPointAdjacentNodes(dx2,dy2,dz2,&node2->children[7],radius,cWidth,F);}}
620 template <
class NodeData,
class Real>
621 template<
class NodeAdjacencyFunction>
622 void OctNode<NodeData,Real>::__ProcessFixedDepthNodeAdjacentNodes(
int dx,
int dy,
int dz,
623 OctNode* node1,
int radius1,
624 OctNode* node2,
int radius2,
int cWidth2,
625 int depth,NodeAdjacencyFunction* F)
627 int cWidth=cWidth2>>1;
628 int radius=radius2>>1;
629 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
637 if(node2->depth()==depth){
638 if(o& 1){F->Function(&node2->children[0],node1);}
639 if(o& 2){F->Function(&node2->children[1],node1);}
640 if(o& 4){F->Function(&node2->children[2],node1);}
641 if(o& 8){F->Function(&node2->children[3],node1);}
642 if(o& 16){F->Function(&node2->children[4],node1);}
643 if(o& 32){F->Function(&node2->children[5],node1);}
644 if(o& 64){F->Function(&node2->children[6],node1);}
645 if(o&128){F->Function(&node2->children[7],node1);}
648 if(o& 1){
if(node2->children[0].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->children[0],radius,cWidth,depth,F);}}
649 if(o& 2){
if(node2->children[1].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->children[1],radius,cWidth,depth,F);}}
650 if(o& 4){
if(node2->children[2].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->children[2],radius,cWidth,depth,F);}}
651 if(o& 8){
if(node2->children[3].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->children[3],radius,cWidth,depth,F);}}
652 if(o& 16){
if(node2->children[4].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->children[4],radius,cWidth,depth,F);}}
653 if(o& 32){
if(node2->children[5].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->children[5],radius,cWidth,depth,F);}}
654 if(o& 64){
if(node2->children[6].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->children[6],radius,cWidth,depth,F);}}
655 if(o&128){
if(node2->children[7].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->children[7],radius,cWidth,depth,F);}}
660 template <
class NodeData,
class Real>
661 template<
class NodeAdjacencyFunction>
662 void OctNode<NodeData,Real>::__ProcessMaxDepthNodeAdjacentNodes(
int dx,
int dy,
int dz,
663 OctNode* node1,
int radius1,
664 OctNode* node2,
int radius2,
int cWidth2,
665 int depth,NodeAdjacencyFunction* F)
667 int cWidth=cWidth2>>1;
668 int radius=radius2>>1;
669 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
677 if(node2->depth()<=depth){
678 if(o& 1){F->Function(&node2->children[0],node1);}
679 if(o& 2){F->Function(&node2->children[1],node1);}
680 if(o& 4){F->Function(&node2->children[2],node1);}
681 if(o& 8){F->Function(&node2->children[3],node1);}
682 if(o& 16){F->Function(&node2->children[4],node1);}
683 if(o& 32){F->Function(&node2->children[5],node1);}
684 if(o& 64){F->Function(&node2->children[6],node1);}
685 if(o&128){F->Function(&node2->children[7],node1);}
687 if(node2->depth()<depth){
688 if(o& 1){
if(node2->children[0].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->children[0],radius,cWidth,depth,F);}}
689 if(o& 2){
if(node2->children[1].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->children[1],radius,cWidth,depth,F);}}
690 if(o& 4){
if(node2->children[2].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->children[2],radius,cWidth,depth,F);}}
691 if(o& 8){
if(node2->children[3].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->children[3],radius,cWidth,depth,F);}}
692 if(o& 16){
if(node2->children[4].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->children[4],radius,cWidth,depth,F);}}
693 if(o& 32){
if(node2->children[5].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->children[5],radius,cWidth,depth,F);}}
694 if(o& 64){
if(node2->children[6].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->children[6],radius,cWidth,depth,F);}}
695 if(o&128){
if(node2->children[7].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->children[7],radius,cWidth,depth,F);}}
700 template <
class NodeData,
class Real>
701 inline int OctNode<NodeData,Real>::ChildOverlap(
int dx,
int dy,
int dz,
int d,
int cRadius2)
708 if(dx<w2 && dx>-w1){test =1;}
709 if(dx<w1 && dx>-w2){test|=2;}
712 if(dz<w2 && dz>-w1){test1 =test;}
713 if(dz<w1 && dz>-w2){test1|=test<<4;}
715 if(!test1){
return 0;}
716 if(dy<w2 && dy>-w1){overlap =test1;}
717 if(dy<w1 && dy>-w2){overlap|=test1<<2;}
721 template <
class NodeData,
class Real>
727 if(!children){
return this;}
728 centerAndWidth(center,width);
731 cIndex=CornerIndex(center,p);
734 if(cIndex&1){center.
coords[0]+=width/2;}
735 else {center.
coords[0]-=width/2;}
736 if(cIndex&2){center.
coords[1]+=width/2;}
737 else {center.
coords[1]-=width/2;}
738 if(cIndex&4){center.
coords[2]+=width/2;}
739 else {center.
coords[2]-=width/2;}
744 template <
class NodeData,
class Real>
748 if(!children){
return this;}
751 if(!i || temp<dist2){
756 return children[nearest].getNearestLeaf(p);
759 template <
class NodeData,
class Real>
761 int o1,o2,i1,i2,j1,j2;
765 if(o1!=o2){
return 0;}
771 case 0: dir[0]=1; dir[1]=2;
break;
772 case 1: dir[0]=0; dir[1]=2;
break;
773 case 2: dir[0]=0; dir[1]=1;
break;
775 int d1,d2,off1[3],off2[3];
778 idx1[0]=off1[dir[0]]+(1<<d1)+i1;
779 idx1[1]=off1[dir[1]]+(1<<d1)+j1;
780 idx2[0]=off2[dir[0]]+(1<<d2)+i2;
781 idx2[1]=off2[dir[1]]+(1<<d2)+j2;
790 if(idx1[0]==idx2[0] && idx1[1]==idx2[1]){
return 1;}
794 template<
class NodeData,
class Real>
803 template <
class NodeData,
class Real>
804 template<
class NodeData2>
811 for(i=0;i<DIMENSION;i++){this->offset[i] = node.offset[i];}
819 template <
class NodeData,
class Real>
824 template<
class NodeData ,
class Real >
829 if( n1->
d!=n2->
d )
return int(n1->
d)-int(n2->
d);
830 else if( n1->
off[0]!=n2->
off[0] )
return int(n1->
off[0]) - int(n2->
off[0]);
831 else if( n1->
off[1]!=n2->
off[1] )
return int(n1->
off[1]) - int(n2->
off[1]);
832 else if( n1->
off[2]!=n2->
off[2] )
return int(n1->
off[2]) - int(n2->
off[2]);
839 long long _p[3] = {p[0],p[1],p[2]};
840 for(
int i=0 ; i<31 ; i++ ) key |= ( ( _p[0] & (1ull<<i) )<<(2*i) ) | ( ( _p[1] & (1ull<<i) )<<(2*i+1) ) | ( ( _p[2] & (1ull<<i) )<<(2*i+2) );
844 template <
class NodeData,
class Real>
849 int d1 , off1[3] , d2 , off2[3];
851 if ( d1>d2 )
return 1;
852 else if( d1<d2 )
return -1;
854 if ( k1>k2 )
return 1;
855 else if( k1<k2 )
return -1;
859 template <
class NodeData,
class Real>
864 if(n1->
d!=n2->
d){
return int(n1->
d)-int(n2->
d);}
870 if(n1->
off[0]!=n2->
off[0]){
return int(n1->
off[0])-int(n2->
off[0]);}
871 if(n1->
off[1]!=n2->
off[1]){
return int(n1->
off[1])-int(n2->
off[1]);}
872 return int(n1->
off[2])-int(n2->
off[2]);
876 template <
class NodeData,
class Real>
881 template <
class NodeData,
class Real>
886 template <
class NodeData,
class Real>
889 Real w=multiplier2+multiplier1*(1<<d);
892 fabs(
Real(offSet2[0]-(offSet1[0]<<d))-w2)>=w ||
893 fabs(
Real(offSet2[1]-(offSet1[1]<<d))-w2)>=w ||
894 fabs(
Real(offSet2[2]-(offSet1[2]<<d))-w2)>=w
899 template <
class NodeData,
class Real>
901 if(c1>=dWidth || c1<=-dWidth || c2>=dWidth || c2<=-dWidth || c3>=dWidth || c3<=-dWidth){
return 0;}
905 template <
class NodeData,
class Real>
907 template <
class NodeData,
class Real>
909 template <
class NodeData,
class Real>
911 if(!parent){
return NULL;}
912 int pIndex=int(
this-(parent->children));
914 if((pIndex & (1<<dir))==(off<<dir)){
return &parent->children[pIndex];}
916 OctNode* temp=parent->__faceNeighbor(dir,off,forceChildren);
917 if(!temp){
return NULL;}
919 if(forceChildren){temp->initChildren();}
922 return &temp->children[pIndex];
926 template <
class NodeData,
class Real>
927 const OctNode<NodeData,Real>* OctNode<NodeData,Real>::__faceNeighbor(
int dir,
int off)
const {
928 if(!parent){
return NULL;}
929 int pIndex=int(
this-(parent->children));
931 if((pIndex & (1<<dir))==(off<<dir)){
return &parent->children[pIndex];}
933 const OctNode* temp=parent->__faceNeighbor(dir,off);
934 if(!temp || !temp->children){
return temp;}
935 else{
return &temp->children[pIndex];}
939 template <
class NodeData,
class Real>
944 case 0: idx[0]=1; idx[1]=2;
break;
945 case 1: idx[0]=0; idx[1]=2;
break;
946 case 2: idx[0]=0; idx[1]=1;
break;
948 return __edgeNeighbor(o,i,idx,forceChildren);
951 template <
class NodeData,
class Real>
956 case 0: idx[0]=1; idx[1]=2;
break;
957 case 1: idx[0]=0; idx[1]=2;
break;
958 case 2: idx[0]=0; idx[1]=1;
break;
960 return __edgeNeighbor(o,i,idx);
963 template <
class NodeData,
class Real>
965 if(!parent){
return NULL;}
966 int pIndex=int(
this-(parent->children));
967 int aIndex,x[DIMENSION];
970 aIndex=(~((i[0] ^ x[idx[0]]) | ((i[1] ^ x[idx[1]])<<1))) & 3;
971 pIndex^=(7 ^ (1<<o));
973 const OctNode* temp=parent->__faceNeighbor(idx[0],i[0]);
974 if(!temp || !temp->children){
return NULL;}
975 else{
return &temp->children[pIndex];}
978 const OctNode* temp=parent->__faceNeighbor(idx[1],i[1]);
979 if(!temp || !temp->children){
return NULL;}
980 else{
return &temp->children[pIndex];}
983 return &parent->children[pIndex];
986 const OctNode* temp=parent->__edgeNeighbor(o,i,idx);
987 if(!temp || !temp->children){
return temp;}
988 else{
return &temp->children[pIndex];}
993 template <
class NodeData,
class Real>
994 OctNode<NodeData,Real>* OctNode<NodeData,Real>::__edgeNeighbor(
int o,
const int i[2],
const int idx[2],
int forceChildren){
995 if(!parent){
return NULL;}
996 int pIndex=int(
this-(parent->children));
997 int aIndex,x[DIMENSION];
1000 aIndex=(~((i[0] ^ x[idx[0]]) | ((i[1] ^ x[idx[1]])<<1))) & 3;
1001 pIndex^=(7 ^ (1<<o));
1003 OctNode* temp=parent->__faceNeighbor(idx[0],i[0],0);
1004 if(!temp || !temp->children){
return NULL;}
1005 else{
return &temp->children[pIndex];}
1007 else if(aIndex==2) {
1008 OctNode* temp=parent->__faceNeighbor(idx[1],i[1],0);
1009 if(!temp || !temp->children){
return NULL;}
1010 else{
return &temp->children[pIndex];}
1012 else if(aIndex==0) {
1013 return &parent->children[pIndex];
1015 else if(aIndex==3) {
1016 OctNode* temp=parent->__edgeNeighbor(o,i,idx,forceChildren);
1017 if(!temp){
return NULL;}
1018 if(!temp->children){
1019 if(forceChildren){temp->initChildren();}
1022 return &temp->children[pIndex];
1027 template <
class NodeData,
class Real>
1029 int pIndex,aIndex=0;
1030 if(!parent){
return NULL;}
1032 pIndex=int(
this-(parent->children));
1033 aIndex=(cornerIndex ^ pIndex);
1036 return &parent->children[pIndex];
1040 if(!temp || !temp->
children){
return temp;}
1041 else{
return &temp->
children[pIndex];}
1044 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(0,cornerIndex & 1);
1045 if(!temp || !temp->
children){
return NULL;}
1046 else{
return & temp->
children[pIndex];}
1049 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(1,(cornerIndex & 2)>>1);
1050 if(!temp || !temp->
children){
return NULL;}
1051 else{
return & temp->
children[pIndex];}
1054 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(2,(cornerIndex & 4)>>2);
1055 if(!temp || !temp->
children){
return NULL;}
1056 else{
return & temp->
children[pIndex];}
1060 if(!temp || !temp->
children){
return NULL;}
1061 else{
return & temp->
children[pIndex];}
1065 if(!temp || !temp->
children){
return NULL;}
1066 else{
return & temp->
children[pIndex];}
1070 if(!temp || !temp->
children){
return NULL;}
1071 else{
return & temp->
children[pIndex];}
1076 template <
class NodeData,
class Real>
1078 int pIndex,aIndex=0;
1079 if(!parent){
return NULL;}
1081 pIndex=int(
this-(parent->children));
1082 aIndex=(cornerIndex ^ pIndex);
1085 return &parent->children[pIndex];
1089 if(!temp){
return NULL;}
1097 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(0,cornerIndex & 1,0);
1098 if(!temp || !temp->
children){
return NULL;}
1099 else{
return & temp->
children[pIndex];}
1102 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(1,(cornerIndex & 2)>>1,0);
1103 if(!temp || !temp->
children){
return NULL;}
1104 else{
return & temp->
children[pIndex];}
1107 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(2,(cornerIndex & 4)>>2,0);
1108 if(!temp || !temp->
children){
return NULL;}
1109 else{
return & temp->
children[pIndex];}
1113 if(!temp || !temp->
children){
return NULL;}
1114 else{
return & temp->
children[pIndex];}
1118 if(!temp || !temp->
children){
return NULL;}
1119 else{
return & temp->
children[pIndex];}
1123 if(!temp || !temp->
children){
return NULL;}
1124 else{
return & temp->
children[pIndex];}
1132 template<
class NodeData,
class Real>
1135 template<
class NodeData,
class Real>
1137 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;}}}
1140 template<
class NodeData,
class Real>
1143 template<
class NodeData,
class Real>
1150 template<
class NodeData,
class Real>
1159 template<
class NodeData ,
class Real >
1162 if( !neighbors[
d].neighbors[1][1][1] || !neighbors[
d].neighbors[1][1][1]->
isInside( p ) )
1166 if( !
d ) neighbors[
d].neighbors[1][1][1] =
root;
1171 int i , j , k , x1 , y1 , z1 , x2 , y2 , z2;
1180 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1225 i=x1<<1 , j=y1<<1 , k=z1<<1;
1233 return neighbors[
d];
1236 template<
class NodeData ,
class Real >
1239 if( !neighbors[
d].neighbors[1][1][1] || !neighbors[
d].neighbors[1][1][1]->
isInside( p ) )
1243 if( !
d ) neighbors[
d].neighbors[1][1][1] =
root;
1246 Neighbors3& temp = getNeighbors(
root , p ,
d-1 );
1248 int i , j , k , x1 , y1 , z1 , x2 , y2 , z2;
1251 temp.neighbors[1][1][1]->centerAndWidth( c , w );
1256 if( !temp.neighbors[1][1][1] || !temp.neighbors[1][1][1]->children )
1260 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1261 neighbors[
d].neighbors[x2+i][y2+j][z2+k] = &temp.neighbors[1][1][1]->children[
Cube::CornerIndex(i,j,k)];
1266 if( temp.neighbors[i][1][1] )
1267 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][y2+j][z2+k] = &temp.neighbors[i][1][1]->children[
Cube::CornerIndex(x2,j,k)];
1269 if( temp.neighbors[1][j][1] )
1270 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[x2+i][j][z2+k] = &temp.neighbors[1][j][1]->children[
Cube::CornerIndex(i,y2,k)];
1272 if( temp.neighbors[1][1][k] )
1273 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[x2+i][y2+j][k] = &temp.neighbors[1][1][k]->children[
Cube::CornerIndex(i,j,z2)];
1277 if( temp.neighbors[i][j][1] )
1278 for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][j][z2+k] = &temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];
1280 if( temp.neighbors[i][1][k] )
1281 for( j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[i][y2+j][k] = &temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];
1283 if( temp.neighbors[1][j][k] )
1284 for( i=0 ; i<2 ; i++ ) neighbors[
d].neighbors[x2+i][j][k] = &temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];
1287 i=x1<<1 , j=y1<<1 , k=z1<<1;
1288 if( temp.neighbors[i][j][k] )
1289 neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1292 return neighbors[
d];
1295 template<
class NodeData ,
class Real >
1298 int d = node->depth();
1299 if( node==neighbors[
d].neighbors[1][1][1] )
1302 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;
1303 if( reset ) neighbors[
d].neighbors[1][1][1] = NULL;
1305 if( node!=neighbors[
d].neighbors[1][1][1] )
1307 neighbors[
d].clear();
1309 if( !node->parent ) neighbors[
d].neighbors[1][1][1] = node;
1312 int i,j,k,x1,y1,z1,x2,y2,z2;
1313 int idx=int(node-node->parent->children);
1319 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1323 Neighbors3& temp=setNeighbors(node->parent);
1327 if(temp.neighbors[i][1][1]){
1328 if(!temp.neighbors[i][1][1]->children){temp.neighbors[i][1][1]->initChildren();}
1329 for(j=0;j<2;j++){
for(k=0;k<2;k++){neighbors[
d].neighbors[i][y2+j][z2+k]=&temp.neighbors[i][1][1]->children[
Cube::CornerIndex(x2,j,k)];}}
1332 if(temp.neighbors[1][j][1]){
1333 if(!temp.neighbors[1][j][1]->children){temp.neighbors[1][j][1]->initChildren();}
1334 for(i=0;i<2;i++){
for(k=0;k<2;k++){neighbors[
d].neighbors[x2+i][j][z2+k]=&temp.neighbors[1][j][1]->children[
Cube::CornerIndex(i,y2,k)];}}
1337 if(temp.neighbors[1][1][k]){
1338 if(!temp.neighbors[1][1][k]->children){temp.neighbors[1][1][k]->initChildren();}
1339 for(i=0;i<2;i++){
for(j=0;j<2;j++){neighbors[
d].neighbors[x2+i][y2+j][k]=&temp.neighbors[1][1][k]->children[
Cube::CornerIndex(i,j,z2)];}}
1344 if(temp.neighbors[i][j][1]){
1345 if(!temp.neighbors[i][j][1]->children){temp.neighbors[i][j][1]->initChildren();}
1346 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1349 if(temp.neighbors[i][1][k]){
1350 if(!temp.neighbors[i][1][k]->children){temp.neighbors[i][1][k]->initChildren();}
1351 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1354 if(temp.neighbors[1][j][k]){
1355 if(!temp.neighbors[1][j][k]->children){temp.neighbors[1][j][k]->initChildren();}
1356 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1360 i=x1<<1; j=y1<<1; k=z1<<1;
1361 if(temp.neighbors[i][j][k]){
1362 if(!temp.neighbors[i][j][k]->children){temp.neighbors[i][j][k]->initChildren();}
1363 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1367 return neighbors[
d];
1372 template<
class NodeData ,
class Real >
1375 int d = node->depth();
1376 if( node==neighbors[
d].neighbors[1][1][1] )
1379 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;
1380 if( reset ) neighbors[
d].neighbors[1][1][1] = NULL;
1382 if( node!=neighbors[
d].neighbors[1][1][1] )
1384 neighbors[
d].clear();
1386 if( !node->parent ) neighbors[
d].neighbors[1][1][1] = node;
1389 int x1,y1,z1,x2,y2,z2;
1390 int idx=int(node-node->parent->children);
1393 for(
int i=0 ; i<2 ; i++ )
1394 for(
int j=0 ; j<2 ; j++ )
1395 for(
int k=0 ; k<2 ; k++ )
1396 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1398 Neighbors3& temp=setNeighbors( node->parent , flags );
1403 if( temp.neighbors[i][1][1] )
1405 if( flags[i][1][1] && !temp.neighbors[i][1][1]->children ) temp.neighbors[i][1][1]->initChildren();
1406 if( temp.neighbors[i][1][1]->children )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][y2+j][z2+k] = &temp.neighbors[i][1][1]->children[
Cube::CornerIndex(x2,j,k)];
1411 if( temp.neighbors[1][j][1] )
1413 if( flags[1][j][1] && !temp.neighbors[1][j][1]->children ) temp.neighbors[1][j][1]->initChildren();
1414 if( temp.neighbors[1][j][1]->children )
for(
int i=0 ; i<2 ; i++ )
for(
int k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[x2+i][j][z2+k] = &temp.neighbors[1][j][1]->children[
Cube::CornerIndex(i,y2,k)];
1419 if( temp.neighbors[1][1][k] )
1421 if( flags[1][1][k] && !temp.neighbors[1][1][k]->children ) temp.neighbors[1][1][k]->initChildren();
1422 if( temp.neighbors[1][1][k]->children )
for(
int i=0 ; i<2 ; i++ )
for(
int j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[x2+i][y2+j][k] = &temp.neighbors[1][1][k]->children[
Cube::CornerIndex(i,j,z2)];
1428 int i=x1<<1 , j=y1<<1;
1429 if( temp.neighbors[i][j][1] )
1431 if( flags[i][j][1] && !temp.neighbors[i][j][1]->children ) temp.neighbors[i][j][1]->initChildren();
1432 if( temp.neighbors[i][j][1]->children )
for(
int k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][j][z2+k] = &temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];
1436 int i=x1<<1 , k=z1<<1;
1437 if( temp.neighbors[i][1][k] )
1439 if( flags[i][1][k] && !temp.neighbors[i][1][k]->children ) temp.neighbors[i][1][k]->initChildren();
1440 if( temp.neighbors[i][1][k]->children )
for(
int j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[i][y2+j][k] = &temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];
1444 int j=y1<<1 , k=z1<<1;
1445 if( temp.neighbors[1][j][k] )
1447 if( flags[1][j][k] && !temp.neighbors[1][j][k]->children ) temp.neighbors[1][j][k]->initChildren();
1448 if( temp.neighbors[1][j][k]->children )
for(
int i=0 ; i<2 ; i++ ) neighbors[
d].neighbors[x2+i][j][k] = &temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];
1454 int i=x1<<1 , j=y1<<1 , k=z1<<1;
1455 if( temp.neighbors[i][j][k] )
1457 if( flags[i][j][k] && !temp.neighbors[i][j][k]->children ) temp.neighbors[i][j][k]->initChildren();
1458 if( temp.neighbors[i][j][k]->children ) neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1463 return neighbors[
d];
1466 template<
class NodeData,
class Real>
1468 int d=node->depth();
1469 if(node!=neighbors[
d].neighbors[1][1][1]){
1470 neighbors[
d].clear();
1472 if(!node->parent){neighbors[
d].neighbors[1][1][1]=node;}
1474 int i,j,k,x1,y1,z1,x2,y2,z2;
1475 int idx=int(node-node->parent->children);
1481 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1485 Neighbors3& temp=getNeighbors(node->parent);
1489 if(temp.neighbors[i][1][1] && temp.neighbors[i][1][1]->children){
1490 for(j=0;j<2;j++){
for(k=0;k<2;k++){neighbors[
d].neighbors[i][y2+j][z2+k]=&temp.neighbors[i][1][1]->children[
Cube::CornerIndex(x2,j,k)];}}
1493 if(temp.neighbors[1][j][1] && temp.neighbors[1][j][1]->children){
1494 for(i=0;i<2;i++){
for(k=0;k<2;k++){neighbors[
d].neighbors[x2+i][j][z2+k]=&temp.neighbors[1][j][1]->children[
Cube::CornerIndex(i,y2,k)];}}
1497 if(temp.neighbors[1][1][k] && temp.neighbors[1][1][k]->children){
1498 for(i=0;i<2;i++){
for(j=0;j<2;j++){neighbors[
d].neighbors[x2+i][y2+j][k]=&temp.neighbors[1][1][k]->children[
Cube::CornerIndex(i,j,z2)];}}
1503 if(temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children){
1504 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1507 if(temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children){
1508 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1511 if(temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children){
1512 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1516 i=x1<<1; j=y1<<1; k=z1<<1;
1517 if(temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children){
1518 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1522 return neighbors[node->depth()];
1528 template<
class NodeData,
class Real>
1531 template<
class NodeData,
class Real>
1533 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;}}}
1536 template<
class NodeData,
class Real>
1539 template<
class NodeData,
class Real>
1545 template<
class NodeData,
class Real>
1553 template<
class NodeData,
class Real>
1556 if(node!=neighbors[
d].neighbors[1][1][1]){
1557 neighbors[
d].clear();
1559 if(!node->
parent){neighbors[
d].neighbors[1][1][1]=node;}
1561 int i,j,k,x1,y1,z1,x2,y2,z2;
1562 int idx=int(node-node->
parent->children);
1572 ConstNeighbors3& temp=getNeighbors(node->
parent);
1576 if(temp.neighbors[i][1][1] && temp.neighbors[i][1][1]->children){
1577 for(j=0;j<2;j++){
for(k=0;k<2;k++){neighbors[
d].neighbors[i][y2+j][z2+k]=&temp.neighbors[i][1][1]->children[
Cube::CornerIndex(x2,j,k)];}}
1580 if(temp.neighbors[1][j][1] && temp.neighbors[1][j][1]->children){
1581 for(i=0;i<2;i++){
for(k=0;k<2;k++){neighbors[
d].neighbors[x2+i][j][z2+k]=&temp.neighbors[1][j][1]->children[
Cube::CornerIndex(i,y2,k)];}}
1584 if(temp.neighbors[1][1][k] && temp.neighbors[1][1][k]->children){
1585 for(i=0;i<2;i++){
for(j=0;j<2;j++){neighbors[
d].neighbors[x2+i][y2+j][k]=&temp.neighbors[1][1][k]->children[
Cube::CornerIndex(i,j,z2)];}}
1590 if(temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children){
1591 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1594 if(temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children){
1595 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1598 if(temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children){
1599 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1603 i=x1<<1; j=y1<<1; k=z1<<1;
1604 if(temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children){
1605 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1609 return neighbors[node->
depth()];
1612 template<
class NodeData,
class Real>
1615 int d=node->depth();
1620 if( node!=neighbors[
d].neighbors[1][1][1] )
1622 neighbors[
d].clear();
1624 if(
d==minDepth ) neighbors[
d].neighbors[1][1][1]=node;
1627 int i,j,k,x1,y1,z1,x2,y2,z2;
1628 int idx = int(node-node->parent->children);
1632 ConstNeighbors3& temp=getNeighbors( node->parent , minDepth );
1635 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1636 neighbors[
d].neighbors[x2+i][y2+j][z2+k] = &node->parent->children[
Cube::CornerIndex(i,j,k) ];
1640 if( temp.neighbors[i][1][1] && temp.neighbors[i][1][1]->children )
1641 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][y2+j][z2+k] = &temp.neighbors[i][1][1]->children[
Cube::CornerIndex(x2,j,k)];
1644 if( temp.neighbors[1][j][1] && temp.neighbors[1][j][1]->children )
1645 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[x2+i][j][z2+k] = &temp.neighbors[1][j][1]->children[
Cube::CornerIndex(i,y2,k)];
1648 if( temp.neighbors[1][1][k] && temp.neighbors[1][1][k]->children )
1649 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[x2+i][y2+j][k] = &temp.neighbors[1][1][k]->children[
Cube::CornerIndex(i,j,z2)];
1653 if( temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children )
1654 for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][j][z2+k] = &temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];
1657 if( temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children )
1658 for( j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[i][y2+j][k] = &temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];
1661 if( temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children )
1662 for( i=0 ; i<2 ; i++ ) neighbors[
d].neighbors[x2+i][j][k] = &temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];
1665 i=x1<<1 , j=y1<<1 , k=z1<<1;
1666 if( temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children )
1667 neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1670 return neighbors[node->depth()];
1677 template<
class NodeData ,
class Real >
1680 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;
1683 template<
class NodeData ,
class Real >
1686 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;
1689 template<
class NodeData ,
class Real >
1696 template<
class NodeData ,
class Real >
1703 template<
class NodeData ,
class Real >
1710 template<
class NodeData ,
class Real >
1717 template<
class NodeData ,
class Real >
1727 template<
class NodeData ,
class Real >
1737 template<
class NodeData ,
class Real >
1741 if( node!=neighbors[
d].neighbors[2][2][2] )
1743 neighbors[
d].clear();
1745 if( !node->
parent ) neighbors[
d].neighbors[2][2][2]=node;
1748 getNeighbors( node->
parent );
1750 int x1 , y1 , z1 , x2 , y2 , z2;
1757 int fx0 = x2+1 , fy0 = y2+1 , fz0 = z2+1;
1758 int cx1 = x1*2+1 , cy1 = y1*2+1 , cz1 = z1*2+1;
1759 int cx2 = x2*2+1 , cy2 = y2*2+1 , cz2 = z2*2+1;
1760 int fx1 = x1*3 , fy1 = y1*3 , fz1 = z1*3;
1761 int fx2 = x2*4 , fy2 = y2*4 , fz2 = z2*4;
1764 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1769 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1772 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1775 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1778 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1781 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1784 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1789 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1792 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1795 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1798 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1801 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1804 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1807 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1810 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1813 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1816 for( k=0 ; k<2 ; k++ )
1819 for( j=0 ; j<2 ; j++ )
1822 for( i=0 ; i<2 ; i++ )
1827 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1830 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1833 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1836 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1839 for( i=0 ; i<2 ; i++ )
1842 for( j=0 ; j<2 ; j++ )
1845 for( k=0 ; k<2 ; k++ )
1851 return neighbors[
d];
1854 template<
class NodeData ,
class Real >
1858 if( node!=neighbors[
d].neighbors[2][2][2] )
1860 neighbors[
d].clear();
1862 if( !node->
parent ) neighbors[
d].neighbors[2][2][2]=node;
1865 setNeighbors( node->
parent , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1867 int x1 , y1 , z1 , x2 , y2 , z2 , ii , jj , kk;
1871 for(
int i=xStart ; i<xEnd ; i++ )
1876 for(
int j=yStart ; j<yEnd ; j++ )
1881 for(
int k=zStart ; k<zEnd ; k++ )
1896 return neighbors[
d];
1899 template<
class NodeData ,
class Real >
1903 if( node!=neighbors[
d].neighbors[2][2][2] )
1905 neighbors[
d].clear();
1907 if(!node->
parent) neighbors[
d].neighbors[2][2][2]=node;
1910 getNeighbors( node->
parent );
1912 int x1,y1,z1,x2,y2,z2,ii,jj,kk;
1916 for(
int i=0;i<5;i++)
1921 for(
int j=0;j<5;j++)
1926 for(
int k=0;k<5;k++)
1938 return neighbors[
d];
1941 template <
class NodeData,
class Real>
1943 FILE* fp=fopen(fileName,
"wb");
1950 template <
class NodeData,
class Real>
1957 template <
class NodeData,
class Real>
1959 FILE* fp=fopen(fileName,
"rb");
1966 template <
class NodeData,
class Real>
1981 template<
class NodeData,
class Real>
1987 template<
class NodeData,
class Real>
const OctNode * neighbors[5][5][5]
OctNode * neighbors[3][3][3]
An exception that is thrown when initialization fails.
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)
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
An exception that is thrown when the arguments number or type is wrong/unhandled. ...
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)