1/* 2Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho 3All rights reserved. 4 5Redistribution and use in source and binary forms, with or without modification, 6are permitted provided that the following conditions are met: 7 8Redistributions of source code must retain the above copyright notice, this list of 9conditions and the following disclaimer. Redistributions in binary form must reproduce 10the above copyright notice, this list of conditions and the following disclaimer 11in the documentation and/or other materials provided with the distribution. 12 13Neither the name of the Johns Hopkins University nor the names of its contributors 14may be used to endorse or promote products derived from this software without specific 15prior written permission. 16 17THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 18EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES 19OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT 20SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 21INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 22TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 23BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 24CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 25ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 26DAMAGE. 27*/ 28 29#include "Octree.h" 30 31#define ITERATION_POWER 1.0/3 32#define MEMORY_ALLOCATOR_BLOCK_SIZE 1<<12 33 34#define READ_SIZE 1024 35 36#define PAD_SIZE (Real(1.0)) 37 38const Real EPSILON=Real(1e-6); 39const Real ROUND_EPS=Real(1e-5); 40 41 42///////////////////// 43// SortedTreeNodes // 44///////////////////// 45SortedTreeNodes::SortedTreeNodes(void){ 46 nodeCount=NULL; 47 treeNodes=NULL; 48 maxDepth=0; 49} 50SortedTreeNodes::~SortedTreeNodes(void){ 51 if(nodeCount){delete[] nodeCount;} 52 if(treeNodes){delete[] treeNodes;} 53 nodeCount=NULL; 54 treeNodes=NULL; 55} 56void SortedTreeNodes::set(TreeOctNode& root,const int& setIndex){ 57 if(nodeCount){delete[] nodeCount;} 58 if(treeNodes){delete[] treeNodes;} 59 maxDepth=root.maxDepth()+1; 60 nodeCount=new int[maxDepth+1]; 61 treeNodes=new TreeOctNode*[root.nodes()]; 62 63 TreeOctNode* temp=root.nextNode(); 64 int i,cnt=0; 65 while(temp){ 66 treeNodes[cnt++]=temp; 67 temp=root.nextNode(temp); 68 } 69 qsort(treeNodes,cnt,sizeof(const TreeOctNode*),TreeOctNode::CompareForwardPointerDepths); 70 for(i=0;i<=maxDepth;i++){nodeCount[i]=0;} 71 for(i=0;i<cnt;i++){ 72 if(setIndex){treeNodes[i]->nodeData.nodeIndex=i;} 73 nodeCount[treeNodes[i]->depth()+1]++; 74 } 75 for(i=1;i<=maxDepth;i++){nodeCount[i]+=nodeCount[i-1];} 76} 77 78 79////////////////// 80// TreeNodeData // 81////////////////// 82int TreeNodeData::UseIndex=1; 83TreeNodeData::TreeNodeData(void){ 84 if(UseIndex){ 85 nodeIndex=-1; 86 centerWeightContribution=0; 87 } 88 else{mcIndex=0;} 89 value=0; 90} 91TreeNodeData::~TreeNodeData(void){;} 92 93 94//////////// 95// Octree // 96//////////// 97template<int Degree> 98double Octree<Degree>::maxMemoryUsage=0; 99 100//template<int Degree> 101//double Octree<Degree>::MemoryUsage(void){ 102// double mem=MemoryInfo::Usage()/(1<<20); 103// if(mem>maxMemoryUsage){maxMemoryUsage=mem;} 104// return mem; 105//} 106 107template<int Degree> 108Octree<Degree>::Octree(void){ 109 radius=0; 110 width=0; 111 postNormalSmooth=0; 112} 113 114template<int Degree> 115void Octree<Degree>::setNodeIndices(TreeOctNode& node,int& idx){ 116 node.nodeData.nodeIndex=idx; 117 idx++; 118 if(node.children){for(int i=0;i<Cube::CORNERS;i++){setNodeIndices(node.children[i],idx);}} 119} 120template<int Degree> 121int Octree<Degree>::NonLinearSplatOrientedPoint(TreeOctNode* node,const Point3D<Real>& position,const Point3D<Real>& normal){ 122 double x,dxdy,dxdydz,dx[DIMENSION][3]; 123 int i,j,k; 124 TreeOctNode::Neighbors& neighbors=neighborKey.setNeighbors(node); 125 double width; 126 Point3D<Real> center; 127 Real w; 128 129 node->centerAndWidth(center,w); 130 width=w; 131 for(int i=0;i<3;i++){ 132 x=(center.coords[i]-position.coords[i]-width)/width; 133 dx[i][0]=1.125+1.500*x+0.500*x*x; 134 x=(center.coords[i]-position.coords[i])/width; 135 dx[i][1]=0.750 - x*x; 136 dx[i][2]=1.0-dx[i][1]-dx[i][0]; 137 } 138 for(i=0;i<3;i++){ 139 for(j=0;j<3;j++){ 140 dxdy=dx[0][i]*dx[1][j]; 141 for(k=0;k<3;k++){ 142 if(neighbors.neighbors[i][j][k]){ 143 dxdydz=dxdy*dx[2][k]; 144 int idx=neighbors.neighbors[i][j][k]->nodeData.nodeIndex; 145 if(idx<0){ 146 Point3D<Real> n; 147 n.coords[0]=n.coords[1]=n.coords[2]=0; 148 idx=neighbors.neighbors[i][j][k]->nodeData.nodeIndex=int(normals->size()); 149 normals->push_back(n); 150 } 151 (*normals)[idx].coords[0]+=Real(normal.coords[0]*dxdydz); 152 (*normals)[idx].coords[1]+=Real(normal.coords[1]*dxdydz); 153 (*normals)[idx].coords[2]+=Real(normal.coords[2]*dxdydz); 154 } 155 } 156 } 157 } 158 return 0; 159} 160template<int Degree> 161void Octree<Degree>::NonLinearSplatOrientedPoint(const Point3D<Real>& position,const Point3D<Real>& normal,const int& splatDepth,const Real& samplesPerNode, 162 const int& minDepth,const int& maxDepth){ 163 double dx; 164 Point3D<Real> n; 165 TreeOctNode* temp; 166 int i,cnt=0; 167 double width; 168 Point3D<Real> myCenter; 169 Real myWidth; 170 myCenter.coords[0]=myCenter.coords[1]=myCenter.coords[2]=Real(0.5); 171 myWidth=Real(1.0); 172 173 temp=&tree; 174 while(temp->depth()<splatDepth){ 175 if(!temp->children){ 176 printf("Octree<Degree>::NonLinearSplatOrientedPoint error\n"); 177 return; 178 } 179 int cIndex=TreeOctNode::CornerIndex(myCenter,position); 180 temp=&temp->children[cIndex]; 181 myWidth/=2; 182 if(cIndex&1){myCenter.coords[0]+=myWidth/2;} 183 else {myCenter.coords[0]-=myWidth/2;} 184 if(cIndex&2){myCenter.coords[1]+=myWidth/2;} 185 else {myCenter.coords[1]-=myWidth/2;} 186 if(cIndex&4){myCenter.coords[2]+=myWidth/2;} 187 else {myCenter.coords[2]-=myWidth/2;} 188 } 189 Real alpha,newDepth; 190 NonLinearGetSampleDepthAndWeight(temp,position,samplesPerNode,newDepth,alpha); 191 192 if(newDepth<minDepth){newDepth=Real(minDepth);} 193 if(newDepth>maxDepth){newDepth=Real(maxDepth);} 194 int topDepth=int(ceil(newDepth)); 195 196 dx=1.0-(topDepth-newDepth); 197 if(topDepth<=minDepth){ 198 topDepth=minDepth; 199 dx=1; 200 } 201 else if(topDepth>maxDepth){ 202 topDepth=maxDepth; 203 dx=1; 204 } 205 while(temp->depth()>topDepth){temp=temp->parent;} 206 while(temp->depth()<topDepth){ 207 if(!temp->children){temp->initChildren();} 208 int cIndex=TreeOctNode::CornerIndex(myCenter,position); 209 temp=&temp->children[cIndex]; 210 myWidth/=2; 211 if(cIndex&1){myCenter.coords[0]+=myWidth/2;} 212 else {myCenter.coords[0]-=myWidth/2;} 213 if(cIndex&2){myCenter.coords[1]+=myWidth/2;} 214 else {myCenter.coords[1]-=myWidth/2;} 215 if(cIndex&4){myCenter.coords[2]+=myWidth/2;} 216 else {myCenter.coords[2]-=myWidth/2;} 217 } 218 width=1.0/(1<<temp->depth()); 219 for(i=0;i<DIMENSION;i++){n.coords[i]=normal.coords[i]*alpha/Real(pow(width,3))*Real(dx);} 220 NonLinearSplatOrientedPoint(temp,position,n); 221 if(fabs(1.0-dx)>EPSILON){ 222 dx=Real(1.0-dx); 223 temp=temp->parent; 224 width=1.0/(1<<temp->depth()); 225 226 for(i=0;i<DIMENSION;i++){n.coords[i]=normal.coords[i]*alpha/Real(pow(width,3))*Real(dx);} 227 NonLinearSplatOrientedPoint(temp,position,n); 228 } 229} 230template<int Degree> 231void Octree<Degree>::NonLinearGetSampleDepthAndWeight(TreeOctNode* node,const Point3D<Real>& position,const Real& samplesPerNode,Real& depth,Real& weight){ 232 TreeOctNode* temp=node; 233 weight=Real(1.0)/NonLinearGetSampleWeight(temp,position); 234 if(weight>=samplesPerNode+1){depth=Real(temp->depth()+log(weight/(samplesPerNode+1))/log(double(1<<(DIMENSION-1))));} 235 else{ 236 Real oldAlpha,newAlpha; 237 oldAlpha=newAlpha=weight; 238 while(newAlpha<(samplesPerNode+1) && temp->parent){ 239 temp=temp->parent; 240 oldAlpha=newAlpha; 241 newAlpha=Real(1.0)/NonLinearGetSampleWeight(temp,position); 242 } 243 depth=Real(temp->depth()+log(newAlpha/(samplesPerNode+1))/log(newAlpha/oldAlpha)); 244 } 245 weight=Real(pow(double(1<<(DIMENSION-1)),-double(depth))); 246} 247 248template<int Degree> 249Real Octree<Degree>::NonLinearGetSampleWeight(TreeOctNode* node,const Point3D<Real>& position){ 250 Real weight=0; 251 double x,dxdy,dx[DIMENSION][3]; 252 int i,j,k; 253 TreeOctNode::Neighbors& neighbors=neighborKey.setNeighbors(node); 254 double width; 255 Point3D<Real> center; 256 Real w; 257 node->centerAndWidth(center,w); 258 width=w; 259 260 for(i=0;i<DIMENSION;i++){ 261 x=(center.coords[i]-position.coords[i]-width)/width; 262 dx[i][0]=1.125+1.500*x+0.500*x*x; 263 x=(center.coords[i]-position.coords[i])/width; 264 dx[i][1]=0.750 - x*x; 265 dx[i][2]=1.0-dx[i][1]-dx[i][0]; 266 } 267 268 for(i=0;i<3;i++){ 269 for(j=0;j<3;j++){ 270 dxdy=dx[0][i]*dx[1][j]; 271 for(k=0;k<3;k++){ 272 if(neighbors.neighbors[i][j][k]){ 273 weight+=Real(dxdy*dx[2][k]*neighbors.neighbors[i][j][k]->nodeData.centerWeightContribution); 274 } 275 } 276 } 277 } 278 return Real(1.0/weight); 279} 280 281template<int Degree> 282int Octree<Degree>::NonLinearUpdateWeightContribution(TreeOctNode* node,const Point3D<Real>& position,const Real& weight){ 283 int i,j,k; 284 TreeOctNode::Neighbors& neighbors=neighborKey.setNeighbors(node); 285 double x,dxdy,dx[DIMENSION][3]; 286 double width; 287 Point3D<Real> center; 288 Real w; 289 node->centerAndWidth(center,w); 290 width=w; 291 292 for(i=0;i<DIMENSION;i++){ 293 x=(center.coords[i]-position.coords[i]-width)/width; 294 dx[i][0]=1.125+1.500*x+0.500*x*x; 295 x=(center.coords[i]-position.coords[i])/width; 296 dx[i][1]=0.750 - x*x; 297 dx[i][2]=1.0-dx[i][1]-dx[i][0]; 298 } 299 300 for(i=0;i<3;i++){ 301 for(j=0;j<3;j++){ 302 dxdy=dx[0][i]*dx[1][j]*weight; 303 for(k=0;k<3;k++){ 304 if(neighbors.neighbors[i][j][k]){neighbors.neighbors[i][j][k]->nodeData.centerWeightContribution+=Real(dxdy*dx[2][k]);} 305 } 306 } 307 } 308 return 0; 309} 310 311template<int Degree> 312int Octree<Degree>::setTree(std::vector<Point3D<Real> > &Pts, std::vector<Point3D<Real> > &Nor, 313 const int& maxDepth, 314 const int& kernelDepth,const Real& samplesPerNode,const Real& scaleFactor,Point3D<Real>& center,Real& scale, 315 const int& resetSamples,const int& useConfidence){ 316 317 Point3D<Real> min,max,position,normal,myCenter; 318 Real myWidth; 319 TreeOctNode* temp; 320 int splatDepth=0; 321 float c[2*DIMENSION]; 322 323 TreeNodeData::UseIndex=1; 324 neighborKey.set(maxDepth); 325 splatDepth=kernelDepth; 326 if(splatDepth<0){splatDepth=0;} 327 328 //DumpOutput("Setting bounding box\n"); 329 // Read through once to get the center and scale 330 int i,cnt=0; 331 for(int pi=0;pi<Pts.size();++pi) 332 { 333 //if(binary){if(fread(c,sizeof(float),2*DIMENSION,fp)!=6){break;}} 334 //else{if(fscanf(fp," %f %f %f %f %f %f ",&c[0],&c[1],&c[2],&c[3],&c[4],&c[5])!=2*DIMENSION){break;}} 335 for(i=0;i<DIMENSION;i++){ 336 c[i]=Pts[pi].coords[i]; 337 if(!cnt || c[i]<min.coords[i]){min.coords[i]=c[i];} 338 if(!cnt || c[i]>max.coords[i]){max.coords[i]=c[i];} 339 } 340 cnt++; 341 } 342 for(i=0;i<DIMENSION;i++){ 343 if(!i || scale<max.coords[i]-min.coords[i]){scale=Real(max.coords[i]-min.coords[i]);} 344 center.coords[i]=Real(max.coords[i]+min.coords[i])/2; 345 } 346 printf("Samples: %d\n",cnt); 347 scale*=scaleFactor; 348 for(i=0;i<DIMENSION;i++){center.coords[i]-=scale/2;} 349 if(splatDepth>0){ 350 printf("Setting sample weights\n"); 351 cnt=0; 352 //fseek(fp,SEEK_SET,0); 353 for(int pi=0;pi<Pts.size();++pi) 354 { 355 //if(binary){if(fread(c,sizeof(float),2*DIMENSION,fp)!=2*DIMENSION){break;}} 356 //else{if(fscanf(fp," %f %f %f %f %f %f ",&c[0],&c[1],&c[2],&c[3],&c[4],&c[5])!=2*DIMENSION){break;}} 357 for(i=0;i<DIMENSION;i++){ 358 c[i]=Pts[pi].coords[i]; 359 c[i+3]=Nor[pi].coords[i]; 360 position.coords[i]=(c[i]-center.coords[i])/scale; 361 normal.coords[i]=c[DIMENSION+i]; 362 } 363 myCenter.coords[0]=myCenter.coords[1]=myCenter.coords[2]=Real(0.5); 364 myWidth=Real(1.0); 365 for(i=0;i<DIMENSION;i++){if(position.coords[i]<myCenter.coords[i]-myWidth/2 || position.coords[i]>myCenter.coords[i]+myWidth/2){break;}} 366 if(i!=DIMENSION){continue;} 367 temp=&tree; 368 int d=0; 369 Real weight=Real(1.0); 370 if(useConfidence){weight=Real(Length(normal));} 371 while(d<splatDepth){ 372 NonLinearUpdateWeightContribution(temp,position,weight); 373 if(!temp->children){temp->initChildren();} 374 int cIndex=TreeOctNode::CornerIndex(myCenter,position); 375 temp=&temp->children[cIndex]; 376 myWidth/=2; 377 if(cIndex&1){myCenter.coords[0]+=myWidth/2;} 378 else {myCenter.coords[0]-=myWidth/2;} 379 if(cIndex&2){myCenter.coords[1]+=myWidth/2;} 380 else {myCenter.coords[1]-=myWidth/2;} 381 if(cIndex&4){myCenter.coords[2]+=myWidth/2;} 382 else {myCenter.coords[2]-=myWidth/2;} 383 d++; 384 } 385 NonLinearUpdateWeightContribution(temp,position,weight); 386 cnt++; 387 } 388 } 389 390 printf("Adding Points and Normals\n"); 391 normals=new std::vector<Point3D<Real> >(); 392 cnt=0; 393 //fseek(fp,SEEK_SET,0); 394 for(int pi=0;pi<Pts.size();++pi) 395 { 396 //if(binary){if(fread(c,sizeof(float),2*DIMENSION,fp)!=2*DIMENSION){break;}} 397 //else{if(fscanf(fp," %f %f %f %f %f %f ",&c[0],&c[1],&c[2],&c[3],&c[4],&c[5])!=2*DIMENSION){break;}} 398 for(i=0;i<DIMENSION;i++){ 399 c[i]=Pts[pi].coords[i]; 400 c[i+3]=Nor[pi].coords[i]; 401 position.coords[i]=(c[i]-center.coords[i])/scale; 402 normal.coords[i]=c[DIMENSION+i]; 403 } 404 myCenter.coords[0]=myCenter.coords[1]=myCenter.coords[2]=Real(0.5); 405 myWidth=Real(1.0); 406 for(i=0;i<DIMENSION;i++){if(position.coords[i]<myCenter.coords[i]-myWidth/2 || position.coords[i]>myCenter.coords[i]+myWidth/2){break;}} 407 if(i!=DIMENSION){continue;} 408 Real l=Real(Length(normal)); 409 if(l<EPSILON){continue;} 410 if(!useConfidence){ 411 normal.coords[0]/=l; 412 normal.coords[1]/=l; 413 normal.coords[2]/=l; 414 } 415 l=Real(2<<maxDepth); 416 normal.coords[0]*=l; 417 normal.coords[1]*=l; 418 normal.coords[2]*=l; 419 420 if(resetSamples && samplesPerNode>0 && splatDepth){ 421 NonLinearSplatOrientedPoint(position,normal,splatDepth,samplesPerNode,1,maxDepth); 422 } 423 else{ 424 Real alpha=1; 425 temp=&tree; 426 if(splatDepth){ 427 int d=0; 428 while(d<splatDepth){ 429 int cIndex=TreeOctNode::CornerIndex(myCenter,position); 430 temp=&temp->children[cIndex]; 431 myWidth/=2; 432 if(cIndex&1){myCenter.coords[0]+=myWidth/2;} 433 else {myCenter.coords[0]-=myWidth/2;} 434 if(cIndex&2){myCenter.coords[1]+=myWidth/2;} 435 else {myCenter.coords[1]-=myWidth/2;} 436 if(cIndex&4){myCenter.coords[2]+=myWidth/2;} 437 else {myCenter.coords[2]-=myWidth/2;} 438 d++; 439 } 440 alpha=NonLinearGetSampleWeight(temp,position); 441 } 442 for(i=0;i<DIMENSION;i++){normal.coords[i]*=alpha;} 443 int d=0; 444 while(d<maxDepth){ 445 if(!temp->children){temp->initChildren();} 446 int cIndex=TreeOctNode::CornerIndex(myCenter,position); 447 temp=&temp->children[cIndex]; 448 myWidth/=2; 449 if(cIndex&1){myCenter.coords[0]+=myWidth/2;} 450 else {myCenter.coords[0]-=myWidth/2;} 451 if(cIndex&2){myCenter.coords[1]+=myWidth/2;} 452 else {myCenter.coords[1]-=myWidth/2;} 453 if(cIndex&4){myCenter.coords[2]+=myWidth/2;} 454 else {myCenter.coords[2]-=myWidth/2;} 455 d++; 456 } 457 NonLinearSplatOrientedPoint(temp,position,normal); 458 } 459 } 460// DumpOutput("Memory Usage: %.3f MB\n",float(MemoryUsage())); 461// fclose(fp); 462 return cnt; 463} 464 465template<int Degree> 466void Octree<Degree>::setFunctionData(const PPolynomial<Degree>& ReconstructionFunction, const int& maxDepth,const int& normalize,const Real& normalSmooth){ 467 468 radius=Real(fabs(ReconstructionFunction.polys[0].start)); 469 width=int(double(radius+0.5-EPSILON)*2); 470 if(normalSmooth>0){postNormalSmooth=normalSmooth;} 471 fData.set(maxDepth,ReconstructionFunction,normalize,1); 472} 473 474template<int Degree> 475void Octree<Degree>::finalize1(const int& refineNeighbors) 476{ 477 TreeOctNode* temp; 478 479 if(refineNeighbors>=0){ 480 RefineFunction rf; 481 temp=tree.nextNode(); 482 while(temp){ 483 if(temp->nodeData.nodeIndex>=0 && Length((*normals)[temp->nodeData.nodeIndex])>EPSILON){ 484 rf.depth=temp->depth()-refineNeighbors; 485 TreeOctNode::ProcessMaxDepthNodeAdjacentNodes(fData.depth,temp,2*width,&tree,1,temp->depth()-refineNeighbors,&rf); 486 } 487 temp=tree.nextNode(temp); 488 } 489 } 490 else if(refineNeighbors==-1234){ 491 temp=tree.nextLeaf(); 492 while(temp){ 493 if(!temp->children && temp->depth()<fData.depth){temp->initChildren();} 494 temp=tree.nextLeaf(temp); 495 } 496 } 497} 498template<int Degree> 499void Octree<Degree>::finalize2(const int& refineNeighbors) 500{ 501 TreeOctNode* temp; 502 503 if(refineNeighbors>=0){ 504 RefineFunction rf; 505 temp=tree.nextNode(); 506 while(temp){ 507 if(fabs(temp->nodeData.value)>EPSILON){ 508 rf.depth=temp->depth()-refineNeighbors; 509 TreeOctNode::ProcessMaxDepthNodeAdjacentNodes(fData.depth,temp,2*width,&tree,1,temp->depth()-refineNeighbors,&rf); 510 } 511 temp=tree.nextNode(temp); 512 } 513 } 514} 515template <int Degree> 516Real Octree<Degree>::GetDivergence(const int idx[DIMENSION],const Point3D<Real>& normal) const 517{ 518 double dot=fData.dotTable[idx[0]]*fData.dotTable[idx[1]]*fData.dotTable[idx[2]]; 519 return Real(dot*(fData.dDotTable[idx[0]]*normal.coords[0]+fData.dDotTable[idx[1]]*normal.coords[1]+fData.dDotTable[idx[2]]*normal.coords[2])); 520} 521template<int Degree> 522Real Octree<Degree>::GetLaplacian(const int idx[DIMENSION]) const 523{ 524 return Real(fData.dotTable[idx[0]]*fData.dotTable[idx[1]]*fData.dotTable[idx[2]]*(fData.d2DotTable[idx[0]]+fData.d2DotTable[idx[1]]+fData.d2DotTable[idx[2]])); 525} 526template<int Degree> 527Real Octree<Degree>::GetDotProduct(const int idx[DIMENSION]) const 528{ 529 return Real(fData.dotTable[idx[0]]*fData.dotTable[idx[1]]*fData.dotTable[idx[2]]); 530} 531 532template<int Degree> 533int Octree<Degree>::GetFixedDepthLaplacian(SparseSymmetricMatrix<float>& matrix,const int& depth,const SortedTreeNodes& sNodes) 534{ 535 LaplacianMatrixFunction mf; 536 mf.ot=this; 537 mf.offset=sNodes.nodeCount[depth]; 538 matrix.Resize(sNodes.nodeCount[depth+1]-sNodes.nodeCount[depth]); 539 mf.rowElements=(MatrixEntry<float>*)malloc(sizeof(MatrixEntry<float>)*matrix.rows); 540 for(int i=sNodes.nodeCount[depth];i<sNodes.nodeCount[depth+1];i++){ 541 mf.elementCount=0; 542 mf.d2=int(sNodes.treeNodes[i]->d); 543 mf.x2=int(sNodes.treeNodes[i]->off[0]); 544 mf.y2=int(sNodes.treeNodes[i]->off[1]); 545 mf.z2=int(sNodes.treeNodes[i]->off[2]); 546 mf.index[0]=mf.x2; 547 mf.index[1]=mf.y2; 548 mf.index[2]=mf.z2; 549 TreeOctNode::ProcessTerminatingNodeAdjacentNodes(fData.depth,sNodes.treeNodes[i],2*width-1,&tree,1,&mf); 550 matrix.SetRowSize(i-sNodes.nodeCount[depth],mf.elementCount); 551 memcpy(matrix.m_ppElements[i-sNodes.nodeCount[depth]],mf.rowElements,sizeof(MatrixEntry<float>)*mf.elementCount); 552 } 553 free(mf.rowElements); 554 return 1; 555} 556template<int Degree> 557int Octree<Degree>::GetRestrictedFixedDepthLaplacian(SparseSymmetricMatrix<float>& matrix,const int& depth,const int* entries,const int& entryCount, 558 const TreeOctNode* rNode,const Real& radius, 559 const SortedTreeNodes& sNodes){ 560 int i; 561 RestrictedLaplacianMatrixFunction mf; 562 Real myRadius=int(2*radius-ROUND_EPS)+ROUND_EPS; 563 mf.ot=this; 564 mf.radius=radius; 565 rNode->depthAndOffset(mf.depth,mf.offset); 566 matrix.Resize(entryCount); 567 mf.rowElements=(MatrixEntry<float>*)malloc(sizeof(MatrixEntry<float>)*matrix.rows); 568 for(i=0;i<entryCount;i++){sNodes.treeNodes[entries[i]]->nodeData.nodeIndex=i;} 569 for(i=0;i<entryCount;i++){ 570 mf.elementCount=0; 571 mf.index[0]=int(sNodes.treeNodes[entries[i]]->off[0]); 572 mf.index[1]=int(sNodes.treeNodes[entries[i]]->off[1]); 573 mf.index[2]=int(sNodes.treeNodes[entries[i]]->off[2]); 574 TreeOctNode::ProcessTerminatingNodeAdjacentNodes(fData.depth,sNodes.treeNodes[entries[i]],2*width-1,&tree,1,&mf); 575 matrix.SetRowSize(i,mf.elementCount); 576 memcpy(matrix.m_ppElements[i],mf.rowElements,sizeof(MatrixEntry<float>)*mf.elementCount); 577 } 578 for(i=0;i<entryCount;i++){sNodes.treeNodes[entries[i]]->nodeData.nodeIndex=entries[i];} 579 free(mf.rowElements); 580 return 1; 581} 582 583 584template<int Degree> 585int Octree<Degree>::LaplacianMatrixIteration(const int& subdivideDepth){ 586 int i,iter=0; 587 SortedTreeNodes sNodes; 588 double t; 589 fData.setDotTables(fData.D2_DOT_FLAG); 590 sNodes.set(tree,1); 591 592 SparseMatrix<float>::SetAllocator(MEMORY_ALLOCATOR_BLOCK_SIZE); 593 594 sNodes.treeNodes[0]->nodeData.value=0; 595 for(i=1;i<sNodes.maxDepth;i++){ 596 printf("Depth: %d/%d\n",i,sNodes.maxDepth-1); 597// t=Time(); 598 if(subdivideDepth>0){iter+=SolveFixedDepthMatrix(i,subdivideDepth,sNodes);} 599 else{iter+=SolveFixedDepthMatrix(i,sNodes);} 600 } 601 SparseMatrix<float>::Allocator.reset(); 602 fData.clearDotTables(fData.DOT_FLAG | fData.D_DOT_FLAG | fData.D2_DOT_FLAG); 603 return iter; 604} 605 606template<int Degree> 607int Octree<Degree>::SolveFixedDepthMatrix(const int& depth,const SortedTreeNodes& sNodes){ 608 int i,iter=0; 609 Vector<double> V,Solution; 610 SparseSymmetricMatrix<Real> matrix; 611 Real myRadius; 612 double gTime,sTime,uTime; 613 Real dx,dy,dz; 614 int x1,x2,y1,y2,z1,z2; 615 Vector<Real> Diagonal; 616 617// gTime=Time(); 618 V.Resize(sNodes.nodeCount[depth+1]-sNodes.nodeCount[depth]); 619 for(i=sNodes.nodeCount[depth];i<sNodes.nodeCount[depth+1];i++){V[i-sNodes.nodeCount[depth]]=sNodes.treeNodes[i]->nodeData.value;} 620 SparseSymmetricMatrix<float>::Allocator.rollBack(); 621 GetFixedDepthLaplacian(matrix,depth,sNodes); 622// gTime=Time()-gTime; 623 printf("\tMatrix entries: %d / %d^2 = %.4f%%\n",matrix.Entries(),matrix.rows,100.0*(matrix.Entries()/double(matrix.rows))/matrix.rows); 624// DumpOutput("\tMemory Usage: %.3f MB\n",float(MemoryUsage())); 625// sTime=Time(); 626 iter+=SparseSymmetricMatrix<Real>::Solve(matrix,V,int(pow(matrix.rows,ITERATION_POWER)),Solution,double(EPSILON),1); 627// sTime=Time()-sTime; 628// uTime=Time(); 629 for(i=sNodes.nodeCount[depth];i<sNodes.nodeCount[depth+1];i++){sNodes.treeNodes[i]->nodeData.value=Real(Solution[i-sNodes.nodeCount[depth]]);} 630 631 myRadius=Real(radius+ROUND_EPS-0.5); 632 myRadius /=(1<<depth); 633 634 if(depth<sNodes.maxDepth-1){ 635 LaplacianProjectionFunction pf; 636 TreeOctNode *node1,*node2; 637 pf.ot=this; 638 int idx1,idx2,off=sNodes.nodeCount[depth]; 639 // First pass: idx2 is the solution coefficient propogated 640 for(i=0;i<matrix.rows;i++){ 641 idx1=i; 642 node1=sNodes.treeNodes[idx1+off]; 643 if(!node1->children){continue;} 644 x1=int(node1->off[0]); 645 y1=int(node1->off[1]); 646 z1=int(node1->off[2]); 647 for(int j=0;j<matrix.rowSizes[i];j++){ 648 idx2=matrix.m_ppElements[i][j].N; 649 node2=sNodes.treeNodes[idx2+off]; 650 x2=int(node2->off[0]); 651 y2=int(node2->off[1]); 652 z2=int(node2->off[2]); 653 pf.value=Solution[idx2]; 654 pf.index[0]=x2; 655 pf.index[1]=y2; 656 pf.index[2]=z2; 657 dx=Real(x2-x1)/(1<<depth); 658 dy=Real(y2-y1)/(1<<depth); 659 dz=Real(z2-z1)/(1<<depth); 660 if(fabs(dx)<myRadius && fabs(dy)<myRadius && fabs(dz)<myRadius){node1->processNodeNodes(node2,&pf,0);} 661 else{TreeOctNode::ProcessNodeAdjacentNodes(fData.depth,node2,width,node1,width,&pf,0);} 662 } 663 } 664 // Second pass: idx1 is the solution coefficient propogated 665 for(i=0;i<matrix.rows;i++){ 666 idx1=i; 667 node1=sNodes.treeNodes[idx1+off]; 668 x1=int(node1->off[0]); 669 y1=int(node1->off[1]); 670 z1=int(node1->off[2]); 671 pf.value=Solution[idx1]; 672 pf.index[0]=x1; 673 pf.index[1]=y1; 674 pf.index[2]=z1; 675 for(int j=0;j<matrix.rowSizes[i];j++){ 676 idx2=matrix.m_ppElements[i][j].N; 677 node2=sNodes.treeNodes[idx2+off]; 678 if(idx1!=idx2 && node2->children){ 679 x2=int(node2->off[0]); 680 y2=int(node2->off[1]); 681 z2=int(node2->off[2]); 682 dx=Real(x1-x2)/(1<<depth); 683 dy=Real(y1-y2)/(1<<depth); 684 dz=Real(z1-z2)/(1<<depth); 685 if(fabs(dx)<myRadius && fabs(dy)<myRadius && fabs(dz)<myRadius){node2->processNodeNodes(node1,&pf,0);} 686 else{TreeOctNode::ProcessNodeAdjacentNodes(fData.depth,node1,width,node2,width,&pf,0);} 687 } 688 } 689 } 690 } 691// uTime=Time()-uTime; 692 printf("\tGot / Solved / Updated in: %6.3f / %6.3f / %6.3f\n",gTime,sTime,uTime); 693 return iter; 694} 695template<int Degree> 696int Octree<Degree>::SolveFixedDepthMatrix(const int& depth,const int& startingDepth,const SortedTreeNodes& sNodes){ 697 int i,j,d,iter=0; 698 SparseSymmetricMatrix<Real> matrix; 699 AdjacencySetFunction asf; 700 AdjacencyCountFunction acf; 701 Vector<Real> Values; 702 Vector<double> SubValues,SubSolution; 703 double gTime,sTime,uTime; 704 Real myRadius,myRadius2; 705 Real dx,dy,dz; 706 Vector<Real> Diagonal; 707 708 if(startingDepth>=depth){return SolveFixedDepthMatrix(depth,sNodes);} 709 710 Values.Resize(sNodes.nodeCount[depth+1]-sNodes.nodeCount[depth]); 711 712 for(i=sNodes.nodeCount[depth];i<sNodes.nodeCount[depth+1];i++){ 713 Values[i-sNodes.nodeCount[depth]]=sNodes.treeNodes[i]->nodeData.value; 714 sNodes.treeNodes[i]->nodeData.value=0; 715 } 716 717 myRadius=2*radius-Real(0.5); 718 myRadius=int(myRadius-ROUND_EPS)+ROUND_EPS; 719 myRadius2=Real(radius+ROUND_EPS-0.5); 720 d=depth-startingDepth; 721 for(i=sNodes.nodeCount[d];i<sNodes.nodeCount[d+1];i++){ 722// gTime=Time(); 723 TreeOctNode* temp; 724 // Get all of the entries associated to the subspace 725 acf.adjacencyCount=0; 726 temp=sNodes.treeNodes[i]->nextNode(); 727 while(temp){ 728 if(temp->depth()==depth){ 729 acf.Function(temp,temp); 730 temp=sNodes.treeNodes[i]->nextBranch(temp); 731 } 732 else{temp=sNodes.treeNodes[i]->nextNode(temp);} 733 } 734 for(j=sNodes.nodeCount[d];j<sNodes.nodeCount[d+1];j++){ 735 if(i==j){continue;} 736 TreeOctNode::ProcessFixedDepthNodeAdjacentNodes(fData.depth,sNodes.treeNodes[i],1,sNodes.treeNodes[j],2*width-1,depth,&acf); 737 } 738 if(!acf.adjacencyCount){continue;} 739 asf.adjacencies=new int[acf.adjacencyCount]; 740 asf.adjacencyCount=0; 741 temp=sNodes.treeNodes[i]->nextNode(); 742 while(temp){ 743 if(temp->depth()==depth){ 744 asf.Function(temp,temp); 745 temp=sNodes.treeNodes[i]->nextBranch(temp); 746 } 747 else{temp=sNodes.treeNodes[i]->nextNode(temp);} 748 } 749 for(j=sNodes.nodeCount[d];j<sNodes.nodeCount[d+1];j++){ 750 if(i==j){continue;} 751 TreeOctNode::ProcessFixedDepthNodeAdjacentNodes(fData.depth,sNodes.treeNodes[i],1,sNodes.treeNodes[j],2*width-1,depth,&asf); 752 } 753 754 printf("\tNodes[%d/%d]: %d\n",i-sNodes.nodeCount[d]+1,sNodes.nodeCount[d+1]-sNodes.nodeCount[d],asf.adjacencyCount); 755 // Get the associated vector 756 SubValues.Resize(asf.adjacencyCount); 757 for(j=0;j<asf.adjacencyCount;j++){SubValues[j]=Values[asf.adjacencies[j]-sNodes.nodeCount[depth]];} 758 SubSolution.Resize(asf.adjacencyCount); 759 for(j=0;j<asf.adjacencyCount;j++){SubSolution[j]=sNodes.treeNodes[asf.adjacencies[j]]->nodeData.value;} 760 // Get the associated matrix 761 SparseSymmetricMatrix<float>::Allocator.rollBack(); 762 GetRestrictedFixedDepthLaplacian(matrix,depth,asf.adjacencies,asf.adjacencyCount,sNodes.treeNodes[i],myRadius,sNodes); 763// gTime=Time()-gTime; 764 printf("\t\tMatrix entries: %d / %d^2 = %.4f%%\n",matrix.Entries(),matrix.rows,100.0*(matrix.Entries()/double(matrix.rows))/matrix.rows); 765// DumpOutput("\t\tMemory Usage: %.3f MB\n",float(MemoryUsage())); 766 767 // Solve the matrix 768// sTime=Time(); 769 iter+=SparseSymmetricMatrix<Real>::Solve(matrix,SubValues,int(pow(matrix.rows,ITERATION_POWER)),SubSolution,double(EPSILON),0); 770// sTime=Time()-sTime; 771 772// uTime=Time(); 773 LaplacianProjectionFunction lpf; 774 lpf.ot=this; 775 776 // Update the solution for all nodes in the sub-tree 777 for(j=0;j<asf.adjacencyCount;j++){ 778 temp=sNodes.treeNodes[asf.adjacencies[j]]; 779 while(temp->depth()>sNodes.treeNodes[i]->depth()){temp=temp->parent;} 780 if(temp->nodeData.nodeIndex>=sNodes.treeNodes[i]->nodeData.nodeIndex){sNodes.treeNodes[asf.adjacencies[j]]->nodeData.value=Real(SubSolution[j]);} 781 } 782// double t=Time(); 783 // Update the values in the next depth 784 int x1,x2,y1,y2,z1,z2; 785 if(depth<sNodes.maxDepth-1){ 786 int idx1,idx2; 787 TreeOctNode *node1,*node2; 788 // First pass: idx2 is the solution coefficient propogated 789 for(j=0;j<matrix.rows;j++){ 790 idx1=asf.adjacencies[j]; 791 node1=sNodes.treeNodes[idx1]; 792 if(!node1->children){continue;} 793 x1=int(node1->off[0]); 794 y1=int(node1->off[1]); 795 z1=int(node1->off[2]); 796 797 for(int k=0;k<matrix.rowSizes[j];k++){ 798 idx2=asf.adjacencies[matrix.m_ppElements[j][k].N]; 799 node2=sNodes.treeNodes[idx2]; 800 temp=node2; 801 while(temp->depth()>d){temp=temp->parent;} 802 if(temp!=sNodes.treeNodes[i]){continue;} 803 lpf.value=Real(SubSolution[matrix.m_ppElements[j][k].N]); 804 x2=int(node2->off[0]); 805 y2=int(node2->off[1]); 806 z2=int(node2->off[2]); 807 lpf.index[0]=x2; 808 lpf.index[1]=y2; 809 lpf.index[2]=z2; 810 dx=Real(x2-x1)/(1<<depth); 811 dy=Real(y2-y1)/(1<<depth); 812 dz=Real(z2-z1)/(1<<depth); 813 if(fabs(dx)<myRadius2 && fabs(dy)<myRadius2 && fabs(dz)<myRadius2){node1->processNodeNodes(node2,&lpf,0);} 814 else{TreeOctNode::ProcessNodeAdjacentNodes(fData.depth,node2,width,node1,width,&lpf,0);} 815 } 816 } 817 // Second pass: idx1 is the solution coefficient propogated 818 for(j=0;j<matrix.rows;j++){ 819 idx1=asf.adjacencies[j]; 820 node1=sNodes.treeNodes[idx1]; 821 temp=node1; 822 while(temp->depth()>d){temp=temp->parent;} 823 if(temp!=sNodes.treeNodes[i]){continue;} 824 x1=int(node1->off[0]); 825 y1=int(node1->off[1]); 826 z1=int(node1->off[2]); 827 828 lpf.value=Real(SubSolution[j]); 829 lpf.index[0]=x1; 830 lpf.index[1]=y1; 831 lpf.index[2]=z1; 832 for(int k=0;k<matrix.rowSizes[j];k++){ 833 idx2=asf.adjacencies[matrix.m_ppElements[j][k].N]; 834 node2=sNodes.treeNodes[idx2]; 835 if(!node2->children){continue;} 836 837 if(idx1!=idx2){ 838 x2=int(node2->off[0]); 839 y2=int(node2->off[1]); 840 z2=int(node2->off[2]); 841 dx=Real(x1-x2)/(1<<depth); 842 dy=Real(y1-y2)/(1<<depth); 843 dz=Real(z1-z2)/(1<<depth); 844 if(fabs(dx)<myRadius2 && fabs(dy)<myRadius2 && fabs(dz)<myRadius2){node2->processNodeNodes(node1,&lpf,0);} 845 else{TreeOctNode::ProcessNodeAdjacentNodes(fData.depth,node1,width,node2,width,&lpf,0);} 846 } 847 } 848 } 849 } 850// uTime=Time()-uTime; 851 printf("\t\tGot / Solved / Updated in: %6.3f / %6.3f / %6.3f\n",gTime,sTime,uTime); 852 delete[] asf.adjacencies; 853 } 854 return iter; 855} 856template<int Degree> 857int Octree<Degree>::HasNormals(TreeOctNode* node,const Real& epsilon){ 858 int hasNormals=0; 859 if(node->nodeData.nodeIndex>=0 && Length((*normals)[node->nodeData.nodeIndex])>epsilon){hasNormals=1;} 860 if(node->children){for(int i=0;i<Cube::CORNERS && !hasNormals;i++){hasNormals|=HasNormals(&node->children[i],epsilon);}} 861 862 return hasNormals; 863} 864template<int Degree> 865void Octree<Degree>::ClipTree(void){ 866 TreeOctNode* temp; 867 temp=tree.nextNode(); 868 while(temp){ 869 if(temp->children){ 870 int hasNormals=0; 871 for(int i=0;i<Cube::CORNERS && !hasNormals;i++){hasNormals=HasNormals(&temp->children[i],EPSILON);} 872 if(!hasNormals){temp->children=NULL;} 873 } 874 temp=tree.nextNode(temp); 875 } 876} 877template<int Degree> 878void Octree<Degree>::SetLaplacianWeights(void){ 879 TreeOctNode* temp; 880 881 fData.setDotTables(fData.DOT_FLAG | fData.D_DOT_FLAG); 882 DivergenceFunction df; 883 df.ot=this; 884 temp=tree.nextNode(); 885 while(temp){ 886 if(temp->nodeData.nodeIndex<0 || Length((*normals)[temp->nodeData.nodeIndex])<=EPSILON){ 887 temp=tree.nextNode(temp); 888 continue; 889 } 890 int d=temp->depth(); 891 df.normal=(*normals)[temp->nodeData.nodeIndex]; 892 df.index[0]=int(temp->off[0]); 893 df.index[1]=int(temp->off[1]); 894 df.index[2]=int(temp->off[2]); 895 TreeOctNode::ProcessNodeAdjacentNodes(fData.depth,temp,width,&tree,width,&df); 896 temp=tree.nextNode(temp); 897 } 898 fData.clearDotTables(fData.D_DOT_FLAG); 899 temp=tree.nextNode(); 900 while(temp){ 901 if(temp->nodeData.nodeIndex<0){temp->nodeData.centerWeightContribution=0;} 902 else{temp->nodeData.centerWeightContribution=Real(Length((*normals)[temp->nodeData.nodeIndex]));} 903 temp=tree.nextNode(temp); 904 } 905// MemoryUsage(); 906 907 delete normals; 908 normals=NULL; 909} 910template<int Degree> 911void Octree<Degree>::DivergenceFunction::Function(TreeOctNode* node1,const TreeOctNode* node2){ 912 Point3D<Real> n=normal; 913 if(FunctionData<Degree,Real>::SymmetricIndex(index[0],int(node1->off[0]),scratch[0])){n.coords[0]=-n.coords[0];} 914 if(FunctionData<Degree,Real>::SymmetricIndex(index[1],int(node1->off[1]),scratch[1])){n.coords[1]=-n.coords[1];} 915 if(FunctionData<Degree,Real>::SymmetricIndex(index[2],int(node1->off[2]),scratch[2])){n.coords[2]=-n.coords[2];} 916 double dot=ot->fData.dotTable[scratch[0]]*ot->fData.dotTable[scratch[1]]*ot->fData.dotTable[scratch[2]]; 917 node1->nodeData.value+=Real(dot*(ot->fData.dDotTable[scratch[0]]*n.coords[0]+ot->fData.dDotTable[scratch[1]]*n.coords[1]+ot->fData.dDotTable[scratch[2]]*n.coords[2])); 918} 919template<int Degree> 920void Octree<Degree>::LaplacianProjectionFunction::Function(TreeOctNode* node1,const TreeOctNode* node2){ 921 scratch[0]=FunctionData<Degree,Real>::SymmetricIndex(index[0],int(node1->off[0])); 922 scratch[1]=FunctionData<Degree,Real>::SymmetricIndex(index[1],int(node1->off[1])); 923 scratch[2]=FunctionData<Degree,Real>::SymmetricIndex(index[2],int(node1->off[2])); 924 node1->nodeData.value-=Real(ot->GetLaplacian(scratch)*value); 925} 926template<int Degree> 927void Octree<Degree>::AdjacencyCountFunction::Function(const TreeOctNode* node1,const TreeOctNode* node2){adjacencyCount++;} 928template<int Degree> 929void Octree<Degree>::AdjacencySetFunction::Function(const TreeOctNode* node1,const TreeOctNode* node2){adjacencies[adjacencyCount++]=node1->nodeData.nodeIndex;} 930template<int Degree> 931void Octree<Degree>::RefineFunction::Function(TreeOctNode* node1,const TreeOctNode* node2){ 932 if(!node1->children && node1->depth()<depth){node1->initChildren();} 933} 934template<int Degree> 935void Octree<Degree>::FaceEdgesFunction::Function(const TreeOctNode* node1,const TreeOctNode* node2){ 936 if(!node1->children && MarchingCubes::HasRoots(node1->nodeData.mcIndex)){ 937 RootInfo ri1,ri2; 938 hash_map<long long,std::pair<RootInfo,int> >::iterator iter; 939 int isoTri[DIMENSION*MarchingCubes::MAX_TRIANGLES]; 940 int count=MarchingCubes::AddTriangleIndices(node1->nodeData.mcIndex,isoTri); 941 942 for(int j=0;j<count;j++){ 943 for(int k=0;k<3;k++){ 944 if(fIndex==Cube::FaceAdjacentToEdges(isoTri[j*3+k],isoTri[j*3+((k+1)%3)])){ 945 if(GetRootIndex(node1,isoTri[j*3+k],maxDepth,ri1) && GetRootIndex(node1,isoTri[j*3+((k+1)%3)],maxDepth,ri2)){ 946 edges->push_back(std::pair<long long,long long>(ri2.key,ri1.key)); 947 iter=vertexCount->find(ri1.key); 948 if(iter==vertexCount->end()){ 949 (*vertexCount)[ri1.key].first=ri1; 950 (*vertexCount)[ri1.key].second=0; 951 } 952 iter=vertexCount->find(ri2.key); 953 if(iter==vertexCount->end()){ 954 (*vertexCount)[ri2.key].first=ri2; 955 (*vertexCount)[ri2.key].second=0; 956 } 957 (*vertexCount)[ri1.key].second--; 958 (*vertexCount)[ri2.key].second++; 959 } 960 else{fprintf(stderr,"Bad Edge 1: %d %d\n",ri1.key,ri2.key);} 961 } 962 } 963 } 964 } 965} 966template<int Degree> 967void Octree<Degree>::PointIndexValueFunction::Function(const TreeOctNode* node){ 968 int idx[DIMENSION]; 969 idx[0]=index[0]+int(node->off[0]); 970 idx[1]=index[1]+int(node->off[1]); 971 idx[2]=index[2]+int(node->off[2]); 972 value+=node->nodeData.value* Real( valueTables[idx[0]]* valueTables[idx[1]]* valueTables[idx[2]]); 973} 974template<int Degree> 975void Octree<Degree>::PointIndexValueAndNormalFunction::Function(const TreeOctNode* node){ 976 int idx[DIMENSION]; 977 idx[0]=index[0]+int(node->off[0]); 978 idx[1]=index[1]+int(node->off[1]); 979 idx[2]=index[2]+int(node->off[2]); 980 value+= node->nodeData.value* Real( valueTables[idx[0]]* valueTables[idx[1]]* valueTables[idx[2]]); 981 normal.coords[0]+= node->nodeData.value* Real(dValueTables[idx[0]]* valueTables[idx[1]]* valueTables[idx[2]]); 982 normal.coords[1]+= node->nodeData.value* Real( valueTables[idx[0]]*dValueTables[idx[1]]* valueTables[idx[2]]); 983 normal.coords[2]+= node->nodeData.value* Real( valueTables[idx[0]]* valueTables[idx[1]]*dValueTables[idx[2]]); 984} 985template<int Degree> 986int Octree<Degree>::LaplacianMatrixFunction::Function(const TreeOctNode* node1,const TreeOctNode* node2){ 987 Real temp; 988 int d1=int(node1->d); 989 int x1,y1,z1; 990 x1=int(node1->off[0]); 991 y1=int(node1->off[1]); 992 z1=int(node1->off[2]); 993 int dDepth=d2-d1; 994 int d; 995 d=(x2>>dDepth)-x1; 996 if(d<0){return 0;} 997 if(!dDepth){ 998 if(!d){ 999 d=y2-y1; 1000 if(d<0){return 0;} 1001 else if(!d){ 1002 d=z2-z1; 1003 if(d<0){return 0;} 1004 } 1005 } 1006 scratch[0]=FunctionData<Degree,Real>::SymmetricIndex(index[0],x1); 1007 scratch[1]=FunctionData<Degree,Real>::SymmetricIndex(index[1],y1); 1008 scratch[2]=FunctionData<Degree,Real>::SymmetricIndex(index[2],z1); 1009 temp=ot->GetLaplacian(scratch); 1010 if(node1==node2){temp/=2;} 1011 if(fabs(temp)>EPSILON){ 1012 rowElements[elementCount].Value=temp; 1013 rowElements[elementCount].N=node1->nodeData.nodeIndex-offset; 1014 elementCount++; 1015 } 1016 return 0; 1017 } 1018 return 1; 1019} 1020 1021template<int Degree> 1022int Octree<Degree>::RestrictedLaplacianMatrixFunction::Function(const TreeOctNode* node1,const TreeOctNode* node2){ 1023 int d1,d2,off1[3],off2[3]; 1024 node1->depthAndOffset(d1,off1); 1025 node2->depthAndOffset(d2,off2); 1026 int dDepth=d2-d1; 1027 int d; 1028 d=(off2[0]>>dDepth)-off1[0]; 1029 if(d<0){return 0;} 1030 1031 if(!dDepth){ 1032 if(!d){ 1033 d=off2[1]-off1[1]; 1034 if(d<0){return 0;} 1035 else if(!d){ 1036 d=off2[2]-off1[2]; 1037 if(d<0){return 0;} 1038 } 1039 } 1040 // Since we are getting the restricted matrix, we don't want to propogate out to terms that don't contribute... 1041 if(!TreeOctNode::Overlap2(depth,offset,0.5,d1,off1,radius)){return 0;} 1042 scratch[0]=FunctionData<Degree,Real>::SymmetricIndex(index[0],BinaryNode<Real>::Index(d1,off1[0])); 1043 scratch[1]=FunctionData<Degree,Real>::SymmetricIndex(index[1],BinaryNode<Real>::Index(d1,off1[1])); 1044 scratch[2]=FunctionData<Degree,Real>::SymmetricIndex(index[2],BinaryNode<Real>::Index(d1,off1[2])); 1045 Real temp=ot->GetLaplacian(scratch); 1046 if(node1==node2){temp/=2;} 1047 if(fabs(temp)>EPSILON){ 1048 rowElements[elementCount].Value=temp; 1049 rowElements[elementCount].N=node1->nodeData.nodeIndex; 1050 elementCount++; 1051 } 1052 return 0; 1053 } 1054 return 1; 1055} 1056 1057template<int Degree> 1058void Octree<Degree>::GetMCIsoTriangles(const Real& isoValue,CoredMeshData* mesh,const int& fullDepthIso,const int& nonLinearFit){ 1059 double t; 1060 TreeOctNode* temp; 1061 1062 hash_map<long long,int> roots; 1063 hash_map<long long,std::pair<Real,Point3D<Real> > > *normalHash=new hash_map<long long,std::pair<Real,Point3D<Real> > >(); 1064 1065 SetIsoSurfaceCorners(isoValue,0,fullDepthIso); 1066 // At the point all of the corner values have been set and all nodes are valid. Now it's just a matter 1067 // of running marching cubes. 1068 1069// t=Time(); 1070 fData.setValueTables(fData.VALUE_FLAG | fData.D_VALUE_FLAG,0,postNormalSmooth); 1071 temp=tree.nextLeaf(); 1072 while(temp){ 1073 SetMCRootPositions(temp,0,isoValue,roots,NULL,*normalHash,NULL,NULL,mesh,nonLinearFit); 1074 temp=tree.nextLeaf(temp); 1075 } 1076// MemoryUsage(); 1077 1078 printf("Normal Size: %.2f MB\n",double(sizeof(Point3D<Real>)*normalHash->size())/1000000); 1079// DumpOutput("Set %d root positions in: %f\n",mesh->inCorePoints.size(),Time()-t); 1080// DumpOutput("Memory Usage: %.3f MB\n",float(MemoryUsage())); 1081 1082 fData.clearValueTables(); 1083 delete normalHash; 1084 1085// DumpOutput("Post deletion size: %.3f MB\n",float(MemoryUsage())); 1086 1087// t=Time(); 1088 1089 // Now get the iso-surfaces, running from finest nodes to coarsest in order to allow for edge propogation from 1090 // finer faces to coarser ones. 1091 temp=tree.nextLeaf(); 1092 while(temp){ 1093 GetMCIsoTriangles(temp,mesh,roots,NULL,NULL,0,0); 1094 temp=tree.nextLeaf(temp); 1095 } 1096// DumpOutput("Added triangles in: %f\n",Time()-t); 1097// DumpOutput("Memory Usage: %.3f MB\n",float(MemoryUsage())); 1098} 1099template<int Degree> 1100void Octree<Degree>::GetMCIsoTriangles(const Real& isoValue,const int& subdivideDepth,CoredMeshData* mesh,const int& fullDepthIso,const int& nonLinearFit){ 1101 TreeOctNode* temp; 1102 hash_map<long long,int> boundaryRoots,*interiorRoots; 1103 hash_map<long long,std::pair<Real,Point3D<Real> > > *boundaryNormalHash,*interiorNormalHash; 1104 std::vector<Point3D<float> >* interiorPoints; 1105 1106 int sDepth; 1107 if(subdivideDepth<=0){sDepth=0;} 1108 else{sDepth=fData.depth-subdivideDepth;} 1109 if(sDepth<0){sDepth=0;} 1110 1111 SetIsoSurfaceCorners(isoValue,sDepth,fullDepthIso); 1112 // At this point all of the corner values have been set and all nodes are valid. Now it's just a matter 1113 // of running marching cubes. 1114 1115 boundaryNormalHash=new hash_map<long long,std::pair<Real,Point3D<Real> > >(); 1116 int offSet=0; 1117 SortedTreeNodes sNodes; 1118 sNodes.set(tree,0); 1119 fData.setValueTables(fData.VALUE_FLAG | fData.D_VALUE_FLAG,0,postNormalSmooth); 1120 1121 // Set the root positions for all leaf nodes below the subdivide threshold 1122 SetBoundaryMCRootPositions(sDepth,isoValue,boundaryRoots,*boundaryNormalHash,mesh,nonLinearFit); 1123 1124 for(int i=sNodes.nodeCount[sDepth];i<sNodes.nodeCount[sDepth+1];i++){ 1125 interiorRoots=new hash_map<long long,int>(); 1126 interiorNormalHash=new hash_map<long long,std::pair<Real,Point3D<Real> > >(); 1127 interiorPoints=new std::vector<Point3D<float> >(); 1128 1129 temp=sNodes.treeNodes[i]->nextLeaf(); 1130 while(temp){ 1131 if(MarchingCubes::HasRoots(temp->nodeData.mcIndex)){ 1132 SetMCRootPositions(temp,sDepth,isoValue,boundaryRoots,interiorRoots,*boundaryNormalHash,interiorNormalHash,interiorPoints,mesh,nonLinearFit); 1133 } 1134 temp=sNodes.treeNodes[i]->nextLeaf(temp); 1135 } 1136 delete interiorNormalHash; 1137 1138 temp=sNodes.treeNodes[i]->nextLeaf(); 1139 while(temp){ 1140 GetMCIsoTriangles(temp,mesh,boundaryRoots,interiorRoots,interiorPoints,offSet,sDepth); 1141 temp=sNodes.treeNodes[i]->nextLeaf(temp); 1142 } 1143 delete interiorRoots; 1144 delete interiorPoints; 1145 offSet=mesh->outOfCorePointCount(); 1146 } 1147 delete boundaryNormalHash; 1148 1149 temp=tree.nextLeaf(); 1150 while(temp){ 1151 if(temp->depth()<sDepth){GetMCIsoTriangles(temp,mesh,boundaryRoots,NULL,NULL,0,0);} 1152 temp=tree.nextLeaf(temp); 1153 } 1154} 1155template<int Degree> 1156Real Octree<Degree>::getCenterValue(const TreeOctNode* node){ 1157 int idx[3]; 1158 Real value=0; 1159 1160 neighborKey2.getNeighbors(node); 1161 VertexData::CenterIndex(node,fData.depth,idx); 1162 idx[0]*=fData.res; 1163 idx[1]*=fData.res; 1164 idx[2]*=fData.res; 1165 for(int i=0;i<=node->depth();i++){ 1166 for(int j=0;j<3;j++){ 1167 for(int k=0;k<3;k++){ 1168 for(int l=0;l<3;l++){ 1169 const TreeOctNode* n=neighborKey2.neighbors[i].neighbors[j][k][l]; 1170 if(n){ 1171 Real temp=n->nodeData.value; 1172 value+=temp*Real( 1173 fData.valueTables[idx[0]+int(n->off[0])]* 1174 fData.valueTables[idx[1]+int(n->off[1])]* 1175 fData.valueTables[idx[2]+int(n->off[2])]); 1176 } 1177 } 1178 } 1179 } 1180 } 1181 if(node->children){ 1182 for(int i=0;i<Cube::CORNERS;i++){ 1183 int ii=Cube::AntipodalCornerIndex(i); 1184 const TreeOctNode* n=&node->children[i]; 1185 while(1){ 1186 value+=n->nodeData.value*Real( 1187 fData.valueTables[idx[0]+int(n->off[0])]* 1188 fData.valueTables[idx[1]+int(n->off[1])]* 1189 fData.valueTables[idx[2]+int(n->off[2])]); 1190 if(n->children){n=&n->children[ii];} 1191 else{break;} 1192 } 1193 } 1194 } 1195 return value; 1196} 1197template<int Degree> 1198Real Octree<Degree>::getCornerValue(const TreeOctNode* node,const int& corner){ 1199 int idx[3]; 1200 Real value=0; 1201 1202 neighborKey2.getNeighbors(node); 1203 VertexData::CornerIndex(node,corner,fData.depth,idx); 1204 idx[0]*=fData.res; 1205 idx[1]*=fData.res; 1206 idx[2]*=fData.res; 1207 for(int i=0;i<=node->depth();i++){ 1208 for(int j=0;j<3;j++){ 1209 for(int k=0;k<3;k++){ 1210 for(int l=0;l<3;l++){ 1211 const TreeOctNode* n=neighborKey2.neighbors[i].neighbors[j][k][l]; 1212 if(n){ 1213 Real temp=n->nodeData.value; 1214 value+=temp*Real( 1215 fData.valueTables[idx[0]+int(n->off[0])]* 1216 fData.valueTables[idx[1]+int(n->off[1])]* 1217 fData.valueTables[idx[2]+int(n->off[2])]); 1218 } 1219 } 1220 } 1221 } 1222 } 1223 int x,y,z,d=node->depth(); 1224 Cube::FactorCornerIndex(corner,x,y,z); 1225 for(int i=0;i<2;i++){ 1226 for(int j=0;j<2;j++){ 1227 for(int k=0;k<2;k++){ 1228 const TreeOctNode* n=neighborKey2.neighbors[d].neighbors[x+i][y+j][z+k]; 1229 if(n){ 1230 int ii=Cube::AntipodalCornerIndex(Cube::CornerIndex(i,j,k)); 1231 while(n->children){ 1232 n=&n->children[ii]; 1233 value+=n->nodeData.value*Real( 1234 fData.valueTables[idx[0]+int(n->off[0])]* 1235 fData.valueTables[idx[1]+int(n->off[1])]* 1236 fData.valueTables[idx[2]+int(n->off[2])]); 1237 } 1238 } 1239 } 1240 } 1241 } 1242 return value; 1243} 1244template<int Degree> 1245void Octree<Degree>::getCornerValueAndNormal(const TreeOctNode* node,const int& corner,Real& value,Point3D<Real>& normal){ 1246 int idx[3],index[3]; 1247 value=normal.coords[0]=normal.coords[1]=normal.coords[2]=0; 1248 1249 neighborKey2.getNeighbors(node); 1250 VertexData::CornerIndex(node,corner,fData.depth,idx); 1251 idx[0]*=fData.res; 1252 idx[1]*=fData.res; 1253 idx[2]*=fData.res; 1254 for(int i=0;i<=node->depth();i++){ 1255 for(int j=0;j<3;j++){ 1256 for(int k=0;k<3;k++){ 1257 for(int l=0;l<3;l++){ 1258 const TreeOctNode* n=neighborKey2.neighbors[i].neighbors[j][k][l]; 1259 if(n){ 1260 Real temp=n->nodeData.value; 1261 index[0]=idx[0]+int(n->off[0]); 1262 index[1]=idx[1]+int(n->off[1]); 1263 index[2]=idx[2]+int(n->off[2]); 1264 value+=temp*Real(fData.valueTables[index[0]]*fData.valueTables[index[1]]*fData.valueTables[index[2]]); 1265 normal.coords[0]+=temp*Real(fData.dValueTables[index[0]]* fData.valueTables[index[1]]* fData.valueTables[index[2]]); 1266 normal.coords[1]+=temp*Real( fData.valueTables[index[0]]*fData.dValueTables[index[1]]* fData.valueTables[index[2]]); 1267 normal.coords[2]+=temp*Real( fData.valueTables[index[0]]* fData.valueTables[index[1]]*fData.dValueTables[index[2]]); 1268 } 1269 } 1270 } 1271 } 1272 } 1273 int x,y,z,d=node->depth(); 1274 Cube::FactorCornerIndex(corner,x,y,z); 1275 for(int i=0;i<2;i++){ 1276 for(int j=0;j<2;j++){ 1277 for(int k=0;k<2;k++){ 1278 const TreeOctNode* n=neighborKey2.neighbors[d].neighbors[x+i][y+j][z+k]; 1279 if(n){ 1280 int ii=Cube::AntipodalCornerIndex(Cube::CornerIndex(i,j,k)); 1281 while(n->children){ 1282 n=&n->children[ii]; 1283 Real temp=n->nodeData.value; 1284 index[0]=idx[0]+int(n->off[0]); 1285 index[1]=idx[1]+int(n->off[1]); 1286 index[2]=idx[2]+int(n->off[2]); 1287 value+=temp*Real(fData.valueTables[index[0]]*fData.valueTables[index[1]]*fData.valueTables[index[2]]); 1288 normal.coords[0]+=temp*Real(fData.dValueTables[index[0]]* fData.valueTables[index[1]]* fData.valueTables[index[2]]); 1289 normal.coords[1]+=temp*Real( fData.valueTables[index[0]]*fData.dValueTables[index[1]]* fData.valueTables[index[2]]); 1290 normal.coords[2]+=temp*Real( fData.valueTables[index[0]]* fData.valueTables[index[1]]*fData.dValueTables[index[2]]); 1291 } 1292 } 1293 } 1294 } 1295 } 1296} 1297template<int Degree> 1298Real Octree<Degree>::GetIsoValue(void){ 1299 if(this->width<=3){ 1300 TreeOctNode* temp; 1301 Real isoValue,weightSum,w; 1302 1303 neighborKey2.set(fData.depth); 1304 fData.setValueTables(fData.VALUE_FLAG,0); 1305 1306 isoValue=weightSum=0; 1307 temp=tree.nextNode(); 1308 while(temp){ 1309 w=temp->nodeData.centerWeightContribution; 1310 if(w>EPSILON){ 1311 isoValue+=getCenterValue(temp)*w; 1312 weightSum+=w; 1313 } 1314 temp=tree.nextNode(temp); 1315 } 1316 return isoValue/weightSum; 1317 } 1318 else{ 1319 const TreeOctNode* temp; 1320 Real isoValue,weightSum,w; 1321 Real myRadius; 1322 PointIndexValueFunction cf; 1323 1324 fData.setValueTables(fData.VALUE_FLAG,0); 1325 cf.valueTables=fData.valueTables; 1326 cf.res2=fData.res2; 1327 myRadius=radius; 1328 isoValue=weightSum=0; 1329 temp=tree.nextNode(); 1330 while(temp){ 1331 w=temp->nodeData.centerWeightContribution; 1332 if(w>EPSILON){ 1333 cf.value=0; 1334 int idx[3]; 1335 VertexData::CenterIndex(temp,fData.depth,idx); 1336 cf.index[0]=idx[0]*fData.res; 1337 cf.index[1]=idx[1]*fData.res; 1338 cf.index[2]=idx[2]*fData.res; 1339 TreeOctNode::ProcessPointAdjacentNodes(fData.depth,idx,&tree,width,&cf); 1340 isoValue+=cf.value*w; 1341 weightSum+=w; 1342 } 1343 temp=tree.nextNode(temp); 1344 } 1345 return isoValue/weightSum; 1346 } 1347} 1348template<int Degree> 1349void Octree<Degree>::SetIsoSurfaceCorners(const Real& isoValue,const int& subdivideDepth,const int& fullDepthIso){ 1350// double t=Time(); 1351 int i,j; 1352 hash_map<long long,Real> values; 1353 Real cornerValues[Cube::CORNERS]; 1354 PointIndexValueFunction cf; 1355 TreeOctNode* temp; 1356 int leafCount=tree.leaves(); 1357 long long key; 1358 SortedTreeNodes *sNodes=new SortedTreeNodes(); 1359 sNodes->set(tree,0); 1360 temp=tree.nextNode(); 1361 while(temp){ 1362 temp->nodeData.mcIndex=0; 1363 temp=tree.nextNode(temp); 1364 } 1365 TreeNodeData::UseIndex=0; 1366 // Start by setting the corner values of all the nodes 1367 cf.valueTables=fData.valueTables; 1368 cf.res2=fData.res2; 1369 for(i=0;i<sNodes->nodeCount[subdivideDepth];i++){ 1370 temp=sNodes->treeNodes[i]; 1371 if(!temp->children){ 1372 for(j=0;j<Cube::CORNERS;j++){ 1373 if(this->width<=3){cornerValues[j]=getCornerValue(temp,j);} 1374 else{ 1375 cf.value=0; 1376 int idx[3]; 1377 VertexData::CornerIndex(temp,j,fData.depth,idx); 1378 cf.index[0]=idx[0]*fData.res; 1379 cf.index[1]=idx[1]*fData.res; 1380 cf.index[2]=idx[2]*fData.res; 1381 TreeOctNode::ProcessPointAdjacentNodes(fData.depth,idx,&tree,width,&cf); 1382 cornerValues[j]=cf.value; 1383 } 1384 } 1385 temp->nodeData.mcIndex=MarchingCubes::GetIndex(cornerValues,isoValue); 1386 1387 if(temp->parent){ 1388 TreeOctNode* parent=temp->parent; 1389 int c=int(temp-temp->parent->children); 1390 int mcid=temp->nodeData.mcIndex&(1<<MarchingCubes::cornerMap[c]); 1391 1392 if(mcid){ 1393 parent->nodeData.mcIndex|=mcid; 1394 while(1){ 1395 if(parent->parent && (parent-parent->parent->children)==c){ 1396 parent->parent->nodeData.mcIndex|=mcid; 1397 parent=parent->parent; 1398 } 1399 else{break;} 1400 } 1401 } 1402 } 1403 } 1404 } 1405 1406// MemoryUsage(); 1407 1408 for(i=sNodes->nodeCount[subdivideDepth];i<sNodes->nodeCount[subdivideDepth+1];i++){ 1409 temp=sNodes->treeNodes[i]->nextLeaf(); 1410 while(temp){ 1411 for(j=0;j<Cube::CORNERS;j++){ 1412 int idx[3]; 1413 key=VertexData::CornerIndex(temp,j,fData.depth,idx); 1414 cf.index[0]=idx[0]*fData.res; 1415 cf.index[1]=idx[1]*fData.res; 1416 cf.index[2]=idx[2]*fData.res; 1417 if(values.find(key)!=values.end()){cornerValues[j]=values[key];} 1418 else{ 1419 if(this->width<=3){values[key]=cornerValues[j]=getCornerValue(temp,j);} 1420 else{ 1421 cf.value=0; 1422 TreeOctNode::ProcessPointAdjacentNodes(fData.depth,idx,&tree,width,&cf); 1423 values[key]=cf.value; 1424 cornerValues[j]=cf.value; 1425 } 1426 } 1427 } 1428 temp->nodeData.mcIndex=MarchingCubes::GetIndex(cornerValues,isoValue); 1429 1430 if(temp->parent){ 1431 TreeOctNode* parent=temp->parent; 1432 int c=int(temp-temp->parent->children); 1433 int mcid=temp->nodeData.mcIndex&(1<<MarchingCubes::cornerMap[c]); 1434 1435 if(mcid){ 1436 parent->nodeData.mcIndex|=mcid; 1437 while(1){ 1438 if(parent->parent && (parent-parent->parent->children)==c){ 1439 parent->parent->nodeData.mcIndex|=mcid; 1440 parent=parent->parent; 1441 } 1442 else{break;} 1443 } 1444 } 1445 } 1446 1447 temp=sNodes->treeNodes[i]->nextLeaf(temp); 1448 } 1449// MemoryUsage(); 1450 values.clear(); 1451 } 1452 delete sNodes; 1453// DumpOutput("Set corner values in: %f\n",Time()-t); 1454// DumpOutput("Memory Usage: %.3f MB\n",float(MemoryUsage())); 1455 1456 if(subdivideDepth){PreValidate(isoValue,fData.depth,subdivideDepth);} 1457} 1458template<int Degree> 1459void Octree<Degree>::Subdivide(TreeOctNode* node,const Real& isoValue,const int& maxDepth){ 1460 int i,j,c[4]; 1461 Real value; 1462 int cornerIndex2[Cube::CORNERS]; 1463 PointIndexValueFunction cf; 1464 cf.valueTables=fData.valueTables; 1465 cf.res2=fData.res2; 1466 node->initChildren(); 1467 // Since we are allocating blocks, it is possible that some of the memory was pre-allocated with 1468 // the wrong initialization 1469 1470 // Now set the corner values for the new children 1471 // Copy old corner values 1472 for(i=0;i<Cube::CORNERS;i++){cornerIndex2[i]=node->nodeData.mcIndex&(1<<MarchingCubes::cornerMap[i]);} 1473 // 8 of 27 corners set 1474 1475 // Set center corner 1476 cf.value=0; 1477 int idx[3]; 1478 VertexData::CenterIndex(node,maxDepth,idx); 1479 cf.index[0]=idx[0]*fData.res; 1480 cf.index[1]=idx[1]*fData.res; 1481 cf.index[2]=idx[2]*fData.res; 1482 if(this->width<=3){value=getCenterValue(node);} 1483 else{ 1484 TreeOctNode::ProcessPointAdjacentNodes(fData.depth,idx,&tree,width,&cf); 1485 value=cf.value; 1486 } 1487 if(value<isoValue){for(i=0;i<Cube::CORNERS;i++){cornerIndex2[i]|=1<<MarchingCubes::cornerMap[Cube::AntipodalCornerIndex(i)];}} 1488 // 9 of 27 set 1489 1490 // Set face corners 1491 for(i=0;i<Cube::NEIGHBORS;i++){ 1492 int dir,offset,e; 1493 Cube::FactorFaceIndex(i,dir,offset); 1494 cf.value=0; 1495 int idx[3]; 1496 VertexData::FaceIndex(node,i,maxDepth,idx); 1497 cf.index[0]=idx[0]*fData.res; 1498 cf.index[1]=idx[1]*fData.res; 1499 cf.index[2]=idx[2]*fData.res; 1500 TreeOctNode::ProcessPointAdjacentNodes(fData.depth,idx,&tree,width,&cf); 1501 value=cf.value; 1502 Cube::FaceCorners(i,c[0],c[1],c[2],c[3]); 1503 e=Cube::EdgeIndex(dir,0,0); 1504 if(value<isoValue){for(j=0;j<4;j++){cornerIndex2[c[j]]|=1<<MarchingCubes::cornerMap[Cube::EdgeReflectCornerIndex(c[j],e)];}} 1505 } 1506 // 15 of 27 set 1507 1508 // Set edge corners 1509 for(i=0;i<Cube::EDGES;i++){ 1510 int o,i1,i2,f; 1511 Cube::FactorEdgeIndex(i,o,i1,i2); 1512 cf.value=0; 1513 int idx[3]; 1514 VertexData::EdgeIndex(node,i,maxDepth,idx); 1515 cf.index[0]=idx[0]*fData.res; 1516 cf.index[1]=idx[1]*fData.res; 1517 cf.index[2]=idx[2]*fData.res; 1518 TreeOctNode::ProcessPointAdjacentNodes(fData.depth,idx,&tree,width,&cf); 1519 value=cf.value; 1520 Cube::EdgeCorners(i,c[0],c[1]); 1521 f=Cube::FaceIndex(o,0); 1522 if(value<isoValue){for(j=0;j<2;j++){cornerIndex2[c[j]]|=1<<MarchingCubes::cornerMap[Cube::FaceReflectCornerIndex(c[j],f)];}} 1523 } 1524 // 27 of 27 set 1525 1526 for(i=0;i<Cube::CORNERS;i++){node->children[i].nodeData.mcIndex=cornerIndex2[i];} 1527} 1528 1529template<int Degree> 1530int Octree<Degree>::InteriorFaceRootCount(const TreeOctNode* node,const int &faceIndex,const int& maxDepth){ 1531 int c1,c2,e1,e2,dir,off,cnt=0; 1532 int corners[Cube::CORNERS/2]; 1533 if(node->children){ 1534 Cube::FaceCorners(faceIndex,corners[0],corners[1],corners[2],corners[3]); 1535 Cube::FactorFaceIndex(faceIndex,dir,off); 1536 c1=corners[0]; 1537 c2=corners[3]; 1538 switch(dir){ 1539 case 0: 1540 e1=Cube::EdgeIndex(1,off,1); 1541 e2=Cube::EdgeIndex(2,off,1); 1542 break; 1543 case 1: 1544 e1=Cube::EdgeIndex(0,off,1); 1545 e2=Cube::EdgeIndex(2,1,off); 1546 break; 1547 case 2: 1548 e1=Cube::EdgeIndex(0,1,off); 1549 e2=Cube::EdgeIndex(1,1,off); 1550 break; 1551 }; 1552 cnt+=EdgeRootCount(&node->children[c1],e1,maxDepth)+EdgeRootCount(&node->children[c1],e2,maxDepth); 1553 switch(dir){ 1554 case 0: 1555 e1=Cube::EdgeIndex(1,off,0); 1556 e2=Cube::EdgeIndex(2,off,0); 1557 break; 1558 case 1: 1559 e1=Cube::EdgeIndex(0,off,0); 1560 e2=Cube::EdgeIndex(2,0,off); 1561 break; 1562 case 2: 1563 e1=Cube::EdgeIndex(0,0,off); 1564 e2=Cube::EdgeIndex(1,0,off); 1565 break; 1566 }; 1567 cnt+=EdgeRootCount(&node->children[c2],e1,maxDepth)+EdgeRootCount(&node->children[c2],e2,maxDepth); 1568 for(int i=0;i<Cube::CORNERS/2;i++){if(node->children[corners[i]].children){cnt+=InteriorFaceRootCount(&node->children[corners[i]],faceIndex,maxDepth);}} 1569 } 1570 return cnt; 1571} 1572 1573template<int Degree> 1574int Octree<Degree>::EdgeRootCount(const TreeOctNode* node,const int& edgeIndex,const int& maxDepth){ 1575 int f1,f2,c1,c2; 1576 const TreeOctNode* temp; 1577 Cube::FacesAdjacentToEdge(edgeIndex,f1,f2); 1578 1579 int eIndex; 1580 const TreeOctNode* finest=node; 1581 eIndex=edgeIndex; 1582 if(node->depth()<maxDepth){ 1583 temp=node->faceNeighbor(f1); 1584 if(temp && temp->children){ 1585 finest=temp; 1586 eIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f1); 1587 } 1588 else{ 1589 temp=node->faceNeighbor(f2); 1590 if(temp && temp->children){ 1591 finest=temp; 1592 eIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f2); 1593 } 1594 else{ 1595 temp=node->edgeNeighbor(edgeIndex); 1596 if(temp && temp->children){ 1597 finest=temp; 1598 eIndex=Cube::EdgeReflectEdgeIndex(edgeIndex); 1599 } 1600 } 1601 } 1602 } 1603 1604 Cube::EdgeCorners(eIndex,c1,c2); 1605 if(finest->children){return EdgeRootCount(&finest->children[c1],eIndex,maxDepth)+EdgeRootCount(&finest->children[c2],eIndex,maxDepth);} 1606 else{return MarchingCubes::HasEdgeRoots(finest->nodeData.mcIndex,eIndex);} 1607} 1608template<int Degree> 1609int Octree<Degree>::IsBoundaryFace(const TreeOctNode* node,const int& faceIndex,const int& subdivideDepth){ 1610 int dir,offset,d,o[3],idx; 1611 1612 if(subdivideDepth<0){return 0;} 1613 if(node->d<=subdivideDepth){return 1;} 1614 Cube::FactorFaceIndex(faceIndex,dir,offset); 1615 node->depthAndOffset(d,o); 1616 1617 idx=(int(o[dir])<<1) + (offset<<1); 1618 return !(idx%(2<<(int(node->d)-subdivideDepth))); 1619} 1620template<int Degree> 1621int Octree<Degree>::IsBoundaryEdge(const TreeOctNode* node,const int& edgeIndex,const int& subdivideDepth){ 1622 int dir,x,y; 1623 Cube::FactorEdgeIndex(edgeIndex,dir,x,y); 1624 return IsBoundaryEdge(node,dir,x,y,subdivideDepth); 1625} 1626template<int Degree> 1627int Octree<Degree>::IsBoundaryEdge(const TreeOctNode* node,const int& dir,const int& x,const int& y,const int& subdivideDepth){ 1628 int d,o[3],idx1,idx2,mask; 1629 1630 if(subdivideDepth<0){return 0;} 1631 if(node->d<=subdivideDepth){return 1;} 1632 node->depthAndOffset(d,o); 1633 1634 switch(dir){ 1635 case 0: 1636 idx1=(int(o[1])<<1) + (x<<1); 1637 idx2=(int(o[2])<<1) + (y<<1); 1638 break; 1639 case 1: 1640 idx1=(int(o[0])<<1) + (x<<1); 1641 idx2=(int(o[2])<<1) + (y<<1); 1642 break; 1643 case 2: 1644 idx1=(int(o[0])<<1) + (x<<1); 1645 idx2=(int(o[1])<<1) + (y<<1); 1646 break; 1647 } 1648 mask=2<<(int(node->d)-subdivideDepth); 1649 return !(idx1%(mask)) || !(idx2%(mask)); 1650} 1651 1652template<int Degree> 1653void Octree<Degree>::PreValidate(TreeOctNode* node,const Real& isoValue,const int& maxDepth,const int& subdivideDepth){ 1654 int sub=0; 1655 if(node->children){printf("Bad Pre-Validate\n");} 1656// if(int(node->d)<subdivideDepth){sub=1;} 1657 for(int i=0;i<Cube::NEIGHBORS && !sub;i++){ 1658 TreeOctNode* neighbor=node->faceNeighbor(i); 1659 if(neighbor && neighbor->children){ 1660 if(IsBoundaryFace(node,i,subdivideDepth) && InteriorFaceRootCount(neighbor,Cube::FaceReflectFaceIndex(i,i),maxDepth)){sub=1;} 1661 } 1662 } 1663 if(sub){ 1664 Subdivide(node,isoValue,maxDepth); 1665 for(int i=0;i<Cube::NEIGHBORS;i++){ 1666 if(IsBoundaryFace(node,i,subdivideDepth) && InteriorFaceRootCount(node,i,maxDepth)){ 1667 TreeOctNode* neighbor=node->faceNeighbor(i); 1668 while(neighbor && !neighbor->children){ 1669 PreValidate(neighbor,isoValue,maxDepth,subdivideDepth); 1670 neighbor=node->faceNeighbor(i); 1671 } 1672 } 1673 } 1674 } 1675} 1676 1677template<int Degree> 1678void Octree<Degree>::PreValidate(const Real& isoValue,const int& maxDepth,const int& subdivideDepth){ 1679 TreeOctNode* temp; 1680 1681 temp=tree.nextLeaf(); 1682 while(temp){ 1683 PreValidate(temp,isoValue,maxDepth,subdivideDepth); 1684 temp=tree.nextLeaf(temp); 1685 } 1686} 1687template<int Degree> 1688void Octree<Degree>::Validate(TreeOctNode* node,const Real& isoValue,const int& maxDepth,const int& fullDepthIso){ 1689 int i,sub=0; 1690 TreeOctNode* treeNode=node; 1691 TreeOctNode* neighbor; 1692 if(node->depth()>=maxDepth || node->children){return;} 1693 1694 // Check if full-depth extraction is enabled and we have an iso-node that is not at maximum depth 1695 if(!sub && fullDepthIso && MarchingCubes::HasRoots(node->nodeData.mcIndex)){sub=1;} 1696 1697 // Check if the node has faces that are ambiguous and are adjacent to finer neighbors 1698 for(i=0;i<Cube::NEIGHBORS && !sub;i++){ 1699 neighbor=treeNode->faceNeighbor(i); 1700 if(neighbor && neighbor->children){if(MarchingCubes::IsAmbiguous(node->nodeData.mcIndex,i)){sub=1;}} 1701 } 1702 1703 // Check if the node has edges with more than one root 1704 for(i=0;i<Cube::EDGES && !sub;i++){if(EdgeRootCount(node,i,maxDepth)>1){sub=1;}} 1705 1706 for(i=0;i<Cube::NEIGHBORS && !sub;i++){ 1707 neighbor=node->faceNeighbor(i); 1708 if( neighbor && neighbor->children && 1709 !MarchingCubes::HasFaceRoots(node->nodeData.mcIndex,i) && 1710 InteriorFaceRootCount(neighbor,Cube::FaceReflectFaceIndex(i,i),maxDepth)){sub=1;} 1711 } 1712 if(sub){ 1713 Subdivide(node,isoValue,maxDepth); 1714 for(i=0;i<Cube::NEIGHBORS;i++){ 1715 neighbor=treeNode->faceNeighbor(i); 1716 if(neighbor && !neighbor->children){Validate(neighbor,isoValue,maxDepth,fullDepthIso);} 1717 } 1718 for(i=0;i<Cube::EDGES;i++){ 1719 neighbor=treeNode->edgeNeighbor(i); 1720 if(neighbor && !neighbor->children){Validate(neighbor,isoValue,maxDepth,fullDepthIso);} 1721 } 1722 for(i=0;i<Cube::CORNERS;i++){if(!node->children[i].children){Validate(&node->children[i],isoValue,maxDepth,fullDepthIso);}} 1723 } 1724} 1725template<int Degree> 1726void Octree<Degree>::Validate(TreeOctNode* node,const Real& isoValue,const int& maxDepth,const int& fullDepthIso,const int& subdivideDepth){ 1727 int i,sub=0; 1728 TreeOctNode* treeNode=node; 1729 TreeOctNode* neighbor; 1730 if(node->depth()>=maxDepth || node->children){return;} 1731 1732 // Check if full-depth extraction is enabled and we have an iso-node that is not at maximum depth 1733 if(!sub && fullDepthIso && MarchingCubes::HasRoots(node->nodeData.mcIndex)){sub=1;} 1734 1735 // Check if the node has faces that are ambiguous and are adjacent to finer neighbors 1736 for(i=0;i<Cube::NEIGHBORS && !sub;i++){ 1737 neighbor=treeNode->faceNeighbor(i); 1738 if(neighbor && neighbor->children){if(MarchingCubes::IsAmbiguous(node->nodeData.mcIndex,i) || IsBoundaryFace(node,i,subdivideDepth)){sub=1;}} 1739 } 1740 1741 // Check if the node has edges with more than one root 1742 for(i=0;i<Cube::EDGES && !sub;i++){if(EdgeRootCount(node,i,maxDepth)>1){sub=1;}} 1743 1744 for(i=0;i<Cube::NEIGHBORS && !sub;i++){ 1745 neighbor=node->faceNeighbor(i); 1746 if( neighbor && neighbor->children && !MarchingCubes::HasFaceRoots(node->nodeData.mcIndex,i) && 1747 InteriorFaceRootCount(neighbor,Cube::FaceReflectFaceIndex(i,i),maxDepth)){sub=1;} 1748 } 1749 if(sub){ 1750 Subdivide(node,isoValue,maxDepth); 1751 for(i=0;i<Cube::NEIGHBORS;i++){ 1752 neighbor=treeNode->faceNeighbor(i); 1753 if(neighbor && !neighbor->children){Validate(neighbor,isoValue,maxDepth,fullDepthIso,subdivideDepth);} 1754 } 1755 for(i=0;i<Cube::EDGES;i++){ 1756 neighbor=treeNode->edgeNeighbor(i); 1757 if(neighbor && !neighbor->children){Validate(neighbor,isoValue,maxDepth,fullDepthIso,subdivideDepth);} 1758 } 1759 for(i=0;i<Cube::CORNERS;i++){if(!node->children[i].children){Validate(&node->children[i],isoValue,maxDepth,fullDepthIso,subdivideDepth);}} 1760 } 1761} 1762////////////////////////////////////////////////////////////////////////////////////// 1763// The assumption made when calling this code is that the edge has at most one root // 1764////////////////////////////////////////////////////////////////////////////////////// 1765template<int Degree> 1766int Octree<Degree>::GetRoot(const RootInfo& ri,const Real& isoValue,Point3D<Real> & position,hash_map<long long,std::pair<Real,Point3D<Real> > >& normalHash,const int& nonLinearFit){ 1767 int c1,c2; 1768 Cube::EdgeCorners(ri.edgeIndex,c1,c2); 1769 if(!MarchingCubes::HasEdgeRoots(ri.node->nodeData.mcIndex,ri.edgeIndex)){return 0;} 1770 1771 long long key; 1772 Point3D<Real> n[2]; 1773 PointIndexValueAndNormalFunction cnf; 1774 cnf.valueTables=fData.valueTables; 1775 cnf.dValueTables=fData.dValueTables; 1776 cnf.res2=fData.res2; 1777 1778 int i,o,i1,i2,rCount=0; 1779 Polynomial<2> P; 1780 std::vector<double> roots; 1781 double x0,x1; 1782 Real center,width; 1783 Real averageRoot=0; 1784 Cube::FactorEdgeIndex(ri.edgeIndex,o,i1,i2); 1785 int idx[3]; 1786 key=VertexData::CornerIndex(ri.node,c1,fData.depth,idx); 1787 cnf.index[0]=idx[0]*fData.res; 1788 cnf.index[1]=idx[1]*fData.res; 1789 cnf.index[2]=idx[2]*fData.res; 1790 1791 if(normalHash.find(key)==normalHash.end()){ 1792 cnf.value=0; 1793 cnf.normal.coords[0]=cnf.normal.coords[1]=cnf.normal.coords[2]=0; 1794 // Careful here as the normal isn't quite accurate... (i.e. postNormalSmooth is ignored) 1795 if(this->width<=3){getCornerValueAndNormal(ri.node,c1,cnf.value,cnf.normal);} 1796 else{TreeOctNode::ProcessPointAdjacentNodes(fData.depth,idx,&tree,this->width,&cnf);} 1797 normalHash[key]=std::pair<Real,Point3D<Real> >(cnf.value,cnf.normal); 1798 } 1799 x0=normalHash[key].first; 1800 n[0]=normalHash[key].second; 1801 1802 key=VertexData::CornerIndex(ri.node,c2,fData.depth,idx); 1803 cnf.index[0]=idx[0]*fData.res; 1804 cnf.index[1]=idx[1]*fData.res; 1805 cnf.index[2]=idx[2]*fData.res; 1806 if(normalHash.find(key)==normalHash.end()){ 1807 cnf.value=0; 1808 cnf.normal.coords[0]=cnf.normal.coords[1]=cnf.normal.coords[2]=0; 1809 if(this->width<=3){getCornerValueAndNormal(ri.node,c2,cnf.value,cnf.normal);} 1810 else{TreeOctNode::ProcessPointAdjacentNodes(fData.depth,idx,&tree,this->width,&cnf);} 1811 normalHash[key]=std::pair<Real,Point3D<Real> >(cnf.value,cnf.normal); 1812 } 1813 x1=normalHash[key].first; 1814 n[1]=normalHash[key].second; 1815 1816 Point3D<Real> c; 1817 ri.node->centerAndWidth(c,width); 1818 center=c.coords[o]; 1819 for(i=0;i<DIMENSION;i++){ 1820 n[0].coords[i]*=width; 1821 n[1].coords[i]*=width; 1822 } 1823 1824 switch(o){ 1825 case 0: 1826 position.coords[1]=c.coords[1]-width/2+width*i1; 1827 position.coords[2]=c.coords[2]-width/2+width*i2; 1828 break; 1829 case 1: 1830 position.coords[0]=c.coords[0]-width/2+width*i1; 1831 position.coords[2]=c.coords[2]-width/2+width*i2; 1832 break; 1833 case 2: 1834 position.coords[0]=c.coords[0]-width/2+width*i1; 1835 position.coords[1]=c.coords[1]-width/2+width*i2; 1836 break; 1837 } 1838 double dx0,dx1; 1839 dx0=n[0].coords[o]; 1840 dx1=n[1].coords[o]; 1841 1842 // The scaling will turn the Hermite Spline into a quadratic 1843 double scl=(x1-x0)/((dx1+dx0)/2); 1844 dx0*=scl; 1845 dx1*=scl; 1846 1847 // Hermite Spline 1848 P.coefficients[0]=x0; 1849 P.coefficients[1]=dx0; 1850 P.coefficients[2]=3*(x1-x0)-dx1-2*dx0; 1851 1852 P.getSolutions(isoValue,roots,EPSILON); 1853 for(i=0;i<int(roots.size());i++){ 1854 if(roots[i]>=0 && roots[i]<=1){ 1855 averageRoot+=Real(roots[i]); 1856 rCount++; 1857 } 1858 } 1859 if(rCount && nonLinearFit) {averageRoot/=rCount;} 1860 else {averageRoot=Real((x0-isoValue)/(x0-x1));} 1861 1862 position.coords[o]=Real(center-width/2+width*averageRoot); 1863 return 1; 1864} 1865 1866template<int Degree> 1867int Octree<Degree>::GetRoot(const RootInfo& ri,const Real& isoValue,const int& maxDepth,Point3D<Real>& position,hash_map<long long,std::pair<Real,Point3D<Real> > >& normals, 1868 Point3D<Real>* normal,const int& nonLinearFit){ 1869 if(!MarchingCubes::HasRoots(ri.node->nodeData.mcIndex)){return 0;} 1870 return GetRoot(ri,isoValue,position,normals,nonLinearFit); 1871} 1872template<int Degree> 1873int Octree<Degree>::GetRootIndex(const TreeOctNode* node,const int& edgeIndex,const int& maxDepth,const int& sDepth,RootInfo& ri){ 1874 int c1,c2,f1,f2; 1875 const TreeOctNode *temp,*finest; 1876 int finestIndex; 1877 1878 Cube::FacesAdjacentToEdge(edgeIndex,f1,f2); 1879 1880 finest=node; 1881 finestIndex=edgeIndex; 1882 if(node->depth()<maxDepth){ 1883 if(IsBoundaryFace(node,f1,sDepth)){temp=NULL;} 1884 else{temp=node->faceNeighbor(f1);} 1885 if(temp && temp->children){ 1886 finest=temp; 1887 finestIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f1); 1888 } 1889 else{ 1890 if(IsBoundaryFace(node,f2,sDepth)){temp=NULL;} 1891 else{temp=node->faceNeighbor(f2);} 1892 if(temp && temp->children){ 1893 finest=temp; 1894 finestIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f2); 1895 } 1896 else{ 1897 if(IsBoundaryEdge(node,edgeIndex,sDepth)){temp=NULL;} 1898 else{temp=node->edgeNeighbor(edgeIndex);} 1899 if(temp && temp->children){ 1900 finest=temp; 1901 finestIndex=Cube::EdgeReflectEdgeIndex(edgeIndex); 1902 } 1903 } 1904 } 1905 } 1906 1907 Cube::EdgeCorners(finestIndex,c1,c2); 1908 if(finest->children){ 1909 if (GetRootIndex(&finest->children[c1],finestIndex,maxDepth,sDepth,ri)) {return 1;} 1910 else if (GetRootIndex(&finest->children[c2],finestIndex,maxDepth,sDepth,ri)) {return 1;} 1911 else {return 0;} 1912 } 1913 else{ 1914 if(!(MarchingCubes::edgeMask[finest->nodeData.mcIndex] & (1<<finestIndex))){return 0;} 1915 1916 int o,i1,i2; 1917 Cube::FactorEdgeIndex(finestIndex,o,i1,i2); 1918 int d,off[3]; 1919 finest->depthAndOffset(d,off); 1920 ri.node=finest; 1921 ri.edgeIndex=finestIndex; 1922 int eIndex[2],offset; 1923 offset=BinaryNode<Real>::Index(d,off[o]); 1924 switch(o){ 1925 case 0: 1926 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i1); 1927 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2); 1928 break; 1929 case 1: 1930 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1); 1931 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2); 1932 break; 1933 case 2: 1934 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1); 1935 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i2); 1936 break; 1937 } 1938 ri.key = (long long)(o) | (long long)(eIndex[0])<<5 | (long long)(eIndex[1])<<25 | (long long)(offset)<<45; 1939 return 1; 1940 } 1941} 1942template<int Degree> 1943int Octree<Degree>::GetRootIndex(const TreeOctNode* node,const int& edgeIndex,const int& maxDepth,RootInfo& ri){ 1944 int c1,c2,f1,f2; 1945 const TreeOctNode *temp,*finest; 1946 int finestIndex; 1947 1948 1949 // The assumption is that the super-edge has a root along it. 1950 if(!(MarchingCubes::edgeMask[node->nodeData.mcIndex] & (1<<edgeIndex))){return 0;} 1951 1952 Cube::FacesAdjacentToEdge(edgeIndex,f1,f2); 1953 1954 finest=node; 1955 finestIndex=edgeIndex; 1956 if(node->depth()<maxDepth){ 1957 temp=node->faceNeighbor(f1); 1958 if(temp && temp->children){ 1959 finest=temp; 1960 finestIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f1); 1961 } 1962 else{ 1963 temp=node->faceNeighbor(f2); 1964 if(temp && temp->children){ 1965 finest=temp; 1966 finestIndex=Cube::FaceReflectEdgeIndex(edgeIndex,f2); 1967 } 1968 else{ 1969 temp=node->edgeNeighbor(edgeIndex); 1970 if(temp && temp->children){ 1971 finest=temp; 1972 finestIndex=Cube::EdgeReflectEdgeIndex(edgeIndex); 1973 } 1974 } 1975 } 1976 } 1977 1978 Cube::EdgeCorners(finestIndex,c1,c2); 1979 if(finest->children){ 1980 if (GetRootIndex(&finest->children[c1],finestIndex,maxDepth,ri)) {return 1;} 1981 else if (GetRootIndex(&finest->children[c2],finestIndex,maxDepth,ri)) {return 1;} 1982 else {return 0;} 1983 } 1984 else{ 1985 int o,i1,i2; 1986 Cube::FactorEdgeIndex(finestIndex,o,i1,i2); 1987 int d,off[3]; 1988 finest->depthAndOffset(d,off); 1989 ri.node=finest; 1990 ri.edgeIndex=finestIndex; 1991 int offset,eIndex[2]; 1992 offset=BinaryNode<Real>::Index(d,off[o]); 1993 switch(o){ 1994 case 0: 1995 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i1); 1996 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2); 1997 break; 1998 case 1: 1999 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1); 2000 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2); 2001 break; 2002 case 2: 2003 eIndex[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1); 2004 eIndex[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i2); 2005 break; 2006 } 2007 ri.key= (long long)(o) | (long long)(eIndex[0])<<5 | (long long)(eIndex[1])<<25 | (long long)(offset)<<45; 2008 return 1; 2009 } 2010} 2011template<int Degree> 2012int Octree<Degree>::GetRootPair(const RootInfo& ri,const int& maxDepth,RootInfo& pair){ 2013 const TreeOctNode* node=ri.node; 2014 int c1,c2,c; 2015 Cube::EdgeCorners(ri.edgeIndex,c1,c2); 2016 while(node->parent){ 2017 c=int(node-node->parent->children); 2018 if(c!=c1 && c!=c2){return 0;} 2019 if(!MarchingCubes::HasEdgeRoots(node->parent->nodeData.mcIndex,ri.edgeIndex)){ 2020 if(c==c1){return GetRootIndex(&node->parent->children[c2],ri.edgeIndex,maxDepth,pair);} 2021 else{return GetRootIndex(&node->parent->children[c1],ri.edgeIndex,maxDepth,pair);} 2022 } 2023 node=node->parent; 2024 } 2025 return 0; 2026 2027} 2028template<int Degree> 2029int Octree<Degree>::GetRootIndex(const long long& key,hash_map<long long,int>& boundaryRoots,hash_map<long long,int>* interiorRoots,CoredPointIndex& index){ 2030 hash_map<long long,int>::iterator rootIter=boundaryRoots.find(key); 2031 if(rootIter!=boundaryRoots.end()){ 2032 index.inCore=1; 2033 index.index=rootIter->second; 2034 return 1; 2035 } 2036 else if(interiorRoots){ 2037 rootIter=interiorRoots->find(key); 2038 if(rootIter!=interiorRoots->end()){ 2039 index.inCore=0; 2040 index.index=rootIter->second; 2041 return 1; 2042 } 2043 } 2044 return 0; 2045} 2046template<int Degree> 2047int Octree<Degree>::SetMCRootPositions(TreeOctNode* node,const int& sDepth,const Real& isoValue, 2048 hash_map<long long,int>& boundaryRoots,hash_map<long long,int>* interiorRoots, 2049 hash_map<long long,std::pair<Real,Point3D<Real> > >& boundaryNormalHash,hash_map<long long,std::pair<Real,Point3D<Real> > >* interiorNormalHash, 2050 std::vector<Point3D<float> >* interiorPositions, 2051 CoredMeshData* mesh,const int& nonLinearFit){ 2052 Point3D<Real> position; 2053 int i,j,k,eIndex; 2054 RootInfo ri; 2055 int count=0; 2056 if(!MarchingCubes::HasRoots(node->nodeData.mcIndex)){return 0;} 2057 for(i=0;i<DIMENSION;i++){ 2058 for(j=0;j<2;j++){ 2059 for(k=0;k<2;k++){ 2060 long long key; 2061 eIndex=Cube::EdgeIndex(i,j,k); 2062 if(GetRootIndex(node,eIndex,fData.depth,ri)){ 2063 key=ri.key; 2064 if(!interiorRoots || IsBoundaryEdge(node,i,j,k,sDepth)){ 2065 if(boundaryRoots.find(key)==boundaryRoots.end()){ 2066 GetRoot(ri,isoValue,fData.depth,position,boundaryNormalHash,NULL,nonLinearFit); 2067 mesh->inCorePoints.push_back(position); 2068 boundaryRoots[key]=int(mesh->inCorePoints.size())-1; 2069 count++; 2070 } 2071 } 2072 else{ 2073 if(interiorRoots->find(key)==interiorRoots->end()){ 2074 GetRoot(ri,isoValue,fData.depth,position,*interiorNormalHash,NULL,nonLinearFit); 2075 (*interiorRoots)[key]=mesh->addOutOfCorePoint(position); 2076 interiorPositions->push_back(position); 2077 count++; 2078 } 2079 } 2080 } 2081 } 2082 } 2083 } 2084 return count; 2085} 2086template<int Degree> 2087int Octree<Degree>::SetBoundaryMCRootPositions(const int& sDepth,const Real& isoValue, 2088 hash_map<long long,int>& boundaryRoots,hash_map<long long,std::pair<Real,Point3D<Real> > >& boundaryNormalHash, 2089 CoredMeshData* mesh,const int& nonLinearFit){ 2090 Point3D<Real> position; 2091 int i,j,k,eIndex,hits; 2092 RootInfo ri; 2093 int count=0; 2094 TreeOctNode* node; 2095 2096 node=tree.nextLeaf(); 2097 while(node){ 2098 if(MarchingCubes::HasRoots(node->nodeData.mcIndex)){ 2099 hits=0; 2100 for(i=0;i<DIMENSION;i++){ 2101 for(j=0;j<2;j++){ 2102 for(k=0;k<2;k++){ 2103 if(IsBoundaryEdge(node,i,j,k,sDepth)){ 2104 hits++; 2105 long long key; 2106 eIndex=Cube::EdgeIndex(i,j,k); 2107 if(GetRootIndex(node,eIndex,fData.depth,ri)){ 2108 key=ri.key; 2109 if(boundaryRoots.find(key)==boundaryRoots.end()){ 2110 GetRoot(ri,isoValue,fData.depth,position,boundaryNormalHash,NULL,nonLinearFit); 2111 mesh->inCorePoints.push_back(position); 2112 boundaryRoots[key]=int(mesh->inCorePoints.size())-1; 2113 count++; 2114 } 2115 } 2116 } 2117 } 2118 } 2119 } 2120 } 2121 if(hits){node=tree.nextLeaf(node);} 2122 else{node=tree.nextBranch(node);} 2123 } 2124 return count; 2125} 2126template<int Degree> 2127void Octree<Degree>::GetMCIsoEdges(TreeOctNode* node,hash_map<long long,int>& boundaryRoots,hash_map<long long,int>* interiorRoots,const int& sDepth, 2128 std::vector<std::pair<long long,long long> >& edges){ 2129 TreeOctNode* temp; 2130 int count=0,tris=0; 2131 int isoTri[DIMENSION*MarchingCubes::MAX_TRIANGLES]; 2132 FaceEdgesFunction fef; 2133 int ref,fIndex; 2134 hash_map<long long,std::pair<RootInfo,int> >::iterator iter; 2135 hash_map<long long,std::pair<RootInfo,int> > vertexCount; 2136 2137 fef.edges=&edges; 2138 fef.maxDepth=fData.depth; 2139 fef.vertexCount=&vertexCount; 2140 count=MarchingCubes::AddTriangleIndices(node->nodeData.mcIndex,isoTri); 2141 for(fIndex=0;fIndex<Cube::NEIGHBORS;fIndex++){ 2142 ref=Cube::FaceReflectFaceIndex(fIndex,fIndex); 2143 fef.fIndex=ref; 2144 temp=node->faceNeighbor(fIndex); 2145 // If the face neighbor exists and has higher resolution than the current node, 2146 // get the iso-curve from the neighbor 2147 if(temp && temp->children && !IsBoundaryFace(node,fIndex,sDepth)){temp->processNodeFaces(temp,&fef,ref);} 2148 // Otherwise, get it from the node 2149 else{ 2150 RootInfo ri1,ri2; 2151 for(int j=0;j<count;j++){ 2152 for(int k=0;k<3;k++){ 2153 if(fIndex==Cube::FaceAdjacentToEdges(isoTri[j*3+k],isoTri[j*3+((k+1)%3)])){ 2154 if(GetRootIndex(node,isoTri[j*3+k],fData.depth,ri1) && GetRootIndex(node,isoTri[j*3+((k+1)%3)],fData.depth,ri2)){ 2155 edges.push_back(std::pair<long long,long long>(ri1.key,ri2.key)); 2156 iter=vertexCount.find(ri1.key); 2157 if(iter==vertexCount.end()){ 2158 vertexCount[ri1.key].first=ri1; 2159 vertexCount[ri1.key].second=0; 2160 } 2161 iter=vertexCount.find(ri2.key); 2162 if(iter==vertexCount.end()){ 2163 vertexCount[ri2.key].first=ri2; 2164 vertexCount[ri2.key].second=0; 2165 } 2166 vertexCount[ri1.key].second++; 2167 vertexCount[ri2.key].second--; 2168 } 2169 else{fprintf(stderr,"Bad Edge 1: %d %d\n",ri1.key,ri2.key);} 2170 } 2171 } 2172 } 2173 } 2174 } 2175 for(int i=0;i<int(edges.size());i++){ 2176 iter=vertexCount.find(edges[i].first); 2177 if(iter==vertexCount.end()){printf("Could not find vertex: %lld\n",edges[i].first);} 2178 else if(vertexCount[edges[i].first].second){ 2179 RootInfo ri; 2180 GetRootPair(vertexCount[edges[i].first].first,fData.depth,ri); 2181 iter=vertexCount.find(ri.key); 2182 if(iter==vertexCount.end()){printf("Vertex pair not in list\n");} 2183 else{ 2184 edges.push_back(std::pair<long long,long long>(ri.key,edges[i].first)); 2185 vertexCount[ri.key].second++; 2186 vertexCount[edges[i].first].second--; 2187 } 2188 } 2189 2190 iter=vertexCount.find(edges[i].second); 2191 if(iter==vertexCount.end()){printf("Could not find vertex: %lld\n",edges[i].second);} 2192 else if(vertexCount[edges[i].second].second){ 2193 RootInfo ri; 2194 GetRootPair(vertexCount[edges[i].second].first,fData.depth,ri); 2195 iter=vertexCount.find(ri.key); 2196 if(iter==vertexCount.end()){printf("Vertex pair not in list\n");} 2197 else{ 2198 edges.push_back(std::pair<long long,long long>(edges[i].second,ri.key)); 2199 vertexCount[edges[i].second].second++; 2200 vertexCount[ri.key].second--; 2201 } 2202 } 2203 } 2204} 2205template<int Degree> 2206int Octree<Degree>::GetMCIsoTriangles(TreeOctNode* node,CoredMeshData* mesh,hash_map<long long,int>& boundaryRoots, 2207 hash_map<long long,int>* interiorRoots,std::vector<Point3D<float> >* interiorPositions,const int& offSet,const int& sDepth) 2208{ 2209 int tris=0; 2210 std::vector<std::pair<long long,long long> > edges; 2211 std::vector<std::vector<std::pair<long long,long long> > > edgeLoops; 2212 GetMCIsoEdges(node,boundaryRoots,interiorRoots,sDepth,edges); 2213 2214 GetEdgeLoops(edges,edgeLoops); 2215 for(int i=0;i<int(edgeLoops.size());i++){ 2216 CoredPointIndex p; 2217 std::vector<CoredPointIndex> edgeIndices; 2218 for(int j=0;j<int(edgeLoops[i].size());j++){ 2219 if(!GetRootIndex(edgeLoops[i][j].first,boundaryRoots,interiorRoots,p)){printf("Bad Point Index\n");} 2220 else{edgeIndices.push_back(p);} 2221 } 2222 tris+=AddTriangles(mesh,edgeIndices,interiorPositions,offSet); 2223 } 2224 return tris; 2225} 2226 2227template<int Degree> 2228int Octree<Degree>::GetEdgeLoops(std::vector<std::pair<long long,long long> >& edges,std::vector<std::vector<std::pair<long long,long long> > >& loops){ 2229 int loopSize=0; 2230 long long frontIdx,backIdx; 2231 std::pair<long long,long long> e,temp; 2232 loops.clear(); 2233 2234 while(edges.size()){ 2235 std::vector<std::pair<long long,long long> > front,back; 2236 e=edges[0]; 2237 loops.resize(loopSize+1); 2238 edges[0]=edges[edges.size()-1]; 2239 edges.pop_back(); 2240 frontIdx=e.second; 2241 backIdx=e.first; 2242 for(int j=int(edges.size())-1;j>=0;j--){ 2243 if(edges[j].first==frontIdx || edges[j].second==frontIdx){ 2244 if(edges[j].first==frontIdx) {temp=edges[j];} 2245 else {temp.first=edges[j].second;temp.second=edges[j].first;} 2246 frontIdx=temp.second; 2247 front.push_back(temp); 2248 edges[j]=edges[edges.size()-1]; 2249 edges.pop_back(); 2250 j=int(edges.size()); 2251 } 2252 else if(edges[j].first==backIdx || edges[j].second==backIdx){ 2253 if(edges[j].second==backIdx) {temp=edges[j];} 2254 else {temp.first=edges[j].second;temp.second=edges[j].first;} 2255 backIdx=temp.first; 2256 back.push_back(temp); 2257 edges[j]=edges[edges.size()-1]; 2258 edges.pop_back(); 2259 j=int(edges.size()); 2260 } 2261 } 2262 for(int j=int(back.size())-1;j>=0;j--){loops[loopSize].push_back(back[j]);} 2263 loops[loopSize].push_back(e); 2264 for(int j=0;j<int(front.size());j++){loops[loopSize].push_back(front[j]);} 2265 loopSize++; 2266 } 2267 return int(loops.size()); 2268} 2269template<int Degree> 2270int Octree<Degree>::AddTriangles(CoredMeshData* mesh,std::vector<CoredPointIndex> edges[3],std::vector<Point3D<float> >* interiorPositions,const int& offSet){ 2271 std::vector<CoredPointIndex> e; 2272 for(int i=0;i<3;i++){for(size_t j=0;j<edges[i].size();j++){e.push_back(edges[i][j]);}} 2273 return AddTriangles(mesh,e,interiorPositions,offSet); 2274} 2275template<int Degree> 2276int Octree<Degree>::AddTriangles(CoredMeshData* mesh,std::vector<CoredPointIndex>& edges,std::vector<Point3D<float> >* interiorPositions,const int& offSet){ 2277 if(edges.size()>3){ 2278 Triangulation<float> t; 2279 2280 // Add the points to the triangulation 2281 for(int i=0;i<int(edges.size());i++){ 2282 Point3D<Real> p; 2283 if(edges[i].inCore) {for(int j=0;j<3;j++){p.coords[j]=mesh->inCorePoints[edges[i].index].coords[j];}} 2284 else {for(int j=0;j<3;j++){p.coords[j]=(*interiorPositions)[edges[i].index-offSet].coords[j];}} 2285 t.points.push_back(p); 2286 } 2287 2288 // Create a fan triangulation 2289 for(int i=1;i<int(edges.size())-1;i++){t.addTriangle(0,i,i+1);} 2290 2291 // Minimize 2292 while(1){ 2293 int i; 2294 for(i=0;i<int(t.edges.size());i++){if(t.flipMinimize(i)){break;}} 2295 if(i==t.edges.size()){break;} 2296 } 2297 // Add the triangles to the mesh 2298 for(int i=0;i<int(t.triangles.size());i++){ 2299 TriangleIndex tri; 2300 int idx[3]; 2301 int inCoreFlag=0; 2302 t.factor(i,idx[0],idx[1],idx[2]); 2303 for(int j=0;j<3;j++){ 2304 tri.idx[j]=edges[idx[j]].index; 2305 if(edges[idx[j]].inCore){inCoreFlag|=CoredMeshData::IN_CORE_FLAG[j];} 2306 } 2307 mesh->addTriangle(tri,inCoreFlag); 2308 } 2309 } 2310 else if(edges.size()==3){ 2311 TriangleIndex tri; 2312 int inCoreFlag=0; 2313 for(int i=0;i<3;i++){ 2314 tri.idx[i]=edges[i].index; 2315 if(edges[i].inCore){inCoreFlag|=CoredMeshData::IN_CORE_FLAG[i];} 2316 } 2317 mesh->addTriangle(tri,inCoreFlag); 2318 } 2319 return int(edges.size())-2; 2320} 2321//////////////// 2322// VertexData // 2323//////////////// 2324long long VertexData::CenterIndex(const TreeOctNode* node,const int& maxDepth){ 2325 int idx[DIMENSION]; 2326 return CenterIndex(node,maxDepth,idx); 2327} 2328long long VertexData::CenterIndex(const TreeOctNode* node,const int& maxDepth,int idx[DIMENSION]){ 2329 int d,o[3]; 2330 node->depthAndOffset(d,o); 2331 for(int i=0;i<DIMENSION;i++){idx[i]=BinaryNode<Real>::CornerIndex(maxDepth+1,d+1,o[i]<<1,1);} 2332 return (long long)(idx[0]) | (long long)(idx[1])<<15 | (long long)(idx[2])<<30; 2333} 2334long long VertexData::CenterIndex(const int& depth,const int offSet[DIMENSION],const int& maxDepth,int idx[DIMENSION]){ 2335 for(int i=0;i<DIMENSION;i++){idx[i]=BinaryNode<Real>::CornerIndex(maxDepth+1,depth+1,offSet[i]<<1,1);} 2336 return (long long)(idx[0]) | (long long)(idx[1])<<15 | (long long)(idx[2])<<30; 2337} 2338long long VertexData::CornerIndex(const TreeOctNode* node,const int& cIndex,const int& maxDepth){ 2339 int idx[DIMENSION]; 2340 return CornerIndex(node,cIndex,maxDepth,idx); 2341} 2342long long VertexData::CornerIndex(const TreeOctNode* node,const int& cIndex,const int& maxDepth,int idx[DIMENSION]){ 2343 int x[DIMENSION]; 2344 Cube::FactorCornerIndex(cIndex,x[0],x[1],x[2]); 2345 int d,o[3]; 2346 node->depthAndOffset(d,o); 2347 for(int i=0;i<DIMENSION;i++){idx[i]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,o[i],x[i]);} 2348 return (long long)(idx[0]) | (long long)(idx[1])<<15 | (long long)(idx[2])<<30; 2349} 2350long long VertexData::CornerIndex(const int& depth,const int offSet[DIMENSION],const int& cIndex,const int& maxDepth,int idx[DIMENSION]){ 2351 int x[DIMENSION]; 2352 Cube::FactorCornerIndex(cIndex,x[0],x[1],x[2]); 2353 for(int i=0;i<DIMENSION;i++){idx[i]=BinaryNode<Real>::CornerIndex(maxDepth+1,depth,offSet[i],x[i]);} 2354 return (long long)(idx[0]) | (long long)(idx[1])<<15 | (long long)(idx[2])<<30; 2355} 2356long long VertexData::FaceIndex(const TreeOctNode* node,const int& fIndex,const int& maxDepth){ 2357 int idx[DIMENSION]; 2358 return FaceIndex(node,fIndex,maxDepth,idx); 2359} 2360long long VertexData::FaceIndex(const TreeOctNode* node,const int& fIndex,const int& maxDepth,int idx[DIMENSION]){ 2361 int dir,offset; 2362 Cube::FactorFaceIndex(fIndex,dir,offset); 2363 int d,o[3]; 2364 node->depthAndOffset(d,o); 2365 for(int i=0;i<DIMENSION;i++){idx[i]=BinaryNode<Real>::CornerIndex(maxDepth+1,d+1,o[i]<<1,1);} 2366 idx[dir]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,o[dir],offset); 2367 return (long long)(idx[0]) | (long long)(idx[1])<<15 | (long long)(idx[2])<<30; 2368} 2369long long VertexData::EdgeIndex(const TreeOctNode* node,const int& eIndex,const int& maxDepth){ 2370 int idx[DIMENSION]; 2371 return EdgeIndex(node,eIndex,maxDepth,idx); 2372} 2373long long VertexData::EdgeIndex(const TreeOctNode* node,const int& eIndex,const int& maxDepth,int idx[DIMENSION]){ 2374 int o,i1,i2; 2375 int d,off[3]; 2376 node->depthAndOffset(d,off); 2377 for(int i=0;i<DIMENSION;i++){idx[i]=BinaryNode<Real>::CornerIndex(maxDepth+1,d+1,off[i]<<1,1);} 2378 Cube::FactorEdgeIndex(eIndex,o,i1,i2); 2379 switch(o){ 2380 case 0: 2381 idx[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i1); 2382 idx[2]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2); 2383 break; 2384 case 1: 2385 idx[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1); 2386 idx[2]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[2],i2); 2387 break; 2388 case 2: 2389 idx[0]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[0],i1); 2390 idx[1]=BinaryNode<Real>::CornerIndex(maxDepth+1,d,off[1],i2); 2391 break; 2392 }; 2393 return (long long)(idx[0]) | (long long)(idx[1])<<15 | (long long)(idx[2])<<30; 2394} 2395