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#include <stdlib.h>
31#include <math.h>
32#include <float.h>
33#ifdef _WIN32
34#include <Windows.h>
35#include <Psapi.h>
36#endif // _WIN32
37#include "MyTime.h"
38#include "MarchingCubes.h"
39#include "Octree.h"
40#include "SparseMatrix.h"
41#include "CmdLineParser.h"
42#include "PPolynomial.h"
43#include "Ply.h"
44#include "MemoryUsage.h"
45#ifdef _OPENMP
46#include "omp.h"
47#endif // _OPENMP
48void DumpOutput( const char* format , ... );
49void DumpOutput2( std::vector< char* >& comments , const char* format , ... );
50#include "MultiGridOctreeData.h"
51
52#define DEFAULT_FULL_DEPTH 5
53
54#define XSTR(x) STR(x)
55#define STR(x) #x
56#if DEFAULT_FULL_DEPTH
57#pragma message ( "[WARNING] Setting default full depth to " XSTR(DEFAULT_FULL_DEPTH) )
58#endif // DEFAULT_FULL_DEPTH
59
60#include <stdarg.h>
61char* outputFile=NULL;
62int echoStdout=0;
63void DumpOutput( const char* format , ... )
64{
65	if( outputFile )
66	{
67		FILE* fp = fopen( outputFile , "a" );
68		va_list args;
69		va_start( args , format );
70		vfprintf( fp , format , args );
71		fclose( fp );
72		va_end( args );
73	}
74	if( echoStdout )
75	{
76		va_list args;
77		va_start( args , format );
78		vprintf( format , args );
79		va_end( args );
80	}
81}
82void DumpOutput2( std::vector< char* >& comments  , const char* format , ... )
83{
84	if( outputFile )
85	{
86		FILE* fp = fopen( outputFile , "a" );
87		va_list args;
88		va_start( args , format );
89		vfprintf( fp , format , args );
90		fclose( fp );
91		va_end( args );
92	}
93	if( echoStdout )
94	{
95		va_list args;
96		va_start( args , format );
97		vprintf( format , args );
98		va_end( args );
99	}
100	comments.push_back( new char[1024] );
101	char* str = comments.back();
102	va_list args;
103	va_start( args , format );
104	vsprintf( str , format , args );
105	va_end( args );
106	if( str[strlen(str)-1]=='\n' ) str[strlen(str)-1] = 0;
107}
108
109
110cmdLineString
111	In( "in" ) ,
112	Out( "out" ) ,
113	VoxelGrid( "voxel" ) ,
114	XForm( "xForm" );
115
116cmdLineReadable
117#ifdef _WIN32
118	Performance( "performance" ) ,
119#endif // _WIN32
120	Complete( "complete" ) ,
121	ShowResidual( "showResidual" ) ,
122	NoComments( "noComments" ) ,
123	PolygonMesh( "polygonMesh" ) ,
124	Confidence( "confidence" ) ,
125	NormalWeights( "nWeights" ) ,
126	NonManifold( "nonManifold" ) ,
127	Dirichlet( "dirichlet" ) ,
128	ASCII( "ascii" ) ,
129	Density( "density" ) ,
130	LinearFit( "linearFit" ) ,
131	PrimalVoxel( "primalVoxel" ) ,
132	Verbose( "verbose" ) ,
133	Double( "double" );
134
135cmdLineInt
136	Degree( "degree" , 2 ) ,
137	Depth( "depth" , 8 ) ,
138	CGDepth( "cgDepth" , 0 ) ,
139	KernelDepth( "kernelDepth" ) ,
140	AdaptiveExponent( "adaptiveExp" , 1 ) ,
141	Iters( "iters" , 8 ) ,
142	VoxelDepth( "voxelDepth" , -1 ) ,
143	FullDepth( "fullDepth" , DEFAULT_FULL_DEPTH ) ,
144	MinDepth( "minDepth" , 0 ) ,
145	MaxSolveDepth( "maxSolveDepth" ) ,
146	Threads( "threads" , omp_get_num_procs() );
147
148cmdLineFloat
149	Color( "color" , 16.f ) ,
150	SamplesPerNode( "samplesPerNode" , 1.5f ) ,
151	Scale( "scale" , 1.1f ) ,
152	CSSolverAccuracy( "cgAccuracy" , float(1e-3) ) ,
153	PointWeight( "pointWeight" , 4.f );
154
155
156cmdLineReadable* params[] =
157{
158	&In , &Degree , &Depth , &Out , &XForm ,
159	&Scale , &Verbose , &CSSolverAccuracy , &NoComments , &Double ,
160	&KernelDepth , &SamplesPerNode , &Confidence , &NormalWeights , &NonManifold , &PolygonMesh , &ASCII , &ShowResidual , &VoxelDepth ,
161	&PointWeight , &VoxelGrid , &Threads , &MaxSolveDepth ,
162	&AdaptiveExponent , &Dirichlet ,
163	&Density ,
164	&FullDepth ,
165	&MinDepth ,
166	&CGDepth , &Iters ,
167	&Complete ,
168	&Color ,
169	&LinearFit ,
170	&PrimalVoxel ,
171#ifdef _WIN32
172	&Performance ,
173#endif // _WIN32
174};
175
176
177void ShowUsage(char* ex)
178{
179	printf( "Usage: %s\n" , ex );
180	printf( "\t --%s  <input points>\n" , In.name );
181
182	printf( "\t[--%s <ouput triangle mesh>]\n" , Out.name );
183	printf( "\t[--%s <ouput voxel grid>]\n" , VoxelGrid.name );
184
185	printf( "\t[--%s <b-spline degree>=%d]\n" , Degree.name , Degree.value );
186
187	printf( "\t[--%s <maximum reconstruction depth>=%d]\n" , Depth.name , Depth.value );
188	printf( "\t\t Running at depth d corresponds to solving on a 2^d x 2^d x 2^d\n" );
189	printf( "\t\t voxel grid.\n" );
190
191	printf( "\t[--%s <full depth>=%d]\n" , FullDepth.name , FullDepth.value );
192	printf( "\t\t This flag specifies the depth up to which the octree should be complete.\n" );
193
194	printf( "\t[--%s <depth at which to extract the voxel grid>=<%s>]\n" , VoxelDepth.name , Depth.name );
195
196	printf( "\t[--%s <conjugate-gradients depth>=%d]\n" , CGDepth.name , CGDepth.value );
197	printf( "\t\t The depth up to which a conjugate-gradients solver should be used.\n");
198
199	printf( "\t[--%s <scale factor>=%f]\n" , Scale.name , Scale.value );
200	printf( "\t\t Specifies the factor of the bounding cube that the input\n" );
201	printf( "\t\t samples should fit into.\n" );
202
203	printf( "\t[--%s <minimum number of samples per node>=%f]\n" , SamplesPerNode.name, SamplesPerNode.value );
204	printf( "\t\t This parameter specifies the minimum number of points that\n" );
205	printf( "\t\t should fall within an octree node.\n" );
206
207	printf( "\t[--%s <interpolation weight>=%f]\n" , PointWeight.name , PointWeight.value );
208	printf( "\t\t This value specifies the weight that point interpolation constraints are\n" );
209	printf( "\t\t given when defining the (screened) Poisson system.\n" );
210
211	printf( "\t[--%s <iterations>=%d]\n" , Iters.name , Iters.value );
212	printf( "\t\t This flag specifies the (maximum if CG) number of solver iterations.\n" );
213
214	printf( "\t[--%s <pull factor>]\n" , Color.name );
215	printf( "\t\t This flag specifies the pull factor for color interpolation\n" );
216
217#ifdef _OPENMP
218	printf( "\t[--%s <num threads>=%d]\n" , Threads.name , Threads.value );
219	printf( "\t\t This parameter specifies the number of threads across which\n" );
220	printf( "\t\t the solver should be parallelized.\n" );
221#endif // _OPENMP
222
223	printf( "\t[--%s]\n" , Confidence.name );
224	printf( "\t\t If this flag is enabled, the size of a sample's normals is\n" );
225	printf( "\t\t used as a confidence value, affecting the sample's\n" );
226	printf( "\t\t constribution to the reconstruction process.\n" );
227
228	printf( "\t[--%s]\n" , NormalWeights.name );
229	printf( "\t\t If this flag is enabled, the size of a sample's normals is\n" );
230	printf( "\t\t used as to modulate the interpolation weight.\n" );
231
232#if 0
233	printf( "\t[--%s]\n" , NonManifold.name );
234	printf( "\t\t If this flag is enabled, the isosurface extraction does not add\n" );
235	printf( "\t\t a planar polygon's barycenter in order to ensure that the output\n" );
236	printf( "\t\t mesh is manifold.\n" );
237#endif
238
239	printf( "\t[--%s]\n" , PolygonMesh.name);
240	printf( "\t\t If this flag is enabled, the isosurface extraction returns polygons\n" );
241	printf( "\t\t rather than triangles.\n" );
242
243#if 0
244	printf( "\t[--%s <minimum depth>=%d]\n" , MinDepth.name , MinDepth.value );
245	printf( "\t\t This flag specifies the coarsest depth at which the system is to be solved.\n" );
246
247	printf( "\t[--%s <cg solver accuracy>=%g]\n" , CSSolverAccuracy.name , CSSolverAccuracy.value );
248	printf( "\t\t This flag specifies the accuracy cut-off to be used for CG.\n" );
249
250	printf( "\t[--%s <adaptive weighting exponent>=%d]\n", AdaptiveExponent.name , AdaptiveExponent.value );
251	printf( "\t\t This flag specifies the exponent scale for the adaptive weighting.\n" );
252
253#ifdef _WIN32
254	printf( "\t[--%s]\n" , Performance.name );
255	printf( "\t\t If this flag is enabled, the running time and peak memory usage\n" );
256	printf( "\t\t is output after the reconstruction.\n" );
257#endif // _WIN32
258#endif
259	printf( "\t[--%s]\n" , Dirichlet.name);
260	printf( "\t\t If this flag is enabled, Dirichlet boundary constraints are used for reconstruction.\n" );
261
262	printf( "\t[--%s]\n" , Density.name );
263	printf( "\t\t If this flag is enabled, the sampling density is written out with the vertices.\n" );
264
265	printf( "\t[--%s]\n" , LinearFit.name );
266	printf( "\t\t If this flag is enabled, the iso-surfacing will be performed using linear fitting.\n" );
267
268	printf( "\t[--%s]\n" , PrimalVoxel.name );
269	printf( "\t\t If this flag is enabled, voxel sampling is performed at corners rather than centers.\n" );
270
271#if 0
272	printf( "\t[--%s]\n" , ASCII.name );
273	printf( "\t\t If this flag is enabled, the output file is written out in ASCII format.\n" );
274
275	printf( "\t[--%s]\n" , NoComments.name );
276	printf( "\t\t If this flag is enabled, the output file will not include comments.\n" );
277#endif
278
279	printf( "\t[--%s]\n" , Double.name );
280	printf( "\t\t If this flag is enabled, the reconstruction will be performed with double-precision floats.\n" );
281
282	printf( "\t[--%s]\n" , Verbose.name );
283	printf( "\t\t If this flag is enabled, the progress of the reconstructor will be output to STDOUT.\n" );
284}
285
286Point3D< unsigned char > ReadASCIIColor( FILE* fp )
287{
288	Point3D< unsigned char > c;
289	if( fscanf( fp , " %c %c %c " , &c[0] , &c[1] , &c[2] )!=3 ) fprintf( stderr , "[ERROR] Failed to read color\n" ) , exit( 0 );
290	return c;
291}
292
293PlyProperty PlyColorProperties[]=
294{
295	{ "r"     , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[0] ) ) , 0 , 0 , 0 , 0 } ,
296	{ "g"     , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[1] ) ) , 0 , 0 , 0 , 0 } ,
297	{ "b"     , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[2] ) ) , 0 , 0 , 0 , 0 } ,
298	{ "red"   , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[0] ) ) , 0 , 0 , 0 , 0 } ,
299	{ "green" , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[1] ) ) , 0 , 0 , 0 , 0 } ,
300	{ "blue"  , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[2] ) ) , 0 , 0 , 0 , 0 }
301};
302
303bool ValidPlyColorProperties( const bool* props ){ return ( props[0] || props[3] ) && ( props[1] || props[4] ) && ( props[2] || props[5] ); }
304
305template< class Real , int Degree , class Vertex >
306int _Execute( int argc , char* argv[] )
307{
308	Reset< Real >();
309	int paramNum = sizeof(params)/sizeof(cmdLineReadable*);
310	std::vector< char* > comments;
311
312	if( Verbose.set ) echoStdout=1;
313
314	XForm4x4< Real > xForm , iXForm;
315	if( XForm.set )
316	{
317		FILE* fp = fopen( XForm.value , "r" );
318		if( !fp )
319		{
320			fprintf( stderr , "[WARNING] Could not read x-form from: %s\n" , XForm.value );
321			xForm = XForm4x4< Real >::Identity();
322		}
323		else
324		{
325			for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ )
326			{
327				float f;
328				if( fscanf( fp , " %f " , &f )!=1 ) fprintf( stderr , "[ERROR] Execute: Failed to read xform\n" ) , exit( 0 );
329				xForm(i,j) = (Real)f;
330			}
331			fclose( fp );
332		}
333	}
334	else xForm = XForm4x4< Real >::Identity();
335	iXForm = xForm.inverse();
336
337	DumpOutput2( comments , "Running Screened Poisson Reconstruction (Version 8.0)\n" );
338	char str[1024];
339	for( int i=0 ; i<paramNum ; i++ )
340		if( params[i]->set )
341		{
342			params[i]->writeValue( str );
343			if( strlen( str ) ) DumpOutput2( comments , "\t--%s %s\n" , params[i]->name , str );
344			else                DumpOutput2( comments , "\t--%s\n" , params[i]->name );
345		}
346
347	double t;
348	double tt=Time();
349	Real isoValue = 0;
350
351	Octree< Real > tree;
352	tree.threads = Threads.value;
353	if( !In.set )
354	{
355		ShowUsage( argv[0] );
356		return 0;
357	}
358	if( !MaxSolveDepth.set ) MaxSolveDepth.value = Depth.value;
359
360	OctNode< TreeNodeData >::SetAllocator( MEMORY_ALLOCATOR_BLOCK_SIZE );
361
362	t=Time();
363	int kernelDepth = KernelDepth.set ? KernelDepth.value : Depth.value-2;
364	if( kernelDepth>Depth.value )
365	{
366		fprintf( stderr,"[WARNING] %s can't be greater than %s: %d <= %d\n" , KernelDepth.name , Depth.name , KernelDepth.value , Depth.value );
367		kernelDepth = Depth.value;
368	}
369
370	double maxMemoryUsage;
371	t=Time() , tree.maxMemoryUsage=0;
372	SparseNodeData< PointData< Real > , 0 >* pointInfo = new SparseNodeData< PointData< Real > , 0 >();
373	SparseNodeData< Point3D< Real > , NORMAL_DEGREE >* normalInfo = new SparseNodeData< Point3D< Real > , NORMAL_DEGREE >();
374	SparseNodeData< Real , WEIGHT_DEGREE >* densityWeights = new SparseNodeData< Real , WEIGHT_DEGREE >();
375	SparseNodeData< Real , NORMAL_DEGREE >* nodeWeights = new SparseNodeData< Real , NORMAL_DEGREE >();
376	int pointCount;
377	typedef typename Octree< Real >::template ProjectiveData< Point3D< Real > > ProjectiveColor;
378	SparseNodeData< ProjectiveColor , DATA_DEGREE >* colorData = NULL;
379
380	char* ext = GetFileExtension( In.value );
381	if( Color.set && Color.value>0 )
382	{
383		colorData = new SparseNodeData< ProjectiveColor , DATA_DEGREE >();
384		OrientedPointStreamWithData< float , Point3D< unsigned char > >* pointStream;
385		if     ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStreamWithData< float , Point3D< unsigned char > >( In.value );
386		else if( !strcasecmp( ext , "ply"   ) ) pointStream = new    PLYOrientedPointStreamWithData< float , Point3D< unsigned char > >( In.value , PlyColorProperties , 6 , ValidPlyColorProperties );
387		else                                    pointStream = new  ASCIIOrientedPointStreamWithData< float , Point3D< unsigned char > >( In.value , ReadASCIIColor );
388		pointCount = tree.template SetTree< float , NORMAL_DEGREE , WEIGHT_DEGREE , DATA_DEGREE , Point3D< unsigned char > >( pointStream , MinDepth.value , Depth.value , FullDepth.value , kernelDepth , Real(SamplesPerNode.value) , Scale.value , Confidence.set , NormalWeights.set , PointWeight.value , AdaptiveExponent.value , *densityWeights , *pointInfo , *normalInfo , *nodeWeights , colorData , xForm , Dirichlet.set , Complete.set );
389		delete pointStream;
390
391		for( const OctNode< TreeNodeData >* n = tree.tree().nextNode() ; n!=NULL ; n=tree.tree().nextNode( n ) )
392		{
393			int idx = colorData->index( n );
394			if( idx>=0 ) colorData->data[idx] *= (Real)pow( Color.value , n->depth() );
395		}
396	}
397	else
398	{
399		OrientedPointStream< float >* pointStream;
400		if     ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStream< float >( In.value );
401		else if( !strcasecmp( ext , "ply"   ) ) pointStream = new    PLYOrientedPointStream< float >( In.value );
402		else                                    pointStream = new  ASCIIOrientedPointStream< float >( In.value );
403		pointCount = tree.template SetTree< float , NORMAL_DEGREE , WEIGHT_DEGREE , DATA_DEGREE , Point3D< unsigned char > >( pointStream , MinDepth.value , Depth.value , FullDepth.value , kernelDepth , Real(SamplesPerNode.value) , Scale.value , Confidence.set , NormalWeights.set , PointWeight.value , AdaptiveExponent.value , *densityWeights , *pointInfo , *normalInfo , *nodeWeights , colorData , xForm , Dirichlet.set , Complete.set );
404		delete pointStream;
405	}
406	delete[] ext;
407	if( !Density.set ) delete densityWeights , densityWeights = NULL;
408	{
409		std::vector< int > indexMap;
410		if( NORMAL_DEGREE>Degree ) tree.template EnableMultigrid< NORMAL_DEGREE >( &indexMap );
411		else                       tree.template EnableMultigrid<        Degree >( &indexMap );
412		if( pointInfo ) pointInfo->remapIndices( indexMap );
413		if( normalInfo ) normalInfo->remapIndices( indexMap );
414		if( densityWeights ) densityWeights->remapIndices( indexMap );
415		if( nodeWeights ) nodeWeights->remapIndices( indexMap );
416		if( colorData ) colorData->remapIndices( indexMap );
417	}
418
419	DumpOutput2( comments , "#             Tree set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
420	DumpOutput( "Input Points: %d\n" , pointCount );
421	DumpOutput( "Leaves/Nodes: %d/%d\n" , tree.leaves() , tree.nodes() );
422	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) );
423
424	maxMemoryUsage = tree.maxMemoryUsage;
425	t=Time() , tree.maxMemoryUsage=0;
426	DenseNodeData< Real , Degree > constraints = tree.template SetLaplacianConstraints< Degree >( *normalInfo );
427	delete normalInfo;
428	DumpOutput2( comments , "#      Constraints set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
429	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) );
430	maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage );
431
432	t=Time() , tree.maxMemoryUsage=0;
433	DenseNodeData< Real , Degree > solution = tree.SolveSystem( *pointInfo , constraints , ShowResidual.set , Iters.value , MaxSolveDepth.value , CGDepth.value , CSSolverAccuracy.value );
434	delete pointInfo;
435	constraints.resize( 0 );
436
437	DumpOutput2( comments , "# Linear system solved in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
438	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) );
439	maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage );
440
441	CoredFileMeshData< Vertex > mesh;
442
443	if( Verbose.set ) tree.maxMemoryUsage=0;
444	t=Time();
445	isoValue = tree.GetIsoValue( solution , *nodeWeights );
446	delete nodeWeights;
447	DumpOutput( "Got average in: %f\n" , Time()-t );
448	DumpOutput( "Iso-Value: %e\n" , isoValue );
449
450	if( VoxelGrid.set )
451	{
452		double t = Time();
453		FILE* fp = fopen( VoxelGrid.value , "wb" );
454		if( !fp ) fprintf( stderr , "Failed to open voxel file for writing: %s\n" , VoxelGrid.value );
455		else
456		{
457			int res = 0;
458			Pointer( Real ) values = tree.Evaluate( solution , res , isoValue , VoxelDepth.value , PrimalVoxel.set );
459			fwrite( &res , sizeof(int) , 1 , fp );
460			if( sizeof(Real)==sizeof(float) ) fwrite( values , sizeof(float) , res*res*res , fp );
461			else
462			{
463				float *fValues = new float[res*res*res];
464				for( int i=0 ; i<res*res*res ; i++ ) fValues[i] = float( values[i] );
465				fwrite( fValues , sizeof(float) , res*res*res , fp );
466				delete[] fValues;
467			}
468			fclose( fp );
469			DeletePointer( values );
470		}
471		DumpOutput( "Got voxel grid in: %f\n" , Time()-t );
472	}
473
474	if( Out.set )
475	{
476		t = Time() , tree.maxMemoryUsage = 0;
477		tree.template GetMCIsoSurface< Degree , WEIGHT_DEGREE , DATA_DEGREE >( densityWeights , colorData , solution , isoValue , mesh , !LinearFit.set , !NonManifold.set , PolygonMesh.set );
478		if( PolygonMesh.set ) DumpOutput2( comments , "#         Got polygons in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
479		else                  DumpOutput2( comments , "#        Got triangles in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
480		maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage );
481		DumpOutput2( comments , "#             Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-tt , maxMemoryUsage );
482
483		if( NoComments.set )
484		{
485			if( ASCII.set ) PlyWritePolygons( Out.value , &mesh , PLY_ASCII         , NULL , 0 , iXForm );
486			else            PlyWritePolygons( Out.value , &mesh , PLY_BINARY_NATIVE , NULL , 0 , iXForm );
487		}
488		else
489		{
490			if( ASCII.set ) PlyWritePolygons( Out.value , &mesh , PLY_ASCII         , &comments[0] , (int)comments.size() , iXForm );
491			else            PlyWritePolygons( Out.value , &mesh , PLY_BINARY_NATIVE , &comments[0] , (int)comments.size() , iXForm );
492		}
493		DumpOutput( "Vertices / Polygons: %d / %d\n" , mesh.outOfCorePointCount()+mesh.inCorePoints.size() , mesh.polygonCount() );
494	}
495	solution.resize( 0 );
496	if( colorData ){ delete colorData ; colorData = NULL; }
497	return 1;
498}
499
500#ifdef _WIN32
501inline double to_seconds( const FILETIME& ft )
502{
503	const double low_to_sec=100e-9; // 100 nanoseconds
504	const double high_to_sec=low_to_sec*4294967296.0;
505	return ft.dwLowDateTime*low_to_sec+ft.dwHighDateTime*high_to_sec;
506}
507#endif // _WIN32
508
509template< class Real , class Vertex >
510int Execute( int argc , char* argv[] )
511{
512	switch( Degree.value )
513	{
514	case 1: return _Execute< Real , 1 , Vertex >( argc , argv );
515	case 2: return _Execute< Real , 2 , Vertex >( argc , argv );
516	case 3: return _Execute< Real , 3 , Vertex >( argc , argv );
517	case 4: return _Execute< Real , 4 , Vertex >( argc , argv );
518	default:
519		fprintf( stderr , "[ERROR] Only B-Splines of degree 1 - 4 are supported" );
520		return EXIT_FAILURE;
521	}
522}
523int main( int argc , char* argv[] )
524{
525#if defined(WIN32) && defined(MAX_MEMORY_GB)
526	if( MAX_MEMORY_GB>0 )
527	{
528		SIZE_T peakMemory = 1;
529		peakMemory <<= 30;
530		peakMemory *= MAX_MEMORY_GB;
531		printf( "Limiting memory usage to %.2f GB\n" , float( peakMemory>>30 ) );
532		HANDLE h = CreateJobObject( NULL , NULL );
533		AssignProcessToJobObject( h , GetCurrentProcess() );
534
535		JOBOBJECT_EXTENDED_LIMIT_INFORMATION jeli = { 0 };
536		jeli.BasicLimitInformation.LimitFlags = JOB_OBJECT_LIMIT_JOB_MEMORY;
537		jeli.JobMemoryLimit = peakMemory;
538		if( !SetInformationJobObject( h , JobObjectExtendedLimitInformation , &jeli , sizeof( jeli ) ) )
539			fprintf( stderr , "Failed to set memory limit\n" );
540	}
541#endif // defined(WIN32) && defined(MAX_MEMORY_GB)
542	double t = Time();
543
544	cmdLineParse( argc-1 , &argv[1] , sizeof(params)/sizeof(cmdLineReadable*) , params , 1 );
545	if( Density.set )
546		if( Color.set )
547			if( Double.set ) Execute< double , PlyColorAndValueVertex< float > >( argc , argv );
548			else             Execute< float  , PlyColorAndValueVertex< float > >( argc , argv );
549		else
550			if( Double.set ) Execute< double , PlyValueVertex< float > >( argc , argv );
551			else             Execute< float  , PlyValueVertex< float > >( argc , argv );
552	else
553		if( Color.set )
554			if( Double.set ) Execute< double , PlyColorVertex< float > >( argc , argv );
555			else             Execute< float  , PlyColorVertex< float > >( argc , argv );
556		else
557			if( Double.set ) Execute< double , PlyVertex< float > >( argc , argv );
558			else             Execute< float  , PlyVertex< float > >( argc , argv );
559#ifdef _WIN32
560	if( Performance.set )
561	{
562		HANDLE cur_thread=GetCurrentThread();
563		FILETIME tcreat, texit, tkernel, tuser;
564		if( GetThreadTimes( cur_thread , &tcreat , &texit , &tkernel , &tuser ) )
565			printf( "Time (Wall/User/Kernel): %.2f / %.2f / %.2f\n" , Time()-t , to_seconds( tuser ) , to_seconds( tkernel ) );
566		else printf( "Time: %.2f\n" , Time()-t );
567		HANDLE h = GetCurrentProcess();
568		PROCESS_MEMORY_COUNTERS pmc;
569		if( GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ) printf( "Peak Memory (MB): %d\n" , pmc.PeakWorkingSetSize>>20 );
570	}
571#endif // _WIN32
572	return EXIT_SUCCESS;
573}
574