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__
59 _InterlockedOr( (
long volatile*)&dest, value );
61 InterlockedOr( (
long volatile*)&dest , value );
63 #else // !_WIN32 || __MINGW32__
66 #endif // _WIN32 && !__MINGW32__
113 if( threads<=0 ) threads = 1;
115 std::vector< std::pair< int , int > > spans( this->maxDepth , std::pair< int , int >( -1 , -1 ) );
116 int minDepth , off[3];
118 cData.
offsets.resize( this->maxDepth , -1 );
122 for(
int d=minDepth ; d<=
maxDepth ; d++ )
124 spans[d] = std::pair< int , int >( start , end+1 );
126 nodeCount += spans[d].second - spans[d].first;
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;
146 for(
int d=minDepth ; d<=
maxDepth ; d++ )
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++ )
152 if( d<maxDepth && node->children )
continue;
153 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node , minDepth );
156 bool cornerOwner =
true;
164 xx += x , yy += y , zz += z;
165 if( neighbors.neighbors[xx][yy][zz] )
166 if( cc<ac || ( d<
maxDepth && neighbors.neighbors[xx][yy][zz]->children ) )
169 neighbors.neighbors[xx][yy][zz]->depthAndOffset( _d , _off );
170 _off[0] >>= (d-minDepth) , _off[1] >>= (d-minDepth) , _off[2] >>= (d-minDepth);
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];
216 unsigned int paralellExceptionCount = 0;
217 #pragma omp parallel for num_threads( threads )
218 for (
int t = 0; t < threads; t++)
220 for (
int d = minDepth; d <=
maxDepth; d++)
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++)
233 ++paralellExceptionCount;
240 idx = rem + offsets[div];
247 if (paralellExceptionCount > 0)
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++ )
260 std::vector< int >& _cornerCount = cornerCount[t];
261 TreeOctNode::ConstNeighborKey3 neighborKey;
264 for(
int i=(range*t)/threads ; i<(range*(t+1))/threads ; i++ )
269 if( d<maxDepth && node->children )
continue;
271 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node , depth );
274 bool cornerOwner =
true;
282 xx += x , yy += y , zz += z;
283 if( neighbors.neighbors[xx][yy][zz] )
284 if( cc<ac || ( d<
maxDepth && neighbors.neighbors[xx][yy][zz]->children ) )
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++ )
298 for(
int t=0 ; t<threads ; t++ ) c += cornerCount[t][i];
299 maxCount = std::max< int >( maxCount , c );
309 if( threads<=0 ) threads = 1;
310 std::vector< std::pair< int , int > > spans( this->maxDepth , std::pair< int , int >( -1 , -1 ) );
312 int minDepth , off[3];
314 eData.
offsets.resize( this->maxDepth , -1 );
318 for(
int d=minDepth ; d<=
maxDepth ; d++ )
320 spans[d] = std::pair< int , int >( start , end+1 );
322 nodeCount += spans[d].second - spans[d].first;
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;
342 for(
int d=minDepth ; d<=
maxDepth ; d++ )
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 );
352 bool edgeOwner =
true;
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;
367 if( neighbors.neighbors[x][y][z] && cc<ac ) { edgeOwner = false ;
break; }
374 int ii , jj , aii , ajj , x , y , z;
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] )
385 eData[ neighbors.neighbors[x][y][z] ][
Cube::EdgeIndex( o , aii , ajj ) ] = count[t]+offset;
394 std::vector< int > offsets( threads+1 );
396 for (
int t = 0; t < threads; t++)
399 offsets[t + 1] = offsets[t] + count[t];
402 unsigned int paralellExceptionCount = 0;
403 #pragma omp parallel for num_threads( threads )
404 for (
int t = 0; t < threads; t++)
406 for (
int d = minDepth; d <=
maxDepth; d++)
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++)
419 ++paralellExceptionCount;
426 idx = rem + offsets[div];
433 if(paralellExceptionCount > 0)
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++ )
447 std::vector< int >& _edgeCount = edgeCount[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 );
460 bool edgeOwner =
true;
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;
475 if( neighbors.neighbors[x][y][z] && cc<ac ) { edgeOwner = false ;
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++ )
485 for(
int t=0 ; t<threads ; t++ ) c += edgeCount[t][i];
486 maxCount = std::max< int >( maxCount , c );
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
576 #endif // SPLAT_ORDER
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];
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
765 int useConfidence ,
Real constraintWeight ,
bool adaptiveWeights )
767 _minDepth = std::min< int >( std::max< int >( 0 , minDepth ) ,
maxDepth );
768 _constrainValues = (constraintWeight>0);
769 double pointWeightSum = 0;
778 splatDepth = kernelDepth;
779 if( splatDepth<0 ) splatDepth = 0;
782 tree.setFullDepth( _minDepth );
784 while (cnt != input_->
size ())
787 p[0] = input_->
points[cnt].x;
788 p[1] = input_->
points[cnt].y;
789 p[2] = input_->
points[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] ) );
800 center = ( max+min ) /2;
802 scale *= scaleFactor;
803 for( i=0 ; i<DIMENSION ; i++ ) center[i] -= scale/2;
807 while (cnt != input_->
size ())
809 position[0] = input_->
points[cnt].x;
810 position[1] = input_->
points[cnt].y;
811 position[2] = input_->
points[cnt].z;
812 normal[0] = input_->
points[cnt].normal_x;
813 normal[1] = input_->
points[cnt].normal_y;
814 normal[2] = input_->
points[cnt].normal_z;
816 for( i=0 ; i<DIMENSION ; i++ ) position[i] = ( position[i]-center[i] ) / scale;
817 myCenter[0] = myCenter[1] = myCenter[2] =
Real(0.5);
819 for( i=0 ; i<DIMENSION ; i++ )
if( position[i]<myCenter[i]-myWidth/2 || position[i]>myCenter[i]+myWidth/2 )
break;
820 if( i!=DIMENSION )
continue;
822 if( useConfidence ) weight =
Real(
Length(normal) );
825 while( d<splatDepth )
827 NonLinearUpdateWeightContribution( temp , position , weight );
832 if(cIndex&1) myCenter[0] += myWidth/2;
833 else myCenter[0] -= myWidth/2;
834 if(cIndex&2) myCenter[1] += myWidth/2;
835 else myCenter[1] -= myWidth/2;
836 if(cIndex&4) myCenter[2] += myWidth/2;
837 else myCenter[2] -= myWidth/2;
840 NonLinearUpdateWeightContribution( temp , position , weight );
845 normals =
new std::vector< Point3D<Real> >();
847 while (cnt != input_->
size ())
849 position[0] = input_->
points[cnt].x;
850 position[1] = input_->
points[cnt].y;
851 position[2] = input_->
points[cnt].z;
852 normal[0] = input_->
points[cnt].normal_x;
853 normal[1] = input_->
points[cnt].normal_y;
854 normal[2] = input_->
points[cnt].normal_z;
856 for( i=0 ; i<DIMENSION ; i++ ) position[i] = ( position[i]-center[i] ) / scale;
857 myCenter[0] = myCenter[1] = myCenter[2] =
Real(0.5);
859 for( i=0 ; i<DIMENSION ; i++ )
if(position[i]<myCenter[i]-myWidth/2 || position[i]>myCenter[i]+myWidth/2)
break;
860 if( i!=DIMENSION )
continue;
862 if( l!=l || l<=
EPSILON )
continue;
863 if( !useConfidence ) normal /= l;
867 if( samplesPerNode>0 && splatDepth )
869 pointWeight = NonLinearSplatOrientedPoint( position , normal , splatDepth , samplesPerNode , _minDepth ,
maxDepth );
878 while( d<splatDepth )
883 if(cIndex&1) myCenter[0]+=myWidth/2;
884 else myCenter[0]-=myWidth/2;
885 if(cIndex&2) myCenter[1]+=myWidth/2;
886 else myCenter[1]-=myWidth/2;
887 if(cIndex&4) myCenter[2]+=myWidth/2;
888 else myCenter[2]-=myWidth/2;
891 alpha = NonLinearGetSampleWeight( temp , position );
893 for( i=0 ; i<DIMENSION ; i++ ) normal[i]*=alpha;
900 if(cIndex&1) myCenter[0]+=myWidth/2;
901 else myCenter[0]-=myWidth/2;
902 if(cIndex&2) myCenter[1]+=myWidth/2;
903 else myCenter[1]-=myWidth/2;
904 if(cIndex&4) myCenter[2]+=myWidth/2;
905 else myCenter[2]-=myWidth/2;
908 NonLinearSplatOrientedPoint( temp , position , normal );
912 pointWeightSum += pointWeight;
913 if( _constrainValues )
917 myCenter[0] = myCenter[1] = myCenter[2] =
Real(0.5);
925 p[0] = p[1] = p[2] = 0;
926 idx = int( _points.size() );
927 _points.push_back( PointData( position*pointWeight , pointWeight ) );
932 _points[idx].weight += pointWeight;
933 _points[idx].position += position * pointWeight;
940 if( cIndex&1 ) myCenter[0] += myWidth/2;
941 else myCenter[0] -= myWidth/2;
942 if( cIndex&2 ) myCenter[1] += myWidth/2;
943 else myCenter[1] -= myWidth/2;
944 if( cIndex&4 ) myCenter[2] += myWidth/2;
945 else myCenter[2] -= myWidth/2;
952 if( _constrainValues )
954 if( n->nodeData.pointIndex!=-1 )
956 int idx = n->nodeData.pointIndex;
957 _points[idx].position /= _points[idx].weight;
958 if( adaptiveWeights ) _points[idx].weight *= (1<<n->d);
959 else _points[idx].weight *= (1<<
maxDepth);
960 _points[idx].weight *=
Real( constraintWeight / pointWeightSum );
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;
972 #endif // FORCE_NEUMANN_FIELD
983 radius = 0.5 + 0.5 * Degree;
984 width=int(
double(radius+0.5-
EPSILON)*2);
985 if( normalSmooth>0 ) postNormalSmooth = normalSmooth;
986 fData.set(
maxDepth ,
true , reflectBoundary );
993 TreeOctNode::NeighborKey5 nKey;
999 int xStart=0 , xEnd=5 , yStart=0 , yEnd=5 , zStart=0 , zEnd=5;
1000 int c = int( node - node->parent->children );
1009 nKey.setNeighbors( node->parent , xStart , xEnd , yStart , yEnd , zStart , zEnd );
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] )
1064 #else // !GRADIENT_DOMAIN_SOLUTION
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] )
1069 #endif // GRADIENT_DOMAIN_SOLUTION
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] ) );
1074 #else // !GRADIENT_DOMAIN_SOLUTION
1075 return -
Real( dot * ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] ) );
1076 #endif // GRADIENT_DOMAIN_SOLUTION
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] )
1094 #else // !GRADIENT_DOMAIN_SOLUTION
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] )
1099 #endif // GRADIENT_DOMAIN_SOLUTION
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 ); }
1131 template<
int Degree >
1132 int Octree< Degree >::GetMatrixRowSize(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 ,
int xStart ,
int xEnd ,
int yStart ,
int yEnd ,
int zStart ,
int zEnd )
const
1135 for(
int x=xStart ; x<=2 ; x++ )
1136 for(
int y=yStart ; y<yEnd ; y++ )
1137 if( x==2 && y>2 )
continue;
1138 else for(
int z=zStart ; z<zEnd ; z++ )
1139 if( x==2 && y==2 && z>2 )
continue;
1140 else if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1144 template<
int Degree >
1145 int Octree< Degree >::SetMatrixRow(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 , MatrixEntry< float >* row ,
int offset ,
const double stencil[5][5][5] )
const
1147 return SetMatrixRow( neighbors5 , row , offset , stencil , 0 , 5 , 0 , 5 , 0 , 5 );
1149 template<
int Degree >
1150 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
1152 bool hasPoints[3][3];
1154 PointInfo samples[3][3][3];
1157 const TreeOctNode* node = neighbors5.neighbors[2][2][2];
1158 int index[] = { int( node->off[0] ) , int( node->off[1] ), int( node->off[2] ) };
1160 if( _constrainValues )
1167 for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
1169 hasPoints[j][k] =
false;
1170 for(
int l=0 ; l<3 ; l++ )
1172 const TreeOctNode* _node = neighbors5.neighbors[j+1][k+1][l+1];
1173 if( _node && _node->nodeData.pointIndex!=-1 )
1175 const PointData& pData = _points[ _node->nodeData.pointIndex ];
1176 PointInfo& pointInfo = samples[j][k][l];
1177 Real weight = pData.weight;
1178 Point3D< Real > p = pData.position;
1179 for(
int s=0 ; s<3 ; s++ )
1181 pointInfo.splineValues[0][s] = float( fData.baseBSplines[ idx[0]+j-s][s]( p[0] ) );
1182 pointInfo.splineValues[1][s] = float( fData.baseBSplines[ idx[1]+k-s][s]( p[1] ) );
1183 pointInfo.splineValues[2][s] = float( fData.baseBSplines[ idx[2]+l-s][s]( p[2] ) );
1185 float value = pointInfo.splineValues[0][j] * pointInfo.splineValues[1][k] * pointInfo.splineValues[2][l];
1186 diagonal += value * value * weight;
1187 pointInfo.weightedValue = value * weight;
1188 for(
int s=0 ; s<3 ; s++ ) pointInfo.splineValues[0][s] *= pointInfo.weightedValue;
1189 hasPoints[j][k] =
true;
1191 else samples[j][k][l].weightedValue = 0;
1198 neighbors5.neighbors[2][2][2]->depthAndOffset( d , off );
1199 int mn = 2 , mx = (1<<d)-2;
1200 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
1201 for(
int x=xStart ; x<=2 ; x++ )
1202 for(
int y=yStart ; y<yEnd ; y++ )
1203 if( x==2 && y>2 )
continue;
1204 else for(
int z=zStart ; z<zEnd ; z++ )
1205 if( x==2 && y==2 && z>2 )
continue;
1206 else if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1208 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
1209 int _index[] = { int( _node->off[0] ) , int( _node->off[1] ), int( _node->off[2] ) };
1211 if( isInterior ) temp =
Real( stencil[x][y][z] );
1212 else temp = GetLaplacian( node , _node );
1213 if( _constrainValues )
1215 int _d[] = { _index[0]-index[0] , _index[1]-index[1] , _index[2]-index[2] };
1216 if( x==2 && y==2 && z==2 ) temp += diagonal;
1217 else temp += GetValue( samples , hasPoints , _d );
1219 if( x==2 && y==2 && z==2 ) temp /= 2;
1222 row[count].N = _node->nodeData.nodeIndex-offset;
1223 row[count].Value = temp;
1229 template<
int Degree >
1230 void Octree< Degree >::SetDivergenceStencil(
int depth , Point3D< double > *stencil ,
bool scatter )
const
1232 int offset[] = { 2 , 2 , 2 };
1235 int index1[3] , index2[3];
1236 if( scatter ) index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1237 else index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1238 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1240 int _offset[] = { x , y , z };
1242 if( scatter ) index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1243 else index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1252 #if GRADIENT_DOMAIN_SOLUTION
1254 fData.Index( index1[0] , index2[0] ) ,
1255 fData.Index( index1[1] , index2[1] ) ,
1256 fData.Index( index1[2] , index2[2] )
1257 #else // !GRADIENT_DOMAIN_SOLUTION
1259 fData.Index( index2[0] , index1[0] ) ,
1260 fData.Index( index2[1] , index1[1] ) ,
1261 fData.Index( index2[2] , index1[2] )
1262 #endif // GRADIENT_DOMAIN_SOLUTION
1265 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1266 #if GRADIENT_DOMAIN_SOLUTION
1267 Point3D<double> temp;
1268 temp[0] = fData.dvDotTable[aSymIndex[0]] *
dot;
1269 temp[1] = fData.dvDotTable[aSymIndex[1]] *
dot;
1270 temp[2] = fData.dvDotTable[aSymIndex[2]] *
dot;
1271 stencil[25*x + 5*y + z] = temp;
1275 #else // !GRADIENT_DOMAIN_SOLUTION
1276 Point3D<double> temp;
1277 temp[0] = -fData.dvDotTable[aSymIndex[0]] *
dot;
1278 temp[1] = -fData.dvDotTable[aSymIndex[1]] *
dot;
1279 temp[2] = -fData.dvDotTable[aSymIndex[2]] *
dot;
1280 stencil[25*x + 5*y + z] = temp;
1284 #endif // GRADIENT_DOMAIN_SOLUTION
1287 template<
int Degree >
1288 void Octree< Degree >::SetLaplacianStencil(
int depth ,
double stencil[5][5][5] )
const
1290 int offset[] = { 2 , 2 , 2 };
1293 int index[] = { int( off[0] ) , int( off[1] ) , int( off[2] ) };
1294 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1296 int _offset[] = { x , y , z };
1299 int _index[] = { int( _off[0] ) , int( _off[1] ) , int( _off[2] ) };
1304 stencil[x][y][z] = GetLaplacian( symIndex );
1307 template<
int Degree >
1308 void Octree< Degree >::SetLaplacianStencils(
int depth , Stencil< double , 5 > stencils[2][2][2] )
const
1310 if( depth<=1 )
return;
1311 for(
int i=0 ; i<2 ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
1314 int offset[] = { 4+i , 4+j , 4+k };
1316 int index[] = { int( off[0] ) , int( off[1] ) , int( off[2] ) };
1317 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1319 int _offset[] = { x , y , z };
1322 int _index[] = { int( _off[0] ) , int( _off[1] ) , int( _off[2] ) };
1327 stencils[i][j][k].values[x][y][z] = GetLaplacian( symIndex );
1331 template<
int Degree >
1332 void Octree< Degree >::SetDivergenceStencils(
int depth , Stencil< Point3D< double > , 5 > stencils[2][2][2] ,
bool scatter )
const
1334 if( depth<=1 )
return;
1335 int index1[3] , index2[3];
1336 for(
int i=0 ; i<2 ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
1339 int offset[] = { 4+i , 4+j , 4+k };
1341 if( scatter ) index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1342 else index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1343 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1345 int _offset[] = { x , y , z };
1347 if( scatter ) index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1348 else index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1358 #if GRADIENT_DOMAIN_SOLUTION
1360 fData.Index( index1[0] , index2[0] ) ,
1361 fData.Index( index1[1] , index2[1] ) ,
1362 fData.Index( index1[2] , index2[2] )
1363 #else // !GRADIENT_DOMAIN_SOLUTION
1365 fData.Index( index2[0] , index1[0] ) ,
1366 fData.Index( index2[1] , index1[1] ) ,
1367 fData.Index( index2[2] , index1[2] )
1368 #endif // GRADIENT_DOMAIN_SOLUTION
1370 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1371 #if GRADIENT_DOMAIN_SOLUTION
1372 stencils[i][j][k].values[x][y][z][0] = fData.dvDotTable[aSymIndex[0]] *
dot;
1373 stencils[i][j][k].values[x][y][z][1] = fData.dvDotTable[aSymIndex[1]] *
dot;
1374 stencils[i][j][k].values[x][y][z][2] = fData.dvDotTable[aSymIndex[2]] *
dot;
1375 #else // !GRADIENT_DOMAIN_SOLUTION
1376 stencils[i][j][k].values[x][y][z][0] = -fData.dvDotTable[aSymIndex[0]] *
dot;
1377 stencils[i][j][k].values[x][y][z][1] = -fData.dvDotTable[aSymIndex[1]] *
dot;
1378 stencils[i][j][k].values[x][y][z][2] = -fData.dvDotTable[aSymIndex[2]] *
dot;
1379 #endif // GRADIENT_DOMAIN_SOLUTION
1383 template<
int Degree >
1384 void Octree< Degree >::SetEvaluationStencils(
int depth , Stencil< double , 3 > stencil1[8] , Stencil< double , 3 > stencil2[8][8] )
const
1389 int off[] = { 2 , 2 , 2 };
1390 for(
int c=0 ; c<8 ; c++ )
1393 idx[0] *= fData.functionCount , idx[1] *= fData.functionCount , idx[2] *= fData.functionCount;
1394 for(
int x=0 ; x<3 ; x++ )
for(
int y=0 ; y<3 ; y++ )
for(
int z=0 ; z<3 ; z++ )
1397 int _offset[] = { x+1 , y+1 , z+1 };
1399 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]) ];
1404 for(
int _c=0 ; _c<8 ; _c++ )
1407 int _cx , _cy , _cz;
1409 int off[] = { 4+_cx , 4+_cy , 4+_cz };
1410 for(
int c=0 ; c<8 ; c++ )
1413 idx[0] *= fData.functionCount , idx[1] *= fData.functionCount , idx[2] *= fData.functionCount;
1414 for(
int x=0 ; x<3 ; x++ )
for(
int y=0 ; y<3 ; y++ )
for(
int z=0 ; z<3 ; z++ )
1417 int _offset[] = { x+1 , y+1 , z+1 };
1419 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]) ];
1424 template<
int Degree >
1425 void Octree< Degree >::UpdateCoarserSupportBounds(
const TreeOctNode* node ,
int& startX ,
int& endX ,
int& startY ,
int& endY ,
int& startZ ,
int& endZ )
1429 int x , y , z , c = int( node - node->parent->children );
1431 if( x==0 ) endX = 4;
1433 if( y==0 ) endY = 4;
1435 if( z==0 ) endZ = 4;
1440 template<
int Degree >
1441 void Octree< Degree >::UpdateConstraintsFromCoarser(
const OctNode< TreeNodeData , Real >::NeighborKey5& neighborKey5 ,
TreeOctNode* node ,
Real* metSolution ,
const Stencil< double , 5 >& lapStencil )
const
1446 node->depthAndOffset( d , off );
1447 int mn = 4 , mx = (1<<d)-4;
1448 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
1451 int depth = node->depth();
1452 if( depth<=_minDepth )
return;
1453 int i = node->nodeData.nodeIndex;
1455 int startX = 0 , endX = 5 , startY = 0 , endY = 5 , startZ = 0 , endZ = 5;
1456 UpdateCoarserSupportBounds( node , startX , endX , startY , endY , startZ , endZ );
1458 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth-1];
1459 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
1460 if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1462 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
1463 Real _solution = metSolution[ _node->nodeData.nodeIndex ];
1465 if( isInterior ) node->nodeData.constraint -=
Real( lapStencil.values[x][y][z] * _solution );
1466 else node->nodeData.constraint -= GetLaplacian( _node , node ) * _solution;
1469 if( _constrainValues )
1472 node->depthAndOffset( d, idx );
1476 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth];
1477 for(
int x=1 ; x<4 ; x++ )
for(
int y=1 ; y<4 ; y++ )
for(
int z=1 ; z<4 ; z++ )
1478 if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.pointIndex!=-1 )
1480 const PointData& pData = _points[ neighbors5.neighbors[x][y][z]->nodeData.pointIndex ];
1481 Real pointValue = pData.value;
1482 Point3D< Real > p = pData.position;
1483 node->nodeData.constraint -=
1485 fData.baseBSplines[idx[0]][x-1]( p[0] ) *
1486 fData.baseBSplines[idx[1]][y-1]( p[1] ) *
1487 fData.baseBSplines[idx[2]][z-1]( p[2] ) *
1498 UpSampleData(
int s ,
double v1 ,
double v2 ) { start = s , v[0] = v1 , v[1] = v2; }
1500 template<
int Degree >
1501 void Octree< Degree >::UpSampleCoarserSolution(
int depth ,
const SortedTreeNodes& sNodes , Vector< Real >& Solution )
const
1503 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1504 Solution.Resize( range );
1505 if( !depth )
return;
1506 else if( depth==1 )
for(
int i=start ; i<end ; i++ ) Solution[i-start] += sNodes.treeNodes[0]->nodeData.solution;
1510 #pragma omp parallel
for num_threads( threads )
1511 for(
int t=0 ; t<threads ; t++ )
1513 TreeOctNode::NeighborKey3 neighborKey;
1514 neighborKey.set( depth );
1515 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1518 UpSampleData usData[3];
1519 sNodes.treeNodes[i]->depthAndOffset( d , off );
1520 for(
int d=0 ; d<3 ; d++ )
1522 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1523 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1524 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1525 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1527 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1528 for(
int ii=0 ; ii<2 ; ii++ )
1530 int _ii = ii + usData[0].start;
1531 double dx = usData[0].v[ii];
1532 for(
int jj=0 ; jj<2 ; jj++ )
1534 int _jj = jj + usData[1].start;
1535 double dxy = dx * usData[1].v[jj];
1536 for(
int kk=0 ; kk<2 ; kk++ )
1538 int _kk = kk + usData[2].start;
1539 double dxyz = dxy * usData[2].v[kk];
1540 Solution[i-start] +=
Real( neighborKey.neighbors[depth-1].neighbors[_ii][_jj][_kk]->nodeData.solution * dxyz );
1548 start = sNodes.nodeCount[depth-1] , end = sNodes.nodeCount[depth] , range = end-start;
1549 #pragma omp parallel for num_threads( threads )
1550 for(
int i=start ; i<end ; i++ ) sNodes.treeNodes[i]->nodeData.solution =
Real( 0. );
1552 template<
int Degree >
1553 void Octree< Degree >::DownSampleFinerConstraints(
int depth , SortedTreeNodes& sNodes )
const
1555 if( !depth )
return;
1556 #pragma omp parallel for num_threads( threads )
1557 for(
int i=sNodes.nodeCount[depth-1] ; i<sNodes.nodeCount[depth] ; i++ )
1558 sNodes.treeNodes[i]->nodeData.constraint =
Real( 0 );
1562 sNodes.treeNodes[0]->nodeData.constraint =
Real( 0 );
1563 for(
int i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ ) sNodes.treeNodes[0]->nodeData.constraint += sNodes.treeNodes[i]->nodeData.constraint;
1566 std::vector< Vector< double > > constraints( threads );
1567 for(
int t=0 ; t<threads ; t++ ) constraints[t].Resize( sNodes.nodeCount[depth] - sNodes.nodeCount[depth-1] ) , constraints[t].SetZero();
1568 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1569 int lStart = sNodes.nodeCount[depth-1] , lEnd = sNodes.nodeCount[depth];
1571 #pragma omp parallel for num_threads( threads )
1572 for(
int t=0 ; t<threads ; t++ )
1574 TreeOctNode::NeighborKey3 neighborKey;
1575 neighborKey.set( depth );
1576 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1579 UpSampleData usData[3];
1580 sNodes.treeNodes[i]->depthAndOffset( d , off );
1581 for(
int d=0 ; d<3 ; d++ )
1583 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1584 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1585 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1586 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1588 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1589 TreeOctNode::Neighbors3& neighbors = neighborKey.neighbors[depth-1];
1590 for(
int ii=0 ; ii<2 ; ii++ )
1592 int _ii = ii + usData[0].start;
1593 double dx = usData[0].v[ii];
1594 for(
int jj=0 ; jj<2 ; jj++ )
1596 int _jj = jj + usData[1].start;
1597 double dxy = dx * usData[1].v[jj];
1598 for(
int kk=0 ; kk<2 ; kk++ )
1600 int _kk = kk + usData[2].start;
1601 double dxyz = dxy * usData[2].v[kk];
1602 constraints[t][neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex-lStart] += sNodes.treeNodes[i]->nodeData.constraint * dxyz;
1608 #pragma omp parallel for num_threads( threads )
1609 for(
int i=lStart ; i<lEnd ; i++ )
1612 for(
int t=0 ; t<threads ; t++ ) cSum += constraints[t][i-lStart];
1613 sNodes.treeNodes[i]->nodeData.constraint += cSum;
1616 template<
int Degree >
1618 void Octree< Degree >::DownSample(
int depth ,
const SortedTreeNodes& sNodes , C* constraints )
const
1620 if( depth==0 )
return;
1623 for(
int i=sNodes.nodeCount[1] ; i<sNodes.nodeCount[2] ; i++ ) constraints[0] += constraints[i];
1626 std::vector< Vector< C > > _constraints( threads );
1627 for(
int t=0 ; t<threads ; t++ ) _constraints[t].Resize( sNodes.nodeCount[depth] - sNodes.nodeCount[depth-1] );
1628 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start , lStart = sNodes.nodeCount[depth-1] , lEnd = sNodes.nodeCount[depth];
1630 #pragma omp parallel for num_threads( threads )
1631 for(
int t=0 ; t<threads ; t++ )
1633 TreeOctNode::NeighborKey3 neighborKey;
1634 neighborKey.set( depth );
1635 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1638 UpSampleData usData[3];
1639 sNodes.treeNodes[i]->depthAndOffset( d , off );
1640 for(
int d=0 ; d<3 ; d++ )
1642 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1643 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1644 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1645 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1647 TreeOctNode::Neighbors3& neighbors = neighborKey.getNeighbors( sNodes.treeNodes[i]->parent );
1648 C c = constraints[i];
1649 for(
int ii=0 ; ii<2 ; ii++ )
1651 int _ii = ii + usData[0].start;
1652 C cx = C( c*usData[0].v[ii] );
1653 for(
int jj=0 ; jj<2 ; jj++ )
1655 int _jj = jj + usData[1].start;
1656 C cxy = C( cx*usData[1].v[jj] );
1657 for(
int kk=0 ; kk<2 ; kk++ )
1659 int _kk = kk + usData[2].start;
1660 if( neighbors.neighbors[_ii][_jj][_kk] )
1661 _constraints[t][neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex-lStart] += C( cxy*usData[2].v[kk] );
1667 #pragma omp parallel for num_threads( threads )
1668 for(
int i=lStart ; i<lEnd ; i++ )
1671 for(
int t=0 ; t<threads ; t++ ) cSum += _constraints[t][i-lStart];
1672 constraints[i] += cSum;
1675 template<
int Degree >
1677 void Octree< Degree >::UpSample(
int depth ,
const SortedTreeNodes& sNodes , C* coefficients )
const
1679 if ( depth==0 )
return;
1682 for(
int i=sNodes.nodeCount[1] ; i<sNodes.nodeCount[2] ; i++ ) coefficients[i] += coefficients[0];
1686 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1688 #pragma omp parallel for num_threads( threads )
1689 for(
int t=0 ; t<threads ; t++ )
1691 TreeOctNode::NeighborKey3 neighborKey;
1692 neighborKey.set( depth-1 );
1693 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1697 UpSampleData usData[3];
1699 for(
int d=0 ; d<3 ; d++ )
1701 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1702 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1703 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1704 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1706 TreeOctNode::Neighbors3& neighbors = neighborKey.getNeighbors( node->parent );
1707 for(
int ii=0 ; ii<2 ; ii++ )
1709 int _ii = ii + usData[0].start;
1710 double dx = usData[0].v[ii];
1711 for(
int jj=0 ; jj<2 ; jj++ )
1713 int _jj = jj + usData[1].start;
1714 double dxy = dx * usData[1].v[jj];
1715 for(
int kk=0 ; kk<2 ; kk++ )
1717 int _kk = kk + usData[2].start;
1718 if( neighbors.neighbors[_ii][_jj][_kk] )
1720 double dxyz = dxy * usData[2].v[kk];
1721 int _i = neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex;
1722 coefficients[i] += coefficients[_i] *
Real( dxyz );
1730 template<
int Degree >
1731 void Octree< Degree >::SetCoarserPointValues(
int depth ,
const SortedTreeNodes& sNodes ,
Real* metSolution )
1733 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1735 #pragma omp parallel for num_threads( threads )
1736 for(
int t=0 ; t<threads ; t++ )
1738 TreeOctNode::NeighborKey3 neighborKey;
1739 neighborKey.set( depth );
1740 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1742 int pIdx = sNodes.treeNodes[i]->nodeData.pointIndex;
1745 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1746 _points[ pIdx ].value = WeightedCoarserFunctionValue( neighborKey , sNodes.treeNodes[i] , metSolution );
1751 template<
int Degree >
1752 Real Octree< Degree >::WeightedCoarserFunctionValue(
const OctNode< TreeNodeData , Real >::NeighborKey3& neighborKey ,
const TreeOctNode* pointNode ,
Real* metSolution )
const
1754 int depth = pointNode->depth();
1755 if( !depth || pointNode->nodeData.pointIndex==-1 )
return Real(0.);
1756 double pointValue = 0;
1758 Real weight = _points[ pointNode->nodeData.pointIndex ].weight;
1759 Point3D< Real > p = _points[ pointNode->nodeData.pointIndex ].position;
1764 const TreeOctNode::Neighbors3& neighbors = neighborKey.neighbors[depth-1];
1765 neighbors.neighbors[1][1][1]->depthAndOffset( d , _idx );
1769 for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
for(
int l=0 ; l<3 ; l++ )
1770 if( neighbors.neighbors[j][k][l] && neighbors.neighbors[j][k][l]->nodeData.nodeIndex>=0 )
1773 const TreeOctNode* basisNode = neighbors.neighbors[j][k][l];
1774 int idx[] = { _idx[0]+j , _idx[1]+k , _idx[2]+l };
1776 fData.baseBSplines[ idx[0] ][2-j]( p[0] ) *
1777 fData.baseBSplines[ idx[1] ][2-k]( p[1] ) *
1778 fData.baseBSplines[ idx[2] ][2-l]( p[2] ) *
1779 metSolution[basisNode->nodeData.nodeIndex];
1782 return Real( pointValue * weight );
1784 template<
int Degree>
1785 int Octree<Degree>::GetFixedDepthLaplacian( SparseSymmetricMatrix< Real >& matrix ,
int depth ,
const SortedTreeNodes& sNodes ,
Real* metSolution )
1787 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1788 double stencil[5][5][5];
1789 SetLaplacianStencil( depth , stencil );
1790 Stencil< double , 5 > stencils[2][2][2];
1791 SetLaplacianStencils( depth , stencils );
1792 matrix.Resize( range );
1793 #pragma omp parallel for num_threads( threads )
1794 for(
int t=0 ; t<threads ; t++ )
1796 TreeOctNode::NeighborKey5 neighborKey5;
1797 neighborKey5.set( depth );
1798 for(
int i=(range*t)/threads ; i<(range*(t+1))/threads ; i++ )
1801 neighborKey5.getNeighbors( node );
1804 int count = GetMatrixRowSize( neighborKey5.neighbors[depth] );
1807 #pragma omp critical (matrix_set_row_size)
1809 matrix.SetRowSize( i , count );
1813 matrix.rowSizes[i] = SetMatrixRow( neighborKey5.neighbors[depth] , matrix[i] , start , stencil );
1819 c = int( node - node->parent->children );
1823 UpdateConstraintsFromCoarser( neighborKey5 , node , metSolution , stencils[x][y][z] );
1828 template<
int Degree>
1829 int Octree<Degree>::GetRestrictedFixedDepthLaplacian( SparseSymmetricMatrix< Real >& matrix,
int depth,
const int* entries,
int entryCount,
1831 const SortedTreeNodes& sNodes ,
Real* metSolution )
1833 for(
int i=0 ; i<entryCount ; i++ ) sNodes.treeNodes[ entries[i] ]->nodeData.nodeIndex = i;
1834 double stencil[5][5][5];
1835 int rDepth , rOff[3];
1836 rNode->depthAndOffset( rDepth , rOff );
1837 matrix.Resize( entryCount );
1838 SetLaplacianStencil( depth , stencil );
1839 Stencil< double , 5 > stencils[2][2][2];
1840 SetLaplacianStencils( depth , stencils );
1841 #pragma omp parallel for num_threads( threads )
1842 for(
int t=0 ; t<threads ; t++ )
1844 TreeOctNode::NeighborKey5 neighborKey5;
1845 neighborKey5.set( depth );
1846 for(
int i=(entryCount*t)/threads ; i<(entryCount*(t+1))/threads ; i++ )
1848 TreeOctNode* node = sNodes.treeNodes[ entries[i] ];
1851 off[0] >>= (depth-rDepth) , off[1] >>= (depth-rDepth) , off[2] >>= (depth-rDepth);
1852 bool isInterior = ( off[0]==rOff[0] && off[1]==rOff[1] && off[2]==rOff[2] );
1854 neighborKey5.getNeighbors( node );
1856 int xStart=0 , xEnd=5 , yStart=0 , yEnd=5 , zStart=0 , zEnd=5;
1857 if( !isInterior ) SetMatrixRowBounds( neighborKey5.neighbors[depth].neighbors[2][2][2] , rDepth , rOff , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1860 int count = GetMatrixRowSize( neighborKey5.neighbors[depth] , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1863 #pragma omp critical (matrix_set_row_size)
1865 matrix.SetRowSize( i , count );
1869 matrix.rowSizes[i] = SetMatrixRow( neighborKey5.neighbors[depth] , matrix[i] , 0 , stencil , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1874 c = int( node - node->parent->children );
1878 UpdateConstraintsFromCoarser( neighborKey5 , node , metSolution , stencils[x][y][z] );
1881 for(
int i=0 ; i<entryCount ; i++ ) sNodes.treeNodes[entries[i]]->nodeData.nodeIndex = entries[i];
1885 template<
int Degree>
1890 fData.setDotTables( fData.DD_DOT_FLAG | fData.DV_DOT_FLAG );
1893 _sNodes.treeNodes[0]->nodeData.solution = 0;
1895 std::vector< Real > metSolution( _sNodes.nodeCount[ _sNodes.maxDepth ] , 0 );
1897 for( i=1 ; i<_sNodes.maxDepth ; i++ )
1899 if( subdivideDepth>0 ) iter += SolveFixedDepthMatrix( i , _sNodes , &metSolution[0] , subdivideDepth , showResidual , minIters , accuracy );
1900 else iter += SolveFixedDepthMatrix( i , _sNodes , &metSolution[0] , showResidual , minIters , accuracy );
1903 fData.clearDotTables( fData.VV_DOT_FLAG | fData.DV_DOT_FLAG | fData.DD_DOT_FLAG );
1908 template<
int Degree>
1914 double systemTime=0. , solveTime=0. , updateTime=0. , evaluateTime = 0.;
1917 if( depth<=_minDepth ) UpSampleCoarserSolution( depth , sNodes , X );
1921 UpSample( depth-1 , sNodes , metSolution );
1924 #pragma omp parallel for num_threads( threads )
1925 for(
int i=_sNodes.nodeCount[depth-1] ; i<_sNodes.nodeCount[depth] ; i++ ) metSolution[i] += _sNodes.treeNodes[i]->nodeData.solution;
1927 if( _constrainValues )
1929 SetCoarserPointValues( depth , sNodes , metSolution );
1934 int maxECount = ( (2*Degree+1)*(2*Degree+1)*(2*Degree+1) + 1 ) / 2;
1935 maxECount = ( ( maxECount + 15 ) / 16 ) * 16;
1943 GetFixedDepthLaplacian( M , depth , sNodes , metSolution );
1955 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;
1956 double bNorm =
B.Norm( 2 ) , rNorm = (
B - M * X ).Norm( 2 );
1957 printf(
"\tResidual: (%d %g) %g -> %g (%f) [%d]\n" , M.
Entries() , sqrt(mNorm) , bNorm , rNorm , rNorm/bNorm , iter );
1965 template<
int Degree>
1966 int Octree<Degree>::SolveFixedDepthMatrix(
int depth ,
const SortedTreeNodes& sNodes ,
Real* metSolution ,
int startingDepth ,
bool showResidual ,
int minIters ,
double accuracy )
1968 if( startingDepth>=depth )
return SolveFixedDepthMatrix( depth , sNodes , metSolution , showResidual , minIters , accuracy );
1970 int i , j , d , tIter=0;
1971 SparseSymmetricMatrix< Real > _M;
1972 Vector< Real >
B , B_ , X_;
1973 AdjacencySetFunction asf;
1974 AdjacencyCountFunction acf;
1975 double systemTime = 0 , solveTime = 0 , memUsage = 0 , evaluateTime = 0 , gTime = 0, sTime = 0;
1976 Real myRadius , myRadius2;
1978 if( depth>_minDepth )
1981 UpSample( depth-1 , sNodes , metSolution );
1984 #pragma omp parallel for num_threads( threads )
1985 for(
int i=_sNodes.nodeCount[depth-1] ; i<_sNodes.nodeCount[depth] ; i++ ) metSolution[i] += _sNodes.treeNodes[i]->nodeData.solution;
1988 if( _constrainValues )
1990 SetCoarserPointValues( depth , sNodes , metSolution );
1992 B.Resize( sNodes.nodeCount[depth+1] - sNodes.nodeCount[depth] );
1995 for( i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ )
1997 B[ i-sNodes.nodeCount[depth] ] = sNodes.treeNodes[i]->nodeData.constraint;
1998 sNodes.treeNodes[i]->nodeData.constraint = 0;
2001 myRadius = 2*radius-
Real(0.5);
2004 d = depth-startingDepth;
2005 std::vector< int > subDimension( sNodes.nodeCount[d+1]-sNodes.nodeCount[d] );
2006 int maxDimension = 0;
2007 for( i=sNodes.nodeCount[d] ; i<sNodes.nodeCount[d+1] ; i++ )
2010 acf.adjacencyCount = 0;
2013 if( temp->depth()==depth ) acf.adjacencyCount++ , temp = sNodes.treeNodes[i]->nextBranch( temp );
2014 else temp = sNodes.treeNodes[i]->nextNode ( temp );
2016 for( j=sNodes.nodeCount[d] ; j<sNodes.nodeCount[d+1] ; j++ )
2018 if( i==j )
continue;
2021 subDimension[i-sNodes.nodeCount[d]] = acf.adjacencyCount;
2022 maxDimension = std::max< int >( maxDimension , subDimension[i-sNodes.nodeCount[d]] );
2024 asf.adjacencies =
new int[maxDimension];
2025 MapReduceVector< Real > mrVector;
2026 mrVector.resize( threads , maxDimension );
2028 for( i=sNodes.nodeCount[d] ; i<sNodes.nodeCount[d+1] ; i++ )
2032 acf.adjacencyCount = subDimension[i-sNodes.nodeCount[d]];
2033 if( !acf.adjacencyCount )
continue;
2036 asf.adjacencyCount = 0;
2039 if( temp->depth()==depth ) asf.adjacencies[ asf.adjacencyCount++ ] = temp->nodeData.nodeIndex , temp = sNodes.treeNodes[i]->nextBranch( temp );
2040 else temp = sNodes.treeNodes[i]->nextNode ( temp );
2042 for( j=sNodes.nodeCount[d] ; j<sNodes.nodeCount[d+1] ; j++ )
2044 if( i==j )
continue;
2049 B_.Resize( asf.adjacencyCount );
2050 for( j=0 ; j<asf.adjacencyCount ; j++ ) B_[j] =
B[ asf.adjacencies[j]-sNodes.nodeCount[depth] ];
2052 X_.Resize( asf.adjacencyCount );
2053 #pragma omp parallel for num_threads( threads ) schedule( static )
2054 for( j=0 ; j<asf.adjacencyCount ; j++ )
2056 X_[j] = sNodes.treeNodes[ asf.adjacencies[j] ]->nodeData.solution;
2060 GetRestrictedFixedDepthLaplacian( _M , depth , asf.adjacencies , asf.adjacencyCount , sNodes.treeNodes[i] , myRadius , sNodes , metSolution );
2061 #pragma omp parallel for num_threads( threads ) schedule( static )
2062 for( j=0 ; j<asf.adjacencyCount ; j++ )
2064 B_[j] += sNodes.treeNodes[asf.adjacencies[j]]->nodeData.constraint;
2065 sNodes.treeNodes[ asf.adjacencies[j] ]->nodeData.constraint = 0;
2075 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;
2076 double bNorm = B_.Norm( 2 ) , rNorm = ( B_ - _M * X_ ).Norm( 2 );
2077 printf(
"\t\tResidual: (%d %g) %g -> %g (%f) [%d]\n" , _M.Entries() , sqrt(mNorm) , bNorm , rNorm , rNorm/bNorm , iter );
2081 for( j=0 ; j<asf.adjacencyCount ; j++ )
2083 TreeOctNode* temp=sNodes.treeNodes[ asf.adjacencies[j] ];
2084 while( temp->depth()>sNodes.treeNodes[i]->depth() ) temp=temp->
parent;
2085 if( temp->nodeData.nodeIndex>=sNodes.treeNodes[i]->nodeData.nodeIndex ) sNodes.treeNodes[ asf.adjacencies[j] ]->
nodeData.solution =
Real( X_[j] );
2087 systemTime += gTime;
2089 memUsage = std::max< double >( MemoryUsage() , memUsage );
2092 delete[] asf.adjacencies;
2095 template<
int Degree>
2099 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;
2100 if( node->children )
for(
int i=0;i<
Cube::CORNERS && !hasNormals;i++) hasNormals |= HasNormals(&node->children[i],epsilon);
2104 template<
int Degree>
2107 int maxDepth = tree.maxDepth();
2109 if( temp->children && temp->d>=_minDepth )
2112 for(
int i=0 ; i<
Cube::CORNERS && !hasNormals ; i++ ) hasNormals = HasNormals( &temp->children[i] ,
EPSILON/(1<<maxDepth) );
2113 if( !hasNormals ) temp->children=NULL;
2117 template<
int Degree>
2125 fData.setDotTables( fData.VV_DOT_FLAG | fData.DV_DOT_FLAG );
2127 int maxDepth = _sNodes.maxDepth-1;
2129 zeroPoint[0] = zeroPoint[1] = zeroPoint[2] = 0;
2130 std::vector< Real > constraints( _sNodes.nodeCount[maxDepth] ,
Real(0) );
2131 std::vector< Point3D< Real > > coefficients( _sNodes.nodeCount[maxDepth] , zeroPoint );
2134 #pragma omp parallel for num_threads( threads )
2135 for(
int i=0 ; i<_sNodes.nodeCount[maxDepth+1] ; i++ ) _sNodes.treeNodes[i]->nodeData.constraint =
Real( 0. );
2138 std::vector< std::vector< Real > > _constraints( threads );
2139 for(
int t=0 ; t<threads ; t++ ) _constraints[t].resize( _sNodes.nodeCount[maxDepth] , 0 );
2141 for(
int d=maxDepth ; d>=0 ; d-- )
2144 SetDivergenceStencil( d , &stencil[0][0][0] ,
false );
2145 Stencil< Point3D< double > , 5 > stencils[2][2][2];
2146 SetDivergenceStencils( d , stencils ,
true );
2147 #pragma omp parallel for num_threads( threads )
2148 for(
int t=0 ; t<threads ; t++ )
2150 TreeOctNode::NeighborKey5 neighborKey5;
2151 neighborKey5.set( fData.depth );
2152 int start = _sNodes.nodeCount[d] , end = _sNodes.nodeCount[d+1] , range = end-start;
2153 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
2156 int startX=0 , endX=5 , startY=0 , endY=5 , startZ=0 , endZ=5;
2157 int depth = node->
depth();
2158 neighborKey5.getNeighbors( node );
2160 bool isInterior , isInterior2;
2164 int mn = 2 , mx = (1<<d)-2;
2165 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
2167 isInterior2 = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
2175 else cx = cy = cz = 0;
2176 Stencil< Point3D< double > , 5 >& _stencil = stencils[cx][cy][cz];
2180 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth];
2183 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2185 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
2189 node->
nodeData.
constraint +=
Real( stencil[x][y][z][0] * _normal[0] + stencil[x][y][z][1] * _normal[1] + stencil[x][y][z][2] * _normal[2] );
2193 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2195 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
2202 UpdateCoarserSupportBounds( neighbors5.neighbors[2][2][2] , startX , endX , startY , endY , startZ , endZ );
2206 if( normal[0]==0 && normal[1]==0 && normal[2]==0 )
continue;
2207 if( depth<maxDepth ) coefficients[i] += normal;
2211 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth-1];
2213 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2214 if( neighbors5.neighbors[x][y][z] )
2216 TreeOctNode* _node = neighbors5.neighbors[x][y][z];
2220 _constraints[t][ _node->
nodeData.
nodeIndex ] +=
Real( div[0] * normal[0] + div[1] * normal[1] + div[2] * normal[2] );
2222 else _constraints[t][ _node->
nodeData.
nodeIndex ] += GetDivergence( node , _node , normal );
2228 #pragma omp parallel for num_threads( threads ) schedule( static )
2229 for(
int i=0 ; i<_sNodes.nodeCount[maxDepth] ; i++ )
2232 for(
int t=0 ; t<threads ; t++ ) cSum += _constraints[t][i];
2233 constraints[i] = cSum;
2236 for(
int d=maxDepth-1 ; d>=0 ; d-- ) DownSample( d , _sNodes , &constraints[0] );
2239 for(
int d=0 ; d<maxDepth ; d++ ) UpSample( d , _sNodes , &coefficients[0] );
2242 #pragma omp parallel for num_threads( threads )
2243 for(
int i=0 ; i<_sNodes.nodeCount[maxDepth] ; i++ ) _sNodes.treeNodes[i]->nodeData.constraint += constraints[i];
2246 for(
int d=0 ; d<=maxDepth ; d++ )
2248 int start = _sNodes.nodeCount[d] , end = _sNodes.nodeCount[d+1] , range = end - start;
2249 Stencil< Point3D< double > , 5 > stencils[2][2][2];
2250 SetDivergenceStencils( d , stencils ,
false );
2251 #pragma omp parallel for num_threads( threads )
2252 for(
int t=0 ; t<threads ; t++ )
2254 TreeOctNode::NeighborKey5 neighborKey5;
2255 neighborKey5.set( maxDepth );
2256 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
2259 int depth = node->
depth();
2260 if( !depth )
continue;
2261 int startX=0 , endX=5 , startY=0 , endY=5 , startZ=0 , endZ=5;
2262 UpdateCoarserSupportBounds( node , startX , endX , startY , endY , startZ , endZ );
2263 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.getNeighbors( node->
parent );
2269 int mn = 4 , mx = (1<<d)-4;
2270 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
2278 else cx = cy = cz = 0;
2279 Stencil< Point3D< double > , 5 >& _stencil = stencils[cx][cy][cz];
2282 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2283 if( neighbors5.neighbors[x][y][z] )
2285 TreeOctNode* _node = neighbors5.neighbors[x][y][z];
2291 constraint +=
Real( div[0] * normal[0] + div[1] * normal[1] + div[2] * normal[2] );
2293 else constraint += GetDivergence( _node , node , coefficients[_i] );
2300 fData.clearDotTables( fData.DV_DOT_FLAG );
2303 #pragma omp parallel for num_threads( threads )
2304 for(
int t=0 ; t<threads ; t++ )
2305 for(
int i=(_sNodes.nodeCount[maxDepth+1]*t)/threads ; i<(_sNodes.nodeCount[maxDepth+1]*(t+1))/threads ; i++ )
2315 template<
int Degree>
2317 template<
int Degree>
2318 void Octree<Degree>::AdjacencySetFunction::Function(
const TreeOctNode* node1,
const TreeOctNode* node2){adjacencies[adjacencyCount++]=node1->nodeData.nodeIndex;}
2319 template<
int Degree>
2322 if( !node1->children && node1->depth()<depth ) node1->initChildren();
2324 template<
int Degree >
2325 void Octree< Degree >::FaceEdgesFunction::Function(
const TreeOctNode* node1 ,
const TreeOctNode* node2 )
2333 for(
int j=0 ; j<count ; j++ )
2334 for(
int k=0 ; k<3 ; k++ )
2336 if( GetRootIndex( node1 , isoTri[j*3+k] , maxDepth , ri1 ) && GetRootIndex( node1 , isoTri[j*3+((k+1)%3)] , maxDepth , ri2 ) )
2338 long long key1=ri1.key , key2=ri2.key;
2339 edges->push_back( std::pair< RootInfo , RootInfo >( ri2 , ri1 ) );
2340 if( vertexCount->count( key1 )==0 )
2342 (*vertexCount)[key1].first = ri1;
2343 (*vertexCount)[key1].second=0;
2345 if( vertexCount->count( key2 )==0 )
2347 (*vertexCount)[key2].first = ri2;
2348 (*vertexCount)[key2].second=0;
2350 (*vertexCount)[key1].second--;
2351 (*vertexCount)[key2].second++;
2353 else fprintf( stderr ,
"Bad Edge 1: %d %d\n" , ri1.key , ri2.key );
2357 template<
int Degree >
2367 bool flags[3][3][3];
2368 int maxDepth = tree.maxDepth();
2371 if( subdivideDepth<=0 ) sDepth = 0;
2372 else sDepth = maxDepth-subdivideDepth;
2373 if( sDepth<=0 )
return;
2377 TreeOctNode::NeighborKey3 nKey;
2378 nKey.set( maxDepth );
2380 if( leaf->depth()>sDepth )
2382 int d , off[3] , _off[3];
2383 leaf->depthAndOffset( d , off );
2384 int res = (1<<d)-1 , _res = ( 1<<(d-sDepth) )-1;
2385 _off[0] = off[0]&_res , _off[1] = off[1]&_res , _off[2] = off[2]&_res;
2386 bool boundary[3][2] =
2388 { (off[0]!=0 && _off[0]==0) , (off[0]!=res && _off[0]==_res) } ,
2389 { (off[1]!=0 && _off[1]==0) , (off[1]!=res && _off[1]==_res) } ,
2390 { (off[2]!=0 && _off[2]==0) , (off[2]!=res && _off[2]==_res) }
2393 if( boundary[0][0] || boundary[0][1] || boundary[1][0] || boundary[1][1] || boundary[2][0] || boundary[2][1] )
2395 TreeOctNode::Neighbors3& neighbors = nKey.getNeighbors( leaf );
2396 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;
2397 int x=0 , y=0 , z=0;
2398 if ( boundary[0][0] && !neighbors.neighbors[0][1][1] ) x = -1;
2399 else if( boundary[0][1] && !neighbors.neighbors[2][1][1] ) x = 1;
2400 if ( boundary[1][0] && !neighbors.neighbors[1][0][1] ) y = -1;
2401 else if( boundary[1][1] && !neighbors.neighbors[1][2][1] ) y = 1;
2402 if ( boundary[2][0] && !neighbors.neighbors[1][1][0] ) z = -1;
2403 else if( boundary[2][1] && !neighbors.neighbors[1][1][2] ) z = 1;
2408 if( x && y && z ) flags[1+x][1+y][1+z] =
true;
2410 if( x && y ) flags[1+x][1+y][1 ] =
true;
2411 if( x && z ) flags[1+x][1 ][1+z] =
true;
2412 if( y && z ) flags[1 ][1+y][1+1] =
true;
2414 if( x ) flags[1+x][1 ][1 ] =
true;
2415 if( y ) flags[1 ][1+y][1 ] =
true;
2416 if( z ) flags[1 ][1 ][1+z] =
true;
2417 nKey.setNeighbors( leaf , flags );
2421 _sNodes.set( tree );
2424 template<
int Degree>
2427 fData.setValueTables( fData.VALUE_FLAG | fData.D_VALUE_FLAG , 0 , postNormalSmooth );
2430 RefineBoundary( subdivideDepth );
2432 RootData rootData , coarseRootData;
2433 std::vector< Point3D< float > >* interiorPoints;
2434 int maxDepth = tree.maxDepth();
2436 int sDepth = subdivideDepth<=0 ? 0 : std::max< int >( 0 , maxDepth-subdivideDepth );
2438 std::vector< Real > metSolution( _sNodes.nodeCount[maxDepth] , 0 );
2439 #pragma omp parallel for num_threads( threads )
2440 for(
int i=_sNodes.nodeCount[_minDepth] ; i<_sNodes.nodeCount[maxDepth] ; i++ ) metSolution[i] = _sNodes.treeNodes[i]->nodeData.solution;
2441 for(
int d=0 ; d<maxDepth ; d++ ) UpSample( d , _sNodes , &metSolution[0] );
2444 #pragma omp parallel for num_threads( threads )
2445 for(
int i=0 ; i<_sNodes.nodeCount[maxDepth+1] ; i++ ) _sNodes.treeNodes[i]->nodeData.mcIndex = 0;
2447 rootData.boundaryValues =
new std::unordered_map<
long long , std::pair<
Real ,
Point3D< Real > > >();
2450 int maxCCount = _sNodes.getMaxCornerCount( &tree , sDepth , maxDepth , threads );
2451 int maxECount = _sNodes.getMaxEdgeCount ( &tree , sDepth , threads );
2452 rootData.cornerValues =
new Real [ maxCCount ];
2454 rootData.interiorRoots =
new int [ maxECount ];
2455 rootData.cornerValuesSet =
new char[ maxCCount ];
2456 rootData.cornerNormalsSet =
new char[ maxCCount ];
2457 rootData.edgesSet =
new char[ maxECount ];
2458 _sNodes.setCornerTable( coarseRootData , &tree , sDepth , threads );
2459 coarseRootData.cornerValues =
new Real[ coarseRootData.cCount ];
2460 coarseRootData.cornerNormals =
new Point3D< Real >[ coarseRootData.cCount ];
2461 coarseRootData.cornerValuesSet =
new char[ coarseRootData.cCount ];
2462 coarseRootData.cornerNormalsSet =
new char[ coarseRootData.cCount ];
2463 memset( coarseRootData.cornerValuesSet , 0 ,
sizeof(
char ) * coarseRootData.cCount );
2464 memset( coarseRootData.cornerNormalsSet , 0 ,
sizeof(
char ) * coarseRootData.cCount );
2467 std::vector< TreeOctNode::ConstNeighborKey3 > nKeys( threads );
2468 for(
int t=0 ; t<threads ; t++ ) nKeys[t].set( maxDepth );
2469 TreeOctNode::ConstNeighborKey3 nKey;
2470 std::vector< TreeOctNode::ConstNeighborKey5 > nKeys5( threads );
2471 for(
int t=0 ; t<threads ; t++ ) nKeys5[t].set( maxDepth );
2472 TreeOctNode::ConstNeighborKey5 nKey5;
2473 nKey5.set( maxDepth );
2474 nKey.set( maxDepth );
2476 for(
int i=_sNodes.nodeCount[sDepth] ; i<_sNodes.nodeCount[sDepth+1] ; i++ )
2478 if( !_sNodes.treeNodes[i]->children )
continue;
2480 _sNodes.setCornerTable( rootData , _sNodes.treeNodes[i] , threads );
2481 _sNodes.setEdgeTable ( rootData , _sNodes.treeNodes[i] , threads );
2482 memset( rootData.cornerValuesSet , 0 ,
sizeof(
char ) * rootData.cCount );
2483 memset( rootData.cornerNormalsSet , 0 ,
sizeof(
char ) * rootData.cCount );
2484 memset( rootData.edgesSet , 0 ,
sizeof(
char ) * rootData.eCount );
2485 interiorPoints =
new std::vector< Point3D< float > >();
2486 for(
int d=maxDepth ; d>sDepth ; d-- )
2488 int leafNodeCount = 0;
2489 std::vector< TreeOctNode* > leafNodes;
2490 for(
TreeOctNode* node=_sNodes.treeNodes[i]->
nextLeaf() ; node ; node=_sNodes.treeNodes[i]->
nextLeaf( node ) )
if( node->d==d ) leafNodeCount++;
2491 leafNodes.reserve( leafNodeCount );
2492 for(
TreeOctNode* node=_sNodes.treeNodes[i]->
nextLeaf() ; node ; node=_sNodes.treeNodes[i]->
nextLeaf( node ) )
if( node->d==d ) leafNodes.push_back( node );
2493 Stencil< double , 3 > stencil1[8] , stencil2[8][8];
2494 SetEvaluationStencils( d , stencil1 , stencil2 );
2497 #pragma omp parallel for num_threads( threads )
2498 for(
int t=0 ; t<threads ; t++ )
for(
int i=(leafNodeCount*t)/threads ; i<(leafNodeCount*(t+1))/threads ; i++ )
2501 SetIsoCorners( isoValue , leaf , rootData , rootData.cornerValuesSet , rootData.cornerValues , nKeys[t] , &metSolution[0] , stencil1 , stencil2 );
2506 int res = 1<<(d-sDepth);
2507 off[0] %= res , off[1] %=res , off[2] %= res;
2509 if( !(off[0]%res) && !(off[1]%res) && !(off[2]%res) )
2512 while( temp->
d!=sDepth ) temp = temp->
parent;
2513 int x = off[0]==0 ? 0 : 1 , y = off[1]==0 ? 0 : 1 , z = off[2]==0 ? 0 : 1;
2515 int idx = coarseRootData.cornerIndices( temp )[ c ];
2516 coarseRootData.cornerValues[ idx ] = rootData.cornerValues[ rootData.cornerIndices( leaf )[c] ];
2517 coarseRootData.cornerValuesSet[ idx ] =
true;
2522 SetMCRootPositions( leaf , sDepth , isoValue , nKeys5[t] , rootData , interiorPoints , mesh , &metSolution[0] , nonLinearFit );
2528 std::vector< Point3D< float > > barycenters;
2529 std::vector< Point3D< float > >* barycenterPtr = addBarycenter ? & barycenters : NULL;
2530 #endif // MISHA_DEBUG
2531 #pragma omp parallel for num_threads( threads )
2532 for(
int t=0 ; t<threads ; t++ )
for(
int i=(leafNodeCount*t)/threads ; i<(leafNodeCount*(t+1))/threads ; i++ )
2537 GetMCIsoTriangles( leaf , mesh , rootData , interiorPoints , offSet , sDepth , polygonMesh , barycenterPtr );
2538 #else // !MISHA_DEBUG
2539 GetMCIsoTriangles( leaf , mesh , rootData , interiorPoints , offSet , sDepth , addBarycenter , polygonMesh );
2540 #endif // MISHA_DEBUG
2543 for(
int i=0 ; i<barycenters.size() ; i++ ) interiorPoints->push_back( barycenters[i] );
2544 #endif // MISHA_DEBUG
2548 delete interiorPoints;
2552 delete[] rootData.cornerValues ,
delete[] rootData.cornerNormals , rootData.cornerValues = NULL , rootData.cornerNormals = NULL;
2553 delete[] rootData.cornerValuesSet ,
delete[] rootData.cornerNormalsSet , rootData.cornerValuesSet = NULL , rootData.cornerNormalsSet = NULL;
2554 delete[] rootData.interiorRoots ; rootData.interiorRoots = NULL;
2555 delete[] rootData.edgesSet ; rootData.edgesSet = NULL;
2556 coarseRootData.interiorRoots = NULL;
2557 coarseRootData.boundaryValues = rootData.boundaryValues;
2558 for(
auto iter=rootData.boundaryRoots.cbegin() ; iter!=rootData.boundaryRoots.cend() ; iter++ )
2559 coarseRootData.boundaryRoots[iter->first] = iter->second;
2561 for(
int d=sDepth ; d>=0 ; d-- )
2563 Stencil< double , 3 > stencil1[8] , stencil2[8][8];
2564 SetEvaluationStencils( d , stencil1 , stencil2 );
2566 std::vector< Point3D< float > > barycenters;
2567 std::vector< Point3D< float > >* barycenterPtr = addBarycenter ? &barycenters : NULL;
2568 #endif // MISHA_DEBUG
2569 for(
int i=_sNodes.nodeCount[d] ; i<_sNodes.nodeCount[d+1] ; i++ )
2575 SetIsoCorners( isoValue , leaf , coarseRootData , coarseRootData.cornerValuesSet , coarseRootData.cornerValues , nKey , &metSolution[0] , stencil1 , stencil2 );
2580 SetMCRootPositions( leaf , 0 , isoValue , nKey5 , coarseRootData , NULL , mesh , &metSolution[0] , nonLinearFit );
2582 GetMCIsoTriangles( leaf , mesh , coarseRootData , NULL , 0 , 0 , polygonMesh , barycenterPtr );
2583 #else // !MISHA_DEBUG
2584 GetMCIsoTriangles( leaf , mesh , coarseRootData , NULL , 0 , 0 , addBarycenter , polygonMesh );
2585 #endif // MISHA_DEBUG
2591 delete[] coarseRootData.cornerValues ,
delete[] coarseRootData.cornerNormals;
2592 delete[] coarseRootData.cornerValuesSet ,
delete[] coarseRootData.cornerNormalsSet;
2593 delete rootData.boundaryValues;
2595 template<
int Degree>
2601 idx[0]*=fData.functionCount;
2602 idx[1]*=fData.functionCount;
2603 idx[2]*=fData.functionCount;
2604 int minDepth = std::max< int >( 0 , std::min< int >( _minDepth , node->
depth()-1 ) );
2605 for(
int i=minDepth ; i<=node->
depth() ; i++ )
2606 for(
int j=0;j<3;j++)
2607 for(
int k=0;k<3;k++)
2608 for(
int l=0;l<3;l++)
2615 fData.valueTables[idx[0]+
int(n->
off[0])]*
2616 fData.valueTables[idx[1]+
int(n->
off[1])]*
2617 fData.valueTables[idx[2]+
int(n->
off[2])]);
2627 fData.valueTables[idx[0]+
int(n->off[0])]*
2628 fData.valueTables[idx[1]+
int(n->off[1])]*
2629 fData.valueTables[idx[2]+
int(n->off[2])]);
2630 if( n->children ) n=&n->children[ii];
2637 template<
int Degree >
2638 Real Octree< Degree >::getCornerValue(
const OctNode< TreeNodeData , Real >::ConstNeighborKey3& neighborKey3 ,
const TreeOctNode* node ,
int corner ,
const Real* metSolution )
2644 idx[0] *= fData.functionCount;
2645 idx[1] *= fData.functionCount;
2646 idx[2] *= fData.functionCount;
2648 int d = node->depth();
2650 int startX = 0 , endX = 3 , startY = 0 , endY = 3 , startZ = 0 , endZ = 3;
2653 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d];
2654 if( cx==0 ) endX = 2;
2656 if( cy==0 ) endY = 2;
2658 if( cz==0 ) endZ = 2;
2660 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2662 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2666 fData.valueTables[ idx[0]+int(n->off[0]) ]*
2667 fData.valueTables[ idx[1]+int(n->off[1]) ]*
2668 fData.valueTables[ idx[2]+int(n->off[2]) ];
2673 if( d>0 && d>_minDepth )
2675 int _corner = int( node - node->parent->children );
2676 int _cx , _cy , _cz;
2678 if( cx!=_cx ) startX = 0 , endX = 3;
2679 if( cy!=_cy ) startY = 0 , endY = 3;
2680 if( cz!=_cz ) startZ = 0 , endZ = 3;
2681 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d-1];
2682 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2684 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2688 fData.valueTables[ idx[0]+int(n->off[0]) ]*
2689 fData.valueTables[ idx[1]+int(n->off[1]) ]*
2690 fData.valueTables[ idx[2]+int(n->off[2]) ];
2691 value += metSolution[ n->nodeData.nodeIndex ] * v;
2695 return Real( value );
2697 template<
int Degree >
2698 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] )
2701 int d = node->depth();
2703 int startX = 0 , endX = 3 , startY = 0 , endY = 3 , startZ = 0 , endZ = 3;
2706 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d];
2707 if( cx==0 ) endX = 2;
2709 if( cy==0 ) endY = 2;
2711 if( cz==0 ) endZ = 2;
2713 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2715 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2719 if( d>0 && d>_minDepth )
2721 int _corner = int( node - node->parent->children );
2722 int _cx , _cy , _cz;
2724 if( cx!=_cx ) startX = 0 , endX = 3;
2725 if( cy!=_cy ) startY = 0 , endY = 3;
2726 if( cz!=_cz ) startZ = 0 , endZ = 3;
2727 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d-1];
2728 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2730 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2731 if( n ) value += metSolution[ n->nodeData.nodeIndex ] * stencil2[x][y][z];
2734 return Real( value );
2736 template<
int Degree >
2737 Point3D< Real > Octree< Degree >::getCornerNormal(
const OctNode< TreeNodeData , Real >::ConstNeighborKey5& neighborKey5 ,
const TreeOctNode* node ,
int corner ,
const Real* metSolution )
2740 Point3D< Real > normal;
2741 normal[0] = normal[1] = normal[2] = 0.;
2744 idx[0] *= fData.functionCount;
2745 idx[1] *= fData.functionCount;
2746 idx[2] *= fData.functionCount;
2748 int d = node->depth();
2751 TreeOctNode::ConstNeighbors5& neighbors = neighborKey5.neighbors[d];
2752 for(
int j=0 ; j<5 ; j++ )
for(
int k=0 ; k<5 ; k++ )
for(
int l=0 ; l<5 ; l++ )
2754 const TreeOctNode* n=neighbors.neighbors[j][k][l];
2757 int _idx[] = { idx[0] + n->
off[0] , idx[1] + n->off[1] , idx[2] + n->off[2] };
2758 double values[] = { fData.valueTables[_idx[0]] , fData.valueTables[_idx[1]] , fData.valueTables[_idx[2]] };
2759 double dValues[] = { fData.dValueTables[_idx[0]] , fData.dValueTables[_idx[1]] , fData.dValueTables[_idx[2]] };
2760 Real solution = n->nodeData.solution;
2761 normal[0] +=
Real( dValues[0] * values[1] * values[2] * solution );
2762 normal[1] +=
Real( values[0] * dValues[1] * values[2] * solution );
2763 normal[2] +=
Real( values[0] * values[1] * dValues[2] * solution );
2767 if( d>0 && d>_minDepth )
2769 TreeOctNode::ConstNeighbors5& neighbors = neighborKey5.neighbors[d-1];
2770 for(
int j=0 ; j<5 ; j++ )
for(
int k=0 ; k<5 ; k++ )
for(
int l=0 ; l<5 ; l++ )
2772 const TreeOctNode* n=neighbors.neighbors[j][k][l];
2775 int _idx[] = { idx[0] + n->
off[0] , idx[1] + n->off[1] , idx[2] + n->off[2] };
2776 double values[] = { fData.valueTables[_idx[0]] , fData.valueTables[_idx[1]] , fData.valueTables[_idx[2]] };
2777 double dValues[] = { fData.dValueTables[_idx[0]] , fData.dValueTables[_idx[1]] , fData.dValueTables[_idx[2]] };
2778 Real solution = metSolution[ n->nodeData.nodeIndex ];
2779 normal[0] +=
Real( dValues[0] * values[1] * values[2] * solution );
2780 normal[1] +=
Real( values[0] * dValues[1] * values[2] * solution );
2781 normal[2] +=
Real( values[0] * values[1] * dValues[2] * solution );
2787 template<
int Degree >
2790 Real isoValue , weightSum;
2792 neighborKey2.set( fData.depth );
2793 fData.setValueTables( fData.VALUE_FLAG , 0 );
2795 isoValue = weightSum = 0;
2796 #pragma omp parallel for num_threads( threads ) reduction( + : isoValue , weightSum )
2797 for(
int t=0 ; t<threads ; t++)
2799 TreeOctNode::ConstNeighborKey3 nKey;
2800 nKey.set( _sNodes.maxDepth-1 );
2801 int nodeCount = _sNodes.nodeCount[ _sNodes.maxDepth ];
2802 for(
int i=(nodeCount*t)/threads ; i<(nodeCount*(t+1))/threads ; i++ )
2805 nKey.getNeighbors( temp );
2809 isoValue += getCenterValue( nKey , temp ) * w;
2814 return isoValue/weightSum;
2817 template<
int Degree >
2818 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] )
2826 int mn = 2 , mx = (1<<d)-2;
2827 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
2828 nKey.getNeighbors( leaf );
2831 int vIndex = cIndices[c];
2832 if( valuesSet[vIndex] ) cornerValues[c] = values[vIndex];
2835 if( isInterior ) cornerValues[c] = getCornerValue( nKey , leaf , c , metSolution , stencil1[c].values , stencil2[
int(leaf - leaf->
parent->
children)][c].values );
2836 else cornerValues[c] = getCornerValue( nKey , leaf , c , metSolution );
2837 values[vIndex] = cornerValues[c];
2838 valuesSet[vIndex] = 1;
2857 if( parent->parent && parent->parent->d>=_minDepth && (parent-parent->parent->children)==c )
2860 parent = parent->parent;
2869 template<
int Degree>
2870 int Octree<Degree>::InteriorFaceRootCount(
const TreeOctNode* node,
const int &faceIndex,
int maxDepth){
2871 int c1,c2,e1,e2,dir,off,cnt=0;
2892 cnt+=EdgeRootCount(&node->children[c1],e1,maxDepth)+EdgeRootCount(&node->children[c1],e2,maxDepth);
2907 cnt+=EdgeRootCount(&node->children[c2],e1,maxDepth)+EdgeRootCount(&node->children[c2],e2,maxDepth);
2908 for(
int i=0;i<
Cube::CORNERS/2;i++){
if(node->children[corners[i]].children){cnt+=InteriorFaceRootCount(&node->children[corners[i]],faceIndex,maxDepth);}}
2913 template<
int Degree>
2914 int Octree<Degree>::EdgeRootCount(
const TreeOctNode* node,
int edgeIndex,
int maxDepth){
2922 if(node->depth()<maxDepth){
2924 if(temp && temp->children){
2929 temp=node->faceNeighbor(f2);
2930 if(temp && temp->children){
2935 temp=node->edgeNeighbor(edgeIndex);
2936 if(temp && temp->children){
2945 if(finest->children)
return EdgeRootCount(&finest->children[c1],eIndex,maxDepth)+EdgeRootCount(&finest->children[c2],eIndex,maxDepth);
2948 template<
int Degree>
2949 int Octree<Degree>::IsBoundaryFace(
const TreeOctNode* node,
int faceIndex,
int subdivideDepth){
2950 int dir,offset,d,o[3],idx;
2952 if(subdivideDepth<0){
return 0;}
2953 if(node->d<=subdivideDepth){
return 1;}
2955 node->depthAndOffset(d,o);
2957 idx=(int(o[dir])<<1) + (offset<<1);
2958 return !(idx%(2<<(int(node->d)-subdivideDepth)));
2960 template<
int Degree>
2961 int Octree<Degree>::IsBoundaryEdge(
const TreeOctNode* node,
int edgeIndex,
int subdivideDepth){
2964 return IsBoundaryEdge(node,dir,x,y,subdivideDepth);
2966 template<
int Degree>
2967 int Octree<Degree>::IsBoundaryEdge(
const TreeOctNode* node ,
int dir ,
int x ,
int y ,
int subdivideDepth )
2969 int d , o[3] , idx1 , idx2 , mask;
2971 if( subdivideDepth<0 )
return 0;
2972 if( node->d<=subdivideDepth )
return 1;
2973 node->depthAndOffset( d , o );
2990 mask = 1<<( int(node->d) - subdivideDepth );
2991 return !(idx1%(mask)) || !(idx2%(mask));
2993 template<
int Degree >
3001 ri.node->centerAndWidth( c , width );
3005 start[0] = c[0] - width/2;
3006 end [0] = c[0] + width/2;
3007 start[1] = end[1] = c[1] - width/2 + width*i1;
3008 start[2] = end[2] = c[2] - width/2 + width*i2;
3011 start[0] = end[0] = c[0] - width/2 + width*i1;
3012 start[1] = c[1] - width/2;
3013 end [1] = c[1] + width/2;
3014 start[2] = end[2] = c[2] - width/2 + width*i2;
3017 start[0] = end[0] = c[0] - width/2 + width*i1;
3018 start[1] = end[1] = c[1] - width/2 + width*i2;
3019 start[2] = c[2] - width/2;
3020 end [2] = c[2] + width/2;
3027 template<
int Degree >
3028 int Octree< Degree >::GetRoot(
const RootInfo& ri ,
Real isoValue , TreeOctNode::ConstNeighborKey5& neighborKey5 , Point3D< Real > & position , RootData& rootData ,
int sDepth ,
const Real* metSolution ,
int nonLinearFit )
3035 long long key1 , key2;
3036 Point3D< Real > n[2];
3038 int i , o , i1 , i2 , rCount=0;
3040 std::vector< double > roots;
3042 Real center , width;
3045 int idx1[3] , idx2[3];
3049 bool isBoundary = ( IsBoundaryEdge( ri.node , ri.edgeIndex , sDepth )!=0 );
3050 bool haveKey1 , haveKey2;
3051 std::pair< Real , Point3D< Real > > keyValue1 , keyValue2;
3054 iter1 = rootData.cornerIndices( ri.node )[ c1 ];
3055 iter2 = rootData.cornerIndices( ri.node )[ c2 ];
3056 keyValue1.first = rootData.cornerValues[iter1];
3057 keyValue2.first = rootData.cornerValues[iter2];
3060 #pragma omp critical (normal_hash_access)
3062 haveKey1 = ( rootData.boundaryValues->count( key1 )>0 );
3063 haveKey2 = ( rootData.boundaryValues->count( key2 )>0 );
3064 if( haveKey1 ) keyValue1 = (*rootData.boundaryValues)[key1];
3065 if( haveKey2 ) keyValue2 = (*rootData.boundaryValues)[key2];
3070 haveKey1 = ( rootData.cornerNormalsSet[ iter1 ] != 0 );
3071 haveKey2 = ( rootData.cornerNormalsSet[ iter2 ] != 0 );
3072 keyValue1.first = rootData.cornerValues[iter1];
3073 keyValue2.first = rootData.cornerValues[iter2];
3074 if( haveKey1 ) keyValue1.second = rootData.cornerNormals[iter1];
3075 if( haveKey2 ) keyValue2.second = rootData.cornerNormals[iter2];
3078 if( !haveKey1 || !haveKey2 ) neighborKey5.getNeighbors( ri.node );
3079 if( !haveKey1 ) keyValue1.second = getCornerNormal( neighborKey5 , ri.node , c1 , metSolution );
3080 x0 = keyValue1.first;
3081 n[0] = keyValue1.second;
3083 if( !haveKey2 ) keyValue2.second = getCornerNormal( neighborKey5 , ri.node , c2 , metSolution );
3084 x1 = keyValue2.first;
3085 n[1] = keyValue2.second;
3087 if( !haveKey1 || !haveKey2 )
3091 #pragma omp critical (normal_hash_access)
3093 if( !haveKey1 ) (*rootData.boundaryValues)[key1] = keyValue1;
3094 if( !haveKey2 ) (*rootData.boundaryValues)[key2] = keyValue2;
3099 if( !haveKey1 ) rootData.cornerNormals[iter1] = keyValue1.second , rootData.cornerNormalsSet[ iter1 ] = 1;
3100 if( !haveKey2 ) rootData.cornerNormals[iter2] = keyValue2.second , rootData.cornerNormalsSet[ iter2 ] = 1;
3105 ri.node->centerAndWidth(c,width);
3107 for( i=0 ; i<DIMENSION ; i++ ) n[0][i] *= width , n[1][i] *= width;
3112 position[1] = c[1]-width/2+width*i1;
3113 position[2] = c[2]-width/2+width*i2;
3116 position[0] = c[0]-width/2+width*i1;
3117 position[2] = c[2]-width/2+width*i2;
3120 position[0] = c[0]-width/2+width*i1;
3121 position[1] = c[1]-width/2+width*i2;
3129 double scl=(x1-x0)/((dx1+dx0)/2);
3134 P.coefficients[0] = x0;
3135 P.coefficients[1] = dx0;
3136 P.coefficients[2] = 3*(x1-x0)-dx1-2*dx0;
3138 P.getSolutions( isoValue , roots ,
EPSILON );
3139 for( i=0 ; i<int(roots.size()) ; i++ )
3140 if( roots[i]>=0 && roots[i]<=1 )
3142 averageRoot +=
Real( roots[i] );
3145 if( rCount && nonLinearFit ) averageRoot /= rCount;
3146 else averageRoot =
Real((x0-isoValue)/(x0-x1));
3147 if( averageRoot<0 || averageRoot>1 )
3149 fprintf( stderr ,
"[WARNING] Bad average root: %f\n" , averageRoot );
3150 fprintf( stderr ,
"\t(%f %f) , (%f %f) (%f)\n" , x0 , x1 , dx0 , dx1 , isoValue );
3151 if( averageRoot<0 ) averageRoot = 0;
3152 if( averageRoot>1 ) averageRoot = 1;
3154 position[o] =
Real(center-width/2+width*averageRoot);
3157 template<
int Degree >
3158 int Octree< Degree >::GetRootIndex(
const TreeOctNode* node ,
int edgeIndex ,
int maxDepth ,
int sDepth,RootInfo& ri )
3167 finestIndex=edgeIndex;
3168 if(node->depth()<maxDepth){
3169 if(IsBoundaryFace(node,f1,sDepth)){temp=NULL;}
3170 else{temp=node->faceNeighbor(f1);}
3171 if(temp && temp->children){
3176 if(IsBoundaryFace(node,f2,sDepth)){temp=NULL;}
3177 else{temp=node->faceNeighbor(f2);}
3178 if(temp && temp->children){
3183 if(IsBoundaryEdge(node,edgeIndex,sDepth)){temp=NULL;}
3184 else{temp=node->edgeNeighbor(edgeIndex);}
3185 if(temp && temp->children){
3194 if(finest->children){
3195 if (GetRootIndex(&finest->children[c1],finestIndex,maxDepth,sDepth,ri)) {
return 1;}
3196 else if (GetRootIndex(&finest->children[c2],finestIndex,maxDepth,sDepth,ri)) {
return 1;}
3199 fprintf( stderr ,
"[WARNING] Couldn't find root index with either child\n" );
3207 fprintf( stderr ,
"[WARNING] Finest node does not have iso-edge\n" );
3214 finest->depthAndOffset(d,off);
3216 ri.edgeIndex=finestIndex;
3217 int eIndex[2],offset;
3234 ri.key = (
long long)(o) | (
long long)(eIndex[0])<<5 | (
long long)(eIndex[1])<<25 | (
long long)(offset)<<45;
3238 template<
int Degree>
3239 int Octree<Degree>::GetRootIndex(
const TreeOctNode* node ,
int edgeIndex ,
int maxDepth , RootInfo& ri )
3252 finestIndex = edgeIndex;
3253 if( node->depth()<maxDepth && !node->children )
3255 temp=node->faceNeighbor( f1 );
3259 temp = node->faceNeighbor( f2 );
3263 temp = node->edgeNeighbor( edgeIndex );
3270 if( finest->children )
3272 if ( GetRootIndex( finest->children + c1 , finestIndex , maxDepth , ri ) )
return 1;
3273 else if( GetRootIndex( finest->children + c2 , finestIndex , maxDepth , ri ) )
return 1;
3276 int d1 , off1[3] , d2 , off2[3];
3277 node->depthAndOffset( d1 , off1 );
3278 finest->depthAndOffset( d2 , off2 );
3279 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 );
3281 for(
int i=0 ; i<8 ; i++ )
if( node->nodeData.mcIndex & (1<<i) ) printf(
"1" );
else printf(
"0" );
3283 for(
int i=0 ; i<8 ; i++ )
if( finest->nodeData.mcIndex & (1<<i) ) printf(
"1" );
else printf(
"0" );
3293 finest->depthAndOffset(d,off);
3295 ri.edgeIndex=finestIndex;
3296 int offset,eIndex[2];
3312 ri.key= (
long long)(o) | (
long long)(eIndex[0])<<5 | (
long long)(eIndex[1])<<25 | (
long long)(offset)<<45;
3316 template<
int Degree>
3317 int Octree<Degree>::GetRootPair(
const RootInfo& ri,
int maxDepth,RootInfo& pair){
3321 while(node->parent){
3322 c=int(node-node->parent->children);
3323 if(c!=c1 && c!=c2){
return 0;}
3325 if(c==c1){
return GetRootIndex(&node->parent->children[c2],ri.edgeIndex,maxDepth,pair);}
3326 else{
return GetRootIndex(&node->parent->children[c1],ri.edgeIndex,maxDepth,pair);}
3333 template<
int Degree>
3334 int Octree< Degree >::GetRootIndex(
const RootInfo& ri , RootData& rootData , CoredPointIndex& index )
3336 long long key = ri.key;
3337 auto rootIter = rootData.boundaryRoots.find( key );
3338 if( rootIter!=rootData.boundaryRoots.end() )
3341 index.index = rootIter->second;
3344 else if( rootData.interiorRoots )
3346 int eIndex = rootData.edgeIndices( ri.node )[ ri.edgeIndex ];
3347 if( rootData.edgesSet[ eIndex ] )
3350 index.index = rootData.interiorRoots[ eIndex ];
3356 template<
int Degree >
3357 int Octree< Degree >::SetMCRootPositions(
TreeOctNode* node ,
int sDepth ,
Real isoValue , TreeOctNode::ConstNeighborKey5& neighborKey5 , RootData& rootData ,
3358 std::vector<
Point3D< float > >* interiorPositions , CoredMeshData* mesh ,
const Real* metSolution ,
int nonLinearFit )
3360 Point3D< Real > position;
3365 for(
int i=0 ; i<DIMENSION ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
3369 if( GetRootIndex( node , eIndex , fData.depth , ri ) )
3372 if( !rootData.interiorRoots || IsBoundaryEdge( node , i , j , k , sDepth ) )
3374 std::unordered_map< long long , int >::iterator iter , end;
3376 #pragma omp critical (boundary_roots_hash_access)
3378 iter = rootData.boundaryRoots.find( key );
3379 end = rootData.boundaryRoots.end();
3384 GetRoot( ri , isoValue , neighborKey5 , position , rootData , sDepth , metSolution , nonLinearFit );
3386 #pragma omp critical (boundary_roots_hash_access)
3388 iter = rootData.boundaryRoots.find( key );
3389 end = rootData.boundaryRoots.end();
3392 mesh->inCorePoints.push_back( position );
3393 rootData.boundaryRoots[key] = int( mesh->inCorePoints.size() ) - 1;
3396 if( iter==end ) count++;
3401 int nodeEdgeIndex = rootData.edgeIndices( ri.node )[ ri.edgeIndex ];
3402 if( !rootData.edgesSet[ nodeEdgeIndex ] )
3405 GetRoot( ri , isoValue , neighborKey5 , position , rootData , sDepth , metSolution , nonLinearFit );
3407 #pragma omp critical (add_point_access)
3409 if( !rootData.edgesSet[ nodeEdgeIndex ] )
3411 rootData.interiorRoots[ nodeEdgeIndex ] = mesh->addOutOfCorePoint( position );
3412 interiorPositions->push_back( position );
3413 rootData.edgesSet[ nodeEdgeIndex ] = 1;
3423 template<
int Degree>
3424 int Octree< Degree >::SetBoundaryMCRootPositions(
int sDepth ,
Real isoValue , RootData& rootData , CoredMeshData* mesh ,
int nonLinearFit )
3426 Point3D< Real > position;
3427 int i,j,k,eIndex,hits=0;
3438 for( i=0 ; i<DIMENSION ; i++ )
3439 for( j=0 ; j<2 ; j++ )
3440 for( k=0 ; k<2 ; k++ )
3441 if( IsBoundaryEdge( node , i , j , k , sDepth ) )
3446 if( GetRootIndex( node , eIndex , fData.depth , ri ) )
3449 if( rootData.boundaryRoots.find(key)==rootData.boundaryRoots.end() )
3451 GetRoot( ri , isoValue , position , rootData , sDepth , nonLinearFit );
3452 mesh->inCorePoints.push_back( position );
3453 rootData.boundaryRoots[key] = int( mesh->inCorePoints.size() )-1;
3459 if( hits ) node=tree.nextLeaf(node);
3460 else node=tree.nextBranch(node);
3464 template<
int Degree>
3465 void Octree< Degree >::GetMCIsoEdges(
TreeOctNode* node ,
int sDepth , std::vector< std::pair< RootInfo , RootInfo > >& edges )
3468 int count=0 , tris=0;
3470 FaceEdgesFunction fef;
3472 std::unordered_map< long long , std::pair< RootInfo , int > > vertexCount;
3475 fef.maxDepth = fData.depth;
3476 fef.vertexCount = &vertexCount;
3482 temp = node->faceNeighbor( fIndex );
3485 if( temp && temp->children && !IsBoundaryFace( node , fIndex , sDepth ) ) temp->processNodeFaces( temp , &fef , ref );
3490 for(
int j=0 ; j<count ; j++ )
3491 for(
int k=0 ; k<3 ; k++ )
3493 if( GetRootIndex( node , isoTri[j*3+k] , fData.depth , ri1 ) && GetRootIndex( node , isoTri[j*3+((k+1)%3)] , fData.depth , ri2 ) )
3495 long long key1 = ri1.key , key2 = ri2.key;
3496 edges.push_back( std::pair< RootInfo , RootInfo >( ri1 , ri2 ) );
3497 if( vertexCount.count( key1 )==0 )
3499 vertexCount[key1].first = ri1;
3500 vertexCount[key1].second = 0;
3502 if( vertexCount.count( key2 )==0 )
3504 vertexCount[key2].first = ri2;
3505 vertexCount[key2].second = 0;
3507 vertexCount[key1].second++;
3508 vertexCount[key2].second--;
3514 fprintf( stderr ,
"Bad Edge 2: %d %d\t%d %d\n" , ri1.key , ri2.key , r1 , r2 );
3518 for(
int i=0 ; i<int(edges.size()) ; i++ )
3520 if( vertexCount.count( edges[i].first.key )==0 ) printf(
"Could not find vertex: %lld\n" , edges[i].first );
3521 else if( vertexCount[ edges[i].first.key ].second )
3524 GetRootPair( vertexCount[edges[i].first.key].first , fData.depth , ri );
3525 long long key = ri.key;
3526 if( vertexCount.count( key )==0 )
3529 node->depthAndOffset( d , off );
3530 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] );
3534 edges.push_back( std::pair< RootInfo , RootInfo >( ri , edges[i].first ) );
3535 vertexCount[ key ].second++;
3536 vertexCount[ edges[i].first.key ].second--;
3540 if( vertexCount.count( edges[i].second.key )==0 ) printf(
"Could not find vertex: %lld\n" , edges[i].second );
3541 else if( vertexCount[edges[i].second.key].second )
3544 GetRootPair( vertexCount[edges[i].second.key].first , fData.depth , ri );
3545 long long key = ri.key;
3546 if( vertexCount.count( key ) )
3549 node->depthAndOffset( d , off );
3550 printf(
"Vertex pair not in list 2\t[%d] (%d %d %d)\n" , d , off[0] , off[1] , off[2] );
3554 edges.push_back( std::pair< RootInfo , RootInfo >( edges[i].second , ri ) );
3555 vertexCount[key].second--;
3556 vertexCount[ edges[i].second.key ].second++;
3561 template<
int Degree>
3563 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 )
3564 #else // !MISHA_DEBUG
3565 int Octree< Degree >::GetMCIsoTriangles(
TreeOctNode* node , CoredMeshData* mesh , RootData& rootData , std::vector<
Point3D< float > >* interiorPositions ,
int offSet ,
int sDepth ,
bool addBarycenter ,
bool polygonMesh )
3566 #endif // MISHA_DEBUG
3569 std::vector< std::pair< RootInfo , RootInfo > > edges;
3570 std::vector< std::vector< std::pair< RootInfo , RootInfo > > > edgeLoops;
3571 GetMCIsoEdges( node , sDepth , edges );
3573 GetEdgeLoops( edges , edgeLoops );
3574 for(
int i=0 ; i<int(edgeLoops.size()) ; i++ )
3577 std::vector<CoredPointIndex> edgeIndices;
3578 for(
int j=0 ; j<int(edgeLoops[i].size()) ; j++ )
3580 if( !GetRootIndex( edgeLoops[i][j].first , rootData , p ) ) printf(
"Bad Point Index\n" );
3581 else edgeIndices.push_back( p );
3584 tris += AddTriangles( mesh , edgeIndices , interiorPositions , offSet , polygonMesh , barycenters );
3585 #else // !MISHA_DEBUG
3586 tris += AddTriangles( mesh , edgeIndices , interiorPositions , offSet , addBarycenter , polygonMesh );
3587 #endif // MISHA_DEBUG
3592 template<
int Degree >
3593 int Octree< Degree >::GetEdgeLoops( std::vector< std::pair< RootInfo , RootInfo > >& edges , std::vector< std::vector< std::pair< RootInfo , RootInfo > > >& loops )
3596 long long frontIdx , backIdx;
3597 std::pair< RootInfo , RootInfo > e , temp;
3600 while( edges.size() )
3602 std::vector< std::pair< RootInfo , RootInfo > > front , back;
3604 loops.resize( loopSize+1 );
3605 edges[0] = edges.back();
3607 frontIdx = e.second.key;
3608 backIdx = e.first.key;
3609 for(
int j=
int(edges.size())-1 ; j>=0 ; j-- )
3611 if( edges[j].first.key==frontIdx || edges[j].second.key==frontIdx )
3613 if( edges[j].first.key==frontIdx ) temp = edges[j];
3614 else temp.first = edges[j].second , temp.second = edges[j].first;
3615 frontIdx = temp.second.key;
3616 front.push_back(temp);
3617 edges[j] = edges.back();
3619 j = int(edges.size());
3621 else if( edges[j].first.key==backIdx || edges[j].second.key==backIdx )
3623 if( edges[j].second.key==backIdx ) temp = edges[j];
3624 else temp.first = edges[j].second , temp.second = edges[j].first;
3625 backIdx = temp.first.key;
3626 back.push_back(temp);
3627 edges[j] = edges.back();
3629 j = int(edges.size());
3632 for(
int j=
int(back.size())-1 ; j>=0 ; j-- ) loops[loopSize].push_back( back[j] );
3633 loops[loopSize].push_back(e);
3634 for(
int j=0 ; j<int(front.size()) ; j++ ) loops[loopSize].push_back( front[j] );
3637 return int(loops.size());
3639 template<
int Degree>
3641 int Octree<Degree>::AddTriangles( CoredMeshData* mesh , std::vector<CoredPointIndex>& edges , std::vector<Point3D<float> >* interiorPositions ,
int offSet ,
bool polygonMesh , std::vector<
Point3D< float > >* barycenters )
3642 #else // !MISHA_DEBUG
3643 int Octree<Degree>::AddTriangles( CoredMeshData* mesh , std::vector<CoredPointIndex>& edges , std::vector<Point3D<float> >* interiorPositions ,
int offSet ,
bool addBarycenter ,
bool polygonMesh )
3644 #endif // MISHA_DEBUG
3646 MinimalAreaTriangulation< float > MAT;
3647 std::vector< Point3D< float > > vertices;
3648 std::vector< TriangleIndex > triangles;
3651 std::vector< CoredVertexIndex > vertices( edges.size() );
3652 for(
int i=0 ; i<int(edges.size()) ; i++ )
3654 vertices[i].idx = edges[i].index;
3655 vertices[i].inCore = (edges[i].inCore!=0);
3657 #pragma omp critical (add_polygon_access)
3659 mesh->addPolygon( vertices );
3663 if( edges.size()>3 )
3665 bool isCoplanar =
false;
3669 #else // !MISHA_DEBUG
3671 #endif // MISHA_DEBUG
3672 for(
int i=0 ; i<int(edges.size()) ; i++ )
3673 for(
int j=0 ; j<i ; j++ )
3674 if( (i+1)%edges.size()!=j && (j+1)%edges.size()!=i )
3676 Point3D< Real > v1 , v2;
3677 if( edges[i].inCore ) v1 = mesh->inCorePoints[ edges[i].index ];
3678 else v1 = (*interiorPositions)[ edges[i].index-offSet ];
3679 if( edges[j].inCore ) v2 = mesh->inCorePoints[ edges[j].index ];
3680 else v2 = (*interiorPositions)[ edges[j].index-offSet ];
3681 for(
int k=0 ; k<3 ; k++ )
if( v1[k]==v2[k] ) isCoplanar =
true;
3686 c[0] = c[1] = c[2] = 0;
3687 for(
int i=0 ; i<int(edges.size()) ; i++ )
3690 if(edges[i].inCore) p = mesh->inCorePoints[edges[i].index ];
3691 else p = (*interiorPositions)[edges[i].index-offSet];
3694 c /=
Real( edges.size() );
3696 #pragma omp critical (add_point_access)
3698 cIdx = mesh->addOutOfCorePoint( c );
3700 barycenters->push_back( c );
3701 #else // !MISHA_DEBUG
3702 interiorPositions->push_back( c );
3703 #endif // MISHA_DEBUG
3705 for(
int i=0 ; i<int(edges.size()) ; i++ )
3707 std::vector< CoredVertexIndex > vertices( 3 );
3708 vertices[0].idx = edges[i ].index;
3709 vertices[1].idx = edges[(i+1)%edges.size()].index;
3710 vertices[2].idx = cIdx;
3711 vertices[0].inCore = (edges[i ].inCore!=0);
3712 vertices[1].inCore = (edges[(i+1)%edges.size()].inCore!=0);
3713 vertices[2].inCore = 0;
3714 #pragma omp critical (add_polygon_access)
3716 mesh->addPolygon( vertices );
3719 return int( edges.size() );
3723 vertices.resize( edges.size() );
3725 for(
int i=0 ; i<int(edges.size()) ; i++ )
3728 if( edges[i].inCore ) p = mesh->inCorePoints[edges[i].index ];
3729 else p = (*interiorPositions)[edges[i].index-offSet];
3732 MAT.GetTriangulation( vertices , triangles );
3733 for(
int i=0 ; i<int(triangles.size()) ; i++ )
3735 std::vector< CoredVertexIndex > _vertices( 3 );
3736 for(
int j=0 ; j<3 ; j++ )
3738 _vertices[j].idx = edges[ triangles[i].idx[j] ].index;
3739 _vertices[j].inCore = (edges[ triangles[i].idx[j] ].inCore!=0);
3741 #pragma omp critical (add_polygon_access)
3743 mesh->addPolygon( _vertices );
3748 else if( edges.size()==3 )
3750 std::vector< CoredVertexIndex > vertices( 3 );
3751 for(
int i=0 ; i<3 ; i++ )
3753 vertices[i].idx = edges[i].index;
3754 vertices[i].inCore = (edges[i].inCore!=0);
3756 #pragma omp critical (add_polygon_access)
3757 mesh->addPolygon( vertices );
3759 return int(edges.size())-2;
3761 template<
int Degree >
3764 if( depth<=0 || depth>tree.maxDepth() ) depth = tree.maxDepth();
3769 Real* values =
new float[ res * res * res ];
3770 memset( values , 0 ,
sizeof(
float ) * res * res * res );
3774 if( n->d>depth )
continue;
3775 if( n->d<_minDepth )
continue;
3776 int d , idx[3] , start[3] , end[3];
3777 n->depthAndOffset( d , idx );
3778 for(
int i=0 ; i<3 ; i++ )
3785 if( !(start[i]&1) ) start[i]++;
3786 if( !( end[i]&1) ) end[i]--;
3788 Real coefficient = n->nodeData.solution;
3789 for(
int x=start[0] ; x<=end[0] ; x+=2 )
3790 for(
int y=start[1] ; y<=end[1] ; y+=2 )
3791 for(
int z=start[2] ; z<=end[2] ; z+=2 )
3793 int xx = (x-1)>>1 , yy = (y-1)>>1 , zz = (z-1)>>1;
3794 values[ zz*res*res + yy*res + xx ] +=
3801 for(
int i=0 ; i<res*res*res ; i++ ) values[i] -= isoValue;
3805 template<
int Degree >
3808 if( depth<=0 || depth>tree.maxDepth() ) depth = tree.maxDepth();
3809 res = 1<<tree.maxDepth();
3810 Real* values =
new float[ res * res * res ];
3811 memset( values , 0 ,
sizeof(
float ) * res * res * res );
3815 if( n->d>depth )
continue;
3816 int d , idx[3] , start[3] , end[3];
3817 n->depthAndOffset( d , idx );
3818 for(
int i=0 ; i<3 ; i++ )
3823 fData.setSampleSpan( idx[i] , start[i] , end[i] );
3825 if( !(start[i]&1) ) start[i]++;
3826 if( !( end[i]&1) ) end[i]--;
3828 for(
int x=start[0] ; x<=end[0] ; x+=2 )
3829 for(
int y=start[1] ; y<=end[1] ; y+=2 )
3830 for(
int z=start[2] ; z<=end[2] ; z+=2 )
3832 int xx = (x-1)>>1 , yy = (y-1)>>1 , zz = (z-1)>>1;
3833 values[ zz*res*res + yy*res + xx ] +=
3834 n->nodeData.centerWeightContribution *
3835 fData.valueTables[ idx[0]+x*fData.functionCount ] *
3836 fData.valueTables[ idx[1]+y*fData.functionCount ] *
3837 fData.valueTables[ idx[2]+z*fData.functionCount ];
3854 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3858 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3870 for(
int i=0 ; i<DIMENSION ; i++ ) idx[i] = BinaryNode<Real>::CornerIndex( maxDepth+1 , d , o[i] , x[i] );
3877 for(
int i=0 ; i<DIMENSION ; i++ ) idx[i] = BinaryNode<Real>::CornerIndex( maxDepth+1 , depth , offSet[i] , x[i] );
3882 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3886 return FaceIndex(node,fIndex,maxDepth,idx);
3895 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3899 return EdgeIndex(node,eIndex,maxDepth,idx);
3921 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;