29#include <unordered_map>
31#include "poisson_exceptions.h"
32#include "octree_poisson.h"
35#if defined WIN32 || defined _WIN32
36 #if !defined __MINGW32__
42#define ITERATION_POWER 1.0/3
43#define MEMORY_ALLOCATOR_BLOCK_SIZE 1<<12
57#if defined _WIN32 && !defined __MINGW32__
113 if( threads<=0 ) threads = 1;
115 std::vector< std::pair< int , int > >
spans( this->maxDepth , std::pair< int , int >( -1 , -1 ) );
118 cData.offsets.resize( this->maxDepth , -1 );
121 int start=
rootNode->nodeData.nodeIndex , end=start;
124 spans[d] = std::pair< int , int >( start , end+1 );
129 while( start< end && !
treeNodes[start]->children ) start++;
130 while( end> start && !
treeNodes[end ]->children ) end--;
131 if( start==end && !
treeNodes[start]->children )
break;
138 std::vector< int > count( threads );
139#pragma omp parallel for num_threads( threads )
140 for(
int t=0 ;
t<threads ;
t++ )
142 TreeOctNode::ConstNeighborKey3 neighborKey;
148 int start =
spans[d].first , end =
spans[d].second , width = end-start;
149 for(
int i=start + (width*
t)/threads ; i<start + (width*(
t+1))/threads ; i++ )
153 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node ,
minDepth );
164 xx += x ,
yy += y ,
zz += z;
165 if( neighbors.neighbors[
xx][
yy][
zz] )
169 neighbors.neighbors[
xx][
yy][
zz]->depthAndOffset(
_d , _off );
171 if( _off[0]==off[0] && _off[1]==off[1] && _off[2]==off[2] )
176 else fprintf(
stderr ,
"[WARNING] How did we leave the subtree?\n" );
185 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.neighbors[d];
191 xx += x ,
yy += y ,
zz += z;
192 if( neighborKey.neighbors[d].neighbors[
xx][
yy][
zz] )
208 std::vector< int > offsets( threads+1 );
210 for (
int t = 0;
t < threads;
t++)
213 offsets[
t + 1] = offsets[
t] + count[
t];
217#pragma omp parallel for num_threads( threads )
218 for (
int t = 0;
t < threads;
t++)
222 int start =
spans[d].first, end =
spans[d].second, width = end - start;
223 for (
int i = start + (width*
t) / threads; i < start + (width*(
t + 1)) / threads; i++)
252 if( threads<=0 ) threads = 1;
254 std::vector< std::vector< int > >
cornerCount( threads );
255 for(
int t=0 ;
t<threads ;
t++ )
cornerCount[
t].resize( res*res*res , 0 );
257#pragma omp parallel for num_threads( threads )
258 for(
int t=0 ;
t<threads ;
t++ )
261 TreeOctNode::ConstNeighborKey3 neighborKey;
264 for(
int i=(range*
t)/threads ; i<(range*(
t+1))/threads ; i++ )
271 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node , depth );
282 xx += x ,
yy += y ,
zz += z;
283 if( neighbors.neighbors[
xx][
yy][
zz] )
290 if(
cornerOwner )
_cornerCount[ ( ( off[0]>>(d-depth) ) * res * res) + ( ( off[1]>>(d-depth) ) * res) + ( off[2]>>(d-depth) ) ]++;
295 for(
int i=0 ; i<res*res*res ; i++ )
309 if( threads<=0 ) threads = 1;
310 std::vector< std::pair< int , int > >
spans( this->maxDepth , std::pair< int , int >( -1 , -1 ) );
314 eData.offsets.resize( this->maxDepth , -1 );
317 int start=
rootNode->nodeData.nodeIndex , end=start;
320 spans[d] = std::pair< int , int >( start , end+1 );
325 while( start< end && !
treeNodes[start]->children ) start++;
326 while( end> start && !
treeNodes[end ]->children ) end--;
327 if( start==end && !
treeNodes[start]->children )
break;
334 std::vector< int > count( threads );
335#pragma omp parallel for num_threads( threads )
336 for(
int t=0 ;
t<threads ;
t++ )
338 TreeOctNode::ConstNeighborKey3 neighborKey;
344 int start =
spans[d].first , end =
spans[d].second , width = end-start;
345 for(
int i=start + (width*
t)/threads ; i<start + (width*(
t+1))/threads ; i++ )
348 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node ,
minDepth );
358 int ii ,
jj , x , y , z;
363 case 0: y =
ii , z =
jj , x = 1 ;
break;
364 case 1: x =
ii , z =
jj , y = 1 ;
break;
365 case 2: x =
ii , y =
jj , z = 1 ;
break;
380 case 0: y =
ii , z =
jj , x = 1 ;
break;
381 case 1: x =
ii , z =
jj , y = 1 ;
break;
382 case 2: x =
ii , y =
jj , z = 1 ;
break;
384 if( neighbors.neighbors[x][y][z] )
394 std::vector< int > offsets( threads+1 );
396 for (
int t = 0;
t < threads;
t++)
399 offsets[
t + 1] = offsets[
t] + count[
t];
403#pragma omp parallel for num_threads( threads )
404 for (
int t = 0;
t < threads;
t++)
408 int start =
spans[d].first, end =
spans[d].second, width = end - start;
409 for (
int i = start + (width*
t) / threads; i < start + (width*(
t + 1)) / threads; i++)
439 if( threads<=0 ) threads = 1;
441 std::vector< std::vector< int > >
edgeCount( threads );
442 for(
int t=0 ;
t<threads ;
t++ )
edgeCount[
t].resize( res*res*res , 0 );
444#pragma omp parallel for num_threads( threads )
445 for(
int t=0 ;
t<threads ;
t++ )
448 TreeOctNode::ConstNeighborKey3 neighborKey;
451 for(
int i=(range*
t)/threads ; i<(range*(
t+1))/threads ; i++ )
454 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node , depth );
466 int ii ,
jj , x , y , z;
471 case 0: y =
ii , z =
jj , x = 1 ;
break;
472 case 1: x =
ii , z =
jj , y = 1 ;
break;
473 case 2: x =
ii , y =
jj , z = 1 ;
break;
477 if(
edgeOwner )
_edgeCount[ ( ( off[0]>>(d-depth) ) * res * res) + ( ( off[1]>>(d-depth) ) * res) + ( off[2]>>(d-depth) ) ]++;
482 for(
int i=0 ; i<res*res*res ; i++ )
502 centerWeightContribution=0;
506 constraint = solution = 0;
523 if(
mem>maxMemoryUsage){maxMemoryUsage=
mem;}
533 postNormalSmooth = 0;
534 _constrainValues =
false;
540 double x ,
dxdy ,
dxdydz , dx[DIMENSION][SPLAT_ORDER+1];
542 TreeOctNode::Neighbors3& neighbors = neighborKey.setNeighbors( node );
548 for(
int i=0 ; i<3 ; i++ )
552 x = (
center[i] - position[i] - width ) / width;
553 dx[i][0] = 1.125+1.500*x+0.500*x*x;
554 x = (
center[i] - position[i] ) / width;
555 dx[i][1] = 0.750 - x*x;
557 dx[i][2] = 1. - dx[i][1] - dx[i][0];
559 x = ( position[i] -
center[i] ) / width;
570 dx[i][1] = 1. - dx[i][0];
575# error Splat order not supported
578 for(
int i=off[0] ; i<=off[0]+SPLAT_ORDER ; i++ )
for(
int j=off[1] ; j<=off[1]+SPLAT_ORDER ; j++ )
580 dxdy = dx[0][i] * dx[1][j];
581 for(
int k=off[2] ; k<=off[2]+SPLAT_ORDER ; k++ )
582 if( neighbors.neighbors[i][j][k] )
584 dxdydz = dxdy * dx[2][k];
586 int idx =_node->
nodeData.normalIndex;
590 n[0] = n[1] = n[2] = 0;
591 _node->nodeData.nodeIndex = 0;
592 idx = _node->nodeData.normalIndex = int(normals->size());
593 normals->push_back(n);
595 (*normals)[idx] += normal *
Real( dxdydz );
601 Real Octree<Degree>::NonLinearSplatOrientedPoint(
const Point3D<Real>& position ,
const Point3D<Real>& normal ,
int splatDepth ,
Real samplesPerNode ,
609 Point3D< Real > myCenter;
611 myCenter[0] = myCenter[1] = myCenter[2] =
Real(0.5);
615 while( temp->depth()<splatDepth )
617 if( !temp->children )
619 fprintf( stderr ,
"Octree<Degree>::NonLinearSplatOrientedPoint error\n" );
623 temp=&temp->children[cIndex];
625 if(cIndex&1) myCenter[0] += myWidth/2;
626 else myCenter[0] -= myWidth/2;
627 if(cIndex&2) myCenter[1] += myWidth/2;
628 else myCenter[1] -= myWidth/2;
629 if(cIndex&4) myCenter[2] += myWidth/2;
630 else myCenter[2] -= myWidth/2;
633 NonLinearGetSampleDepthAndWeight( temp , position , samplesPerNode , newDepth , alpha );
635 if( newDepth<minDepth ) newDepth=
Real(minDepth);
637 int topDepth=int(ceil(newDepth));
639 dx = 1.0-(topDepth-newDepth);
640 if( topDepth<=minDepth )
650 while( temp->depth()>topDepth ) temp=temp->parent;
651 while( temp->depth()<topDepth )
653 if(!temp->children) temp->initChildren();
655 temp=&temp->children[cIndex];
657 if(cIndex&1) myCenter[0] += myWidth/2;
658 else myCenter[0] -= myWidth/2;
659 if(cIndex&2) myCenter[1] += myWidth/2;
660 else myCenter[1] -= myWidth/2;
661 if(cIndex&4) myCenter[2] += myWidth/2;
662 else myCenter[2] -= myWidth/2;
664 width = 1.0 / ( 1<<temp->depth() );
665 n = normal * alpha /
Real( pow( width , 3 ) ) *
Real( dx );
666 NonLinearSplatOrientedPoint( temp , position , n );
671 width = 1.0 / ( 1<<temp->depth() );
673 n = normal * alpha /
Real( pow( width , 3 ) ) *
Real( dx );
674 NonLinearSplatOrientedPoint( temp , position , n );
679 void Octree<Degree>::NonLinearGetSampleDepthAndWeight(
TreeOctNode* node,
const Point3D<Real>& position,
Real samplesPerNode,
Real& depth,
Real& weight){
681 weight =
Real(1.0)/NonLinearGetSampleWeight(temp,position);
682 if( weight>=samplesPerNode ) depth=
Real( temp->depth() + log( weight / samplesPerNode ) / log(
double(1<<(DIMENSION-1))) );
685 Real oldAlpha,newAlpha;
686 oldAlpha=newAlpha=weight;
687 while( newAlpha<samplesPerNode && temp->parent )
691 newAlpha=
Real(1.0)/NonLinearGetSampleWeight(temp,position);
693 depth =
Real( temp->depth() + log( newAlpha / samplesPerNode ) / log( newAlpha / oldAlpha ) );
695 weight=
Real(pow(
double(1<<(DIMENSION-1)),-
double(depth)));
699 Real Octree<Degree>::NonLinearGetSampleWeight(
TreeOctNode* node ,
const Point3D<Real>& position )
702 double x,dxdy,dx[DIMENSION][3];
703 TreeOctNode::Neighbors3& neighbors=neighborKey.setNeighbors( node );
705 Point3D<Real> center;
707 node->centerAndWidth(center,w);
710 for(
int i=0 ; i<DIMENSION ; i++ )
712 x = ( center[i] - position[i] - width ) / width;
713 dx[i][0] = 1.125 + 1.500*x + 0.500*x*x;
714 x = ( center[i] - position[i] ) / width;
715 dx[i][1] = 0.750 - x*x;
717 dx[i][2] = 1.0 - dx[i][1] - dx[i][0];
720 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
722 dxdy = dx[0][i] * dx[1][j];
723 for(
int k=0 ; k<3 ; k++ )
if( neighbors.neighbors[i][j][k] )
724 weight +=
Real( dxdy * dx[2][k] * neighbors.neighbors[i][j][k]->nodeData.centerWeightContribution );
726 return Real( 1.0 / weight );
730 int Octree<Degree>::NonLinearUpdateWeightContribution(
TreeOctNode* node ,
const Point3D<Real>& position ,
Real weight )
732 TreeOctNode::Neighbors3& neighbors = neighborKey.setNeighbors( node );
733 double x,dxdy,dx[DIMENSION][3];
735 Point3D<Real> center;
737 node->centerAndWidth( center , w );
739 const double SAMPLE_SCALE = 1. / ( 0.125 * 0.125 + 0.75 * 0.75 + 0.125 * 0.125 );
741 for(
int i=0 ; i<DIMENSION ; i++ )
743 x = ( center[i] - position[i] - width ) / width;
744 dx[i][0] = 1.125 + 1.500*x + 0.500*x*x;
745 x = ( center[i] - position[i] ) / width;
746 dx[i][1] = 0.750 - x*x;
747 dx[i][2] = 1. - dx[i][1] - dx[i][0];
750 dx[i][0] *= SAMPLE_SCALE;
753 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
755 dxdy = dx[0][i] * dx[1][j] * weight;
756 for(
int k=0 ; k<3 ; k++ )
if( neighbors.neighbors[i][j][k] )
757 neighbors.neighbors[i][j][k]->nodeData.centerWeightContribution +=
Real( dxdy * dx[2][k] );
762 template<
int Degree >
template<
typename Po
intNT>
int
782 tree.setFullDepth( _minDepth );
784 while (
cnt != input_->size ())
787 p[0] = (*input_)[
cnt].x;
788 p[1] = (*input_)[
cnt].y;
789 p[2] = (*input_)[
cnt].z;
791 for (i = 0; i < DIMENSION; i++)
793 if( !
cnt || p[i]<min[i] ) min[i] = p[i];
794 if( !
cnt || p[i]>max[i] ) max[i] = p[i];
799 scale = std::max< Real >( max[0]-min[0] , std::max< Real >( max[1]-min[1] , max[2]-min[2] ) );
803 for( i=0 ; i<DIMENSION ; i++ )
center[i] -= scale/2;
807 while (
cnt != input_->size ())
809 position[0] = (*input_)[
cnt].x;
810 position[1] = (*input_)[
cnt].y;
811 position[2] = (*input_)[
cnt].z;
812 normal[0] = (*input_)[
cnt].normal_x;
813 normal[1] = (*input_)[
cnt].normal_y;
814 normal[2] = (*input_)[
cnt].normal_z;
817 for( i=0 ; i<DIMENSION ; i++ ) position[i] = ( position[i]-
center[i] ) / scale;
821 if( i!=DIMENSION )
continue;
828 NonLinearUpdateWeightContribution(
temp , position , weight );
829 if( !
temp->children )
temp->initChildren();
841 NonLinearUpdateWeightContribution(
temp , position , weight );
845 normals =
new std::vector< Point3D<Real> >();
847 while (
cnt != input_->size ())
849 position[0] = (*input_)[
cnt].x;
850 position[1] = (*input_)[
cnt].y;
851 position[2] = (*input_)[
cnt].z;
852 normal[0] = (*input_)[
cnt].normal_x;
853 normal[1] = (*input_)[
cnt].normal_y;
854 normal[2] = (*input_)[
cnt].normal_z;
856 for( i=0 ; i<DIMENSION ; i++ ) position[i] = ( position[i]-
center[i] ) / scale;
860 if( i!=DIMENSION )
continue;
862 if( l!=l || l<=
EPSILON )
continue;
891 alpha = NonLinearGetSampleWeight(
temp , position );
893 for( i=0 ; i<DIMENSION ; i++ ) normal[i]*=
alpha;
896 if(!
temp->children){
temp->initChildren();}
908 NonLinearSplatOrientedPoint(
temp , position , normal );
913 if( _constrainValues )
921 int idx =
temp->nodeData.pointIndex;
925 p[0] = p[1] = p[2] = 0;
926 idx =
int( _points.size() );
928 temp->nodeData.pointIndex = idx;
937 if( !
temp->children )
break;
952 if( _constrainValues )
954 if( n->nodeData.pointIndex!=-1 )
956 int idx = n->nodeData.pointIndex;
957 _points[idx].position /= _points[idx].weight;
959 else _points[idx].weight *= (1<<
maxDepth);
962#if FORCE_NEUMANN_FIELD
965 int d , off[3] , res;
966 node->depthAndOffset( d , off );
968 if( node->nodeData.normalIndex<0 )
continue;
970 for(
int d=0 ; d<3 ; d++ )
if( off[d]==0 || off[d]==res-1 ) normal[d] = 0;
983 radius = 0.5 + 0.5 *
Degree;
986 fData.set(
maxDepth ,
true , reflectBoundary );
993 TreeOctNode::NeighborKey5
nKey;
1000 int c =
int( node - node->parent->children );
1011 _sNodes.set( tree );
1014 template<
int Degree >
1018 const int min[] = { std::max<int>( 0 , d[0]+0 ) , std::max<int>( 0 , d[1]+0 ) , std::max<int>( 0 , d[2]+0 ) };
1019 const int max[] = { std::min<int>( 2 , d[0]+2 ) , std::min<int>( 2 , d[1]+2 ) , std::min<int>( 2 , d[2]+2 ) };
1020 for(
int i=min[0] ; i<=max[0] ; i++ )
for(
int j=min[1] ; j<=max[1] ; j++ )
1022 if( !hasPoints[i][j] )
continue;
1023 const PointInfo* pInfo = points[i][j];
1026 for(
int k=min[2] ; k<=max[2] ; k++ )
1027 if( pInfo[k].weightedValue )
1028 v += pInfo[k].splineValues[0][ii] * pInfo[k].splineValues[1][jj] * pInfo[k].splineValues[2][-d[2]+k];
1032 template<
int Degree>
1033 Real Octree<Degree>::GetLaplacian(
const int idx[DIMENSION] )
const
1035 return Real( fData.vvDotTable[idx[0]] * fData.vvDotTable[idx[1]] * fData.vvDotTable[idx[2]] * (fData.ddDotTable[idx[0]]+fData.ddDotTable[idx[1]]+fData.ddDotTable[idx[2]] ) );
1037 template<
int Degree >
1046 return GetLaplacian( symIndex );
1048 template<
int Degree >
1049 Real Octree< Degree >::GetDivergence(
const TreeOctNode* node1 ,
const TreeOctNode* node2 ,
const Point3D< Real >& normal1 )
const
1059 #if GRADIENT_DOMAIN_SOLUTION
1061 fData.Index( node2->off[0] , node1->off[0] ) ,
1062 fData.Index( node2->off[1] , node1->off[1] ) ,
1063 fData.Index( node2->off[2] , node1->off[2] )
1066 fData.Index( node1->off[0] , node2->off[0] ) ,
1067 fData.Index( node1->off[1] , node2->off[1] ) ,
1068 fData.Index( node1->off[2] , node2->off[2] )
1071 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1072#if GRADIENT_DOMAIN_SOLUTION
1073 return Real( dot * ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] ) );
1075 return -
Real( dot * ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] ) );
1078 template<
int Degree >
1079 Real Octree< Degree >::GetDivergenceMinusLaplacian(
const TreeOctNode* node1 ,
const TreeOctNode* node2 ,
Real value1 ,
const Point3D<Real>& normal1 )
const
1089 #if GRADIENT_DOMAIN_SOLUTION
1091 fData.Index( node2->off[0] , node1->off[0] ) ,
1092 fData.Index( node2->off[1] , node1->off[1] ) ,
1093 fData.Index( node2->off[2] , node1->off[2] )
1096 fData.Index( node1->off[0] , node2->off[0] ) ,
1097 fData.Index( node1->off[1] , node2->off[1] ) ,
1098 fData.Index( node1->off[2] , node2->off[2] )
1101 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1104 #
if GRADIENT_DOMAIN_SOLUTION
1105 - ( fData.ddDotTable[ symIndex[0]] + fData.ddDotTable[ symIndex[1]] + fData.ddDotTable[ symIndex[2]] ) * value1
1106 + ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] )
1108 - ( fData.ddDotTable[ symIndex[0]] + fData.ddDotTable[ symIndex[1]] + fData.ddDotTable[ symIndex[2]] ) * value1
1109 - ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] )
1114 template<
int Degree >
1115 void Octree< Degree >::SetMatrixRowBounds(
const TreeOctNode* node ,
int rDepth ,
const int rOff[3] ,
int& xStart ,
int& xEnd ,
int& yStart ,
int& yEnd ,
int& zStart ,
int& zEnd )
const
1117 xStart = 0 , yStart = 0 , zStart = 0 , xEnd = 5 , yEnd = 5 , zEnd = 5;
1119 node->depthAndOffset( depth , off );
1120 int width = 1 << ( depth-rDepth );
1121 off[0] -= rOff[0] << ( depth-rDepth ) , off[1] -= rOff[1] << ( depth-rDepth ) , off[2] -= rOff[2] << ( depth-rDepth );
1122 if( off[0]<0 ) xStart = -off[0];
1123 if( off[1]<0 ) yStart = -off[1];
1124 if( off[2]<0 ) zStart = -off[2];
1125 if( off[0]>=width ) xEnd = 4 - ( off[0]-width );
1126 if( off[1]>=width ) yEnd = 4 - ( off[1]-width );
1127 if( off[2]>=width ) zEnd = 4 - ( off[2]-width );
1129 template<
int Degree >
1130 int Octree< Degree >::GetMatrixRowSize(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 )
const {
return GetMatrixRowSize( neighbors5 , 0 , 5 , 0 , 5 , 0 , 5 ); }
1132 template<
int Degree >
1133 int Octree< Degree >::GetMatrixRowSize(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 ,
int xStart ,
int xEnd ,
int yStart ,
int yEnd ,
int zStart ,
int zEnd )
const
1136 for(
int x=xStart ; x<=2 ; x++ )
1137 for(
int y=yStart ; y<yEnd ; y++ )
1138 if( x==2 && y>2 )
continue;
1139 else for(
int z=zStart ; z<zEnd ; z++ )
1140 if( x==2 && y==2 && z>2 )
continue;
1141 else if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1146 template<
int Degree >
1147 int Octree< Degree >::SetMatrixRow(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 , MatrixEntry< float >* row ,
int offset ,
const double stencil[5][5][5] )
const
1149 return SetMatrixRow( neighbors5 , row , offset , stencil , 0 , 5 , 0 , 5 , 0 , 5 );
1152 template<
int Degree >
1153 int Octree< Degree >::SetMatrixRow(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 , MatrixEntry< float >* row ,
int offset ,
const double stencil[5][5][5] ,
int xStart ,
int xEnd ,
int yStart ,
int yEnd ,
int zStart ,
int zEnd )
const
1155 bool hasPoints[3][3];
1157 PointInfo samples[3][3][3];
1160 const TreeOctNode* node = neighbors5.neighbors[2][2][2];
1161 int index[] = { int( node->off[0] ) , int( node->off[1] ), int( node->off[2] ) };
1163 if( _constrainValues )
1166 node->depthAndOffset( d , idx );
1170 for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
1172 hasPoints[j][k] =
false;
1173 for(
int l=0 ; l<3 ; l++ )
1175 const TreeOctNode* _node = neighbors5.neighbors[j+1][k+1][l+1];
1176 if( _node && _node->nodeData.pointIndex!=-1 )
1178 const PointData& pData = _points[ _node->nodeData.pointIndex ];
1179 PointInfo& pointInfo = samples[j][k][l];
1180 Real weight = pData.weight;
1181 Point3D< Real > p = pData.position;
1182 for(
int s=0 ; s<3 ; s++ )
1184 pointInfo.splineValues[0][s] = float( fData.baseBSplines[ idx[0]+j-s][s]( p[0] ) );
1185 pointInfo.splineValues[1][s] = float( fData.baseBSplines[ idx[1]+k-s][s]( p[1] ) );
1186 pointInfo.splineValues[2][s] = float( fData.baseBSplines[ idx[2]+l-s][s]( p[2] ) );
1188 float value = pointInfo.splineValues[0][j] * pointInfo.splineValues[1][k] * pointInfo.splineValues[2][l];
1189 diagonal += value * value * weight;
1190 pointInfo.weightedValue = value * weight;
1191 for(
int s=0 ; s<3 ; s++ ) pointInfo.splineValues[0][s] *= pointInfo.weightedValue;
1192 hasPoints[j][k] =
true;
1194 else samples[j][k][l].weightedValue = 0;
1201 neighbors5.neighbors[2][2][2]->depthAndOffset( d , off );
1202 int mn = 2 , mx = (1<<d)-2;
1203 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
1204 for(
int x=xStart ; x<=2 ; x++ )
1205 for(
int y=yStart ; y<yEnd ; y++ )
1206 if( x==2 && y>2 )
continue;
1207 else for(
int z=zStart ; z<zEnd ; z++ )
1208 if( x==2 && y==2 && z>2 )
continue;
1209 else if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1211 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
1212 int _index[] = { int( _node->off[0] ) , int( _node->off[1] ), int( _node->off[2] ) };
1214 if( isInterior ) temp =
Real( stencil[x][y][z] );
1215 else temp = GetLaplacian( node , _node );
1216 if( _constrainValues )
1218 int _d[] = { _index[0]-index[0] , _index[1]-index[1] , _index[2]-index[2] };
1219 if( x==2 && y==2 && z==2 ) temp += diagonal;
1220 else temp += GetValue( samples , hasPoints , _d );
1222 if( x==2 && y==2 && z==2 ) temp /= 2;
1225 row[count].N = _node->nodeData.nodeIndex-offset;
1226 row[count].Value = temp;
1233 template<
int Degree >
1234 void Octree< Degree >::SetDivergenceStencil(
int depth , Point3D< double > *stencil ,
bool scatter )
const
1236 int offset[] = { 2 , 2 , 2 };
1239 int index1[3] , index2[3];
1240 if( scatter ) index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1241 else index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1242 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1244 int _offset[] = { x , y , z };
1246 if( scatter ) index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1247 else index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1256 #if GRADIENT_DOMAIN_SOLUTION
1258 fData.Index( index1[0] , index2[0] ) ,
1259 fData.Index( index1[1] , index2[1] ) ,
1260 fData.Index( index1[2] , index2[2] )
1263 fData.Index( index2[0] , index1[0] ) ,
1264 fData.Index( index2[1] , index1[1] ) ,
1265 fData.Index( index2[2] , index1[2] )
1269 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1270#if GRADIENT_DOMAIN_SOLUTION
1271 Point3D<double> temp;
1272 temp[0] = fData.dvDotTable[aSymIndex[0]] *
dot;
1273 temp[1] = fData.dvDotTable[aSymIndex[1]] *
dot;
1274 temp[2] = fData.dvDotTable[aSymIndex[2]] *
dot;
1275 stencil[25*x + 5*y + z] = temp;
1280 Point3D<double> temp;
1281 temp[0] = -fData.dvDotTable[aSymIndex[0]] *
dot;
1282 temp[1] = -fData.dvDotTable[aSymIndex[1]] *
dot;
1283 temp[2] = -fData.dvDotTable[aSymIndex[2]] *
dot;
1284 stencil[25*x + 5*y + z] = temp;
1292 template<
int Degree >
1293 void Octree< Degree >::SetLaplacianStencil(
int depth ,
double stencil[5][5][5] )
const
1295 int offset[] = { 2 , 2 , 2 };
1298 int index[] = { int( off[0] ) , int( off[1] ) , int( off[2] ) };
1299 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1301 int _offset[] = { x , y , z };
1304 int _index[] = { int( _off[0] ) , int( _off[1] ) , int( _off[2] ) };
1309 stencil[x][y][z] = GetLaplacian( symIndex );
1313 template<
int Degree >
1314 void Octree< Degree >::SetLaplacianStencils(
int depth , Stencil< double , 5 > stencils[2][2][2] )
const
1316 if( depth<=1 )
return;
1317 for(
int i=0 ; i<2 ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
1320 int offset[] = { 4+i , 4+j , 4+k };
1322 int index[] = { int( off[0] ) , int( off[1] ) , int( off[2] ) };
1323 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1325 int _offset[] = { x , y , z };
1328 int _index[] = { int( _off[0] ) , int( _off[1] ) , int( _off[2] ) };
1333 stencils[i][j][k].values[x][y][z] = GetLaplacian( symIndex );
1338 template<
int Degree >
1339 void Octree< Degree >::SetDivergenceStencils(
int depth , Stencil< Point3D< double > , 5 > stencils[2][2][2] ,
bool scatter )
const
1341 if( depth<=1 )
return;
1342 int index1[3] , index2[3];
1343 for(
int i=0 ; i<2 ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
1346 int offset[] = { 4+i , 4+j , 4+k };
1348 if( scatter ) index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1349 else index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1350 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1352 int _offset[] = { x , y , z };
1354 if( scatter ) index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1355 else index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1365 #if GRADIENT_DOMAIN_SOLUTION
1367 fData.Index( index1[0] , index2[0] ) ,
1368 fData.Index( index1[1] , index2[1] ) ,
1369 fData.Index( index1[2] , index2[2] )
1372 fData.Index( index2[0] , index1[0] ) ,
1373 fData.Index( index2[1] , index1[1] ) ,
1374 fData.Index( index2[2] , index1[2] )
1377 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1378#if GRADIENT_DOMAIN_SOLUTION
1379 stencils[i][j][k].values[x][y][z][0] = fData.dvDotTable[aSymIndex[0]] *
dot;
1380 stencils[i][j][k].values[x][y][z][1] = fData.dvDotTable[aSymIndex[1]] *
dot;
1381 stencils[i][j][k].values[x][y][z][2] = fData.dvDotTable[aSymIndex[2]] *
dot;
1383 stencils[i][j][k].values[x][y][z][0] = -fData.dvDotTable[aSymIndex[0]] *
dot;
1384 stencils[i][j][k].values[x][y][z][1] = -fData.dvDotTable[aSymIndex[1]] *
dot;
1385 stencils[i][j][k].values[x][y][z][2] = -fData.dvDotTable[aSymIndex[2]] *
dot;
1391 template<
int Degree >
1392 void Octree< Degree >::SetEvaluationStencils(
int depth , Stencil< double , 3 > stencil1[8] , Stencil< double , 3 > stencil2[8][8] )
const
1397 int off[] = { 2 , 2 , 2 };
1398 for(
int c=0 ; c<8 ; c++ )
1401 idx[0] *= fData.functionCount , idx[1] *= fData.functionCount , idx[2] *= fData.functionCount;
1402 for(
int x=0 ; x<3 ; x++ )
for(
int y=0 ; y<3 ; y++ )
for(
int z=0 ; z<3 ; z++ )
1405 int _offset[] = { x+1 , y+1 , z+1 };
1407 stencil1[c].values[x][y][z] = fData.valueTables[ idx[0]+int(_off[0]) ] * fData.valueTables[ idx[1]+int(_off[1]) ] * fData.valueTables[ idx[2]+int(_off[2]) ];
1412 for(
int _c=0 ; _c<8 ; _c++ )
1415 int _cx , _cy , _cz;
1417 int off[] = { 4+_cx , 4+_cy , 4+_cz };
1418 for(
int c=0 ; c<8 ; c++ )
1421 idx[0] *= fData.functionCount , idx[1] *= fData.functionCount , idx[2] *= fData.functionCount;
1422 for(
int x=0 ; x<3 ; x++ )
for(
int y=0 ; y<3 ; y++ )
for(
int z=0 ; z<3 ; z++ )
1425 int _offset[] = { x+1 , y+1 , z+1 };
1427 stencil2[_c][c].values[x][y][z] = fData.valueTables[ idx[0]+int(_off[0]) ] * fData.valueTables[ idx[1]+int(_off[1]) ] * fData.valueTables[ idx[2]+int(_off[2]) ];
1433 template<
int Degree >
1434 void Octree< Degree >::UpdateCoarserSupportBounds(
const TreeOctNode* node ,
int& startX ,
int& endX ,
int& startY ,
int& endY ,
int& startZ ,
int& endZ )
1438 int x , y , z , c = int( node - node->parent->children );
1440 if( x==0 ) endX = 4;
1442 if( y==0 ) endY = 4;
1444 if( z==0 ) endZ = 4;
1449 template<
int Degree >
1450 void Octree< Degree >::UpdateConstraintsFromCoarser(
const OctNode< TreeNodeData , Real >::NeighborKey5& neighborKey5 ,
TreeOctNode* node ,
Real* metSolution ,
const Stencil< double , 5 >& lapStencil )
const
1455 node->depthAndOffset( d , off );
1456 int mn = 4 , mx = (1<<d)-4;
1457 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
1460 int depth = node->depth();
1461 if( depth<=_minDepth )
return;
1462 int i = node->nodeData.nodeIndex;
1464 int startX = 0 , endX = 5 , startY = 0 , endY = 5 , startZ = 0 , endZ = 5;
1465 UpdateCoarserSupportBounds( node , startX , endX , startY , endY , startZ , endZ );
1467 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth-1];
1468 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
1469 if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1471 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
1472 Real _solution = metSolution[ _node->nodeData.nodeIndex ];
1474 if( isInterior ) node->nodeData.constraint -=
Real( lapStencil.values[x][y][z] * _solution );
1475 else node->nodeData.constraint -= GetLaplacian( _node , node ) * _solution;
1478 if( _constrainValues )
1481 node->depthAndOffset( d, idx );
1485 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth];
1486 for(
int x=1 ; x<4 ; x++ )
for(
int y=1 ; y<4 ; y++ )
for(
int z=1 ; z<4 ; z++ )
1487 if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.pointIndex!=-1 )
1489 const PointData& pData = _points[ neighbors5.neighbors[x][y][z]->nodeData.pointIndex ];
1490 Real pointValue = pData.value;
1491 Point3D< Real > p = pData.position;
1492 node->nodeData.constraint -=
1494 fData.baseBSplines[idx[0]][x-1]( p[0] ) *
1495 fData.baseBSplines[idx[1]][y-1]( p[1] ) *
1496 fData.baseBSplines[idx[2]][z-1]( p[2] ) *
1510 template<
int Degree >
1513 int start =
sNodes.nodeCount[depth] , end =
sNodes.nodeCount[depth+1] , range = end-start;
1515 if( !depth )
return;
1516 else if( depth==1 )
for(
int i=start ; i<end ; i++ )
Solution[i-start] +=
sNodes.treeNodes[0]->nodeData.solution;
1520#pragma omp parallel for num_threads( threads )
1521 for(
int t=0 ;
t<threads ;
t++ )
1523 TreeOctNode::NeighborKey3 neighborKey;
1524 neighborKey.set( depth );
1525 for(
int i=start+(range*
t)/threads ; i<start+(range*(
t+1))/threads ; i++ )
1529 sNodes.treeNodes[i]->depthAndOffset( d , off );
1530 for(
int d=0 ; d<3 ; d++ )
1532 if ( off[d] ==0 )
usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1533 else if( off[d]+1==(1<<depth) )
usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1534 else if( off[d]%2 )
usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1535 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1537 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1538 for(
int ii=0 ; ii<2 ; ii++ )
1540 int _ii = ii + usData[0].start;
1541 double dx = usData[0].v[ii];
1542 for(
int jj=0 ; jj<2 ; jj++ )
1544 int _jj = jj + usData[1].start;
1545 double dxy = dx * usData[1].v[jj];
1546 for(
int kk=0 ; kk<2 ; kk++ )
1548 int _kk = kk + usData[2].start;
1549 double dxyz = dxy * usData[2].v[kk];
1550 Solution[i-start] +=
Real( neighborKey.neighbors[depth-1].neighbors[_ii][_jj][_kk]->nodeData.solution * dxyz );
1558 start = sNodes.nodeCount[depth-1] , end = sNodes.nodeCount[depth] , range = end-start;
1559#pragma omp parallel for num_threads( threads )
1560 for(
int i=start ; i<end ; i++ ) sNodes.treeNodes[i]->nodeData.solution =
Real( 0. );
1562 template<
int Degree >
1563 void Octree< Degree >::DownSampleFinerConstraints(
int depth ,
SortedTreeNodes& sNodes )
const
1565 if( !depth )
return;
1566#pragma omp parallel for num_threads( threads )
1567 for(
int i=sNodes.nodeCount[depth-1] ; i<sNodes.nodeCount[depth] ; i++ )
1568 sNodes.treeNodes[i]->nodeData.constraint =
Real( 0 );
1572 sNodes.treeNodes[0]->nodeData.constraint =
Real( 0 );
1573 for(
int i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ ) sNodes.treeNodes[0]->nodeData.constraint += sNodes.treeNodes[i]->nodeData.constraint;
1576 std::vector< Vector< double > > constraints( threads );
1577 for(
int t=0 ; t<threads ; t++ ) constraints[t].Resize( sNodes.nodeCount[depth] - sNodes.nodeCount[depth-1] ) , constraints[t].SetZero();
1578 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1579 int lStart = sNodes.nodeCount[depth-1] , lEnd = sNodes.nodeCount[depth];
1581#pragma omp parallel for num_threads( threads )
1582 for(
int t=0 ; t<threads ; t++ )
1584 TreeOctNode::NeighborKey3 neighborKey;
1585 neighborKey.set( depth );
1586 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1589 UpSampleData usData[3];
1590 sNodes.treeNodes[i]->depthAndOffset( d , off );
1591 for(
int d=0 ; d<3 ; d++ )
1593 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1594 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1595 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1596 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1598 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1599 TreeOctNode::Neighbors3& neighbors = neighborKey.neighbors[depth-1];
1600 for(
int ii=0 ; ii<2 ; ii++ )
1602 int _ii = ii + usData[0].start;
1603 double dx = usData[0].v[ii];
1604 for(
int jj=0 ; jj<2 ; jj++ )
1606 int _jj = jj + usData[1].start;
1607 double dxy = dx * usData[1].v[jj];
1608 for(
int kk=0 ; kk<2 ; kk++ )
1610 int _kk = kk + usData[2].start;
1611 double dxyz = dxy * usData[2].v[kk];
1612 constraints[t][neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex-lStart] += sNodes.treeNodes[i]->nodeData.constraint * dxyz;
1618#pragma omp parallel for num_threads( threads )
1619 for(
int i=lStart ; i<lEnd ; i++ )
1622 for(
int t=0 ; t<threads ; t++ ) cSum += constraints[t][i-lStart];
1623 sNodes.treeNodes[i]->nodeData.constraint += cSum;
1626 template<
int Degree >
1628 void Octree< Degree >::DownSample(
int depth ,
const SortedTreeNodes& sNodes , C* constraints )
const
1630 if( depth==0 )
return;
1633 for(
int i=sNodes.nodeCount[1] ; i<sNodes.nodeCount[2] ; i++ ) constraints[0] += constraints[i];
1636 std::vector< Vector< C > > _constraints( threads );
1637 for(
int t=0 ; t<threads ; t++ ) _constraints[t].Resize( sNodes.nodeCount[depth] - sNodes.nodeCount[depth-1] );
1638 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start , lStart = sNodes.nodeCount[depth-1] , lEnd = sNodes.nodeCount[depth];
1640#pragma omp parallel for num_threads( threads )
1641 for(
int t=0 ; t<threads ; t++ )
1643 TreeOctNode::NeighborKey3 neighborKey;
1644 neighborKey.set( depth );
1645 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1648 UpSampleData usData[3];
1649 sNodes.treeNodes[i]->depthAndOffset( d , off );
1650 for(
int d=0 ; d<3 ; d++ )
1652 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1653 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1654 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1655 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1657 TreeOctNode::Neighbors3& neighbors = neighborKey.getNeighbors( sNodes.treeNodes[i]->parent );
1658 C c = constraints[i];
1659 for(
int ii=0 ; ii<2 ; ii++ )
1661 int _ii = ii + usData[0].start;
1662 C cx = C( c*usData[0].v[ii] );
1663 for(
int jj=0 ; jj<2 ; jj++ )
1665 int _jj = jj + usData[1].start;
1666 C cxy = C( cx*usData[1].v[jj] );
1667 for(
int kk=0 ; kk<2 ; kk++ )
1669 int _kk = kk + usData[2].start;
1670 if( neighbors.neighbors[_ii][_jj][_kk] )
1671 _constraints[t][neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex-lStart] += C( cxy*usData[2].v[kk] );
1677#pragma omp parallel for num_threads( threads )
1678 for(
int i=lStart ; i<lEnd ; i++ )
1681 for(
int t=0 ; t<threads ; t++ ) cSum += _constraints[t][i-lStart];
1682 constraints[i] += cSum;
1686 template<
int Degree >
1688 void Octree< Degree >::UpSample(
int depth ,
const SortedTreeNodes& sNodes , C* coefficients )
const
1690 if ( depth==0 )
return;
1693 for(
int i=sNodes.nodeCount[1] ; i<sNodes.nodeCount[2] ; i++ ) coefficients[i] += coefficients[0];
1697 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1699#pragma omp parallel for num_threads( threads )
1700 for(
int t=0 ; t<threads ; t++ )
1702 TreeOctNode::NeighborKey3 neighborKey;
1703 neighborKey.set( depth-1 );
1704 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1708 UpSampleData usData[3];
1710 for(
int d=0 ; d<3 ; d++ )
1712 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1713 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1714 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1715 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1717 TreeOctNode::Neighbors3& neighbors = neighborKey.getNeighbors( node->parent );
1718 for(
int ii=0 ; ii<2 ; ii++ )
1720 int _ii = ii + usData[0].start;
1721 double dx = usData[0].v[ii];
1722 for(
int jj=0 ; jj<2 ; jj++ )
1724 int _jj = jj + usData[1].start;
1725 double dxy = dx * usData[1].v[jj];
1726 for(
int kk=0 ; kk<2 ; kk++ )
1728 int _kk = kk + usData[2].start;
1729 if( neighbors.neighbors[_ii][_jj][_kk] )
1731 double dxyz = dxy * usData[2].v[kk];
1732 int _i = neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex;
1733 coefficients[i] += coefficients[_i] *
Real( dxyz );
1742 template<
int Degree >
1743 void Octree< Degree >::SetCoarserPointValues(
int depth ,
const SortedTreeNodes& sNodes ,
Real* metSolution )
1745 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1747#pragma omp parallel for num_threads( threads )
1748 for(
int t=0 ; t<threads ; t++ )
1750 TreeOctNode::NeighborKey3 neighborKey;
1751 neighborKey.set( depth );
1752 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1754 int pIdx = sNodes.treeNodes[i]->nodeData.pointIndex;
1757 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1758 _points[ pIdx ].value = WeightedCoarserFunctionValue( neighborKey , sNodes.treeNodes[i] , metSolution );
1764 template<
int Degree >
1765 Real Octree< Degree >::WeightedCoarserFunctionValue(
const OctNode< TreeNodeData , Real >::NeighborKey3& neighborKey ,
const TreeOctNode* pointNode ,
Real* metSolution )
const
1767 int depth = pointNode->depth();
1768 if( !depth || pointNode->nodeData.pointIndex==-1 )
return Real(0.);
1769 double pointValue = 0;
1771 Real weight = _points[ pointNode->nodeData.pointIndex ].weight;
1772 Point3D< Real > p = _points[ pointNode->nodeData.pointIndex ].position;
1777 const TreeOctNode::Neighbors3& neighbors = neighborKey.neighbors[depth-1];
1778 neighbors.neighbors[1][1][1]->depthAndOffset( d , _idx );
1782 for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
for(
int l=0 ; l<3 ; l++ )
1783 if( neighbors.neighbors[j][k][l] && neighbors.neighbors[j][k][l]->nodeData.nodeIndex>=0 )
1786 const TreeOctNode* basisNode = neighbors.neighbors[j][k][l];
1787 int idx[] = { _idx[0]+j , _idx[1]+k , _idx[2]+l };
1789 fData.baseBSplines[ idx[0] ][2-j]( p[0] ) *
1790 fData.baseBSplines[ idx[1] ][2-k]( p[1] ) *
1791 fData.baseBSplines[ idx[2] ][2-l]( p[2] ) *
1792 metSolution[basisNode->nodeData.nodeIndex];
1795 return Real( pointValue * weight );
1798 template<
int Degree>
1799 int Octree<Degree>::GetFixedDepthLaplacian( SparseSymmetricMatrix< Real >& matrix ,
int depth ,
const SortedTreeNodes& sNodes ,
Real* metSolution )
1801 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1802 double stencil[5][5][5];
1803 SetLaplacianStencil( depth , stencil );
1804 Stencil< double , 5 > stencils[2][2][2];
1805 SetLaplacianStencils( depth , stencils );
1806 matrix.Resize( range );
1807#pragma omp parallel for num_threads( threads )
1808 for(
int t=0 ; t<threads ; t++ )
1810 TreeOctNode::NeighborKey5 neighborKey5;
1811 neighborKey5.set( depth );
1812 for(
int i=(range*t)/threads ; i<(range*(t+1))/threads ; i++ )
1815 neighborKey5.getNeighbors( node );
1818 int count = GetMatrixRowSize( neighborKey5.neighbors[depth] );
1821#pragma omp critical (matrix_set_row_size)
1823 matrix.SetRowSize( i , count );
1827 matrix.rowSizes[i] = SetMatrixRow( neighborKey5.neighbors[depth] , matrix[i] , start , stencil );
1833 c = int( node - node->parent->children );
1837 UpdateConstraintsFromCoarser( neighborKey5 , node , metSolution , stencils[x][y][z] );
1843 template<
int Degree>
1844 int Octree<Degree>::GetRestrictedFixedDepthLaplacian( SparseSymmetricMatrix< Real >& matrix,
int depth,
const int* entries,
int entryCount,
1848 for(
int i=0 ; i<entryCount ; i++ ) sNodes.treeNodes[ entries[i] ]->nodeData.nodeIndex = i;
1849 double stencil[5][5][5];
1850 int rDepth , rOff[3];
1851 rNode->depthAndOffset( rDepth , rOff );
1852 matrix.Resize( entryCount );
1853 SetLaplacianStencil( depth , stencil );
1854 Stencil< double , 5 > stencils[2][2][2];
1855 SetLaplacianStencils( depth , stencils );
1856#pragma omp parallel for num_threads( threads )
1857 for(
int t=0 ; t<threads ; t++ )
1859 TreeOctNode::NeighborKey5 neighborKey5;
1860 neighborKey5.set( depth );
1861 for(
int i=(entryCount*t)/threads ; i<(entryCount*(t+1))/threads ; i++ )
1863 TreeOctNode* node = sNodes.treeNodes[ entries[i] ];
1866 off[0] >>= (depth-rDepth) , off[1] >>= (depth-rDepth) , off[2] >>= (depth-rDepth);
1867 bool isInterior = ( off[0]==rOff[0] && off[1]==rOff[1] && off[2]==rOff[2] );
1869 neighborKey5.getNeighbors( node );
1871 int xStart=0 , xEnd=5 , yStart=0 , yEnd=5 , zStart=0 , zEnd=5;
1872 if( !isInterior ) SetMatrixRowBounds( neighborKey5.neighbors[depth].neighbors[2][2][2] , rDepth , rOff , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1875 int count = GetMatrixRowSize( neighborKey5.neighbors[depth] , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1878#pragma omp critical (matrix_set_row_size)
1880 matrix.SetRowSize( i , count );
1884 matrix.rowSizes[i] = SetMatrixRow( neighborKey5.neighbors[depth] , matrix[i] , 0 , stencil , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1889 c = int( node - node->parent->children );
1893 UpdateConstraintsFromCoarser( neighborKey5 , node , metSolution , stencils[x][y][z] );
1896 for(
int i=0 ; i<entryCount ; i++ ) sNodes.treeNodes[entries[i]]->nodeData.nodeIndex = entries[i];
1900 template<
int Degree>
1905 fData.setDotTables( fData.DD_DOT_FLAG | fData.DV_DOT_FLAG );
1908 _sNodes.treeNodes[0]->nodeData.solution = 0;
1910 std::vector< Real >
metSolution( _sNodes.nodeCount[ _sNodes.maxDepth ] , 0 );
1912 for( i=1 ; i<_sNodes.maxDepth ; i++ )
1918 fData.clearDotTables( fData.VV_DOT_FLAG | fData.DV_DOT_FLAG | fData.DD_DOT_FLAG );
1923 template<
int Degree>
1931 X.Resize(
sNodes.nodeCount[depth+1]-
sNodes.nodeCount[depth] );
1932 if( depth<=_minDepth ) UpSampleCoarserSolution( depth ,
sNodes , X );
1939#pragma omp parallel for num_threads( threads )
1940 for(
int i=_sNodes.nodeCount[depth-1] ; i<_sNodes.nodeCount[depth] ; i++ )
metSolution[i] += _sNodes.treeNodes[i]->nodeData.solution;
1942 if( _constrainValues )
1949 int maxECount = ( (2*Degree+1)*(2*Degree+1)*(2*Degree+1) + 1 ) / 2;
1950 maxECount = ( ( maxECount + 15 ) / 16 ) * 16;
1952 for(
int i=0 ; i<M.rows ; i++ ) M.SetRowSize( i , maxECount );
1958 GetFixedDepthLaplacian( M , depth , sNodes , metSolution );
1965 iter +=
SparseSymmetricMatrix< Real >::Solve( M ,
B , std::max< int >(
int( pow( M.rows , ITERATION_POWER ) ) , minIters ) , X ,
Real(accuracy) , 0 , threads , (depth<=_minDepth) && !_constrainValues );
1970 for(
int i=0 ; i<M.rows ; i++ )
for(
int j=0 ; j<M.rowSizes[i] ; j++ ) mNorm += M[i][j].Value * M[i][j].Value;
1971 double bNorm =
B.Norm( 2 ) , rNorm = (
B - M * X ).Norm( 2 );
1972 printf(
"\tResidual: (%d %g) %g -> %g (%f) [%d]\n" , M.Entries() , sqrt(mNorm) , bNorm , rNorm , rNorm/bNorm , iter );
1981 template<
int Degree>
1982 int Octree<Degree>::SolveFixedDepthMatrix(
int depth ,
const SortedTreeNodes& sNodes ,
Real* metSolution ,
int startingDepth ,
bool showResidual ,
int minIters ,
double accuracy )
1984 if( startingDepth>=depth )
return SolveFixedDepthMatrix( depth , sNodes , metSolution , showResidual , minIters , accuracy );
1986 int i , j , d , tIter=0;
1987 SparseSymmetricMatrix< Real > _M;
1988 Vector< Real >
B , B_ , X_;
1989 AdjacencySetFunction asf;
1990 AdjacencyCountFunction acf;
1991 double systemTime = 0 , solveTime = 0 , memUsage = 0 , evaluateTime = 0 , gTime = 0, sTime = 0;
1992 Real myRadius , myRadius2;
1994 if( depth>_minDepth )
1997 UpSample( depth-1 , sNodes , metSolution );
2000#pragma omp parallel for num_threads( threads )
2001 for(
int i=_sNodes.nodeCount[depth-1] ; i<_sNodes.nodeCount[depth] ; i++ ) metSolution[i] += _sNodes.treeNodes[i]->nodeData.solution;
2004 if( _constrainValues )
2006 SetCoarserPointValues( depth , sNodes , metSolution );
2008 B.Resize( sNodes.nodeCount[depth+1] - sNodes.nodeCount[depth] );
2011 for( i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ )
2013 B[ i-sNodes.nodeCount[depth] ] = sNodes.treeNodes[i]->nodeData.constraint;
2014 sNodes.treeNodes[i]->nodeData.constraint = 0;
2017 myRadius = 2*radius-
Real(0.5);
2020 d = depth-startingDepth;
2021 std::vector< int > subDimension( sNodes.nodeCount[d+1]-sNodes.nodeCount[d] );
2022 int maxDimension = 0;
2023 for( i=sNodes.nodeCount[d] ; i<sNodes.nodeCount[d+1] ; i++ )
2026 acf.adjacencyCount = 0;
2029 if( temp->depth()==depth ) acf.adjacencyCount++ , temp = sNodes.treeNodes[i]->nextBranch( temp );
2030 else temp = sNodes.treeNodes[i]->nextNode ( temp );
2032 for( j=sNodes.nodeCount[d] ; j<sNodes.nodeCount[d+1] ; j++ )
2034 if( i==j )
continue;
2037 subDimension[i-sNodes.nodeCount[d]] = acf.adjacencyCount;
2038 maxDimension = std::max< int >( maxDimension , subDimension[i-sNodes.nodeCount[d]] );
2040 asf.adjacencies =
new int[maxDimension];
2041 MapReduceVector< Real > mrVector;
2042 mrVector.resize( threads , maxDimension );
2044 for( i=sNodes.nodeCount[d] ; i<sNodes.nodeCount[d+1] ; i++ )
2048 acf.adjacencyCount = subDimension[i-sNodes.nodeCount[d]];
2049 if( !acf.adjacencyCount )
continue;
2052 asf.adjacencyCount = 0;
2055 if( temp->depth()==depth ) asf.adjacencies[ asf.adjacencyCount++ ] = temp->nodeData.nodeIndex , temp = sNodes.treeNodes[i]->nextBranch( temp );
2056 else temp = sNodes.treeNodes[i]->nextNode ( temp );
2058 for( j=sNodes.nodeCount[d] ; j<sNodes.nodeCount[d+1] ; j++ )
2060 if( i==j )
continue;
2065 B_.Resize( asf.adjacencyCount );
2066 for( j=0 ; j<asf.adjacencyCount ; j++ ) B_[j] =
B[ asf.adjacencies[j]-sNodes.nodeCount[depth] ];
2068 X_.Resize( asf.adjacencyCount );
2069#pragma omp parallel for num_threads( threads ) schedule( static )
2070 for( j=0 ; j<asf.adjacencyCount ; j++ )
2072 X_[j] = sNodes.treeNodes[ asf.adjacencies[j] ]->nodeData.solution;
2076 GetRestrictedFixedDepthLaplacian( _M , depth , asf.adjacencies , asf.adjacencyCount , sNodes.treeNodes[i] , myRadius , sNodes , metSolution );
2077#pragma omp parallel for num_threads( threads ) schedule( static )
2078 for( j=0 ; j<asf.adjacencyCount ; j++ )
2080 B_[j] += sNodes.treeNodes[asf.adjacencies[j]]->nodeData.constraint;
2081 sNodes.treeNodes[ asf.adjacencies[j] ]->nodeData.constraint = 0;
2091 for(
int i=0 ; i<_M.rows ; i++ )
for(
int j=0 ; j<_M.rowSizes[i] ; j++ ) mNorm += _M[i][j].Value * _M[i][j].Value;
2092 double bNorm = B_.Norm( 2 ) , rNorm = ( B_ - _M * X_ ).Norm( 2 );
2093 printf(
"\t\tResidual: (%d %g) %g -> %g (%f) [%d]\n" , _M.Entries() , sqrt(mNorm) , bNorm , rNorm , rNorm/bNorm , iter );
2097 for( j=0 ; j<asf.adjacencyCount ; j++ )
2099 TreeOctNode* temp=sNodes.treeNodes[ asf.adjacencies[j] ];
2100 while( temp->depth()>sNodes.treeNodes[i]->depth() ) temp=temp->
parent;
2101 if( temp->nodeData.nodeIndex>=sNodes.treeNodes[i]->nodeData.nodeIndex ) sNodes.treeNodes[ asf.adjacencies[j] ]->
nodeData.solution =
Real( X_[j] );
2103 systemTime += gTime;
2105 memUsage = std::max< double >( MemoryUsage() , memUsage );
2108 delete[] asf.adjacencies;
2112 template<
int Degree>
2116 if( node->nodeData.normalIndex>=0 && ( (*normals)[node->nodeData.normalIndex][0]!=0 || (*normals)[node->nodeData.normalIndex][1]!=0 || (*normals)[node->nodeData.normalIndex][2]!=0 ) ) hasNormals=1;
2117 if( node->children )
for(
int i=0;i<
Cube::CORNERS && !hasNormals;i++) hasNormals |= HasNormals(&node->children[i],epsilon);
2122 template<
int Degree>
2127 if(
temp->children &&
temp->d>=_minDepth )
2136 template<
int Degree>
2144 fData.setDotTables( fData.VV_DOT_FLAG | fData.DV_DOT_FLAG );
2150 std::vector< Point3D< Real > > coefficients( _sNodes.nodeCount[
maxDepth] ,
zeroPoint );
2153#pragma omp parallel for num_threads( threads )
2154 for(
int i=0 ; i<_sNodes.nodeCount[
maxDepth+1] ; i++ ) _sNodes.treeNodes[i]->nodeData.constraint =
Real( 0. );
2157 std::vector< std::vector< Real > >
_constraints( threads );
2163 SetDivergenceStencil( d , &
stencil[0][0][0] ,
false );
2165 SetDivergenceStencils( d ,
stencils ,
true );
2166#pragma omp parallel for num_threads( threads )
2167 for(
int t=0 ;
t<threads ;
t++ )
2171 int start = _sNodes.nodeCount[d] , end = _sNodes.nodeCount[d+1] , range = end-start;
2172 for(
int i=start+(range*
t)/threads ; i<start+(range*(
t+1))/threads ; i++ )
2176 int depth = node->
depth();
2183 int mn = 2 ,
mx = (1<<d)-2;
2194 else cx = cy =
cz = 0;
2205 if(
_node &&
_node->nodeData.normalIndex>=0 )
2215 if(
_node &&
_node->nodeData.normalIndex>=0 )
2225 if( normal[0]==0 && normal[1]==0 && normal[2]==0 )
continue;
2226 if( depth<
maxDepth ) coefficients[i] += normal;
2247#pragma omp parallel for num_threads( threads ) schedule( static )
2248 for(
int i=0 ; i<_sNodes.nodeCount[
maxDepth] ; i++ )
2258 for(
int d=0 ; d<
maxDepth ; d++ ) UpSample( d , _sNodes , &coefficients[0] );
2261#pragma omp parallel for num_threads( threads )
2262 for(
int i=0 ; i<_sNodes.nodeCount[
maxDepth] ; i++ ) _sNodes.treeNodes[i]->nodeData.constraint +=
constraints[i];
2267 int start = _sNodes.nodeCount[d] , end = _sNodes.nodeCount[d+1] , range = end - start;
2269 SetDivergenceStencils( d ,
stencils ,
false );
2270#pragma omp parallel for num_threads( threads )
2271 for(
int t=0 ;
t<threads ;
t++ )
2275 for(
int i=start+(range*
t)/threads ; i<start+(range*(
t+1))/threads ; i++ )
2278 int depth = node->
depth();
2279 if( !depth )
continue;
2288 int mn = 4 ,
mx = (1<<d)-4;
2297 else cx = cy =
cz = 0;
2305 int _i =
_node->nodeData.nodeIndex;
2310 constraint +=
Real(
div[0] * normal[0] +
div[1] * normal[1] +
div[2] * normal[2] );
2312 else constraint += GetDivergence(
_node , node , coefficients[
_i] );
2314 node->
nodeData.constraint += constraint;
2319 fData.clearDotTables( fData.DV_DOT_FLAG );
2322#pragma omp parallel for num_threads( threads )
2323 for(
int t=0 ;
t<threads ;
t++ )
2324 for(
int i=(_sNodes.nodeCount[
maxDepth+1]*
t)/threads ; i<(_sNodes.nodeCount[
maxDepth+1]*(
t+1))/threads ; i++ )
2327 if(
temp->nodeData.nodeIndex<0 ||
temp->nodeData.normalIndex<0 )
temp->nodeData.centerWeightContribution = 0;
2328 else temp->nodeData.centerWeightContribution =
Real(
Length((*normals)[
temp->nodeData.normalIndex]) );
2335 template<
int Degree>
2338 template<
int Degree>
2341 template<
int Degree>
2344 if( !node1->children && node1->depth()<depth ) node1->initChildren();
2347 template<
int Degree >
2348 void Octree< Degree >::FaceEdgesFunction::Function(
const TreeOctNode* node1 ,
const TreeOctNode* node2 )
2356 for(
int j=0 ; j<count ; j++ )
2357 for(
int k=0 ; k<3 ; k++ )
2359 if( GetRootIndex( node1 , isoTri[j*3+k] ,
maxDepth , ri1 ) && GetRootIndex( node1 , isoTri[j*3+((k+1)%3)] ,
maxDepth , ri2 ) )
2361 long long key1=ri1.key , key2=ri2.key;
2362 edges->push_back( std::pair< RootInfo , RootInfo >( ri2 , ri1 ) );
2363 if( vertexCount->count( key1 )==0 )
2365 (*vertexCount)[key1].first = ri1;
2366 (*vertexCount)[key1].second=0;
2368 if( vertexCount->count( key2 )==0 )
2370 (*vertexCount)[key2].first = ri2;
2371 (*vertexCount)[key2].second=0;
2373 (*vertexCount)[key1].second--;
2374 (*vertexCount)[key2].second++;
2376 else fprintf( stderr ,
"Bad Edge 1: %d %d\n" , ri1.key , ri2.key );
2380 template<
int Degree >
2390 bool flags[3][3][3];
2400 TreeOctNode::NeighborKey3
nKey;
2405 int d , off[3] , _off[3];
2406 leaf->depthAndOffset( d , off );
2407 int res = (1<<d)-1 ,
_res = ( 1<<(d-
sDepth) )-1;
2408 _off[0] = off[0]&
_res , _off[1] = off[1]&
_res , _off[2] = off[2]&
_res;
2411 { (off[0]!=0 && _off[0]==0) , (off[0]!=res && _off[0]==
_res) } ,
2412 { (off[1]!=0 && _off[1]==0) , (off[1]!=res && _off[1]==
_res) } ,
2413 { (off[2]!=0 && _off[2]==0) , (off[2]!=res && _off[2]==
_res) }
2418 TreeOctNode::Neighbors3& neighbors =
nKey.getNeighbors(
leaf );
2419 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ ) flags[i][j][k] =
false;
2420 int x=0 , y=0 , z=0;
2421 if (
boundary[0][0] && !neighbors.neighbors[0][1][1] ) x = -1;
2422 else if(
boundary[0][1] && !neighbors.neighbors[2][1][1] ) x = 1;
2423 if (
boundary[1][0] && !neighbors.neighbors[1][0][1] ) y = -1;
2424 else if(
boundary[1][1] && !neighbors.neighbors[1][2][1] ) y = 1;
2425 if (
boundary[2][0] && !neighbors.neighbors[1][1][0] ) z = -1;
2426 else if(
boundary[2][1] && !neighbors.neighbors[1][1][2] ) z = 1;
2431 if( x && y && z ) flags[1+x][1+y][1+z] =
true;
2433 if( x && y ) flags[1+x][1+y][1 ] =
true;
2434 if( x && z ) flags[1+x][1 ][1+z] =
true;
2435 if( y && z ) flags[1 ][1+y][1+1] =
true;
2437 if( x ) flags[1+x][1 ][1 ] =
true;
2438 if( y ) flags[1 ][1+y][1 ] =
true;
2439 if( z ) flags[1 ][1 ][1+z] =
true;
2444 _sNodes.set( tree );
2448 template<
int Degree>
2451 fData.setValueTables( fData.VALUE_FLAG | fData.D_VALUE_FLAG , 0 , postNormalSmooth );
2463#pragma omp parallel for num_threads( threads )
2464 for(
int i=_sNodes.nodeCount[_minDepth] ; i<_sNodes.nodeCount[
maxDepth] ; i++ )
metSolution[i] = _sNodes.treeNodes[i]->nodeData.solution;
2468#pragma omp parallel for num_threads( threads )
2469 for(
int i=0 ; i<_sNodes.nodeCount[
maxDepth+1] ; i++ ) _sNodes.treeNodes[i]->nodeData.mcIndex = 0;
2471 rootData.boundaryValues =
new std::unordered_map< long long , std::pair< Real , Point3D< Real > > >();
2491 std::vector< TreeOctNode::ConstNeighborKey3 >
nKeys( threads );
2493 TreeOctNode::ConstNeighborKey3
nKey;
2494 std::vector< TreeOctNode::ConstNeighborKey5 >
nKeys5( threads );
2496 TreeOctNode::ConstNeighborKey5
nKey5;
2500 for(
int i=_sNodes.nodeCount[
sDepth] ; i<_sNodes.nodeCount[
sDepth+1] ; i++ )
2502 if( !_sNodes.treeNodes[i]->children )
continue;
2504 _sNodes.setCornerTable(
rootData , _sNodes.treeNodes[i] , threads );
2505 _sNodes.setEdgeTable (
rootData , _sNodes.treeNodes[i] , threads );
2516 for(
TreeOctNode* node=_sNodes.treeNodes[i]->
nextLeaf() ; node ; node=_sNodes.treeNodes[i]->nextLeaf( node ) )
if( node->d==d )
leafNodes.push_back( node );
2521#pragma omp parallel for num_threads( threads )
2529 leaf->depthAndOffset( d , off );
2531 off[0] %= res , off[1] %=res , off[2] %= res;
2533 if( !(off[0]%res) && !(off[1]%res) && !(off[2]%res) )
2537 int x = off[0]==0 ? 0 : 1 , y = off[1]==0 ? 0 : 1 , z = off[2]==0 ? 0 : 1;
2555#pragma omp parallel for num_threads( threads )
2582 for(
auto iter=
rootData.boundaryRoots.cbegin() ; iter!=
rootData.boundaryRoots.cend() ; iter++ )
2585 for(
int d=
sDepth ; d>=0 ; d-- )
2593 for(
int i=_sNodes.nodeCount[d] ; i<_sNodes.nodeCount[d+1] ; i++ )
2596 if(
leaf->children )
continue;
2620 template<
int Degree>
2626 idx[0]*=fData.functionCount;
2627 idx[1]*=fData.functionCount;
2628 idx[2]*=fData.functionCount;
2629 int minDepth = std::max< int >( 0 , std::min< int >( _minDepth , node->
depth()-1 ) );
2631 for(
int j=0;j<3;j++)
2632 for(
int k=0;k<3;k++)
2633 for(
int l=0;l<3;l++)
2640 fData.valueTables[idx[0]+
int(n->
off[0])]*
2641 fData.valueTables[idx[1]+
int(n->
off[1])]*
2642 fData.valueTables[idx[2]+
int(n->
off[2])]);
2652 fData.valueTables[idx[0]+
int(n->off[0])]*
2653 fData.valueTables[idx[1]+
int(n->off[1])]*
2654 fData.valueTables[idx[2]+
int(n->off[2])]);
2655 if( n->children ) n=&n->children[ii];
2663 template<
int Degree >
2664 Real Octree< Degree >::getCornerValue(
const OctNode< TreeNodeData , Real >::ConstNeighborKey3& neighborKey3 ,
const TreeOctNode* node ,
int corner ,
const Real* metSolution )
2670 idx[0] *= fData.functionCount;
2671 idx[1] *= fData.functionCount;
2672 idx[2] *= fData.functionCount;
2674 int d = node->depth();
2676 int startX = 0 , endX = 3 , startY = 0 , endY = 3 , startZ = 0 , endZ = 3;
2679 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d];
2680 if( cx==0 ) endX = 2;
2682 if( cy==0 ) endY = 2;
2684 if( cz==0 ) endZ = 2;
2686 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2688 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2692 fData.valueTables[ idx[0]+int(n->off[0]) ]*
2693 fData.valueTables[ idx[1]+int(n->off[1]) ]*
2694 fData.valueTables[ idx[2]+int(n->off[2]) ];
2695 value += n->nodeData.solution * v;
2699 if( d>0 && d>_minDepth )
2701 int _corner = int( node - node->parent->children );
2702 int _cx , _cy , _cz;
2704 if( cx!=_cx ) startX = 0 , endX = 3;
2705 if( cy!=_cy ) startY = 0 , endY = 3;
2706 if( cz!=_cz ) startZ = 0 , endZ = 3;
2707 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d-1];
2708 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2710 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2714 fData.valueTables[ idx[0]+int(n->off[0]) ]*
2715 fData.valueTables[ idx[1]+int(n->off[1]) ]*
2716 fData.valueTables[ idx[2]+int(n->off[2]) ];
2717 value += metSolution[ n->nodeData.nodeIndex ] * v;
2721 return Real( value );
2724 template<
int Degree >
2725 Real Octree< Degree >::getCornerValue(
const OctNode< TreeNodeData , Real >::ConstNeighborKey3& neighborKey3 ,
const TreeOctNode* node ,
int corner ,
const Real* metSolution ,
const double stencil1[3][3][3] ,
const double stencil2[3][3][3] )
2728 int d = node->depth();
2730 int startX = 0 , endX = 3 , startY = 0 , endY = 3 , startZ = 0 , endZ = 3;
2733 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d];
2734 if( cx==0 ) endX = 2;
2736 if( cy==0 ) endY = 2;
2738 if( cz==0 ) endZ = 2;
2740 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2742 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2743 if( n ) value += n->
nodeData.solution * stencil1[x][y][z];
2746 if( d>0 && d>_minDepth )
2748 int _corner = int( node - node->parent->children );
2749 int _cx , _cy , _cz;
2751 if( cx!=_cx ) startX = 0 , endX = 3;
2752 if( cy!=_cy ) startY = 0 , endY = 3;
2753 if( cz!=_cz ) startZ = 0 , endZ = 3;
2754 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d-1];
2755 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2757 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2758 if( n ) value += metSolution[ n->nodeData.nodeIndex ] * stencil2[x][y][z];
2761 return Real( value );
2764 template<
int Degree >
2765 Point3D< Real > Octree< Degree >::getCornerNormal(
const OctNode< TreeNodeData , Real >::ConstNeighborKey5& neighborKey5 ,
const TreeOctNode* node ,
int corner ,
const Real* metSolution )
2768 Point3D< Real > normal;
2769 normal[0] = normal[1] = normal[2] = 0.;
2772 idx[0] *= fData.functionCount;
2773 idx[1] *= fData.functionCount;
2774 idx[2] *= fData.functionCount;
2776 int d = node->depth();
2779 TreeOctNode::ConstNeighbors5& neighbors = neighborKey5.neighbors[d];
2780 for(
int j=0 ; j<5 ; j++ )
for(
int k=0 ; k<5 ; k++ )
for(
int l=0 ; l<5 ; l++ )
2782 const TreeOctNode* n=neighbors.neighbors[j][k][l];
2785 int _idx[] = { idx[0] + n->
off[0] , idx[1] + n->off[1] , idx[2] + n->off[2] };
2786 double values[] = { fData.valueTables[_idx[0]] , fData.valueTables[_idx[1]] , fData.valueTables[_idx[2]] };
2787 double dValues[] = { fData.dValueTables[_idx[0]] , fData.dValueTables[_idx[1]] , fData.dValueTables[_idx[2]] };
2788 Real solution = n->nodeData.solution;
2789 normal[0] +=
Real( dValues[0] * values[1] * values[2] * solution );
2790 normal[1] +=
Real( values[0] * dValues[1] * values[2] * solution );
2791 normal[2] +=
Real( values[0] * values[1] * dValues[2] * solution );
2795 if( d>0 && d>_minDepth )
2797 TreeOctNode::ConstNeighbors5& neighbors = neighborKey5.neighbors[d-1];
2798 for(
int j=0 ; j<5 ; j++ )
for(
int k=0 ; k<5 ; k++ )
for(
int l=0 ; l<5 ; l++ )
2800 const TreeOctNode* n=neighbors.neighbors[j][k][l];
2803 int _idx[] = { idx[0] + n->
off[0] , idx[1] + n->off[1] , idx[2] + n->off[2] };
2804 double values[] = { fData.valueTables[_idx[0]] , fData.valueTables[_idx[1]] , fData.valueTables[_idx[2]] };
2805 double dValues[] = { fData.dValueTables[_idx[0]] , fData.dValueTables[_idx[1]] , fData.dValueTables[_idx[2]] };
2806 Real solution = metSolution[ n->nodeData.nodeIndex ];
2807 normal[0] +=
Real( dValues[0] * values[1] * values[2] * solution );
2808 normal[1] +=
Real( values[0] * dValues[1] * values[2] * solution );
2809 normal[2] +=
Real( values[0] * values[1] * dValues[2] * solution );
2816 template<
int Degree >
2821 neighborKey2.set( fData.depth );
2822 fData.setValueTables( fData.VALUE_FLAG , 0 );
2825#pragma omp parallel for num_threads( threads ) reduction( + : isoValue , weightSum )
2826 for(
int t=0 ;
t<threads ;
t++)
2828 TreeOctNode::ConstNeighborKey3
nKey;
2829 nKey.set( _sNodes.maxDepth-1 );
2830 int nodeCount = _sNodes.nodeCount[ _sNodes.maxDepth ];
2835 Real w =
temp->nodeData.centerWeightContribution;
2846 template<
int Degree >
2847 void Octree< Degree >::SetIsoCorners(
Real isoValue ,
TreeOctNode*
leaf ,
SortedTreeNodes::CornerTableData&
cData ,
char*
valuesSet ,
Real* values , TreeOctNode::ConstNeighborKey3&
nKey ,
const Real*
metSolution ,
const Stencil< double , 3 > stencil1[8] ,
const Stencil< double , 3 > stencil2[8][8] )
2854 leaf->depthAndOffset( d , off );
2855 int mn = 2 ,
mx = (1<<d)-2;
2866 values[
vIndex] = cornerValues[
c];
2886 if( parent->parent && parent->parent->d>=_minDepth && (parent-parent->parent->children)==c )
2889 parent = parent->parent;
2897 template<
int Degree>
2898 int Octree<Degree>::InteriorFaceRootCount(
const TreeOctNode* node,
const int &faceIndex,
int maxDepth){
2899 int c1,c2,e1,e2,dir,off,cnt=0;
2920 cnt+=EdgeRootCount(&node->children[c1],e1,
maxDepth)+EdgeRootCount(&node->children[c1],e2,
maxDepth);
2935 cnt+=EdgeRootCount(&node->children[c2],e1,
maxDepth)+EdgeRootCount(&node->children[c2],e2,
maxDepth);
2936 for(
int i=0;i<
Cube::CORNERS/2;i++){
if(node->children[corners[i]].children){cnt+=InteriorFaceRootCount(&node->children[corners[i]],faceIndex,
maxDepth);}}
2941 template<
int Degree>
2952 if(temp && temp->children){
2957 temp=node->faceNeighbor(f2);
2958 if(temp && temp->children){
2963 temp=node->edgeNeighbor(edgeIndex);
2964 if(temp && temp->children){
2973 if(finest->children)
return EdgeRootCount(&finest->children[c1],eIndex,
maxDepth)+EdgeRootCount(&finest->children[c2],eIndex,
maxDepth);
2977 template<
int Degree>
2978 int Octree<Degree>::IsBoundaryFace(
const TreeOctNode* node,
int faceIndex,
int subdivideDepth){
2979 int dir,offset,d,o[3],idx;
2981 if(subdivideDepth<0){
return 0;}
2982 if(node->d<=subdivideDepth){
return 1;}
2984 node->depthAndOffset(d,o);
2986 idx=(int(o[dir])<<1) + (offset<<1);
2987 return !(idx%(2<<(int(node->d)-subdivideDepth)));
2990 template<
int Degree>
2991 int Octree<Degree>::IsBoundaryEdge(
const TreeOctNode* node,
int edgeIndex,
int subdivideDepth){
2994 return IsBoundaryEdge(node,dir,x,y,subdivideDepth);
2997 template<
int Degree>
2998 int Octree<Degree>::IsBoundaryEdge(
const TreeOctNode* node ,
int dir ,
int x ,
int y ,
int subdivideDepth )
3000 int d , o[3] , idx1 , idx2 , mask;
3002 if( subdivideDepth<0 )
return 0;
3003 if( node->d<=subdivideDepth )
return 1;
3004 node->depthAndOffset( d , o );
3021 mask = 1<<( int(node->d) - subdivideDepth );
3022 return !(idx1%(mask)) || !(idx2%(mask));
3025 template<
int Degree >
3033 ri.node->centerAndWidth( c , width );
3037 start[0] = c[0] - width/2;
3038 end [0] = c[0] + width/2;
3039 start[1] = end[1] = c[1] - width/2 + width*i1;
3040 start[2] = end[2] = c[2] - width/2 + width*i2;
3043 start[0] = end[0] = c[0] - width/2 + width*i1;
3044 start[1] = c[1] - width/2;
3045 end [1] = c[1] + width/2;
3046 start[2] = end[2] = c[2] - width/2 + width*i2;
3049 start[0] = end[0] = c[0] - width/2 + width*i1;
3050 start[1] = end[1] = c[1] - width/2 + width*i2;
3051 start[2] = c[2] - width/2;
3052 end [2] = c[2] + width/2;
3060 template<
int Degree >
3061 int Octree< Degree >::GetRoot(
const RootInfo& ri ,
Real isoValue , TreeOctNode::ConstNeighborKey5& neighborKey5 , Point3D< Real > & position , RootData& rootData ,
int sDepth ,
const Real* metSolution ,
int nonLinearFit )
3068 long long key1 , key2;
3069 Point3D< Real > n[2];
3071 int i , o , i1 , i2 , rCount=0;
3073 std::vector< double > roots;
3075 Real center , width;
3078 int idx1[3] , idx2[3];
3082 bool isBoundary = ( IsBoundaryEdge( ri.node , ri.edgeIndex , sDepth )!=0 );
3083 bool haveKey1 , haveKey2;
3084 std::pair< Real , Point3D< Real > > keyValue1 , keyValue2;
3087 iter1 = rootData.cornerIndices( ri.node )[ c1 ];
3088 iter2 = rootData.cornerIndices( ri.node )[ c2 ];
3089 keyValue1.first = rootData.cornerValues[iter1];
3090 keyValue2.first = rootData.cornerValues[iter2];
3093#pragma omp critical (normal_hash_access)
3095 haveKey1 = ( rootData.boundaryValues->count( key1 )>0 );
3096 haveKey2 = ( rootData.boundaryValues->count( key2 )>0 );
3097 if( haveKey1 ) keyValue1 = (*rootData.boundaryValues)[key1];
3098 if( haveKey2 ) keyValue2 = (*rootData.boundaryValues)[key2];
3103 haveKey1 = ( rootData.cornerNormalsSet[ iter1 ] != 0 );
3104 haveKey2 = ( rootData.cornerNormalsSet[ iter2 ] != 0 );
3105 keyValue1.first = rootData.cornerValues[iter1];
3106 keyValue2.first = rootData.cornerValues[iter2];
3107 if( haveKey1 ) keyValue1.second = rootData.cornerNormals[iter1];
3108 if( haveKey2 ) keyValue2.second = rootData.cornerNormals[iter2];
3111 if( !haveKey1 || !haveKey2 ) neighborKey5.getNeighbors( ri.node );
3112 if( !haveKey1 ) keyValue1.second = getCornerNormal( neighborKey5 , ri.node , c1 , metSolution );
3113 x0 = keyValue1.first;
3114 n[0] = keyValue1.second;
3116 if( !haveKey2 ) keyValue2.second = getCornerNormal( neighborKey5 , ri.node , c2 , metSolution );
3117 x1 = keyValue2.first;
3118 n[1] = keyValue2.second;
3120 if( !haveKey1 || !haveKey2 )
3124#pragma omp critical (normal_hash_access)
3126 if( !haveKey1 ) (*rootData.boundaryValues)[key1] = keyValue1;
3127 if( !haveKey2 ) (*rootData.boundaryValues)[key2] = keyValue2;
3132 if( !haveKey1 ) rootData.cornerNormals[iter1] = keyValue1.second , rootData.cornerNormalsSet[ iter1 ] = 1;
3133 if( !haveKey2 ) rootData.cornerNormals[iter2] = keyValue2.second , rootData.cornerNormalsSet[ iter2 ] = 1;
3138 ri.node->centerAndWidth(c,width);
3140 for( i=0 ; i<DIMENSION ; i++ ) n[0][i] *= width , n[1][i] *= width;
3145 position[1] = c[1]-width/2+width*i1;
3146 position[2] = c[2]-width/2+width*i2;
3149 position[0] = c[0]-width/2+width*i1;
3150 position[2] = c[2]-width/2+width*i2;
3153 position[0] = c[0]-width/2+width*i1;
3154 position[1] = c[1]-width/2+width*i2;
3162 double scl=(x1-x0)/((dx1+dx0)/2);
3167 P.coefficients[0] = x0;
3168 P.coefficients[1] = dx0;
3169 P.coefficients[2] = 3*(x1-x0)-dx1-2*dx0;
3171 P.getSolutions( isoValue , roots ,
EPSILON );
3172 for( i=0 ; i<int(roots.size()) ; i++ )
3173 if( roots[i]>=0 && roots[i]<=1 )
3175 averageRoot +=
Real( roots[i] );
3178 if( rCount && nonLinearFit ) averageRoot /= rCount;
3179 else averageRoot =
Real((x0-isoValue)/(x0-x1));
3180 if( averageRoot<0 || averageRoot>1 )
3182 fprintf( stderr ,
"[WARNING] Bad average root: %f\n" , averageRoot );
3183 fprintf( stderr ,
"\t(%f %f) , (%f %f) (%f)\n" , x0 , x1 , dx0 , dx1 , isoValue );
3184 if( averageRoot<0 ) averageRoot = 0;
3185 if( averageRoot>1 ) averageRoot = 1;
3187 position[o] =
Real(center-width/2+width*averageRoot);
3191 template<
int Degree >
3192 int Octree< Degree >::GetRootIndex(
const TreeOctNode* node ,
int edgeIndex ,
int maxDepth ,
int sDepth,RootInfo& ri )
3201 finestIndex=edgeIndex;
3203 if(IsBoundaryFace(node,f1,sDepth)){temp=NULL;}
3204 else{temp=node->faceNeighbor(f1);}
3205 if(temp && temp->children){
3210 if(IsBoundaryFace(node,f2,sDepth)){temp=NULL;}
3211 else{temp=node->faceNeighbor(f2);}
3212 if(temp && temp->children){
3217 if(IsBoundaryEdge(node,edgeIndex,sDepth)){temp=NULL;}
3218 else{temp=node->edgeNeighbor(edgeIndex);}
3219 if(temp && temp->children){
3228 if(finest->children){
3229 if (GetRootIndex(&finest->children[c1],finestIndex,
maxDepth,sDepth,ri)) {
return 1;}
3230 else if (GetRootIndex(&finest->children[c2],finestIndex,
maxDepth,sDepth,ri)) {
return 1;}
3233 fprintf( stderr ,
"[WARNING] Couldn't find root index with either child\n" );
3241 fprintf( stderr ,
"[WARNING] Finest node does not have iso-edge\n" );
3248 finest->depthAndOffset(d,off);
3250 ri.edgeIndex=finestIndex;
3251 int eIndex[2],offset;
3268 ri.key = (
long long)(o) | (
long long)(eIndex[0])<<5 | (
long long)(eIndex[1])<<25 | (
long long)(offset)<<45;
3273 template<
int Degree>
3274 int Octree<Degree>::GetRootIndex(
const TreeOctNode* node ,
int edgeIndex ,
int maxDepth , RootInfo& ri )
3287 finestIndex = edgeIndex;
3288 if( node->depth()<
maxDepth && !node->children )
3290 temp=node->faceNeighbor( f1 );
3294 temp = node->faceNeighbor( f2 );
3298 temp = node->edgeNeighbor( edgeIndex );
3305 if( finest->children )
3307 if ( GetRootIndex( finest->children + c1 , finestIndex ,
maxDepth , ri ) )
return 1;
3308 else if( GetRootIndex( finest->children + c2 , finestIndex ,
maxDepth , ri ) )
return 1;
3311 int d1 , off1[3] , d2 , off2[3];
3312 node->depthAndOffset( d1 , off1 );
3313 finest->depthAndOffset( d2 , off2 );
3314 fprintf( stderr ,
"[WARNING] Couldn't find root index with either child [%d] (%d %d %d) -> [%d] (%d %d %d) (%d %d)\n" , d1 , off1[0] , off1[1] , off1[2] , d2 , off2[0] , off2[1] , off2[2] , node->children!=NULL , finest->children!=NULL );
3316 for(
int i=0 ; i<8 ; i++ )
if( node->nodeData.mcIndex & (1<<i) ) printf(
"1" );
else printf(
"0" );
3318 for(
int i=0 ; i<8 ; i++ )
if( finest->nodeData.mcIndex & (1<<i) ) printf(
"1" );
else printf(
"0" );
3328 finest->depthAndOffset(d,off);
3330 ri.edgeIndex=finestIndex;
3331 int offset,eIndex[2];
3347 ri.key= (
long long)(o) | (
long long)(eIndex[0])<<5 | (
long long)(eIndex[1])<<25 | (
long long)(offset)<<45;
3352 template<
int Degree>
3353 int Octree<Degree>::GetRootPair(
const RootInfo& ri,
int maxDepth,RootInfo& pair){
3357 while(node->parent){
3358 c=int(node-node->parent->children);
3359 if(c!=c1 && c!=c2){
return 0;}
3361 if(c==c1){
return GetRootIndex(&node->parent->children[c2],ri.edgeIndex,
maxDepth,pair);}
3362 else{
return GetRootIndex(&node->parent->children[c1],ri.edgeIndex,
maxDepth,pair);}
3369 template<
int Degree>
3370 int Octree< Degree >::GetRootIndex(
const RootInfo& ri , RootData& rootData , CoredPointIndex& index )
3372 long long key = ri.key;
3373 auto rootIter = rootData.boundaryRoots.find( key );
3374 if( rootIter!=rootData.boundaryRoots.end() )
3377 index.index = rootIter->second;
3380 else if( rootData.interiorRoots )
3382 int eIndex = rootData.edgeIndices( ri.node )[ ri.edgeIndex ];
3383 if( rootData.edgesSet[ eIndex ] )
3386 index.index = rootData.interiorRoots[ eIndex ];
3393 template<
int Degree >
3394 int Octree< Degree >::SetMCRootPositions(
TreeOctNode* node ,
int sDepth ,
Real isoValue , TreeOctNode::ConstNeighborKey5& neighborKey5 , RootData& rootData ,
3395 std::vector<
Point3D< float > >* interiorPositions , CoredMeshData* mesh ,
const Real* metSolution ,
int nonLinearFit )
3397 Point3D< Real > position;
3402 for(
int i=0 ; i<DIMENSION ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
3406 if( GetRootIndex( node , eIndex , fData.depth , ri ) )
3409 if( !rootData.interiorRoots || IsBoundaryEdge( node , i , j , k , sDepth ) )
3411 std::unordered_map< long long , int >::iterator iter , end;
3413#pragma omp critical (boundary_roots_hash_access)
3415 iter = rootData.boundaryRoots.find( key );
3416 end = rootData.boundaryRoots.end();
3421 GetRoot( ri , isoValue , neighborKey5 , position , rootData , sDepth , metSolution , nonLinearFit );
3423#pragma omp critical (boundary_roots_hash_access)
3425 iter = rootData.boundaryRoots.find( key );
3426 end = rootData.boundaryRoots.end();
3429 mesh->inCorePoints.push_back( position );
3430 rootData.boundaryRoots[key] = int( mesh->inCorePoints.size() ) - 1;
3433 if( iter==end ) count++;
3438 int nodeEdgeIndex = rootData.edgeIndices( ri.node )[ ri.edgeIndex ];
3439 if( !rootData.edgesSet[ nodeEdgeIndex ] )
3442 GetRoot( ri , isoValue , neighborKey5 , position , rootData , sDepth , metSolution , nonLinearFit );
3444#pragma omp critical (add_point_access)
3446 if( !rootData.edgesSet[ nodeEdgeIndex ] )
3448 rootData.interiorRoots[ nodeEdgeIndex ] = mesh->addOutOfCorePoint( position );
3449 interiorPositions->push_back( position );
3450 rootData.edgesSet[ nodeEdgeIndex ] = 1;
3461 template<
int Degree>
3462 int Octree< Degree >::SetBoundaryMCRootPositions(
int sDepth ,
Real isoValue , RootData& rootData , CoredMeshData* mesh ,
int nonLinearFit )
3464 Point3D< Real > position;
3465 int i,j,k,eIndex,hits=0;
3476 for( i=0 ; i<DIMENSION ; i++ )
3477 for( j=0 ; j<2 ; j++ )
3478 for( k=0 ; k<2 ; k++ )
3479 if( IsBoundaryEdge( node , i , j , k , sDepth ) )
3484 if( GetRootIndex( node , eIndex , fData.depth , ri ) )
3487 if( rootData.boundaryRoots.find(key)==rootData.boundaryRoots.end() )
3489 GetRoot( ri , isoValue , position , rootData , sDepth , nonLinearFit );
3490 mesh->inCorePoints.push_back( position );
3491 rootData.boundaryRoots[key] = int( mesh->inCorePoints.size() )-1;
3497 if( hits ) node=tree.nextLeaf(node);
3498 else node=tree.nextBranch(node);
3503 template<
int Degree>
3504 void Octree< Degree >::GetMCIsoEdges(
TreeOctNode* node ,
int sDepth , std::vector< std::pair< RootInfo , RootInfo > >& edges )
3507 int count=0 , tris=0;
3509 FaceEdgesFunction fef;
3511 std::unordered_map< long long , std::pair< RootInfo , int > > vertexCount;
3514 fef.maxDepth = fData.depth;
3515 fef.vertexCount = &vertexCount;
3521 temp = node->faceNeighbor( fIndex );
3524 if( temp && temp->children && !IsBoundaryFace( node , fIndex , sDepth ) ) temp->processNodeFaces( temp , &fef , ref );
3529 for(
int j=0 ; j<count ; j++ )
3530 for(
int k=0 ; k<3 ; k++ )
3532 if( GetRootIndex( node , isoTri[j*3+k] , fData.depth , ri1 ) && GetRootIndex( node , isoTri[j*3+((k+1)%3)] , fData.depth , ri2 ) )
3534 long long key1 = ri1.key , key2 = ri2.key;
3535 edges.push_back( std::pair< RootInfo , RootInfo >( ri1 , ri2 ) );
3536 if( vertexCount.count( key1 )==0 )
3538 vertexCount[key1].first = ri1;
3539 vertexCount[key1].second = 0;
3541 if( vertexCount.count( key2 )==0 )
3543 vertexCount[key2].first = ri2;
3544 vertexCount[key2].second = 0;
3546 vertexCount[key1].second++;
3547 vertexCount[key2].second--;
3553 fprintf( stderr ,
"Bad Edge 2: %d %d\t%d %d\n" , ri1.key , ri2.key , r1 , r2 );
3557 for(
int i=0 ; i<int(edges.size()) ; i++ )
3559 if( vertexCount.count( edges[i].first.key )==0 ) printf(
"Could not find vertex: %lld\n" , edges[i].first );
3560 else if( vertexCount[ edges[i].first.key ].second )
3563 GetRootPair( vertexCount[edges[i].first.key].first , fData.depth , ri );
3564 long long key = ri.key;
3565 if( vertexCount.count( key )==0 )
3568 node->depthAndOffset( d , off );
3569 printf(
"Vertex pair not in list 1 (%lld) %d\t[%d] (%d %d %d)\n" , key , IsBoundaryEdge( ri.node , ri.edgeIndex , sDepth ) , d , off[0] , off[1] , off[2] );
3573 edges.push_back( std::pair< RootInfo , RootInfo >( ri , edges[i].first ) );
3574 vertexCount[ key ].second++;
3575 vertexCount[ edges[i].first.key ].second--;
3579 if( vertexCount.count( edges[i].second.key )==0 ) printf(
"Could not find vertex: %lld\n" , edges[i].second );
3580 else if( vertexCount[edges[i].second.key].second )
3583 GetRootPair( vertexCount[edges[i].second.key].first , fData.depth , ri );
3584 long long key = ri.key;
3585 if( vertexCount.count( key ) )
3588 node->depthAndOffset( d , off );
3589 printf(
"Vertex pair not in list 2\t[%d] (%d %d %d)\n" , d , off[0] , off[1] , off[2] );
3593 edges.push_back( std::pair< RootInfo , RootInfo >( edges[i].second , ri ) );
3594 vertexCount[key].second--;
3595 vertexCount[ edges[i].second.key ].second++;
3601 template<
int Degree>
3603 int Octree< Degree >::GetMCIsoTriangles(
TreeOctNode* node , CoredMeshData* mesh , RootData& rootData , std::vector<
Point3D< float > >* interiorPositions ,
int offSet ,
int sDepth ,
bool polygonMesh , std::vector<
Point3D< float > >* barycenters )
3605 int Octree< Degree >::GetMCIsoTriangles(
TreeOctNode* node , CoredMeshData* mesh , RootData& rootData , std::vector<
Point3D< float > >* interiorPositions ,
int offSet ,
int sDepth ,
bool addBarycenter ,
bool polygonMesh )
3609 std::vector< std::pair< RootInfo , RootInfo > > edges;
3610 std::vector< std::vector< std::pair< RootInfo , RootInfo > > > edgeLoops;
3611 GetMCIsoEdges( node , sDepth , edges );
3613 GetEdgeLoops( edges , edgeLoops );
3614 for(
int i=0 ; i<int(edgeLoops.size()) ; i++ )
3617 std::vector<CoredPointIndex> edgeIndices;
3618 for(
int j=0 ; j<int(edgeLoops[i].size()) ; j++ )
3620 if( !GetRootIndex( edgeLoops[i][j].first , rootData , p ) ) printf(
"Bad Point Index\n" );
3621 else edgeIndices.push_back( p );
3624 tris += AddTriangles( mesh , edgeIndices , interiorPositions , offSet , polygonMesh , barycenters );
3626 tris += AddTriangles( mesh , edgeIndices , interiorPositions , offSet , addBarycenter , polygonMesh );
3632 template<
int Degree >
3633 int Octree< Degree >::GetEdgeLoops( std::vector< std::pair< RootInfo , RootInfo > >& edges , std::vector< std::vector< std::pair< RootInfo , RootInfo > > >& loops )
3636 long long frontIdx , backIdx;
3637 std::pair< RootInfo , RootInfo > e , temp;
3640 while( edges.size() )
3642 std::vector< std::pair< RootInfo , RootInfo > > front , back;
3644 loops.resize( loopSize+1 );
3645 edges[0] = edges.back();
3647 frontIdx = e.second.key;
3648 backIdx = e.first.key;
3649 for(
int j=
int(edges.size())-1 ; j>=0 ; j-- )
3651 if( edges[j].first.key==frontIdx || edges[j].second.key==frontIdx )
3653 if( edges[j].first.key==frontIdx ) temp = edges[j];
3654 else temp.first = edges[j].second , temp.second = edges[j].first;
3655 frontIdx = temp.second.key;
3656 front.push_back(temp);
3657 edges[j] = edges.back();
3659 j = int(edges.size());
3661 else if( edges[j].first.key==backIdx || edges[j].second.key==backIdx )
3663 if( edges[j].second.key==backIdx ) temp = edges[j];
3664 else temp.first = edges[j].second , temp.second = edges[j].first;
3665 backIdx = temp.first.key;
3666 back.push_back(temp);
3667 edges[j] = edges.back();
3669 j = int(edges.size());
3672 for(
int j=
int(back.size())-1 ; j>=0 ; j-- ) loops[loopSize].push_back( back[j] );
3673 loops[loopSize].push_back(e);
3674 for(
int j=0 ; j<int(front.size()) ; j++ ) loops[loopSize].push_back( front[j] );
3677 return int(loops.size());
3680 template<
int Degree>
3682 int Octree<Degree>::AddTriangles( CoredMeshData* mesh , std::vector<CoredPointIndex>& edges , std::vector<Point3D<float> >* interiorPositions ,
int offSet ,
bool polygonMesh , std::vector<
Point3D< float > >* barycenters )
3684 int Octree<Degree>::AddTriangles( CoredMeshData* mesh , std::vector<CoredPointIndex>& edges , std::vector<Point3D<float> >* interiorPositions ,
int offSet ,
bool addBarycenter ,
bool polygonMesh )
3687 MinimalAreaTriangulation< float > MAT;
3688 std::vector< Point3D< float > > vertices;
3689 std::vector< TriangleIndex > triangles;
3692 std::vector< CoredVertexIndex > vertices( edges.size() );
3693 for(
int i=0 ; i<int(edges.size()) ; i++ )
3695 vertices[i].idx = edges[i].index;
3696 vertices[i].inCore = (edges[i].inCore!=0);
3698#pragma omp critical (add_polygon_access)
3700 mesh->addPolygon( vertices );
3704 if( edges.size()>3 )
3706 bool isCoplanar =
false;
3713 for(
int i=0 ; i<int(edges.size()) ; i++ )
3714 for(
int j=0 ; j<i ; j++ )
3715 if( (i+1)%edges.size()!=j && (j+1)%edges.size()!=i )
3717 Point3D< Real > v1 , v2;
3718 if( edges[i].inCore ) v1 = mesh->inCorePoints[ edges[i].index ];
3719 else v1 = (*interiorPositions)[ edges[i].index-offSet ];
3720 if( edges[j].inCore ) v2 = mesh->inCorePoints[ edges[j].index ];
3721 else v2 = (*interiorPositions)[ edges[j].index-offSet ];
3722 for(
int k=0 ; k<3 ; k++ )
if( v1[k]==v2[k] ) isCoplanar =
true;
3727 c[0] = c[1] = c[2] = 0;
3728 for(
int i=0 ; i<int(edges.size()) ; i++ )
3731 if(edges[i].inCore) p = mesh->inCorePoints[edges[i].index ];
3732 else p = (*interiorPositions)[edges[i].index-offSet];
3735 c /=
Real( edges.size() );
3737#pragma omp critical (add_point_access)
3739 cIdx = mesh->addOutOfCorePoint( c );
3741 barycenters->push_back( c );
3743 interiorPositions->push_back( c );
3746 for(
int i=0 ; i<int(edges.size()) ; i++ )
3748 std::vector< CoredVertexIndex > vertices( 3 );
3749 vertices[0].idx = edges[i ].index;
3750 vertices[1].idx = edges[(i+1)%edges.size()].index;
3751 vertices[2].idx = cIdx;
3752 vertices[0].inCore = (edges[i ].inCore!=0);
3753 vertices[1].inCore = (edges[(i+1)%edges.size()].inCore!=0);
3754 vertices[2].inCore = 0;
3755#pragma omp critical (add_polygon_access)
3757 mesh->addPolygon( vertices );
3760 return int( edges.size() );
3764 vertices.resize( edges.size() );
3766 for(
int i=0 ; i<int(edges.size()) ; i++ )
3769 if( edges[i].inCore ) p = mesh->inCorePoints[edges[i].index ];
3770 else p = (*interiorPositions)[edges[i].index-offSet];
3773 MAT.GetTriangulation( vertices , triangles );
3774 for(
int i=0 ; i<int(triangles.size()) ; i++ )
3776 std::vector< CoredVertexIndex > _vertices( 3 );
3777 for(
int j=0 ; j<3 ; j++ )
3779 _vertices[j].idx = edges[ triangles[i].idx[j] ].index;
3780 _vertices[j].inCore = (edges[ triangles[i].idx[j] ].inCore!=0);
3782#pragma omp critical (add_polygon_access)
3784 mesh->addPolygon( _vertices );
3789 else if( edges.size()==3 )
3791 std::vector< CoredVertexIndex > vertices( 3 );
3792 for(
int i=0 ; i<3 ; i++ )
3794 vertices[i].idx = edges[i].index;
3795 vertices[i].inCore = (edges[i].inCore!=0);
3797#pragma omp critical (add_polygon_access)
3798 mesh->addPolygon( vertices );
3800 return int(edges.size())-2;
3803 template<
int Degree >
3809 fData.setValueTables( fData.VALUE_FLAG );
3811 Real* values =
new float[ res * res * res ];
3812 memset( values , 0 ,
sizeof(
float ) * res * res * res );
3816 if( n->d>depth )
continue;
3817 if( n->d<_minDepth )
continue;
3818 int d , idx[3] , start[3] , end[3];
3819 n->depthAndOffset( d , idx );
3820 for(
int i=0 ; i<3 ; i++ )
3825 fData.setSampleSpan( idx[i] , start[i] , end[i] );
3827 if( !(start[i]&1) ) start[i]++;
3828 if( !( end[i]&1) ) end[i]--;
3831 for(
int x=start[0] ; x<=end[0] ; x+=2 )
3832 for(
int y=start[1] ; y<=end[1] ; y+=2 )
3833 for(
int z=start[2] ; z<=end[2] ; z+=2 )
3835 int xx = (x-1)>>1 ,
yy = (y-1)>>1 ,
zz = (z-1)>>1;
3836 values[
zz*res*res +
yy*res +
xx ] +=
3838 fData.valueTables[ idx[0]+x*fData.functionCount ] *
3839 fData.valueTables[ idx[1]+y*fData.functionCount ] *
3840 fData.valueTables[ idx[2]+z*fData.functionCount ];
3843 for(
int i=0 ; i<res*res*res ; i++ ) values[i] -=
isoValue;
3848 template<
int Degree >
3852 res = 1<<tree.maxDepth();
3853 Real* values =
new float[ res * res * res ];
3854 memset( values , 0 ,
sizeof(
float ) * res * res * res );
3858 if( n->d>depth )
continue;
3859 int d , idx[3] , start[3] , end[3];
3860 n->depthAndOffset( d , idx );
3861 for(
int i=0 ; i<3 ; i++ )
3866 fData.setSampleSpan( idx[i] , start[i] , end[i] );
3868 if( !(start[i]&1) ) start[i]++;
3869 if( !( end[i]&1) ) end[i]--;
3871 for(
int x=start[0] ; x<=end[0] ; x+=2 )
3872 for(
int y=start[1] ; y<=end[1] ; y+=2 )
3873 for(
int z=start[2] ; z<=end[2] ; z+=2 )
3875 int xx = (x-1)>>1 ,
yy = (y-1)>>1 ,
zz = (z-1)>>1;
3876 values[
zz*res*res +
yy*res +
xx ] +=
3877 n->nodeData.centerWeightContribution *
3878 fData.valueTables[ idx[0]+x*fData.functionCount ] *
3879 fData.valueTables[ idx[1]+y*fData.functionCount ] *
3880 fData.valueTables[ idx[2]+z*fData.functionCount ];
3891 return CenterIndex(node,
maxDepth,idx);
3898 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3903 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3918 return CornerIndexKey( idx );
3926 return CornerIndexKey( idx );
3931 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3936 return FaceIndex(node,fIndex,
maxDepth,idx);
3946 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3974 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
Iterator class for point clouds with or without given indices.
std::size_t size() const
Size of the range the iterator is going through.
shared_ptr< const PointCloud< PointT > > ConstPtr
static int SymmetricIndex(int i1, int i2)
static int Index(int depth, int offSet)
static int CornerIndex(int depth, int offSet)
static int CenterIndex(int depth, int offSet)
static int FaceReflectEdgeIndex(int idx, int faceIndex)
static int FaceAdjacentToEdges(int eIndex1, int eIndex2)
static void FactorEdgeIndex(int idx, int &orientation, int &i, int &j)
static void FactorCornerIndex(int idx, int &x, int &y, int &z)
static void FacesAdjacentToEdge(int eIndex, int &f1Index, int &f2Index)
static int FaceReflectFaceIndex(int idx, int faceIndex)
static void FaceCorners(int idx, int &c1, int &c2, int &c3, int &c4)
static int AntipodalCornerIndex(int idx)
static int CornerIndex(int x, int y, int z)
static int EdgeIndex(int orientation, int i, int j)
static int EdgeReflectEdgeIndex(int edgeIndex)
static void FactorFaceIndex(int idx, int &x, int &y, int &z)
static void EdgeCorners(int idx, int &c1, int &c2)
static const int * cornerMap()
static int AddTriangleIndices(int mcIndex, int *triangles)
static int HasEdgeRoots(int mcIndex, int edgeIndex)
static int HasRoots(const double v[Cube::CORNERS], double isoValue)
static int GetIndex(const double values[Cube::CORNERS], double iso)
static const int * edgeMask()
ConstNeighbors3 * neighbors
const OctNode * neighbors[3][3][3]
void depthAndOffset(int &depth, int offset[DIMENSION]) const
static int CornerIndex(const Point3D< Real > ¢er, const Point3D< Real > &p)
const OctNode * nextNode(const OctNode *currentNode=NULL) const
void centerAndWidth(Point3D< Real > ¢er, Real &width) const
const OctNode * nextLeaf(const OctNode *currentLeaf=NULL) const
static void ProcessFixedDepthNodeAdjacentNodes(int maxDepth, OctNode *node1, int width1, OctNode *node2, int width2, int depth, NodeAdjacencyFunction *F, int processCurrent=1)
static void Index(int depth, const int offset[3], short &d, short off[DIMENSION])
OctNode * faceNeighbor(int faceIndex, int forceChildren=0)
void RefineBoundary(int subdivisionDepth)
void setBSplineData(int maxDepth, Real normalSmooth=-1, bool reflectBoundary=false)
static double MemoryUsage(void)
int LaplacianMatrixIteration(int subdivideDepth, bool showResidual, int minIters, double accuracy)
Real * GetSolutionGrid(int &res, float isoValue=0.f, int depth=-1)
void SetLaplacianConstraints(void)
Real * GetWeightGrid(int &res, int depth=-1)
int setTree(typename pcl::PointCloud< PointNT >::ConstPtr input_, int maxDepth, int minDepth, int kernelDepth, Real samplesPerNode, Real scaleFactor, Point3D< Real > ¢er, Real &scale, int useConfidence, Real constraintWeight, bool adaptiveWeights)
An exception that is thrown when something goes wrong inside an openMP for loop.
void setEdgeTable(EdgeTableData &eData, const TreeOctNode *rootNode, int depth, int threads)
int getMaxEdgeCount(const TreeOctNode *rootNode, int depth, int threads) const
int getMaxCornerCount(const TreeOctNode *rootNode, int depth, int maxDepth, int threads) const
void set(TreeOctNode &root)
void setCornerTable(CornerTableData &cData, const TreeOctNode *rootNode, int depth, int threads) const
static void SetAllocator(int blockSize)
static Allocator< MatrixEntry< T > > internalAllocator
static int Solve(const SparseSymmetricMatrix< T > &M, const Vector< T2 > &b, int iters, Vector< T2 > &solution, T2 eps=1e-8, int reset=1, int threads=0, bool addDCTerm=false, bool solveNormal=false)
static int CornerIndex(int x, int y)
static int AntipodalCornerIndex(int idx)
static void FactorCornerIndex(int idx, int &x, int &y)
static long long CenterIndex(int depth, const int offSet[DIMENSION], int maxDepth, int index[DIMENSION])
static long long CornerIndex(int depth, const int offSet[DIMENSION], int cIndex, int maxDepth, int index[DIMENSION])
static long long FaceIndex(const TreeOctNode *node, int fIndex, int maxDepth, int index[DIMENSION])
static long long CornerIndexKey(const int index[DIMENSION])
static long long EdgeIndex(const TreeOctNode *node, int eIndex, int maxDepth, int index[DIMENSION])
__device__ __forceinline__ float dot(const float3 &v1, const float3 &v2)
pcl::poisson::OctNode< class TreeNodeData, Real > TreeOctNode
void atomicOr(volatile int &dest, int value)
const Real MATRIX_ENTRY_EPSILON
double Length(const Point3D< Real > &p)
std::vector< int > offsets
CornerIndices & operator[](const TreeOctNode *node)
std::vector< CornerIndices > cTable
CornerIndices & cornerIndices(const TreeOctNode *node)
EdgeIndices & operator[](const TreeOctNode *node)
EdgeIndices & edgeIndices(const TreeOctNode *node)
UpSampleData(int s, double v1, double v2)