32 #include "poisson_exceptions.h"
50 template<
class NodeData,
class Real>
int OctNode<NodeData,Real>::UseAlloc=0;
53 template<
class NodeData,
class Real>
59 internalAllocator.set(blockSize);
64 template<
class NodeData,
class Real>
67 template <
class NodeData,
class Real>
70 d=off[0]=off[1]=off[2]=0;
73 template <
class NodeData,
class Real>
75 if(!UseAlloc){
delete[] children;}
79 template <
class NodeData,
class Real>
83 if( !children ) initChildren();
84 for(
int i=0 ; i<8 ; i++ ) children[i].setFullDepth( maxDepth-1 );
88 template <
class NodeData,
class Real>
92 if(UseAlloc){children=internalAllocator.newElements(8);}
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);
120 template <
class NodeData,
class Real>
123 off[0]=short((1<<depth)+offset[0]-1);
124 off[1]=short((1<<depth)+offset[1]-1);
125 off[2]=short((1<<depth)+offset[2]-1);
128 template<
class NodeData,
class Real>
131 offset[0]=(int(off[0])+1)&(~(1<<depth));
132 offset[1]=(int(off[1])+1)&(~(1<<depth));
133 offset[2]=(int(off[2])+1)&(~(1<<depth));
136 template<
class NodeData,
class Real>
138 template<
class NodeData,
class Real>
140 depth=int(index&DepthMask);
141 offset[0]=(int((index>>OffsetShift1)&OffsetMask)+1)&(~(1<<depth));
142 offset[1]=(int((index>>OffsetShift2)&OffsetMask)+1)&(~(1<<depth));
143 offset[2]=(int((index>>OffsetShift3)&OffsetMask)+1)&(~(1<<depth));
146 template<
class NodeData,
class Real>
148 template <
class NodeData,
class Real>
152 offset[0]=(int(off[0])+1)&(~(1<<depth));
153 offset[1]=(int(off[1])+1)&(~(1<<depth));
154 offset[2]=(int(off[2])+1)&(~(1<<depth));
155 width=
Real(1.0/(1<<depth));
156 for(
int dim=0;dim<DIMENSION;dim++){center.
coords[dim]=
Real(0.5+offset[dim])*width;}
159 template<
class NodeData ,
class Real >
164 centerAndWidth( c , w );
166 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);
169 template <
class NodeData,
class Real>
172 depth=index&DepthMask;
173 offset[0]=(int((index>>OffsetShift1)&OffsetMask)+1)&(~(1<<depth));
174 offset[1]=(int((index>>OffsetShift2)&OffsetMask)+1)&(~(1<<depth));
175 offset[2]=(int((index>>OffsetShift3)&OffsetMask)+1)&(~(1<<depth));
176 width=
Real(1.0/(1<<depth));
177 for(
int dim=0;dim<DIMENSION;dim++){center.
coords[dim]=
Real(0.5+offset[dim])*width;}
180 template <
class NodeData,
class Real>
182 if(!children){
return 0;}
186 d=children[i].maxDepth();
193 template <
class NodeData,
class Real>
195 if(!children){
return 1;}
203 template <
class NodeData,
class Real>
205 if(!children){
return 1;}
213 template<
class NodeData,
class Real>
215 if(depth()>maxDepth){
return 0;}
216 if(!children){
return 1;}
219 for(
int i=0;i<
Cube::CORNERS;i++){c+=children[i].maxDepthLeaves(maxDepth);}
224 template <
class NodeData,
class Real>
231 template <
class NodeData,
class Real>
234 if( !current->
parent || current==
this )
return NULL;
236 else return current+1;
239 template <
class NodeData,
class Real>
241 if(!current->
parent || current==
this){
return NULL;}
243 else{
return current+1;}
246 template<
class NodeData ,
class Real >
249 if( !current->
parent || current==
this )
return NULL;
251 else return current-1;
254 template<
class NodeData ,
class Real >
257 if( !current->
parent || current==
this )
return NULL;
259 else return current-1;
262 template <
class NodeData,
class Real>
270 const OctNode* temp=nextBranch(current);
271 if(!temp){
return NULL;}
275 template <
class NodeData,
class Real>
283 OctNode* temp=nextBranch(current);
284 if(!temp){
return NULL;}
288 template <
class NodeData,
class Real>
291 if( !current )
return this;
293 else return nextBranch(current);
296 template <
class NodeData,
class Real>
299 if( !current )
return this;
301 else return nextBranch( current );
304 template <
class NodeData,
class Real>
308 centerAndWidth(center,width);
309 for(
int dim=0;dim<DIMENSION;dim++){
310 printf(
"%[%f,%f]",center.
coords[dim]-width/2,center.
coords[dim]+width/2);
311 if(dim<DIMENSION-1){printf(
"x");}
316 template <
class NodeData,
class Real>
319 template <
class NodeData,
class Real>
320 template<
class NodeAdjacencyFunction>
322 if(processCurrent){F->Function(
this,node);}
323 if(children){__processNodeNodes(node,F);}
326 template <
class NodeData,
class Real>
327 template<
class NodeAdjacencyFunction>
329 if(processCurrent){F->Function(
this,node);}
333 __processNodeFaces(node,F,c1,c2,c3,c4);
337 template <
class NodeData,
class Real>
338 template<
class NodeAdjacencyFunction>
340 if(processCurrent){F->Function(
this,node);}
344 __processNodeEdges(node,F,c1,c2);
348 template <
class NodeData,
class Real>
349 template<
class NodeAdjacencyFunction>
351 if(processCurrent){F->Function(
this,node);}
355 F->Function(temp,node);
359 template <
class NodeData,
class Real>
360 template<
class NodeAdjacencyFunction>
362 F->Function(&children[0],node);
363 F->Function(&children[1],node);
364 F->Function(&children[2],node);
365 F->Function(&children[3],node);
366 F->Function(&children[4],node);
367 F->Function(&children[5],node);
368 F->Function(&children[6],node);
369 F->Function(&children[7],node);
370 if(children[0].children){children[0].__processNodeNodes(node,F);}
371 if(children[1].children){children[1].__processNodeNodes(node,F);}
372 if(children[2].children){children[2].__processNodeNodes(node,F);}
373 if(children[3].children){children[3].__processNodeNodes(node,F);}
374 if(children[4].children){children[4].__processNodeNodes(node,F);}
375 if(children[5].children){children[5].__processNodeNodes(node,F);}
376 if(children[6].children){children[6].__processNodeNodes(node,F);}
377 if(children[7].children){children[7].__processNodeNodes(node,F);}
380 template <
class NodeData,
class Real>
381 template<
class NodeAdjacencyFunction>
382 void OctNode<NodeData,Real>::__processNodeEdges(OctNode* node,NodeAdjacencyFunction* F,
int cIndex1,
int cIndex2){
383 F->Function(&children[cIndex1],node);
384 F->Function(&children[cIndex2],node);
385 if(children[cIndex1].children){children[cIndex1].__processNodeEdges(node,F,cIndex1,cIndex2);}
386 if(children[cIndex2].children){children[cIndex2].__processNodeEdges(node,F,cIndex1,cIndex2);}
389 template <
class NodeData,
class Real>
390 template<
class NodeAdjacencyFunction>
391 void OctNode<NodeData,Real>::__processNodeFaces(OctNode* node,NodeAdjacencyFunction* F,
int cIndex1,
int cIndex2,
int cIndex3,
int cIndex4){
392 F->Function(&children[cIndex1],node);
393 F->Function(&children[cIndex2],node);
394 F->Function(&children[cIndex3],node);
395 F->Function(&children[cIndex4],node);
396 if(children[cIndex1].children){children[cIndex1].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
397 if(children[cIndex2].children){children[cIndex2].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
398 if(children[cIndex3].children){children[cIndex3].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
399 if(children[cIndex4].children){children[cIndex4].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
402 template<
class NodeData,
class Real>
403 template<
class NodeAdjacencyFunction>
405 int c1[3],c2[3],w1,w2;
408 w1=node1->
width(maxDepth+1);
409 w2=node2->
width(maxDepth+1);
411 ProcessNodeAdjacentNodes(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 NodeAdjacencyFunction>
419 NodeAdjacencyFunction* F,
int processCurrent){
420 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
421 if(processCurrent){F->Function(node2,node1);}
423 __ProcessNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,F);
426 template<
class NodeData,
class Real>
427 template<
class TerminatingNodeAdjacencyFunction>
429 int c1[3],c2[3],w1,w2;
432 w1=node1->
width(maxDepth+1);
433 w2=node2->
width(maxDepth+1);
435 ProcessTerminatingNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,F,processCurrent);
438 template<
class NodeData,
class Real>
439 template<
class TerminatingNodeAdjacencyFunction>
443 TerminatingNodeAdjacencyFunction* F,
int processCurrent)
445 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
446 if(processCurrent){F->Function(node2,node1);}
448 __ProcessTerminatingNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,F);
451 template<
class NodeData,
class Real>
452 template<
class Po
intAdjacencyFunction>
457 w2 = node2->
width( maxDepth+1 );
458 ProcessPointAdjacentNodes( c1[0]-c2[0] , c1[1]-c2[1] , c1[2]-c2[2] , node2 , (width2*w2)>>1 , w2 , F , processCurrent );
461 template<
class NodeData,
class Real>
462 template<
class Po
intAdjacencyFunction>
465 PointAdjacencyFunction* F,
int processCurrent)
467 if( !Overlap(dx,dy,dz,radius2) )
return;
468 if( processCurrent ) F->Function(node2);
470 __ProcessPointAdjacentNodes( -dx , -dy , -dz , node2 , radius2 , width2>>1 , F );
473 template<
class NodeData,
class Real>
474 template<
class NodeAdjacencyFunction>
478 int depth,NodeAdjacencyFunction* F,
int processCurrent)
480 int c1[3],c2[3],w1,w2;
483 w1=node1->
width(maxDepth+1);
484 w2=node2->
width(maxDepth+1);
486 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);
489 template<
class NodeData,
class Real>
490 template<
class NodeAdjacencyFunction>
494 int depth,NodeAdjacencyFunction* F,
int processCurrent)
496 int d=node2->
depth();
498 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
499 if(d==depth){
if(processCurrent){F->Function(node2,node1);}}
502 __ProcessFixedDepthNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,depth-1,F);
506 template<
class NodeData,
class Real>
507 template<
class NodeAdjacencyFunction>
511 int depth,NodeAdjacencyFunction* F,
int processCurrent)
513 int c1[3],c2[3],w1,w2;
516 w1=node1->
width(maxDepth+1);
517 w2=node2->
width(maxDepth+1);
518 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);
521 template<
class NodeData,
class Real>
522 template<
class NodeAdjacencyFunction>
526 int depth,NodeAdjacencyFunction* F,
int processCurrent)
528 int d=node2->
depth();
530 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
531 if(processCurrent){F->Function(node2,node1);}
532 if(d<depth && node2->children){__ProcessMaxDepthNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2>>1,depth-1,F);}
535 template <
class NodeData,
class Real>
536 template<
class NodeAdjacencyFunction>
539 OctNode* node2,
int radius2,
int cWidth2,
540 NodeAdjacencyFunction* F)
542 int cWidth=cWidth2>>1;
543 int radius=radius2>>1;
544 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
552 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);}}
553 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);}}
554 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);}}
555 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);}}
556 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);}}
557 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);}}
558 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);}}
559 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);}}
563 template <
class NodeData,
class Real>
564 template<
class TerminatingNodeAdjacencyFunction>
565 void OctNode<NodeData,Real>::__ProcessTerminatingNodeAdjacentNodes(
int dx,
int dy,
int dz,
566 OctNode* node1,
int radius1,
567 OctNode* node2,
int radius2,
int cWidth2,
568 TerminatingNodeAdjacencyFunction* F)
570 int cWidth=cWidth2>>1;
571 int radius=radius2>>1;
572 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
580 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);}}
581 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);}}
582 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);}}
583 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);}}
584 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);}}
585 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);}}
586 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);}}
587 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);}}
591 template <
class NodeData,
class Real>
592 template<
class Po
intAdjacencyFunction>
593 void OctNode<NodeData,Real>::__ProcessPointAdjacentNodes(
int dx,
int dy,
int dz,
594 OctNode* node2,
int radius2,
int cWidth2,
595 PointAdjacencyFunction* F)
597 int cWidth=cWidth2>>1;
598 int radius=radius2>>1;
599 int o=ChildOverlap(dx,dy,dz,radius,cWidth);
608 if(o& 1){F->Function(&node2->children[0]);
if(node2->children[0].children){__ProcessPointAdjacentNodes(dx1,dy1,dz1,&node2->children[0],radius,cWidth,F);}}
609 if(o& 2){F->Function(&node2->children[1]);
if(node2->children[1].children){__ProcessPointAdjacentNodes(dx2,dy1,dz1,&node2->children[1],radius,cWidth,F);}}
610 if(o& 4){F->Function(&node2->children[2]);
if(node2->children[2].children){__ProcessPointAdjacentNodes(dx1,dy2,dz1,&node2->children[2],radius,cWidth,F);}}
611 if(o& 8){F->Function(&node2->children[3]);
if(node2->children[3].children){__ProcessPointAdjacentNodes(dx2,dy2,dz1,&node2->children[3],radius,cWidth,F);}}
612 if(o& 16){F->Function(&node2->children[4]);
if(node2->children[4].children){__ProcessPointAdjacentNodes(dx1,dy1,dz2,&node2->children[4],radius,cWidth,F);}}
613 if(o& 32){F->Function(&node2->children[5]);
if(node2->children[5].children){__ProcessPointAdjacentNodes(dx2,dy1,dz2,&node2->children[5],radius,cWidth,F);}}
614 if(o& 64){F->Function(&node2->children[6]);
if(node2->children[6].children){__ProcessPointAdjacentNodes(dx1,dy2,dz2,&node2->children[6],radius,cWidth,F);}}
615 if(o&128){F->Function(&node2->children[7]);
if(node2->children[7].children){__ProcessPointAdjacentNodes(dx2,dy2,dz2,&node2->children[7],radius,cWidth,F);}}
619 template <
class NodeData,
class Real>
620 template<
class NodeAdjacencyFunction>
621 void OctNode<NodeData,Real>::__ProcessFixedDepthNodeAdjacentNodes(
int dx,
int dy,
int dz,
622 OctNode* node1,
int radius1,
623 OctNode* node2,
int radius2,
int cWidth2,
624 int depth,NodeAdjacencyFunction* F)
626 int cWidth=cWidth2>>1;
627 int radius=radius2>>1;
628 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
636 if(node2->depth()==depth){
637 if(o& 1){F->Function(&node2->children[0],node1);}
638 if(o& 2){F->Function(&node2->children[1],node1);}
639 if(o& 4){F->Function(&node2->children[2],node1);}
640 if(o& 8){F->Function(&node2->children[3],node1);}
641 if(o& 16){F->Function(&node2->children[4],node1);}
642 if(o& 32){F->Function(&node2->children[5],node1);}
643 if(o& 64){F->Function(&node2->children[6],node1);}
644 if(o&128){F->Function(&node2->children[7],node1);}
647 if(o& 1){
if(node2->children[0].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->children[0],radius,cWidth,depth,F);}}
648 if(o& 2){
if(node2->children[1].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->children[1],radius,cWidth,depth,F);}}
649 if(o& 4){
if(node2->children[2].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->children[2],radius,cWidth,depth,F);}}
650 if(o& 8){
if(node2->children[3].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->children[3],radius,cWidth,depth,F);}}
651 if(o& 16){
if(node2->children[4].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->children[4],radius,cWidth,depth,F);}}
652 if(o& 32){
if(node2->children[5].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->children[5],radius,cWidth,depth,F);}}
653 if(o& 64){
if(node2->children[6].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->children[6],radius,cWidth,depth,F);}}
654 if(o&128){
if(node2->children[7].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->children[7],radius,cWidth,depth,F);}}
659 template <
class NodeData,
class Real>
660 template<
class NodeAdjacencyFunction>
661 void OctNode<NodeData,Real>::__ProcessMaxDepthNodeAdjacentNodes(
int dx,
int dy,
int dz,
662 OctNode* node1,
int radius1,
663 OctNode* node2,
int radius2,
int cWidth2,
664 int depth,NodeAdjacencyFunction* F)
666 int cWidth=cWidth2>>1;
667 int radius=radius2>>1;
668 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
676 if(node2->depth()<=depth){
677 if(o& 1){F->Function(&node2->children[0],node1);}
678 if(o& 2){F->Function(&node2->children[1],node1);}
679 if(o& 4){F->Function(&node2->children[2],node1);}
680 if(o& 8){F->Function(&node2->children[3],node1);}
681 if(o& 16){F->Function(&node2->children[4],node1);}
682 if(o& 32){F->Function(&node2->children[5],node1);}
683 if(o& 64){F->Function(&node2->children[6],node1);}
684 if(o&128){F->Function(&node2->children[7],node1);}
686 if(node2->depth()<depth){
687 if(o& 1){
if(node2->children[0].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->children[0],radius,cWidth,depth,F);}}
688 if(o& 2){
if(node2->children[1].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->children[1],radius,cWidth,depth,F);}}
689 if(o& 4){
if(node2->children[2].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->children[2],radius,cWidth,depth,F);}}
690 if(o& 8){
if(node2->children[3].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->children[3],radius,cWidth,depth,F);}}
691 if(o& 16){
if(node2->children[4].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->children[4],radius,cWidth,depth,F);}}
692 if(o& 32){
if(node2->children[5].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->children[5],radius,cWidth,depth,F);}}
693 if(o& 64){
if(node2->children[6].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->children[6],radius,cWidth,depth,F);}}
694 if(o&128){
if(node2->children[7].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->children[7],radius,cWidth,depth,F);}}
699 template <
class NodeData,
class Real>
700 inline int OctNode<NodeData,Real>::ChildOverlap(
int dx,
int dy,
int dz,
int d,
int cRadius2)
707 if(dx<w2 && dx>-w1){test =1;}
708 if(dx<w1 && dx>-w2){test|=2;}
711 if(dz<w2 && dz>-w1){test1 =test;}
712 if(dz<w1 && dz>-w2){test1|=test<<4;}
714 if(!test1){
return 0;}
715 if(dy<w2 && dy>-w1){overlap =test1;}
716 if(dy<w1 && dy>-w2){overlap|=test1<<2;}
720 template <
class NodeData,
class Real>
726 if(!children){
return this;}
727 centerAndWidth(center,width);
730 cIndex=CornerIndex(center,p);
733 if(cIndex&1){center.
coords[0]+=width/2;}
734 else {center.
coords[0]-=width/2;}
735 if(cIndex&2){center.
coords[1]+=width/2;}
736 else {center.
coords[1]-=width/2;}
737 if(cIndex&4){center.
coords[2]+=width/2;}
738 else {center.
coords[2]-=width/2;}
743 template <
class NodeData,
class Real>
747 if(!children){
return this;}
750 if(!i || temp<dist2){
755 return children[nearest].getNearestLeaf(p);
758 template <
class NodeData,
class Real>
760 int o1,o2,i1,i2,j1,j2;
764 if(o1!=o2){
return 0;}
770 case 0: dir[0]=1; dir[1]=2;
break;
771 case 1: dir[0]=0; dir[1]=2;
break;
772 case 2: dir[0]=0; dir[1]=1;
break;
774 int d1,d2,off1[3],off2[3];
777 idx1[0]=off1[dir[0]]+(1<<d1)+i1;
778 idx1[1]=off1[dir[1]]+(1<<d1)+j1;
779 idx2[0]=off2[dir[0]]+(1<<d2)+i2;
780 idx2[1]=off2[dir[1]]+(1<<d2)+j2;
789 if(idx1[0]==idx2[0] && idx1[1]==idx2[1]){
return 1;}
793 template<
class NodeData,
class Real>
802 template <
class NodeData,
class Real>
803 template<
class NodeData2>
810 for(i=0;i<DIMENSION;i++){this->offset[i] = node.offset[i];}
818 template <
class NodeData,
class Real>
823 template<
class NodeData ,
class Real >
828 if( n1->
d!=n2->
d )
return int(n1->
d)-int(n2->
d);
829 else if( n1->
off[0]!=n2->
off[0] )
return int(n1->
off[0]) - int(n2->
off[0]);
830 else if( n1->
off[1]!=n2->
off[1] )
return int(n1->
off[1]) - int(n2->
off[1]);
831 else if( n1->
off[2]!=n2->
off[2] )
return int(n1->
off[2]) - int(n2->
off[2]);
838 long long _p[3] = {p[0],p[1],p[2]};
839 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) );
843 template <
class NodeData,
class Real>
848 int d1 , off1[3] , d2 , off2[3];
850 if ( d1>d2 )
return 1;
851 else if( d1<d2 )
return -1;
853 if ( k1>k2 )
return 1;
854 else if( k1<k2 )
return -1;
858 template <
class NodeData,
class Real>
863 if(n1->
d!=n2->
d){
return int(n1->
d)-int(n2->
d);}
869 if(n1->
off[0]!=n2->
off[0]){
return int(n1->
off[0])-int(n2->
off[0]);}
870 if(n1->
off[1]!=n2->
off[1]){
return int(n1->
off[1])-int(n2->
off[1]);}
871 return int(n1->
off[2])-int(n2->
off[2]);
875 template <
class NodeData,
class Real>
880 template <
class NodeData,
class Real>
885 template <
class NodeData,
class Real>
888 Real w=multiplier2+multiplier1*(1<<d);
891 fabs(
Real(offSet2[0]-(offSet1[0]<<d))-w2)>=w ||
892 fabs(
Real(offSet2[1]-(offSet1[1]<<d))-w2)>=w ||
893 fabs(
Real(offSet2[2]-(offSet1[2]<<d))-w2)>=w
898 template <
class NodeData,
class Real>
900 if(c1>=dWidth || c1<=-dWidth || c2>=dWidth || c2<=-dWidth || c3>=dWidth || c3<=-dWidth){
return 0;}
904 template <
class NodeData,
class Real>
906 template <
class NodeData,
class Real>
908 template <
class NodeData,
class Real>
910 if(!parent){
return NULL;}
911 int pIndex=int(
this-(parent->children));
913 if((pIndex & (1<<dir))==(off<<dir)){
return &parent->children[pIndex];}
915 OctNode* temp=parent->__faceNeighbor(dir,off,forceChildren);
916 if(!temp){
return NULL;}
918 if(forceChildren){temp->initChildren();}
921 return &temp->children[pIndex];
925 template <
class NodeData,
class Real>
926 const OctNode<NodeData,Real>* OctNode<NodeData,Real>::__faceNeighbor(
int dir,
int off)
const {
927 if(!parent){
return NULL;}
928 int pIndex=int(
this-(parent->children));
930 if((pIndex & (1<<dir))==(off<<dir)){
return &parent->children[pIndex];}
932 const OctNode* temp=parent->__faceNeighbor(dir,off);
933 if(!temp || !temp->children){
return temp;}
934 else{
return &temp->children[pIndex];}
938 template <
class NodeData,
class Real>
943 case 0: idx[0]=1; idx[1]=2;
break;
944 case 1: idx[0]=0; idx[1]=2;
break;
945 case 2: idx[0]=0; idx[1]=1;
break;
947 return __edgeNeighbor(o,i,idx,forceChildren);
950 template <
class NodeData,
class Real>
955 case 0: idx[0]=1; idx[1]=2;
break;
956 case 1: idx[0]=0; idx[1]=2;
break;
957 case 2: idx[0]=0; idx[1]=1;
break;
959 return __edgeNeighbor(o,i,idx);
962 template <
class NodeData,
class Real>
964 if(!parent){
return NULL;}
965 int pIndex=int(
this-(parent->children));
966 int aIndex,x[DIMENSION];
969 aIndex=(~((i[0] ^ x[idx[0]]) | ((i[1] ^ x[idx[1]])<<1))) & 3;
970 pIndex^=(7 ^ (1<<o));
972 const OctNode* temp=parent->__faceNeighbor(idx[0],i[0]);
973 if(!temp || !temp->children){
return NULL;}
974 else{
return &temp->children[pIndex];}
977 const OctNode* temp=parent->__faceNeighbor(idx[1],i[1]);
978 if(!temp || !temp->children){
return NULL;}
979 else{
return &temp->children[pIndex];}
982 return &parent->children[pIndex];
985 const OctNode* temp=parent->__edgeNeighbor(o,i,idx);
986 if(!temp || !temp->children){
return temp;}
987 else{
return &temp->children[pIndex];}
992 template <
class NodeData,
class Real>
993 OctNode<NodeData,Real>* OctNode<NodeData,Real>::__edgeNeighbor(
int o,
const int i[2],
const int idx[2],
int forceChildren){
994 if(!parent){
return NULL;}
995 int pIndex=int(
this-(parent->children));
996 int aIndex,x[DIMENSION];
999 aIndex=(~((i[0] ^ x[idx[0]]) | ((i[1] ^ x[idx[1]])<<1))) & 3;
1000 pIndex^=(7 ^ (1<<o));
1002 OctNode* temp=parent->__faceNeighbor(idx[0],i[0],0);
1003 if(!temp || !temp->children){
return NULL;}
1004 else{
return &temp->children[pIndex];}
1006 else if(aIndex==2) {
1007 OctNode* temp=parent->__faceNeighbor(idx[1],i[1],0);
1008 if(!temp || !temp->children){
return NULL;}
1009 else{
return &temp->children[pIndex];}
1011 else if(aIndex==0) {
1012 return &parent->children[pIndex];
1014 else if(aIndex==3) {
1015 OctNode* temp=parent->__edgeNeighbor(o,i,idx,forceChildren);
1016 if(!temp){
return NULL;}
1017 if(!temp->children){
1018 if(forceChildren){temp->initChildren();}
1021 return &temp->children[pIndex];
1026 template <
class NodeData,
class Real>
1028 int pIndex,aIndex=0;
1029 if(!parent){
return NULL;}
1031 pIndex=int(
this-(parent->children));
1032 aIndex=(cornerIndex ^ pIndex);
1035 return &parent->children[pIndex];
1039 if(!temp || !temp->
children){
return temp;}
1040 else{
return &temp->
children[pIndex];}
1043 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(0,cornerIndex & 1);
1044 if(!temp || !temp->
children){
return NULL;}
1045 else{
return & temp->
children[pIndex];}
1048 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(1,(cornerIndex & 2)>>1);
1049 if(!temp || !temp->
children){
return NULL;}
1050 else{
return & temp->
children[pIndex];}
1053 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(2,(cornerIndex & 4)>>2);
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];}
1075 template <
class NodeData,
class Real>
1077 int pIndex,aIndex=0;
1078 if(!parent){
return NULL;}
1080 pIndex=int(
this-(parent->children));
1081 aIndex=(cornerIndex ^ pIndex);
1084 return &parent->children[pIndex];
1088 if(!temp){
return NULL;}
1096 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(0,cornerIndex & 1,0);
1097 if(!temp || !temp->
children){
return NULL;}
1098 else{
return & temp->
children[pIndex];}
1101 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(1,(cornerIndex & 2)>>1,0);
1102 if(!temp || !temp->
children){
return NULL;}
1103 else{
return & temp->
children[pIndex];}
1106 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(2,(cornerIndex & 4)>>2,0);
1107 if(!temp || !temp->
children){
return NULL;}
1108 else{
return & temp->
children[pIndex];}
1112 if(!temp || !temp->
children){
return NULL;}
1113 else{
return & temp->
children[pIndex];}
1117 if(!temp || !temp->
children){
return NULL;}
1118 else{
return & temp->
children[pIndex];}
1122 if(!temp || !temp->
children){
return NULL;}
1123 else{
return & temp->
children[pIndex];}
1131 template<
class NodeData,
class Real>
1134 template<
class NodeData,
class Real>
1136 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;}}}
1139 template<
class NodeData,
class Real>
1142 template<
class NodeData,
class Real>
1149 template<
class NodeData,
class Real>
1158 template<
class NodeData ,
class Real >
1161 if( !neighbors[
d].neighbors[1][1][1] || !neighbors[
d].neighbors[1][1][1]->
isInside( p ) )
1165 if( !
d ) neighbors[
d].neighbors[1][1][1] =
root;
1170 int i , j , k , x1 , y1 , z1 , x2 , y2 , z2;
1179 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1224 i=x1<<1 , j=y1<<1 , k=z1<<1;
1232 return neighbors[
d];
1235 template<
class NodeData ,
class Real >
1238 if( !neighbors[
d].neighbors[1][1][1] || !neighbors[
d].neighbors[1][1][1]->
isInside( p ) )
1242 if( !
d ) neighbors[
d].neighbors[1][1][1] =
root;
1245 Neighbors3& temp = getNeighbors(
root , p ,
d-1 );
1247 int i , j , k , x1 , y1 , z1 , x2 , y2 , z2;
1250 temp.neighbors[1][1][1]->centerAndWidth( c , w );
1255 if( !temp.neighbors[1][1][1] || !temp.neighbors[1][1][1]->children )
1259 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1260 neighbors[
d].neighbors[x2+i][y2+j][z2+k] = &temp.neighbors[1][1][1]->children[
Cube::CornerIndex(i,j,k)];
1265 if( temp.neighbors[i][1][1] )
1266 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)];
1268 if( temp.neighbors[1][j][1] )
1269 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)];
1271 if( temp.neighbors[1][1][k] )
1272 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)];
1276 if( temp.neighbors[i][j][1] )
1277 for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][j][z2+k] = &temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];
1279 if( temp.neighbors[i][1][k] )
1280 for( j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[i][y2+j][k] = &temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];
1282 if( temp.neighbors[1][j][k] )
1283 for( i=0 ; i<2 ; i++ ) neighbors[
d].neighbors[x2+i][j][k] = &temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];
1286 i=x1<<1 , j=y1<<1 , k=z1<<1;
1287 if( temp.neighbors[i][j][k] )
1288 neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1291 return neighbors[
d];
1294 template<
class NodeData ,
class Real >
1297 int d = node->depth();
1298 if( node==neighbors[
d].neighbors[1][1][1] )
1301 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;
1302 if( reset ) neighbors[
d].neighbors[1][1][1] = NULL;
1304 if( node!=neighbors[
d].neighbors[1][1][1] )
1306 neighbors[
d].clear();
1308 if( !node->parent ) neighbors[
d].neighbors[1][1][1] = node;
1311 int i,j,k,x1,y1,z1,x2,y2,z2;
1312 int idx=int(node-node->parent->children);
1318 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1322 Neighbors3& temp=setNeighbors(node->parent);
1326 if(temp.neighbors[i][1][1]){
1327 if(!temp.neighbors[i][1][1]->children){temp.neighbors[i][1][1]->initChildren();}
1328 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)];}}
1331 if(temp.neighbors[1][j][1]){
1332 if(!temp.neighbors[1][j][1]->children){temp.neighbors[1][j][1]->initChildren();}
1333 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)];}}
1336 if(temp.neighbors[1][1][k]){
1337 if(!temp.neighbors[1][1][k]->children){temp.neighbors[1][1][k]->initChildren();}
1338 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)];}}
1343 if(temp.neighbors[i][j][1]){
1344 if(!temp.neighbors[i][j][1]->children){temp.neighbors[i][j][1]->initChildren();}
1345 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1348 if(temp.neighbors[i][1][k]){
1349 if(!temp.neighbors[i][1][k]->children){temp.neighbors[i][1][k]->initChildren();}
1350 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1353 if(temp.neighbors[1][j][k]){
1354 if(!temp.neighbors[1][j][k]->children){temp.neighbors[1][j][k]->initChildren();}
1355 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1359 i=x1<<1; j=y1<<1; k=z1<<1;
1360 if(temp.neighbors[i][j][k]){
1361 if(!temp.neighbors[i][j][k]->children){temp.neighbors[i][j][k]->initChildren();}
1362 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1366 return neighbors[
d];
1371 template<
class NodeData ,
class Real >
1374 int d = node->depth();
1375 if( node==neighbors[
d].neighbors[1][1][1] )
1378 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;
1379 if( reset ) neighbors[
d].neighbors[1][1][1] = NULL;
1381 if( node!=neighbors[
d].neighbors[1][1][1] )
1383 neighbors[
d].clear();
1385 if( !node->parent ) neighbors[
d].neighbors[1][1][1] = node;
1388 int x1,y1,z1,x2,y2,z2;
1389 int idx=int(node-node->parent->children);
1392 for(
int i=0 ; i<2 ; i++ )
1393 for(
int j=0 ; j<2 ; j++ )
1394 for(
int k=0 ; k<2 ; k++ )
1395 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1397 Neighbors3& temp=setNeighbors( node->parent , flags );
1402 if( temp.neighbors[i][1][1] )
1404 if( flags[i][1][1] && !temp.neighbors[i][1][1]->children ) temp.neighbors[i][1][1]->initChildren();
1405 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)];
1410 if( temp.neighbors[1][j][1] )
1412 if( flags[1][j][1] && !temp.neighbors[1][j][1]->children ) temp.neighbors[1][j][1]->initChildren();
1413 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)];
1418 if( temp.neighbors[1][1][k] )
1420 if( flags[1][1][k] && !temp.neighbors[1][1][k]->children ) temp.neighbors[1][1][k]->initChildren();
1421 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)];
1427 int i=x1<<1 , j=y1<<1;
1428 if( temp.neighbors[i][j][1] )
1430 if( flags[i][j][1] && !temp.neighbors[i][j][1]->children ) temp.neighbors[i][j][1]->initChildren();
1431 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)];
1435 int i=x1<<1 , k=z1<<1;
1436 if( temp.neighbors[i][1][k] )
1438 if( flags[i][1][k] && !temp.neighbors[i][1][k]->children ) temp.neighbors[i][1][k]->initChildren();
1439 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)];
1443 int j=y1<<1 , k=z1<<1;
1444 if( temp.neighbors[1][j][k] )
1446 if( flags[1][j][k] && !temp.neighbors[1][j][k]->children ) temp.neighbors[1][j][k]->initChildren();
1447 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)];
1453 int i=x1<<1 , j=y1<<1 , k=z1<<1;
1454 if( temp.neighbors[i][j][k] )
1456 if( flags[i][j][k] && !temp.neighbors[i][j][k]->children ) temp.neighbors[i][j][k]->initChildren();
1457 if( temp.neighbors[i][j][k]->children ) neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1462 return neighbors[
d];
1465 template<
class NodeData,
class Real>
1467 int d=node->depth();
1468 if(node!=neighbors[
d].neighbors[1][1][1]){
1469 neighbors[
d].clear();
1471 if(!node->parent){neighbors[
d].neighbors[1][1][1]=node;}
1473 int i,j,k,x1,y1,z1,x2,y2,z2;
1474 int idx=int(node-node->parent->children);
1480 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1484 Neighbors3& temp=getNeighbors(node->parent);
1488 if(temp.neighbors[i][1][1] && temp.neighbors[i][1][1]->children){
1489 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)];}}
1492 if(temp.neighbors[1][j][1] && temp.neighbors[1][j][1]->children){
1493 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)];}}
1496 if(temp.neighbors[1][1][k] && temp.neighbors[1][1][k]->children){
1497 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)];}}
1502 if(temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children){
1503 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1506 if(temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children){
1507 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1510 if(temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children){
1511 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1515 i=x1<<1; j=y1<<1; k=z1<<1;
1516 if(temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children){
1517 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1521 return neighbors[node->depth()];
1527 template<
class NodeData,
class Real>
1530 template<
class NodeData,
class Real>
1532 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;}}}
1535 template<
class NodeData,
class Real>
1538 template<
class NodeData,
class Real>
1544 template<
class NodeData,
class Real>
1552 template<
class NodeData,
class Real>
1555 if(node!=neighbors[
d].neighbors[1][1][1]){
1556 neighbors[
d].clear();
1558 if(!node->
parent){neighbors[
d].neighbors[1][1][1]=node;}
1560 int i,j,k,x1,y1,z1,x2,y2,z2;
1561 int idx=int(node-node->
parent->children);
1571 ConstNeighbors3& temp=getNeighbors(node->
parent);
1575 if(temp.neighbors[i][1][1] && temp.neighbors[i][1][1]->children){
1576 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)];}}
1579 if(temp.neighbors[1][j][1] && temp.neighbors[1][j][1]->children){
1580 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)];}}
1583 if(temp.neighbors[1][1][k] && temp.neighbors[1][1][k]->children){
1584 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)];}}
1589 if(temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children){
1590 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1593 if(temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children){
1594 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1597 if(temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children){
1598 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1602 i=x1<<1; j=y1<<1; k=z1<<1;
1603 if(temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children){
1604 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1608 return neighbors[node->
depth()];
1611 template<
class NodeData,
class Real>
1614 int d=node->depth();
1619 if( node!=neighbors[
d].neighbors[1][1][1] )
1621 neighbors[
d].clear();
1623 if(
d==minDepth ) neighbors[
d].neighbors[1][1][1]=node;
1626 int i,j,k,x1,y1,z1,x2,y2,z2;
1627 int idx = int(node-node->parent->children);
1631 ConstNeighbors3& temp=getNeighbors( node->parent , minDepth );
1634 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1635 neighbors[
d].neighbors[x2+i][y2+j][z2+k] = &node->parent->children[
Cube::CornerIndex(i,j,k) ];
1639 if( temp.neighbors[i][1][1] && temp.neighbors[i][1][1]->children )
1640 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)];
1643 if( temp.neighbors[1][j][1] && temp.neighbors[1][j][1]->children )
1644 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)];
1647 if( temp.neighbors[1][1][k] && temp.neighbors[1][1][k]->children )
1648 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)];
1652 if( temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children )
1653 for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][j][z2+k] = &temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];
1656 if( temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children )
1657 for( j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[i][y2+j][k] = &temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];
1660 if( temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children )
1661 for( i=0 ; i<2 ; i++ ) neighbors[
d].neighbors[x2+i][j][k] = &temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];
1664 i=x1<<1 , j=y1<<1 , k=z1<<1;
1665 if( temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children )
1666 neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1669 return neighbors[node->depth()];
1676 template<
class NodeData ,
class Real >
1679 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;
1682 template<
class NodeData ,
class Real >
1685 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;
1688 template<
class NodeData ,
class Real >
1695 template<
class NodeData ,
class Real >
1702 template<
class NodeData ,
class Real >
1709 template<
class NodeData ,
class Real >
1716 template<
class NodeData ,
class Real >
1726 template<
class NodeData ,
class Real >
1736 template<
class NodeData ,
class Real >
1740 if( node!=neighbors[
d].neighbors[2][2][2] )
1742 neighbors[
d].clear();
1744 if( !node->
parent ) neighbors[
d].neighbors[2][2][2]=node;
1747 getNeighbors( node->
parent );
1749 int x1 , y1 , z1 , x2 , y2 , z2;
1756 int fx0 = x2+1 , fy0 = y2+1 , fz0 = z2+1;
1757 int cx1 = x1*2+1 , cy1 = y1*2+1 , cz1 = z1*2+1;
1758 int cx2 = x2*2+1 , cy2 = y2*2+1 , cz2 = z2*2+1;
1759 int fx1 = x1*3 , fy1 = y1*3 , fz1 = z1*3;
1760 int fx2 = x2*4 , fy2 = y2*4 , fz2 = z2*4;
1763 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1768 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1771 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1774 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1777 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1780 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1783 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1788 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1791 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1794 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1797 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1800 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1803 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1806 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1809 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1812 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1815 for( k=0 ; k<2 ; k++ )
1818 for( j=0 ; j<2 ; j++ )
1821 for( i=0 ; i<2 ; i++ )
1826 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1829 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1832 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1835 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1838 for( i=0 ; i<2 ; i++ )
1841 for( j=0 ; j<2 ; j++ )
1844 for( k=0 ; k<2 ; k++ )
1850 return neighbors[
d];
1853 template<
class NodeData ,
class Real >
1857 if( node!=neighbors[
d].neighbors[2][2][2] )
1859 neighbors[
d].clear();
1861 if( !node->
parent ) neighbors[
d].neighbors[2][2][2]=node;
1864 setNeighbors( node->
parent , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1866 int x1 , y1 , z1 , x2 , y2 , z2 , ii , jj , kk;
1870 for(
int i=xStart ; i<xEnd ; i++ )
1875 for(
int j=yStart ; j<yEnd ; j++ )
1880 for(
int k=zStart ; k<zEnd ; k++ )
1895 return neighbors[
d];
1898 template<
class NodeData ,
class Real >
1902 if( node!=neighbors[
d].neighbors[2][2][2] )
1904 neighbors[
d].clear();
1906 if(!node->
parent) neighbors[
d].neighbors[2][2][2]=node;
1909 getNeighbors( node->
parent );
1911 int x1,y1,z1,x2,y2,z2,ii,jj,kk;
1915 for(
int i=0;i<5;i++)
1920 for(
int j=0;j<5;j++)
1925 for(
int k=0;k<5;k++)
1937 return neighbors[
d];
1940 template <
class NodeData,
class Real>
1942 FILE* fp=fopen(fileName,
"wb");
1949 template <
class NodeData,
class Real>
1956 template <
class NodeData,
class Real>
1958 FILE* fp=fopen(fileName,
"rb");
1965 template <
class NodeData,
class Real>
1980 template<
class NodeData,
class Real>
1986 template<
class NodeData,
class Real>