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 <stdio.h> 30 31template<class Real> 32Real Random(void){return Real(rand())/RAND_MAX;} 33 34template<class Real> 35Point3D<Real> RandomBallPoint(void){ 36 Point3D<Real> p; 37 while(1){ 38 p.coords[0]=Real(1.0-2.0*Random<Real>()); 39 p.coords[1]=Real(1.0-2.0*Random<Real>()); 40 p.coords[2]=Real(1.0-2.0*Random<Real>()); 41 double l=SquareLength(p); 42 if(l<=1){return p;} 43 } 44} 45template<class Real> 46Point3D<Real> RandomSpherePoint(void){ 47 Point3D<Real> p=RandomBallPoint<Real>(); 48 Real l=Real(Length(p)); 49 p.coords[0]/=l; 50 p.coords[1]/=l; 51 p.coords[2]/=l; 52 return p; 53} 54 55template<class Real> 56double SquareLength(const Point3D<Real>& p){return p.coords[0]*p.coords[0]+p.coords[1]*p.coords[1]+p.coords[2]*p.coords[2];} 57 58template<class Real> 59double Length(const Point3D<Real>& p){return sqrt(SquareLength(p));} 60 61template<class Real> 62double SquareDistance(const Point3D<Real>& p1,const Point3D<Real>& p2){ 63 return (p1.coords[0]-p2.coords[0])*(p1.coords[0]-p2.coords[0])+(p1.coords[1]-p2.coords[1])*(p1.coords[1]-p2.coords[1])+(p1.coords[2]-p2.coords[2])*(p1.coords[2]-p2.coords[2]); 64} 65 66template<class Real> 67double Distance(const Point3D<Real>& p1,const Point3D<Real>& p2){return sqrt(SquareDistance(p1,p2));} 68 69template <class Real> 70void CrossProduct(const Point3D<Real>& p1,const Point3D<Real>& p2,Point3D<Real>& p){ 71 p.coords[0]= p1.coords[1]*p2.coords[2]-p1.coords[2]*p2.coords[1]; 72 p.coords[1]=-p1.coords[0]*p2.coords[2]+p1.coords[2]*p2.coords[0]; 73 p.coords[2]= p1.coords[0]*p2.coords[1]-p1.coords[1]*p2.coords[0]; 74} 75template<class Real> 76void EdgeCollapse(const Real& edgeRatio,std::vector<TriangleIndex>& triangles,std::vector< Point3D<Real> >& positions,std::vector< Point3D<Real> >* normals){ 77 int i,j,*remapTable,*pointCount,idx[3]; 78 Point3D<Real> p[3],q[2],c; 79 double d[3],a; 80 double Ratio=12.0/sqrt(3.0); // (Sum of Squares Length / Area) for and equilateral triangle 81 82 remapTable=new int[positions.size()]; 83 pointCount=new int[positions.size()]; 84 for(i=0;i<int(positions.size());i++){ 85 remapTable[i]=i; 86 pointCount[i]=1; 87 } 88 for(i=int(triangles.size()-1);i>=0;i--){ 89 for(j=0;j<3;j++){ 90 idx[j]=triangles[i].idx[j]; 91 while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];} 92 } 93 if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){ 94 triangles[i]=triangles[triangles.size()-1]; 95 triangles.pop_back(); 96 continue; 97 } 98 for(j=0;j<3;j++){ 99 p[j].coords[0]=positions[idx[j]].coords[0]/pointCount[idx[j]]; 100 p[j].coords[1]=positions[idx[j]].coords[1]/pointCount[idx[j]]; 101 p[j].coords[2]=positions[idx[j]].coords[2]/pointCount[idx[j]]; 102 } 103 for(j=0;j<3;j++){ 104 q[0].coords[j]=p[1].coords[j]-p[0].coords[j]; 105 q[1].coords[j]=p[2].coords[j]-p[0].coords[j]; 106 d[j]=SquareDistance(p[j],p[(j+1)%3]); 107 } 108 CrossProduct(q[0],q[1],c); 109 a=Length(c)/2; 110 111 if((d[0]+d[1]+d[2])*edgeRatio > a*Ratio){ 112 // Find the smallest edge 113 j=0; 114 if(d[1]<d[j]){j=1;} 115 if(d[2]<d[j]){j=2;} 116 117 int idx1,idx2; 118 if(idx[j]<idx[(j+1)%3]){ 119 idx1=idx[j]; 120 idx2=idx[(j+1)%3]; 121 } 122 else{ 123 idx2=idx[j]; 124 idx1=idx[(j+1)%3]; 125 } 126 positions[idx1].coords[0]+=positions[idx2].coords[0]; 127 positions[idx1].coords[1]+=positions[idx2].coords[1]; 128 positions[idx1].coords[2]+=positions[idx2].coords[2]; 129 if(normals){ 130 (*normals)[idx1].coords[0]+=(*normals)[idx2].coords[0]; 131 (*normals)[idx1].coords[1]+=(*normals)[idx2].coords[1]; 132 (*normals)[idx1].coords[2]+=(*normals)[idx2].coords[2]; 133 } 134 pointCount[idx1]+=pointCount[idx2]; 135 remapTable[idx2]=idx1; 136 triangles[i]=triangles[triangles.size()-1]; 137 triangles.pop_back(); 138 } 139 } 140 int pCount=0; 141 for(i=0;i<int(positions.size());i++){ 142 for(j=0;j<3;j++){positions[i].coords[j]/=pointCount[i];} 143 if(normals){ 144 Real l=Real(Length((*normals)[i])); 145 for(j=0;j<3;j++){(*normals)[i].coords[j]/=l;} 146 } 147 if(remapTable[i]==i){ // If vertex i is being used 148 positions[pCount]=positions[i]; 149 if(normals){(*normals)[pCount]=(*normals)[i];} 150 pointCount[i]=pCount; 151 pCount++; 152 } 153 } 154 positions.resize(pCount); 155 for(i=int(triangles.size()-1);i>=0;i--){ 156 for(j=0;j<3;j++){ 157 idx[j]=triangles[i].idx[j]; 158 while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];} 159 triangles[i].idx[j]=pointCount[idx[j]]; 160 } 161 if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){ 162 triangles[i]=triangles[triangles.size()-1]; 163 triangles.pop_back(); 164 } 165 } 166 167 delete[] pointCount; 168 delete[] remapTable; 169} 170template<class Real> 171void TriangleCollapse(const Real& edgeRatio,std::vector<TriangleIndex>& triangles,std::vector< Point3D<Real> >& positions,std::vector< Point3D<Real> >* normals){ 172 int i,j,*remapTable,*pointCount,idx[3]; 173 Point3D<Real> p[3],q[2],c; 174 double d[3],a; 175 double Ratio=12.0/sqrt(3.0); // (Sum of Squares Length / Area) for and equilateral triangle 176 177 remapTable=new int[positions.size()]; 178 pointCount=new int[positions.size()]; 179 for(i=0;i<int(positions.size());i++){ 180 remapTable[i]=i; 181 pointCount[i]=1; 182 } 183 for(i=int(triangles.size()-1);i>=0;i--){ 184 for(j=0;j<3;j++){ 185 idx[j]=triangles[i].idx[j]; 186 while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];} 187 } 188 if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){ 189 triangles[i]=triangles[triangles.size()-1]; 190 triangles.pop_back(); 191 continue; 192 } 193 for(j=0;j<3;j++){ 194 p[j].coords[0]=positions[idx[j]].coords[0]/pointCount[idx[j]]; 195 p[j].coords[1]=positions[idx[j]].coords[1]/pointCount[idx[j]]; 196 p[j].coords[2]=positions[idx[j]].coords[2]/pointCount[idx[j]]; 197 } 198 for(j=0;j<3;j++){ 199 q[0].coords[j]=p[1].coords[j]-p[0].coords[j]; 200 q[1].coords[j]=p[2].coords[j]-p[0].coords[j]; 201 d[j]=SquareDistance(p[j],p[(j+1)%3]); 202 } 203 CrossProduct(q[0],q[1],c); 204 a=Length(c)/2; 205 206 if((d[0]+d[1]+d[2])*edgeRatio > a*Ratio){ 207 // Find the smallest edge 208 j=0; 209 if(d[1]<d[j]){j=1;} 210 if(d[2]<d[j]){j=2;} 211 212 int idx1,idx2,idx3; 213 if(idx[0]<idx[1]){ 214 if(idx[0]<idx[2]){ 215 idx1=idx[0]; 216 idx2=idx[2]; 217 idx3=idx[1]; 218 } 219 else{ 220 idx1=idx[2]; 221 idx2=idx[0]; 222 idx3=idx[1]; 223 } 224 } 225 else{ 226 if(idx[1]<idx[2]){ 227 idx1=idx[1]; 228 idx2=idx[2]; 229 idx3=idx[0]; 230 } 231 else{ 232 idx1=idx[2]; 233 idx2=idx[1]; 234 idx3=idx[0]; 235 } 236 } 237 positions[idx1].coords[0]+=positions[idx2].coords[0]+positions[idx3].coords[0]; 238 positions[idx1].coords[1]+=positions[idx2].coords[1]+positions[idx3].coords[1]; 239 positions[idx1].coords[2]+=positions[idx2].coords[2]+positions[idx3].coords[2]; 240 if(normals){ 241 (*normals)[idx1].coords[0]+=(*normals)[idx2].coords[0]+(*normals)[idx3].coords[0]; 242 (*normals)[idx1].coords[1]+=(*normals)[idx2].coords[1]+(*normals)[idx3].coords[1]; 243 (*normals)[idx1].coords[2]+=(*normals)[idx2].coords[2]+(*normals)[idx3].coords[2]; 244 } 245 pointCount[idx1]+=pointCount[idx2]+pointCount[idx3]; 246 remapTable[idx2]=idx1; 247 remapTable[idx3]=idx1; 248 triangles[i]=triangles[triangles.size()-1]; 249 triangles.pop_back(); 250 } 251 } 252 int pCount=0; 253 for(i=0;i<int(positions.size());i++){ 254 for(j=0;j<3;j++){positions[i].coords[j]/=pointCount[i];} 255 if(normals){ 256 Real l=Real(Length((*normals)[i])); 257 for(j=0;j<3;j++){(*normals)[i].coords[j]/=l;} 258 } 259 if(remapTable[i]==i){ // If vertex i is being used 260 positions[pCount]=positions[i]; 261 if(normals){(*normals)[pCount]=(*normals)[i];} 262 pointCount[i]=pCount; 263 pCount++; 264 } 265 } 266 positions.resize(pCount); 267 for(i=int(triangles.size()-1);i>=0;i--){ 268 for(j=0;j<3;j++){ 269 idx[j]=triangles[i].idx[j]; 270 while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];} 271 triangles[i].idx[j]=pointCount[idx[j]]; 272 } 273 if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){ 274 triangles[i]=triangles[triangles.size()-1]; 275 triangles.pop_back(); 276 } 277 } 278 delete[] pointCount; 279 delete[] remapTable; 280} 281 282/////////////////// 283// Triangulation // 284/////////////////// 285template<class Real> 286long long Triangulation<Real>::EdgeIndex( int p1 , int p2 ) 287{ 288 if(p1>p2) {return ((long long)(p1)<<32) | ((long long)(p2));} 289 else {return ((long long)(p2)<<32) | ((long long)(p1));} 290} 291 292template<class Real> 293int Triangulation<Real>::factor(int tIndex,int& p1,int& p2,int & p3){ 294 if(triangles[tIndex].eIndex[0]<0 || triangles[tIndex].eIndex[1]<0 || triangles[tIndex].eIndex[2]<0){return 0;} 295 if(edges[triangles[tIndex].eIndex[0]].tIndex[0]==tIndex){p1=edges[triangles[tIndex].eIndex[0]].pIndex[0];} 296 else {p1=edges[triangles[tIndex].eIndex[0]].pIndex[1];} 297 if(edges[triangles[tIndex].eIndex[1]].tIndex[0]==tIndex){p2=edges[triangles[tIndex].eIndex[1]].pIndex[0];} 298 else {p2=edges[triangles[tIndex].eIndex[1]].pIndex[1];} 299 if(edges[triangles[tIndex].eIndex[2]].tIndex[0]==tIndex){p3=edges[triangles[tIndex].eIndex[2]].pIndex[0];} 300 else {p3=edges[triangles[tIndex].eIndex[2]].pIndex[1];} 301 return 1; 302} 303template<class Real> 304double Triangulation<Real>::area(int p1,int p2,int p3){ 305 Point3D<Real> q1,q2,q; 306 for(int i=0;i<3;i++){ 307 q1.coords[i]=points[p2].coords[i]-points[p1].coords[i]; 308 q2.coords[i]=points[p3].coords[i]-points[p1].coords[i]; 309 } 310 CrossProduct(q1,q2,q); 311 return Length(q); 312} 313template<class Real> 314double Triangulation<Real>::area(int tIndex){ 315 int p1,p2,p3; 316 factor(tIndex,p1,p2,p3); 317 return area(p1,p2,p3); 318} 319template<class Real> 320double Triangulation<Real>::area(void){ 321 double a=0; 322 for(int i=0;i<int(triangles.size());i++){a+=area(i);} 323 return a; 324} 325template<class Real> 326int Triangulation<Real>::addTriangle(int p1,int p2,int p3){ 327 hash_map<long long,int>::iterator iter; 328 int tIdx,eIdx,p[3]; 329 p[0]=p1; 330 p[1]=p2; 331 p[2]=p3; 332 triangles.push_back(TriangulationTriangle()); 333 tIdx=int(triangles.size())-1; 334 335 for(int i=0;i<3;i++) 336 { 337 long long e = EdgeIndex(p[i],p[(i+1)%3]); 338 iter=edgeMap.find(e); 339 if(iter==edgeMap.end()) 340 { 341 TriangulationEdge edge; 342 edge.pIndex[0]=p[i]; 343 edge.pIndex[1]=p[(i+1)%3]; 344 edges.push_back(edge); 345 eIdx=int(edges.size())-1; 346 edgeMap[e]=eIdx; 347 edges[eIdx].tIndex[0]=tIdx; 348 } 349 else{ 350 eIdx=edgeMap[e]; 351 if(edges[eIdx].pIndex[0]==p[i]){ 352 if(edges[eIdx].tIndex[0]<0){edges[eIdx].tIndex[0]=tIdx;} 353 else{printf("Edge Triangle in use 1\n");return 0;} 354 } 355 else{ 356 if(edges[eIdx].tIndex[1]<0){edges[eIdx].tIndex[1]=tIdx;} 357 else{printf("Edge Triangle in use 2\n");return 0;} 358 } 359 360 } 361 triangles[tIdx].eIndex[i]=eIdx; 362 } 363 return tIdx; 364} 365template<class Real> 366int Triangulation<Real>::flipMinimize(int eIndex){ 367 double oldArea,newArea; 368 int oldP[3],oldQ[3],newP[3],newQ[3]; 369 TriangulationEdge newEdge; 370 371 if(edges[eIndex].tIndex[0]<0 || edges[eIndex].tIndex[1]<0){return 0;} 372 373 if(!factor(edges[eIndex].tIndex[0],oldP[0],oldP[1],oldP[2])){return 0;} 374 if(!factor(edges[eIndex].tIndex[1],oldQ[0],oldQ[1],oldQ[2])){return 0;} 375 376 oldArea=area(oldP[0],oldP[1],oldP[2])+area(oldQ[0],oldQ[1],oldQ[2]); 377 int idxP,idxQ; 378 for(idxP=0;idxP<3;idxP++){ 379 int i; 380 for(i=0;i<3;i++){if(oldP[idxP]==oldQ[i]){break;}} 381 if(i==3){break;} 382 } 383 for(idxQ=0;idxQ<3;idxQ++){ 384 int i; 385 for(i=0;i<3;i++){if(oldP[i]==oldQ[idxQ]){break;}} 386 if(i==3){break;} 387 } 388 if(idxP==3 || idxQ==3){return 0;} 389 newP[0]=oldP[idxP]; 390 newP[1]=oldP[(idxP+1)%3]; 391 newP[2]=oldQ[idxQ]; 392 newQ[0]=oldQ[idxQ]; 393 newQ[1]=oldP[(idxP+2)%3]; 394 newQ[2]=oldP[idxP]; 395 396 newArea=area(newP[0],newP[1],newP[2])+area(newQ[0],newQ[1],newQ[2]); 397 if(oldArea<=newArea){return 0;} 398 399 // Remove the entry in the hash_table for the old edge 400 edgeMap.erase(EdgeIndex(edges[eIndex].pIndex[0],edges[eIndex].pIndex[1])); 401 // Set the new edge so that the zero-side is newQ 402 edges[eIndex].pIndex[0]=newP[0]; 403 edges[eIndex].pIndex[1]=newQ[0]; 404 // Insert the entry into the hash_table for the new edge 405 edgeMap[EdgeIndex(newP[0],newQ[0])]=eIndex; 406 // Update the triangle information 407 for(int i=0;i<3;i++){ 408 int idx; 409 idx=edgeMap[EdgeIndex(newQ[i],newQ[(i+1)%3])]; 410 triangles[edges[eIndex].tIndex[0]].eIndex[i]=idx; 411 if(idx!=eIndex){ 412 if(edges[idx].tIndex[0]==edges[eIndex].tIndex[1]){edges[idx].tIndex[0]=edges[eIndex].tIndex[0];} 413 if(edges[idx].tIndex[1]==edges[eIndex].tIndex[1]){edges[idx].tIndex[1]=edges[eIndex].tIndex[0];} 414 } 415 416 idx=edgeMap[EdgeIndex(newP[i],newP[(i+1)%3])]; 417 triangles[edges[eIndex].tIndex[1]].eIndex[i]=idx; 418 if(idx!=eIndex){ 419 if(edges[idx].tIndex[0]==edges[eIndex].tIndex[0]){edges[idx].tIndex[0]=edges[eIndex].tIndex[1];} 420 if(edges[idx].tIndex[1]==edges[eIndex].tIndex[0]){edges[idx].tIndex[1]=edges[eIndex].tIndex[1];} 421 } 422 } 423 return 1; 424} 425///////////////////////// 426// CoredVectorMeshData // 427///////////////////////// 428template< class Vertex > 429CoredVectorMeshData< Vertex >::CoredVectorMeshData( void ) { oocPointIndex = polygonIndex = 0; } 430template< class Vertex > 431void CoredVectorMeshData< Vertex >::resetIterator ( void ) { oocPointIndex = polygonIndex = 0; } 432template< class Vertex > 433int CoredVectorMeshData< Vertex >::addOutOfCorePoint( const Vertex& p ) 434{ 435 oocPoints.push_back(p); 436 return int(oocPoints.size())-1; 437} 438template< class Vertex > 439int CoredVectorMeshData< Vertex >::addOutOfCorePoint_s( const Vertex& p ) 440{ 441 size_t sz; 442#pragma omp critical (CoredVectorMeshData_addOutOfCorePoint_s ) 443 { 444 sz = oocPoints.size(); 445 oocPoints.push_back(p); 446 } 447 return (int)sz; 448} 449template< class Vertex > 450int CoredVectorMeshData< Vertex >::addPolygon_s( const std::vector< int >& polygon ) 451{ 452 size_t sz; 453#pragma omp critical (CoredVectorMeshData_addPolygon_s) 454 { 455 sz = polygon.size(); 456 polygons.push_back( polygon ); 457 } 458 return (int)sz; 459} 460template< class Vertex > 461int CoredVectorMeshData< Vertex >::addPolygon_s( const std::vector< CoredVertexIndex >& vertices ) 462{ 463 std::vector< int > polygon( vertices.size() ); 464 for( int i=0 ; i<(int)vertices.size() ; i++ ) 465 if( vertices[i].inCore ) polygon[i] = vertices[i].idx; 466 else polygon[i] = -vertices[i].idx-1; 467 return addPolygon_s( polygon ); 468} 469template< class Vertex > 470int CoredVectorMeshData< Vertex >::nextOutOfCorePoint( Vertex& p ) 471{ 472 if( oocPointIndex<int(oocPoints.size()) ) 473 { 474 p=oocPoints[oocPointIndex++]; 475 return 1; 476 } 477 else{return 0;} 478} 479template< class Vertex > 480int CoredVectorMeshData< Vertex >::nextPolygon( std::vector< CoredVertexIndex >& vertices ) 481{ 482 if( polygonIndex<int( polygons.size() ) ) 483 { 484 std::vector< int >& polygon = polygons[ polygonIndex++ ]; 485 vertices.resize( polygon.size() ); 486 for( int i=0 ; i<int(polygon.size()) ; i++ ) 487 if( polygon[i]<0 ) vertices[i].idx = -polygon[i]-1 , vertices[i].inCore = false; 488 else vertices[i].idx = polygon[i] , vertices[i].inCore = true; 489 return 1; 490 } 491 else return 0; 492} 493template< class Vertex > 494int CoredVectorMeshData< Vertex >::outOfCorePointCount(void){return int(oocPoints.size());} 495template< class Vertex > 496int CoredVectorMeshData< Vertex >::polygonCount( void ) { return int( polygons.size() ); } 497 498/////////////////////// 499// CoredFileMeshData // 500/////////////////////// 501template< class Vertex > 502CoredFileMeshData< Vertex >::CoredFileMeshData( void ) 503{ 504 oocPoints = polygons = 0; 505 506 oocPointFile = new BufferedReadWriteFile(); 507 polygonFile = new BufferedReadWriteFile(); 508} 509template< class Vertex > 510CoredFileMeshData< Vertex >::~CoredFileMeshData( void ) 511{ 512 delete oocPointFile; 513 delete polygonFile; 514} 515template< class Vertex > 516void CoredFileMeshData< Vertex >::resetIterator ( void ) 517{ 518 oocPointFile->reset(); 519 polygonFile->reset(); 520} 521template< class Vertex > 522int CoredFileMeshData< Vertex >::addOutOfCorePoint( const Vertex& p ) 523{ 524 oocPointFile->write( &p , sizeof( Vertex ) ); 525 oocPoints++; 526 return oocPoints-1; 527} 528template< class Vertex > 529int CoredFileMeshData< Vertex >::addOutOfCorePoint_s( const Vertex& p ) 530{ 531 int sz; 532#pragma omp critical (CoredFileMeshData_addOutOfCorePoint_s) 533 { 534 sz = oocPoints; 535 oocPointFile->write( &p , sizeof( Vertex ) ); 536 oocPoints++; 537 } 538 return sz; 539} 540template< class Vertex > 541int CoredFileMeshData< Vertex >::addPolygon_s( const std::vector< int >& vertices ) 542{ 543 int sz , vSize = (int)vertices.size(); 544#pragma omp critical (CoredFileMeshData_addPolygon_s ) 545 { 546 sz = polygons; 547 polygonFile->write( &vSize , sizeof(int) ); 548 polygonFile->write( &vertices[0] , sizeof(int) * vSize ); 549 polygons++; 550 } 551 return sz; 552} 553template< class Vertex > 554int CoredFileMeshData< Vertex >::addPolygon_s( const std::vector< CoredVertexIndex >& vertices ) 555{ 556 std::vector< int > polygon( vertices.size() ); 557 for( int i=0 ; i<(int)vertices.size() ; i++ ) 558 if( vertices[i].inCore ) polygon[i] = vertices[i].idx; 559 else polygon[i] = -vertices[i].idx-1; 560 return addPolygon_s( polygon ); 561} 562template< class Vertex > 563int CoredFileMeshData< Vertex >::nextOutOfCorePoint( Vertex& p ) 564{ 565 if( oocPointFile->read( &p , sizeof( Vertex ) ) ) return 1; 566 else return 0; 567} 568template< class Vertex > 569int CoredFileMeshData< Vertex >::nextPolygon( std::vector< CoredVertexIndex >& vertices ) 570{ 571 int pSize; 572 if( polygonFile->read( &pSize , sizeof(int) ) ) 573 { 574 std::vector< int > polygon( pSize ); 575 if( polygonFile->read( &polygon[0] , sizeof(int)*pSize ) ) 576 { 577 vertices.resize( pSize ); 578 for( int i=0 ; i<int(polygon.size()) ; i++ ) 579 if( polygon[i]<0 ) vertices[i].idx = -polygon[i]-1 , vertices[i].inCore = false; 580 else vertices[i].idx = polygon[i] , vertices[i].inCore = true; 581 return 1; 582 } 583 return 0; 584 } 585 else return 0; 586} 587template< class Vertex > 588int CoredFileMeshData< Vertex >::outOfCorePointCount( void ){ return oocPoints; } 589template< class Vertex > 590int CoredFileMeshData< Vertex >::polygonCount( void ) { return polygons; } 591