From 85e597e11b37f647114e7066b5180fafcdc24d11 Mon Sep 17 00:00:00 2001
From: mkazhdan <misha@cs.jhu.edu>
Date: Tue, 21 Jul 2015 10:07:22 -0400
Subject: [PATCH] Version 7.0 (Continued clean-up)

Extended SurfaceTrimmer to support colors
Consolidated color-wrapper
---
 Src/InterpolateColors.V1.cpp           | 540 -------------------------
 Src/InterpolateColors.cpp              | 407 -------------------
 Src/MultiGridOctreeData.IsoSurface.inl |  14 +-
 Src/MultiGridOctreeData.h              |   7 +-
 Src/MultiGridOctreeData.inl            |   5 +-
 Src/Ply.h                              | 300 ++++++++++----
 Src/PointStream.inl                    |  12 +-
 Src/PoissonRecon.cpp                   |  21 +-
 Src/SurfaceTrimmer.cpp                 | 291 ++++++-------
 9 files changed, 383 insertions(+), 1214 deletions(-)
 delete mode 100644 Src/InterpolateColors.V1.cpp
 delete mode 100644 Src/InterpolateColors.cpp

diff --git a/Src/InterpolateColors.V1.cpp b/Src/InterpolateColors.V1.cpp
deleted file mode 100644
index 86084fa..0000000
--- a/Src/InterpolateColors.V1.cpp
+++ /dev/null
@@ -1,540 +0,0 @@
-
-/*
-Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without modification,
-are permitted provided that the following conditions are met:
-
-Redistributions of source code must retain the above copyright notice, this list of
-conditions and the following disclaimer. Redistributions in binary form must reproduce
-the above copyright notice, this list of conditions and the following disclaimer
-in the documentation and/or other materials provided with the distribution. 
-
-Neither the name of the Johns Hopkins University nor the names of its contributors
-may be used to endorse or promote products derived from this software without specific
-prior written permission. 
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
-EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES 
-OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
-SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
-INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
-TO, PROCUREMENT OF SUBSTITUTE  GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
-BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
-ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
-DAMAGE.
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <float.h>
-#ifdef _WIN32
-#include <Windows.h>
-#include <Psapi.h>
-#endif // _WIN32
-#include "MyTime.h"
-#include "MarchingCubes.h"
-#include "Octree.h"
-#include "SparseMatrix.h"
-#include "CmdLineParser.h"
-#include "PPolynomial.h"
-#include "Ply.h"
-#include "MemoryUsage.h"
-#ifdef _OPENMP
-#include "omp.h"
-#endif // _OPENMP
-void DumpOutput( const char* format , ... );
-#include "MultiGridOctreeData.h"
-void DumpOutput2( std::vector< char* >& comments , const char* format , ... );
-
-#define DEFAULT_FULL_DEPTH 5
-
-#define XSTR(x) STR(x)
-#define STR(x) #x
-#if DEFAULT_FULL_DEPTH
-#pragma message ( "[WARNING] Setting default full depth to " XSTR(DEFAULT_FULL_DEPTH) )
-#endif // DEFAULT_FULL_DEPTH
-
-#include <stdarg.h>
-char* outputFile=NULL;
-int echoStdout=0;
-void DumpOutput( const char* format , ... )
-{
-	if( outputFile )
-	{
-		FILE* fp = fopen( outputFile , "a" );
-		va_list args;
-		va_start( args , format );
-		vfprintf( fp , format , args );
-		fclose( fp );
-		va_end( args );
-	}
-	if( echoStdout )
-	{
-		va_list args;
-		va_start( args , format );
-		vprintf( format , args );
-		va_end( args );
-	}
-}
-void DumpOutput2( std::vector< char* >& comments , const char* format , ... )
-{
-	if( outputFile )
-	{
-		FILE* fp = fopen( outputFile , "a" );
-		va_list args;
-		va_start( args , format );
-		vfprintf( fp , format , args );
-		fclose( fp );
-		va_end( args );
-	}
-	if( echoStdout )
-	{
-		va_list args;
-		va_start( args , format );
-		vprintf( format , args );
-		va_end( args );
-	}
-	comments.push_back( new char[1024] );
-	char* str = comments.back();
-	va_list args;
-	va_start( args , format );
-	vsprintf( str , format , args );
-	va_end( args );
-	if( str[strlen(str)-1]=='\n' ) str[strlen(str)-1] = 0;
-}
-
-
-cmdLineString
-	InPoints( "inPoints" ) ,
-	InMesh( "inMesh" ) ,
-	Out( "out" ) ,
-	VoxelGrid( "voxel" ) ,
-	XForm( "xForm" );
-
-cmdLineReadable
-#ifdef _WIN32
-	Performance( "performance" ) ,
-#endif // _WIN32
-	NoComments( "noComments" ) ,
-	Confidence( "confidence" ) ,
-	NormalWeights( "nWeights" ) ,
-	Verbose( "verbose" ) ,
-	Double( "double" );
-
-cmdLineInt
-	Depth( "depth" , 8 ) ,
-	KernelDepth( "kernelDepth" ) ,
-	VoxelDepth( "voxelDepth" , -1 ) ,
-	FullDepth( "fullDepth" , DEFAULT_FULL_DEPTH ) ,
-	MinDepth( "minDepth" , 0 ) ,
-	BoundaryType( "boundary" , 1 ) ,
-	Threads( "threads" , omp_get_num_procs() );
-
-cmdLineFloat
-	SamplesPerNode( "samplesPerNode" , 1.f ) ,
-	Pull( "pull" , 1.f ) ,
-	Scale( "scale" , 1.1f );
-
-
-cmdLineReadable* params[] =
-{
-	&InPoints , &InMesh , &Depth , &Out , &XForm ,
-	&Pull , &Scale , &Verbose , &NoComments , &Double ,
-	&KernelDepth , &SamplesPerNode , &Confidence , &NormalWeights , &VoxelDepth ,
-	&VoxelGrid , &Threads ,
-	&BoundaryType ,
-	&FullDepth ,
-	&MinDepth ,
-#ifdef _WIN32
-	&Performance ,
-#endif // _WIN32
-};
-
-
-void ShowUsage(char* ex)
-{
-	printf( "Usage: %s\n" , ex );
-	printf( "\t --%s  <input points>\n" , InPoints.name );
-	printf( "\t --%s  <input mesh>\n" , InMesh.name );
-
-	printf( "\t[--%s <ouput triangle mesh>]\n" , Out.name );
-	printf( "\t[--%s <ouput voxel grid>]\n" , VoxelGrid.name );
-
-	printf( "\t[--%s <maximum reconstruction depth>=%d]\n" , Depth.name , Depth.value );
-	printf( "\t\t Running at depth d corresponds to solving on a 2^d x 2^d x 2^d\n" );
-	printf( "\t\t voxel grid.\n" );
-
-	printf( "\t[--%s <full depth>=%d]\n" , FullDepth.name , FullDepth.value );
-	printf( "\t\t This flag specifies the depth up to which the octree should be complete.\n" );
-
-	printf( "\t[--%s <depth at which to extract the voxel grid>=<%s>]\n" , VoxelDepth.name , Depth.name );
-
-	printf( "\t[--%s <scale factor>=%f]\n" , Scale.name , Scale.value );
-	printf( "\t\t Specifies the factor of the bounding cube that the input\n" );
-	printf( "\t\t samples should fit into.\n" );
-
-	printf( "\t[--%s <minimum number of samples per node>=%f]\n" , SamplesPerNode.name, SamplesPerNode.value );
-	printf( "\t\t This parameter specifies the minimum number of points that\n" );
-	printf( "\t\t should fall within an octree node.\n" );
-
-	printf( "\t[--%s <pull cut off>]\n" , Pull.name );
-	printf( "\t\t This parameter specifies the thresholld below which coarser colors functions should be used.\n" );
-
-#ifdef _OPENMP
-	printf( "\t[--%s <num threads>=%d]\n" , Threads.name , Threads.value );
-	printf( "\t\t This parameter specifies the number of threads across which\n" );
-	printf( "\t\t the solver should be parallelized.\n" );
-#endif // _OPENMP
-
-	printf( "\t[--%s]\n" , Confidence.name );
-	printf( "\t\t If this flag is enabled, the size of a sample's normals is\n" );
-	printf( "\t\t used as a confidence value, affecting the sample's\n" );
-	printf( "\t\t constribution to the reconstruction process.\n" );
-
-	printf( "\t[--%s]\n" , NormalWeights.name );
-	printf( "\t\t If this flag is enabled, the size of a sample's normals is\n" );
-	printf( "\t\t used as to modulate the interpolation weight.\n" );
-
-#if 0
-	printf( "\t[--%s <minimum depth>=%d]\n" , MinDepth.name , MinDepth.value );
-	printf( "\t\t This flag specifies the coarsest depth at which the system is to be solved.\n" );
-
-#ifdef _WIN32
-	printf( "\t[--%s]\n" , Performance.name );
-	printf( "\t\t If this flag is enabled, the running time and peak memory usage\n" );
-	printf( "\t\t is output after the reconstruction.\n" );
-#endif // _WIN32
-
-	printf( "\t[--%s]\n" , NoComments.name );
-	printf( "\t\t If this flag is enabled, the output file will not include comments.\n" );
-#endif
-	
-	printf( "\t[--%s]\n" , Double.name );
-	printf( "\t\t If this flag is enabled, the reconstruction will be performed with double-precision floats.\n" );
-
-	printf( "\t[--%s]\n" , Verbose.name );
-	printf( "\t\t If this flag is enabled, the progress of the reconstructor will be output to STDOUT.\n" );
-}
-Point3D< float > ReadASCIIColor( FILE* fp )
-{
-	unsigned char c[3];
-	if( fscanf( fp , " %c %c %c " , c+0 , c+1 , c+2 )!=3 ) fprintf( stderr , "[ERROR] Failed to read color\n" ) , exit( 0 );
-	return Point3D< float >( (float)c[0] , (float)c[1] , (float)c[2] );
-}
-
-PlyProperty PlyColorProperties[]=
-{
-	{ "r"     , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[0] ) ) , 0 , 0 , 0 , 0 },
-	{ "g"     , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[1] ) ) , 0 , 0 , 0 , 0 },
-	{ "b"     , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[2] ) ) , 0 , 0 , 0 , 0 },
-	{ "red"   , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[0] ) ) , 0 , 0 , 0 , 0 },
-	{ "green" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[1] ) ) , 0 , 0 , 0 , 0 },
-	{ "blue"  , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[2] ) ) , 0 , 0 , 0 , 0 },
-};
-bool ValidPlyColorProperties( const bool* props ){ return ( props[0] || props[3] ) && ( props[1] || props[4] ) && ( props[2] || props[5] ); }
-
-template< class Real , class Vertex >
-int Execute( int argc , char* argv[] )
-{
-	Reset< Real >();
-	int paramNum = sizeof(params)/sizeof(cmdLineReadable*);
-	std::vector< char* > comments;
-
-	if( Verbose.set ) echoStdout=1;
-
-	XForm4x4< Real > xForm , iXForm;
-	if( XForm.set )
-	{
-		FILE* fp = fopen( XForm.value , "r" );
-		if( !fp )
-		{
-			fprintf( stderr , "[WARNING] Could not read x-form from: %s\n" , XForm.value );
-			xForm = XForm4x4< Real >::Identity();
-		}
-		else
-		{
-			for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ )
-			{
-				float f;
-				fscanf( fp , " %f " , &f );
-				xForm(i,j) = (Real)f;
-			}
-			fclose( fp );
-		}
-	}
-	else xForm = XForm4x4< Real >::Identity();
-	iXForm = xForm.inverse();
-
-	DumpOutput2( comments , "Running Color Interpolator (Version 1.0)\n" );
-	char str[1024];
-	for( int i=0 ; i<paramNum ; i++ )
-		if( params[i]->set )
-		{
-			params[i]->writeValue( str );
-			if( strlen( str ) ) DumpOutput2( comments , "\t--%s %s\n" , params[i]->name , str );
-			else                DumpOutput2( comments , "\t--%s\n" , params[i]->name );
-		}
-
-	double t;
-	double tt=Time();
-
-	Octree< Real > tree;
-	tree.threads = Threads.value;
-	if( !InPoints.set || !InMesh.set )
-	{
-		ShowUsage( argv[0] );
-		return 0;
-	}
-	
-	OctNode< TreeNodeData >::SetAllocator( MEMORY_ALLOCATOR_BLOCK_SIZE );
-
-	t=Time();
-	int kernelDepth = KernelDepth.set ?  KernelDepth.value : Depth.value-2;
-	if( kernelDepth>Depth.value )
-	{
-		fprintf( stderr,"[ERROR] %s can't be greater than %s: %d <= %d\n" , KernelDepth.name , Depth.name , KernelDepth.value , Depth.value );
-		return EXIT_FAILURE;
-	}
-
-	double maxMemoryUsage;
-	t=Time() , tree.maxMemoryUsage=0;
-	typedef Octree< Real >::ProjectiveData< Point3D< float > > ProjectiveColor;
-	typename Octree< Real >::SparseNodeData< ProjectiveColor > colorData;
-	std::vector< Real > kernelDensityWeights;
-	OrientedPointStreamWithData< float , Point3D< float > >* pointStream;
-
-	char* ext = GetFileExtension( InPoints.value );
-	if     ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStreamWithData< float , Point3D< float > >( InPoints.value );
-	else if( !strcasecmp( ext , "ply"   ) ) pointStream = new    PLYOrientedPointStreamWithData< float , Point3D< float > >( InPoints.value , PlyColorProperties , 6 , ValidPlyColorProperties );
-	else                                    pointStream = new  ASCIIOrientedPointStreamWithData< float , Point3D< float > >( InPoints.value , ReadASCIIColor );
-	delete[] ext;
-	int pointCount = tree.template SetTree< float , Point3D< float > >( pointStream , MinDepth.value , Depth.value , FullDepth.value , kernelDepth , Real(SamplesPerNode.value) , Scale.value , Confidence.set , NormalWeights.set , kernelDensityWeights , colorData , xForm , BoundaryType.value );
-
-	DumpOutput2( comments , "#             Tree set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
-	DumpOutput( "Input Points: %d\n" , pointCount );
-	DumpOutput( "Depth/Nodes/Leaves: %d/%d/%d\n" , tree.tree.maxDepth() , tree.tree.nodes() , tree.tree.leaves() );
-	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) );
-
-	maxMemoryUsage = tree.maxMemoryUsage;
-	t=Time() , tree.maxMemoryUsage=0;
-	if( InMesh.set && Out.set )
-	{
-		BSplineData< 2 > fData;
-		fData.set( tree.tree.maxDepth() , BoundaryType.value );
-
-		std::vector< PlyColorVertex< float > > vertices;
-		std::vector< std::vector< int > > polygons;
-		int file_type;
-		PlyReadPolygons( InMesh.value , vertices , polygons , PlyVertex< float >::Properties , PlyVertex< float >::Components , file_type );
-
-#if 1
-#if 1
-		for( const OctNode< TreeNodeData >* n = tree.tree.nextNode() ; n!=NULL ; n=tree.tree.nextNode( n ) )
-		{
-			int idx = colorData.index( n );
-			if( idx>=0 ) colorData.data[idx] *= (Real)pow( Pull.value , n->depth() );
-		}
-		std::vector< ProjectiveColor > values( tree.tree.maxDepth()+1 );
-#pragma omp parallel for num_threads( Threads.value )
-		for( int i=0 ; i<vertices.size() ; i++ )
-		{
-			Point3D< Real > p = xForm * Point3D< Real >( vertices[i].point );
-			ProjectiveColor c = tree.Evaluate( colorData , p , &fData );
-			Point3D< float > _c = Point3D< float >( c );
-			for( int j=0 ; j<3 ; j++ ) vertices[i].color[j] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( _c[j]+0.5 ) ) );
-		}
-
-#else
-		int depth = tree.tree.maxDepth();
-		Pointer( Pointer( ProjectiveColor ) ) colors = AllocPointer< Pointer( ProjectiveColor ) >( depth+1 );
-		Pointer( int ) nodeCount = AllocPointer< int >( depth+1 );
-		for( int d=0 ; d<=depth ; d++ ) nodeCount[d] = 0;
-		for( const OctNode< TreeNodeData >* n = tree.tree.nextNode() ; n!=NULL ; n=tree.tree.nextNode( n ) ) nodeCount[ n->depth() ]++;
-		for( int d=1 ; d<=depth ; d++ ) nodeCount[d] += nodeCount[d-1];
-		for( int d=0 ; d<=depth ; d++ ) colors[d] = NewPointer< ProjectiveColor >( nodeCount[d] );
-
-		// Copy the finest level functions
-		for( const OctNode< TreeNodeData >* n = tree.tree.nextNode() ; n!=NULL ; n=tree.tree.nextNode( n ) )
-		{
-			int idx = colorData.index( n );
-			if( idx>=0 ) colors[depth][ n->nodeData.nodeIndex ] = colorData.data[idx];
-		}
-		// Down-sample to the coarser levels
-		for( int d=depth-1 ; d>=0 ; d-- ) for( const OctNode< TreeNodeData >* n=tree.tree.nextNode() ; n!=NULL ; )
-			if( n->depth()==d )
-			{
-				colors[d][ n->nodeData.nodeIndex ] = colors[d+1][ n->nodeData.nodeIndex ];
-				if( n->children ) for( int c=0 ; c<Cube::CORNERS ; c++ ) colors[d][ n->nodeData.nodeIndex ] += colors[d+1][ n->children[c].nodeData.nodeIndex ] / (Real)4;
-				n = tree.tree.nextBranch( n );
-			}
-			else n = tree.tree.nextNode( n );
-
-		double et = Time();
-		std::vector< typename OctNode< TreeNodeData >::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , Threads.value ) );
-		for( int i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( depth+1 );
-#pragma omp parallel for num_threads( Threads.value )
-		for( int i=0 ; i<vertices.size() ; i++ )
-		{
-			typename OctNode< TreeNodeData >::ConstNeighborKey3& neighborKey = neighborKeys[ omp_get_thread_num() ];
-			Point3D< Real > p = xForm * Point3D< Real >( vertices[i].point );
-			Real _depth , _weight;
-			tree.GetSampleDepthAndWeight( GetPointer( kernelDensityWeights ) , p , neighborKey , _depth , _weight );
-if( i==100 )
-{
-	printf( "Depth / Weight: %f / %f\n" , _depth , _weight );
-	for( int d=0 ; d<=depth ; d++ )
-	{
-		ProjectiveColor c = tree.Evaluate( colors[d] , p , &fData , d );
-		printf( "\t%d] %f\n" , d , c.w );
-	}
-}
-_depth -= 0.01;
-			ProjectiveColor c;
-			if( _depth>=depth ) c = tree.Evaluate( colors[depth] , p , &fData , depth );
-			else
-			{
-				int d = (int)floor( _depth );
-				Real dx = _depth - d;
-				c = tree.Evaluate( colors[d] , p , &fData , d ) * Real(1.-dx) + tree.Evaluate( colors[d+1] , p , &fData , d+1 ) * Real(dx);
-			}
-
-			Point3D< float > _c = Point3D< float >( c );
-			for( int j=0 ; j<3 ; j++ ) vertices[i].color[j] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( _c[j]+0.5 ) ) );
-		}
-		et = Time()-et;
-
-		for( int d=0 ; d<=depth ; d++ ) FreePointer( colors[d] );
-		DeletePointer( colors );
-		FreePointer( nodeCount );
-#endif
-#else
-		if( Pull.value>0 )
-		{
-			// Transform the sparse color representation into a dense one
-			Pointer( ProjectiveColor ) colors = new ProjectiveColor[ tree.tree.nodes() ];
-			for( const OctNode< TreeNodeData >* n = tree.tree.nextNode() ; n!=NULL ; n=tree.tree.nextNode( n ) )
-			{
-				int idx = colorData.index( n );
-				if( idx>=0 ) colors[ n->nodeData.nodeIndex ] = colorData.data[idx];
-			}
-
-			std::vector< ProjectiveColor > _colors( vertices.size() );
-			int depth = tree.tree.maxDepth();
-
-			// Iterate from fine to coarse, sampling color data as needed
-			for( int d=depth ; d>=1 ; d-- )
-			{
-
-				// Accumulate the color contribution if the weight is to low
-				double et = Time();
-#pragma omp parallel for num_threads( Threads.value )
-				for( int i=0 ; i<vertices.size() ; i++ ) if( _colors[i].w<Pull.value ) 
-				{
-					Point3D< Real > p = xForm * Point3D< Real >( vertices[i].point );
-					_colors[i] += tree.Evaluate( colors , p , &fData , d );
-				}
-				et = Time()-et;
-
-				// Pull the signal to the next coarser level
-				double ct = Time();
-				const OctNode< TreeNodeData >* n = tree.tree.nextNode();
-				while( n )
-				{
-					if( n->depth()==d-1 )
-					{
-						if( n->children ) for( int c=0 ; c<Cube::CORNERS ; c++ ) colors[ n->nodeData.nodeIndex ] += colors[ n->children[c].nodeData.nodeIndex ];
-						n = tree.tree.nextBranch( n );
-					}
-					else n = tree.tree.nextNode( n );
-				}
-				ct = Time()-ct;
-
-				if( Verbose.set ) printf( "Depth[%d] Eval./Compress: %.2f / %.2f (s)\n" , d , et , ct );
-			}
-			tree.MemoryUsage();
-			delete[] colors;
-
-#pragma omp parallel for num_threads( Threads.value )
-			for( int i=0 ; i<vertices.size() ; i++ )
-			{
-				Point3D< float > c = Point3D< float >( _colors[i] );
-				for( int j=0 ; j<3 ; j++ ) vertices[i].color[j] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( c[j]+0.5 ) ) );
-			}
-		}
-		else
-		{
-#pragma omp parallel for num_threads( Threads.value )
-			for( int i=0 ; i<vertices.size() ; i++ )
-			{
-				Point3D< Real > p = xForm * Point3D< Real >( vertices[i].point );
-				ProjectiveColor c = tree.Evaluate( colorData , p , &fData );
-				Point3D< float > _c = Point3D< float >( c );
-				for( int j=0 ; j<3 ; j++ ) vertices[i].color[j] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( _c[j]+0.5 ) ) );
-			}
-		}
-#endif
-
-		DumpOutput2( comments , "#       Sampled colors in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.MemoryUsage() );
-
-		PlyWritePolygons( Out.value , vertices , polygons , PlyColorVertex< float >::Properties , PlyColorVertex< float >::Components , file_type );
-		if( NoComments.set ) PlyWritePolygons( Out.value , vertices , polygons , PlyColorVertex< float >::Properties , PlyColorVertex< float >::Components , file_type , NULL , 0 );
-		else                 PlyWritePolygons( Out.value , vertices , polygons , PlyColorVertex< float >::Properties , PlyColorVertex< float >::Components , file_type , &comments[0] , (int)comments.size() );
-	}
-
-	for( int i=0 ; i<comments.size() ; i++ ) delete[] comments[i];
-	
-	return 1;
-}
-
-#ifdef _WIN32
-inline double to_seconds( const FILETIME& ft )
-{
-	const double low_to_sec=100e-9; // 100 nanoseconds
-	const double high_to_sec=low_to_sec*4294967296.0;
-	return ft.dwLowDateTime*low_to_sec+ft.dwHighDateTime*high_to_sec;
-}
-#endif // _WIN32
-
-int main( int argc , char* argv[] )
-{
-#if defined(WIN32) && defined(MAX_MEMORY_GB)
-	if( MAX_MEMORY_GB>0 )
-	{
-		SIZE_T peakMemory = 1;
-		peakMemory <<= 30;
-		peakMemory *= MAX_MEMORY_GB;
-		printf( "Limiting memory usage to %.2f GB\n" , float( peakMemory>>30 ) );
-		HANDLE h = CreateJobObject( NULL , NULL );
-		AssignProcessToJobObject( h , GetCurrentProcess() );
-
-		JOBOBJECT_EXTENDED_LIMIT_INFORMATION jeli = { 0 };
-		jeli.BasicLimitInformation.LimitFlags = JOB_OBJECT_LIMIT_JOB_MEMORY;
-		jeli.JobMemoryLimit = peakMemory;
-		if( !SetInformationJobObject( h , JobObjectExtendedLimitInformation , &jeli , sizeof( jeli ) ) )
-			fprintf( stderr , "Failed to set memory limit\n" );
-	}
-#endif // defined(WIN32) && defined(MAX_MEMORY_GB)
-	double t = Time();
-
-	cmdLineParse( argc-1 , &argv[1] , sizeof(params)/sizeof(cmdLineReadable*) , params , 1 );
-	if( Double.set ) Execute< double , PlyVertex< float > >( argc , argv );
-	else             Execute< float  , PlyVertex< float > >( argc , argv );
-#ifdef _WIN32
-	if( Performance.set )
-	{
-		HANDLE cur_thread=GetCurrentThread();
-		FILETIME tcreat, texit, tkernel, tuser;
-		if( GetThreadTimes( cur_thread , &tcreat , &texit , &tkernel , &tuser ) )
-			printf( "Time (Wall/User/Kernel): %.2f / %.2f / %.2f\n" , Time()-t , to_seconds( tuser ) , to_seconds( tkernel ) );
-		else printf( "Time: %.2f\n" , Time()-t );
-		HANDLE h = GetCurrentProcess();
-		PROCESS_MEMORY_COUNTERS pmc;
-		if( GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ) printf( "Peak Memory (MB): %d\n" , pmc.PeakWorkingSetSize>>20 );
-	}
-#endif // _WIN32
-	return EXIT_SUCCESS;
-}
diff --git a/Src/InterpolateColors.cpp b/Src/InterpolateColors.cpp
deleted file mode 100644
index 353c3f0..0000000
--- a/Src/InterpolateColors.cpp
+++ /dev/null
@@ -1,407 +0,0 @@
-
-/*
-Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without modification,
-are permitted provided that the following conditions are met:
-
-Redistributions of source code must retain the above copyright notice, this list of
-conditions and the following disclaimer. Redistributions in binary form must reproduce
-the above copyright notice, this list of conditions and the following disclaimer
-in the documentation and/or other materials provided with the distribution. 
-
-Neither the name of the Johns Hopkins University nor the names of its contributors
-may be used to endorse or promote products derived from this software without specific
-prior written permission. 
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
-EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES 
-OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
-SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
-INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
-TO, PROCUREMENT OF SUBSTITUTE  GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
-BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
-ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
-DAMAGE.
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <float.h>
-#ifdef _WIN32
-#include <Windows.h>
-#include <Psapi.h>
-#endif // _WIN32
-#include "MyTime.h"
-#include "MarchingCubes.h"
-#include "Octree.h"
-#include "SparseMatrix.h"
-#include "CmdLineParser.h"
-#include "PPolynomial.h"
-#include "Ply.h"
-#include "MemoryUsage.h"
-#ifdef _OPENMP
-#include "omp.h"
-#endif // _OPENMP
-void DumpOutput( const char* format , ... );
-#include "MultiGridOctreeData.h"
-void DumpOutput2( std::vector< char* >& comments , const char* format , ... );
-
-#define DEFAULT_FULL_DEPTH 5
-
-#define XSTR(x) STR(x)
-#define STR(x) #x
-#if DEFAULT_FULL_DEPTH
-#pragma message ( "[WARNING] Setting default full depth to " XSTR(DEFAULT_FULL_DEPTH) )
-#endif // DEFAULT_FULL_DEPTH
-
-#include <stdarg.h>
-char* outputFile=NULL;
-int echoStdout=0;
-void DumpOutput( const char* format , ... )
-{
-	if( outputFile )
-	{
-		FILE* fp = fopen( outputFile , "a" );
-		va_list args;
-		va_start( args , format );
-		vfprintf( fp , format , args );
-		fclose( fp );
-		va_end( args );
-	}
-	if( echoStdout )
-	{
-		va_list args;
-		va_start( args , format );
-		vprintf( format , args );
-		va_end( args );
-	}
-}
-void DumpOutput2( std::vector< char* >& comments , const char* format , ... )
-{
-	if( outputFile )
-	{
-		FILE* fp = fopen( outputFile , "a" );
-		va_list args;
-		va_start( args , format );
-		vfprintf( fp , format , args );
-		fclose( fp );
-		va_end( args );
-	}
-	if( echoStdout )
-	{
-		va_list args;
-		va_start( args , format );
-		vprintf( format , args );
-		va_end( args );
-	}
-	comments.push_back( new char[1024] );
-	char* str = comments.back();
-	va_list args;
-	va_start( args , format );
-	vsprintf( str , format , args );
-	va_end( args );
-	if( str[strlen(str)-1]=='\n' ) str[strlen(str)-1] = 0;
-}
-
-
-cmdLineString
-	InPoints( "inPoints" ) ,
-	InMesh( "inMesh" ) ,
-	Out( "out" ) ,
-	VoxelGrid( "voxel" ) ,
-	XForm( "xForm" );
-
-cmdLineReadable
-#ifdef _WIN32
-	Performance( "performance" ) ,
-#endif // _WIN32
-	NoComments( "noComments" ) ,
-	Confidence( "confidence" ) ,
-	NormalWeights( "nWeights" ) ,
-	Verbose( "verbose" ) ,
-	Double( "double" );
-
-cmdLineInt
-	Depth( "depth" , 8 ) ,
-	KernelDepth( "kernelDepth" ) ,
-	VoxelDepth( "voxelDepth" , -1 ) ,
-	FullDepth( "fullDepth" , DEFAULT_FULL_DEPTH ) ,
-	MinDepth( "minDepth" , 0 ) ,
-	BoundaryType( "boundary" , 1 ) ,
-	Threads( "threads" , omp_get_num_procs() );
-
-cmdLineFloat
-	SamplesPerNode( "samplesPerNode" , 1.f ) ,
-	Pull( "pull" , 1.f ) ,
-	Scale( "scale" , 1.1f );
-
-
-cmdLineReadable* params[] =
-{
-	&InPoints , &InMesh , &Depth , &Out , &XForm ,
-	&Pull , &Scale , &Verbose , &NoComments , &Double ,
-	&KernelDepth , &SamplesPerNode , &Confidence , &NormalWeights , &VoxelDepth ,
-	&VoxelGrid , &Threads ,
-	&BoundaryType ,
-	&FullDepth ,
-	&MinDepth ,
-#ifdef _WIN32
-	&Performance ,
-#endif // _WIN32
-};
-
-
-void ShowUsage(char* ex)
-{
-	printf( "Usage: %s\n" , ex );
-	printf( "\t --%s  <input points>\n" , InPoints.name );
-	printf( "\t --%s  <input mesh>\n" , InMesh.name );
-
-	printf( "\t[--%s <ouput triangle mesh>]\n" , Out.name );
-	printf( "\t[--%s <ouput voxel grid>]\n" , VoxelGrid.name );
-
-	printf( "\t[--%s <maximum reconstruction depth>=%d]\n" , Depth.name , Depth.value );
-	printf( "\t\t Running at depth d corresponds to solving on a 2^d x 2^d x 2^d\n" );
-	printf( "\t\t voxel grid.\n" );
-
-	printf( "\t[--%s <full depth>=%d]\n" , FullDepth.name , FullDepth.value );
-	printf( "\t\t This flag specifies the depth up to which the octree should be complete.\n" );
-
-	printf( "\t[--%s <depth at which to extract the voxel grid>=<%s>]\n" , VoxelDepth.name , Depth.name );
-
-	printf( "\t[--%s <scale factor>=%f]\n" , Scale.name , Scale.value );
-	printf( "\t\t Specifies the factor of the bounding cube that the input\n" );
-	printf( "\t\t samples should fit into.\n" );
-
-	printf( "\t[--%s <minimum number of samples per node>=%f]\n" , SamplesPerNode.name, SamplesPerNode.value );
-	printf( "\t\t This parameter specifies the minimum number of points that\n" );
-	printf( "\t\t should fall within an octree node.\n" );
-
-	printf( "\t[--%s <pull cut off>]\n" , Pull.name );
-	printf( "\t\t This parameter specifies the thresholld below which coarser colors functions should be used.\n" );
-
-#ifdef _OPENMP
-	printf( "\t[--%s <num threads>=%d]\n" , Threads.name , Threads.value );
-	printf( "\t\t This parameter specifies the number of threads across which\n" );
-	printf( "\t\t the solver should be parallelized.\n" );
-#endif // _OPENMP
-
-	printf( "\t[--%s]\n" , Confidence.name );
-	printf( "\t\t If this flag is enabled, the size of a sample's normals is\n" );
-	printf( "\t\t used as a confidence value, affecting the sample's\n" );
-	printf( "\t\t constribution to the reconstruction process.\n" );
-
-	printf( "\t[--%s]\n" , NormalWeights.name );
-	printf( "\t\t If this flag is enabled, the size of a sample's normals is\n" );
-	printf( "\t\t used as to modulate the interpolation weight.\n" );
-
-#if 0
-	printf( "\t[--%s <minimum depth>=%d]\n" , MinDepth.name , MinDepth.value );
-	printf( "\t\t This flag specifies the coarsest depth at which the system is to be solved.\n" );
-
-#ifdef _WIN32
-	printf( "\t[--%s]\n" , Performance.name );
-	printf( "\t\t If this flag is enabled, the running time and peak memory usage\n" );
-	printf( "\t\t is output after the reconstruction.\n" );
-#endif // _WIN32
-
-	printf( "\t[--%s]\n" , NoComments.name );
-	printf( "\t\t If this flag is enabled, the output file will not include comments.\n" );
-#endif
-	
-	printf( "\t[--%s]\n" , Double.name );
-	printf( "\t\t If this flag is enabled, the reconstruction will be performed with double-precision floats.\n" );
-
-	printf( "\t[--%s]\n" , Verbose.name );
-	printf( "\t\t If this flag is enabled, the progress of the reconstructor will be output to STDOUT.\n" );
-}
-Point3D< float > ReadASCIIColor( FILE* fp )
-{
-	unsigned char c[3];
-	if( fscanf( fp , " %c %c %c " , c+0 , c+1 , c+2 )!=3 ) fprintf( stderr , "[ERROR] Failed to read color\n" ) , exit( 0 );
-	return Point3D< float >( (float)c[0] , (float)c[1] , (float)c[2] );
-}
-
-PlyProperty PlyColorProperties[]=
-{
-	{ "r"     , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[0] ) ) , 0 , 0 , 0 , 0 },
-	{ "g"     , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[1] ) ) , 0 , 0 , 0 , 0 },
-	{ "b"     , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[2] ) ) , 0 , 0 , 0 , 0 },
-	{ "red"   , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[0] ) ) , 0 , 0 , 0 , 0 },
-	{ "green" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[1] ) ) , 0 , 0 , 0 , 0 },
-	{ "blue"  , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[2] ) ) , 0 , 0 , 0 , 0 },
-};
-bool ValidPlyColorProperties( const bool* props ){ return ( props[0] || props[3] ) && ( props[1] || props[4] ) && ( props[2] || props[5] ); }
-
-template< class Real , class Vertex >
-int Execute( int argc , char* argv[] )
-{
-	Reset< Real >();
-	int paramNum = sizeof(params)/sizeof(cmdLineReadable*);
-	std::vector< char* > comments;
-
-	if( Verbose.set ) echoStdout=1;
-
-	XForm4x4< Real > xForm , iXForm;
-	if( XForm.set )
-	{
-		FILE* fp = fopen( XForm.value , "r" );
-		if( !fp )
-		{
-			fprintf( stderr , "[WARNING] Could not read x-form from: %s\n" , XForm.value );
-			xForm = XForm4x4< Real >::Identity();
-		}
-		else
-		{
-			for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ )
-			{
-				float f;
-				fscanf( fp , " %f " , &f );
-				xForm(i,j) = (Real)f;
-			}
-			fclose( fp );
-		}
-	}
-	else xForm = XForm4x4< Real >::Identity();
-	iXForm = xForm.inverse();
-
-	DumpOutput2( comments , "Running Color Interpolator (Version 1.0)\n" );
-	char str[1024];
-	for( int i=0 ; i<paramNum ; i++ )
-		if( params[i]->set )
-		{
-			params[i]->writeValue( str );
-			if( strlen( str ) ) DumpOutput2( comments , "\t--%s %s\n" , params[i]->name , str );
-			else                DumpOutput2( comments , "\t--%s\n" , params[i]->name );
-		}
-
-	double t;
-	double tt=Time();
-
-	Octree< Real > tree;
-	tree.threads = Threads.value;
-	if( !InPoints.set || !InMesh.set )
-	{
-		ShowUsage( argv[0] );
-		return 0;
-	}
-	
-	OctNode< TreeNodeData >::SetAllocator( MEMORY_ALLOCATOR_BLOCK_SIZE );
-
-	t=Time();
-	int kernelDepth = KernelDepth.set ?  KernelDepth.value : Depth.value-2;
-	if( kernelDepth>Depth.value )
-	{
-		fprintf( stderr,"[ERROR] %s can't be greater than %s: %d <= %d\n" , KernelDepth.name , Depth.name , KernelDepth.value , Depth.value );
-		return EXIT_FAILURE;
-	}
-
-	double maxMemoryUsage;
-	t=Time() , tree.maxMemoryUsage=0;
-	typedef Octree< Real >::ProjectiveData< Point3D< float > > ProjectiveColor;
-	typename Octree< Real >::SparseNodeData< ProjectiveColor > colorData;
-	std::vector< Real > kernelDensityWeights;
-	OrientedPointStreamWithData< float , Point3D< float > >* pointStream;
-
-	char* ext = GetFileExtension( InPoints.value );
-	if     ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStreamWithData< float , Point3D< float > >( InPoints.value );
-	else if( !strcasecmp( ext , "ply"   ) ) pointStream = new    PLYOrientedPointStreamWithData< float , Point3D< float > >( InPoints.value , PlyColorProperties , 6 , ValidPlyColorProperties );
-	else                                    pointStream = new  ASCIIOrientedPointStreamWithData< float , Point3D< float > >( InPoints.value , ReadASCIIColor );
-	delete[] ext;
-	int pointCount = tree.template SetTree< float , Point3D< float > >( pointStream , MinDepth.value , Depth.value , FullDepth.value , kernelDepth , Real(SamplesPerNode.value) , Scale.value , Confidence.set , NormalWeights.set , kernelDensityWeights , colorData , xForm , BoundaryType.value );
-
-	DumpOutput2( comments , "#             Tree set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
-	DumpOutput( "Input Points: %d\n" , pointCount );
-	DumpOutput( "Depth/Nodes/Leaves: %d/%d/%d\n" , tree.tree.maxDepth() , tree.tree.nodes() , tree.tree.leaves() );
-	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) );
-
-	maxMemoryUsage = tree.maxMemoryUsage;
-	t=Time() , tree.maxMemoryUsage=0;
-	if( InMesh.set && Out.set )
-	{
-		BSplineData< 2 > fData;
-		fData.set( tree.tree.maxDepth() , BoundaryType.value );
-
-		std::vector< PlyColorVertex< float > > vertices;
-		std::vector< std::vector< int > > polygons;
-		int file_type;
-		PlyReadPolygons( InMesh.value , vertices , polygons , PlyVertex< float >::Properties , PlyVertex< float >::Components , file_type );
-
-		for( const OctNode< TreeNodeData >* n = tree.tree.nextNode() ; n!=NULL ; n=tree.tree.nextNode( n ) )
-		{
-			int idx = colorData.index( n );
-			if( idx>=0 ) colorData.data[idx] *= (Real)pow( Pull.value , n->depth() );
-		}
-		std::vector< ProjectiveColor > values( tree.tree.maxDepth()+1 );
-#pragma omp parallel for num_threads( Threads.value )
-		for( int i=0 ; i<vertices.size() ; i++ )
-		{
-			Point3D< Real > p = xForm * Point3D< Real >( vertices[i].point );
-			ProjectiveColor c = tree.Evaluate( colorData , p , &fData );
-			Point3D< float > _c = Point3D< float >( c );
-			for( int j=0 ; j<3 ; j++ ) vertices[i].color[j] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( _c[j]+0.5 ) ) );
-		}
-		DumpOutput2( comments , "#       Sampled colors in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.MemoryUsage() );
-
-		PlyWritePolygons( Out.value , vertices , polygons , PlyColorVertex< float >::Properties , PlyColorVertex< float >::Components , file_type );
-		if( NoComments.set ) PlyWritePolygons( Out.value , vertices , polygons , PlyColorVertex< float >::Properties , PlyColorVertex< float >::Components , file_type , NULL , 0 );
-		else                 PlyWritePolygons( Out.value , vertices , polygons , PlyColorVertex< float >::Properties , PlyColorVertex< float >::Components , file_type , &comments[0] , (int)comments.size() );
-	}
-
-	for( int i=0 ; i<comments.size() ; i++ ) delete[] comments[i];
-	
-	return 1;
-}
-
-#ifdef _WIN32
-inline double to_seconds( const FILETIME& ft )
-{
-	const double low_to_sec=100e-9; // 100 nanoseconds
-	const double high_to_sec=low_to_sec*4294967296.0;
-	return ft.dwLowDateTime*low_to_sec+ft.dwHighDateTime*high_to_sec;
-}
-#endif // _WIN32
-
-int main( int argc , char* argv[] )
-{
-#if defined(WIN32) && defined(MAX_MEMORY_GB)
-	if( MAX_MEMORY_GB>0 )
-	{
-		SIZE_T peakMemory = 1;
-		peakMemory <<= 30;
-		peakMemory *= MAX_MEMORY_GB;
-		printf( "Limiting memory usage to %.2f GB\n" , float( peakMemory>>30 ) );
-		HANDLE h = CreateJobObject( NULL , NULL );
-		AssignProcessToJobObject( h , GetCurrentProcess() );
-
-		JOBOBJECT_EXTENDED_LIMIT_INFORMATION jeli = { 0 };
-		jeli.BasicLimitInformation.LimitFlags = JOB_OBJECT_LIMIT_JOB_MEMORY;
-		jeli.JobMemoryLimit = peakMemory;
-		if( !SetInformationJobObject( h , JobObjectExtendedLimitInformation , &jeli , sizeof( jeli ) ) )
-			fprintf( stderr , "Failed to set memory limit\n" );
-	}
-#endif // defined(WIN32) && defined(MAX_MEMORY_GB)
-	double t = Time();
-
-	cmdLineParse( argc-1 , &argv[1] , sizeof(params)/sizeof(cmdLineReadable*) , params , 1 );
-	if( Double.set ) Execute< double , PlyVertex< float > >( argc , argv );
-	else             Execute< float  , PlyVertex< float > >( argc , argv );
-#ifdef _WIN32
-	if( Performance.set )
-	{
-		HANDLE cur_thread=GetCurrentThread();
-		FILETIME tcreat, texit, tkernel, tuser;
-		if( GetThreadTimes( cur_thread , &tcreat , &texit , &tkernel , &tuser ) )
-			printf( "Time (Wall/User/Kernel): %.2f / %.2f / %.2f\n" , Time()-t , to_seconds( tuser ) , to_seconds( tkernel ) );
-		else printf( "Time: %.2f\n" , Time()-t );
-		HANDLE h = GetCurrentProcess();
-		PROCESS_MEMORY_COUNTERS pmc;
-		if( GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ) printf( "Peak Memory (MB): %d\n" , pmc.PeakWorkingSetSize>>20 );
-	}
-#endif // _WIN32
-	return EXIT_SUCCESS;
-}
diff --git a/Src/MultiGridOctreeData.IsoSurface.inl b/Src/MultiGridOctreeData.IsoSurface.inl
index 930688c..124ddd1 100644
--- a/Src/MultiGridOctreeData.IsoSurface.inl
+++ b/Src/MultiGridOctreeData.IsoSurface.inl
@@ -136,7 +136,7 @@ void Octree< Real >::XSliceValues< Vertex >::reset( void )
 }
 
 template< class Real >
-template< class Vertex , class _Vertex >
+template< class Vertex >
 void Octree< Real >::GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , ConstPointer( Real ) solution , Real isoValue , CoredMeshData< Vertex >& mesh , bool nonLinearFit , bool addBarycenter , bool polygonMesh )
 {
 	typename BSplineData< 2 >::template CornerEvaluator< 2 > evaluator;
@@ -205,7 +205,7 @@ void Octree< Real >::GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights
 			SetXSliceIsoEdges( d , o-1 , slabValues , threads );
 
 			// Add the triangles
-			SetIsoSurface< Vertex , _Vertex >( d , o-1 , slabValues[d].sliceValues(o-1) , slabValues[d].sliceValues(o) , slabValues[d].xSliceValues(o-1) , mesh , polygonMesh , addBarycenter , vertexOffset , threads );
+			SetIsoSurface( d , o-1 , slabValues[d].sliceValues(o-1) , slabValues[d].sliceValues(o) , slabValues[d].xSliceValues(o-1) , mesh , polygonMesh , addBarycenter , vertexOffset , threads );
 
 			if( o&1 ) break;
 		}
@@ -791,7 +791,7 @@ void Octree< Real >::SetXSliceIsoEdges( int depth , int slab , std::vector< Slab
 	}
 }
 template< class Real >
-template< class Vertex , class _Vertex >
+template< class Vertex >
 void Octree< Real >::SetIsoSurface( int depth , int offset , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , const XSliceValues< Vertex >& xValues , CoredMeshData< Vertex >& mesh , bool polygonMesh , bool addBarycenter , int& vOffset , int threads )
 {
 	std::vector< std::pair< int , Vertex > > polygon;
@@ -902,7 +902,7 @@ void Octree< Real >::SetIsoSurface( int depth , int offset , const SliceValues<
 						else if( ( iter=xValues.edgeVertexMap.find( key ) )!=xValues.edgeVertexMap.end() ) polygon[k] = iter->second;
 						else fprintf( stderr , "[ERROR] Couldn't find vertex in edge map\n" ) , exit( 0 );
 					}
-					AddIsoPolygons< Vertex , _Vertex >( mesh , polygon , polygonMesh , addBarycenter , vOffset );
+					AddIsoPolygons( mesh , polygon , polygonMesh , addBarycenter , vOffset );
 				}
 			}
 		}
@@ -1061,7 +1061,7 @@ bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , c
 }
 
 template< class Real >
-template< class Vertex , class _Vertex>
+template< class Vertex >
 int Octree< Real >::AddIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector< std::pair< int , Vertex > >& polygon , bool polygonMesh , bool addBarycenter , int& vOffset )
 {
 	if( polygonMesh )
@@ -1087,9 +1087,9 @@ int Octree< Real >::AddIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector<
 		if( isCoplanar )
 		{
 			Vertex c;
-			_Vertex _c;
+			typename Vertex::Wrapper _c;
 			_c *= 0;
-			for( int i=0 ; i<(int)polygon.size() ; i++ ) _c += _Vertex( polygon[i].second );
+			for( int i=0 ; i<(int)polygon.size() ; i++ ) _c += typename Vertex::Wrapper( polygon[i].second );
 			_c /= Real( polygon.size() );
 			c = Vertex( _c );
 			int cIdx;
diff --git a/Src/MultiGridOctreeData.h b/Src/MultiGridOctreeData.h
index f3ed4dc..fe02875 100644
--- a/Src/MultiGridOctreeData.h
+++ b/Src/MultiGridOctreeData.h
@@ -30,6 +30,7 @@ DAMAGE.
 #define MULTI_GRID_OCTREE_DATA_INCLUDED
 
 #define NEW_CODE 1
+#define NEW_NEW_CODE 1
 
 //#define MAX_MEMORY_GB 15
 #define MAX_MEMORY_GB 0
@@ -408,10 +409,10 @@ protected:
 	template< class Vertex >
 	void SetXSliceIsoEdges( int depth , int slice , std::vector< SlabValues< Vertex > >& slabValues , int threads );
 
-	template< class Vertex , class _Vertex >
+	template< class Vertex >
 	void SetIsoSurface( int depth , int offset , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , const XSliceValues< Vertex >& xValues , CoredMeshData< Vertex >& mesh , bool polygonMesh , bool addBarycenter , int& vOffset , int threads );
 
-	template< class Vertex , class _Vertex >
+	template< class Vertex >
 	static int AddIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector< std::pair< int , Vertex > >& polygon , bool polygonMesh , bool addBarycenter , int& vOffset );
 
 	template< class Vertex >
@@ -474,7 +475,7 @@ public:
 	Pointer( Real ) SolveSystem( SparseNodeData< PointData >& pointInfo , Pointer( Real ) constraints , bool showResidual , int iters , int maxSolveDepth , int cgDepth=0 , double cgAccuracy=0 );
 
 	Real GetIsoValue( ConstPointer( Real ) solution , const std::vector< Real >& centerWeights );
-	template< class Vertex , class _Vertex >
+	template< class Vertex >
 	void GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , ConstPointer( Real ) solution , Real isoValue , CoredMeshData< Vertex >& mesh , bool nonLinearFit=true , bool addBarycenter=false , bool polygonMesh=false );
 };
 template< class Real >
diff --git a/Src/MultiGridOctreeData.inl b/Src/MultiGridOctreeData.inl
index 70a6579..ad02751 100644
--- a/Src/MultiGridOctreeData.inl
+++ b/Src/MultiGridOctreeData.inl
@@ -574,7 +574,6 @@ int Octree< Real >::SetTree( OrientedPointStream< PointReal >* pointStream , int
 #endif
 
 	_fData.set( maxDepth , _boundaryType );
-
 	_minDepth = minDepth;
 	_fullDepth = fullDepth;
 	double pointWeightSum = 0;
@@ -2751,6 +2750,8 @@ template< class V >
 V Octree< Real >::getCornerValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey , const TreeOctNode* node , int corner , ConstPointer( V ) solution , ConstPointer( V ) metSolution , const typename BSplineData< 2 >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 > stencils[8] , bool isInterior ) const
 {
 	V value(0);
+	// [NOTE] For screening, we force the samples to have value 0.0, so Dirichlet boundary constraints are implemented
+	// by having the points on the boundary have value -0.5. Hence the offsetting.
 	if( _boundaryType==-1 ) value = -0.5;
 	int d , off[3];
 	node->depthAndOffset( d , off );
@@ -3142,6 +3143,8 @@ Pointer( V ) Octree< Real >::Evaluate( ConstPointer( V ) coefficients , int& res
 						vTables.valueTable[ idx[2] + z*vTables.functionCount ];
 				}
 	}
+	// [NOTE] For screening, we force the samples to have value 0.0, so Dirichlet boundary constraints are implemented
+	// by having the points on the boundary have value -0.5. Hence the offsetting.
 	if( _boundaryType==-1 ) for( int i=0 ; i<res*res*res ; i++ ) values[i] -= Real(0.5);
 	for( int i=0 ; i<res*res*res ; i++ ) values[i] -= isoValue;
 
diff --git a/Src/Ply.h b/Src/Ply.h
index 5a1b8f0..699381f 100644
--- a/Src/Ply.h
+++ b/Src/Ply.h
@@ -33,6 +33,8 @@
 #ifndef __PLY_H__
 #define __PLY_H__
 
+#define USE_PLY_WRAPPER 1
+
 #ifndef WIN32
 #define _strdup strdup
 #endif
@@ -239,12 +241,22 @@ static PlyProperty face_props[] =
 	{ _strdup( "vertex_indices" ) , PLY_INT , PLY_INT , offsetof( PlyFace , vertices ) , 1 , PLY_UCHAR, PLY_UCHAR , offsetof(PlyFace,nr_vertices) },
 };
 
+
+///////////////////
+// PlyVertexType //
+///////////////////
+
+// The "Wrapper" class indicates the class to cast to/from in order to support linear operations.
 template< class Real >
 class PlyVertex
 {
 public:
-	const static int Components=3;
-	static PlyProperty Properties[];
+	typedef PlyVertex Wrapper;
+
+	const static int ReadComponents=3;
+	const static int WriteComponents=3;
+	static PlyProperty ReadProperties[];
+	static PlyProperty WriteProperties[];
 
 	Point3D< Real > point;
 
@@ -260,7 +272,13 @@ public:
 	template< class _Real > PlyVertex& operator /= ( _Real s ) { point /= s ; return *this; }
 };
 template< class Real , class _Real > PlyVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyVertex< Real > v ) { return PlyVertex< Real >( xForm * v.point ); }
-template< class Real > PlyProperty PlyVertex< Real >::Properties[]=
+template< class Real > PlyProperty PlyVertex< Real >::ReadProperties[]=
+{
+	{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 }
+};
+template< class Real > PlyProperty PlyVertex< Real >::WriteProperties[]=
 {
 	{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
 	{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
@@ -270,8 +288,12 @@ template< class Real >
 class PlyValueVertex
 {
 public:
-	const static int Components=4;
-	static PlyProperty Properties[];
+	typedef PlyValueVertex Wrapper;
+
+	const static int ReadComponents=4;
+	const static int WriteComponents=4;
+	static PlyProperty ReadProperties[];
+	static PlyProperty WriteProperties[];
 
 	Point3D<Real> point;
 	Real value;
@@ -288,8 +310,14 @@ public:
 	template< class _Real > PlyValueVertex& operator /= ( _Real s ) { point /= s , value /= Real(s) ; return *this; }
 };
 template< class Real , class _Real > PlyValueVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyValueVertex< Real > v ) { return PlyValueVertex< Real >( xForm * v.point , v.value ); }
-template< class Real >
-PlyProperty PlyValueVertex< Real >::Properties[]=
+template< class Real > PlyProperty PlyValueVertex< Real >::ReadProperties[]=
+{
+	{ _strdup( "x"     ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "y"     ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "z"     ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "value" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , value           ) ) , 0 , 0 , 0 , 0 }
+};
+template< class Real > PlyProperty PlyValueVertex< Real >::WriteProperties[]=
 {
 	{ _strdup( "x"     ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
 	{ _strdup( "y"     ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
@@ -300,8 +328,12 @@ template< class Real >
 class PlyOrientedVertex
 {
 public:
-	const static int Components=6;
-	static PlyProperty Properties[];
+	typedef PlyOrientedVertex Wrapper;
+
+	const static int ReadComponents=6;
+	const static int WriteComponents=6;
+	static PlyProperty ReadProperties[];
+	static PlyProperty WriteProperties[];
 
 	Point3D<Real> point , normal;
 
@@ -317,8 +349,16 @@ public:
 	template< class _Real > PlyOrientedVertex& operator /= ( _Real s ) { point /= s , normal /= s ; return *this; }
 };
 template< class Real , class _Real > PlyOrientedVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyOrientedVertex< Real > v ) { return PlyOrientedVertex< Real >( xForm * v.point , xForm.inverse().transpose() * v.normal ); }
-template< class Real >
-PlyProperty PlyOrientedVertex< Real >::Properties[]=
+template< class Real > PlyProperty PlyOrientedVertex< Real >::ReadProperties[]=
+{
+	{ _strdup( "x"  ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex ,  point.coords[0] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "y"  ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex ,  point.coords[1] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "z"  ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex ,  point.coords[2] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "nx" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[0] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "ny" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[1] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "nz" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[2] ) ) , 0 , 0 , 0 , 0 }
+};
+template< class Real > PlyProperty PlyOrientedVertex< Real >::WriteProperties[]=
 {
 	{ _strdup( "x"  ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex ,  point.coords[0] ) ) , 0 , 0 , 0 , 0 },
 	{ _strdup( "y"  ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex ,  point.coords[1] ) ) , 0 , 0 , 0 , 0 },
@@ -331,8 +371,36 @@ template< class Real >
 class PlyColorVertex
 {
 public:
-	const static int Components=6;
-	static PlyProperty Properties[];
+	struct _PlyColorVertex
+	{
+		Point3D< Real > point , color;
+		_PlyColorVertex( void ) { ; }
+		_PlyColorVertex( Point3D< Real > p , Point3D< Real > c ) : point(p) , color(c) { ; }
+		_PlyColorVertex( PlyColorVertex< Real > p ){ point = p.point ; for( int c=0 ; c<3 ; c++ ) color[c] = (Real) p.color[c]; }
+		operator PlyColorVertex< Real > ()
+		{
+			PlyColorVertex< Real > p;
+			p.point = point;
+			for( int c=0 ; c<3 ; c++ ) p.color[c] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[c]+0.5 ) ) );
+			return p;
+		}
+
+	  	_PlyColorVertex operator + ( _PlyColorVertex p ) const { return _PlyColorVertex( point+p.point , color+p.color ); }
+		_PlyColorVertex operator - ( _PlyColorVertex p ) const { return _PlyColorVertex( point-p.value , color-p.color ); }
+		template< class _Real > _PlyColorVertex operator * ( _Real s ) const { return _PlyColorVertex( point*s , color*s ); }
+		template< class _Real > _PlyColorVertex operator / ( _Real s ) const { return _PlyColorVertex( point/s , color/s ); }
+		_PlyColorVertex& operator += ( _PlyColorVertex p ) { point += p.point , color += p.color ; return *this; }
+		_PlyColorVertex& operator -= ( _PlyColorVertex p ) { point -= p.point , color -= p.color ; return *this; }
+		template< class _Real > _PlyColorVertex& operator *= ( _Real s ) { point *= s , color *= s ; return *this; }
+		template< class _Real > _PlyColorVertex& operator /= ( _Real s ) { point /= s , color /= s ; return *this; }
+	};
+
+	typedef _PlyColorVertex Wrapper;
+
+	const static int ReadComponents=9;
+	const static int WriteComponents=6;
+	static PlyProperty ReadProperties[];
+	static PlyProperty WriteProperties[];
 
 	Point3D< Real > point;
 	unsigned char color[3];
@@ -345,8 +413,19 @@ public:
 };
 template< class Real , class _Real > PlyColorVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyColorVertex< Real > v ) { return PlyColorVertex< Real >( xForm * v.point , v.color ); }
 
-template< class Real >
-PlyProperty PlyColorVertex< Real >::Properties[]=
+template< class Real > PlyProperty PlyColorVertex< Real >::ReadProperties[]=
+{
+	{ _strdup( "x"     ) , PLYType<          Real >() , PLYType<          Real >(), int( offsetof( PlyColorVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "y"     ) , PLYType<          Real >() , PLYType<          Real >(), int( offsetof( PlyColorVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "z"     ) , PLYType<          Real >() , PLYType<          Real >(), int( offsetof( PlyColorVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "red"   ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex ,        color[0] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex ,        color[1] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "blue"  ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex ,        color[2] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "r"     ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex ,        color[0] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "g"     ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex ,        color[1] ) ) , 0 , 0 , 0 , 0 },
+	{ _strdup( "b"     ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex ,        color[2] ) ) , 0 , 0 , 0 , 0 }
+};
+template< class Real > PlyProperty PlyColorVertex< Real >::WriteProperties[]=
 {
 	{ _strdup( "x"     ) , PLYType<          Real >() , PLYType<          Real >(), int( offsetof( PlyColorVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
 	{ _strdup( "y"     ) , PLYType<          Real >() , PLYType<          Real >(), int( offsetof( PlyColorVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
@@ -359,8 +438,38 @@ template< class Real >
 class PlyColorAndValueVertex
 {
 public:
-	const static int Components=7;
-	static PlyProperty Properties[];
+	struct _PlyColorAndValueVertex
+	{
+		Point3D< Real > point , color;
+		Real value;
+		_PlyColorAndValueVertex( void ) : value(0) { ; }
+		_PlyColorAndValueVertex( Point3D< Real > p , Point3D< Real > c , Real v ) : point(p) , color(c) , value(v) { ; }
+		_PlyColorAndValueVertex( PlyColorAndValueVertex< Real > p ){ point = p.point ; for( int c=0 ; c<3 ; c++ ) color[c] = (Real) p.color[c] ; value = p.value; }
+		operator PlyColorAndValueVertex< Real > ()
+		{
+			PlyColorAndValueVertex< Real > p;
+			p.point = point;
+			for( int c=0 ; c<3 ; c++ ) p.color[c] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[c]+0.5 ) ) );
+			p.value = value;
+			return p;
+		}
+
+	  	_PlyColorAndValueVertex operator + ( _PlyColorAndValueVertex p ) const { return _PlyColorAndValueVertex( point+p.point , color+p.color , value+p.value ); }
+		_PlyColorAndValueVertex operator - ( _PlyColorAndValueVertex p ) const { return _PlyColorAndValueVertex( point-p.value , color-p.color , value+p.value ); }
+		template< class _Real > _PlyColorAndValueVertex operator * ( _Real s ) const { return _PlyColorAndValueVertex( point*s , color*s , value*s ); }
+		template< class _Real > _PlyColorAndValueVertex operator / ( _Real s ) const { return _PlyColorAndValueVertex( point/s , color/s , value/s ); }
+		_PlyColorAndValueVertex& operator += ( _PlyColorAndValueVertex p ) { point += p.point , color += p.color , value += p.value ; return *this; }
+		_PlyColorAndValueVertex& operator -= ( _PlyColorAndValueVertex p ) { point -= p.point , color -= p.color , value -= p.value ; return *this; }
+		template< class _Real > _PlyColorAndValueVertex& operator *= ( _Real s ) { point *= s , color *= s , value *= (Real)s ; return *this; }
+		template< class _Real > _PlyColorAndValueVertex& operator /= ( _Real s ) { point /= s , color /= s , value /= (Real)s ; return *this; }
+	};
+
+	typedef _PlyColorAndValueVertex Wrapper;
+
+	const static int ReadComponents=10;
+	const static int WriteComponents=7;
+	static PlyProperty ReadProperties[];
+	static PlyProperty WriteProperties[];
 
 	Point3D< Real > point;
 	unsigned char color[3];
@@ -373,8 +482,20 @@ public:
 	PlyColorAndValueVertex( const Point3D< Real >& p , const unsigned char c[3] , Real v) { point = p , color[0] = c[0] , color[1] = c[1] , color[2] = c[2] , value = v; }
 };
 template< class Real , class _Real > PlyColorAndValueVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyColorAndValueVertex< Real > v ) { return PlyColorAndValueVertex< Real >( xForm * v.point , v.color , v.value ); }
-template< class Real >
-PlyProperty PlyColorAndValueVertex< Real >::Properties[]=
+template< class Real > PlyProperty PlyColorAndValueVertex< Real >::ReadProperties[]=
+{
+	{ _strdup( "x"     ) , PLYType<          Real >() , PLYType<          Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "y"     ) , PLYType<          Real >() , PLYType<          Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "z"     ) , PLYType<          Real >() , PLYType<          Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "value" ) , PLYType<          Real >() , PLYType<          Real >() , int( offsetof( PlyColorAndValueVertex ,        value    ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "red"   ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex ,        color[0] ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex ,        color[1] ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "blue"  ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex ,        color[2] ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "r"     ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex ,        color[0] ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "g"     ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex ,        color[1] ) ) , 0 , 0 , 0 , 0 } ,
+	{ _strdup( "b"     ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex ,        color[2] ) ) , 0 , 0 , 0 , 0 }
+};
+template< class Real > PlyProperty PlyColorAndValueVertex< Real >::WriteProperties[]=
 {
 	{ _strdup( "x"     ) , PLYType<          Real >() , PLYType<          Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 } ,
 	{ _strdup( "y"     ) , PLYType<          Real >() , PLYType<          Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 } ,
@@ -385,63 +506,100 @@ PlyProperty PlyColorAndValueVertex< Real >::Properties[]=
 	{ _strdup( "blue"  ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex ,        color[2] ) ) , 0 , 0 , 0 , 0 }
 };
 
-template< class Real >
-struct _PlyColorVertex
-{
-	Point3D< Real > point , color;
-	_PlyColorVertex( void ) { ; }
-	_PlyColorVertex( Point3D< Real > p , Point3D< Real > c ) : point(p) , color(c) { ; }
-	_PlyColorVertex( PlyColorVertex< Real > p ){ point = p.point ; for( int c=0 ; c<3 ; c++ ) color[c] = (Real) p.color[c]; }
-	operator PlyColorVertex< Real > ()
-	{
-		PlyColorVertex< Real > p;
-		p.point = point;
-		for( int c=0 ; c<3 ; c++ ) p.color[c] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[c]+0.5 ) ) );
-		return p;
-	}
-
-  	_PlyColorVertex operator + ( _PlyColorVertex p ) const { return _PlyColorVertex( point+p.point , color+p.color ); }
-	_PlyColorVertex operator - ( _PlyColorVertex p ) const { return _PlyColorVertex( point-p.value , color-p.color ); }
-	template< class _Real > _PlyColorVertex operator * ( _Real s ) const { return _PlyColorVertex( point*s , color*s ); }
-	template< class _Real > _PlyColorVertex operator / ( _Real s ) const { return _PlyColorVertex( point/s , color/s ); }
-	_PlyColorVertex& operator += ( _PlyColorVertex p ) { point += p.point , color += p.color ; return *this; }
-	_PlyColorVertex& operator -= ( _PlyColorVertex p ) { point -= p.point , color -= p.color ; return *this; }
-	template< class _Real > _PlyColorVertex& operator *= ( _Real s ) { point *= s , color *= s ; return *this; }
-	template< class _Real > _PlyColorVertex& operator /= ( _Real s ) { point /= s , color /= s ; return *this; }
-};
-template< class Real >
-struct _PlyColorAndValueVertex
-{
-	Point3D< Real > point , color;
-	Real value;
-	_PlyColorAndValueVertex( void ) : value(0) { ; }
-	_PlyColorAndValueVertex( Point3D< Real > p , Point3D< Real > c , Real v ) : point(p) , color(c) , value(v) { ; }
-	_PlyColorAndValueVertex( PlyColorAndValueVertex< Real > p ){ point = p.point ; for( int c=0 ; c<3 ; c++ ) color[c] = (Real) p.color[c] ; value = p.value; }
-	operator PlyColorAndValueVertex< Real > ()
-	{
-		PlyColorAndValueVertex< Real > p;
-		p.point = point;
-		for( int c=0 ; c<3 ; c++ ) p.color[c] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[c]+0.5 ) ) );
-		p.value = value;
-		return p;
-	}
-
-  	_PlyColorAndValueVertex operator + ( _PlyColorAndValueVertex p ) const { return _PlyColorAndValueVertex( point+p.point , color+p.color , value+p.value ); }
-	_PlyColorAndValueVertex operator - ( _PlyColorAndValueVertex p ) const { return _PlyColorAndValueVertex( point-p.value , color-p.color , value+p.value ); }
-	template< class _Real > _PlyColorAndValueVertex operator * ( _Real s ) const { return _PlyColorAndValueVertex( point*s , color*s , value*s ); }
-	template< class _Real > _PlyColorAndValueVertex operator / ( _Real s ) const { return _PlyColorAndValueVertex( point/s , color/s , value/s ); }
-	_PlyColorAndValueVertex& operator += ( _PlyColorAndValueVertex p ) { point += p.point , color += p.color , value += p.value ; return *this; }
-	_PlyColorAndValueVertex& operator -= ( _PlyColorAndValueVertex p ) { point -= p.point , color -= p.color , value -= p.value ; return *this; }
-	template< class _Real > _PlyColorAndValueVertex& operator *= ( _Real s ) { point *= s , color *= s , value *= (Real)s ; return *this; }
-	template< class _Real > _PlyColorAndValueVertex& operator /= ( _Real s ) { point /= s , color /= s , value /= (Real)s ; return *this; }
-};
-
 template< class Vertex , class Real >
 int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >*  mesh , int file_type , const Point3D< float >& translate , float scale , char** comments=NULL , int commentNum=0 , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() );
 
 template< class Vertex , class Real >
 int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >*  mesh , int file_type , char** comments=NULL , int commentNum=0 , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() );
 
+inline bool PlyReadHeader( char* fileName , PlyProperty* properties , int propertyNum , bool* readFlags , int& file_type )
+{
+	int nr_elems;
+	char **elist;
+	float version;
+	PlyFile* ply;
+	char* elem_name;
+	int num_elems;
+	int nr_props;
+	PlyProperty** plist;
+
+	ply = ply_open_for_reading( fileName , &nr_elems , &elist , &file_type , &version );
+	if( !ply ) return false;
+
+	for( int i=0 ; i<nr_elems ; i++ )
+	{
+		elem_name = elist[i];
+		plist = ply_get_element_description( ply , elem_name , &num_elems , &nr_props );
+		if( !plist )
+		{
+			for( int i=0 ; i<nr_elems ; i++ )
+			{
+				free( ply->elems[i]->name );
+				free( ply->elems[i]->store_prop );
+				for( int j=0 ; j<ply->elems[i]->nprops ; j++ )
+				{
+					free( ply->elems[i]->props[j]->name );
+					free( ply->elems[i]->props[j] );
+				}
+				free( ply->elems[i]->props );
+			}
+			for( int i=0 ; i<nr_elems ; i++ ) free( ply->elems[i] );
+			free( ply->elems );
+			for( int i=0 ; i<ply->num_comments ; i++ ) free( ply->comments[i] );
+			free( ply->comments );
+			for( int i=0 ; i<ply->num_obj_info ; i++ ) free( ply->obj_info[i] );
+			free( ply->obj_info );
+			ply_free_other_elements( ply->other_elems );
+			
+			for( int i=0 ; i<nr_elems ; i++ ) free( elist[i] );
+			free( elist );
+			ply_close( ply );
+			return 0;
+		}		
+		if( equal_strings( "vertex" , elem_name ) )
+			for( int i=0 ; i<propertyNum ; i++ )
+				if( readFlags ) readFlags[i] = ply_get_property( ply , elem_name , &properties[i] )!=0;
+
+		for( int j=0 ; j<nr_props ; j++ )
+		{
+			free( plist[j]->name );
+			free( plist[j] );
+		}
+		free( plist );
+	}  // for each type of element
+	
+	for( int i=0 ; i<nr_elems ; i++ )
+	{
+		free( ply->elems[i]->name );
+		free( ply->elems[i]->store_prop );
+		for( int j=0 ; j<ply->elems[i]->nprops ; j++ )
+		{
+			free( ply->elems[i]->props[j]->name );
+			free( ply->elems[i]->props[j] );
+		}
+		if( ply->elems[i]->props && ply->elems[i]->nprops ) free(ply->elems[i]->props);
+	}
+	for( int i=0 ; i<nr_elems ; i++ ) free(ply->elems[i]);
+	free( ply->elems) ;
+	for( int i=0 ; i<ply->num_comments ; i++ ) free( ply->comments[i] );
+	free( ply->comments );
+	for( int i=0 ; i<ply->num_obj_info ; i++ ) free( ply->obj_info[i] );
+	free( ply->obj_info );
+	ply_free_other_elements(ply->other_elems);
+	
+	
+	for( int i=0 ; i<nr_elems ; i++ ) free( elist[i] );
+	free( elist );
+	ply_close( ply );
+	return true;
+}
+inline bool PlyReadHeader( char* fileName , PlyProperty* properties , int propertyNum , bool* readFlags )
+{
+	int file_type;
+	return PlyReadHeader( fileName , properties , propertyNum , readFlags , file_type );
+}
+
+
 template<class Vertex>
 int PlyReadPolygons(char* fileName,
 					std::vector<Vertex>& vertices,std::vector<std::vector<int> >& polygons,
@@ -709,7 +867,7 @@ int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_
 	// describe vertex and face properties
 	//
 	ply_element_count( ply , "vertex" , nr_vertices );
-	for( int i=0 ; i<Vertex::Components ; i++ ) ply_describe_property( ply , "vertex" , &Vertex::Properties[i] );
+	for( int i=0 ; i<Vertex::WriteComponents ; i++ ) ply_describe_property( ply , "vertex" , &Vertex::WriteProperties[i] );
 	
 	ply_element_count( ply , "face" , nr_faces );
 	ply_describe_property( ply , "face" , &face_props[0] );
diff --git a/Src/PointStream.inl b/Src/PointStream.inl
index f59685c..4883d22 100644
--- a/Src/PointStream.inl
+++ b/Src/PointStream.inl
@@ -153,10 +153,10 @@ void PLYOrientedPointStream< Real >::reset( void )
 		{
 			foundVertices = true;
 			_pCount = num_elems , _pIdx = 0;
-			for( int i=0 ; i<PlyOrientedVertex< Real >::Components ; i++ ) 
-				if( !ply_get_property( _ply , elem_name , &(PlyOrientedVertex< Real >::Properties[i]) ) )
+			for( int i=0 ; i<PlyOrientedVertex< Real >::ReadComponents ; i++ ) 
+				if( !ply_get_property( _ply , elem_name , &(PlyOrientedVertex< Real >::ReadProperties[i]) ) )
 				{
-					fprintf( stderr , "[ERROR] Failed to find property in ply file: %s\n" , PlyOrientedVertex< Real >::Properties[i].name );
+					fprintf( stderr , "[ERROR] Failed to find property in ply file: %s\n" , PlyOrientedVertex< Real >::ReadProperties[i].name );
 					exit( 0 );
 				}
 		}
@@ -337,10 +337,10 @@ void PLYOrientedPointStreamWithData< Real , Data >::reset( void )
 		{
 			foundVertices = true;
 			_pCount = num_elems , _pIdx = 0;
-			for( int i=0 ; i<PlyOrientedVertex< Real >::Components ; i++ ) 
-				if( !ply_get_property( _ply , elem_name , &(PlyOrientedVertex< Real >::Properties[i]) ) )
+			for( int i=0 ; i<PlyOrientedVertex< Real >::ReadComponents ; i++ ) 
+				if( !ply_get_property( _ply , elem_name , &(PlyOrientedVertex< Real >::ReadProperties[i]) ) )
 				{
-					fprintf( stderr , "[ERROR] Failed to find property in ply file: %s\n" , PlyOrientedVertex< Real >::Properties[i].name );
+					fprintf( stderr , "[ERROR] Failed to find property in ply file: %s\n" , PlyOrientedVertex< Real >::ReadProperties[i].name );
 					exit( 0 );
 				}
 			if( _validationFunction )
diff --git a/Src/PoissonRecon.cpp b/Src/PoissonRecon.cpp
index af995ff..dba820f 100644
--- a/Src/PoissonRecon.cpp
+++ b/Src/PoissonRecon.cpp
@@ -287,7 +287,7 @@ PlyProperty PlyColorProperties[]=
 
 bool ValidPlyColorProperties( const bool* props ){ return ( props[0] || props[3] ) && ( props[1] || props[4] ) && ( props[2] || props[5] ); }
 
-template< class Real , class Vertex , class _Vertex >
+template< class Real , class Vertex >
 int Execute( int argc , char* argv[] )
 {
 	Reset< Real >();
@@ -407,6 +407,7 @@ int Execute( int argc , char* argv[] )
 	Pointer( Real ) solution = tree.SolveSystem( *pointInfo , constraints , ShowResidual.set , Iters.value , MaxSolveDepth.value , CGDepth.value , CSSolverAccuracy.value );
 	delete pointInfo;
 	FreePointer( constraints );
+
 	DumpOutput2( comments , "# Linear system solved in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
 	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) );
 	maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage );
@@ -447,7 +448,7 @@ int Execute( int argc , char* argv[] )
 	if( Out.set )
 	{
 		t = Time() , tree.maxMemoryUsage = 0;
-		tree.template GetMCIsoSurface< Vertex , _Vertex >( kernelDensityWeights ? GetPointer( *kernelDensityWeights ) : NullPointer( Real ) , Color.set ? &colorData : NULL , solution , isoValue , mesh , true , !NonManifold.set , PolygonMesh.set );
+		tree.GetMCIsoSurface( kernelDensityWeights ? GetPointer( *kernelDensityWeights ) : NullPointer( Real ) , Color.set ? &colorData : NULL , solution , isoValue , mesh , true , !NonManifold.set , PolygonMesh.set );
 		if( PolygonMesh.set ) DumpOutput2( comments , "#         Got polygons in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
 		else                  DumpOutput2( comments , "#        Got triangles in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
 		maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage );
@@ -502,18 +503,18 @@ int main( int argc , char* argv[] )
 	cmdLineParse( argc-1 , &argv[1] , sizeof(params)/sizeof(cmdLineReadable*) , params , 1 );
 	if( Density.set )
 		if( Color.set )
-			if( Double.set ) Execute< double , PlyColorAndValueVertex< float > , _PlyColorAndValueVertex< float > >( argc , argv );
-			else             Execute< float  , PlyColorAndValueVertex< float > , _PlyColorAndValueVertex< float > >( argc , argv );
+			if( Double.set ) Execute< double , PlyColorAndValueVertex< float > >( argc , argv );
+			else             Execute< float  , PlyColorAndValueVertex< float > >( argc , argv );
 		else
-			if( Double.set ) Execute< double , PlyValueVertex< float > , PlyValueVertex< float > >( argc , argv );
-			else             Execute< float  , PlyValueVertex< float > , PlyValueVertex< float > >( argc , argv );
+			if( Double.set ) Execute< double , PlyValueVertex< float > >( argc , argv );
+			else             Execute< float  , PlyValueVertex< float > >( argc , argv );
 	else
 		if( Color.set )
-			if( Double.set ) Execute< double , PlyColorVertex< float > , _PlyColorVertex< float > >( argc , argv );
-			else             Execute< float  , PlyColorVertex< float > , _PlyColorVertex< float > >( argc , argv );
+			if( Double.set ) Execute< double , PlyColorVertex< float > >( argc , argv );
+			else             Execute< float  , PlyColorVertex< float > >( argc , argv );
 		else
-			if( Double.set ) Execute< double , PlyVertex< float > , PlyVertex< float > >( argc , argv );
-			else             Execute< float  , PlyVertex< float > , PlyVertex< float > >( argc , argv );
+			if( Double.set ) Execute< double , PlyVertex< float > >( argc , argv );
+			else             Execute< float  , PlyVertex< float > >( argc , argv );
 #ifdef _WIN32
 	if( Performance.set )
 	{
diff --git a/Src/SurfaceTrimmer.cpp b/Src/SurfaceTrimmer.cpp
index 9243f35..cc7a557 100644
--- a/Src/SurfaceTrimmer.cpp
+++ b/Src/SurfaceTrimmer.cpp
@@ -37,17 +37,14 @@ DAMAGE.
 #include "MAT.h"
 #include "MyTime.h"
 
-#define FOR_RELEASE 1
-
 cmdLineString In( "in" ) , Out( "out" );
 cmdLineInt Smooth( "smooth" , 5 );
 cmdLineFloat Trim( "trim" ) , IslandAreaRatio( "aRatio" , 0.001f );
-cmdLineFloatArray< 2 > ColorRange( "color" );
 cmdLineReadable PolygonMesh( "polygonMesh" );
 
 cmdLineReadable* params[] =
 {
-	&In , &Out , &Trim , &PolygonMesh , &ColorRange , &Smooth , &IslandAreaRatio
+	&In , &Out , &Trim , &PolygonMesh , &Smooth , &IslandAreaRatio
 };
 
 void ShowUsage( char* ex )
@@ -59,9 +56,6 @@ void ShowUsage( char* ex )
 	printf( "\t[--%s <trimming value>]\n" , Trim.name );
 	printf( "\t[--%s <relative area of islands>=%f]\n" , IslandAreaRatio.name , IslandAreaRatio.value );
 	printf( "\t[--%s]\n" , PolygonMesh.name );
-#if !FOR_RELEASE
-	printf( "\t[--%s <color range>]\n" , ColorRange.name );
-#endif // !FOR_RELEASE
 }
 
 long long EdgeKey( int key1 , int key2 )
@@ -69,33 +63,18 @@ long long EdgeKey( int key1 , int key2 )
 	if( key1<key2 ) return ( ( (long long)key1 )<<32 ) | ( (long long)key2 );
 	else            return ( ( (long long)key2 )<<32 ) | ( (long long)key1 );
 }
-template< class Real >
-PlyValueVertex< Real > InterpolateVertices( const PlyValueVertex< Real >& v1 , const PlyValueVertex< Real >& v2 , const float& value )
-{
-	if( v1.value==v2.value ) return (v1+v2)/Real(2.);
-	PlyValueVertex< Real > v;
 
-	Real dx = (v1.value-value)/(v1.value-v2.value);
-	for( int i=0 ; i<3 ; i++ ) v.point.coords[i]=v1.point.coords[i]*(1.f-dx)+v2.point.coords[i]*dx;
-	v.value=v1.value*(1.f-dx)+v2.value*dx;
-	return v;
-}
-
-template< class Real >
-void ColorVertices( const std::vector< PlyValueVertex< Real > >& inVertices , std::vector< PlyColorVertex< Real > >& outVertices , float min , float max )
+template< class Real , class Vertex >
+Vertex InterpolateVertices( const Vertex& v1 , const Vertex& v2 , Real value )
 {
-	outVertices.resize( inVertices.size() );
-	for( size_t i=0 ; i<inVertices.size() ; i++ )
-	{
-		outVertices[i].point = inVertices[i].point;
-		float temp = ( inVertices[i].value-min ) / ( max - min );
-		temp = std::max< float >( 0.f , std::min< float >( 1.f , temp ) );
-		temp *= 255;
-		outVertices[i].color[0] = outVertices[i].color[1] = outVertices[i].color[2] = (int)temp;
-	}
+	typename Vertex::Wrapper _v1(v1) , _v2(v2);
+	if( _v1.value==_v2.value ) return Vertex( (_v1+_v2)/Real(2.) );
+
+	Real dx = ( _v1.value-value ) / ( _v1.value-_v2.value );
+	return Vertex( _v1*(1.f-dx) + _v2*dx );
 }
-template< class Real >
-void SmoothValues( std::vector< PlyValueVertex< Real > >& vertices , const std::vector< std::vector< int > >& polygons )
+template< class Real , class Vertex >
+void SmoothValues( std::vector< Vertex >& vertices , const std::vector< std::vector< int > >& polygons )
 {
 	std::vector< int > count( vertices.size() );
 	std::vector< Real > sums( vertices.size() , 0 );
@@ -112,26 +91,11 @@ void SmoothValues( std::vector< PlyValueVertex< Real > >& vertices , const std::
 	}
 	for( size_t i=0 ; i<vertices.size() ; i++ ) vertices[i].value = ( sums[i] + vertices[i].value ) / ( count[i] + 1 );
 }
-template< class Real >
-void SmoothValues( std::vector< PlyValueVertex< Real > >& vertices , const std::vector< std::vector< int > >& polygons , Real min , Real max )
-{
-	std::vector< int > count( vertices.size() );
-	std::vector< Real > sums( vertices.size() , 0 );
-	for( int i=0 ; i<polygons.size() ; i++ ) for( int j=0 ; j<polygons[i].size() ; j++ )
-	{
-		int j1 = j , j2 = (j+1)%polygons[i].size();
-		int v1 = polygons[i][j1] , v2 = polygons[i][j2];
-		if( vertices[v1].value>min && vertices[v1].value<max ) count[v1]++ , sums[v1] += vertices[v2].value;
-		if( vertices[v2].value>min && vertices[v2].value<max ) count[v2]++ , sums[v2] += vertices[v1].value;
-	}
-	for( int i=0 ; i<vertices.size() ; i++ ) vertices[i].value = ( sums[i] + vertices[i].value ) / ( count[i] + 1 );
-}
-
-template< class Real >
+template< class Real , class Vertex >
 void SplitPolygon
 	(
 	const std::vector< int >& polygon ,
-	std::vector< PlyValueVertex< Real > >& vertices , 
+	std::vector< Vertex >& vertices , 
 	std::vector< std::vector< int > >* ltPolygons , std::vector< std::vector< int > >* gtPolygons ,
 	std::vector< bool >* ltFlags , std::vector< bool >* gtFlags ,
 	hash_map< long long , int >& vertexTable ,
@@ -195,8 +159,8 @@ void SplitPolygon
 		}
 	}
 }
-template< class Real >
-void Triangulate( const std::vector< PlyValueVertex< Real > >& vertices , const std::vector< std::vector< int > >& polygons , std::vector< std::vector< int > >& triangles )
+template< class Real , class Vertex >
+void Triangulate( const std::vector< Vertex >& vertices , const std::vector< std::vector< int > >& polygons , std::vector< std::vector< int > >& triangles )
 {
 	triangles.clear();
 	for( size_t i=0 ; i<polygons.size() ; i++ )
@@ -219,6 +183,22 @@ void Triangulate( const std::vector< PlyValueVertex< Real > >& vertices , const
 		}
 		else if( polygons[i].size()==3 ) triangles.push_back( polygons[i] );
 }
+template< class Real , class Vertex >
+double PolygonArea( const std::vector< Vertex >& vertices , const std::vector< int >& polygon )
+{
+	if( polygon.size()<3 ) return 0.;
+	else if( polygon.size()==3 ) return TriangleArea( vertices[polygon[0]].point , vertices[polygon[1]].point , vertices[polygon[2]].point );
+	else
+	{
+		Point3D< Real > center;
+		for( size_t i=0 ; i<polygon.size() ; i++ ) center += vertices[ polygon[i] ].point;
+		center /= Real( polygon.size() );
+		double area = 0;
+		for( size_t i=0 ; i<polygon.size() ; i++ ) area += TriangleArea( center , vertices[ polygon[i] ].point , vertices[ polygon[ (i+1)%polygon.size() ] ].point );
+		return area;
+	}
+}
+
 template< class Vertex >
 void RemoveHangingVertices( std::vector< Vertex >& vertices , std::vector< std::vector< int > >& polygons )
 {
@@ -288,139 +268,112 @@ double TriangleArea( Point3D< Real > v1 , Point3D< Real > v2 , Point3D< Real > v
 	Point3D< Real > n = CrossProduct( v2-v1 , v3-v1 );
 	return sqrt( n[0]*n[0] + n[1]*n[1] + n[2]*n[2] ) / 2.;
 }
-template< class Real >
-double PolygonArea( const std::vector< PlyValueVertex< Real > >& vertices , const std::vector< int >& polygon )
+template< class Vertex >
+int Execute( void )
 {
-	if( polygon.size()<3 ) return 0.;
-	else if( polygon.size()==3 ) return TriangleArea( vertices[polygon[0]].point , vertices[polygon[1]].point , vertices[polygon[2]].point );
-	else
+	float min , max;
+	int paramNum = sizeof(params)/sizeof(cmdLineReadable*);
+	std::vector< Vertex > vertices;
+	std::vector< std::vector< int > > polygons;
+
+	int ft , commentNum = paramNum+2;
+	char** comments;
+	PlyReadPolygons( In.value , vertices , polygons , Vertex::ReadProperties , Vertex::ReadComponents , ft , &comments , &commentNum );
+	for( int i=0 ; i<Smooth.value ; i++ ) SmoothValues< float , Vertex >( vertices , polygons );
+	min = max = vertices[0].value;
+	for( size_t i=0 ; i<vertices.size() ; i++ ) min = std::min< float >( min , vertices[i].value ) , max = std::max< float >( max , vertices[i].value );
+	printf( "Value Range: [%f,%f]\n" , min , max );
+
+
+	hash_map< long long , int > vertexTable;
+	std::vector< std::vector< int > > ltPolygons , gtPolygons;
+	std::vector< bool > ltFlags , gtFlags;
+
+	for( int i=0 ; i<paramNum+2 ; i++ ) comments[i+commentNum]=new char[1024];
+	sprintf( comments[commentNum++] , "Running Surface Trimmer (V5)" );
+	if(              In.set ) sprintf(comments[commentNum++],"\t--%s %s" , In.name , In.value );
+	if(             Out.set ) sprintf(comments[commentNum++],"\t--%s %s" , Out.name , Out.value );
+	if(            Trim.set ) sprintf(comments[commentNum++],"\t--%s %f" , Trim.name , Trim.value );
+	if(          Smooth.set ) sprintf(comments[commentNum++],"\t--%s %d" , Smooth.name , Smooth.value );
+	if( IslandAreaRatio.set ) sprintf(comments[commentNum++],"\t--%s %f" , IslandAreaRatio.name , IslandAreaRatio.value );
+	if(     PolygonMesh.set ) sprintf(comments[commentNum++],"\t--%s" , PolygonMesh.name );
+
+	double t=Time();
+	for( size_t i=0 ; i<polygons.size() ; i++ ) SplitPolygon( polygons[i] , vertices , &ltPolygons , &gtPolygons , &ltFlags , &gtFlags , vertexTable , Trim.value );
+	if( IslandAreaRatio.value>0 )
 	{
-		Point3D< Real > center;
-		for( size_t i=0 ; i<polygon.size() ; i++ ) center += vertices[ polygon[i] ].point;
-		center /= Real( polygon.size() );
-		double area = 0;
-		for( size_t i=0 ; i<polygon.size() ; i++ ) area += TriangleArea( center , vertices[ polygon[i] ].point , vertices[ polygon[ (i+1)%polygon.size() ] ].point );
-		return area;
+		std::vector< std::vector< int > > _ltPolygons , _gtPolygons;
+		std::vector< std::vector< int > > ltComponents , gtComponents;
+		SetConnectedComponents( ltPolygons , ltComponents );
+		SetConnectedComponents( gtPolygons , gtComponents );
+		std::vector< double > ltAreas( ltComponents.size() , 0. ) , gtAreas( gtComponents.size() , 0. );
+		std::vector< bool > ltComponentFlags( ltComponents.size() , false ) , gtComponentFlags( gtComponents.size() , false );
+		double area = 0.;
+		for( size_t i=0 ; i<ltComponents.size() ; i++ )
+		{
+			for( size_t j=0 ; j<ltComponents[i].size() ; j++ )
+			{
+				ltAreas[i] += PolygonArea< float , Vertex >( vertices , ltPolygons[ ltComponents[i][j] ] );
+				ltComponentFlags[i] = ( ltComponentFlags[i] || ltFlags[ ltComponents[i][j] ] );
+			}
+			area += ltAreas[i];
+		}
+		for( size_t i=0 ; i<gtComponents.size() ; i++ )
+		{
+			for( size_t j=0 ; j<gtComponents[i].size() ; j++ )
+			{
+				gtAreas[i] += PolygonArea< float , Vertex >( vertices , gtPolygons[ gtComponents[i][j] ] );
+				gtComponentFlags[i] = ( gtComponentFlags[i] || gtFlags[ gtComponents[i][j] ] );
+			}
+			area += gtAreas[i];
+		}
+		for( size_t i=0 ; i<ltComponents.size() ; i++ )
+		{
+			if( ltAreas[i]<area*IslandAreaRatio.value && ltComponentFlags[i] ) for( size_t j=0 ; j<ltComponents[i].size() ; j++ ) _gtPolygons.push_back( ltPolygons[ ltComponents[i][j] ] );
+			else                                                               for( size_t j=0 ; j<ltComponents[i].size() ; j++ ) _ltPolygons.push_back( ltPolygons[ ltComponents[i][j] ] );
+		}
+		for( size_t i=0 ; i<gtComponents.size() ; i++ )
+		{
+			if( gtAreas[i]<area*IslandAreaRatio.value && gtComponentFlags[i] ) for( size_t j=0 ; j<gtComponents[i].size() ; j++ ) _ltPolygons.push_back( gtPolygons[ gtComponents[i][j] ] );
+			else                                                               for( size_t j=0 ; j<gtComponents[i].size() ; j++ ) _gtPolygons.push_back( gtPolygons[ gtComponents[i][j] ] );
+		}
+		ltPolygons = _ltPolygons , gtPolygons = _gtPolygons;
 	}
+	if( !PolygonMesh.set )
+	{
+		{
+			std::vector< std::vector< int > > polys = ltPolygons;
+			Triangulate< float , Vertex >( vertices , ltPolygons , polys ) , ltPolygons = polys;
+		}
+		{
+			std::vector< std::vector< int > > polys = gtPolygons;
+			Triangulate< float  , Vertex >( vertices , gtPolygons , polys ) , gtPolygons = polys;
+		}
+	}
+
+	RemoveHangingVertices( vertices , gtPolygons );
+	sprintf( comments[commentNum++] , "#Trimmed In: %9.1f (s)" , Time()-t );
+	if( Out.set ) PlyWritePolygons( Out.value , vertices , gtPolygons , Vertex::WriteProperties , Vertex::WriteComponents , ft , comments , commentNum );
+
+	return EXIT_SUCCESS;
 }
 int main( int argc , char* argv[] )
 {
 	int paramNum = sizeof(params)/sizeof(cmdLineReadable*);
 	cmdLineParse( argc-1 , &argv[1] , paramNum , params , 0 );
 
-#if FOR_RELEASE
 	if( !In.set || !Trim.set )
 	{
 		ShowUsage( argv[0] );
 		return EXIT_FAILURE;
 	}
-#else // !FOR_RELEASE
-	if( !In.set )
-	{
-		ShowUsage( argv[0] );
-		return EXIT_FAILURE;
-	}
-#endif // FOR_RELEASE
-	float min , max;
-	std::vector< PlyValueVertex< float > > vertices;
-	std::vector< std::vector< int > > polygons;
+	bool readFlags[ PlyColorAndValueVertex< float >::ReadComponents ];
+	if( !PlyReadHeader( In.value , PlyColorAndValueVertex< float >::ReadProperties , PlyColorAndValueVertex< float >::ReadComponents , readFlags ) ) fprintf( stderr , "[ERROR] Failed to read ply header: %s\n" , In.value ) , exit( 0 );
 
-	int ft , commentNum = paramNum+2;
-	char** comments;
-	bool readFlags[ PlyValueVertex< float >::Components ];
-	PlyReadPolygons( In.value , vertices , polygons , PlyValueVertex< float >::Properties , PlyValueVertex< float >::Components , ft , &comments , &commentNum , readFlags );
-	if( !readFlags[3] ){ fprintf( stderr , "[ERROR] vertices do not have value flag\n" ) ; return EXIT_FAILURE; }
-#if 0
-	if( Trim.set ) for( int i=0 ; i<Smooth.value ; i++ ) SmoothValues( vertices , polygons , Trim.value-0.5f , Trim.value+0.5f );
-	else           for( int i=0 ; i<Smooth.value ; i++ ) SmoothValues( vertices , polygons );
-#else
-	for( int i=0 ; i<Smooth.value ; i++ ) SmoothValues( vertices , polygons );
-#endif 
-	min = max = vertices[0].value;
-	for( size_t i=0 ; i<vertices.size() ; i++ ) min = std::min< float >( min , vertices[i].value ) , max = std::max< float >( max , vertices[i].value );
-	printf( "Value Range: [%f,%f]\n" , min , max );
+	bool hasValue = readFlags[3];
+	bool hasColor = ( readFlags[4] || readFlags[7] ) && ( readFlags[5] || readFlags[8] ) && ( readFlags[6] || readFlags[9] );
 
-
-	if( Trim.set )
-	{
-		hash_map< long long , int > vertexTable;
-		std::vector< std::vector< int > > ltPolygons , gtPolygons;
-		std::vector< bool > ltFlags , gtFlags;
-
-		for( int i=0 ; i<paramNum+2 ; i++ ) comments[i+commentNum]=new char[1024];
-		sprintf( comments[commentNum++] , "Running Surface Trimmer (V5)" );
-		if(              In.set ) sprintf(comments[commentNum++],"\t--%s %s" , In.name , In.value );
-		if(             Out.set ) sprintf(comments[commentNum++],"\t--%s %s" , Out.name , Out.value );
-		if(            Trim.set ) sprintf(comments[commentNum++],"\t--%s %f" , Trim.name , Trim.value );
-		if(          Smooth.set ) sprintf(comments[commentNum++],"\t--%s %d" , Smooth.name , Smooth.value );
-		if( IslandAreaRatio.set ) sprintf(comments[commentNum++],"\t--%s %f" , IslandAreaRatio.name , IslandAreaRatio.value );
-		if(     PolygonMesh.set ) sprintf(comments[commentNum++],"\t--%s" , PolygonMesh.name );
-
-		double t=Time();
-		for( size_t i=0 ; i<polygons.size() ; i++ ) SplitPolygon( polygons[i] , vertices , &ltPolygons , &gtPolygons , &ltFlags , &gtFlags , vertexTable , Trim.value );
-		if( IslandAreaRatio.value>0 )
-		{
-			std::vector< std::vector< int > > _ltPolygons , _gtPolygons;
-			std::vector< std::vector< int > > ltComponents , gtComponents;
-			SetConnectedComponents( ltPolygons , ltComponents );
-			SetConnectedComponents( gtPolygons , gtComponents );
-			std::vector< double > ltAreas( ltComponents.size() , 0. ) , gtAreas( gtComponents.size() , 0. );
-			std::vector< bool > ltComponentFlags( ltComponents.size() , false ) , gtComponentFlags( gtComponents.size() , false );
-			double area = 0.;
-			for( size_t i=0 ; i<ltComponents.size() ; i++ )
-			{
-				for( size_t j=0 ; j<ltComponents[i].size() ; j++ )
-				{
-					ltAreas[i] += PolygonArea( vertices , ltPolygons[ ltComponents[i][j] ] );
-					ltComponentFlags[i] = ( ltComponentFlags[i] || ltFlags[ ltComponents[i][j] ] );
-				}
-				area += ltAreas[i];
-			}
-			for( size_t i=0 ; i<gtComponents.size() ; i++ )
-			{
-				for( size_t j=0 ; j<gtComponents[i].size() ; j++ )
-				{
-					gtAreas[i] += PolygonArea( vertices , gtPolygons[ gtComponents[i][j] ] );
-					gtComponentFlags[i] = ( gtComponentFlags[i] || gtFlags[ gtComponents[i][j] ] );
-				}
-				area += gtAreas[i];
-			}
-			for( size_t i=0 ; i<ltComponents.size() ; i++ )
-			{
-				if( ltAreas[i]<area*IslandAreaRatio.value && ltComponentFlags[i] ) for( size_t j=0 ; j<ltComponents[i].size() ; j++ ) _gtPolygons.push_back( ltPolygons[ ltComponents[i][j] ] );
-				else                                                               for( size_t j=0 ; j<ltComponents[i].size() ; j++ ) _ltPolygons.push_back( ltPolygons[ ltComponents[i][j] ] );
-			}
-			for( size_t i=0 ; i<gtComponents.size() ; i++ )
-			{
-				if( gtAreas[i]<area*IslandAreaRatio.value && gtComponentFlags[i] ) for( size_t j=0 ; j<gtComponents[i].size() ; j++ ) _ltPolygons.push_back( gtPolygons[ gtComponents[i][j] ] );
-				else                                                               for( size_t j=0 ; j<gtComponents[i].size() ; j++ ) _gtPolygons.push_back( gtPolygons[ gtComponents[i][j] ] );
-			}
-			ltPolygons = _ltPolygons , gtPolygons = _gtPolygons;
-		}
-		if( !PolygonMesh.set )
-		{
-			{
-				std::vector< std::vector< int > > polys = ltPolygons;
-				Triangulate( vertices , ltPolygons , polys ) , ltPolygons = polys;
-			}
-			{
-				std::vector< std::vector< int > > polys = gtPolygons;
-				Triangulate( vertices , gtPolygons , polys ) , gtPolygons = polys;
-			}
-		}
-
-		RemoveHangingVertices( vertices , gtPolygons );
-		sprintf( comments[commentNum++] , "#Trimmed In: %9.1f (s)" , Time()-t );
-		if( Out.set ) PlyWritePolygons( Out.value , vertices , gtPolygons , PlyValueVertex< float >::Properties , PlyValueVertex< float >::Components , ft , comments , commentNum );
-	}
-	else
-	{
-		if( ColorRange.set ) min = ColorRange.values[0] , max = ColorRange.values[1];
-		std::vector< PlyColorVertex< float > > outVertices;
-		ColorVertices( vertices , outVertices , min , max );
-		if( Out.set ) PlyWritePolygons( Out.value , outVertices , polygons , PlyColorVertex< float >::Properties , PlyColorVertex< float >::Components , ft , comments , commentNum );
-	}
-
-	return EXIT_SUCCESS;
+	if( !hasValue ) fprintf( stderr , "[ERROR] Ply file does not contain values\n" ) , exit( 0 );
+	if( hasColor ) return Execute< PlyColorAndValueVertex< float > >();
+	else           return Execute< PlyValueVertex< float > >();
 }
-