33 #include "poisson_exceptions.h"
52 template<
class NodeData,
class Real>
int OctNode<NodeData,Real>::UseAlloc=0;
55 template<
class NodeData,
class Real>
61 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){
if(children){
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);}
94 if(children){
delete[] children;}
102 depthAndOffset(d,off);
107 children[idx].parent=
this;
108 children[idx].children=NULL;
110 off2[0]=(off[0]<<1)+i;
111 off2[1]=(off[1]<<1)+j;
112 off2[2]=(off[2]<<1)+k;
113 Index(d+1,off2,children[idx].d,children[idx].off);
119 template <
class NodeData,
class Real>
122 off[0]=short((1<<depth)+offset[0]-1);
123 off[1]=short((1<<depth)+offset[1]-1);
124 off[2]=short((1<<depth)+offset[2]-1);
127 template<
class NodeData,
class Real>
130 offset[0]=(int(off[0])+1)&(~(1<<depth));
131 offset[1]=(int(off[1])+1)&(~(1<<depth));
132 offset[2]=(int(off[2])+1)&(~(1<<depth));
134 template<
class NodeData,
class Real>
136 template<
class NodeData,
class Real>
138 depth=int(index&DepthMask);
139 offset[0]=(int((index>>OffsetShift1)&OffsetMask)+1)&(~(1<<depth));
140 offset[1]=(int((index>>OffsetShift2)&OffsetMask)+1)&(~(1<<depth));
141 offset[2]=(int((index>>OffsetShift3)&OffsetMask)+1)&(~(1<<depth));
143 template<
class NodeData,
class Real>
145 template <
class NodeData,
class Real>
149 offset[0]=(int(off[0])+1)&(~(1<<depth));
150 offset[1]=(int(off[1])+1)&(~(1<<depth));
151 offset[2]=(int(off[2])+1)&(~(1<<depth));
152 width=
Real(1.0/(1<<depth));
153 for(
int dim=0;dim<DIMENSION;dim++){center.
coords[dim]=
Real(0.5+offset[dim])*width;}
155 template<
class NodeData ,
class Real >
160 centerAndWidth( c , w );
162 return (c[0]-w)<p[0] && p[0]<=(c[0]+w) && (c[1]-w)<p[1] && p[1]<=(c[1]+w) && (c[2]-w)<p[2] && p[2]<=(c[2]+w);
164 template <
class NodeData,
class Real>
167 depth=index&DepthMask;
168 offset[0]=(int((index>>OffsetShift1)&OffsetMask)+1)&(~(1<<depth));
169 offset[1]=(int((index>>OffsetShift2)&OffsetMask)+1)&(~(1<<depth));
170 offset[2]=(int((index>>OffsetShift3)&OffsetMask)+1)&(~(1<<depth));
171 width=
Real(1.0/(1<<depth));
172 for(
int dim=0;dim<DIMENSION;dim++){center.
coords[dim]=
Real(0.5+offset[dim])*width;}
175 template <
class NodeData,
class Real>
177 if(!children){
return 0;}
181 d=children[i].maxDepth();
187 template <
class NodeData,
class Real>
189 if(!children){
return 1;}
196 template <
class NodeData,
class Real>
198 if(!children){
return 1;}
205 template<
class NodeData,
class Real>
207 if(depth()>maxDepth){
return 0;}
208 if(!children){
return 1;}
211 for(
int i=0;i<
Cube::CORNERS;i++){c+=children[i].maxDepthLeaves(maxDepth);}
215 template <
class NodeData,
class Real>
223 template <
class NodeData,
class Real>
226 if( !current->
parent || current==
this )
return NULL;
228 else return current+1;
230 template <
class NodeData,
class Real>
232 if(!current->
parent || current==
this){
return NULL;}
234 else{
return current+1;}
236 template<
class NodeData ,
class Real >
239 if( !current->
parent || current==
this )
return NULL;
241 else return current-1;
243 template<
class NodeData ,
class Real >
246 if( !current->
parent || current==
this )
return NULL;
248 else return current-1;
250 template <
class NodeData,
class Real>
258 const OctNode* temp=nextBranch(current);
259 if(!temp){
return NULL;}
262 template <
class NodeData,
class Real>
270 OctNode* temp=nextBranch(current);
271 if(!temp){
return NULL;}
275 template <
class NodeData,
class Real>
278 if( !current )
return this;
280 else return nextBranch(current);
282 template <
class NodeData,
class Real>
285 if( !current )
return this;
287 else return nextBranch( current );
290 template <
class NodeData,
class Real>
294 centerAndWidth(center,width);
295 for(
int dim=0;dim<DIMENSION;dim++){
296 printf(
"%[%f,%f]",center.
coords[dim]-width/2,center.
coords[dim]+width/2);
297 if(dim<DIMENSION-1){printf(
"x");}
302 template <
class NodeData,
class Real>
305 template <
class NodeData,
class Real>
306 template<
class NodeAdjacencyFunction>
308 if(processCurrent){F->Function(
this,node);}
309 if(children){__processNodeNodes(node,F);}
311 template <
class NodeData,
class Real>
312 template<
class NodeAdjacencyFunction>
314 if(processCurrent){F->Function(
this,node);}
318 __processNodeFaces(node,F,c1,c2,c3,c4);
321 template <
class NodeData,
class Real>
322 template<
class NodeAdjacencyFunction>
324 if(processCurrent){F->Function(
this,node);}
328 __processNodeEdges(node,F,c1,c2);
331 template <
class NodeData,
class Real>
332 template<
class NodeAdjacencyFunction>
334 if(processCurrent){F->Function(
this,node);}
338 F->Function(temp,node);
341 template <
class NodeData,
class Real>
342 template<
class NodeAdjacencyFunction>
344 F->Function(&children[0],node);
345 F->Function(&children[1],node);
346 F->Function(&children[2],node);
347 F->Function(&children[3],node);
348 F->Function(&children[4],node);
349 F->Function(&children[5],node);
350 F->Function(&children[6],node);
351 F->Function(&children[7],node);
352 if(children[0].children){children[0].__processNodeNodes(node,F);}
353 if(children[1].children){children[1].__processNodeNodes(node,F);}
354 if(children[2].children){children[2].__processNodeNodes(node,F);}
355 if(children[3].children){children[3].__processNodeNodes(node,F);}
356 if(children[4].children){children[4].__processNodeNodes(node,F);}
357 if(children[5].children){children[5].__processNodeNodes(node,F);}
358 if(children[6].children){children[6].__processNodeNodes(node,F);}
359 if(children[7].children){children[7].__processNodeNodes(node,F);}
361 template <
class NodeData,
class Real>
362 template<
class NodeAdjacencyFunction>
363 void OctNode<NodeData,Real>::__processNodeEdges(OctNode* node,NodeAdjacencyFunction* F,
int cIndex1,
int cIndex2){
364 F->Function(&children[cIndex1],node);
365 F->Function(&children[cIndex2],node);
366 if(children[cIndex1].children){children[cIndex1].__processNodeEdges(node,F,cIndex1,cIndex2);}
367 if(children[cIndex2].children){children[cIndex2].__processNodeEdges(node,F,cIndex1,cIndex2);}
369 template <
class NodeData,
class Real>
370 template<
class NodeAdjacencyFunction>
371 void OctNode<NodeData,Real>::__processNodeFaces(OctNode* node,NodeAdjacencyFunction* F,
int cIndex1,
int cIndex2,
int cIndex3,
int cIndex4){
372 F->Function(&children[cIndex1],node);
373 F->Function(&children[cIndex2],node);
374 F->Function(&children[cIndex3],node);
375 F->Function(&children[cIndex4],node);
376 if(children[cIndex1].children){children[cIndex1].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
377 if(children[cIndex2].children){children[cIndex2].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
378 if(children[cIndex3].children){children[cIndex3].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
379 if(children[cIndex4].children){children[cIndex4].__processNodeFaces(node,F,cIndex1,cIndex2,cIndex3,cIndex4);}
381 template<
class NodeData,
class Real>
382 template<
class NodeAdjacencyFunction>
384 int c1[3],c2[3],w1,w2;
387 w1=node1->
width(maxDepth+1);
388 w2=node2->
width(maxDepth+1);
390 ProcessNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,F,processCurrent);
392 template<
class NodeData,
class Real>
393 template<
class NodeAdjacencyFunction>
397 NodeAdjacencyFunction* F,
int processCurrent){
398 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
399 if(processCurrent){F->Function(node2,node1);}
401 __ProcessNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,F);
403 template<
class NodeData,
class Real>
404 template<
class TerminatingNodeAdjacencyFunction>
406 int c1[3],c2[3],w1,w2;
409 w1=node1->
width(maxDepth+1);
410 w2=node2->
width(maxDepth+1);
412 ProcessTerminatingNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,F,processCurrent);
414 template<
class NodeData,
class Real>
415 template<
class TerminatingNodeAdjacencyFunction>
419 TerminatingNodeAdjacencyFunction* F,
int processCurrent)
421 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
422 if(processCurrent){F->Function(node2,node1);}
424 __ProcessTerminatingNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,F);
426 template<
class NodeData,
class Real>
427 template<
class Po
intAdjacencyFunction>
432 w2 = node2->
width( maxDepth+1 );
433 ProcessPointAdjacentNodes( c1[0]-c2[0] , c1[1]-c2[1] , c1[2]-c2[2] , node2 , (width2*w2)>>1 , w2 , F , processCurrent );
435 template<
class NodeData,
class Real>
436 template<
class Po
intAdjacencyFunction>
439 PointAdjacencyFunction* F,
int processCurrent)
441 if( !Overlap(dx,dy,dz,radius2) )
return;
442 if( processCurrent ) F->Function(node2);
444 __ProcessPointAdjacentNodes( -dx , -dy , -dz , node2 , radius2 , width2>>1 , F );
446 template<
class NodeData,
class Real>
447 template<
class NodeAdjacencyFunction>
451 int depth,NodeAdjacencyFunction* F,
int processCurrent)
453 int c1[3],c2[3],w1,w2;
456 w1=node1->
width(maxDepth+1);
457 w2=node2->
width(maxDepth+1);
459 ProcessFixedDepthNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,depth,F,processCurrent);
461 template<
class NodeData,
class Real>
462 template<
class NodeAdjacencyFunction>
466 int depth,NodeAdjacencyFunction* F,
int processCurrent)
468 int d=node2->
depth();
470 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
471 if(d==depth){
if(processCurrent){F->Function(node2,node1);}}
474 __ProcessFixedDepthNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2/2,depth-1,F);
477 template<
class NodeData,
class Real>
478 template<
class NodeAdjacencyFunction>
482 int depth,NodeAdjacencyFunction* F,
int processCurrent)
484 int c1[3],c2[3],w1,w2;
487 w1=node1->
width(maxDepth+1);
488 w2=node2->
width(maxDepth+1);
489 ProcessMaxDepthNodeAdjacentNodes(c1[0]-c2[0],c1[1]-c2[1],c1[2]-c2[2],node1,(width1*w1)>>1,node2,(width2*w2)>>1,w2,depth,F,processCurrent);
491 template<
class NodeData,
class Real>
492 template<
class NodeAdjacencyFunction>
496 int depth,NodeAdjacencyFunction* F,
int processCurrent)
498 int d=node2->
depth();
500 if(!Overlap(dx,dy,dz,radius1+radius2)){
return;}
501 if(processCurrent){F->Function(node2,node1);}
502 if(d<depth && node2->children){__ProcessMaxDepthNodeAdjacentNodes(-dx,-dy,-dz,node1,radius1,node2,radius2,width2>>1,depth-1,F);}
504 template <
class NodeData,
class Real>
505 template<
class NodeAdjacencyFunction>
508 OctNode* node2,
int radius2,
int cWidth2,
509 NodeAdjacencyFunction* F)
511 int cWidth=cWidth2>>1;
512 int radius=radius2>>1;
513 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
521 if(o& 1){F->Function(&node2->
children[0],node1);
if(node2->
children[0].
children){__ProcessNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->
children[0],radius,cWidth,F);}}
522 if(o& 2){F->Function(&node2->
children[1],node1);
if(node2->
children[1].
children){__ProcessNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->
children[1],radius,cWidth,F);}}
523 if(o& 4){F->Function(&node2->
children[2],node1);
if(node2->
children[2].
children){__ProcessNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->
children[2],radius,cWidth,F);}}
524 if(o& 8){F->Function(&node2->
children[3],node1);
if(node2->
children[3].
children){__ProcessNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->
children[3],radius,cWidth,F);}}
525 if(o& 16){F->Function(&node2->
children[4],node1);
if(node2->
children[4].
children){__ProcessNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->
children[4],radius,cWidth,F);}}
526 if(o& 32){F->Function(&node2->
children[5],node1);
if(node2->
children[5].
children){__ProcessNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->
children[5],radius,cWidth,F);}}
527 if(o& 64){F->Function(&node2->
children[6],node1);
if(node2->
children[6].
children){__ProcessNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->
children[6],radius,cWidth,F);}}
528 if(o&128){F->Function(&node2->
children[7],node1);
if(node2->
children[7].
children){__ProcessNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->
children[7],radius,cWidth,F);}}
531 template <
class NodeData,
class Real>
532 template<
class TerminatingNodeAdjacencyFunction>
533 void OctNode<NodeData,Real>::__ProcessTerminatingNodeAdjacentNodes(
int dx,
int dy,
int dz,
534 OctNode* node1,
int radius1,
535 OctNode* node2,
int radius2,
int cWidth2,
536 TerminatingNodeAdjacencyFunction* F)
538 int cWidth=cWidth2>>1;
539 int radius=radius2>>1;
540 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
548 if(o& 1){
if(F->Function(&node2->children[0],node1) && node2->children[0].children){__ProcessTerminatingNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->children[0],radius,cWidth,F);}}
549 if(o& 2){
if(F->Function(&node2->children[1],node1) && node2->children[1].children){__ProcessTerminatingNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->children[1],radius,cWidth,F);}}
550 if(o& 4){
if(F->Function(&node2->children[2],node1) && node2->children[2].children){__ProcessTerminatingNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->children[2],radius,cWidth,F);}}
551 if(o& 8){
if(F->Function(&node2->children[3],node1) && node2->children[3].children){__ProcessTerminatingNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->children[3],radius,cWidth,F);}}
552 if(o& 16){
if(F->Function(&node2->children[4],node1) && node2->children[4].children){__ProcessTerminatingNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->children[4],radius,cWidth,F);}}
553 if(o& 32){
if(F->Function(&node2->children[5],node1) && node2->children[5].children){__ProcessTerminatingNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->children[5],radius,cWidth,F);}}
554 if(o& 64){
if(F->Function(&node2->children[6],node1) && node2->children[6].children){__ProcessTerminatingNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->children[6],radius,cWidth,F);}}
555 if(o&128){
if(F->Function(&node2->children[7],node1) && node2->children[7].children){__ProcessTerminatingNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->children[7],radius,cWidth,F);}}
558 template <
class NodeData,
class Real>
559 template<
class Po
intAdjacencyFunction>
560 void OctNode<NodeData,Real>::__ProcessPointAdjacentNodes(
int dx,
int dy,
int dz,
561 OctNode* node2,
int radius2,
int cWidth2,
562 PointAdjacencyFunction* F)
564 int cWidth=cWidth2>>1;
565 int radius=radius2>>1;
566 int o=ChildOverlap(dx,dy,dz,radius,cWidth);
575 if(o& 1){F->Function(&node2->children[0]);
if(node2->children[0].children){__ProcessPointAdjacentNodes(dx1,dy1,dz1,&node2->children[0],radius,cWidth,F);}}
576 if(o& 2){F->Function(&node2->children[1]);
if(node2->children[1].children){__ProcessPointAdjacentNodes(dx2,dy1,dz1,&node2->children[1],radius,cWidth,F);}}
577 if(o& 4){F->Function(&node2->children[2]);
if(node2->children[2].children){__ProcessPointAdjacentNodes(dx1,dy2,dz1,&node2->children[2],radius,cWidth,F);}}
578 if(o& 8){F->Function(&node2->children[3]);
if(node2->children[3].children){__ProcessPointAdjacentNodes(dx2,dy2,dz1,&node2->children[3],radius,cWidth,F);}}
579 if(o& 16){F->Function(&node2->children[4]);
if(node2->children[4].children){__ProcessPointAdjacentNodes(dx1,dy1,dz2,&node2->children[4],radius,cWidth,F);}}
580 if(o& 32){F->Function(&node2->children[5]);
if(node2->children[5].children){__ProcessPointAdjacentNodes(dx2,dy1,dz2,&node2->children[5],radius,cWidth,F);}}
581 if(o& 64){F->Function(&node2->children[6]);
if(node2->children[6].children){__ProcessPointAdjacentNodes(dx1,dy2,dz2,&node2->children[6],radius,cWidth,F);}}
582 if(o&128){F->Function(&node2->children[7]);
if(node2->children[7].children){__ProcessPointAdjacentNodes(dx2,dy2,dz2,&node2->children[7],radius,cWidth,F);}}
585 template <
class NodeData,
class Real>
586 template<
class NodeAdjacencyFunction>
587 void OctNode<NodeData,Real>::__ProcessFixedDepthNodeAdjacentNodes(
int dx,
int dy,
int dz,
588 OctNode* node1,
int radius1,
589 OctNode* node2,
int radius2,
int cWidth2,
590 int depth,NodeAdjacencyFunction* F)
592 int cWidth=cWidth2>>1;
593 int radius=radius2>>1;
594 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
602 if(node2->depth()==depth){
603 if(o& 1){F->Function(&node2->children[0],node1);}
604 if(o& 2){F->Function(&node2->children[1],node1);}
605 if(o& 4){F->Function(&node2->children[2],node1);}
606 if(o& 8){F->Function(&node2->children[3],node1);}
607 if(o& 16){F->Function(&node2->children[4],node1);}
608 if(o& 32){F->Function(&node2->children[5],node1);}
609 if(o& 64){F->Function(&node2->children[6],node1);}
610 if(o&128){F->Function(&node2->children[7],node1);}
613 if(o& 1){
if(node2->children[0].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->children[0],radius,cWidth,depth,F);}}
614 if(o& 2){
if(node2->children[1].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->children[1],radius,cWidth,depth,F);}}
615 if(o& 4){
if(node2->children[2].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->children[2],radius,cWidth,depth,F);}}
616 if(o& 8){
if(node2->children[3].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->children[3],radius,cWidth,depth,F);}}
617 if(o& 16){
if(node2->children[4].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->children[4],radius,cWidth,depth,F);}}
618 if(o& 32){
if(node2->children[5].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->children[5],radius,cWidth,depth,F);}}
619 if(o& 64){
if(node2->children[6].children){__ProcessFixedDepthNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->children[6],radius,cWidth,depth,F);}}
620 if(o&128){
if(node2->children[7].children){__ProcessFixedDepthNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->children[7],radius,cWidth,depth,F);}}
624 template <
class NodeData,
class Real>
625 template<
class NodeAdjacencyFunction>
626 void OctNode<NodeData,Real>::__ProcessMaxDepthNodeAdjacentNodes(
int dx,
int dy,
int dz,
627 OctNode* node1,
int radius1,
628 OctNode* node2,
int radius2,
int cWidth2,
629 int depth,NodeAdjacencyFunction* F)
631 int cWidth=cWidth2>>1;
632 int radius=radius2>>1;
633 int o=ChildOverlap(dx,dy,dz,radius1+radius,cWidth);
641 if(node2->depth()<=depth){
642 if(o& 1){F->Function(&node2->children[0],node1);}
643 if(o& 2){F->Function(&node2->children[1],node1);}
644 if(o& 4){F->Function(&node2->children[2],node1);}
645 if(o& 8){F->Function(&node2->children[3],node1);}
646 if(o& 16){F->Function(&node2->children[4],node1);}
647 if(o& 32){F->Function(&node2->children[5],node1);}
648 if(o& 64){F->Function(&node2->children[6],node1);}
649 if(o&128){F->Function(&node2->children[7],node1);}
651 if(node2->depth()<depth){
652 if(o& 1){
if(node2->children[0].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy1,dz1,node1,radius1,&node2->children[0],radius,cWidth,depth,F);}}
653 if(o& 2){
if(node2->children[1].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy1,dz1,node1,radius1,&node2->children[1],radius,cWidth,depth,F);}}
654 if(o& 4){
if(node2->children[2].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy2,dz1,node1,radius1,&node2->children[2],radius,cWidth,depth,F);}}
655 if(o& 8){
if(node2->children[3].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy2,dz1,node1,radius1,&node2->children[3],radius,cWidth,depth,F);}}
656 if(o& 16){
if(node2->children[4].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy1,dz2,node1,radius1,&node2->children[4],radius,cWidth,depth,F);}}
657 if(o& 32){
if(node2->children[5].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy1,dz2,node1,radius1,&node2->children[5],radius,cWidth,depth,F);}}
658 if(o& 64){
if(node2->children[6].children){__ProcessMaxDepthNodeAdjacentNodes(dx1,dy2,dz2,node1,radius1,&node2->children[6],radius,cWidth,depth,F);}}
659 if(o&128){
if(node2->children[7].children){__ProcessMaxDepthNodeAdjacentNodes(dx2,dy2,dz2,node1,radius1,&node2->children[7],radius,cWidth,depth,F);}}
663 template <
class NodeData,
class Real>
664 inline int OctNode<NodeData,Real>::ChildOverlap(
int dx,
int dy,
int dz,
int d,
int cRadius2)
671 if(dx<w2 && dx>-w1){test =1;}
672 if(dx<w1 && dx>-w2){test|=2;}
675 if(dz<w2 && dz>-w1){test1 =test;}
676 if(dz<w1 && dz>-w2){test1|=test<<4;}
678 if(!test1){
return 0;}
679 if(dy<w2 && dy>-w1){overlap =test1;}
680 if(dy<w1 && dy>-w2){overlap|=test1<<2;}
684 template <
class NodeData,
class Real>
690 if(!children){
return this;}
691 centerAndWidth(center,width);
694 cIndex=CornerIndex(center,p);
697 if(cIndex&1){center.
coords[0]+=width/2;}
698 else {center.
coords[0]-=width/2;}
699 if(cIndex&2){center.
coords[1]+=width/2;}
700 else {center.
coords[1]-=width/2;}
701 if(cIndex&4){center.
coords[2]+=width/2;}
702 else {center.
coords[2]-=width/2;}
706 template <
class NodeData,
class Real>
710 if(!children){
return this;}
713 if(!i || temp<dist2){
718 return children[nearest].getNearestLeaf(p);
721 template <
class NodeData,
class Real>
723 int o1,o2,i1,i2,j1,j2;
727 if(o1!=o2){
return 0;}
733 case 0: dir[0]=1; dir[1]=2;
break;
734 case 1: dir[0]=0; dir[1]=2;
break;
735 case 2: dir[0]=0; dir[1]=1;
break;
737 int d1,d2,off1[3],off2[3];
740 idx1[0]=off1[dir[0]]+(1<<d1)+i1;
741 idx1[1]=off1[dir[1]]+(1<<d1)+j1;
742 idx2[0]=off2[dir[0]]+(1<<d2)+i2;
743 idx2[1]=off2[dir[1]]+(1<<d2)+j2;
752 if(idx1[0]==idx2[0] && idx1[1]==idx2[1]){
return 1;}
755 template<
class NodeData,
class Real>
763 template <
class NodeData,
class Real>
764 template<
class NodeData2>
767 if(children){
delete[] children;}
771 for(i=0;i<DIMENSION;i++){this->offset[i] = node.offset[i];}
778 template <
class NodeData,
class Real>
782 template<
class NodeData ,
class Real >
787 if( n1->
d!=n2->
d )
return int(n1->
d)-int(n2->
d);
788 else if( n1->
off[0]!=n2->
off[0] )
return int(n1->
off[0]) - int(n2->
off[0]);
789 else if( n1->
off[1]!=n2->
off[1] )
return int(n1->
off[1]) - int(n2->
off[1]);
790 else if( n1->
off[2]!=n2->
off[2] )
return int(n1->
off[2]) - int(n2->
off[2]);
797 long long _p[3] = {p[0],p[1],p[2]};
798 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) );
801 template <
class NodeData,
class Real>
806 int d1 , off1[3] , d2 , off2[3];
808 if ( d1>d2 )
return 1;
809 else if( d1<d2 )
return -1;
811 if ( k1>k2 )
return 1;
812 else if( k1<k2 )
return -1;
816 template <
class NodeData,
class Real>
821 if(n1->
d!=n2->
d){
return int(n1->
d)-int(n2->
d);}
827 if(n1->
off[0]!=n2->
off[0]){
return int(n1->
off[0])-int(n2->
off[0]);}
828 if(n1->
off[1]!=n2->
off[1]){
return int(n1->
off[1])-int(n2->
off[1]);}
829 return int(n1->
off[2])-int(n2->
off[2]);
832 template <
class NodeData,
class Real>
836 template <
class NodeData,
class Real>
840 template <
class NodeData,
class Real>
843 Real w=multiplier2+multiplier1*(1<<d);
846 fabs(
Real(offSet2[0]-(offSet1[0]<<d))-w2)>=w ||
847 fabs(
Real(offSet2[1]-(offSet1[1]<<d))-w2)>=w ||
848 fabs(
Real(offSet2[2]-(offSet1[2]<<d))-w2)>=w
852 template <
class NodeData,
class Real>
854 if(c1>=dWidth || c1<=-dWidth || c2>=dWidth || c2<=-dWidth || c3>=dWidth || c3<=-dWidth){
return 0;}
857 template <
class NodeData,
class Real>
859 template <
class NodeData,
class Real>
861 template <
class NodeData,
class Real>
863 if(!parent){
return NULL;}
864 int pIndex=int(
this-(parent->children));
866 if((pIndex & (1<<dir))==(off<<dir)){
return &parent->children[pIndex];}
868 OctNode* temp=parent->__faceNeighbor(dir,off,forceChildren);
869 if(!temp){
return NULL;}
871 if(forceChildren){temp->initChildren();}
874 return &temp->children[pIndex];
877 template <
class NodeData,
class Real>
878 const OctNode<NodeData,Real>* OctNode<NodeData,Real>::__faceNeighbor(
int dir,
int off)
const {
879 if(!parent){
return NULL;}
880 int pIndex=int(
this-(parent->children));
882 if((pIndex & (1<<dir))==(off<<dir)){
return &parent->children[pIndex];}
884 const OctNode* temp=parent->__faceNeighbor(dir,off);
885 if(!temp || !temp->children){
return temp;}
886 else{
return &temp->children[pIndex];}
890 template <
class NodeData,
class Real>
895 case 0: idx[0]=1; idx[1]=2;
break;
896 case 1: idx[0]=0; idx[1]=2;
break;
897 case 2: idx[0]=0; idx[1]=1;
break;
899 return __edgeNeighbor(o,i,idx,forceChildren);
901 template <
class NodeData,
class Real>
906 case 0: idx[0]=1; idx[1]=2;
break;
907 case 1: idx[0]=0; idx[1]=2;
break;
908 case 2: idx[0]=0; idx[1]=1;
break;
910 return __edgeNeighbor(o,i,idx);
912 template <
class NodeData,
class Real>
914 if(!parent){
return NULL;}
915 int pIndex=int(
this-(parent->children));
916 int aIndex,x[DIMENSION];
919 aIndex=(~((i[0] ^ x[idx[0]]) | ((i[1] ^ x[idx[1]])<<1))) & 3;
920 pIndex^=(7 ^ (1<<o));
922 const OctNode* temp=parent->__faceNeighbor(idx[0],i[0]);
923 if(!temp || !temp->children){
return NULL;}
924 else{
return &temp->children[pIndex];}
927 const OctNode* temp=parent->__faceNeighbor(idx[1],i[1]);
928 if(!temp || !temp->children){
return NULL;}
929 else{
return &temp->children[pIndex];}
932 return &parent->children[pIndex];
935 const OctNode* temp=parent->__edgeNeighbor(o,i,idx);
936 if(!temp || !temp->children){
return temp;}
937 else{
return &temp->children[pIndex];}
941 template <
class NodeData,
class Real>
942 OctNode<NodeData,Real>* OctNode<NodeData,Real>::__edgeNeighbor(
int o,
const int i[2],
const int idx[2],
int forceChildren){
943 if(!parent){
return NULL;}
944 int pIndex=int(
this-(parent->children));
945 int aIndex,x[DIMENSION];
948 aIndex=(~((i[0] ^ x[idx[0]]) | ((i[1] ^ x[idx[1]])<<1))) & 3;
949 pIndex^=(7 ^ (1<<o));
951 OctNode* temp=parent->__faceNeighbor(idx[0],i[0],0);
952 if(!temp || !temp->children){
return NULL;}
953 else{
return &temp->children[pIndex];}
956 OctNode* temp=parent->__faceNeighbor(idx[1],i[1],0);
957 if(!temp || !temp->children){
return NULL;}
958 else{
return &temp->children[pIndex];}
961 return &parent->children[pIndex];
964 OctNode* temp=parent->__edgeNeighbor(o,i,idx,forceChildren);
965 if(!temp){
return NULL;}
967 if(forceChildren){temp->initChildren();}
970 return &temp->children[pIndex];
975 template <
class NodeData,
class Real>
978 if(!parent){
return NULL;}
980 pIndex=int(
this-(parent->children));
981 aIndex=(cornerIndex ^ pIndex);
984 return &parent->children[pIndex];
988 if(!temp || !temp->
children){
return temp;}
989 else{
return &temp->
children[pIndex];}
992 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(0,cornerIndex & 1);
993 if(!temp || !temp->
children){
return NULL;}
994 else{
return & temp->
children[pIndex];}
997 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(1,(cornerIndex & 2)>>1);
998 if(!temp || !temp->
children){
return NULL;}
999 else{
return & temp->
children[pIndex];}
1002 const OctNode* temp=((
const OctNode*)parent)->__faceNeighbor(2,(cornerIndex & 4)>>2);
1003 if(!temp || !temp->
children){
return NULL;}
1004 else{
return & temp->
children[pIndex];}
1008 if(!temp || !temp->
children){
return NULL;}
1009 else{
return & temp->
children[pIndex];}
1013 if(!temp || !temp->
children){
return NULL;}
1014 else{
return & temp->
children[pIndex];}
1018 if(!temp || !temp->
children){
return NULL;}
1019 else{
return & temp->
children[pIndex];}
1023 template <
class NodeData,
class Real>
1025 int pIndex,aIndex=0;
1026 if(!parent){
return NULL;}
1028 pIndex=int(
this-(parent->children));
1029 aIndex=(cornerIndex ^ pIndex);
1032 return &parent->children[pIndex];
1036 if(!temp){
return NULL;}
1044 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(0,cornerIndex & 1,0);
1045 if(!temp || !temp->
children){
return NULL;}
1046 else{
return & temp->
children[pIndex];}
1049 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(1,(cornerIndex & 2)>>1,0);
1050 if(!temp || !temp->
children){
return NULL;}
1051 else{
return & temp->
children[pIndex];}
1054 OctNode* temp=((
OctNode*)parent)->__faceNeighbor(2,(cornerIndex & 4)>>2,0);
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];}
1078 template<
class NodeData,
class Real>
1080 template<
class NodeData,
class Real>
1082 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;}}}
1084 template<
class NodeData,
class Real>
1086 template<
class NodeData,
class Real>
1089 if( neighbors )
delete[] neighbors;
1093 template<
class NodeData,
class Real>
1096 if( neighbors )
delete[] neighbors;
1101 template<
class NodeData ,
class Real >
1104 if( !neighbors[
d].neighbors[1][1][1] || !neighbors[
d].neighbors[1][1][1]->
isInside( p ) )
1108 if( !
d ) neighbors[
d].neighbors[1][1][1] =
root;
1113 int i , j , k , x1 , y1 , z1 , x2 , y2 , z2;
1122 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1167 i=x1<<1 , j=y1<<1 , k=z1<<1;
1175 return neighbors[
d];
1177 template<
class NodeData ,
class Real >
1180 if( !neighbors[
d].neighbors[1][1][1] || !neighbors[
d].neighbors[1][1][1]->
isInside( p ) )
1184 if( !
d ) neighbors[
d].neighbors[1][1][1] =
root;
1187 Neighbors3& temp = getNeighbors(
root , p ,
d-1 );
1189 int i , j , k , x1 , y1 , z1 , x2 , y2 , z2;
1192 temp.neighbors[1][1][1]->centerAndWidth( c , w );
1197 if( !temp.neighbors[1][1][1] || !temp.neighbors[1][1][1]->children )
1201 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1202 neighbors[
d].neighbors[x2+i][y2+j][z2+k] = &temp.neighbors[1][1][1]->children[
Cube::CornerIndex(i,j,k)];
1207 if( temp.neighbors[i][1][1] )
1208 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)];
1210 if( temp.neighbors[1][j][1] )
1211 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)];
1213 if( temp.neighbors[1][1][k] )
1214 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)];
1218 if( temp.neighbors[i][j][1] )
1219 for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][j][z2+k] = &temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];
1221 if( temp.neighbors[i][1][k] )
1222 for( j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[i][y2+j][k] = &temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];
1224 if( temp.neighbors[1][j][k] )
1225 for( i=0 ; i<2 ; i++ ) neighbors[
d].neighbors[x2+i][j][k] = &temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];
1228 i=x1<<1 , j=y1<<1 , k=z1<<1;
1229 if( temp.neighbors[i][j][k] )
1230 neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1233 return neighbors[
d];
1236 template<
class NodeData ,
class Real >
1239 int d = node->depth();
1240 if( node==neighbors[
d].neighbors[1][1][1] )
1243 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
if( !neighbors[
d].neighbors[i][j][k] ) reset =
true;
1244 if( reset ) neighbors[
d].neighbors[1][1][1] = NULL;
1246 if( node!=neighbors[
d].neighbors[1][1][1] )
1248 neighbors[
d].clear();
1250 if( !node->parent ) neighbors[
d].neighbors[1][1][1] = node;
1253 int i,j,k,x1,y1,z1,x2,y2,z2;
1254 int idx=int(node-node->parent->children);
1260 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1264 Neighbors3& temp=setNeighbors(node->parent);
1268 if(temp.neighbors[i][1][1]){
1269 if(!temp.neighbors[i][1][1]->children){temp.neighbors[i][1][1]->initChildren();}
1270 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)];}}
1273 if(temp.neighbors[1][j][1]){
1274 if(!temp.neighbors[1][j][1]->children){temp.neighbors[1][j][1]->initChildren();}
1275 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)];}}
1278 if(temp.neighbors[1][1][k]){
1279 if(!temp.neighbors[1][1][k]->children){temp.neighbors[1][1][k]->initChildren();}
1280 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)];}}
1285 if(temp.neighbors[i][j][1]){
1286 if(!temp.neighbors[i][j][1]->children){temp.neighbors[i][j][1]->initChildren();}
1287 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1290 if(temp.neighbors[i][1][k]){
1291 if(!temp.neighbors[i][1][k]->children){temp.neighbors[i][1][k]->initChildren();}
1292 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1295 if(temp.neighbors[1][j][k]){
1296 if(!temp.neighbors[1][j][k]->children){temp.neighbors[1][j][k]->initChildren();}
1297 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1301 i=x1<<1; j=y1<<1; k=z1<<1;
1302 if(temp.neighbors[i][j][k]){
1303 if(!temp.neighbors[i][j][k]->children){temp.neighbors[i][j][k]->initChildren();}
1304 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1308 return neighbors[
d];
1312 template<
class NodeData ,
class Real >
1315 int d = node->depth();
1316 if( node==neighbors[
d].neighbors[1][1][1] )
1319 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
if( flags[i][j][k] && !neighbors[
d].neighbors[i][j][k] ) reset =
true;
1320 if( reset ) neighbors[
d].neighbors[1][1][1] = NULL;
1322 if( node!=neighbors[
d].neighbors[1][1][1] )
1324 neighbors[
d].clear();
1326 if( !node->parent ) neighbors[
d].neighbors[1][1][1] = node;
1329 int x1,y1,z1,x2,y2,z2;
1330 int idx=int(node-node->parent->children);
1333 for(
int i=0 ; i<2 ; i++ )
1334 for(
int j=0 ; j<2 ; j++ )
1335 for(
int k=0 ; k<2 ; k++ )
1336 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1338 Neighbors3& temp=setNeighbors( node->parent , flags );
1343 if( temp.neighbors[i][1][1] )
1345 if( flags[i][1][1] && !temp.neighbors[i][1][1]->children ) temp.neighbors[i][1][1]->initChildren();
1346 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)];
1351 if( temp.neighbors[1][j][1] )
1353 if( flags[1][j][1] && !temp.neighbors[1][j][1]->children ) temp.neighbors[1][j][1]->initChildren();
1354 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)];
1359 if( temp.neighbors[1][1][k] )
1361 if( flags[1][1][k] && !temp.neighbors[1][1][k]->children ) temp.neighbors[1][1][k]->initChildren();
1362 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)];
1368 int i=x1<<1 , j=y1<<1;
1369 if( temp.neighbors[i][j][1] )
1371 if( flags[i][j][1] && !temp.neighbors[i][j][1]->children ) temp.neighbors[i][j][1]->initChildren();
1372 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)];
1376 int i=x1<<1 , k=z1<<1;
1377 if( temp.neighbors[i][1][k] )
1379 if( flags[i][1][k] && !temp.neighbors[i][1][k]->children ) temp.neighbors[i][1][k]->initChildren();
1380 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)];
1384 int j=y1<<1 , k=z1<<1;
1385 if( temp.neighbors[1][j][k] )
1387 if( flags[1][j][k] && !temp.neighbors[1][j][k]->children ) temp.neighbors[1][j][k]->initChildren();
1388 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)];
1394 int i=x1<<1 , j=y1<<1 , k=z1<<1;
1395 if( temp.neighbors[i][j][k] )
1397 if( flags[i][j][k] && !temp.neighbors[i][j][k]->children ) temp.neighbors[i][j][k]->initChildren();
1398 if( temp.neighbors[i][j][k]->children ) neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1403 return neighbors[
d];
1406 template<
class NodeData,
class Real>
1408 int d=node->depth();
1409 if(node!=neighbors[
d].neighbors[1][1][1]){
1410 neighbors[
d].clear();
1412 if(!node->parent){neighbors[
d].neighbors[1][1][1]=node;}
1414 int i,j,k,x1,y1,z1,x2,y2,z2;
1415 int idx=int(node-node->parent->children);
1421 neighbors[
d].neighbors[x2+i][y2+j][z2+k]=&node->parent->children[
Cube::CornerIndex(i,j,k)];
1425 Neighbors3& temp=getNeighbors(node->parent);
1429 if(temp.neighbors[i][1][1] && temp.neighbors[i][1][1]->children){
1430 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)];}}
1433 if(temp.neighbors[1][j][1] && temp.neighbors[1][j][1]->children){
1434 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)];}}
1437 if(temp.neighbors[1][1][k] && temp.neighbors[1][1][k]->children){
1438 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)];}}
1443 if(temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children){
1444 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1447 if(temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children){
1448 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1451 if(temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children){
1452 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1456 i=x1<<1; j=y1<<1; k=z1<<1;
1457 if(temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children){
1458 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1462 return neighbors[node->depth()];
1468 template<
class NodeData,
class Real>
1470 template<
class NodeData,
class Real>
1472 for(
int i=0;i<3;i++){
for(
int j=0;j<3;j++){
for(
int k=0;k<3;k++){neighbors[i][j][k]=NULL;}}}
1474 template<
class NodeData,
class Real>
1476 template<
class NodeData,
class Real>
1478 if(neighbors){
delete[] neighbors;}
1482 template<
class NodeData,
class Real>
1484 if(neighbors){
delete[] neighbors;}
1489 template<
class NodeData,
class Real>
1492 if(node!=neighbors[
d].neighbors[1][1][1]){
1493 neighbors[
d].clear();
1495 if(!node->
parent){neighbors[
d].neighbors[1][1][1]=node;}
1497 int i,j,k,x1,y1,z1,x2,y2,z2;
1498 int idx=int(node-node->
parent->children);
1508 ConstNeighbors3& temp=getNeighbors(node->
parent);
1512 if(temp.neighbors[i][1][1] && temp.neighbors[i][1][1]->children){
1513 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)];}}
1516 if(temp.neighbors[1][j][1] && temp.neighbors[1][j][1]->children){
1517 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)];}}
1520 if(temp.neighbors[1][1][k] && temp.neighbors[1][1][k]->children){
1521 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)];}}
1526 if(temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children){
1527 for(k=0;k<2;k++){neighbors[
d].neighbors[i][j][z2+k]=&temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];}
1530 if(temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children){
1531 for(j=0;j<2;j++){neighbors[
d].neighbors[i][y2+j][k]=&temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];}
1534 if(temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children){
1535 for(i=0;i<2;i++){neighbors[
d].neighbors[x2+i][j][k]=&temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];}
1539 i=x1<<1; j=y1<<1; k=z1<<1;
1540 if(temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children){
1541 neighbors[
d].neighbors[i][j][k]=&temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1545 return neighbors[node->
depth()];
1547 template<
class NodeData,
class Real>
1550 int d=node->depth();
1555 if( node!=neighbors[
d].neighbors[1][1][1] )
1557 neighbors[
d].clear();
1559 if(
d==minDepth ) neighbors[
d].neighbors[1][1][1]=node;
1562 int i,j,k,x1,y1,z1,x2,y2,z2;
1563 int idx = int(node-node->parent->children);
1567 ConstNeighbors3& temp=getNeighbors( node->parent , minDepth );
1570 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1571 neighbors[
d].neighbors[x2+i][y2+j][z2+k] = &node->parent->children[
Cube::CornerIndex(i,j,k) ];
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)];
1588 if( temp.neighbors[i][j][1] && temp.neighbors[i][j][1]->children )
1589 for( k=0 ; k<2 ; k++ ) neighbors[
d].neighbors[i][j][z2+k] = &temp.neighbors[i][j][1]->children[
Cube::CornerIndex(x2,y2,k)];
1592 if( temp.neighbors[i][1][k] && temp.neighbors[i][1][k]->children )
1593 for( j=0 ; j<2 ; j++ ) neighbors[
d].neighbors[i][y2+j][k] = &temp.neighbors[i][1][k]->children[
Cube::CornerIndex(x2,j,z2)];
1596 if( temp.neighbors[1][j][k] && temp.neighbors[1][j][k]->children )
1597 for( i=0 ; i<2 ; i++ ) neighbors[
d].neighbors[x2+i][j][k] = &temp.neighbors[1][j][k]->children[
Cube::CornerIndex(i,y2,z2)];
1600 i=x1<<1 , j=y1<<1 , k=z1<<1;
1601 if( temp.neighbors[i][j][k] && temp.neighbors[i][j][k]->children )
1602 neighbors[
d].neighbors[i][j][k] = &temp.neighbors[i][j][k]->children[
Cube::CornerIndex(x2,y2,z2)];
1605 return neighbors[node->depth()];
1610 template<
class NodeData ,
class Real >
1613 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;
1615 template<
class NodeData ,
class Real >
1618 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;
1620 template<
class NodeData ,
class Real >
1626 template<
class NodeData ,
class Real >
1632 template<
class NodeData ,
class Real >
1635 if( neighbors )
delete[] neighbors;
1638 template<
class NodeData ,
class Real >
1641 if( neighbors )
delete[] neighbors;
1644 template<
class NodeData ,
class Real >
1647 if( neighbors )
delete[] neighbors;
1653 template<
class NodeData ,
class Real >
1656 if( neighbors )
delete[] neighbors;
1662 template<
class NodeData ,
class Real >
1666 if( node!=neighbors[
d].neighbors[2][2][2] )
1668 neighbors[
d].clear();
1670 if( !node->
parent ) neighbors[
d].neighbors[2][2][2]=node;
1673 getNeighbors( node->
parent );
1675 int x1 , y1 , z1 , x2 , y2 , z2;
1682 int fx0 = x2+1 , fy0 = y2+1 , fz0 = z2+1;
1683 int cx1 = x1*2+1 , cy1 = y1*2+1 , cz1 = z1*2+1;
1684 int cx2 = x2*2+1 , cy2 = y2*2+1 , cz2 = z2*2+1;
1685 int fx1 = x1*3 , fy1 = y1*3 , fz1 = z1*3;
1686 int fx2 = x2*4 , fy2 = y2*4 , fz2 = z2*4;
1689 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1694 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1697 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1700 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1703 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1706 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1709 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1714 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1717 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1720 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1723 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1726 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1729 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1732 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1735 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1738 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1741 for( k=0 ; k<2 ; k++ )
1744 for( j=0 ; j<2 ; j++ )
1747 for( i=0 ; i<2 ; i++ )
1752 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1755 for( i=0 ; i<2 ; i++ )
for( j=0 ; j<2 ; j++ )
1758 for( i=0 ; i<2 ; i++ )
for( k=0 ; k<2 ; k++ )
1761 for( j=0 ; j<2 ; j++ )
for( k=0 ; k<2 ; k++ )
1764 for( i=0 ; i<2 ; i++ )
1767 for( j=0 ; j<2 ; j++ )
1770 for( k=0 ; k<2 ; k++ )
1776 return neighbors[
d];
1778 template<
class NodeData ,
class Real >
1782 if( node!=neighbors[
d].neighbors[2][2][2] )
1784 neighbors[
d].clear();
1786 if( !node->
parent ) neighbors[
d].neighbors[2][2][2]=node;
1789 setNeighbors( node->
parent , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1791 int x1 , y1 , z1 , x2 , y2 , z2 , ii , jj , kk;
1795 for(
int i=xStart ; i<xEnd ; i++ )
1800 for(
int j=yStart ; j<yEnd ; j++ )
1805 for(
int k=zStart ; k<zEnd ; k++ )
1820 return neighbors[
d];
1822 template<
class NodeData ,
class Real >
1826 if( node!=neighbors[
d].neighbors[2][2][2] )
1828 neighbors[
d].clear();
1830 if(!node->
parent) neighbors[
d].neighbors[2][2][2]=node;
1833 getNeighbors( node->
parent );
1835 int x1,y1,z1,x2,y2,z2,ii,jj,kk;
1839 for(
int i=0;i<5;i++)
1844 for(
int j=0;j<5;j++)
1849 for(
int k=0;k<5;k++)
1861 return neighbors[
d];
1865 template <
class NodeData,
class Real>
1867 FILE* fp=fopen(fileName,
"wb");
1873 template <
class NodeData,
class Real>
1879 template <
class NodeData,
class Real>
1881 FILE* fp=fopen(fileName,
"rb");
1887 template <
class NodeData,
class Real>
1901 template<
class NodeData,
class Real>
1906 template<
class NodeData,
class Real>