29 #include "octree_poisson.h"
32 #if defined WIN32 || defined _WIN32
33 #if !defined __MINGW32__
40 #define ITERATION_POWER 1.0/3
41 #define MEMORY_ALLOCATOR_BLOCK_SIZE 1<<12
47 using namespace __gnu_cxx;
60 #if defined _WIN32 && !defined __MINGW32__
61 using stdext::hash_map;
69 #if defined _WIN32 && !defined __MINGW32__
71 _InterlockedOr( (
long volatile*)&dest, value );
73 InterlockedOr( (
long volatile*)&dest , value );
75 #else // !_WIN32 || __MINGW32__
78 #endif // _WIN32 && !__MINGW32__
85 SortedTreeNodes::SortedTreeNodes(
void)
91 SortedTreeNodes::~SortedTreeNodes(
void){
92 if( nodeCount )
delete[] nodeCount;
93 if( treeNodes )
delete[] treeNodes;
100 if( nodeCount )
delete[] nodeCount;
101 if( treeNodes )
delete[] treeNodes;
103 nodeCount =
new int[ maxDepth+1 ];
106 nodeCount[0] = 0 , nodeCount[1] = 1;
107 treeNodes[0] = &root;
108 for(
int d=1 ; d<maxDepth ; d++ )
110 nodeCount[d+1] = nodeCount[d];
111 for(
int i=nodeCount[d-1] ; i<nodeCount[d] ; i++ )
114 if( temp->
children )
for(
int c=0 ; c<8 ; c++ ) treeNodes[ nodeCount[d+1]++ ] = temp->
children + c;
117 for(
int i=0 ; i<nodeCount[maxDepth] ; i++ ) treeNodes[i]->nodeData.nodeIndex = i;
125 if( threads<=0 ) threads = 1;
127 std::vector< std::pair< int , int > > spans( this->maxDepth , std::pair< int , int >( -1 , -1 ) );
128 int minDepth , off[3];
130 cData.
offsets.resize( this->maxDepth , -1 );
134 for(
int d=minDepth ; d<=maxDepth ; d++ )
136 spans[d] = std::pair< int , int >( start , end+1 );
137 cData.
offsets[d] = nodeCount - spans[d].first;
138 nodeCount += spans[d].second - spans[d].first;
141 while( start< end && !treeNodes[start]->children ) start++;
142 while( end> start && !treeNodes[end ]->children ) end--;
143 if( start==end && !treeNodes[start]->children )
break;
144 start = treeNodes[start]->children[0].nodeData.nodeIndex;
145 end = treeNodes[end ]->children[7].nodeData.nodeIndex;
149 cData.
cTable.resize( nodeCount );
150 std::vector< int > count( threads );
151 #pragma omp parallel for num_threads( threads )
152 for(
int t=0 ; t<threads ; t++ )
154 TreeOctNode::ConstNeighborKey3 neighborKey;
155 neighborKey.set( maxDepth );
156 int offset = nodeCount * t * Cube::CORNERS;
158 for(
int d=minDepth ; d<=maxDepth ; d++ )
160 int start = spans[d].first , end = spans[d].second , width = end-start;
161 for(
int i=start + (width*t)/threads ; i<start + (width*(t+1))/threads ; i++ )
164 if( d<maxDepth && node->children )
continue;
165 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node , minDepth );
166 for(
int c=0 ; c<Cube::CORNERS ; c++ )
168 bool cornerOwner =
true;
170 int ac = Cube::AntipodalCornerIndex( c );
171 Cube::FactorCornerIndex( c , x , y , z );
172 for(
int cc=0 ; cc<Cube::CORNERS ; cc++ )
175 Cube::FactorCornerIndex( cc , xx , yy , zz );
176 xx += x , yy += y , zz += z;
177 if( neighbors.neighbors[xx][yy][zz] )
178 if( cc<ac || ( d<maxDepth && neighbors.neighbors[xx][yy][zz]->children ) )
181 neighbors.neighbors[xx][yy][zz]->depthAndOffset( _d , _off );
182 _off[0] >>= (d-minDepth) , _off[1] >>= (d-minDepth) , _off[2] >>= (d-minDepth);
183 if( _off[0]==off[0] && _off[1]==off[1] && _off[2]==off[2] )
188 else fprintf( stderr ,
"[WARNING] How did we leave the subtree?\n" );
197 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.neighbors[d];
199 for(
int cc=0 ; cc<Cube::CORNERS ; cc++ )
202 Cube::FactorCornerIndex( cc , xx , yy , zz );
203 xx += x , yy += y , zz += z;
204 if( neighborKey.neighbors[d].neighbors[xx][yy][zz] )
205 cData[ neighbors.neighbors[xx][yy][zz] ][ Cube::AntipodalCornerIndex(cc) ] = count[t] + offset;
220 std::vector< int > offsets( threads+1 );
222 for(
int t=0 ; t<threads ; t++ ) cData.
cCount += count[t] , offsets[t+1] = offsets[t] + count[t];
223 #pragma omp parallel
for num_threads( threads )
224 for(
int t=0 ; t<threads ; t++ )
225 for(
int d=minDepth ; d<=maxDepth ; d++ )
227 int start = spans[d].first , end = spans[d].second , width = end - start;
228 for(
int i=start + (width*t)/threads ; i<start+(width*(t+1))/threads ; i++ )
229 for(
int c=0 ; c<Cube::CORNERS ; c++ )
231 int& idx = cData[ treeNodes[i] ][c];
234 fprintf( stderr ,
"[ERROR] Found unindexed corner nodes[%d][%d] = %d (%d,%d)\n" , treeNodes[i]->nodeData.nodeIndex , c , idx , minDepth , maxDepth );
236 treeNodes[i]->depthAndOffset( _d , _off );
237 printf(
"(%d [%d %d %d) <-> (%d [%d %d %d])\n" , minDepth , off[0] , off[1] , off[2] , _d , _off[0] , _off[1] , _off[2] );
238 printf(
"[%d %d]\n" , spans[d].first , spans[d].second );
243 int div = idx / ( nodeCount*Cube::CORNERS );
244 int rem = idx % ( nodeCount*Cube::CORNERS );
245 idx = rem + offsets[div];
250 int SortedTreeNodes::getMaxCornerCount(
const TreeOctNode* rootNode ,
int depth ,
int maxDepth ,
int threads )
const
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;
262 neighborKey.set( maxDepth );
263 int start = nodeCount[depth] , end = nodeCount[maxDepth+1] , range = end-start;
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 );
272 for(
int c=0 ; c<Cube::CORNERS ; c++ )
274 bool cornerOwner =
true;
276 int ac = Cube::AntipodalCornerIndex( c );
277 Cube::FactorCornerIndex( c , x , y , z );
278 for(
int cc=0 ; cc<Cube::CORNERS ; cc++ )
281 Cube::FactorCornerIndex( cc , xx , yy , zz );
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 );
321 eData.
offsets[d] = nodeCount - spans[d].first;
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;
328 start = treeNodes[start]->children[0].nodeData.nodeIndex;
329 end = treeNodes[end ]->children[7].nodeData.nodeIndex;
333 eData.
eTable.resize( nodeCount );
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;
339 neighborKey.set( maxDepth );
340 int offset = nodeCount * t * Cube::EDGES;
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 );
350 for(
int e=0 ; e<Cube::EDGES ; e++ )
352 bool edgeOwner =
true;
354 Cube::FactorEdgeIndex( e , o , i , j );
355 int ac = Square::AntipodalCornerIndex( Square::CornerIndex( i , j ) );
356 for(
int cc=0 ; cc<Square::CORNERS ; cc++ )
358 int ii , jj , x , y , z;
359 Square::FactorCornerIndex( cc , ii , jj );
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; }
372 for(
int cc=0 ; cc<Square::CORNERS ; cc++ )
374 int ii , jj , aii , ajj , x , y , z;
375 Square::FactorCornerIndex( cc , ii , jj );
376 Square::FactorCornerIndex( Square::AntipodalCornerIndex( cc ) , aii , ajj );
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++ ) eData.
eCount += count[t] , offsets[t+1] = offsets[t] + count[t];
397 #pragma omp parallel
for num_threads( threads )
398 for(
int t=0 ; t<threads ; t++ )
399 for(
int d=minDepth ; d<=maxDepth ; d++ )
401 int start = spans[d].first , end = spans[d].second , width = end - start;
402 for(
int i=start + (width*t)/threads ; i<start+(width*(t+1))/threads ; i++ )
403 for(
int e=0 ; e<Cube::EDGES ; e++ )
405 int& idx = eData[ treeNodes[i] ][e];
406 if( idx<0 ) fprintf( stderr ,
"[ERROR] Found unindexed edge %d (%d,%d)\n" , idx , minDepth , maxDepth ) , exit( 0 );
409 int div = idx / ( nodeCount*Cube::EDGES );
410 int rem = idx % ( nodeCount*Cube::EDGES );
411 idx = rem + offsets[div];
416 int SortedTreeNodes::getMaxEdgeCount(
const TreeOctNode* rootNode ,
int depth ,
int threads )
const
418 if( threads<=0 ) threads = 1;
420 std::vector< std::vector< int > > edgeCount( threads );
421 for(
int t=0 ; t<threads ; t++ ) edgeCount[t].resize( res*res*res , 0 );
423 #pragma omp parallel for num_threads( threads )
424 for(
int t=0 ; t<threads ; t++ )
426 std::vector< int >& _edgeCount = edgeCount[t];
427 TreeOctNode::ConstNeighborKey3 neighborKey;
428 neighborKey.set( maxDepth-1 );
429 int start = nodeCount[depth] , end = nodeCount[maxDepth] , range = end-start;
430 for(
int i=(range*t)/threads ; i<(range*(t+1))/threads ; i++ )
433 const TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node , depth );
437 for(
int e=0 ; e<Cube::EDGES ; e++ )
439 bool edgeOwner =
true;
441 Cube::FactorEdgeIndex( e , o , i , j );
442 int ac = Square::AntipodalCornerIndex( Square::CornerIndex( i , j ) );
443 for(
int cc=0 ; cc<Square::CORNERS ; cc++ )
445 int ii , jj , x , y , z;
446 Square::FactorCornerIndex( cc , ii , jj );
450 case 0: y = ii , z = jj , x = 1 ;
break;
451 case 1: x = ii , z = jj , y = 1 ;
break;
452 case 2: x = ii , y = jj , z = 1 ;
break;
454 if( neighbors.neighbors[x][y][z] && cc<ac ) { edgeOwner = false ;
break; }
456 if( edgeOwner ) _edgeCount[ ( ( off[0]>>(d-depth) ) * res * res) + ( ( off[1]>>(d-depth) ) * res) + ( off[2]>>(d-depth) ) ]++;
461 for(
int i=0 ; i<res*res*res ; i++ )
464 for(
int t=0 ; t<threads ; t++ ) c += edgeCount[t][i];
465 maxCount = std::max< int >( maxCount , c );
475 int TreeNodeData::UseIndex=1;
476 TreeNodeData::TreeNodeData(
void )
481 centerWeightContribution=0;
485 constraint = solution = 0;
488 TreeNodeData::~TreeNodeData(
void ) { }
502 if(mem>maxMemoryUsage){maxMemoryUsage=mem;}
512 postNormalSmooth = 0;
513 _constrainValues =
false;
519 double x , dxdy , dxdydz , dx[DIMENSION][SPLAT_ORDER+1];
521 TreeOctNode::Neighbors3& neighbors = neighborKey.setNeighbors( node );
527 for(
int i=0 ; i<3 ; i++ )
531 x = ( center[i] - position[i] - width ) / width;
532 dx[i][0] = 1.125+1.500*x+0.500*x*x;
533 x = ( center[i] - position[i] ) / width;
534 dx[i][1] = 0.750 - x*x;
536 dx[i][2] = 1. - dx[i][1] - dx[i][0];
538 x = ( position[i] - center[i] ) / width;
549 dx[i][1] = 1. - dx[i][0];
554 # error Splat order not supported
555 #endif // SPLAT_ORDER
557 for(
int i=off[0] ; i<=off[0]+SPLAT_ORDER ; i++ )
for(
int j=off[1] ; j<=off[1]+SPLAT_ORDER ; j++ )
559 dxdy = dx[0][i] * dx[1][j];
560 for(
int k=off[2] ; k<=off[2]+SPLAT_ORDER ; k++ )
561 if( neighbors.neighbors[i][j][k] )
563 dxdydz = dxdy * dx[2][k];
569 n[0] = n[1] = n[2] = 0;
570 _node->nodeData.nodeIndex = 0;
571 idx = _node->nodeData.normalIndex = int(normals->size());
572 normals->push_back(n);
574 (*normals)[idx] += normal *
Real( dxdydz );
580 Real Octree<Degree>::NonLinearSplatOrientedPoint(
const Point3D<Real>& position ,
const Point3D<Real>& normal ,
int splatDepth , Real samplesPerNode ,
581 int minDepth ,
int maxDepth )
588 Point3D< Real > myCenter;
590 myCenter[0] = myCenter[1] = myCenter[2] =
Real(0.5);
594 while( temp->depth()<splatDepth )
596 if( !temp->children )
598 fprintf( stderr ,
"Octree<Degree>::NonLinearSplatOrientedPoint error\n" );
601 int cIndex=TreeOctNode::CornerIndex(myCenter,position);
602 temp=&temp->children[cIndex];
604 if(cIndex&1) myCenter[0] += myWidth/2;
605 else myCenter[0] -= myWidth/2;
606 if(cIndex&2) myCenter[1] += myWidth/2;
607 else myCenter[1] -= myWidth/2;
608 if(cIndex&4) myCenter[2] += myWidth/2;
609 else myCenter[2] -= myWidth/2;
612 NonLinearGetSampleDepthAndWeight( temp , position , samplesPerNode , newDepth , alpha );
614 if( newDepth<minDepth ) newDepth=
Real(minDepth);
615 if( newDepth>maxDepth ) newDepth=
Real(maxDepth);
616 int topDepth=int(ceil(newDepth));
618 dx = 1.0-(topDepth-newDepth);
619 if( topDepth<=minDepth )
624 else if( topDepth>maxDepth )
629 while( temp->depth()>topDepth ) temp=temp->parent;
630 while( temp->depth()<topDepth )
632 if(!temp->children) temp->initChildren();
633 int cIndex=TreeOctNode::CornerIndex(myCenter,position);
634 temp=&temp->children[cIndex];
636 if(cIndex&1) myCenter[0] += myWidth/2;
637 else myCenter[0] -= myWidth/2;
638 if(cIndex&2) myCenter[1] += myWidth/2;
639 else myCenter[1] -= myWidth/2;
640 if(cIndex&4) myCenter[2] += myWidth/2;
641 else myCenter[2] -= myWidth/2;
643 width = 1.0 / ( 1<<temp->depth() );
644 n = normal * alpha /
Real( pow( width , 3 ) ) *
Real( dx );
645 NonLinearSplatOrientedPoint( temp , position , n );
646 if( fabs(1.0-dx) > EPSILON )
650 width = 1.0 / ( 1<<temp->depth() );
652 n = normal * alpha /
Real( pow( width , 3 ) ) *
Real( dx );
653 NonLinearSplatOrientedPoint( temp , position , n );
658 void Octree<Degree>::NonLinearGetSampleDepthAndWeight(TreeOctNode* node,
const Point3D<Real>& position,Real samplesPerNode,Real& depth,Real& weight){
660 weight =
Real(1.0)/NonLinearGetSampleWeight(temp,position);
661 if( weight>=samplesPerNode ) depth=
Real( temp->depth() + log( weight / samplesPerNode ) / log(
double(1<<(DIMENSION-1))) );
664 Real oldAlpha,newAlpha;
665 oldAlpha=newAlpha=weight;
666 while( newAlpha<samplesPerNode && temp->parent )
670 newAlpha=
Real(1.0)/NonLinearGetSampleWeight(temp,position);
672 depth =
Real( temp->depth() + log( newAlpha / samplesPerNode ) / log( newAlpha / oldAlpha ) );
674 weight=
Real(pow(
double(1<<(DIMENSION-1)),-
double(depth)));
678 Real Octree<Degree>::NonLinearGetSampleWeight( TreeOctNode* node ,
const Point3D<Real>& position )
681 double x,dxdy,dx[DIMENSION][3];
682 TreeOctNode::Neighbors3& neighbors=neighborKey.setNeighbors( node );
684 Point3D<Real> center;
686 node->centerAndWidth(center,w);
689 for(
int i=0 ; i<DIMENSION ; i++ )
691 x = ( center[i] - position[i] - width ) / width;
692 dx[i][0] = 1.125 + 1.500*x + 0.500*x*x;
693 x = ( center[i] - position[i] ) / width;
694 dx[i][1] = 0.750 - x*x;
696 dx[i][2] = 1.0 - dx[i][1] - dx[i][0];
699 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
701 dxdy = dx[0][i] * dx[1][j];
702 for(
int k=0 ; k<3 ; k++ )
if( neighbors.neighbors[i][j][k] )
703 weight +=
Real( dxdy * dx[2][k] * neighbors.neighbors[i][j][k]->nodeData.centerWeightContribution );
705 return Real( 1.0 / weight );
709 int Octree<Degree>::NonLinearUpdateWeightContribution( TreeOctNode* node ,
const Point3D<Real>& position , Real weight )
711 TreeOctNode::Neighbors3& neighbors = neighborKey.setNeighbors( node );
712 double x,dxdy,dx[DIMENSION][3];
714 Point3D<Real> center;
716 node->centerAndWidth( center , w );
718 const double SAMPLE_SCALE = 1. / ( 0.125 * 0.125 + 0.75 * 0.75 + 0.125 * 0.125 );
720 for(
int i=0 ; i<DIMENSION ; i++ )
722 x = ( center[i] - position[i] - width ) / width;
723 dx[i][0] = 1.125 + 1.500*x + 0.500*x*x;
724 x = ( center[i] - position[i] ) / width;
725 dx[i][1] = 0.750 - x*x;
726 dx[i][2] = 1. - dx[i][1] - dx[i][0];
729 dx[i][0] *= SAMPLE_SCALE;
732 for(
int i=0 ; i<3 ; i++ )
for(
int j=0 ; j<3 ; j++ )
734 dxdy = dx[0][i] * dx[1][j] * weight;
735 for(
int k=0 ; k<3 ; k++ )
if( neighbors.neighbors[i][j][k] )
736 neighbors.neighbors[i][j][k]->nodeData.centerWeightContribution +=
Real( dxdy * dx[2][k] );
741 template<
int Degree >
template<
typename Po
intNT>
int
744 int useConfidence ,
Real constraintWeight ,
bool adaptiveWeights )
746 _minDepth = std::min< int >( std::max< int >( 0 , minDepth ) , maxDepth );
747 _constrainValues = (constraintWeight>0);
748 double pointWeightSum = 0;
755 TreeNodeData::UseIndex = 1;
756 neighborKey.set( maxDepth );
757 splatDepth = kernelDepth;
758 if( splatDepth<0 ) splatDepth = 0;
763 while (cnt != input_->size ())
766 p[0] = input_->points[cnt].x;
767 p[1] = input_->points[cnt].y;
768 p[2] = input_->points[cnt].z;
770 for (i = 0; i < DIMENSION; i++)
772 if( !cnt || p[i]<min[i] ) min[i] = p[i];
773 if( !cnt || p[i]>max[i] ) max[i] = p[i];
778 scale = std::max< Real >( max[0]-min[0] , std::max< Real >( max[1]-min[1] , max[2]-min[2] ) );
779 center = ( max+min ) /2;
781 scale *= scaleFactor;
782 for( i=0 ; i<DIMENSION ; i++ ) center[i] -= scale/2;
786 while (cnt != input_->size ())
788 position[0] = input_->points[cnt].x;
789 position[1] = input_->points[cnt].y;
790 position[2] = input_->points[cnt].z;
791 normal[0] = input_->points[cnt].normal_x;
792 normal[1] = input_->points[cnt].normal_y;
793 normal[2] = input_->points[cnt].normal_z;
795 for( i=0 ; i<DIMENSION ; i++ ) position[i] = ( position[i]-center[i] ) / scale;
796 myCenter[0] = myCenter[1] = myCenter[2] =
Real(0.5);
798 for( i=0 ; i<DIMENSION ; i++ )
if( position[i]<myCenter[i]-myWidth/2 || position[i]>myCenter[i]+myWidth/2 )
break;
799 if( i!=DIMENSION )
continue;
801 if( useConfidence ) weight =
Real(
Length(normal) );
804 while( d<splatDepth )
806 NonLinearUpdateWeightContribution( temp , position , weight );
808 int cIndex=TreeOctNode::CornerIndex(myCenter,position);
811 if(cIndex&1) myCenter[0] += myWidth/2;
812 else myCenter[0] -= myWidth/2;
813 if(cIndex&2) myCenter[1] += myWidth/2;
814 else myCenter[1] -= myWidth/2;
815 if(cIndex&4) myCenter[2] += myWidth/2;
816 else myCenter[2] -= myWidth/2;
819 NonLinearUpdateWeightContribution( temp , position , weight );
824 normals =
new std::vector< Point3D<Real> >();
826 while (cnt != input_->size ())
828 position[0] = input_->points[cnt].x;
829 position[1] = input_->points[cnt].y;
830 position[2] = input_->points[cnt].z;
831 normal[0] = input_->points[cnt].normal_x;
832 normal[1] = input_->points[cnt].normal_y;
833 normal[2] = input_->points[cnt].normal_z;
835 for( i=0 ; i<DIMENSION ; i++ ) position[i] = ( position[i]-center[i] ) / scale;
836 myCenter[0] = myCenter[1] = myCenter[2] =
Real(0.5);
838 for( i=0 ; i<DIMENSION ; i++ )
if(position[i]<myCenter[i]-myWidth/2 || position[i]>myCenter[i]+myWidth/2)
break;
839 if( i!=DIMENSION )
continue;
841 if( l!=l || l<=
EPSILON )
continue;
842 if( !useConfidence ) normal /= l;
846 if( samplesPerNode>0 && splatDepth )
848 pointWeight = NonLinearSplatOrientedPoint( position , normal , splatDepth , samplesPerNode , _minDepth , maxDepth );
857 while( d<splatDepth )
859 int cIndex=TreeOctNode::CornerIndex(myCenter,position);
862 if(cIndex&1) myCenter[0]+=myWidth/2;
863 else myCenter[0]-=myWidth/2;
864 if(cIndex&2) myCenter[1]+=myWidth/2;
865 else myCenter[1]-=myWidth/2;
866 if(cIndex&4) myCenter[2]+=myWidth/2;
867 else myCenter[2]-=myWidth/2;
870 alpha = NonLinearGetSampleWeight( temp , position );
872 for( i=0 ; i<DIMENSION ; i++ ) normal[i]*=alpha;
876 int cIndex=TreeOctNode::CornerIndex(myCenter,position);
879 if(cIndex&1) myCenter[0]+=myWidth/2;
880 else myCenter[0]-=myWidth/2;
881 if(cIndex&2) myCenter[1]+=myWidth/2;
882 else myCenter[1]-=myWidth/2;
883 if(cIndex&4) myCenter[2]+=myWidth/2;
884 else myCenter[2]-=myWidth/2;
887 NonLinearSplatOrientedPoint( temp , position , normal );
891 pointWeightSum += pointWeight;
892 if( _constrainValues )
896 myCenter[0] = myCenter[1] = myCenter[2] =
Real(0.5);
904 p[0] = p[1] = p[2] = 0;
905 idx = int( _points.size() );
906 _points.push_back( PointData( position*pointWeight , pointWeight ) );
911 _points[idx].weight += pointWeight;
912 _points[idx].position += position * pointWeight;
915 int cIndex = TreeOctNode::CornerIndex( myCenter , position );
919 if( cIndex&1 ) myCenter[0] += myWidth/2;
920 else myCenter[0] -= myWidth/2;
921 if( cIndex&2 ) myCenter[1] += myWidth/2;
922 else myCenter[1] -= myWidth/2;
923 if( cIndex&4 ) myCenter[2] += myWidth/2;
924 else myCenter[2] -= myWidth/2;
931 if( _constrainValues )
933 if( n->nodeData.pointIndex!=-1 )
935 int idx = n->nodeData.pointIndex;
936 _points[idx].position /= _points[idx].weight;
937 if( adaptiveWeights ) _points[idx].weight *= (1<<n->d);
938 else _points[idx].weight *= (1<<maxDepth);
939 _points[idx].weight *=
Real( constraintWeight / pointWeightSum );
941 #if FORCE_NEUMANN_FIELD
944 int d , off[3] , res;
945 node->depthAndOffset( d , off );
947 if( node->nodeData.normalIndex<0 )
continue;
949 for(
int d=0 ; d<3 ; d++ )
if( off[d]==0 || off[d]==res-1 ) normal[d] = 0;
951 #endif // FORCE_NEUMANN_FIELD
962 radius = 0.5 + 0.5 * Degree;
963 width=int(
double(radius+0.5-
EPSILON)*2);
964 if( normalSmooth>0 ) postNormalSmooth = normalSmooth;
965 fData.set( maxDepth ,
true , reflectBoundary );
971 int maxDepth = tree.maxDepth( );
972 TreeOctNode::NeighborKey5 nKey;
973 nKey.set( maxDepth );
974 for(
int d=maxDepth ; d>0 ; d-- )
978 int xStart=0 , xEnd=5 , yStart=0 , yEnd=5 , zStart=0 , zEnd=5;
979 int c = int( node - node->parent->children );
981 Cube::FactorCornerIndex( c , x , y , z );
988 nKey.setNeighbors( node->parent , xStart , xEnd , yStart , yEnd , zStart , zEnd );
993 template<
int Degree >
997 const int min[] = { std::max<int>( 0 , d[0]+0 ) , std::max<int>( 0 , d[1]+0 ) , std::max<int>( 0 , d[2]+0 ) };
998 const int max[] = { std::min<int>( 2 , d[0]+2 ) , std::min<int>( 2 , d[1]+2 ) , std::min<int>( 2 , d[2]+2 ) };
999 for(
int i=min[0] ; i<=max[0] ; i++ )
for(
int j=min[1] ; j<=max[1] ; j++ )
1001 if( !hasPoints[i][j] )
continue;
1002 const PointInfo* pInfo = points[i][j];
1005 for(
int k=min[2] ; k<=max[2] ; k++ )
1006 if( pInfo[k].weightedValue )
1007 v += pInfo[k].splineValues[0][ii] * pInfo[k].splineValues[1][jj] * pInfo[k].splineValues[2][-d[2]+k];
1011 template<
int Degree>
1012 Real Octree<Degree>::GetLaplacian(
const int idx[DIMENSION] )
const
1014 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]] ) );
1016 template<
int Degree >
1017 Real Octree< Degree >::GetLaplacian(
const TreeOctNode* node1 ,
const TreeOctNode* node2 )
const
1021 BSplineData< Degree , Real >::SymmetricIndex( node1->off[0] , node2->off[0] ) ,
1022 BSplineData< Degree , Real >::SymmetricIndex( node1->off[1] , node2->off[1] ) ,
1023 BSplineData< Degree , Real >::SymmetricIndex( node1->off[2] , node2->off[2] )
1025 return GetLaplacian( symIndex );
1027 template<
int Degree >
1028 Real Octree< Degree >::GetDivergence(
const TreeOctNode* node1 ,
const TreeOctNode* node2 ,
const Point3D< Real >& normal1 )
const
1032 BSplineData< Degree , Real >::SymmetricIndex( node1->off[0] , node2->off[0] ) ,
1033 BSplineData< Degree , Real >::SymmetricIndex( node1->off[1] , node2->off[1] ) ,
1034 BSplineData< Degree , Real >::SymmetricIndex( node1->off[2] , node2->off[2] ) ,
1038 #if GRADIENT_DOMAIN_SOLUTION
1040 fData.Index( node2->off[0] , node1->off[0] ) ,
1041 fData.Index( node2->off[1] , node1->off[1] ) ,
1042 fData.Index( node2->off[2] , node1->off[2] )
1043 #else // !GRADIENT_DOMAIN_SOLUTION
1045 fData.Index( node1->off[0] , node2->off[0] ) ,
1046 fData.Index( node1->off[1] , node2->off[1] ) ,
1047 fData.Index( node1->off[2] , node2->off[2] )
1048 #endif // GRADIENT_DOMAIN_SOLUTION
1050 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1051 #if GRADIENT_DOMAIN_SOLUTION
1052 return Real( dot * ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] ) );
1053 #else // !GRADIENT_DOMAIN_SOLUTION
1054 return -
Real( dot * ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] ) );
1055 #endif // GRADIENT_DOMAIN_SOLUTION
1057 template<
int Degree >
1058 Real Octree< Degree >::GetDivergenceMinusLaplacian(
const TreeOctNode* node1 ,
const TreeOctNode* node2 , Real value1 ,
const Point3D<Real>& normal1 )
const
1062 BSplineData< Degree , Real >::SymmetricIndex( node1->off[0] , node2->off[0] ) ,
1063 BSplineData< Degree , Real >::SymmetricIndex( node1->off[1] , node2->off[1] ) ,
1064 BSplineData< Degree , Real >::SymmetricIndex( node1->off[2] , node2->off[2] )
1068 #if GRADIENT_DOMAIN_SOLUTION
1070 fData.Index( node2->off[0] , node1->off[0] ) ,
1071 fData.Index( node2->off[1] , node1->off[1] ) ,
1072 fData.Index( node2->off[2] , node1->off[2] )
1073 #else // !GRADIENT_DOMAIN_SOLUTION
1075 fData.Index( node1->off[0] , node2->off[0] ) ,
1076 fData.Index( node1->off[1] , node2->off[1] ) ,
1077 fData.Index( node1->off[2] , node2->off[2] )
1078 #endif // GRADIENT_DOMAIN_SOLUTION
1080 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1083 #
if GRADIENT_DOMAIN_SOLUTION
1084 - ( fData.ddDotTable[ symIndex[0]] + fData.ddDotTable[ symIndex[1]] + fData.ddDotTable[ symIndex[2]] ) * value1
1085 + ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] )
1087 - ( fData.ddDotTable[ symIndex[0]] + fData.ddDotTable[ symIndex[1]] + fData.ddDotTable[ symIndex[2]] ) * value1
1088 - ( fData.dvDotTable[aSymIndex[0]]*normal1[0] + fData.dvDotTable[aSymIndex[1]]*normal1[1] + fData.dvDotTable[aSymIndex[2]]*normal1[2] )
1093 template<
int Degree >
1094 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
1096 xStart = 0 , yStart = 0 , zStart = 0 , xEnd = 5 , yEnd = 5 , zEnd = 5;
1098 node->depthAndOffset( depth , off );
1099 int width = 1 << ( depth-rDepth );
1100 off[0] -= rOff[0] << ( depth-rDepth ) , off[1] -= rOff[1] << ( depth-rDepth ) , off[2] -= rOff[2] << ( depth-rDepth );
1101 if( off[0]<0 ) xStart = -off[0];
1102 if( off[1]<0 ) yStart = -off[1];
1103 if( off[2]<0 ) zStart = -off[2];
1104 if( off[0]>=width ) xEnd = 4 - ( off[0]-width );
1105 if( off[1]>=width ) yEnd = 4 - ( off[1]-width );
1106 if( off[2]>=width ) zEnd = 4 - ( off[2]-width );
1108 template<
int Degree >
1109 int Octree< Degree >::GetMatrixRowSize(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 )
const {
return GetMatrixRowSize( neighbors5 , 0 , 5 , 0 , 5 , 0 , 5 ); }
1110 template<
int Degree >
1111 int Octree< Degree >::GetMatrixRowSize(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 ,
int xStart ,
int xEnd ,
int yStart ,
int yEnd ,
int zStart ,
int zEnd )
const
1114 for(
int x=xStart ; x<=2 ; x++ )
1115 for(
int y=yStart ; y<yEnd ; y++ )
1116 if( x==2 && y>2 )
continue;
1117 else for(
int z=zStart ; z<zEnd ; z++ )
1118 if( x==2 && y==2 && z>2 )
continue;
1119 else if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1123 template<
int Degree >
1124 int Octree< Degree >::SetMatrixRow(
const OctNode< TreeNodeData , Real >::Neighbors5& neighbors5 , MatrixEntry< float >* row ,
int offset ,
const double stencil[5][5][5] )
const
1126 return SetMatrixRow( neighbors5 , row , offset , stencil , 0 , 5 , 0 , 5 , 0 , 5 );
1128 template<
int Degree >
1129 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
1131 bool hasPoints[3][3];
1133 PointInfo samples[3][3][3];
1136 const TreeOctNode* node = neighbors5.neighbors[2][2][2];
1137 int index[] = { int( node->off[0] ) , int( node->off[1] ), int( node->off[2] ) };
1139 if( _constrainValues )
1143 idx[0] = BinaryNode< double >::CenterIndex( d , idx[0] );
1144 idx[1] = BinaryNode< double >::CenterIndex( d , idx[1] );
1145 idx[2] = BinaryNode< double >::CenterIndex( d , idx[2] );
1146 for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
1148 hasPoints[j][k] =
false;
1149 for(
int l=0 ; l<3 ; l++ )
1151 const TreeOctNode* _node = neighbors5.neighbors[j+1][k+1][l+1];
1152 if( _node && _node->nodeData.pointIndex!=-1 )
1154 const PointData& pData = _points[ _node->nodeData.pointIndex ];
1155 PointInfo& pointInfo = samples[j][k][l];
1156 Real weight = pData.weight;
1157 Point3D< Real > p = pData.position;
1158 for(
int s=0 ; s<3 ; s++ )
1160 pointInfo.splineValues[0][s] = float( fData.baseBSplines[ idx[0]+j-s][s]( p[0] ) );
1161 pointInfo.splineValues[1][s] = float( fData.baseBSplines[ idx[1]+k-s][s]( p[1] ) );
1162 pointInfo.splineValues[2][s] = float( fData.baseBSplines[ idx[2]+l-s][s]( p[2] ) );
1164 float value = pointInfo.splineValues[0][j] * pointInfo.splineValues[1][k] * pointInfo.splineValues[2][l];
1165 diagonal += value * value * weight;
1166 pointInfo.weightedValue = value * weight;
1167 for(
int s=0 ; s<3 ; s++ ) pointInfo.splineValues[0][s] *= pointInfo.weightedValue;
1168 hasPoints[j][k] =
true;
1170 else samples[j][k][l].weightedValue = 0;
1177 neighbors5.neighbors[2][2][2]->depthAndOffset( d , off );
1178 int mn = 2 , mx = (1<<d)-2;
1179 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
1180 for(
int x=xStart ; x<=2 ; x++ )
1181 for(
int y=yStart ; y<yEnd ; y++ )
1182 if( x==2 && y>2 )
continue;
1183 else for(
int z=zStart ; z<zEnd ; z++ )
1184 if( x==2 && y==2 && z>2 )
continue;
1185 else if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1187 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
1188 int _index[] = { int( _node->off[0] ) , int( _node->off[1] ), int( _node->off[2] ) };
1190 if( isInterior ) temp =
Real( stencil[x][y][z] );
1191 else temp = GetLaplacian( node , _node );
1192 if( _constrainValues )
1194 int _d[] = { _index[0]-index[0] , _index[1]-index[1] , _index[2]-index[2] };
1195 if( x==2 && y==2 && z==2 ) temp += diagonal;
1196 else temp += GetValue( samples , hasPoints , _d );
1198 if( x==2 && y==2 && z==2 ) temp /= 2;
1199 if( fabs(temp)>MATRIX_ENTRY_EPSILON )
1201 row[count].N = _node->nodeData.nodeIndex-offset;
1202 row[count].Value = temp;
1208 template<
int Degree >
1209 void Octree< Degree >::SetDivergenceStencil(
int depth , Point3D< double > *stencil ,
bool scatter )
const
1211 int offset[] = { 2 , 2 , 2 };
1213 TreeOctNode::Index( depth , offset , d , off );
1214 int index1[3] , index2[3];
1215 if( scatter ) index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1216 else index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1217 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1219 int _offset[] = { x , y , z };
1220 TreeOctNode::Index( depth , _offset , d , off );
1221 if( scatter ) index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1222 else index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1225 BSplineData< Degree , Real >::SymmetricIndex( index1[0] , index2[0] ) ,
1226 BSplineData< Degree , Real >::SymmetricIndex( index1[1] , index2[1] ) ,
1227 BSplineData< Degree , Real >::SymmetricIndex( index1[2] , index2[2] ) ,
1231 #if GRADIENT_DOMAIN_SOLUTION
1233 fData.Index( index1[0] , index2[0] ) ,
1234 fData.Index( index1[1] , index2[1] ) ,
1235 fData.Index( index1[2] , index2[2] )
1236 #else // !GRADIENT_DOMAIN_SOLUTION
1238 fData.Index( index2[0] , index1[0] ) ,
1239 fData.Index( index2[1] , index1[1] ) ,
1240 fData.Index( index2[2] , index1[2] )
1241 #endif // GRADIENT_DOMAIN_SOLUTION
1244 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1245 #if GRADIENT_DOMAIN_SOLUTION
1246 Point3D<double> temp;
1247 temp[0] = fData.dvDotTable[aSymIndex[0]] * dot;
1248 temp[1] = fData.dvDotTable[aSymIndex[1]] * dot;
1249 temp[2] = fData.dvDotTable[aSymIndex[2]] * dot;
1250 stencil[25*x + 5*y + z] = temp;
1254 #else // !GRADIENT_DOMAIN_SOLUTION
1255 Point3D<double> temp;
1256 temp[0] = -fData.dvDotTable[aSymIndex[0]] * dot;
1257 temp[1] = -fData.dvDotTable[aSymIndex[1]] * dot;
1258 temp[2] = -fData.dvDotTable[aSymIndex[2]] * dot;
1259 stencil[25*x + 5*y + z] = temp;
1263 #endif // GRADIENT_DOMAIN_SOLUTION
1266 template<
int Degree >
1267 void Octree< Degree >::SetLaplacianStencil(
int depth ,
double stencil[5][5][5] )
const
1269 int offset[] = { 2 , 2 , 2 };
1271 TreeOctNode::Index( depth , offset , d , off );
1272 int index[] = { int( off[0] ) , int( off[1] ) , int( off[2] ) };
1273 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1275 int _offset[] = { x , y , z };
1277 TreeOctNode::Index( depth , _offset , _d , _off );
1278 int _index[] = { int( _off[0] ) , int( _off[1] ) , int( _off[2] ) };
1280 symIndex[0] = BSplineData< Degree , Real >::SymmetricIndex( index[0] , _index[0] );
1281 symIndex[1] = BSplineData< Degree , Real >::SymmetricIndex( index[1] , _index[1] );
1282 symIndex[2] = BSplineData< Degree , Real >::SymmetricIndex( index[2] , _index[2] );
1283 stencil[x][y][z] = GetLaplacian( symIndex );
1286 template<
int Degree >
1287 void Octree< Degree >::SetLaplacianStencils(
int depth , Stencil< double , 5 > stencils[2][2][2] )
const
1289 if( depth<=1 )
return;
1290 for(
int i=0 ; i<2 ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
1293 int offset[] = { 4+i , 4+j , 4+k };
1294 TreeOctNode::Index( depth , offset , d , off );
1295 int index[] = { int( off[0] ) , int( off[1] ) , int( off[2] ) };
1296 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1298 int _offset[] = { x , y , z };
1300 TreeOctNode::Index( depth-1 , _offset , _d , _off );
1301 int _index[] = { int( _off[0] ) , int( _off[1] ) , int( _off[2] ) };
1303 symIndex[0] = BSplineData< Degree , Real >::SymmetricIndex( index[0] , _index[0] );
1304 symIndex[1] = BSplineData< Degree , Real >::SymmetricIndex( index[1] , _index[1] );
1305 symIndex[2] = BSplineData< Degree , Real >::SymmetricIndex( index[2] , _index[2] );
1306 stencils[i][j][k].values[x][y][z] = GetLaplacian( symIndex );
1310 template<
int Degree >
1311 void Octree< Degree >::SetDivergenceStencils(
int depth , Stencil< Point3D< double > , 5 > stencils[2][2][2] ,
bool scatter )
const
1313 if( depth<=1 )
return;
1314 int index1[3] , index2[3];
1315 for(
int i=0 ; i<2 ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
1318 int offset[] = { 4+i , 4+j , 4+k };
1319 TreeOctNode::Index( depth , offset , d , off );
1320 if( scatter ) index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1321 else index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1322 for(
int x=0 ; x<5 ; x++ )
for(
int y=0 ; y<5 ; y++ )
for(
int z=0 ; z<5 ; z++ )
1324 int _offset[] = { x , y , z };
1325 TreeOctNode::Index( depth-1 , _offset , d , off );
1326 if( scatter ) index1[0] = int( off[0] ) , index1[1] = int( off[1] ) , index1[2] = int( off[2] );
1327 else index2[0] = int( off[0] ) , index2[1] = int( off[1] ) , index2[2] = int( off[2] );
1331 BSplineData< Degree , Real >::SymmetricIndex( index1[0] , index2[0] ) ,
1332 BSplineData< Degree , Real >::SymmetricIndex( index1[1] , index2[1] ) ,
1333 BSplineData< Degree , Real >::SymmetricIndex( index1[2] , index2[2] ) ,
1337 #if GRADIENT_DOMAIN_SOLUTION
1339 fData.Index( index1[0] , index2[0] ) ,
1340 fData.Index( index1[1] , index2[1] ) ,
1341 fData.Index( index1[2] , index2[2] )
1342 #else // !GRADIENT_DOMAIN_SOLUTION
1344 fData.Index( index2[0] , index1[0] ) ,
1345 fData.Index( index2[1] , index1[1] ) ,
1346 fData.Index( index2[2] , index1[2] )
1347 #endif // GRADIENT_DOMAIN_SOLUTION
1349 double dot = fData.vvDotTable[symIndex[0]] * fData.vvDotTable[symIndex[1]] * fData.vvDotTable[symIndex[2]];
1350 #if GRADIENT_DOMAIN_SOLUTION
1351 stencils[i][j][k].values[x][y][z][0] = fData.dvDotTable[aSymIndex[0]] * dot;
1352 stencils[i][j][k].values[x][y][z][1] = fData.dvDotTable[aSymIndex[1]] * dot;
1353 stencils[i][j][k].values[x][y][z][2] = fData.dvDotTable[aSymIndex[2]] * dot;
1354 #else // !GRADIENT_DOMAIN_SOLUTION
1355 stencils[i][j][k].values[x][y][z][0] = -fData.dvDotTable[aSymIndex[0]] * dot;
1356 stencils[i][j][k].values[x][y][z][1] = -fData.dvDotTable[aSymIndex[1]] * dot;
1357 stencils[i][j][k].values[x][y][z][2] = -fData.dvDotTable[aSymIndex[2]] * dot;
1358 #endif // GRADIENT_DOMAIN_SOLUTION
1362 template<
int Degree >
1363 void Octree< Degree >::SetEvaluationStencils(
int depth , Stencil< double , 3 > stencil1[8] , Stencil< double , 3 > stencil2[8][8] )
const
1368 int off[] = { 2 , 2 , 2 };
1369 for(
int c=0 ; c<8 ; c++ )
1371 VertexData::CornerIndex( depth , off , c , fData.depth , idx );
1372 idx[0] *= fData.functionCount , idx[1] *= fData.functionCount , idx[2] *= fData.functionCount;
1373 for(
int x=0 ; x<3 ; x++ )
for(
int y=0 ; y<3 ; y++ )
for(
int z=0 ; z<3 ; z++ )
1376 int _offset[] = { x+1 , y+1 , z+1 };
1377 TreeOctNode::Index( depth , _offset , _d , _off );
1378 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]) ];
1383 for(
int _c=0 ; _c<8 ; _c++ )
1386 int _cx , _cy , _cz;
1387 Cube::FactorCornerIndex( _c , _cx , _cy , _cz );
1388 int off[] = { 4+_cx , 4+_cy , 4+_cz };
1389 for(
int c=0 ; c<8 ; c++ )
1391 VertexData::CornerIndex( depth , off , c , fData.depth , idx );
1392 idx[0] *= fData.functionCount , idx[1] *= fData.functionCount , idx[2] *= fData.functionCount;
1393 for(
int x=0 ; x<3 ; x++ )
for(
int y=0 ; y<3 ; y++ )
for(
int z=0 ; z<3 ; z++ )
1396 int _offset[] = { x+1 , y+1 , z+1 };
1397 TreeOctNode::Index( depth-1 , _offset , _d , _off );
1398 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]) ];
1403 template<
int Degree >
1404 void Octree< Degree >::UpdateCoarserSupportBounds(
const TreeOctNode* node ,
int& startX ,
int& endX ,
int& startY ,
int& endY ,
int& startZ ,
int& endZ )
1408 int x , y , z , c = int( node - node->parent->children );
1409 Cube::FactorCornerIndex( c , x , y , z );
1410 if( x==0 ) endX = 4;
1412 if( y==0 ) endY = 4;
1414 if( z==0 ) endZ = 4;
1419 template<
int Degree >
1420 void Octree< Degree >::UpdateConstraintsFromCoarser(
const OctNode< TreeNodeData , Real >::NeighborKey5& neighborKey5 , TreeOctNode* node , Real* metSolution ,
const Stencil< double , 5 >& lapStencil )
const
1425 node->depthAndOffset( d , off );
1426 int mn = 4 , mx = (1<<d)-4;
1427 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
1430 int depth = node->depth();
1431 if( depth<=_minDepth )
return;
1432 int i = node->nodeData.nodeIndex;
1434 int startX = 0 , endX = 5 , startY = 0 , endY = 5 , startZ = 0 , endZ = 5;
1435 UpdateCoarserSupportBounds( node , startX , endX , startY , endY , startZ , endZ );
1437 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth-1];
1438 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
1439 if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.nodeIndex>=0 )
1441 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
1442 Real _solution = metSolution[ _node->nodeData.nodeIndex ];
1444 if( isInterior ) node->nodeData.constraint -=
Real( lapStencil.values[x][y][z] * _solution );
1445 else node->nodeData.constraint -= GetLaplacian( _node , node ) * _solution;
1448 if( _constrainValues )
1451 node->depthAndOffset( d, idx );
1452 idx[0] = BinaryNode< double >::CenterIndex( d , idx[0] );
1453 idx[1] = BinaryNode< double >::CenterIndex( d , idx[1] );
1454 idx[2] = BinaryNode< double >::CenterIndex( d , idx[2] );
1455 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth];
1456 for(
int x=1 ; x<4 ; x++ )
for(
int y=1 ; y<4 ; y++ )
for(
int z=1 ; z<4 ; z++ )
1457 if( neighbors5.neighbors[x][y][z] && neighbors5.neighbors[x][y][z]->nodeData.pointIndex!=-1 )
1459 const PointData& pData = _points[ neighbors5.neighbors[x][y][z]->nodeData.pointIndex ];
1460 Real pointValue = pData.value;
1461 Point3D< Real > p = pData.position;
1462 node->nodeData.constraint -=
1464 fData.baseBSplines[idx[0]][x-1]( p[0] ) *
1465 fData.baseBSplines[idx[1]][y-1]( p[1] ) *
1466 fData.baseBSplines[idx[2]][z-1]( p[2] ) *
1477 UpSampleData(
int s ,
double v1 ,
double v2 ) { start = s , v[0] = v1 , v[1] = v2; }
1479 template<
int Degree >
1480 void Octree< Degree >::UpSampleCoarserSolution(
int depth ,
const SortedTreeNodes& sNodes , Vector< Real >& Solution )
const
1482 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1483 Solution.Resize( range );
1484 if( !depth )
return;
1485 else if( depth==1 )
for(
int i=start ; i<end ; i++ ) Solution[i-start] += sNodes.treeNodes[0]->nodeData.solution;
1489 #pragma omp parallel
for num_threads( threads )
1490 for(
int t=0 ; t<threads ; t++ )
1492 TreeOctNode::NeighborKey3 neighborKey;
1493 neighborKey.set( depth );
1494 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1497 UpSampleData usData[3];
1498 sNodes.treeNodes[i]->depthAndOffset( d , off );
1499 for(
int d=0 ; d<3 ; d++ )
1501 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1502 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1503 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1504 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1506 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1507 for(
int ii=0 ; ii<2 ; ii++ )
1509 int _ii = ii + usData[0].start;
1510 double dx = usData[0].v[ii];
1511 for(
int jj=0 ; jj<2 ; jj++ )
1513 int _jj = jj + usData[1].start;
1514 double dxy = dx * usData[1].v[jj];
1515 for(
int kk=0 ; kk<2 ; kk++ )
1517 int _kk = kk + usData[2].start;
1518 double dxyz = dxy * usData[2].v[kk];
1519 Solution[i-start] += Real( neighborKey.neighbors[depth-1].neighbors[_ii][_jj][_kk]->nodeData.solution * dxyz );
1527 start = sNodes.nodeCount[depth-1] , end = sNodes.nodeCount[depth] , range = end-start;
1528 #pragma omp parallel for num_threads( threads )
1529 for(
int i=start ; i<end ; i++ ) sNodes.treeNodes[i]->nodeData.solution =
Real( 0. );
1531 template<
int Degree >
1532 void Octree< Degree >::DownSampleFinerConstraints(
int depth , SortedTreeNodes& sNodes )
const
1534 if( !depth )
return;
1535 #pragma omp parallel for num_threads( threads )
1536 for(
int i=sNodes.nodeCount[depth-1] ; i<sNodes.nodeCount[depth] ; i++ )
1537 sNodes.treeNodes[i]->nodeData.constraint =
Real( 0 );
1541 sNodes.treeNodes[0]->nodeData.constraint =
Real( 0 );
1542 for(
int i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ ) sNodes.treeNodes[0]->nodeData.constraint += sNodes.treeNodes[i]->nodeData.constraint;
1545 std::vector< Vector< double > > constraints( threads );
1546 for(
int t=0 ; t<threads ; t++ ) constraints[t].Resize( sNodes.nodeCount[depth] - sNodes.nodeCount[depth-1] ) , constraints[t].SetZero();
1547 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1548 int lStart = sNodes.nodeCount[depth-1] , lEnd = sNodes.nodeCount[depth];
1550 #pragma omp parallel for num_threads( threads )
1551 for(
int t=0 ; t<threads ; t++ )
1553 TreeOctNode::NeighborKey3 neighborKey;
1554 neighborKey.set( depth );
1555 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1558 UpSampleData usData[3];
1559 sNodes.treeNodes[i]->depthAndOffset( d , off );
1560 for(
int d=0 ; d<3 ; d++ )
1562 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1563 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1564 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1565 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1567 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1568 TreeOctNode::Neighbors3& neighbors = neighborKey.neighbors[depth-1];
1569 for(
int ii=0 ; ii<2 ; ii++ )
1571 int _ii = ii + usData[0].start;
1572 double dx = usData[0].v[ii];
1573 for(
int jj=0 ; jj<2 ; jj++ )
1575 int _jj = jj + usData[1].start;
1576 double dxy = dx * usData[1].v[jj];
1577 for(
int kk=0 ; kk<2 ; kk++ )
1579 int _kk = kk + usData[2].start;
1580 double dxyz = dxy * usData[2].v[kk];
1581 constraints[t][neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex-lStart] += sNodes.treeNodes[i]->nodeData.constraint * dxyz;
1587 #pragma omp parallel for num_threads( threads )
1588 for(
int i=lStart ; i<lEnd ; i++ )
1591 for(
int t=0 ; t<threads ; t++ ) cSum += constraints[t][i-lStart];
1592 sNodes.treeNodes[i]->nodeData.constraint += cSum;
1595 template<
int Degree >
1597 void Octree< Degree >::DownSample(
int depth ,
const SortedTreeNodes& sNodes , C* constraints )
const
1599 if( depth==0 )
return;
1602 for(
int i=sNodes.nodeCount[1] ; i<sNodes.nodeCount[2] ; i++ ) constraints[0] += constraints[i];
1605 std::vector< Vector< C > > _constraints( threads );
1606 for(
int t=0 ; t<threads ; t++ ) _constraints[t].Resize( sNodes.nodeCount[depth] - sNodes.nodeCount[depth-1] );
1607 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start , lStart = sNodes.nodeCount[depth-1] , lEnd = sNodes.nodeCount[depth];
1609 #pragma omp parallel for num_threads( threads )
1610 for(
int t=0 ; t<threads ; t++ )
1612 TreeOctNode::NeighborKey3 neighborKey;
1613 neighborKey.set( depth );
1614 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1617 UpSampleData usData[3];
1618 sNodes.treeNodes[i]->depthAndOffset( d , off );
1619 for(
int d=0 ; d<3 ; d++ )
1621 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1622 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1623 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1624 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1626 TreeOctNode::Neighbors3& neighbors = neighborKey.getNeighbors( sNodes.treeNodes[i]->parent );
1627 C c = constraints[i];
1628 for(
int ii=0 ; ii<2 ; ii++ )
1630 int _ii = ii + usData[0].start;
1631 C cx = C( c*usData[0].v[ii] );
1632 for(
int jj=0 ; jj<2 ; jj++ )
1634 int _jj = jj + usData[1].start;
1635 C cxy = C( cx*usData[1].v[jj] );
1636 for(
int kk=0 ; kk<2 ; kk++ )
1638 int _kk = kk + usData[2].start;
1639 if( neighbors.neighbors[_ii][_jj][_kk] )
1640 _constraints[t][neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex-lStart] += C( cxy*usData[2].v[kk] );
1646 #pragma omp parallel for num_threads( threads )
1647 for(
int i=lStart ; i<lEnd ; i++ )
1650 for(
int t=0 ; t<threads ; t++ ) cSum += _constraints[t][i-lStart];
1651 constraints[i] += cSum;
1654 template<
int Degree >
1656 void Octree< Degree >::UpSample(
int depth ,
const SortedTreeNodes& sNodes , C* coefficients )
const
1658 if ( depth==0 )
return;
1661 for(
int i=sNodes.nodeCount[1] ; i<sNodes.nodeCount[2] ; i++ ) coefficients[i] += coefficients[0];
1665 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1667 #pragma omp parallel for num_threads( threads )
1668 for(
int t=0 ; t<threads ; t++ )
1670 TreeOctNode::NeighborKey3 neighborKey;
1671 neighborKey.set( depth-1 );
1672 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1676 UpSampleData usData[3];
1678 for(
int d=0 ; d<3 ; d++ )
1680 if ( off[d] ==0 ) usData[d] = UpSampleData( 1 , 1.00 , 0.00 );
1681 else if( off[d]+1==(1<<depth) ) usData[d] = UpSampleData( 0 , 0.00 , 1.00 );
1682 else if( off[d]%2 ) usData[d] = UpSampleData( 1 , 0.75 , 0.25 );
1683 else usData[d] = UpSampleData( 0 , 0.25 , 0.75 );
1685 TreeOctNode::Neighbors3& neighbors = neighborKey.getNeighbors( node->parent );
1686 for(
int ii=0 ; ii<2 ; ii++ )
1688 int _ii = ii + usData[0].start;
1689 double dx = usData[0].v[ii];
1690 for(
int jj=0 ; jj<2 ; jj++ )
1692 int _jj = jj + usData[1].start;
1693 double dxy = dx * usData[1].v[jj];
1694 for(
int kk=0 ; kk<2 ; kk++ )
1696 int _kk = kk + usData[2].start;
1697 if( neighbors.neighbors[_ii][_jj][_kk] )
1699 double dxyz = dxy * usData[2].v[kk];
1700 int _i = neighbors.neighbors[_ii][_jj][_kk]->nodeData.nodeIndex;
1701 coefficients[i] += coefficients[_i] *
Real( dxyz );
1709 template<
int Degree >
1710 void Octree< Degree >::SetCoarserPointValues(
int depth ,
const SortedTreeNodes& sNodes , Real* metSolution )
1712 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1714 #pragma omp parallel for num_threads( threads )
1715 for(
int t=0 ; t<threads ; t++ )
1717 TreeOctNode::NeighborKey3 neighborKey;
1718 neighborKey.set( depth );
1719 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
1721 int pIdx = sNodes.treeNodes[i]->nodeData.pointIndex;
1724 neighborKey.getNeighbors( sNodes.treeNodes[i] );
1725 _points[ pIdx ].value = WeightedCoarserFunctionValue( neighborKey , sNodes.treeNodes[i] , metSolution );
1730 template<
int Degree >
1731 Real Octree< Degree >::WeightedCoarserFunctionValue(
const OctNode< TreeNodeData , Real >::NeighborKey3& neighborKey ,
const TreeOctNode* pointNode , Real* metSolution )
const
1733 int depth = pointNode->depth();
1734 if( !depth || pointNode->nodeData.pointIndex==-1 )
return Real(0.);
1735 double pointValue = 0;
1737 Real weight = _points[ pointNode->nodeData.pointIndex ].weight;
1738 Point3D< Real > p = _points[ pointNode->nodeData.pointIndex ].position;
1743 const TreeOctNode::Neighbors3& neighbors = neighborKey.neighbors[depth-1];
1744 neighbors.neighbors[1][1][1]->depthAndOffset( d , _idx );
1745 _idx[0] = BinaryNode< double >::CenterIndex( d , _idx[0]-1 );
1746 _idx[1] = BinaryNode< double >::CenterIndex( d , _idx[1]-1 );
1747 _idx[2] = BinaryNode< double >::CenterIndex( d , _idx[2]-1 );
1748 for(
int j=0 ; j<3 ; j++ )
for(
int k=0 ; k<3 ; k++ )
for(
int l=0 ; l<3 ; l++ )
1749 if( neighbors.neighbors[j][k][l] && neighbors.neighbors[j][k][l]->nodeData.nodeIndex>=0 )
1752 const TreeOctNode* basisNode = neighbors.neighbors[j][k][l];
1753 int idx[] = { _idx[0]+j , _idx[1]+k , _idx[2]+l };
1755 fData.baseBSplines[ idx[0] ][2-j]( p[0] ) *
1756 fData.baseBSplines[ idx[1] ][2-k]( p[1] ) *
1757 fData.baseBSplines[ idx[2] ][2-l]( p[2] ) *
1758 metSolution[basisNode->nodeData.nodeIndex];
1761 return Real( pointValue * weight );
1763 template<
int Degree>
1764 int Octree<Degree>::GetFixedDepthLaplacian( SparseSymmetricMatrix< Real >& matrix ,
int depth ,
const SortedTreeNodes& sNodes , Real* metSolution )
1766 int start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start;
1767 double stencil[5][5][5];
1768 SetLaplacianStencil( depth , stencil );
1769 Stencil< double , 5 > stencils[2][2][2];
1770 SetLaplacianStencils( depth , stencils );
1771 matrix.Resize( range );
1772 #pragma omp parallel for num_threads( threads )
1773 for(
int t=0 ; t<threads ; t++ )
1775 TreeOctNode::NeighborKey5 neighborKey5;
1776 neighborKey5.set( depth );
1777 for(
int i=(range*t)/threads ; i<(range*(t+1))/threads ; i++ )
1780 neighborKey5.getNeighbors( node );
1783 int count = GetMatrixRowSize( neighborKey5.neighbors[depth] );
1786 #pragma omp critical (matrix_set_row_size)
1788 matrix.SetRowSize( i , count );
1792 matrix.rowSizes[i] = SetMatrixRow( neighborKey5.neighbors[depth] , matrix[i] , start , stencil );
1798 c = int( node - node->parent->children );
1799 Cube::FactorCornerIndex( c , x , y , z );
1802 UpdateConstraintsFromCoarser( neighborKey5 , node , metSolution , stencils[x][y][z] );
1807 template<
int Degree>
1808 int Octree<Degree>::GetRestrictedFixedDepthLaplacian( SparseSymmetricMatrix< Real >& matrix,
int depth,
const int* entries,
int entryCount,
1809 const TreeOctNode* rNode , Real radius ,
1810 const SortedTreeNodes& sNodes , Real* metSolution )
1812 for(
int i=0 ; i<entryCount ; i++ ) sNodes.treeNodes[ entries[i] ]->nodeData.nodeIndex = i;
1813 double stencil[5][5][5];
1814 int rDepth , rOff[3];
1815 rNode->depthAndOffset( rDepth , rOff );
1816 matrix.Resize( entryCount );
1817 SetLaplacianStencil( depth , stencil );
1818 Stencil< double , 5 > stencils[2][2][2];
1819 SetLaplacianStencils( depth , stencils );
1820 #pragma omp parallel for num_threads( threads )
1821 for(
int t=0 ; t<threads ; t++ )
1823 TreeOctNode::NeighborKey5 neighborKey5;
1824 neighborKey5.set( depth );
1825 for(
int i=(entryCount*t)/threads ; i<(entryCount*(t+1))/threads ; i++ )
1827 TreeOctNode* node = sNodes.treeNodes[ entries[i] ];
1830 off[0] >>= (depth-rDepth) , off[1] >>= (depth-rDepth) , off[2] >>= (depth-rDepth);
1831 bool isInterior = ( off[0]==rOff[0] && off[1]==rOff[1] && off[2]==rOff[2] );
1833 neighborKey5.getNeighbors( node );
1835 int xStart=0 , xEnd=5 , yStart=0 , yEnd=5 , zStart=0 , zEnd=5;
1836 if( !isInterior ) SetMatrixRowBounds( neighborKey5.neighbors[depth].neighbors[2][2][2] , rDepth , rOff , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1839 int count = GetMatrixRowSize( neighborKey5.neighbors[depth] , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1842 #pragma omp critical (matrix_set_row_size)
1844 matrix.SetRowSize( i , count );
1848 matrix.rowSizes[i] = SetMatrixRow( neighborKey5.neighbors[depth] , matrix[i] , 0 , stencil , xStart , xEnd , yStart , yEnd , zStart , zEnd );
1853 c = int( node - node->parent->children );
1854 Cube::FactorCornerIndex( c , x , y , z );
1857 UpdateConstraintsFromCoarser( neighborKey5 , node , metSolution , stencils[x][y][z] );
1860 for(
int i=0 ; i<entryCount ; i++ ) sNodes.treeNodes[entries[i]]->nodeData.nodeIndex = entries[i];
1864 template<
int Degree>
1869 fData.setDotTables( fData.DD_DOT_FLAG | fData.DV_DOT_FLAG );
1872 _sNodes.treeNodes[0]->nodeData.solution = 0;
1874 std::vector< Real > metSolution( _sNodes.nodeCount[ _sNodes.maxDepth ] , 0 );
1876 for( i=1 ; i<_sNodes.maxDepth ; i++ )
1878 if( subdivideDepth>0 ) iter += SolveFixedDepthMatrix( i , _sNodes , &metSolution[0] , subdivideDepth , showResidual , minIters , accuracy );
1879 else iter += SolveFixedDepthMatrix( i , _sNodes , &metSolution[0] , showResidual , minIters , accuracy );
1882 fData.clearDotTables( fData.VV_DOT_FLAG | fData.DV_DOT_FLAG | fData.DD_DOT_FLAG );
1887 template<
int Degree>
1893 double systemTime=0. , solveTime=0. , updateTime=0. , evaluateTime = 0.;
1896 if( depth<=_minDepth ) UpSampleCoarserSolution( depth , sNodes , X );
1900 UpSample( depth-1 , sNodes , metSolution );
1903 #pragma omp parallel for num_threads( threads )
1904 for(
int i=_sNodes.nodeCount[depth-1] ; i<_sNodes.nodeCount[depth] ; i++ ) metSolution[i] += _sNodes.treeNodes[i]->nodeData.solution;
1906 if( _constrainValues )
1908 SetCoarserPointValues( depth , sNodes , metSolution );
1911 SparseSymmetricMatrix< Real >::internalAllocator.rollBack();
1913 int maxECount = ( (2*Degree+1)*(2*Degree+1)*(2*Degree+1) + 1 ) / 2;
1914 maxECount = ( ( maxECount + 15 ) / 16 ) * 16;
1921 SparseSymmetricMatrix< Real >::internalAllocator.rollBack();
1922 GetFixedDepthLaplacian( M , depth , sNodes , metSolution );
1929 iter += SparseSymmetricMatrix< Real >::Solve( M , B , std::max< int >(
int( pow( M.
rows , ITERATION_POWER ) ) , minIters ) , X ,
Real(accuracy) , 0 , threads , (depth<=_minDepth) && !_constrainValues );
1934 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;
1935 double bNorm = B.
Norm( 2 ) , rNorm = ( B - M * X ).Norm( 2 );
1936 printf(
"\tResidual: (%d %g) %g -> %g (%f) [%d]\n" , M.
Entries() , sqrt(mNorm) , bNorm , rNorm , rNorm/bNorm , iter );
1944 template<
int Degree>
1945 int Octree<Degree>::SolveFixedDepthMatrix(
int depth ,
const SortedTreeNodes& sNodes , Real* metSolution ,
int startingDepth ,
bool showResidual ,
int minIters ,
double accuracy )
1947 if( startingDepth>=depth )
return SolveFixedDepthMatrix( depth , sNodes , metSolution , showResidual , minIters , accuracy );
1949 int i , j , d , tIter=0;
1950 SparseSymmetricMatrix< Real > _M;
1951 Vector< Real >
B , _B , _X;
1952 AdjacencySetFunction asf;
1953 AdjacencyCountFunction acf;
1954 double systemTime = 0 , solveTime = 0 , memUsage = 0 , evaluateTime = 0 , gTime = 0, sTime = 0;
1955 Real myRadius , myRadius2;
1957 if( depth>_minDepth )
1960 UpSample( depth-1 , sNodes , metSolution );
1963 #pragma omp parallel for num_threads( threads )
1964 for(
int i=_sNodes.nodeCount[depth-1] ; i<_sNodes.nodeCount[depth] ; i++ ) metSolution[i] += _sNodes.treeNodes[i]->nodeData.solution;
1967 if( _constrainValues )
1969 SetCoarserPointValues( depth , sNodes , metSolution );
1971 B.Resize( sNodes.nodeCount[depth+1] - sNodes.nodeCount[depth] );
1974 for( i=sNodes.nodeCount[depth] ; i<sNodes.nodeCount[depth+1] ; i++ )
1976 B[ i-sNodes.nodeCount[depth] ] = sNodes.treeNodes[i]->nodeData.constraint;
1977 sNodes.treeNodes[i]->nodeData.constraint = 0;
1980 myRadius = 2*radius-
Real(0.5);
1981 myRadius = int(myRadius-ROUND_EPS)+
ROUND_EPS;
1982 myRadius2 =
Real(radius+ROUND_EPS-0.5);
1983 d = depth-startingDepth;
1984 std::vector< int > subDimension( sNodes.nodeCount[d+1]-sNodes.nodeCount[d] );
1985 int maxDimension = 0;
1986 for( i=sNodes.nodeCount[d] ; i<sNodes.nodeCount[d+1] ; i++ )
1989 acf.adjacencyCount = 0;
1990 for( TreeOctNode* temp=sNodes.treeNodes[i]->nextNode() ; temp ; )
1992 if( temp->depth()==depth ) acf.adjacencyCount++ , temp = sNodes.treeNodes[i]->nextBranch( temp );
1993 else temp = sNodes.treeNodes[i]->nextNode ( temp );
1995 for( j=sNodes.nodeCount[d] ; j<sNodes.nodeCount[d+1] ; j++ )
1997 if( i==j )
continue;
1998 TreeOctNode::ProcessFixedDepthNodeAdjacentNodes( fData.depth , sNodes.treeNodes[i] , 1 , sNodes.treeNodes[j] , 2*width-1 , depth , &acf );
2000 subDimension[i-sNodes.nodeCount[d]] = acf.adjacencyCount;
2001 maxDimension = std::max< int >( maxDimension , subDimension[i-sNodes.nodeCount[d]] );
2003 asf.adjacencies =
new int[maxDimension];
2004 MapReduceVector< Real > mrVector;
2005 mrVector.resize( threads , maxDimension );
2007 for( i=sNodes.nodeCount[d] ; i<sNodes.nodeCount[d+1] ; i++ )
2011 acf.adjacencyCount = subDimension[i-sNodes.nodeCount[d]];
2012 if( !acf.adjacencyCount )
continue;
2015 asf.adjacencyCount = 0;
2016 for( TreeOctNode* temp=sNodes.treeNodes[i]->nextNode() ; temp ; )
2018 if( temp->depth()==depth ) asf.adjacencies[ asf.adjacencyCount++ ] = temp->nodeData.nodeIndex , temp = sNodes.treeNodes[i]->nextBranch( temp );
2019 else temp = sNodes.treeNodes[i]->nextNode ( temp );
2021 for( j=sNodes.nodeCount[d] ; j<sNodes.nodeCount[d+1] ; j++ )
2023 if( i==j )
continue;
2024 TreeOctNode::ProcessFixedDepthNodeAdjacentNodes( fData.depth , sNodes.treeNodes[i] , 1 , sNodes.treeNodes[j] , 2*width-1 , depth , &asf );
2028 _B.Resize( asf.adjacencyCount );
2029 for( j=0 ; j<asf.adjacencyCount ; j++ ) _B[j] = B[ asf.adjacencies[j]-sNodes.nodeCount[depth] ];
2031 _X.Resize( asf.adjacencyCount );
2032 #pragma omp parallel for num_threads( threads ) schedule( static )
2033 for( j=0 ; j<asf.adjacencyCount ; j++ )
2035 _X[j] = sNodes.treeNodes[ asf.adjacencies[j] ]->nodeData.solution;
2038 SparseSymmetricMatrix< Real >::internalAllocator.rollBack();
2039 GetRestrictedFixedDepthLaplacian( _M , depth , asf.adjacencies , asf.adjacencyCount , sNodes.treeNodes[i] , myRadius , sNodes , metSolution );
2040 #pragma omp parallel for num_threads( threads ) schedule( static )
2041 for( j=0 ; j<asf.adjacencyCount ; j++ )
2043 _B[j] += sNodes.treeNodes[asf.adjacencies[j]]->nodeData.constraint;
2044 sNodes.treeNodes[ asf.adjacencies[j] ]->nodeData.constraint = 0;
2049 iter += SparseSymmetricMatrix< Real >::Solve( _M , _B , std::max< int >(
int( pow( _M.rows , ITERATION_POWER ) ) , minIters ) , _X , mrVector ,
Real(accuracy) , 0 );
2054 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;
2055 double bNorm = _B.Norm( 2 ) , rNorm = ( _B - _M * _X ).Norm( 2 );
2056 printf(
"\t\tResidual: (%d %g) %g -> %g (%f) [%d]\n" , _M.Entries() , sqrt(mNorm) , bNorm , rNorm , rNorm/bNorm , iter );
2060 for( j=0 ; j<asf.adjacencyCount ; j++ )
2062 TreeOctNode* temp=sNodes.treeNodes[ asf.adjacencies[j] ];
2063 while( temp->depth()>sNodes.treeNodes[i]->depth() ) temp=temp->
parent;
2064 if( temp->nodeData.nodeIndex>=sNodes.treeNodes[i]->nodeData.nodeIndex ) sNodes.treeNodes[ asf.adjacencies[j] ]->
nodeData.solution =
Real( _X[j] );
2066 systemTime += gTime;
2068 memUsage = std::max< double >( MemoryUsage() , memUsage );
2071 delete[] asf.adjacencies;
2074 template<
int Degree>
2075 int Octree<Degree>::HasNormals(TreeOctNode* node,Real epsilon)
2078 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;
2079 if( node->children )
for(
int i=0;i<Cube::CORNERS && !hasNormals;i++) hasNormals |= HasNormals(&node->children[i],epsilon);
2083 template<
int Degree>
2086 int maxDepth = tree.maxDepth();
2088 if( temp->children && temp->d>=_minDepth )
2091 for(
int i=0 ; i<Cube::CORNERS && !hasNormals ; i++ ) hasNormals = HasNormals( &temp->children[i] ,
EPSILON/(1<<maxDepth) );
2092 if( !hasNormals ) temp->children=NULL;
2096 template<
int Degree>
2104 fData.setDotTables( fData.VV_DOT_FLAG | fData.DV_DOT_FLAG );
2106 int maxDepth = _sNodes.maxDepth-1;
2108 zeroPoint[0] = zeroPoint[1] = zeroPoint[2] = 0;
2109 std::vector< Real > constraints( _sNodes.nodeCount[maxDepth] ,
Real(0) );
2110 std::vector< Point3D< Real > > coefficients( _sNodes.nodeCount[maxDepth] , zeroPoint );
2113 #pragma omp parallel for num_threads( threads )
2114 for(
int i=0 ; i<_sNodes.nodeCount[maxDepth+1] ; i++ ) _sNodes.treeNodes[i]->nodeData.constraint =
Real( 0. );
2117 std::vector< std::vector< Real > > _constraints( threads );
2118 for(
int t=0 ; t<threads ; t++ ) _constraints[t].resize( _sNodes.nodeCount[maxDepth] , 0 );
2120 for(
int d=maxDepth ; d>=0 ; d-- )
2123 SetDivergenceStencil( d , &stencil[0][0][0] ,
false );
2124 Stencil< Point3D< double > , 5 > stencils[2][2][2];
2125 SetDivergenceStencils( d , stencils ,
true );
2126 #pragma omp parallel for num_threads( threads )
2127 for(
int t=0 ; t<threads ; t++ )
2129 TreeOctNode::NeighborKey5 neighborKey5;
2130 neighborKey5.set( fData.depth );
2131 int start = _sNodes.nodeCount[d] , end = _sNodes.nodeCount[d+1] , range = end-start;
2132 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
2135 int startX=0 , endX=5 , startY=0 , endY=5 , startZ=0 , endZ=5;
2136 int depth = node->
depth();
2137 neighborKey5.getNeighbors( node );
2139 bool isInterior , isInterior2;
2143 int mn = 2 , mx = (1<<d)-2;
2144 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
2146 isInterior2 = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
2152 Cube::FactorCornerIndex( c , cx , cy , cz );
2154 else cx = cy = cz = 0;
2155 Stencil< Point3D< double > , 5 >& _stencil = stencils[cx][cy][cz];
2159 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth];
2162 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2164 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
2168 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] );
2172 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2174 const TreeOctNode* _node = neighbors5.neighbors[x][y][z];
2181 UpdateCoarserSupportBounds( neighbors5.neighbors[2][2][2] , startX , endX , startY , endY , startZ , endZ );
2185 if( normal[0]==0 && normal[1]==0 && normal[2]==0 )
continue;
2186 if( depth<maxDepth ) coefficients[i] += normal;
2190 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.neighbors[depth-1];
2192 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2193 if( neighbors5.neighbors[x][y][z] )
2195 TreeOctNode* _node = neighbors5.neighbors[x][y][z];
2199 _constraints[t][ _node->
nodeData.
nodeIndex ] +=
Real( div[0] * normal[0] + div[1] * normal[1] + div[2] * normal[2] );
2201 else _constraints[t][ _node->
nodeData.
nodeIndex ] += GetDivergence( node , _node , normal );
2207 #pragma omp parallel for num_threads( threads ) schedule( static )
2208 for(
int i=0 ; i<_sNodes.nodeCount[maxDepth] ; i++ )
2211 for(
int t=0 ; t<threads ; t++ ) cSum += _constraints[t][i];
2212 constraints[i] = cSum;
2215 for(
int d=maxDepth-1 ; d>=0 ; d-- ) DownSample( d , _sNodes , &constraints[0] );
2218 for(
int d=0 ; d<maxDepth ; d++ ) UpSample( d , _sNodes , &coefficients[0] );
2221 #pragma omp parallel for num_threads( threads )
2222 for(
int i=0 ; i<_sNodes.nodeCount[maxDepth] ; i++ ) _sNodes.treeNodes[i]->nodeData.constraint += constraints[i];
2225 for(
int d=0 ; d<=maxDepth ; d++ )
2227 int start = _sNodes.nodeCount[d] , end = _sNodes.nodeCount[d+1] , range = end - start;
2228 Stencil< Point3D< double > , 5 > stencils[2][2][2];
2229 SetDivergenceStencils( d , stencils ,
false );
2230 #pragma omp parallel for num_threads( threads )
2231 for(
int t=0 ; t<threads ; t++ )
2233 TreeOctNode::NeighborKey5 neighborKey5;
2234 neighborKey5.set( maxDepth );
2235 for(
int i=start+(range*t)/threads ; i<start+(range*(t+1))/threads ; i++ )
2238 int depth = node->
depth();
2239 if( !depth )
continue;
2240 int startX=0 , endX=5 , startY=0 , endY=5 , startZ=0 , endZ=5;
2241 UpdateCoarserSupportBounds( node , startX , endX , startY , endY , startZ , endZ );
2242 const TreeOctNode::Neighbors5& neighbors5 = neighborKey5.getNeighbors( node->
parent );
2248 int mn = 4 , mx = (1<<d)-4;
2249 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
2255 Cube::FactorCornerIndex( c , cx , cy , cz );
2257 else cx = cy = cz = 0;
2258 Stencil< Point3D< double > , 5 >& _stencil = stencils[cx][cy][cz];
2261 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2262 if( neighbors5.neighbors[x][y][z] )
2264 TreeOctNode* _node = neighbors5.neighbors[x][y][z];
2270 constraint +=
Real( div[0] * normal[0] + div[1] * normal[1] + div[2] * normal[2] );
2272 else constraint += GetDivergence( _node , node , coefficients[_i] );
2279 fData.clearDotTables( fData.DV_DOT_FLAG );
2282 #pragma omp parallel for num_threads( threads )
2283 for(
int t=0 ; t<threads ; t++ )
2284 for(
int i=(_sNodes.nodeCount[maxDepth+1]*t)/threads ; i<(_sNodes.nodeCount[maxDepth+1]*(t+1))/threads ; i++ )
2294 template<
int Degree>
2296 template<
int Degree>
2297 void Octree<Degree>::AdjacencySetFunction::Function(
const TreeOctNode* node1,
const TreeOctNode* node2){adjacencies[adjacencyCount++]=node1->nodeData.nodeIndex;}
2298 template<
int Degree>
2299 void Octree<Degree>::RefineFunction::Function( TreeOctNode* node1 ,
const TreeOctNode* node2 )
2301 if( !node1->children && node1->depth()<depth ) node1->initChildren();
2303 template<
int Degree >
2304 void Octree< Degree >::FaceEdgesFunction::Function(
const TreeOctNode* node1 ,
const TreeOctNode* node2 )
2306 if( !node1->children && MarchingCubes::HasRoots( node1->nodeData.mcIndex ) )
2309 hash_map< long long , std::pair< RootInfo , int > >::iterator iter;
2310 int isoTri[DIMENSION*MarchingCubes::MAX_TRIANGLES];
2311 int count=MarchingCubes::AddTriangleIndices( node1->nodeData.mcIndex , isoTri );
2313 for(
int j=0 ; j<count ; j++ )
2314 for(
int k=0 ; k<3 ; k++ )
2315 if( fIndex==Cube::FaceAdjacentToEdges( isoTri[j*3+k] , isoTri[j*3+((k+1)%3)] ) )
2316 if( GetRootIndex( node1 , isoTri[j*3+k] , maxDepth , ri1 ) && GetRootIndex( node1 , isoTri[j*3+((k+1)%3)] , maxDepth , ri2 ) )
2318 long long key1=ri1.key , key2=ri2.key;
2319 edges->push_back( std::pair< RootInfo , RootInfo >( ri2 , ri1 ) );
2320 iter = vertexCount->find( key1 );
2321 if( iter==vertexCount->end() )
2323 (*vertexCount)[key1].first = ri1;
2324 (*vertexCount)[key1].second=0;
2326 iter=vertexCount->find(key2);
2327 if( iter==vertexCount->end() )
2329 (*vertexCount)[key2].first = ri2;
2330 (*vertexCount)[key2].second=0;
2332 (*vertexCount)[key1].second--;
2333 (*vertexCount)[key2].second++;
2335 else fprintf( stderr ,
"Bad Edge 1: %d %d\n" , ri1.key , ri2.key );
2339 template<
int Degree >
2349 bool flags[3][3][3];
2350 int maxDepth = tree.maxDepth();
2353 if( subdivideDepth<=0 ) sDepth = 0;
2354 else sDepth = maxDepth-subdivideDepth;
2355 if( sDepth<=0 )
return;
2359 TreeOctNode::NeighborKey3 nKey;
2360 nKey.set( maxDepth );
2362 if( leaf->depth()>sDepth )
2364 int d , off[3] , _off[3];
2365 leaf->depthAndOffset( d , off );
2366 int res = (1<<d)-1 , _res = ( 1<<(d-sDepth) )-1;
2367 _off[0] = off[0]&_res , _off[1] = off[1]&_res , _off[2] = off[2]&_res;
2368 bool boundary[3][2] =
2370 { (off[0]!=0 && _off[0]==0) , (off[0]!=res && _off[0]==_res) } ,
2371 { (off[1]!=0 && _off[1]==0) , (off[1]!=res && _off[1]==_res) } ,
2372 { (off[2]!=0 && _off[2]==0) , (off[2]!=res && _off[2]==_res) }
2375 if( boundary[0][0] || boundary[0][1] || boundary[1][0] || boundary[1][1] || boundary[2][0] || boundary[2][1] )
2377 TreeOctNode::Neighbors3& neighbors = nKey.getNeighbors( leaf );
2378 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;
2379 int x=0 , y=0 , z=0;
2380 if ( boundary[0][0] && !neighbors.neighbors[0][1][1] ) x = -1;
2381 else if( boundary[0][1] && !neighbors.neighbors[2][1][1] ) x = 1;
2382 if ( boundary[1][0] && !neighbors.neighbors[1][0][1] ) y = -1;
2383 else if( boundary[1][1] && !neighbors.neighbors[1][2][1] ) y = 1;
2384 if ( boundary[2][0] && !neighbors.neighbors[1][1][0] ) z = -1;
2385 else if( boundary[2][1] && !neighbors.neighbors[1][1][2] ) z = 1;
2390 if( x && y && z ) flags[1+x][1+y][1+z] =
true;
2392 if( x && y ) flags[1+x][1+y][1 ] =
true;
2393 if( x && z ) flags[1+x][1 ][1+z] =
true;
2394 if( y && z ) flags[1 ][1+y][1+1] =
true;
2396 if( x ) flags[1+x][1 ][1 ] =
true;
2397 if( y ) flags[1 ][1+y][1 ] =
true;
2398 if( z ) flags[1 ][1 ][1+z] =
true;
2399 nKey.setNeighbors( leaf , flags );
2403 _sNodes.set( tree );
2406 template<
int Degree>
2409 fData.setValueTables( fData.VALUE_FLAG | fData.D_VALUE_FLAG , 0 , postNormalSmooth );
2412 RefineBoundary( subdivideDepth );
2414 RootData rootData , coarseRootData;
2415 std::vector< Point3D< float > >* interiorPoints;
2416 int maxDepth = tree.maxDepth();
2418 int sDepth = subdivideDepth<=0 ? 0 : std::max< int >( 0 , maxDepth-subdivideDepth );
2420 std::vector< Real > metSolution( _sNodes.nodeCount[maxDepth] , 0 );
2421 #pragma omp parallel for num_threads( threads )
2422 for(
int i=_sNodes.nodeCount[_minDepth] ; i<_sNodes.nodeCount[maxDepth] ; i++ ) metSolution[i] = _sNodes.treeNodes[i]->nodeData.solution;
2423 for(
int d=0 ; d<maxDepth ; d++ ) UpSample( d , _sNodes , &metSolution[0] );
2426 #pragma omp parallel for num_threads( threads )
2427 for(
int i=0 ; i<_sNodes.nodeCount[maxDepth+1] ; i++ ) _sNodes.treeNodes[i]->nodeData.mcIndex = 0;
2429 rootData.boundaryValues =
new hash_map<
long long , std::pair<
Real ,
Point3D< Real > > >();
2432 int maxCCount = _sNodes.getMaxCornerCount( &tree , sDepth , maxDepth , threads );
2433 int maxECount = _sNodes.getMaxEdgeCount ( &tree , sDepth , threads );
2434 rootData.cornerValues =
new Real [ maxCCount ];
2436 rootData.interiorRoots =
new int [ maxECount ];
2437 rootData.cornerValuesSet =
new char[ maxCCount ];
2438 rootData.cornerNormalsSet =
new char[ maxCCount ];
2439 rootData.edgesSet =
new char[ maxECount ];
2440 _sNodes.setCornerTable( coarseRootData , &tree , sDepth , threads );
2441 coarseRootData.cornerValues =
new Real[ coarseRootData.cCount ];
2442 coarseRootData.cornerNormals =
new Point3D< Real >[ coarseRootData.cCount ];
2443 coarseRootData.cornerValuesSet =
new char[ coarseRootData.cCount ];
2444 coarseRootData.cornerNormalsSet =
new char[ coarseRootData.cCount ];
2445 memset( coarseRootData.cornerValuesSet , 0 ,
sizeof(
char ) * coarseRootData.cCount );
2446 memset( coarseRootData.cornerNormalsSet , 0 ,
sizeof(
char ) * coarseRootData.cCount );
2449 std::vector< TreeOctNode::ConstNeighborKey3 > nKeys( threads );
2450 for(
int t=0 ; t<threads ; t++ ) nKeys[t].set( maxDepth );
2451 TreeOctNode::ConstNeighborKey3 nKey;
2452 std::vector< TreeOctNode::ConstNeighborKey5 > nKeys5( threads );
2453 for(
int t=0 ; t<threads ; t++ ) nKeys5[t].set( maxDepth );
2454 TreeOctNode::ConstNeighborKey5 nKey5;
2455 nKey5.set( maxDepth );
2456 nKey.set( maxDepth );
2458 for(
int i=_sNodes.nodeCount[sDepth] ; i<_sNodes.nodeCount[sDepth+1] ; i++ )
2460 if( !_sNodes.treeNodes[i]->children )
continue;
2462 _sNodes.setCornerTable( rootData , _sNodes.treeNodes[i] , threads );
2463 _sNodes.setEdgeTable ( rootData , _sNodes.treeNodes[i] , threads );
2464 memset( rootData.cornerValuesSet , 0 ,
sizeof(
char ) * rootData.cCount );
2465 memset( rootData.cornerNormalsSet , 0 ,
sizeof(
char ) * rootData.cCount );
2466 memset( rootData.edgesSet , 0 ,
sizeof(
char ) * rootData.eCount );
2467 interiorPoints =
new std::vector< Point3D< float > >();
2468 for(
int d=maxDepth ; d>sDepth ; d-- )
2470 int leafNodeCount = 0;
2471 std::vector< TreeOctNode* > leafNodes;
2472 for(
TreeOctNode* node=_sNodes.treeNodes[i]->
nextLeaf() ; node ; node=_sNodes.treeNodes[i]->
nextLeaf( node ) )
if( node->d==d ) leafNodeCount++;
2473 leafNodes.reserve( leafNodeCount );
2474 for(
TreeOctNode* node=_sNodes.treeNodes[i]->
nextLeaf() ; node ; node=_sNodes.treeNodes[i]->
nextLeaf( node ) )
if( node->d==d ) leafNodes.push_back( node );
2475 Stencil< double , 3 > stencil1[8] , stencil2[8][8];
2476 SetEvaluationStencils( d , stencil1 , stencil2 );
2479 #pragma omp parallel for num_threads( threads )
2480 for(
int t=0 ; t<threads ; t++ )
for(
int i=(leafNodeCount*t)/threads ; i<(leafNodeCount*(t+1))/threads ; i++ )
2483 SetIsoCorners( isoValue , leaf , rootData , rootData.cornerValuesSet , rootData.cornerValues , nKeys[t] , &metSolution[0] , stencil1 , stencil2 );
2488 int res = 1<<(d-sDepth);
2489 off[0] %= res , off[1] %=res , off[2] %= res;
2491 if( !(off[0]%res) && !(off[1]%res) && !(off[2]%res) )
2494 while( temp->
d!=sDepth ) temp = temp->
parent;
2495 int x = off[0]==0 ? 0 : 1 , y = off[1]==0 ? 0 : 1 , z = off[2]==0 ? 0 : 1;
2496 int c = Cube::CornerIndex( x , y , z );
2497 int idx = coarseRootData.cornerIndices( temp )[ c ];
2498 coarseRootData.cornerValues[ idx ] = rootData.cornerValues[ rootData.cornerIndices( leaf )[c] ];
2499 coarseRootData.cornerValuesSet[ idx ] =
true;
2504 SetMCRootPositions( leaf , sDepth , isoValue , nKeys5[t] , rootData , interiorPoints , mesh , &metSolution[0] , nonLinearFit );
2510 std::vector< Point3D< float > > barycenters;
2511 std::vector< Point3D< float > >* barycenterPtr = addBarycenter ? & barycenters : NULL;
2512 #endif // MISHA_DEBUG
2513 #pragma omp parallel for num_threads( threads )
2514 for(
int t=0 ; t<threads ; t++ )
for(
int i=(leafNodeCount*t)/threads ; i<(leafNodeCount*(t+1))/threads ; i++ )
2519 GetMCIsoTriangles( leaf , mesh , rootData , interiorPoints , offSet , sDepth , polygonMesh , barycenterPtr );
2520 #else // !MISHA_DEBUG
2521 GetMCIsoTriangles( leaf , mesh , rootData , interiorPoints , offSet , sDepth , addBarycenter , polygonMesh );
2522 #endif // MISHA_DEBUG
2525 for(
int i=0 ; i<barycenters.size() ; i++ ) interiorPoints->push_back( barycenters[i] );
2526 #endif // MISHA_DEBUG
2530 delete interiorPoints;
2534 delete[] rootData.cornerValues ,
delete[] rootData.cornerNormals , rootData.cornerValues = NULL , rootData.cornerNormals = NULL;
2535 delete[] rootData.cornerValuesSet ,
delete[] rootData.cornerNormalsSet , rootData.cornerValuesSet = NULL , rootData.cornerNormalsSet = NULL;
2536 delete[] rootData.interiorRoots ; rootData.interiorRoots = NULL;
2537 delete[] rootData.edgesSet ; rootData.edgesSet = NULL;
2538 coarseRootData.interiorRoots = NULL;
2539 coarseRootData.boundaryValues = rootData.boundaryValues;
2540 for( poisson::hash_map< long long , int >::iterator iter=rootData.boundaryRoots.begin() ; iter!=rootData.boundaryRoots.end() ; iter++ )
2541 coarseRootData.boundaryRoots[iter->first] = iter->second;
2543 for(
int d=sDepth ; d>=0 ; d-- )
2545 Stencil< double , 3 > stencil1[8] , stencil2[8][8];
2546 SetEvaluationStencils( d , stencil1 , stencil2 );
2548 std::vector< Point3D< float > > barycenters;
2549 std::vector< Point3D< float > >* barycenterPtr = addBarycenter ? &barycenters : NULL;
2550 #endif // MISHA_DEBUG
2551 for(
int i=_sNodes.nodeCount[d] ; i<_sNodes.nodeCount[d+1] ; i++ )
2557 SetIsoCorners( isoValue , leaf , coarseRootData , coarseRootData.cornerValuesSet , coarseRootData.cornerValues , nKey , &metSolution[0] , stencil1 , stencil2 );
2562 SetMCRootPositions( leaf , 0 , isoValue , nKey5 , coarseRootData , NULL , mesh , &metSolution[0] , nonLinearFit );
2564 GetMCIsoTriangles( leaf , mesh , coarseRootData , NULL , 0 , 0 , polygonMesh , barycenterPtr );
2565 #else // !MISHA_DEBUG
2566 GetMCIsoTriangles( leaf , mesh , coarseRootData , NULL , 0 , 0 , addBarycenter , polygonMesh );
2567 #endif // MISHA_DEBUG
2573 delete[] coarseRootData.cornerValues ,
delete[] coarseRootData.cornerNormals;
2574 delete[] coarseRootData.cornerValuesSet ,
delete[] coarseRootData.cornerNormalsSet;
2575 delete rootData.boundaryValues;
2577 template<
int Degree>
2582 VertexData::CenterIndex(node,fData.
depth,idx);
2583 idx[0]*=fData.functionCount;
2584 idx[1]*=fData.functionCount;
2585 idx[2]*=fData.functionCount;
2586 int minDepth = std::max< int >( 0 , std::min< int >( _minDepth , node->
depth()-1 ) );
2587 for(
int i=minDepth ; i<=node->
depth() ; i++ )
2588 for(
int j=0;j<3;j++)
2589 for(
int k=0;k<3;k++)
2590 for(
int l=0;l<3;l++)
2597 fData.valueTables[idx[0]+
int(n->
off[0])]*
2598 fData.valueTables[idx[1]+
int(n->
off[1])]*
2599 fData.valueTables[idx[2]+
int(n->
off[2])]);
2604 for(
int i=0;i<Cube::CORNERS;i++){
2605 int ii=Cube::AntipodalCornerIndex(i);
2606 const TreeOctNode* n=&node->
children[i];
2609 fData.valueTables[idx[0]+
int(n->off[0])]*
2610 fData.valueTables[idx[1]+
int(n->off[1])]*
2611 fData.valueTables[idx[2]+
int(n->off[2])]);
2612 if( n->children ) n=&n->children[ii];
2619 template<
int Degree >
2620 Real Octree< Degree >::getCornerValue(
const OctNode< TreeNodeData , Real >::ConstNeighborKey3& neighborKey3 ,
const TreeOctNode* node ,
int corner ,
const Real* metSolution )
2625 VertexData::CornerIndex( node , corner , fData.depth , idx );
2626 idx[0] *= fData.functionCount;
2627 idx[1] *= fData.functionCount;
2628 idx[2] *= fData.functionCount;
2630 int d = node->depth();
2632 int startX = 0 , endX = 3 , startY = 0 , endY = 3 , startZ = 0 , endZ = 3;
2633 Cube::FactorCornerIndex( corner , cx , cy , cz );
2635 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d];
2636 if( cx==0 ) endX = 2;
2638 if( cy==0 ) endY = 2;
2640 if( cz==0 ) endZ = 2;
2642 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2644 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2648 fData.valueTables[ idx[0]+int(n->off[0]) ]*
2649 fData.valueTables[ idx[1]+int(n->off[1]) ]*
2650 fData.valueTables[ idx[2]+int(n->off[2]) ];
2655 if( d>0 && d>_minDepth )
2657 int _corner = int( node - node->parent->children );
2658 int _cx , _cy , _cz;
2659 Cube::FactorCornerIndex( _corner , _cx , _cy , _cz );
2660 if( cx!=_cx ) startX = 0 , endX = 3;
2661 if( cy!=_cy ) startY = 0 , endY = 3;
2662 if( cz!=_cz ) startZ = 0 , endZ = 3;
2663 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d-1];
2664 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2666 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2670 fData.valueTables[ idx[0]+int(n->off[0]) ]*
2671 fData.valueTables[ idx[1]+int(n->off[1]) ]*
2672 fData.valueTables[ idx[2]+int(n->off[2]) ];
2673 value += metSolution[ n->nodeData.nodeIndex ] * v;
2677 return Real( value );
2679 template<
int Degree >
2680 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] )
2683 int d = node->depth();
2685 int startX = 0 , endX = 3 , startY = 0 , endY = 3 , startZ = 0 , endZ = 3;
2686 Cube::FactorCornerIndex( corner , cx , cy , cz );
2688 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d];
2689 if( cx==0 ) endX = 2;
2691 if( cy==0 ) endY = 2;
2693 if( cz==0 ) endZ = 2;
2695 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2697 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2701 if( d>0 && d>_minDepth )
2703 int _corner = int( node - node->parent->children );
2704 int _cx , _cy , _cz;
2705 Cube::FactorCornerIndex( _corner , _cx , _cy , _cz );
2706 if( cx!=_cx ) startX = 0 , endX = 3;
2707 if( cy!=_cy ) startY = 0 , endY = 3;
2708 if( cz!=_cz ) startZ = 0 , endZ = 3;
2709 TreeOctNode::ConstNeighbors3& neighbors = neighborKey3.neighbors[d-1];
2710 for(
int x=startX ; x<endX ; x++ )
for(
int y=startY ; y<endY ; y++ )
for(
int z=startZ ; z<endZ ; z++ )
2712 const TreeOctNode* n=neighbors.neighbors[x][y][z];
2713 if( n ) value += metSolution[ n->nodeData.nodeIndex ] * stencil2[x][y][z];
2716 return Real( value );
2718 template<
int Degree >
2719 Point3D< Real > Octree< Degree >::getCornerNormal(
const OctNode< TreeNodeData , Real >::ConstNeighborKey5& neighborKey5 ,
const TreeOctNode* node ,
int corner ,
const Real* metSolution )
2722 Point3D< Real > normal;
2723 normal[0] = normal[1] = normal[2] = 0.;
2725 VertexData::CornerIndex( node , corner , fData.depth , idx );
2726 idx[0] *= fData.functionCount;
2727 idx[1] *= fData.functionCount;
2728 idx[2] *= fData.functionCount;
2730 int d = node->depth();
2733 TreeOctNode::ConstNeighbors5& neighbors = neighborKey5.neighbors[d];
2734 for(
int j=0 ; j<5 ; j++ )
for(
int k=0 ; k<5 ; k++ )
for(
int l=0 ; l<5 ; l++ )
2736 const TreeOctNode* n=neighbors.neighbors[j][k][l];
2739 int _idx[] = { idx[0] + n->
off[0] , idx[1] + n->off[1] , idx[2] + n->off[2] };
2740 double values[] = { fData.valueTables[_idx[0]] , fData.valueTables[_idx[1]] , fData.valueTables[_idx[2]] };
2741 double dValues[] = { fData.dValueTables[_idx[0]] , fData.dValueTables[_idx[1]] , fData.dValueTables[_idx[2]] };
2742 Real solution = n->nodeData.solution;
2743 normal[0] +=
Real( dValues[0] * values[1] * values[2] * solution );
2744 normal[1] +=
Real( values[0] * dValues[1] * values[2] * solution );
2745 normal[2] +=
Real( values[0] * values[1] * dValues[2] * solution );
2749 if( d>0 && d>_minDepth )
2751 TreeOctNode::ConstNeighbors5& neighbors = neighborKey5.neighbors[d-1];
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 = metSolution[ n->nodeData.nodeIndex ];
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 );
2769 template<
int Degree >
2772 Real isoValue , weightSum;
2774 neighborKey2.set( fData.depth );
2775 fData.setValueTables( fData.VALUE_FLAG , 0 );
2777 isoValue = weightSum = 0;
2778 #pragma omp parallel for num_threads( threads ) reduction( + : isoValue , weightSum )
2779 for(
int t=0 ; t<threads ; t++)
2781 TreeOctNode::ConstNeighborKey3 nKey;
2782 nKey.set( _sNodes.maxDepth-1 );
2783 int nodeCount = _sNodes.nodeCount[ _sNodes.maxDepth ];
2784 for(
int i=(nodeCount*t)/threads ; i<(nodeCount*(t+1))/threads ; i++ )
2787 nKey.getNeighbors( temp );
2791 isoValue += getCenterValue( nKey , temp ) * w;
2796 return isoValue/weightSum;
2799 template<
int Degree >
2800 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] )
2802 Real cornerValues[ Cube::CORNERS ];
2808 int mn = 2 , mx = (1<<d)-2;
2809 isInterior = ( off[0]>=mn && off[0]<mx && off[1]>=mn && off[1]<mx && off[2]>=mn && off[2]<mx );
2810 nKey.getNeighbors( leaf );
2811 for(
int c=0 ; c<Cube::CORNERS ; c++ )
2813 int vIndex = cIndices[c];
2814 if( valuesSet[vIndex] ) cornerValues[c] = values[vIndex];
2817 if( isInterior ) cornerValues[c] = getCornerValue( nKey , leaf , c , metSolution , stencil1[c].values , stencil2[
int(leaf - leaf->
parent->
children)][c].values );
2818 else cornerValues[c] = getCornerValue( nKey , leaf , c , metSolution );
2819 values[vIndex] = cornerValues[c];
2820 valuesSet[vIndex] = 1;
2824 leaf->
nodeData.
mcIndex = MarchingCubes::GetIndex( cornerValues , isoValue );
2829 TreeOctNode* parent = leaf->
parent;
2831 int mcid = leaf->
nodeData.
mcIndex & (1<<MarchingCubes::cornerMap()[c]);
2835 poisson::atomicOr(parent->nodeData.mcIndex, mcid);
2839 if( parent->parent && parent->parent->d>=_minDepth && (parent-parent->parent->children)==c )
2841 poisson::atomicOr(parent->parent->nodeData.mcIndex, mcid);
2842 parent = parent->parent;
2851 template<
int Degree>
2852 int Octree<Degree>::InteriorFaceRootCount(
const TreeOctNode* node,
const int &faceIndex,
int maxDepth){
2853 int c1,c2,e1,e2,dir,off,cnt=0;
2854 int corners[Cube::CORNERS/2];
2856 Cube::FaceCorners(faceIndex,corners[0],corners[1],corners[2],corners[3]);
2857 Cube::FactorFaceIndex(faceIndex,dir,off);
2862 e1=Cube::EdgeIndex(1,off,1);
2863 e2=Cube::EdgeIndex(2,off,1);
2866 e1=Cube::EdgeIndex(0,off,1);
2867 e2=Cube::EdgeIndex(2,1,off);
2870 e1=Cube::EdgeIndex(0,1,off);
2871 e2=Cube::EdgeIndex(1,1,off);
2874 cnt+=EdgeRootCount(&node->children[c1],e1,maxDepth)+EdgeRootCount(&node->children[c1],e2,maxDepth);
2877 e1=Cube::EdgeIndex(1,off,0);
2878 e2=Cube::EdgeIndex(2,off,0);
2881 e1=Cube::EdgeIndex(0,off,0);
2882 e2=Cube::EdgeIndex(2,0,off);
2885 e1=Cube::EdgeIndex(0,0,off);
2886 e2=Cube::EdgeIndex(1,0,off);
2889 cnt+=EdgeRootCount(&node->children[c2],e1,maxDepth)+EdgeRootCount(&node->children[c2],e2,maxDepth);
2890 for(
int i=0;i<Cube::CORNERS/2;i++){
if(node->children[corners[i]].children){cnt+=InteriorFaceRootCount(&node->children[corners[i]],faceIndex,maxDepth);}}
2895 template<
int Degree>
2896 int Octree<Degree>::EdgeRootCount(
const TreeOctNode* node,
int edgeIndex,
int maxDepth){
2899 Cube::FacesAdjacentToEdge(edgeIndex,f1,f2);
2904 if(node->depth()<maxDepth){
2906 if(temp && temp->children){
2908 eIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f1);
2911 temp=node->faceNeighbor(f2);
2912 if(temp && temp->children){
2914 eIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f2);
2917 temp=node->edgeNeighbor(edgeIndex);
2918 if(temp && temp->children){
2920 eIndex=Cube::EdgeReflectEdgeIndex(edgeIndex);
2926 Cube::EdgeCorners(eIndex,c1,c2);
2927 if(finest->children)
return EdgeRootCount(&finest->children[c1],eIndex,maxDepth)+EdgeRootCount(&finest->children[c2],eIndex,maxDepth);
2928 else return MarchingCubes::HasEdgeRoots(finest->nodeData.mcIndex,eIndex);
2930 template<
int Degree>
2931 int Octree<Degree>::IsBoundaryFace(
const TreeOctNode* node,
int faceIndex,
int subdivideDepth){
2932 int dir,offset,d,o[3],idx;
2934 if(subdivideDepth<0){
return 0;}
2935 if(node->d<=subdivideDepth){
return 1;}
2936 Cube::FactorFaceIndex(faceIndex,dir,offset);
2937 node->depthAndOffset(d,o);
2939 idx=(int(o[dir])<<1) + (offset<<1);
2940 return !(idx%(2<<(int(node->d)-subdivideDepth)));
2942 template<
int Degree>
2943 int Octree<Degree>::IsBoundaryEdge(
const TreeOctNode* node,
int edgeIndex,
int subdivideDepth){
2945 Cube::FactorEdgeIndex(edgeIndex,dir,x,y);
2946 return IsBoundaryEdge(node,dir,x,y,subdivideDepth);
2948 template<
int Degree>
2949 int Octree<Degree>::IsBoundaryEdge(
const TreeOctNode* node ,
int dir ,
int x ,
int y ,
int subdivideDepth )
2951 int d , o[3] , idx1 , idx2 , mask;
2953 if( subdivideDepth<0 )
return 0;
2954 if( node->d<=subdivideDepth )
return 1;
2955 node->depthAndOffset( d , o );
2972 mask = 1<<( int(node->d) - subdivideDepth );
2973 return !(idx1%(mask)) || !(idx2%(mask));
2975 template<
int Degree >
2982 Cube::FactorEdgeIndex( ri.edgeIndex , o , i1 , i2 );
2983 ri.node->centerAndWidth( c , width );
2987 start[0] = c[0] - width/2;
2988 end [0] = c[0] + width/2;
2989 start[1] = end[1] = c[1] - width/2 + width*i1;
2990 start[2] = end[2] = c[2] - width/2 + width*i2;
2993 start[0] = end[0] = c[0] - width/2 + width*i1;
2994 start[1] = c[1] - width/2;
2995 end [1] = c[1] + width/2;
2996 start[2] = end[2] = c[2] - width/2 + width*i2;
2999 start[0] = end[0] = c[0] - width/2 + width*i1;
3000 start[1] = end[1] = c[1] - width/2 + width*i2;
3001 start[2] = c[2] - width/2;
3002 end [2] = c[2] + width/2;
3009 template<
int Degree >
3010 int Octree< Degree >::GetRoot(
const RootInfo& ri , Real isoValue , TreeOctNode::ConstNeighborKey5& neighborKey5 , Point3D< Real > & position , RootData& rootData ,
int sDepth ,
const Real* metSolution ,
int nonLinearFit )
3012 if( !MarchingCubes::HasRoots( ri.node->nodeData.mcIndex ) )
return 0;
3014 Cube::EdgeCorners( ri.edgeIndex , c1 , c2 );
3015 if( !MarchingCubes::HasEdgeRoots( ri.node->nodeData.mcIndex , ri.edgeIndex ) )
return 0;
3017 long long key1 , key2;
3018 Point3D< Real > n[2];
3020 int i , o , i1 , i2 , rCount=0;
3022 std::vector< double > roots;
3024 Real center , width;
3026 Cube::FactorEdgeIndex( ri.edgeIndex , o , i1 , i2 );
3027 int idx1[3] , idx2[3];
3028 key1 = VertexData::CornerIndex( ri.node , c1 , fData.depth , idx1 );
3029 key2 = VertexData::CornerIndex( ri.node , c2 , fData.depth , idx2 );
3031 bool isBoundary = ( IsBoundaryEdge( ri.node , ri.edgeIndex , sDepth )!=0 );
3032 bool haveKey1 , haveKey2;
3033 std::pair< Real , Point3D< Real > > keyValue1 , keyValue2;
3036 iter1 = rootData.cornerIndices( ri.node )[ c1 ];
3037 iter2 = rootData.cornerIndices( ri.node )[ c2 ];
3038 keyValue1.first = rootData.cornerValues[iter1];
3039 keyValue2.first = rootData.cornerValues[iter2];
3042 #pragma omp critical (normal_hash_access)
3044 haveKey1 = ( rootData.boundaryValues->find( key1 )!=rootData.boundaryValues->end() );
3045 haveKey2 = ( rootData.boundaryValues->find( key2 )!=rootData.boundaryValues->end() );
3046 if( haveKey1 ) keyValue1 = (*rootData.boundaryValues)[key1];
3047 if( haveKey2 ) keyValue2 = (*rootData.boundaryValues)[key2];
3052 haveKey1 = ( rootData.cornerNormalsSet[ iter1 ] != 0 );
3053 haveKey2 = ( rootData.cornerNormalsSet[ iter2 ] != 0 );
3054 keyValue1.first = rootData.cornerValues[iter1];
3055 keyValue2.first = rootData.cornerValues[iter2];
3056 if( haveKey1 ) keyValue1.second = rootData.cornerNormals[iter1];
3057 if( haveKey2 ) keyValue2.second = rootData.cornerNormals[iter2];
3060 if( !haveKey1 || !haveKey2 ) neighborKey5.getNeighbors( ri.node );
3061 if( !haveKey1 ) keyValue1.second = getCornerNormal( neighborKey5 , ri.node , c1 , metSolution );
3062 x0 = keyValue1.first;
3063 n[0] = keyValue1.second;
3065 if( !haveKey2 ) keyValue2.second = getCornerNormal( neighborKey5 , ri.node , c2 , metSolution );
3066 x1 = keyValue2.first;
3067 n[1] = keyValue2.second;
3069 if( !haveKey1 || !haveKey2 )
3073 #pragma omp critical (normal_hash_access)
3075 if( !haveKey1 ) (*rootData.boundaryValues)[key1] = keyValue1;
3076 if( !haveKey2 ) (*rootData.boundaryValues)[key2] = keyValue2;
3081 if( !haveKey1 ) rootData.cornerNormals[iter1] = keyValue1.second , rootData.cornerNormalsSet[ iter1 ] = 1;
3082 if( !haveKey2 ) rootData.cornerNormals[iter2] = keyValue2.second , rootData.cornerNormalsSet[ iter2 ] = 1;
3087 ri.node->centerAndWidth(c,width);
3089 for( i=0 ; i<DIMENSION ; i++ ) n[0][i] *= width , n[1][i] *= width;
3094 position[1] = c[1]-width/2+width*i1;
3095 position[2] = c[2]-width/2+width*i2;
3098 position[0] = c[0]-width/2+width*i1;
3099 position[2] = c[2]-width/2+width*i2;
3102 position[0] = c[0]-width/2+width*i1;
3103 position[1] = c[1]-width/2+width*i2;
3111 double scl=(x1-x0)/((dx1+dx0)/2);
3116 P.coefficients[0] = x0;
3117 P.coefficients[1] = dx0;
3118 P.coefficients[2] = 3*(x1-x0)-dx1-2*dx0;
3120 P.getSolutions( isoValue , roots , EPSILON );
3121 for( i=0 ; i<int(roots.size()) ; i++ )
3122 if( roots[i]>=0 && roots[i]<=1 )
3124 averageRoot +=
Real( roots[i] );
3127 if( rCount && nonLinearFit ) averageRoot /= rCount;
3128 else averageRoot =
Real((x0-isoValue)/(x0-x1));
3129 if( averageRoot<0 || averageRoot>1 )
3131 fprintf( stderr ,
"[WARNING] Bad average root: %f\n" , averageRoot );
3132 fprintf( stderr ,
"\t(%f %f) , (%f %f) (%f)\n" , x0 , x1 , dx0 , dx1 , isoValue );
3133 if( averageRoot<0 ) averageRoot = 0;
3134 if( averageRoot>1 ) averageRoot = 1;
3136 position[o] =
Real(center-width/2+width*averageRoot);
3139 template<
int Degree >
3140 int Octree< Degree >::GetRootIndex(
const TreeOctNode* node ,
int edgeIndex ,
int maxDepth ,
int sDepth,RootInfo& ri )
3146 Cube::FacesAdjacentToEdge(edgeIndex,f1,f2);
3149 finestIndex=edgeIndex;
3150 if(node->depth()<maxDepth){
3151 if(IsBoundaryFace(node,f1,sDepth)){temp=NULL;}
3152 else{temp=node->faceNeighbor(f1);}
3153 if(temp && temp->children){
3155 finestIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f1);
3158 if(IsBoundaryFace(node,f2,sDepth)){temp=NULL;}
3159 else{temp=node->faceNeighbor(f2);}
3160 if(temp && temp->children){
3162 finestIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f2);
3165 if(IsBoundaryEdge(node,edgeIndex,sDepth)){temp=NULL;}
3166 else{temp=node->edgeNeighbor(edgeIndex);}
3167 if(temp && temp->children){
3169 finestIndex=Cube::EdgeReflectEdgeIndex(edgeIndex);
3175 Cube::EdgeCorners(finestIndex,c1,c2);
3176 if(finest->children){
3177 if (GetRootIndex(&finest->children[c1],finestIndex,maxDepth,sDepth,ri)) {
return 1;}
3178 else if (GetRootIndex(&finest->children[c2],finestIndex,maxDepth,sDepth,ri)) {
return 1;}
3181 fprintf( stderr ,
"[WARNING] Couldn't find root index with either child\n" );
3187 if( !(MarchingCubes::edgeMask()[finest->nodeData.mcIndex] & (1<<finestIndex)) )
3189 fprintf( stderr ,
"[WARNING] Finest node does not have iso-edge\n" );
3194 Cube::FactorEdgeIndex(finestIndex,o,i1,i2);
3196 finest->depthAndOffset(d,off);
3198 ri.edgeIndex=finestIndex;
3199 int eIndex[2],offset;
3200 offset=BinaryNode<Real>::Index( d , off[o] );
3204 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i1);
3205 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2);
3208 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1);
3209 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2);
3212 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1);
3213 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i2);
3216 ri.key = (
long long)(o) | (
long long)(eIndex[0])<<5 | (
long long)(eIndex[1])<<25 | (
long long)(offset)<<45;
3220 template<
int Degree>
3221 int Octree<Degree>::GetRootIndex(
const TreeOctNode* node ,
int edgeIndex ,
int maxDepth , RootInfo& ri )
3229 if(!(MarchingCubes::edgeMask()[node->nodeData.mcIndex] & (1<<edgeIndex))){
return 0;}
3231 Cube::FacesAdjacentToEdge(edgeIndex,f1,f2);
3234 finestIndex = edgeIndex;
3235 if( node->depth()<maxDepth && !node->children )
3237 temp=node->faceNeighbor( f1 );
3238 if( temp && temp->children ) finest = temp , finestIndex = Cube::FaceReflectEdgeIndex( edgeIndex , f1 );
3241 temp = node->faceNeighbor( f2 );
3242 if( temp && temp->children ) finest = temp , finestIndex = Cube::FaceReflectEdgeIndex( edgeIndex , f2 );
3245 temp = node->edgeNeighbor( edgeIndex );
3246 if( temp && temp->children ) finest = temp , finestIndex = Cube::EdgeReflectEdgeIndex( edgeIndex );
3251 Cube::EdgeCorners( finestIndex , c1 , c2 );
3252 if( finest->children )
3254 if ( GetRootIndex( finest->children + c1 , finestIndex , maxDepth , ri ) )
return 1;
3255 else if( GetRootIndex( finest->children + c2 , finestIndex , maxDepth , ri ) )
return 1;
3258 int d1 , off1[3] , d2 , off2[3];
3259 node->depthAndOffset( d1 , off1 );
3260 finest->depthAndOffset( d2 , off2 );
3261 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 );
3263 for(
int i=0 ; i<8 ; i++ )
if( node->nodeData.mcIndex & (1<<i) ) printf(
"1" );
else printf(
"0" );
3265 for(
int i=0 ; i<8 ; i++ )
if( finest->nodeData.mcIndex & (1<<i) ) printf(
"1" );
else printf(
"0" );
3273 Cube::FactorEdgeIndex(finestIndex,o,i1,i2);
3275 finest->depthAndOffset(d,off);
3277 ri.edgeIndex=finestIndex;
3278 int offset,eIndex[2];
3279 offset = BinaryNode< Real >::CenterIndex( d , off[o] );
3282 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i1);
3283 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2);
3286 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1);
3287 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2);
3290 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1);
3291 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i2);
3294 ri.key= (
long long)(o) | (
long long)(eIndex[0])<<5 | (
long long)(eIndex[1])<<25 | (
long long)(offset)<<45;
3298 template<
int Degree>
3299 int Octree<Degree>::GetRootPair(
const RootInfo& ri,
int maxDepth,RootInfo& pair){
3302 Cube::EdgeCorners(ri.edgeIndex,c1,c2);
3303 while(node->parent){
3304 c=int(node-node->parent->children);
3305 if(c!=c1 && c!=c2){
return 0;}
3306 if(!MarchingCubes::HasEdgeRoots(node->parent->nodeData.mcIndex,ri.edgeIndex)){
3307 if(c==c1){
return GetRootIndex(&node->parent->children[c2],ri.edgeIndex,maxDepth,pair);}
3308 else{
return GetRootIndex(&node->parent->children[c1],ri.edgeIndex,maxDepth,pair);}
3315 template<
int Degree>
3316 int Octree< Degree >::GetRootIndex(
const RootInfo& ri , RootData& rootData , CoredPointIndex& index )
3318 long long key = ri.key;
3319 hash_map< long long , int >::iterator rootIter;
3320 rootIter = rootData.boundaryRoots.find( key );
3321 if( rootIter!=rootData.boundaryRoots.end() )
3324 index.index = rootIter->second;
3327 else if( rootData.interiorRoots )
3329 int eIndex = rootData.edgeIndices( ri.node )[ ri.edgeIndex ];
3330 if( rootData.edgesSet[ eIndex ] )
3333 index.index = rootData.interiorRoots[ eIndex ];
3339 template<
int Degree >
3340 int Octree< Degree >::SetMCRootPositions( TreeOctNode* node ,
int sDepth , Real isoValue , TreeOctNode::ConstNeighborKey5& neighborKey5 , RootData& rootData ,
3341 std::vector<
Point3D< float > >* interiorPositions , CoredMeshData* mesh ,
const Real* metSolution ,
int nonLinearFit )
3343 Point3D< Real > position;
3347 if( !MarchingCubes::HasRoots( node->nodeData.mcIndex ) )
return 0;
3348 for(
int i=0 ; i<DIMENSION ; i++ )
for(
int j=0 ; j<2 ; j++ )
for(
int k=0 ; k<2 ; k++ )
3351 eIndex = Cube::EdgeIndex( i , j , k );
3352 if( GetRootIndex( node , eIndex , fData.depth , ri ) )
3355 if( !rootData.interiorRoots || IsBoundaryEdge( node , i , j , k , sDepth ) )
3357 poisson::hash_map< long long , int >::iterator iter , end;
3359 #pragma omp critical (boundary_roots_hash_access)
3361 iter = rootData.boundaryRoots.find( key );
3362 end = rootData.boundaryRoots.end();
3367 GetRoot( ri , isoValue , neighborKey5 , position , rootData , sDepth , metSolution , nonLinearFit );
3369 #pragma omp critical (boundary_roots_hash_access)
3371 iter = rootData.boundaryRoots.find( key );
3372 end = rootData.boundaryRoots.end();
3375 mesh->inCorePoints.push_back( position );
3376 rootData.boundaryRoots[key] = int( mesh->inCorePoints.size() ) - 1;
3379 if( iter==end ) count++;
3384 int nodeEdgeIndex = rootData.edgeIndices( ri.node )[ ri.edgeIndex ];
3385 if( !rootData.edgesSet[ nodeEdgeIndex ] )
3388 GetRoot( ri , isoValue , neighborKey5 , position , rootData , sDepth , metSolution , nonLinearFit );
3390 #pragma omp critical (add_point_access)
3392 if( !rootData.edgesSet[ nodeEdgeIndex ] )
3394 rootData.interiorRoots[ nodeEdgeIndex ] = mesh->addOutOfCorePoint( position );
3395 interiorPositions->push_back( position );
3396 rootData.edgesSet[ nodeEdgeIndex ] = 1;
3406 template<
int Degree>
3407 int Octree< Degree >::SetBoundaryMCRootPositions(
int sDepth , Real isoValue , RootData& rootData , CoredMeshData* mesh ,
int nonLinearFit )
3409 Point3D< Real > position;
3410 int i,j,k,eIndex,hits=0;
3418 if( MarchingCubes::HasRoots( node->nodeData.mcIndex ) )
3421 for( i=0 ; i<DIMENSION ; i++ )
3422 for( j=0 ; j<2 ; j++ )
3423 for( k=0 ; k<2 ; k++ )
3424 if( IsBoundaryEdge( node , i , j , k , sDepth ) )
3428 eIndex = Cube::EdgeIndex( i , j , k );
3429 if( GetRootIndex( node , eIndex , fData.depth , ri ) )
3432 if( rootData.boundaryRoots.find(key)==rootData.boundaryRoots.end() )
3434 GetRoot( ri , isoValue , position , rootData , sDepth , nonLinearFit );
3435 mesh->inCorePoints.push_back( position );
3436 rootData.boundaryRoots[key] = int( mesh->inCorePoints.size() )-1;
3442 if( hits ) node=tree.nextLeaf(node);
3443 else node=tree.nextBranch(node);
3447 template<
int Degree>
3448 void Octree< Degree >::GetMCIsoEdges( TreeOctNode* node ,
int sDepth , std::vector< std::pair< RootInfo , RootInfo > >& edges )
3451 int count=0 , tris=0;
3452 int isoTri[ DIMENSION * MarchingCubes::MAX_TRIANGLES ];
3453 FaceEdgesFunction fef;
3455 hash_map< long long , std::pair< RootInfo , int > >::iterator iter;
3456 hash_map< long long , std::pair< RootInfo , int > > vertexCount;
3460 fef.vertexCount = &vertexCount;
3461 count = MarchingCubes::AddTriangleIndices( node->nodeData.mcIndex , isoTri );
3462 for( fIndex=0 ; fIndex<Cube::NEIGHBORS ; fIndex++ )
3464 ref = Cube::FaceReflectFaceIndex( fIndex , fIndex );
3466 temp = node->faceNeighbor( fIndex );
3469 if( temp && temp->children && !IsBoundaryFace( node , fIndex , sDepth ) ) temp->processNodeFaces( temp , &fef , ref );
3474 for(
int j=0 ; j<count ; j++ )
3475 for(
int k=0 ; k<3 ; k++ )
3476 if( fIndex==Cube::FaceAdjacentToEdges( isoTri[j*3+k] , isoTri[j*3+((k+1)%3)] ) )
3477 if( GetRootIndex( node , isoTri[j*3+k] , fData.depth , ri1 ) && GetRootIndex( node , isoTri[j*3+((k+1)%3)] , fData.depth , ri2 ) )
3479 long long key1 = ri1.key , key2 = ri2.key;
3480 edges.push_back( std::pair< RootInfo , RootInfo >( ri1 , ri2 ) );
3481 iter=vertexCount.find( key1 );
3482 if( iter==vertexCount.end() )
3484 vertexCount[key1].first = ri1;
3485 vertexCount[key1].second = 0;
3487 iter=vertexCount.find( key2 );
3488 if( iter==vertexCount.end() )
3490 vertexCount[key2].first = ri2;
3491 vertexCount[key2].second = 0;
3493 vertexCount[key1].second++;
3494 vertexCount[key2].second--;
3498 int r1 = MarchingCubes::HasEdgeRoots( node->nodeData.mcIndex , isoTri[j*3+k] );
3499 int r2 = MarchingCubes::HasEdgeRoots( node->nodeData.mcIndex , isoTri[j*3+((k+1)%3)] );
3500 fprintf( stderr ,
"Bad Edge 2: %d %d\t%d %d\n" , ri1.key , ri2.key , r1 , r2 );
3504 for(
int i=0 ; i<int(edges.size()) ; i++ )
3506 iter = vertexCount.find( edges[i].first.key );
3507 if( iter==vertexCount.end() ) printf(
"Could not find vertex: %lld\n" , edges[i].first );
3508 else if( vertexCount[ edges[i].first.key ].second )
3511 GetRootPair( vertexCount[edges[i].first.key].first , fData.depth , ri );
3512 long long key = ri.key;
3513 iter = vertexCount.find( key );
3514 if( iter==vertexCount.end() )
3517 node->depthAndOffset( d , off );
3518 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] );
3522 edges.push_back( std::pair< RootInfo , RootInfo >( ri , edges[i].first ) );
3523 vertexCount[ key ].second++;
3524 vertexCount[ edges[i].first.key ].second--;
3528 iter = vertexCount.find( edges[i].second.key );
3529 if( iter==vertexCount.end() ) printf(
"Could not find vertex: %lld\n" , edges[i].second );
3530 else if( vertexCount[edges[i].second.key].second )
3533 GetRootPair( vertexCount[edges[i].second.key].first , fData.depth , ri );
3534 long long key = ri.key;
3535 iter=vertexCount.find( key );
3536 if( iter==vertexCount.end() )
3539 node->depthAndOffset( d , off );
3540 printf(
"Vertex pair not in list 2\t[%d] (%d %d %d)\n" , d , off[0] , off[1] , off[2] );
3544 edges.push_back( std::pair< RootInfo , RootInfo >( edges[i].second , ri ) );
3545 vertexCount[key].second--;
3546 vertexCount[ edges[i].second.key ].second++;
3551 template<
int Degree>
3553 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 )
3554 #else // !MISHA_DEBUG
3555 int Octree< Degree >::GetMCIsoTriangles( TreeOctNode* node , CoredMeshData* mesh , RootData& rootData , std::vector<
Point3D< float > >* interiorPositions ,
int offSet ,
int sDepth ,
bool addBarycenter ,
bool polygonMesh )
3556 #endif // MISHA_DEBUG
3559 std::vector< std::pair< RootInfo , RootInfo > > edges;
3560 std::vector< std::vector< std::pair< RootInfo , RootInfo > > > edgeLoops;
3561 GetMCIsoEdges( node , sDepth , edges );
3563 GetEdgeLoops( edges , edgeLoops );
3564 for(
int i=0 ; i<int(edgeLoops.size()) ; i++ )
3567 std::vector<CoredPointIndex> edgeIndices;
3568 for(
int j=0 ; j<int(edgeLoops[i].size()) ; j++ )
3570 if( !GetRootIndex( edgeLoops[i][j].first , rootData , p ) ) printf(
"Bad Point Index\n" );
3571 else edgeIndices.push_back( p );
3574 tris += AddTriangles( mesh , edgeIndices , interiorPositions , offSet , polygonMesh , barycenters );
3575 #else // !MISHA_DEBUG
3576 tris += AddTriangles( mesh , edgeIndices , interiorPositions , offSet , addBarycenter , polygonMesh );
3577 #endif // MISHA_DEBUG
3582 template<
int Degree >
3583 int Octree< Degree >::GetEdgeLoops( std::vector< std::pair< RootInfo , RootInfo > >& edges , std::vector< std::vector< std::pair< RootInfo , RootInfo > > >& loops )
3586 long long frontIdx , backIdx;
3587 std::pair< RootInfo , RootInfo > e , temp;
3590 while( edges.size() )
3592 std::vector< std::pair< RootInfo , RootInfo > > front , back;
3594 loops.resize( loopSize+1 );
3595 edges[0] = edges.back();
3597 frontIdx = e.second.key;
3598 backIdx = e.first.key;
3599 for(
int j=
int(edges.size())-1 ; j>=0 ; j-- )
3601 if( edges[j].first.key==frontIdx || edges[j].second.key==frontIdx )
3603 if( edges[j].first.key==frontIdx ) temp = edges[j];
3604 else temp.first = edges[j].second , temp.second = edges[j].first;
3605 frontIdx = temp.second.key;
3606 front.push_back(temp);
3607 edges[j] = edges.back();
3609 j = int(edges.size());
3611 else if( edges[j].first.key==backIdx || edges[j].second.key==backIdx )
3613 if( edges[j].second.key==backIdx ) temp = edges[j];
3614 else temp.first = edges[j].second , temp.second = edges[j].first;
3615 backIdx = temp.first.key;
3616 back.push_back(temp);
3617 edges[j] = edges.back();
3619 j = int(edges.size());
3622 for(
int j=
int(back.size())-1 ; j>=0 ; j-- ) loops[loopSize].push_back( back[j] );
3623 loops[loopSize].push_back(e);
3624 for(
int j=0 ; j<int(front.size()) ; j++ ) loops[loopSize].push_back( front[j] );
3627 return int(loops.size());
3629 template<
int Degree>
3631 int Octree<Degree>::AddTriangles( CoredMeshData* mesh , std::vector<CoredPointIndex>& edges , std::vector<Point3D<float> >* interiorPositions ,
int offSet ,
bool polygonMesh , std::vector<
Point3D< float > >* barycenters )
3632 #else // !MISHA_DEBUG
3633 int Octree<Degree>::AddTriangles( CoredMeshData* mesh , std::vector<CoredPointIndex>& edges , std::vector<Point3D<float> >* interiorPositions ,
int offSet ,
bool addBarycenter ,
bool polygonMesh )
3634 #endif // MISHA_DEBUG
3636 MinimalAreaTriangulation< float > MAT;
3637 std::vector< Point3D< float > > vertices;
3638 std::vector< TriangleIndex > triangles;
3641 std::vector< CoredVertexIndex > vertices( edges.size() );
3642 for(
int i=0 ; i<int(edges.size()) ; i++ )
3644 vertices[i].idx = edges[i].index;
3645 vertices[i].inCore = (edges[i].inCore!=0);
3647 #pragma omp critical (add_polygon_access)
3649 mesh->addPolygon( vertices );
3653 if( edges.size()>3 )
3655 bool isCoplanar =
false;
3659 #else // !MISHA_DEBUG
3661 #endif // MISHA_DEBUG
3662 for(
int i=0 ; i<int(edges.size()) ; i++ )
3663 for(
int j=0 ; j<i ; j++ )
3664 if( (i+1)%edges.size()!=j && (j+1)%edges.size()!=i )
3666 Point3D< Real > v1 , v2;
3667 if( edges[i].inCore ) v1 = mesh->inCorePoints[ edges[i].index ];
3668 else v1 = (*interiorPositions)[ edges[i].index-offSet ];
3669 if( edges[j].inCore ) v2 = mesh->inCorePoints[ edges[j].index ];
3670 else v2 = (*interiorPositions)[ edges[j].index-offSet ];
3671 for(
int k=0 ; k<3 ; k++ )
if( v1[k]==v2[k] ) isCoplanar =
true;
3676 c[0] = c[1] = c[2] = 0;
3677 for(
int i=0 ; i<int(edges.size()) ; i++ )
3680 if(edges[i].inCore) p = mesh->inCorePoints[edges[i].index ];
3681 else p = (*interiorPositions)[edges[i].index-offSet];
3684 c /=
Real( edges.size() );
3686 #pragma omp critical (add_point_access)
3688 cIdx = mesh->addOutOfCorePoint( c );
3690 barycenters->push_back( c );
3691 #else // !MISHA_DEBUG
3692 interiorPositions->push_back( c );
3693 #endif // MISHA_DEBUG
3695 for(
int i=0 ; i<int(edges.size()) ; i++ )
3697 std::vector< CoredVertexIndex > vertices( 3 );
3698 vertices[0].idx = edges[i ].index;
3699 vertices[1].idx = edges[(i+1)%edges.size()].index;
3700 vertices[2].idx = cIdx;
3701 vertices[0].inCore = (edges[i ].inCore!=0);
3702 vertices[1].inCore = (edges[(i+1)%edges.size()].inCore!=0);
3703 vertices[2].inCore = 0;
3704 #pragma omp critical (add_polygon_access)
3706 mesh->addPolygon( vertices );
3709 return int( edges.size() );
3713 vertices.resize( edges.size() );
3715 for(
int i=0 ; i<int(edges.size()) ; i++ )
3718 if( edges[i].inCore ) p = mesh->inCorePoints[edges[i].index ];
3719 else p = (*interiorPositions)[edges[i].index-offSet];
3722 MAT.GetTriangulation( vertices , triangles );
3723 for(
int i=0 ; i<int(triangles.size()) ; i++ )
3725 std::vector< CoredVertexIndex > _vertices( 3 );
3726 for(
int j=0 ; j<3 ; j++ )
3728 _vertices[j].idx = edges[ triangles[i].idx[j] ].index;
3729 _vertices[j].inCore = (edges[ triangles[i].idx[j] ].inCore!=0);
3731 #pragma omp critical (add_polygon_access)
3733 mesh->addPolygon( _vertices );
3738 else if( edges.size()==3 )
3740 std::vector< CoredVertexIndex > vertices( 3 );
3741 for(
int i=0 ; i<3 ; i++ )
3743 vertices[i].idx = edges[i].index;
3744 vertices[i].inCore = (edges[i].inCore!=0);
3746 #pragma omp critical (add_polygon_access)
3747 mesh->addPolygon( vertices );
3749 return int(edges.size())-2;
3751 template<
int Degree >
3754 if( depth<=0 || depth>tree.maxDepth() ) depth = tree.maxDepth();
3759 Real* values =
new float[ res * res * res ];
3760 memset( values , 0 ,
sizeof(
float ) * res * res * res );
3764 if( n->d>depth )
continue;
3765 if( n->d<_minDepth )
continue;
3766 int d , idx[3] , start[3] , end[3];
3767 n->depthAndOffset( d , idx );
3768 for(
int i=0 ; i<3 ; i++ )
3775 if( !(start[i]&1) ) start[i]++;
3776 if( !( end[i]&1) ) end[i]--;
3778 Real coefficient = n->nodeData.solution;
3779 for(
int x=start[0] ; x<=end[0] ; x+=2 )
3780 for(
int y=start[1] ; y<=end[1] ; y+=2 )
3781 for(
int z=start[2] ; z<=end[2] ; z+=2 )
3783 int xx = (x-1)>>1 , yy = (y-1)>>1 , zz = (z-1)>>1;
3784 values[ zz*res*res + yy*res + xx ] +=
3791 for(
int i=0 ; i<res*res*res ; i++ ) values[i] -= isoValue;
3795 template<
int Degree >
3798 if( depth<=0 || depth>tree.maxDepth() ) depth = tree.maxDepth();
3799 res = 1<<tree.maxDepth();
3800 Real* values =
new float[ res * res * res ];
3801 memset( values , 0 ,
sizeof(
float ) * res * res * res );
3805 if( n->d>depth )
continue;
3806 int d , idx[3] , start[3] , end[3];
3807 n->depthAndOffset( d , idx );
3808 for(
int i=0 ; i<3 ; i++ )
3813 fData.setSampleSpan( idx[i] , start[i] , end[i] );
3815 if( !(start[i]&1) ) start[i]++;
3816 if( !( end[i]&1) ) end[i]--;
3818 for(
int x=start[0] ; x<=end[0] ; x+=2 )
3819 for(
int y=start[1] ; y<=end[1] ; y+=2 )
3820 for(
int z=start[2] ; z<=end[2] ; z+=2 )
3822 int xx = (x-1)>>1 , yy = (y-1)>>1 , zz = (z-1)>>1;
3823 values[ zz*res*res + yy*res + xx ] +=
3824 n->nodeData.centerWeightContribution *
3825 fData.valueTables[ idx[0]+x*fData.functionCount ] *
3826 fData.valueTables[ idx[1]+y*fData.functionCount ] *
3827 fData.valueTables[ idx[2]+z*fData.functionCount ];
3836 long long VertexData::CenterIndex(
const TreeOctNode* node,
int maxDepth){
3838 return CenterIndex(node,maxDepth,idx);
3840 long long VertexData::CenterIndex(
const TreeOctNode* node,
int maxDepth,
int idx[DIMENSION]){
3844 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3846 long long VertexData::CenterIndex(
int depth,
const int offSet[DIMENSION],
int maxDepth,
int idx[DIMENSION]){
3848 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3850 long long VertexData::CornerIndex(
const TreeOctNode* node,
int cIndex,
int maxDepth){
3852 return CornerIndex(node,cIndex,maxDepth,idx);
3854 long long VertexData::CornerIndex(
const TreeOctNode* node ,
int cIndex ,
int maxDepth ,
int idx[DIMENSION] )
3857 Cube::FactorCornerIndex( cIndex , x[0] , x[1] , x[2] );
3860 for(
int i=0 ; i<DIMENSION ; i++ ) idx[i] = BinaryNode<Real>::CornerIndex( maxDepth+1 , d , o[i] , x[i] );
3861 return CornerIndexKey( idx );
3863 long long VertexData::CornerIndex(
int depth ,
const int offSet[DIMENSION] ,
int cIndex ,
int maxDepth ,
int idx[DIMENSION] )
3866 Cube::FactorCornerIndex( cIndex , x[0] , x[1] , x[2] );
3867 for(
int i=0 ; i<DIMENSION ; i++ ) idx[i] = BinaryNode<Real>::CornerIndex( maxDepth+1 , depth , offSet[i] , x[i] );
3868 return CornerIndexKey( idx );
3870 long long VertexData::CornerIndexKey(
const int idx[DIMENSION] )
3872 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3874 long long VertexData::FaceIndex(
const TreeOctNode* node,
int fIndex,
int maxDepth){
3876 return FaceIndex(node,fIndex,maxDepth,idx);
3878 long long VertexData::FaceIndex(
const TreeOctNode* node,
int fIndex,
int maxDepth,
int idx[DIMENSION]){
3880 Cube::FactorFaceIndex(fIndex,dir,offset);
3885 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
3887 long long VertexData::EdgeIndex(
const TreeOctNode* node,
int eIndex,
int maxDepth){
3889 return EdgeIndex(node,eIndex,maxDepth,idx);
3891 long long VertexData::EdgeIndex(
const TreeOctNode* node,
int eIndex,
int maxDepth,
int idx[DIMENSION]){
3896 Cube::FactorEdgeIndex(eIndex,o,i1,i2);
3911 return (
long long)(idx[0]) | (
long long)(idx[1])<<15 | (
long long)(idx[2])<<30;
const OctNode * nextLeaf(const OctNode *currentLeaf=NULL) const
void SetRowSize(int row, int count)
std::vector< int > offsets
OctNode * faceNeighbor(int faceIndex, int forceChildren=0)
std::vector< EdgeIndices > eTable
double Length(const Point3D< Real > &p)
UpSampleData(int s, double v1, double v2)
void atomicOr(volatile int &dest, int value)
ConstNeighbors3 * neighbors
const Real MATRIX_ENTRY_EPSILON
static const int VALUE_FLAG
std::vector< CornerIndices > cTable
const OctNode * nextNode(const OctNode *currentNode=NULL) const
void set(int maxDepth, bool useDotRatios=true, bool reflectBoundary=false)
void centerAndWidth(Point3D< Real > ¢er, Real &width) const
const OctNode * neighbors[3][3][3]
void setSampleSpan(int idx, int &start, int &end, double smooth=0) const
Real centerWeightContribution
virtual void setValueTables(int flags, double smooth=0)
void setFullDepth(int maxDepth)
virtual int outOfCorePointCount(void)=0
pcl::poisson::OctNode< class TreeNodeData, Real > TreeOctNode
void depthAndOffset(int &depth, int offset[DIMENSION]) const
std::vector< int > offsets