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