diff --git a/.gitignore b/.gitignore index 4bccb39..bcc3d09 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,4 @@ *.sdf -*.suo -*.filters /Obj/ /Debug/ /Bin/ diff --git a/PoissonRecon.v12.suo b/PoissonRecon.v12.suo new file mode 100644 index 0000000..5f8f226 Binary files /dev/null and b/PoissonRecon.v12.suo differ diff --git a/PoissonRecon.vcxproj.filters b/PoissonRecon.vcxproj.filters new file mode 100644 index 0000000..852bbac --- /dev/null +++ b/PoissonRecon.vcxproj.filters @@ -0,0 +1,143 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + inc;inl + + + + + Source Files + + + Source Files + + + Source Files + + + Source Files + + + Source Files + + + Source Files + + + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + Include Files + + + \ No newline at end of file diff --git a/Src/Array.h b/Src/Array.h index 5added0..5138513 100644 --- a/Src/Array.h +++ b/Src/Array.h @@ -62,6 +62,7 @@ void aligned_free( void* mem ) { free( ( ( void** )mem )[-1] ); } #include "Array.inl" #define Pointer( ... ) Array< __VA_ARGS__ > #define ConstPointer( ... ) ConstArray< __VA_ARGS__ > +#define NullPointer( ... ) Array< __VA_ARGS__ >() template< class C > void FreePointer( Array< C >& a ){ a.Free( ); } template< class C > void AlignedFreePointer( Array< C >& a ){ a.Free( ); } template< class C > void VFreePointer( Array< C >& a ){ a.Free( ); } @@ -72,7 +73,7 @@ template< class C > Array< C > AllocPointer( size_t size template< class C > Array< C > AlignedAllocPointer( size_t size , size_t alignment , const char* name=NULL ){ return Array< C >::AlignedAlloc( size , alignment , false , name ); } template< class C > Array< C > ReAllocPointer( Array< C >& a , size_t size , const char* name=NULL ){ return Array< C >::ReAlloc ( a , size , false , name ); } -template< class C > Array< C > NullPointer( void ){ return Array< C >( ); } +//template< class C > Array< C > NullPointer( void ){ return Array< C >( ); } template< class C > C* PointerAddress( Array< C >& a ) { return a.pointer(); } template< class C > const C* PointerAddress( ConstArray< C >& a ) { return a.pointer(); } @@ -84,6 +85,7 @@ template< class C > ConstArray< C > GetPointer( const std::vector< C >& v ){ ret #else // !ARRAY_DEBUG #define Pointer( ... ) __VA_ARGS__* #define ConstPointer( ... ) const __VA_ARGS__* +#define NullPointer( ... ) NULL #define FreePointer( ... ) { if( __VA_ARGS__ ) free( __VA_ARGS__ ) , __VA_ARGS__ = NULL; } #define AlignedFreePointer( ... ) { if( __VA_ARGS__ ) aligned_free( __VA_ARGS__ ) , __VA_ARGS__ = NULL; } @@ -94,7 +96,7 @@ template< class C > C* AllocPointer( size_t size , template< class C > C* AlignedAllocPointer( size_t size , size_t alignment , const char* name=NULL ){ return (C*)aligned_malloc( sizeof(C) * size , alignment ); } template< class C > C* ReAllocPointer( C* c , size_t size , const char* name=NULL ){ return (C*) realloc( c , sizeof(C) * size ); } -template< class C > C* NullPointer( void ){ return NULL; } +//template< class C > C* NullPointer( void ){ return NULL; } template< class C > C* PointerAddress( C* c ){ return c; } template< class C > const C* PointerAddress( const C* c ){ return c; } diff --git a/Src/BSplineData.inl b/Src/BSplineData.inl index ef9cb7b..062858b 100644 --- a/Src/BSplineData.inl +++ b/Src/BSplineData.inl @@ -183,7 +183,7 @@ double BSplineData< Degree >::CornerEvaluator< Radius >::value( int depth , int template< int Degree > void BSplineData< Degree >::set( int maxDepth , int boundaryType ) { - _boundaryType = boundaryType; + _boundaryType = boundaryType>0 ? 1 : ( boundaryType<0 ? -1 : 0 ); depth = maxDepth; // [Warning] This assumes that the functions spacing is dual @@ -378,9 +378,9 @@ template< int Degree > template< class Real > BSplineData< Degree >::DotTables< Real >::DotTables( void ) { - vvDotTable = NullPointer< Real >(); - dvDotTable = NullPointer< Real >(); - ddDotTable = NullPointer< Real >(); + vvDotTable = NullPointer( Real ); + dvDotTable = NullPointer( Real ); + ddDotTable = NullPointer( Real ); } template< int Degree > template< class Real > @@ -533,8 +533,8 @@ template< int Degree > template< class Real > BSplineData< Degree >::ValueTables< Real >::ValueTables( void ) { - valueTable = NullPointer< Real >(); - dValueTable = NullPointer< Real >(); + valueTable = NullPointer( Real ); + dValueTable = NullPointer( Real ); } template< int Degree > template< class Real > diff --git a/Src/Geometry.h b/Src/Geometry.h index 867e993..0d73eda 100644 --- a/Src/Geometry.h +++ b/Src/Geometry.h @@ -42,6 +42,7 @@ struct Point3D { Real coords[3]; Point3D( void ) { coords[0] = coords[1] = coords[2] = Real(0); } + Point3D( Real v ) { coords[0] = coords[1] = coords[2] = v; } template< class _Real > Point3D( _Real v0 , _Real v1 , _Real v2 ){ coords[0] = Real(v0) , coords[1] = Real(v1) , coords[2] = Real(v2); } template< class _Real > Point3D( const Point3D< _Real >& p ){ coords[0] = Real( p[0] ) , coords[1] = Real( p[1] ) , coords[2] = Real( p[2] ); } inline Real& operator[] ( int i ) { return coords[i]; } diff --git a/Src/InterpolateColors.V1.cpp b/Src/InterpolateColors.V1.cpp new file mode 100644 index 0000000..86084fa --- /dev/null +++ b/Src/InterpolateColors.V1.cpp @@ -0,0 +1,540 @@ + +/* +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 +#include +#include +#include +#ifdef _WIN32 +#include +#include +#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 +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 \n" , InPoints.name ); + printf( "\t --%s \n" , InMesh.name ); + + printf( "\t[--%s ]\n" , Out.name ); + printf( "\t[--%s ]\n" , VoxelGrid.name ); + + printf( "\t[--%s =%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 =%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 =<%s>]\n" , VoxelDepth.name , Depth.name ); + + printf( "\t[--%s =%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 =%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 ]\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 =%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 =%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 ; iset ) + { + 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 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 ; cnodeData.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::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 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 ; cnodeData.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 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 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 ; i0 ) + { + 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 new file mode 100644 index 0000000..353c3f0 --- /dev/null +++ b/Src/InterpolateColors.cpp @@ -0,0 +1,407 @@ + +/* +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 +#include +#include +#include +#ifdef _WIN32 +#include +#include +#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 +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 \n" , InPoints.name ); + printf( "\t --%s \n" , InMesh.name ); + + printf( "\t[--%s ]\n" , Out.name ); + printf( "\t[--%s ]\n" , VoxelGrid.name ); + + printf( "\t[--%s =%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 =%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 =<%s>]\n" , VoxelDepth.name , Depth.name ); + + printf( "\t[--%s =%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 =%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 ]\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 =%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 =%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 ; iset ) + { + 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 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 ; i0 ) + { + 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 8a90e42..77d7eb5 100644 --- a/Src/MultiGridOctreeData.IsoSurface.inl +++ b/Src/MultiGridOctreeData.IsoSurface.inl @@ -38,10 +38,10 @@ template< class Vertex > Octree< Real >::SliceValues< Vertex >::SliceValues( void ) { _oldCCount = _oldECount = _oldFCount = _oldNCount = 0; - cornerValues = NullPointer< Real >() ; cornerNormals = NullPointer< Point3D< Real > >() ; cornerSet = NullPointer< char >(); - edgeKeys = NullPointer< long long >() ; edgeSet = NullPointer< char >(); - faceEdges = NullPointer< FaceEdges >() ; faceSet = NullPointer< char >(); - mcIndices = NullPointer< char >(); + cornerValues = NullPointer( Real ) ; cornerNormals = NullPointer( Point3D< Real > ) ; cornerSet = NullPointer( char ); + edgeKeys = NullPointer( long long ) ; edgeSet = NullPointer( char ); + faceEdges = NullPointer( FaceEdges ) ; faceSet = NullPointer( char ); + mcIndices = NullPointer( char ); } template< class Real > template< class Vertex > @@ -100,8 +100,8 @@ template< class Vertex > Octree< Real >::XSliceValues< Vertex >::XSliceValues( void ) { _oldECount = _oldFCount = 0; - edgeKeys = NullPointer< long long >() ; edgeSet = NullPointer< char >(); - faceEdges = NullPointer< FaceEdges >() ; faceSet = NullPointer< char >(); + edgeKeys = NullPointer( long long ) ; edgeSet = NullPointer( char ); + faceEdges = NullPointer( FaceEdges ) ; faceSet = NullPointer( char ); } template< class Real > template< class Vertex > @@ -136,8 +136,8 @@ void Octree< Real >::XSliceValues< Vertex >::reset( void ) } template< class Real > -template< class Vertex > -void Octree< Real >::GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights , ConstPointer( Real ) solution , Real isoValue , CoredMeshData< Vertex >& mesh , bool nonLinearFit , bool addBarycenter , bool polygonMesh ) +template< class Vertex , 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; _fData.setCornerEvaluator( evaluator , 0 , 0 , _boundaryType==0 ); @@ -179,7 +179,7 @@ void Octree< Real >::GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights // Copy edges from finer if( d::GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights // Set the slice values/vertices SetSliceIsoCorners( solution , GetPointer( coarseSolution ) , isoValue , d , o , slabValues , evaluator , vStencils[d].stencil , vStencils[d].stencils , nStencils[d].stencil , nStencils[d].stencils , threads ); - SetSliceIsoVertices( kernelDensityWeights , isoValue , d , o , vertexOffset , mesh , slabValues , threads ); + SetSliceIsoVertices( kernelDensityWeights , colorData , isoValue , d , o , vertexOffset , mesh , slabValues , threads ); SetSliceIsoEdges( d , o , slabValues , threads ); // Set the cross-slice edges - SetXSliceIsoVertices( kernelDensityWeights , isoValue , d , o-1 , vertexOffset , mesh , slabValues , threads ); + SetXSliceIsoVertices( kernelDensityWeights , colorData , isoValue , d , o-1 , vertexOffset , mesh , slabValues , threads ); SetXSliceIsoEdges( d , o-1 , slabValues , threads ); // Add the triangles - SetIsoSurface( d , o-1 , slabValues[d].sliceValues(o-1) , slabValues[d].sliceValues(o) , slabValues[d].xSliceValues(o-1) , mesh , polygonMesh , addBarycenter , vertexOffset , threads ); + 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 ); if( o&1 ) break; } @@ -353,14 +353,14 @@ void Octree< Real >::SetSliceIsoCorners( ConstPointer( Real ) solution , ConstPo template< class Real > template< class Vertex > -void Octree< Real >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , Real isoValue , int depth , int slice , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +void Octree< Real >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , int depth , int slice , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) { - if( slice>0 ) SetSliceIsoVertices( kernelDensityWeight , isoValue , depth , slice , 1 , vOffset , mesh , slabValues , threads ); - if( slice<(1<0 ) SetSliceIsoVertices( kernelDensityWeight , colorData , isoValue , depth , slice , 1 , vOffset , mesh , slabValues , threads ); + if( slice<(1< template< class Vertex > -void Octree< Real >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , Real isoValue , int depth , int slice , int z , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +void Octree< Real >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , int depth , int slice , int z , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) { typename Octree< Real >::template SliceValues< Vertex >& sValues = slabValues[depth].sliceValues( slice ); std::vector< typename TreeOctNode::ConstNeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) ); @@ -387,7 +387,7 @@ void Octree< Real >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeig int o , y; Square::FactorEdgeIndex( e , o , y ); long long key = VertexData::EdgeIndex( leaf , Cube::EdgeIndex( o , y , z ) , _sNodes.maxDepth ); - GetIsoVertex( kernelDensityWeight , isoValue , neighborKey , leaf , e , z , sValues , vertex ); + GetIsoVertex( kernelDensityWeight , colorData , isoValue , neighborKey , leaf , e , z , sValues , vertex ); vertex.point = vertex.point * _scale + _center; bool stillOwner = false; std::pair< int , Vertex > hashed_vertex; @@ -449,7 +449,7 @@ void Octree< Real >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeig } template< class Real > template< class Vertex > -void Octree< Real >::SetXSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , Real isoValue , int depth , int slab , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +void Octree< Real >::SetXSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , int depth , int slab , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) { typename Octree< Real >::template SliceValues< Vertex >& bValues = slabValues[depth].sliceValues ( slab ); typename Octree< Real >::template SliceValues< Vertex >& fValues = slabValues[depth].sliceValues ( slab+1 ); @@ -480,7 +480,7 @@ void Octree< Real >::SetXSliceIsoVertices( ConstPointer( Real ) kernelDensityWei { Vertex vertex; long long key = VertexData::EdgeIndex( leaf , e , _sNodes.maxDepth ); - GetIsoVertex( kernelDensityWeight , isoValue , neighborKey , leaf , c , bValues , fValues , vertex ); + GetIsoVertex( kernelDensityWeight , colorData , isoValue , neighborKey , leaf , c , bValues , fValues , vertex ); vertex.point = vertex.point * _scale + _center; bool stillOwner = false; std::pair< int , Vertex > hashed_vertex; @@ -791,7 +791,7 @@ void Octree< Real >::SetXSliceIsoEdges( int depth , int slab , std::vector< Slab } } template< class Real > -template< class Vertex > +template< class Vertex , 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,19 +902,26 @@ 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( mesh , polygon , polygonMesh , addBarycenter , vOffset ); + AddIsoPolygons< Vertex , _Vertex >( mesh , polygon , polygonMesh , addBarycenter , vOffset ); } } } } } -template< class Real > void SetIsoVertexValue( PlyVertex< float >& vertex , Real value ){ ; } -template< class Real > void SetIsoVertexValue( PlyValueVertex< float >& vertex , Real value ){ vertex.value = float(value); } -template< class Real > void SetIsoVertexValue( PlyVertex< double >& vertex , Real value ){ ; } -template< class Real > void SetIsoVertexValue( PlyValueVertex< double >& vertex , Real value ){ vertex.value = double(value); } +template< class Real > void SetColor( Point3D< Real >& color , unsigned char c[3] ){ for( int i=0 ; i<3 ; i++ ) c[i] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[i]+0.5 ) ) ); } + +template< class Real > void SetIsoVertex( PlyVertex< float >& vertex , Point3D< Real > color , Real value ){ ; } +template< class Real > void SetIsoVertex( PlyColorVertex< float >& vertex , Point3D< Real > color , Real value ){ SetColor( color , vertex.color ); } +template< class Real > void SetIsoVertex( PlyValueVertex< float >& vertex , Point3D< Real > color , Real value ){ vertex.value = float(value); } +template< class Real > void SetIsoVertex( PlyColorAndValueVertex< float >& vertex , Point3D< Real > color , Real value ){ SetColor( color , vertex.color ) , vertex.value = float(value); } +template< class Real > void SetIsoVertex( PlyVertex< double >& vertex , Point3D< Real > color , Real value ){ ; } +template< class Real > void SetIsoVertex( PlyColorVertex< double >& vertex , Point3D< Real > color , Real value ){ SetColor( color , vertex.color ); } +template< class Real > void SetIsoVertex( PlyValueVertex< double >& vertex , Point3D< Real > color , Real value ){ vertex.value = double(value); } +template< class Real > void SetIsoVertex( PlyColorAndValueVertex< double >& vertex , Point3D< Real > color , Real value ){ SetColor( color , vertex.color ) , vertex.value = double(value); } + template< class Real > template< class Vertex > -bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int edgeIndex , int z , const SliceValues< Vertex >& sValues , Vertex& vertex ) +bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int edgeIndex , int z , const SliceValues< Vertex >& sValues , Vertex& vertex ) { Point3D< Real > position; int c0 , c1; @@ -948,11 +955,6 @@ bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , R } double dx0 = n0[o] , dx1 = n1[o]; - { - double dAverage = ( dx0 + dx1 ) / 2.; - dx0 = dx0 * (1.-_normalSmooth) + dAverage * _normalSmooth; - dx1 = dx1 * (1.-_normalSmooth) + dAverage * _normalSmooth; - } // The scaling will turn the Hermite Spline into a quadratic double scl = (x1-x0) / ( (dx1+dx0 ) / 2 ); @@ -979,19 +981,22 @@ bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , R } position[o] = Real( center - width/2 + width*averageRoot ); vertex.point = position; + Point3D< Real > color; + Real depth(0); if( kernelDensityWeights ) { - Real depth , weight; + Real weight; const TreeOctNode* temp = node; while( temp->depth()>_splatDepth ) temp=temp->parent; - GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey3 , _samplesPerNode , depth , weight ); - SetIsoVertexValue( vertex , depth ); + GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey3 , depth , weight ); } + if( colorData ) color = Point3D< Real >( _Evaluate( *colorData , position , neighborKey3 ) ); + SetIsoVertex( vertex , color , depth ); return true; } template< class Real > template< class Vertex > -bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int cornerIndex , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , Vertex& vertex ) +bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int cornerIndex , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , Vertex& vertex ) { Point3D< Real > position; @@ -1015,11 +1020,6 @@ bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , R position[1] = c[1] - width/2 + width*y; double dx0 = n0[2] , dx1 = n1[2]; - { - double dAverage = ( dx0 + dx1 ) / 2.; - dx0 = dx0 * (1.-_normalSmooth) + dAverage * _normalSmooth; - dx1 = dx1 * (1.-_normalSmooth) + dAverage * _normalSmooth; - } // The scaling will turn the Hermite Spline into a quadratic double scl = (x1-x0) / ( (dx1+dx0 ) / 2 ); @@ -1046,19 +1046,22 @@ bool Octree< Real >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , R } position[2] = Real( center - width/2 + width*averageRoot ); vertex.point = position; + Point3D< Real > color; + Real depth(0); if( kernelDensityWeights ) { - Real depth , weight; + Real weight; const TreeOctNode* temp = node; while( temp->depth()>_splatDepth ) temp=temp->parent; - GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey3 , _samplesPerNode , depth , weight ); - SetIsoVertexValue( vertex , depth ); + GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey3 , depth , weight ); } + if( colorData ) color = Point3D< Real >( _Evaluate( *colorData , position , neighborKey3 ) ); + SetIsoVertex( vertex , color , depth ); return true; } template< class Real > -template< class Vertex > +template< class Vertex , class _Vertex> int Octree< Real >::AddIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector< std::pair< int , Vertex > >& polygon , bool polygonMesh , bool addBarycenter , int& vOffset ) { if( polygonMesh ) @@ -1084,9 +1087,11 @@ int Octree< Real >::AddIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector< if( isCoplanar ) { Vertex c; - c *= 0; - for( int i=0 ; i<(int)polygon.size() ; i++ ) c += polygon[i].second; - c /= Real( polygon.size() ); + _Vertex _c; + _c *= 0; + for( int i=0 ; i<(int)polygon.size() ; i++ ) _c += _Vertex( polygon[i].second ); + _c /= Real( polygon.size() ); + c = Vertex( _c ); int cIdx; #pragma omp critical (add_barycenter_point_access) { diff --git a/Src/MultiGridOctreeData.SortedTreeNodes.inl b/Src/MultiGridOctreeData.SortedTreeNodes.inl index e02fdfd..0fe89e1 100644 --- a/Src/MultiGridOctreeData.SortedTreeNodes.inl +++ b/Src/MultiGridOctreeData.SortedTreeNodes.inl @@ -32,9 +32,9 @@ DAMAGE. SortedTreeNodes::SortedTreeNodes( void ) { nodeCount = NULL; - treeNodes = NullPointer< TreeOctNode* >(); + treeNodes = NullPointer( TreeOctNode* ); maxDepth = 0; - sliceOffsets = NullPointer< Pointer( int ) >(); + sliceOffsets = NullPointer( Pointer( int ) ); } SortedTreeNodes::~SortedTreeNodes( void ) { diff --git a/Src/MultiGridOctreeData.h b/Src/MultiGridOctreeData.h index 762de97..8862a62 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 @@ -195,29 +196,39 @@ template< class Real > class Octree { typedef OctNode< TreeNodeData > TreeOctNode; - struct _PointData +public: + template< class V > + struct ProjectiveData + { + V v; + Real w; + ProjectiveData( V vv=V(0) , Real ww=Real(0) ) : v(vv) , w(ww) { } + operator V (){ return w!=0 ? v/w : v*w; } + ProjectiveData& operator += ( const ProjectiveData& p ){ v += p.v , w += p.w ; return *this; } + ProjectiveData& operator -= ( const ProjectiveData& p ){ v -= p.v , w -= p.w ; return *this; } + ProjectiveData& operator *= ( Real s ){ v *= s , w *= s ; return *this; } + ProjectiveData& operator /= ( Real s ){ v /= s , w /= s ; return *this; } + ProjectiveData operator + ( const ProjectiveData& p ) const { return ProjectiveData( v+p.v , w+p.w ); } + ProjectiveData operator - ( const ProjectiveData& p ) const { return ProjectiveData( v-p.v , w-p.w ); } + ProjectiveData operator * ( Real s ) const { return ProjectiveData( v*s , w*s ); } + ProjectiveData operator / ( Real s ) const { return ProjectiveData( v/s , w/s ); } + }; + struct PointData { Point3D< Real > position; Real weightedCoarserValue; Real weight; - _PointData( Point3D< Real > p=Point3D< Real >() , Real w=0 ) { position = p , weight = w , weightedCoarserValue = Real(0); } + PointData( Point3D< Real > p=Point3D< Real >() , Real w=0 ) { position = p , weight = w , weightedCoarserValue = Real(0); } }; -public: - struct NormalInfo + template< class Data > + struct SparseNodeData { - std::vector< int > normalIndices; - std::vector< Point3D< Real > > normals; - int normalIndex( const TreeOctNode* node ) const { return node->nodeData.nodeIndex>=normalIndices.size() ? -1 : normalIndices[ node->nodeData.nodeIndex ]; } - }; - struct PointInfo - { - std::vector< int > pointIndices; - std::vector< _PointData > points; - int pointIndex( const TreeOctNode* node ) const { return node->nodeData.nodeIndex>=pointIndices.size() ? -1 : pointIndices[ node->nodeData.nodeIndex ]; } + std::vector< int > indices; + std::vector< Data > data; + int index( const TreeOctNode* node ) const { return node->nodeData.nodeIndex>=indices.size() ? -1 : indices[ node->nodeData.nodeIndex ]; } }; protected: SortedTreeNodes _sNodes; - Real _samplesPerNode; int _splatDepth; int _minDepth; int _fullDepth; @@ -226,7 +237,6 @@ protected: Real _scale; Point3D< Real > _center; std::vector< int > _pointCount; - Real _normalSmooth; BSplineData< 2 > _fData; bool _InBounds( Point3D< Real > ) const; @@ -255,11 +265,11 @@ protected: }; void _setMultiColorIndices( int start , int end , std::vector< std::vector< int > >& indices ) const; - int _SolveSystemGS( PointInfo& pointInfo , int depth , const typename BSplineData< 2 >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual=false , double* bNorm2=NULL , double* inRNorm2=NULL , double* outRNorm2=NULL , bool forceSilent=false ); - int _SolveSystemCG( PointInfo& pointInfo , int depth , const typename BSplineData< 2 >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual=false , double* bNorm2=NULL , double* inRNorm2=NULL , double* outRNorm2=NULL , double accuracy=0 ); + int _SolveSystemGS( SparseNodeData< PointData >& pointInfo , int depth , const typename BSplineData< 2 >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual=false , double* bNorm2=NULL , double* inRNorm2=NULL , double* outRNorm2=NULL , bool forceSilent=false ); + int _SolveSystemCG( SparseNodeData< PointData >& pointInfo , int depth , const typename BSplineData< 2 >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual=false , double* bNorm2=NULL , double* inRNorm2=NULL , double* outRNorm2=NULL , double accuracy=0 ); int GetMatrixRowSize( const typename TreeOctNode::Neighbors5& neighbors5 , bool symmetric ) const; - int SetMatrixRow( const PointInfo& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , Pointer( MatrixEntry< Real > ) row , int offset , const typename BSplineData< 2 >::Integrator& integrator , const Stencil< double , 5 >& stencil , bool symmetric ) const; + int SetMatrixRow( const SparseNodeData< PointData >& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , Pointer( MatrixEntry< Real > ) row , int offset , const typename BSplineData< 2 >::Integrator& integrator , const Stencil< double , 5 >& stencil , bool symmetric ) const; void SetDivergenceStencil ( int depth , const typename BSplineData< 2 >::Integrator& integrator , Stencil< Point3D< double > , 5 >& stencil , bool scatter ) const; void SetDivergenceStencils( int depth , const typename BSplineData< 2 >::Integrator& integrator , Stencil< Point3D< double > , 5 > stencil[2][2][2] , bool scatter ) const; @@ -276,34 +286,45 @@ protected: static void UpdateCoarserSupportBounds( const TreeOctNode* node , int& startX , int& endX , int& startY , int& endY , int& startZ , int& endZ ); - void UpdateConstraintsFromCoarser( const PointInfo& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , const typename TreeOctNode::Neighbors5& pNeighbors5 , TreeOctNode* node , Pointer( Real ) constraints , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::Integrator& integrator , const Stencil< double , 5 >& stencil ) const; + void UpdateConstraintsFromCoarser( const SparseNodeData< PointData >& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , const typename TreeOctNode::Neighbors5& pNeighbors5 , TreeOctNode* node , Pointer( Real ) constraints , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::Integrator& integrator , const Stencil< double , 5 >& stencil ) const; // Updates the constraints @(depth-1) based on the solution coefficients @(depth) void UpdateConstraintsFromFiner( const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) fineSolution , Pointer( Real ) coarseConstraints ) const; // Evaluate the points @(depth) using coefficients @(depth-1) - void SetPointValuesFromCoarser( PointInfo& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) coarseCoefficients ); + void SetPointValuesFromCoarser( SparseNodeData< PointData >& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) coarseCoefficients ); // Evalutes the solution @(depth) at the points @(depth-1) and updates the met constraints @(depth-1) - void SetPointConstraintsFromFiner( const PointInfo& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) finerCoefficients , Pointer( Real ) metConstraints) const; - Real _WeightedCoarserFunctionValue( const _PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) coarseCoefficients ) const; - Real _WeightedFinerFunctionValue ( const _PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) finerCoefficients ) const; + void SetPointConstraintsFromFiner ( const SparseNodeData< PointData >& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) finerCoefficients , Pointer( Real ) metConstraints ) const; + Real _WeightedCoarserFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) coarseCoefficients ) const; + Real _WeightedFinerFunctionValue ( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) finerCoefficients ) const; // Down samples constraints @(depth) to constraints @(depth-1) template< class C > void DownSample( int depth , const SortedTreeNodes& sNodes , ConstPointer( C ) fineConstraints , Pointer( C ) coarseConstraints ) const; // Up samples solution @(depth-1) to solution @(depth) template< class C > void UpSample ( int depth , const SortedTreeNodes& sNodes , ConstPointer( C ) coarseCoefficients , Pointer( C ) fineCoefficients ) const; - int GetSliceMatrixAndUpdateConstraints( const PointInfo& pointInfo , SparseMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine , int nStart , int nEnd ); - int GetMatrixAndUpdateConstraints( const PointInfo& pointInfo , SparseSymmetricMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine ); + int GetSliceMatrixAndUpdateConstraints( const SparseNodeData< PointData >& pointInfo , SparseMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine , int nStart , int nEnd ); + int GetMatrixAndUpdateConstraints( const SparseNodeData< PointData >& pointInfo , SparseSymmetricMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine ); int UpdateWeightContribution( std::vector< Real >& kernelDensityWeights , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , Real weight=Real(1.0) ); - Real GetSampleWeight( ConstPointer( Real ) kernelDensityWeight , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , int splatDepth ); - Real GetSampleWeight( ConstPointer( Real ) kernelDensityWeight , const TreeOctNode* node , const Point3D& position , typename TreeOctNode::ConstNeighborKey3& neighborKey ); - void GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeight , const TreeOctNode* node , const Point3D& position , typename TreeOctNode::ConstNeighborKey3& neighborKey , Real samplesPerNode , Real& depth , Real& weight ); - Real GetSampleWeight( ConstPointer( Real ) kernelDensityWeight , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey ); - void GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeight , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , Real samplesPerNode , Real& depth , Real& weight ); - int SplatOrientedPoint( ConstPointer( Real ) kernelDensityWeights , TreeOctNode* node , const Point3D& point , const Point3D< Real >& normal , NormalInfo& normalInfo , typename TreeOctNode::NeighborKey3& neighborKey ); - Real SplatOrientedPoint( ConstPointer( Real ) kernelDensityWeights , const Point3D& point , const Point3D& normal , NormalInfo& normalInfo , typename TreeOctNode::NeighborKey3& neighborKey , int kernelDepth , Real samplesPerNode , int minDepth , int maxDepth ); + Real GetSamplesPerNode( ConstPointer( Real ) kernelDensityWeight , const TreeOctNode* node , const Point3D< Real >& position , typename TreeOctNode::ConstNeighborKey3& neighborKey ); + Real GetSamplesPerNode( ConstPointer( Real ) kernelDensityWeight , TreeOctNode* node , const Point3D< Real >& position , typename TreeOctNode::NeighborKey3& neighborKey ); +public: + void GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeight , const Point3D< Real >& position , typename TreeOctNode::NeighborKey3& neighborKey , Real& depth , Real& weight ); + void GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeight , const Point3D< Real >& position , typename TreeOctNode::ConstNeighborKey3& neighborKey , Real& depth , Real& weight ); +protected: + void GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeight , const TreeOctNode* node , const Point3D< Real >& position , typename TreeOctNode::ConstNeighborKey3& neighborKey , Real& depth , Real& weight ); + void GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeight , TreeOctNode* node , const Point3D< Real >& position , typename TreeOctNode::NeighborKey3& neighborKey , Real& depth , Real& weight ); + template< class V > + int SplatPointData( TreeOctNode* node , const Point3D& point , const V& v , SparseNodeData< V >& data , typename TreeOctNode::NeighborKey3& neighborKey ); + template< class V > + int SplatPointData( TreeOctNode* node , const Point3D& point , const V& v , SparseNodeData< V >& data , typename TreeOctNode::ConstNeighborKey3& neighborKey ); + template< class V > + Real SplatPointData( ConstPointer( Real ) kernelDensityWeights , const Point3D< Real >& point , const V& v , SparseNodeData< V >& data , typename TreeOctNode::NeighborKey3& neighborKey , int minDepth , int maxDepth , int dim=DIMENSION ); + template< class V > + void MultiSplatPointData( ConstPointer( Real ) kernelDensityWeights , const Point3D< Real >& point , const V& v , SparseNodeData< V >& data , typename TreeOctNode::NeighborKey3& neighborKey , int maxDepth , int dim=DIMENSION ); + template< class V > + void MultiSplatPointData( ConstPointer( Real ) kernelDensityWeights , const Point3D< Real >& point , const V& v , SparseNodeData< V >& data , typename TreeOctNode::ConstNeighborKey3& neighborKey , int dim=DIMENSION ); - int HasNormals( TreeOctNode* node , const NormalInfo& normalInfo ); + int HasNormals( TreeOctNode* node , const SparseNodeData< Point3D< Real > >& normalInfo ); /////////////////////////// // Iso-Surfacing Methods // @@ -370,11 +391,11 @@ protected: template< class Vertex > void SetSliceIsoCorners( ConstPointer( Real ) solution , ConstPointer( Real ) coarseSolution , Real isoValue , int depth , int slice , int z , std::vector< SlabValues< Vertex > >& sValues , const typename BSplineData< 2 >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 > stencil[8] , const Stencil< double , 3 > stencils[8][8] , const Stencil< Point3D< double > , 3 > nStencil[8] , const Stencil< Point3D< double > , 3 > nStencils[8][8] , int threads ); template< class Vertex > - void SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , Real isoValue , int depth , int slice , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); + void SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , int depth , int slice , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); template< class Vertex > - void SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , Real isoValue , int depth , int slice , int z , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); + void SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , int depth , int slice , int z , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); template< class Vertex > - void SetXSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , Real isoValue , int depth , int slab , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); + void SetXSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , int depth , int slab , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); template< class Vertex > void CopyFinerSliceIsoEdgeKeys( int depth , int slice , std::vector< SlabValues< Vertex > >& sValues , int threads ); template< class Vertex > @@ -388,33 +409,35 @@ protected: template< class Vertex > void SetXSliceIsoEdges( int depth , int slice , std::vector< SlabValues< Vertex > >& slabValues , int threads ); - template< class Vertex > + template< class Vertex , 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 > + template< class Vertex , class _Vertex > static int AddIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector< std::pair< int , Vertex > >& polygon , bool polygonMesh , bool addBarycenter , int& vOffset ); template< class Vertex > - bool GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int edgeIndex , int z , const SliceValues< Vertex >& sValues , Vertex& vertex ); + bool GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int edgeIndex , int z , const SliceValues< Vertex >& sValues , Vertex& vertex ); template< class Vertex > - bool GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int cornerIndex , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , Vertex& vertex ); + bool GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > > >* colorData , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int cornerIndex , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , Vertex& vertex ); //////////////////////// // Evaluation Methods // //////////////////////// - Real getCornerValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int corner , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 > stencils[8] , bool isInterior ) const; + template< class V > + V getCornerValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey3 , 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; Point3D< Real > getCornerNormal( const typename TreeOctNode::ConstNeighbors5& neighbors5 , const typename TreeOctNode::ConstNeighbors5& pNeighbors5 , const TreeOctNode* node , int corner , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::template CornerEvaluator< 2 >& evaluator , const Stencil< Point3D< double > , 5 >& nStencil , const Stencil< Point3D< double > , 5 > nStencils[8] , bool isInterior ) const; std::pair< Real , Point3D< Real > > getCornerValueAndNormal( const typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int corner , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 >& vStencil , const Stencil< double , 3 > vStencils[8] , const Stencil< Point3D< double > , 3 >& nStencil , const Stencil< Point3D< double > , 3 > nStencils[8] , bool isInterior ) const; - Real getCenterValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::template CenterEvaluator< 1 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 >& pStencil , bool isInterior ) const; + template< class V > + V getCenterValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( V ) solution , ConstPointer( V ) metSolution , const typename BSplineData< 2 >::template CenterEvaluator< 1 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 >& pStencil , bool isInterior ) const; static bool _IsInset( const TreeOctNode* node ); static bool _IsInsetSupported( const TreeOctNode* node ); void refineBoundary( std::vector< int >* map ); public: - int threads; static double maxMemoryUsage; + int threads; TreeOctNode tree; static double MemoryUsage( void ); @@ -422,26 +445,45 @@ public: void MakeComplete( std::vector< int >* map=NULL ); void Finalize( std::vector< int >* map=NULL ); - void ClipTree( const NormalInfo& normalInfo ); + void ClipTree( const SparseNodeData< Point3D< Real > >& normalInfo ); - Real Evaluate( ConstPointer( Real ) coefficients , Point3D< Real > p , const BSplineData< 2 >* fData=NULL ) const; - Pointer( Real ) Evaluate( ConstPointer( Real ) coefficients , int& res , Real isoValue=0.f , int depth=-1 ); +protected: + template< class V > V _Evaluate( ConstPointer( V ) coefficients , Point3D< Real > p , typename TreeOctNode::ConstNeighborKey3& neighborKey3 ) const; + template< class V > V _Evaluate( const SparseNodeData< V >& coefficients , Point3D< Real > p , typename TreeOctNode::ConstNeighborKey3& neighborKey3 ) const; +public: + template< class V > V Evaluate( ConstPointer( V ) coefficients , Point3D< Real > p , const BSplineData< 2 >* fData = NULL , int depth=-1 ) const; + template< class V > V Evaluate( const SparseNodeData< V >& coefficients , Point3D< Real > p , const BSplineData< 2 >* fData = NULL , int depth=-1 ) const; + template< class V > Pointer( V ) Evaluate( ConstPointer( V ) coefficients , int& res , Real isoValue=0.f , int depth=-1 ); + // After calling set tree, the indices of the octree node will be stored by depth, and within depth they will be sorted by z-coordinate template< class PointReal > - int SetTree( PointStream< PointReal >* pointStream , int minDepth , int maxDepth , int fullDepth , int splatDepth , Real samplesPerNode , + int SetTree( OrientedPointStream< PointReal >* pointStream , int minDepth , int maxDepth , int fullDepth , int splatDepth , Real samplesPerNode , Real scaleFactor , bool useConfidence , bool useNormalWeight , Real constraintWeight , int adaptiveExponent , - PointInfo& pointInfo , NormalInfo& normalInfo , std::vector< Real >& kernelDensityWeights , std::vector< Real >& centerWeights , -#if 1 - // Bug fix courtesy of asmaloney - int boundaryType=BSplineElements< 2 >::NONE , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() , bool makeComplete=false ); -#else - int boundaryType=BSplineElements< 2 >::NONE , XForm4x4< Real > xForm=XForm4x4< Real >::Identity , bool makeComplete=false ); -#endif - Pointer( Real ) SetLaplacianConstraints( const NormalInfo& normalInfo ); - Pointer( Real ) SolveSystem( PointInfo& pointInfo , Pointer( Real ) constraints , bool showResidual , int iters , int maxSolveDepth , int cgDepth=0 , double cgAccuracy=0 ); + std::vector< Real >& kernelDensityWeights , + SparseNodeData< PointData >& pointInfo , SparseNodeData< Point3D< Real > >& normalInfo , std::vector< Real >& centerWeights , + XForm4x4< Real >& xForm , int boundaryType=BSplineElements< 2 >::NONE , bool makeComplete=false ); + + template< class PointReal , class Data , class _Data > + int SetTree( OrientedPointStreamWithData< PointReal , Data >* pointStream , int minDepth , int maxDepth , int fullDepth , int splatDepth , Real samplesPerNode , + Real scaleFactor , bool useConfidence , bool useNormalWeight , + std::vector< Real >& kernelDensityWeights , + SparseNodeData< ProjectiveData< _Data > >& dataValues , + XForm4x4< Real >& xForm , int boundaryType=BSplineElements< 2 >::NONE); + + template< class PointReal , class Data , class _Data > + int SetTree( OrientedPointStreamWithData< PointReal , Data >* pointStream , int minDepth , int maxDepth , int fullDepth , int splatDepth , Real samplesPerNode , + Real scaleFactor , bool useConfidence , bool useNormalWeight , + Real constraintWeight , int adaptiveExponent , + std::vector< Real >& kernelDensityWeights , + SparseNodeData< PointData >& pointInfo , SparseNodeData< Point3D< Real > >& normalInfo , std::vector< Real >& centerWeights , + SparseNodeData< ProjectiveData< _Data > >& dataValues , + XForm4x4< Real >& xForm , int boundaryType=BSplineElements< 2 >::NONE , bool makeComplete=false ); + + Pointer( Real ) SetLaplacianConstraints( const SparseNodeData< Point3D< Real > >& normalInfo ); + 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 > - void GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights , ConstPointer( Real ) solution , Real isoValue , CoredMeshData< Vertex >& mesh , bool nonLinearFit=true , bool addBarycenter=false , bool polygonMesh=false ); + template< class Vertex , 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 > void Reset( void ) diff --git a/Src/MultiGridOctreeData.inl b/Src/MultiGridOctreeData.inl index fc8d35b..8a007ca 100644 --- a/Src/MultiGridOctreeData.inl +++ b/Src/MultiGridOctreeData.inl @@ -67,7 +67,6 @@ template< class Real > Octree< Real >::Octree( void ) { threads = 1; - _normalSmooth = 0; _constrainValues = false; } @@ -88,7 +87,73 @@ bool Octree< Real >::_IsInsetSupported( const TreeOctNode* node ) return ( off[0]>=o && off[0]=o && off[1]=o && off[2] -int Octree< Real >::SplatOrientedPoint( ConstPointer( Real ) kernelDensityWeights , TreeOctNode* node , const Point3D& position , const Point3D& normal , NormalInfo& normalInfo , typename TreeOctNode::NeighborKey3& neighborKey ) +template< class V > +int Octree< Real >::SplatPointData( TreeOctNode* node , const Point3D< Real >& position , const V& v , SparseNodeData< V >& dataInfo , typename TreeOctNode::ConstNeighborKey3& neighborKey ) +{ + double x , dxdy , dxdydz , dx[DIMENSION][SPLAT_ORDER+1]; + double width; + int off[3]; + typename TreeOctNode::ConstNeighbors3& neighbors = neighborKey.getNeighbors( node ); + Point3D center; + Real w; + node->centerAndWidth( center , w ); + width=w; + for( int i=0 ; i<3 ; i++ ) + { +#if SPLAT_ORDER==2 + off[i] = 0; + x = ( center[i] - position[i] - width ) / width; + dx[i][0] = 1.125+1.500*x+0.500*x*x; + x = ( center[i] - position[i] ) / width; + dx[i][1] = 0.750 - x*x; + + dx[i][2] = 1. - dx[i][1] - dx[i][0]; +#elif SPLAT_ORDER==1 + x = ( position[i] - center[i] ) / width; + if( x<0 ) + { + off[i] = 0; + dx[i][0] = -x; + } + else + { + off[i] = 1; + dx[i][0] = 1. - x; + } + dx[i][1] = 1. - dx[i][0]; +#elif SPLAT_ORDER==0 + off[i] = 1; + dx[i][0] = 1.; +#else +# error Splat order not supported +#endif // SPLAT_ORDER + } + for( int i=off[0] ; i<=off[0]+SPLAT_ORDER ; i++ ) for( int j=off[1] ; j<=off[1]+SPLAT_ORDER ; j++ ) + { + dxdy = dx[0][i] * dx[1][j]; + for( int k=off[2] ; k<=off[2]+SPLAT_ORDER ; k++ ) + if( neighbors.neighbors[i][j][k] ) + { + dxdydz = dxdy * dx[2][k]; + const TreeOctNode* _node = neighbors.neighbors[i][j][k]; + if( _node ) + { + if( dataInfo.indices.size()nodeData.nodeIndex ] = (int)dataInfo.data.size(); + dataInfo.data.push_back( v * Real(dxdydz) ); + } + else dataInfo.data[idx] += v * Real( dxdydz ); + } + } + } + return 0; +} +template< class Real > +template< class V > +int Octree< Real >::SplatPointData( TreeOctNode* node , const Point3D< Real >& position , const V& v , SparseNodeData< V >& dataInfo , typename TreeOctNode::NeighborKey3& neighborKey ) { double x , dxdy , dxdydz , dx[DIMENSION][SPLAT_ORDER+1]; double width; @@ -136,170 +201,242 @@ int Octree< Real >::SplatOrientedPoint( ConstPointer( Real ) kernelDensityWeight { dxdydz = dxdy * dx[2][k]; TreeOctNode* _node = neighbors.neighbors[i][j][k]; - if( normalInfo.normalIndices.size() n; - n[0] = n[1] = n[2] = 0; - idx = normalInfo.normalIndices[ _node->nodeData.nodeIndex ] = (int)normalInfo.normals.size(); - normalInfo.normals.push_back( n ); + dataInfo.indices[ _node->nodeData.nodeIndex ] = (int)dataInfo.data.size(); + dataInfo.data.push_back( v * Real(dxdydz) ); } - normalInfo.normals[idx] += normal * Real( dxdydz ); + else dataInfo.data[idx] += v * Real( dxdydz ); } } return 0; } template< class Real > -Real Octree< Real >::SplatOrientedPoint( ConstPointer( Real ) kernelDensityWeights , const Point3D& position , const Point3D& normal , NormalInfo& normalInfo , typename TreeOctNode::NeighborKey3& neighborKey , int splatDepth , Real samplesPerNode , int minDepth , int maxDepth ) +template< class V > +Real Octree< Real >::SplatPointData( ConstPointer( Real ) kernelDensityWeights , const Point3D< Real >& position , const V& v , SparseNodeData< V >& dataInfo , typename TreeOctNode::NeighborKey3& neighborKey , int minDepth , int maxDepth , int dim ) { double dx; - Point3D n; + V _v; TreeOctNode* temp; int cnt=0; double width; - Point3D< Real > myCenter; - Real myWidth; - myCenter[0] = myCenter[1] = myCenter[2] = Real(0.5); - myWidth = Real(1.0); + Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 ); + Real myWidth = (Real)1.; temp = &tree; - while( temp->depth()depth()<_splatDepth ) { - if( !temp->children ) - { - fprintf( stderr , "Octree::SplatOrientedPoint error\n" ); - return -1; - } - int cIndex=TreeOctNode::CornerIndex(myCenter,position); + if( !temp->children ) fprintf( stderr , "[ERROR] Octree::SplatPointData\n" ) , exit( 0 ); + int cIndex=TreeOctNode::CornerIndex( myCenter , position ); temp=&temp->children[cIndex]; myWidth/=2; - if(cIndex&1) myCenter[0] += myWidth/2; - else myCenter[0] -= myWidth/2; - if(cIndex&2) myCenter[1] += myWidth/2; - else myCenter[1] -= myWidth/2; - if(cIndex&4) myCenter[2] += myWidth/2; - else myCenter[2] -= myWidth/2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; } Real weight , depth; - GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey , samplesPerNode , depth , weight ); + GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey , depth , weight ); - if( depthmaxDepth ) depth=Real(maxDepth); + if( depthmaxDepth ) depth = Real(maxDepth); int topDepth=int(ceil(depth)); dx = 1.0-(topDepth-depth); if( topDepth<=minDepth ) { - topDepth=minDepth; + topDepth = minDepth; dx=1; } else if( topDepth>maxDepth ) { - topDepth=maxDepth; + topDepth = maxDepth; dx=1; } while( temp->depth()>topDepth ) temp=temp->parent; while( temp->depth()children) temp->initChildren(); - int cIndex=TreeOctNode::CornerIndex(myCenter,position); + int cIndex=TreeOctNode::CornerIndex( myCenter , position ); temp=&temp->children[cIndex]; myWidth/=2; - if(cIndex&1) myCenter[0] += myWidth/2; - else myCenter[0] -= myWidth/2; - if(cIndex&2) myCenter[1] += myWidth/2; - else myCenter[1] -= myWidth/2; - if(cIndex&4) myCenter[2] += myWidth/2; - else myCenter[2] -= myWidth/2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; } width = 1.0 / ( 1<depth() ); - n = normal * weight / Real( pow( width , 3 ) ) * Real( dx ); - SplatOrientedPoint( kernelDensityWeights , temp , position , n , normalInfo , neighborKey ); + _v = v * weight / Real( pow( width , dim ) ) * Real( dx ); + SplatPointData( temp , position , _v , dataInfo , neighborKey ); if( fabs(1.0-dx) > EPSILON ) { dx = Real(1.0-dx); temp = temp->parent; width = 1.0 / ( 1<depth() ); - n = normal * weight / Real( pow( width , 3 ) ) * Real( dx ); - SplatOrientedPoint( kernelDensityWeights , temp , position , n , normalInfo , neighborKey ); + _v = v * weight / Real( pow( width , dim ) ) * Real( dx ); + SplatPointData( temp , position , _v , dataInfo , neighborKey ); } return weight; } +template< class Real > +template< class V > +void Octree< Real >::MultiSplatPointData( ConstPointer( Real ) kernelDensityWeights , const Point3D< Real >& position , const V& v , SparseNodeData< V >& dataInfo , typename TreeOctNode::NeighborKey3& neighborKey , int maxDepth , int dim ) +{ + Real _depth , weight; + if( kernelDensityWeights ) GetSampleDepthAndWeight( kernelDensityWeights , position , neighborKey , _depth , weight ); + else weight = (Real)1. , _depth = (Real)maxDepth; + int depth = std::min< int >( maxDepth , (int)ceil( _depth ) ); + V _v = v * weight; -template< class Real > -void Octree< Real >::GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeights , const TreeOctNode* node , const Point3D& position , typename TreeOctNode::ConstNeighborKey3& neighborKey , Real samplesPerNode , Real& depth , Real& weight ) -{ - const TreeOctNode* temp=node; - weight = Real(1.0)/GetSampleWeight( kernelDensityWeights , temp , position , neighborKey ); - if( weight>=samplesPerNode ) depth = Real( temp->depth() + log( weight / samplesPerNode ) / log(double(1<<(DIMENSION-1))) ); - else - { - Real oldWeight , newWeight; - oldWeight = newWeight = weight; - while( newWeightparent ) - { - temp=temp->parent; - oldWeight = newWeight; - newWeight = Real(1.0)/GetSampleWeight( kernelDensityWeights , temp , position, neighborKey ); - } - depth = Real( temp->depth() + log( newWeight / samplesPerNode ) / log( newWeight / oldWeight ) ); - } - weight = Real( pow( double(1<<(DIMENSION-1)) , -double(depth) ) ); -} -template< class Real > -void Octree< Real >::GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeights , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , Real samplesPerNode , Real& depth , Real& weight ) -{ - TreeOctNode* temp=node; - weight = Real(1.0)/GetSampleWeight( kernelDensityWeights , temp , position , neighborKey ); - if( weight>=samplesPerNode ) depth = Real( temp->depth() + log( weight / samplesPerNode ) / log(double(1<<(DIMENSION-1))) ); - else - { - Real oldWeight , newWeight; - oldWeight = newWeight = weight; - while( newWeightparent ) - { - temp=temp->parent; - oldWeight = newWeight; - newWeight = Real(1.0)/GetSampleWeight( kernelDensityWeights , temp , position, neighborKey ); - } - depth = Real( temp->depth() + log( newWeight / samplesPerNode ) / log( newWeight / oldWeight ) ); - } - weight = Real( pow( double(1<<(DIMENSION-1)) , -double(depth) ) ); -} -template< class Real > -Real Octree< Real >::GetSampleWeight( ConstPointer( Real ) kernelDensityWeights , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , int splatDepth ) -{ - Point3D< Real > myCenter; - Real myWidth; - myCenter[0] = myCenter[1] = myCenter[2] = Real(0.5); - myWidth = Real(1.0); + Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 ); + Real myWidth = (Real)1.; TreeOctNode* temp = &tree; - int d = 0; - while( ddepth()<=depth ) { - if( !temp->children ) + SplatPointData( temp , position , _v * Real( pow( 1<depth() , dim ) ) , dataInfo , neighborKey ); + if( temp->depth()children ) temp->initChildren(); + int cIndex = TreeOctNode::CornerIndex( myCenter , position ); + temp = &temp->children[cIndex]; + myWidth /= 2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; } + else break; + } +} +template< class Real > +template< class V > +void Octree< Real >::MultiSplatPointData( ConstPointer( Real ) kernelDensityWeights , const Point3D< Real >& position , const V& v , SparseNodeData< V >& dataInfo , typename TreeOctNode::ConstNeighborKey3& neighborKey , int dim ) +{ + Real depth , weight; + if( kernelDensityWeights ) GetSampleDepthAndWeight( kernelDensityWeights , position , neighborKey , depth , weight ); + else weight = (Real)1.; + V _v = v * weight; + + Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 ); + Real myWidth = (Real)1.; + + TreeOctNode* temp = &tree; + while( true ) + { + SplatPointData( temp , position , _v * Real( pow( 1<depth() , dim ) ) , dataInfo , neighborKey ); + if( !temp->children ) break; int cIndex = TreeOctNode::CornerIndex( myCenter , position ); temp = &temp->children[cIndex]; myWidth /= 2; if( cIndex&1 ) myCenter[0] += myWidth/2; - else myCenter[0] -= myWidth/2; + else myCenter[0] -= myWidth/2; if( cIndex&2 ) myCenter[1] += myWidth/2; - else myCenter[1] -= myWidth/2; + else myCenter[1] -= myWidth/2; if( cIndex&4 ) myCenter[2] += myWidth/2; else myCenter[2] -= myWidth/2; - d++; } - return GetSampleWeight( kernelDensityWeights , temp , position , neighborKey ); +} + +template< class Real > +void Octree< Real >::GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeights , const Point3D< Real >& position , typename TreeOctNode::ConstNeighborKey3& neighborKey , Real& depth , Real& weight ) +{ + TreeOctNode* temp; + Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 ); + Real myWidth = Real( 1. ); + + // Get the finest node with depth less than or equal to the splat depth that contains the point + temp = &tree; + while( temp->depth()<_splatDepth ) + { + if( !temp->children ) break;// fprintf( stderr , "[ERROR] Octree::GetSampleDepthAndWeight\n" ) , exit( 0 ); + int cIndex = TreeOctNode::CornerIndex( myCenter , position ); + temp = &temp->children[cIndex]; + myWidth /= 2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + } + return GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey , depth , weight ); } template< class Real > -Real Octree< Real >::GetSampleWeight( ConstPointer( Real ) kernelDensityWeights , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey ) +void Octree< Real >::GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeights , const Point3D< Real >& position , typename TreeOctNode::NeighborKey3& neighborKey , Real& depth , Real& weight ) +{ + TreeOctNode* temp; + Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 ); + Real myWidth = Real( 1. ); + + // Get the finest node with depth less than or equal to the splat depth that contains the point + temp = &tree; + while( temp->depth()<_splatDepth ) + { + if( !temp->children ) break;// fprintf( stderr , "[ERROR] Octree::GetSampleDepthAndWeight\n" ) , exit( 0 ); + int cIndex = TreeOctNode::CornerIndex( myCenter , position ); + temp = &temp->children[cIndex]; + myWidth /= 2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + } + return GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey , depth , weight ); +} +template< class Real > +void Octree< Real >::GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeights , const TreeOctNode* node , const Point3D& position , typename TreeOctNode::ConstNeighborKey3& neighborKey , Real& depth , Real& weight ) +{ + const TreeOctNode* temp=node; + weight = GetSamplesPerNode( kernelDensityWeights , temp , position , neighborKey ); + if( weight>=(Real)1. ) depth = Real( temp->depth() + log( weight ) / log(double(1<<(DIMENSION-1))) ); + else + { + Real oldWeight , newWeight; + oldWeight = newWeight = weight; + while( newWeight<(Real)1. && temp->parent ) + { + temp=temp->parent; + oldWeight = newWeight; + newWeight = GetSamplesPerNode( kernelDensityWeights , temp , position, neighborKey ); + } + depth = Real( temp->depth() + log( newWeight ) / log( newWeight / oldWeight ) ); + } + weight = Real( pow( double(1<<(DIMENSION-1)) , -double(depth) ) ); +} +template< class Real > +void Octree< Real >::GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeights , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , Real& depth , Real& weight ) +{ + TreeOctNode* temp=node; + weight = GetSamplesPerNode( kernelDensityWeights , temp , position , neighborKey ); + if( weight>=(Real)1. ) depth = Real( temp->depth() + log( weight ) / log(double(1<<(DIMENSION-1))) ); + else + { + Real oldWeight , newWeight; + oldWeight = newWeight = weight; + while( newWeight<(Real)1. && temp->parent ) + { + temp=temp->parent; + oldWeight = newWeight; + newWeight = GetSamplesPerNode( kernelDensityWeights , temp , position, neighborKey ); + } + depth = Real( temp->depth() + log( newWeight ) / log( newWeight / oldWeight ) ); + } + weight = Real( pow( double(1<<(DIMENSION-1)) , -double(depth) ) ); +} +template< class Real > +Real Octree< Real >::GetSamplesPerNode( ConstPointer( Real ) kernelDensityWeights , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey ) { Real weight=0; double x , dxdy , dx[DIMENSION][3]; @@ -307,8 +444,8 @@ Real Octree< Real >::GetSampleWeight( ConstPointer( Real ) kernelDensityWeights typename TreeOctNode::Neighbors3& neighbors = neighborKey.setNeighbors( node ); Point3D center; Real w; - node->centerAndWidth(center,w); - width=w; + node->centerAndWidth( center , w ); + width = w; for( int i=0 ; i::GetSampleWeight( ConstPointer( Real ) kernelDensityWeights for( int k=0 ; k<3 ; k++ ) if( neighbors.neighbors[i][j][k] ) weight += Real( dxdy * dx[2][k] * kernelDensityWeights[ neighbors.neighbors[i][j][k]->nodeData.nodeIndex ] ); } - return Real( 1.0 / weight ); + return weight; } template< class Real > -Real Octree< Real >::GetSampleWeight( ConstPointer( Real ) kernelDensityWeights , const TreeOctNode* node , const Point3D& position , typename TreeOctNode::ConstNeighborKey3& neighborKey ) +Real Octree< Real >::GetSamplesPerNode( ConstPointer( Real ) kernelDensityWeights , const TreeOctNode* node , const Point3D& position , typename TreeOctNode::ConstNeighborKey3& neighborKey ) { Real weight=0; double x,dxdy,dx[DIMENSION][3]; @@ -356,7 +493,7 @@ Real Octree< Real >::GetSampleWeight( ConstPointer( Real ) kernelDensityWeights for( int k=0 ; k<3 ; k++ ) if( neighbors.neighbors[i][j][k] ) weight += Real( dxdy * dx[2][k] * kernelDensityWeights[ neighbors.neighbors[i][j][k]->nodeData.nodeIndex ] ); } - return Real( 1.0 / weight ); + return weight; } template< class Real > int Octree< Real >::UpdateWeightContribution( std::vector< Real >& kernelDensityWeights , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , Real weight ) @@ -367,8 +504,7 @@ int Octree< Real >::UpdateWeightContribution( std::vector< Real >& kernelDensity Point3D< Real > center; Real w; node->centerAndWidth( center , w ); - width=w; - const double SAMPLE_SCALE = 1. / ( 0.125 * 0.125 + 0.75 * 0.75 + 0.125 * 0.125 ); + width = w; for( int i=0 ; i::UpdateWeightContribution( std::vector< Real >& kernelDensity dx[i][0] = 1.125 + 1.500*x + 0.500*x*x; dx[i][1] = -0.25 - 2.*x - x*x; dx[i][2] = 1. - dx[i][1] - dx[i][0]; - // Note that we are splatting along a co-dimension one manifold, so uniform point samples - // do not generate a unit sample weight. - dx[i][0] *= SAMPLE_SCALE; } + + // Suppose that the samples are uniformly placed along the middle of the three slices. + // Then splatting the points we get coefficients: + // 0.125 / 0.75 / 0.125 across the three slices. + // Sampling at the center slice we get: + // 0.125^2 + 0.75^2 + 0.125^2 = 19/32 + const double SAMPLE_SCALE = 1. / ( 0.125 * 0.125 + 0.75 * 0.75 + 0.125 * 0.125 ); + weight *= (Real)SAMPLE_SCALE; + for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) { dxdy = dx[0][i] * dx[1][j] * weight; @@ -397,18 +539,17 @@ bool Octree< Real >::_InBounds( Point3D< Real > p ) const } template< class Real > template< class PointReal > -int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDepth , int maxDepth , int fullDepth , +int Octree< Real >::SetTree( OrientedPointStream< PointReal >* pointStream , int minDepth , int maxDepth , int fullDepth , int splatDepth , Real samplesPerNode , Real scaleFactor , bool useConfidence , bool useNormalWeights , Real constraintWeight , int adaptiveExponent , - PointInfo& pointInfo , NormalInfo& normalInfo , std::vector< Real >& kernelDensityWeights , std::vector< Real >& centerWeights , - int boundaryType , XForm4x4< Real > xForm , bool makeComplete ) + std::vector< Real >& kernelDensityWeights , SparseNodeData< PointData >& pointInfo , SparseNodeData< Point3D< Real > >& normalInfo , std::vector< Real >& centerWeights , + XForm4x4< Real >& xForm , int boundaryType , bool makeComplete ) { if( splatDepth<0 ) splatDepth = 0; _boundaryType = boundaryType; if ( _boundaryType<0 ) _boundaryType = -1; else if( _boundaryType>0 ) _boundaryType = 1; - _samplesPerNode = samplesPerNode; _splatDepth = splatDepth; _constrainValues = (constraintWeight>0); @@ -451,10 +592,10 @@ int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDept { double t = Time(); Point3D< Real > p; - Point3D< PointReal > _p , _n; - while( pointStream->nextPoint( _p , _n ) ) + OrientedPoint3D< PointReal > _p; + while( pointStream->nextPoint( _p ) ) { - p = xForm * Point3D< Real >(_p); + p = xForm * Point3D< Real >(_p.p); for( i=0 ; i::SetTree( PointStream< PointReal >* pointStream , int minDept _scale *= scaleFactor; for( i=0 ; i trans = XForm4x4< Real >::Identity() , scale = XForm4x4< Real >::Identity(); + for( int i=0 ; i<3 ; i++ ) scale(i,i) = (Real)(1./_scale ) , trans(3,i) = -_center[i]; + xForm = scale * trans * xForm; + } + if( splatDepth>0 ) { double t = Time(); cnt = 0; pointStream->reset(); Point3D< Real > p , n; - Point3D< PointReal > _p , _n; - while( pointStream->nextPoint( _p , _n ) ) + OrientedPoint3D< PointReal > _p; + while( pointStream->nextPoint( _p ) ) { - p = xForm * Point3D< Real >(_p) , n = xFormN * Point3D< Real >(_n); - p = ( p - _center ) / _scale; + p = xForm * Point3D< Real >(_p.p) , n = xFormN * Point3D< Real >(_p.n); if( !_InBounds(p) ) continue; myCenter = Point3D< Real >( Real(0.5) , Real(0.5) , Real(0.5) ); myWidth = Real(1.0); Real weight=Real( 1. ); if( useConfidence ) weight = Real( Length(n) ); + if( samplesPerNode>0 ) weight /= (Real)samplesPerNode; temp = &tree; int d=0; while( d::SetTree( PointStream< PointReal >* pointStream , int minDept } kernelDensityWeights.resize( TreeNodeData::NodeCount , 0 ); - std::vector< _PointData >& points = pointInfo.points; + std::vector< PointData >& points = pointInfo.data; cnt = 0; pointStream->reset(); Point3D< Real > p , n; - Point3D< PointReal > _p , _n; - while( pointStream->nextPoint( _p , _n ) ) + OrientedPoint3D< PointReal > _p; + while( pointStream->nextPoint( _p ) ) { - p = xForm * Point3D< Real >(_p) , n = xFormN * Point3D< Real >(_n); + p = xForm * Point3D< Real >(_p.p) , n = xFormN * Point3D< Real >(_p.n); n *= Real(-1.); - p = ( p - _center ) / _scale; if( !_InBounds(p) ) continue; myCenter = Point3D< Real >( Real(0.5) , Real(0.5) , Real(0.5) ); myWidth = Real(1.0); Real normalLength = Real( Length( n ) ); - if( normalLength!=normalLength || normalLength<=EPSILON ) continue; + if( isnan( normalLength ) || !isfinite( normalLength ) || normalLength<=EPSILON ) continue; if( !useConfidence ) n /= normalLength; Real pointWeight = Real(1.f); - if( samplesPerNode>0 && splatDepth ) pointWeight = SplatOrientedPoint( GetPointer( kernelDensityWeights ) , p , n , normalInfo , neighborKey , splatDepth , samplesPerNode , _minDepth , maxDepth ); + if( samplesPerNode>0 && splatDepth ) pointWeight = SplatPointData( GetPointer( kernelDensityWeights ) , p , n , normalInfo , neighborKey , _minDepth , maxDepth ); else { temp = &tree; @@ -548,7 +696,7 @@ int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDept else myCenter[2] -= myWidth/2; d++; } - pointWeight = GetSampleWeight( GetPointer( kernelDensityWeights ) , temp , p , neighborKey ); + pointWeight = (Real)1.0/GetSamplesPerNode( GetPointer( kernelDensityWeights ) , temp , p , neighborKey ); } for( i=0 ; i::SetTree( PointStream< PointReal >* pointStream , int minDept else myCenter[2] -= myWidth/2; d++; } - SplatOrientedPoint( GetPointer( kernelDensityWeights ) , temp , p , n , normalInfo , neighborKey ); + SplatPointData( temp , p , n , normalInfo , neighborKey ); } pointWeightSum += pointWeight; if( _constrainValues ) @@ -577,14 +725,14 @@ int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDept myWidth = Real(1.0); while( 1 ) { - if( pointInfo.pointIndices.size()nodeData.nodeIndex ] = idx; + points.push_back( PointData( p*pointScreeningWeight , pointScreeningWeight ) ); + pointInfo.indices[ temp->nodeData.nodeIndex ] = idx; } else { @@ -616,9 +764,9 @@ int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDept if( _constrainValues ) // Set the average position and scale the weights for( TreeOctNode* node=tree.nextNode() ; node ; node=tree.nextNode(node) ) - if( pointInfo.pointIndex( node )!=-1 ) + if( pointInfo.index( node )!=-1 ) { - int idx = pointInfo.pointIndex( node ); + int idx = pointInfo.index( node ); points[idx].position /= points[idx].weight; int e = ( _boundaryType==0 ? node->depth()-1 : node->depth() ) * adaptiveExponent - ( _boundaryType==0 ? maxDepth-1 : maxDepth ) * (adaptiveExponent-1); if( e<0 ) points[idx].weight /= Real( 1<<(-e) ); @@ -632,9 +780,9 @@ int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDept int d , off[3] , res; node->depthAndOffset( d , off ); res = 1<& normal = normalInfo.normals[ idx ]; + Point3D< Real >& normal = normalInfo.data[ idx ]; for( int d=0 ; d<3 ; d++ ) if( off[d]==0 || off[d]==res-1 ) normal[d] = 0; } #endif // FORCE_NEUMANN_FIELD @@ -643,9 +791,9 @@ int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDept // Set the point weights for evaluating the iso-value for( TreeOctNode* node=tree.nextNode() ; node ; node=tree.nextNode(node) ) { - int idx = normalInfo.normalIndex( node ); + int idx = normalInfo.index( node ); if( idx<0 ) centerWeights[ node->nodeData.nodeIndex ] = 0; - else centerWeights[ node->nodeData.nodeIndex ] = Real( Length( normalInfo.normals[ idx ] ) ); + else centerWeights[ node->nodeData.nodeIndex ] = Real( Length( normalInfo.data[ idx ] ) ); } MemoryUsage(); { @@ -654,18 +802,18 @@ int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDept else ClipTree( normalInfo ) , Finalize( &indexMap ); { - std::vector< int > temp = pointInfo.pointIndices; - pointInfo.pointIndices.resize( indexMap.size() ); + std::vector< int > temp = pointInfo.indices; + pointInfo.indices.resize( indexMap.size() ); for( int i=0 ; i temp = normalInfo.normalIndices; - normalInfo.normalIndices.resize( indexMap.size() ); + std::vector< int > temp = normalInfo.indices; + normalInfo.indices.resize( indexMap.size() ); for( int i=0 ; i temp = centerWeights; @@ -685,6 +833,501 @@ int Octree< Real >::SetTree( PointStream< PointReal >* pointStream , int minDept return cnt; } template< class Real > +template< class PointReal , class Data , class _Data > +int Octree< Real >::SetTree( OrientedPointStreamWithData< PointReal , Data >* pointStream , int minDepth , int maxDepth , int fullDepth , + int splatDepth , Real samplesPerNode , Real scaleFactor , + bool useConfidence , bool useNormalWeights , Real constraintWeight , int adaptiveExponent , + std::vector< Real >& kernelDensityWeights , SparseNodeData< PointData >& pointInfo , SparseNodeData< Point3D< Real > >& normalInfo , std::vector< Real >& centerWeights , + SparseNodeData< ProjectiveData< _Data > >& dataValues , + XForm4x4< Real >& xForm , int boundaryType , bool makeComplete ) +{ + if( splatDepth<0 ) splatDepth = 0; + + _boundaryType = boundaryType; + if ( _boundaryType<0 ) _boundaryType = -1; + else if( _boundaryType>0 ) _boundaryType = 1; + + _splatDepth = splatDepth; + _constrainValues = (constraintWeight>0); + + XForm3x3< Real > xFormN; + for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) xFormN(i,j) = xForm(i,j); + xFormN = xFormN.transpose().inverse(); + minDepth = std::min< int >( minDepth , maxDepth ); // minDepth <= maxDepth + fullDepth = std::max< int >( minDepth , std::min< int >( fullDepth , maxDepth ) ); // minDepth <= fullDepth <= maxDepth + // If _boundaryType==0, points are scaled to be in the [0.25,0.75]^3 cube so all depths have to be offset by + // and the minDepth has to be 2. + if( _boundaryType==0 ) + { + minDepth++ , maxDepth++ , fullDepth++; + if( splatDepth ) splatDepth++; + minDepth = std::max< int >( minDepth , 2 ); + } + // Otherwise the points are in the [0,1]^3 cube. + // However, for Neumann constraints, the function at depth 0 is constant so the system matrix is zero if there + // is no screening. +#if 0 + else if( _boundaryType==1 && !_constrainValues ) minDepth = std::max< int >( minDepth , 1 ); +#endif + + _fData.set( maxDepth , _boundaryType ); + + _minDepth = minDepth; + _fullDepth = fullDepth; + double pointWeightSum = 0; + Point3D< Real > min , max , myCenter; + Real myWidth; + int i , cnt=0; + TreeOctNode* temp; + + typename TreeOctNode::NeighborKey3 neighborKey; + neighborKey.set( maxDepth ); + + tree.setFullDepth( _fullDepth ); + + // Read through once to get the center and scale + { + double t = Time(); + Point3D< Real > p; + OrientedPoint3D< PointReal > _p; + while( pointStream->nextPoint( _p ) ) + { + p = xForm * Point3D< Real >(_p.p); + for( i=0 ; imax[i] ) max[i] = p[i]; + } + cnt++; + } + + if( _boundaryType==0 ) _scale = std::max< Real >( max[0]-min[0] , std::max< Real >( max[1]-min[1] , max[2]-min[2] ) ) * 2; + else _scale = std::max< Real >( max[0]-min[0] , std::max< Real >( max[1]-min[1] , max[2]-min[2] ) ); + _center = ( max+min ) /2; + } + + _scale *= scaleFactor; + for( i=0 ; i trans = XForm4x4< Real >::Identity() , scale = XForm4x4< Real >::Identity(); + for( int i=0 ; i<3 ; i++ ) scale(i,i) = (Real)(1./_scale ) , trans(3,i) = -_center[i]; + xForm = scale * trans * xForm; + } + + if( splatDepth>0 ) + { + double t = Time(); + cnt = 0; + pointStream->reset(); + Point3D< Real > p , n; + OrientedPoint3D< PointReal > _p; + while( pointStream->nextPoint( _p ) ) + { + p = xForm * Point3D< Real >(_p.p) , n = xFormN * Point3D< Real >(_p.n); + if( !_InBounds(p) ) continue; + myCenter = Point3D< Real >( Real(0.5) , Real(0.5) , Real(0.5) ); + myWidth = Real(1.0); + Real weight=Real( 1. ); + if( useConfidence ) weight = Real( Length(n) ); + if( samplesPerNode>0 ) weight /= (Real)samplesPerNode; + temp = &tree; + int d=0; + while( dchildren ) temp->initChildren(); + int cIndex=TreeOctNode::CornerIndex( myCenter , p ); + temp = temp->children + cIndex; + myWidth/=2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + d++; + } + UpdateWeightContribution( kernelDensityWeights , temp , p , neighborKey , weight ); + cnt++; + } + } + kernelDensityWeights.resize( TreeNodeData::NodeCount , 0 ); + + std::vector< PointData >& points = pointInfo.data; + + cnt = 0; + pointStream->reset(); + Point3D< Real > p , n; + OrientedPoint3D< PointReal > _p; + Data _d; + while( pointStream->nextPoint( _p , _d ) ) + { + p = xForm * Point3D< Real >(_p.p) , n = xFormN * Point3D< Real >(_p.n); + n *= Real(-1.); + if( !_InBounds(p) ) continue; + myCenter = Point3D< Real >( Real(0.5) , Real(0.5) , Real(0.5) ); + myWidth = Real(1.0); + Real normalLength = Real( Length( n ) ); + if( normalLength!=normalLength || normalLength<=EPSILON ) continue; + if( !useConfidence ) n /= normalLength; + + Real pointWeight = Real(1.f); + if( samplesPerNode>0 && splatDepth ) + { + MultiSplatPointData( GetPointer( kernelDensityWeights ) , p , ProjectiveData< _Data >( _Data( _d ) , (Real)1. ) , dataValues , neighborKey , maxDepth , 2 ); + pointWeight = SplatPointData( GetPointer( kernelDensityWeights ) , p , n , normalInfo , neighborKey , _minDepth , maxDepth , 3 ); + } + else + { + MultiSplatPointData( NullPointer( Real ) , p , ProjectiveData< _Data >( _Data( _d ) , (Real)1. ) , dataValues , neighborKey , maxDepth , 2 ); + + temp = &tree; + int d=0; + if( splatDepth ) + { + while( dchildren[cIndex]; + myWidth /= 2; + if(cIndex&1) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if(cIndex&2) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if(cIndex&4) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + d++; + } + pointWeight = (Real)1.0/GetSamplesPerNode( GetPointer( kernelDensityWeights ) , temp , p , neighborKey ); + } + for( i=0 ; ichildren ) temp->initChildren(); + int cIndex=TreeOctNode::CornerIndex( myCenter , p ); + temp=&temp->children[cIndex]; + myWidth/=2; + if(cIndex&1) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if(cIndex&2) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if(cIndex&4) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + d++; + } + SplatPointData( temp , p , n , normalInfo , neighborKey ); + } + pointWeightSum += pointWeight; + if( _constrainValues ) + { + Real pointScreeningWeight = useNormalWeights ? Real( normalLength ) : Real(1.f); + int d = 0; + TreeOctNode* temp = &tree; + myCenter = Point3D< Real >( Real(0.5) , Real(0.5) , Real(0.5) ); + myWidth = Real(1.0); + while( 1 ) + { + if( pointInfo.indices.size()nodeData.nodeIndex ] = idx; + } + points[idx].addPoint( p*pointScreeningWeight , myCenter , myWidth , pointScreeningWeight ); +#else // !POINT_DATA_RES + if( idx==-1 ) + { + idx = (int)points.size(); + points.push_back( PointData( p*pointScreeningWeight , pointScreeningWeight ) ); + pointInfo.indices[ temp->nodeData.nodeIndex ] = idx; + } + else + { + points[idx].weight += pointScreeningWeight; + points[idx].position += p*pointScreeningWeight; + } +#endif // POINT_DATA_RES + + int cIndex = TreeOctNode::CornerIndex( myCenter , p ); + if( !temp->children ) break; + temp = &temp->children[cIndex]; + myWidth /= 2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + d++; + } + } + cnt++; + } + + if( _boundaryType==0 ) pointWeightSum *= Real(4.); + constraintWeight *= Real( pointWeightSum ); + constraintWeight /= cnt; + + MemoryUsage( ); + if( _constrainValues ) + // Set the average position and scale the weights + for( TreeOctNode* node=tree.nextNode() ; node ; node=tree.nextNode(node) ) + if( pointInfo.index( node )!=-1 ) + { + int idx = pointInfo.index( node ); +#if POINT_DATA_RES + for( int c=0 ; c<_PointData::SAMPLES ; c++ ) + if( points[idx].weights[c] ) + { + points[idx].positions[c] /= points[idx].weights[c]; + int e = ( _boundaryType==0 ? node->depth()-1 : node->depth() ) * adaptiveExponent - ( _boundaryType==0 ? maxDepth-1 : maxDepth ) * (adaptiveExponent-1); + if( e<0 ) points[idx].weights[c] /= Real( 1<<(-e) ); + else points[idx].weights[c] *= Real( 1<< e ); + points[idx].weights[c] *= Real( constraintWeight ); + } +#else // !POINT_DATA_RES + points[idx].position /= points[idx].weight; + int e = ( _boundaryType==0 ? node->depth()-1 : node->depth() ) * adaptiveExponent - ( _boundaryType==0 ? maxDepth-1 : maxDepth ) * (adaptiveExponent-1); + if( e<0 ) points[idx].weight /= Real( 1<<(-e) ); + else points[idx].weight *= Real( 1<< e ); + points[idx].weight *= Real( constraintWeight ); +#endif // POINT_DATA_RES + } +#if FORCE_NEUMANN_FIELD + if( _boundaryType==1 ) + for( TreeOctNode* node=tree.nextNode() ; node ; node=tree.nextNode( node ) ) + { + int d , off[3] , res; + node->depthAndOffset( d , off ); + res = 1<& normal = normalInfo.data[ idx ]; + for( int d=0 ; d<3 ; d++ ) if( off[d]==0 || off[d]==res-1 ) normal[d] = 0; + } +#endif // FORCE_NEUMANN_FIELD + centerWeights.resize( tree.nodes() , 0 ); + kernelDensityWeights.resize( tree.nodes() , 0 ); + // Set the point weights for evaluating the iso-value + for( TreeOctNode* node=tree.nextNode() ; node ; node=tree.nextNode(node) ) + { + int idx = normalInfo.index( node ); + if( idx<0 ) centerWeights[ node->nodeData.nodeIndex ] = 0; + else centerWeights[ node->nodeData.nodeIndex ] = Real( Length( normalInfo.data[ idx ] ) ); + } + MemoryUsage(); + { + std::vector< int > indexMap; + if( makeComplete ) MakeComplete( &indexMap ); + else ClipTree( normalInfo ) , Finalize( &indexMap ); + + { + std::vector< int > temp = pointInfo.indices; + pointInfo.indices.resize( indexMap.size() ); + for( int i=0 ; i temp = normalInfo.indices; + normalInfo.indices.resize( indexMap.size() ); + for( int i=0 ; i temp = centerWeights; + centerWeights.resize( indexMap.size() ); + for( int i=0 ; i temp = kernelDensityWeights; + kernelDensityWeights.resize( indexMap.size() ); + for( int i=0 ; i temp = dataValues.indices; + dataValues.indices.resize( indexMap.size() ); + for( int i=0 ; i +template< class PointReal , class Data , class _Data > +int Octree< Real >::SetTree( OrientedPointStreamWithData< PointReal , Data >* pointStream , int minDepth , int maxDepth , int fullDepth , + int splatDepth , Real samplesPerNode , Real scaleFactor , + bool useConfidence , bool useNormalWeights , + std::vector< Real >& kernelDensityWeights , SparseNodeData< ProjectiveData< _Data > >& dataValues , + XForm4x4< Real >& xForm , int boundaryType ) +{ + if( splatDepth<0 ) splatDepth = 0; + + _boundaryType = boundaryType; + if ( _boundaryType<0 ) _boundaryType = -1; + else if( _boundaryType>0 ) _boundaryType = 1; + _splatDepth = splatDepth; + + XForm3x3< Real > xFormN; + for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) xFormN(i,j) = xForm(i,j); + xFormN = xFormN.transpose().inverse(); + minDepth = std::min< int >( minDepth , maxDepth ); // minDepth <= maxDepth + fullDepth = std::max< int >( minDepth , std::min< int >( fullDepth , maxDepth ) ); // minDepth <= fullDepth <= maxDepth + // If _boundaryType==0, points are scaled to be in the [0.25,0.75]^3 cube so all depths have to be offset by + // and the minDepth has to be 2. + if( _boundaryType==0 ) + { + minDepth++ , maxDepth++ , fullDepth++; + if( splatDepth ) splatDepth++; + minDepth = std::max< int >( minDepth , 2 ); + } + // Otherwise the points are in the [0,1]^3 cube. + // However, for Neumann constraints, the function at depth 0 is constant so the system matrix is zero if there + // is no screening. +#if 0 + else if( _boundaryType==1 ) minDepth = std::max< int >( minDepth , 1 ); +#endif + + _fData.set( maxDepth , _boundaryType ); + + _minDepth = minDepth; + _fullDepth = fullDepth; + Point3D< Real > min , max , myCenter; + Real myWidth; + int i , cnt=0; + TreeOctNode* temp; + + typename TreeOctNode::NeighborKey3 neighborKey; + neighborKey.set( maxDepth ); + + tree.setFullDepth( _fullDepth ); + + // Read through once to get the center and scale + { + double t = Time(); + Point3D< Real > p; + OrientedPoint3D< PointReal > _p; + while( pointStream->nextPoint( _p ) ) + { + p = xForm * Point3D< Real >(_p.p); + for( i=0 ; imax[i] ) max[i] = p[i]; + } + cnt++; + } + + if( _boundaryType==0 ) _scale = std::max< Real >( max[0]-min[0] , std::max< Real >( max[1]-min[1] , max[2]-min[2] ) ) * 2; + else _scale = std::max< Real >( max[0]-min[0] , std::max< Real >( max[1]-min[1] , max[2]-min[2] ) ); + _center = ( max+min ) /2; + } + + _scale *= scaleFactor; + for( i=0 ; i trans = XForm4x4< Real >::Identity() , scale = XForm4x4< Real >::Identity(); + for( int i=0 ; i<3 ; i++ ) scale(i,i) = (Real)(1./_scale ) , trans(3,i) = -_center[i]; + xForm = scale * trans * xForm; + } + + if( splatDepth>0 ) + { + double t = Time(); + cnt = 0; + pointStream->reset(); + Point3D< Real > p , n; + OrientedPoint3D< PointReal > _p; + while( pointStream->nextPoint( _p ) ) + { + p = xForm * Point3D< Real >(_p.p) , n = xFormN * Point3D< Real >(_p.n); + if( !_InBounds(p) ) continue; + myCenter = Point3D< Real >( Real(0.5) , Real(0.5) , Real(0.5) ); + myWidth = Real(1.0); + Real weight=Real( 1. ); + if( useConfidence ) weight = Real( Length(n) ); + if( samplesPerNode ) weight /= (Real)samplesPerNode; + temp = &tree; + int d=0; + while( dchildren ) temp->initChildren(); + int cIndex=TreeOctNode::CornerIndex( myCenter , p ); + temp = temp->children + cIndex; + myWidth/=2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + d++; + } + UpdateWeightContribution( kernelDensityWeights , temp , p , neighborKey , weight ); + cnt++; + } + } + kernelDensityWeights.resize( TreeNodeData::NodeCount , 0 ); + + cnt = 0; + pointStream->reset(); + Point3D< Real > p; + OrientedPoint3D< PointReal > _p; + Data _d; + while( pointStream->nextPoint( _p , _d ) ) + { + p = xForm * Point3D< Real >(_p.p); + if( !_InBounds(p) ) continue; + myCenter = Point3D< Real >( Real(0.5) , Real(0.5) , Real(0.5) ); + myWidth = Real(1.0); + + if( samplesPerNode>0 && splatDepth ) MultiSplatPointData( GetPointer( kernelDensityWeights ) , p , ProjectiveData< _Data >( _Data( _d ) , (Real)1. ) , dataValues , neighborKey , maxDepth , 2 ); + else MultiSplatPointData( NullPointer( Real ) , p , ProjectiveData< _Data >( _Data( _d ) , (Real)1. ) , dataValues , neighborKey , maxDepth , 2 ); + cnt++; + } + + MemoryUsage( ); + kernelDensityWeights.resize( tree.nodes() , 0 ); + // Set the point weights for evaluating the iso-value + MemoryUsage(); + { + std::vector< int > indexMap; + Finalize( &indexMap ); + + { + std::vector< int > temp = dataValues.indices; + dataValues.indices.resize( indexMap.size() ); + for( int i=0 ; i temp = kernelDensityWeights; + kernelDensityWeights.resize( indexMap.size() ); + for( int i=0 ; i void Octree< Real >::MakeComplete( std::vector< int >* map ) { tree.setFullDepth( tree.maxDepth() ); @@ -692,7 +1335,7 @@ void Octree< Real >::MakeComplete( std::vector< int >* map ) MemoryUsage(); } template< class Real > -void Octree< Real >::ClipTree( const NormalInfo& normalInfo ) +void Octree< Real >::ClipTree( const SparseNodeData< Point3D< Real > >& normalInfo ) { int maxDepth = tree.maxDepth(); for( TreeOctNode* temp=tree.nextNode() ; temp ; temp=tree.nextNode(temp) ) @@ -825,9 +1468,9 @@ int Octree< Real >::GetMatrixRowSize( const typename TreeOctNode::Neighbors5& ne } template< class Real > -int Octree< Real >::SetMatrixRow( const PointInfo& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , Pointer( MatrixEntry< Real > ) row , int offset , const typename BSplineData< 2 >::Integrator& integrator , const Stencil< double , 5 >& stencil , bool symmetric ) const +int Octree< Real >::SetMatrixRow( const SparseNodeData< PointData >& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , Pointer( MatrixEntry< Real > ) row , int offset , const typename BSplineData< 2 >::Integrator& integrator , const Stencil< double , 5 >& stencil , bool symmetric ) const { - const std::vector< _PointData >& points = pointInfo.points; + const std::vector< PointData >& points = pointInfo.data; bool hasYZPoints[3] , hasZPoints[3][3]; Real diagonal = 0; Real splineValues[3*3*3*3*3]; @@ -857,9 +1500,9 @@ int Octree< Real >::SetMatrixRow( const PointInfo& pointInfo , const typename Tr for( int l=0 ; l<3 ; l++ ) { const TreeOctNode* _node = neighbors5.neighbors[j+1][k+1][l+1]; - if( _node && pointInfo.pointIndex( _node )!=-1 ) + if( _node && pointInfo.index( _node )!=-1 ) { - const _PointData& pData = points[ pointInfo.pointIndex( _node ) ]; + const PointData& pData = points[ pointInfo.index( _node ) ]; Real* _splineValues = splineValues + 3*3*(3*(3*j+k)+l); Real weight = pData.weight; Point3D< Real > p = pData.position; @@ -895,7 +1538,7 @@ int Octree< Real >::SetMatrixRow( const PointInfo& pointInfo , const typename Tr for( int k=0 ; k<3 ; k++ ) { const TreeOctNode* _node = neighbors5.neighbors[i+1][j+1][k+1]; - if( _node && pointInfo.pointIndex( _node )!=-1 ) + if( _node && pointInfo.index( _node )!=-1 ) { const Real* _splineValuesX = splineValues + 3*(3*(3*(3*i+j)+k)+0)+2; const Real* _splineValuesY = splineValues + 3*(3*(3*(3*i+j)+k)+1)+2; @@ -1238,9 +1881,9 @@ void Octree< Real >::UpdateConstraintsFromFiner( const typename BSplineData< 2 > } template< class Real > -void Octree< Real >::UpdateConstraintsFromCoarser( const PointInfo& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , const typename TreeOctNode::Neighbors5& pNeighbors5 , TreeOctNode* node , Pointer( Real ) constraints , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::Integrator& integrator , const Stencil< double , 5 >& lapStencil ) const +void Octree< Real >::UpdateConstraintsFromCoarser( const SparseNodeData< PointData >& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , const typename TreeOctNode::Neighbors5& pNeighbors5 , TreeOctNode* node , Pointer( Real ) constraints , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::Integrator& integrator , const Stencil< double , 5 >& lapStencil ) const { - const std::vector< _PointData >& points = pointInfo.points; + const std::vector< PointData >& points = pointInfo.data; if( node->depth()<=_minDepth ) return; bool isInterior; int d , off[3]; @@ -1277,9 +1920,9 @@ void Octree< Real >::UpdateConstraintsFromCoarser( const PointInfo& pointInfo , node->centerIndex( idx ); // Evaluate the current node's basis function at adjacent points for( int x=1 ; x<4 ; x++ ) for( int y=1 ; y<4 ; y++ ) for( int z=1 ; z<4 ; z++ ) - if( neighbors5.neighbors[x][y][z] && pointInfo.pointIndex( neighbors5.neighbors[x][y][z] )!=-1 ) + if( neighbors5.neighbors[x][y][z] && pointInfo.index( neighbors5.neighbors[x][y][z] )!=-1 ) { - const _PointData& pData = points[ pointInfo.pointIndex( neighbors5.neighbors[x][y][z] ) ]; + const PointData& pData = points[ pointInfo.index( neighbors5.neighbors[x][y][z] ) ]; Real weightedPointValue = pData.weightedCoarserValue; Point3D< Real > p = pData.position; constraint += @@ -1399,9 +2042,9 @@ void Octree< Real >::UpSample( int depth , const SortedTreeNodes& sNodes , Const } // At each point @( depth ), evaluate the met solution @( depth-1 ) template< class Real > -void Octree< Real >::SetPointValuesFromCoarser( PointInfo& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) coarseCoefficients ) +void Octree< Real >::SetPointValuesFromCoarser( SparseNodeData< PointData >& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) coarseCoefficients ) { - std::vector< _PointData >& points = pointInfo.points; + std::vector< PointData >& points = pointInfo.data; // For every node at the current depth std::vector< typename TreeOctNode::NeighborKey3 > neighborKeys( std::max< int >( 1 , threads ) ); for( int i=0 ; i::SetPointValuesFromCoarser( PointInfo& pointInfo , int depth for( int i=sNodes.nodeCount[depth] ; i::SetPointValuesFromCoarser( PointInfo& pointInfo , int depth } } template< class Real > -Real Octree< Real >::_WeightedCoarserFunctionValue( const _PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey , const TreeOctNode* pointNode , ConstPointer( Real ) coarseCoefficients ) const +Real Octree< Real >::_WeightedCoarserFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey , const TreeOctNode* pointNode , ConstPointer( Real ) coarseCoefficients ) const { double pointValue = 0; int depth = pointNode->depth(); @@ -1476,9 +2119,9 @@ Real Octree< Real >::_WeightedCoarserFunctionValue( const _PointData& pointData return Real( pointValue * weight ); } template< class Real > -void Octree< Real >::SetPointConstraintsFromFiner( const PointInfo& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) finerCoefficients , Pointer( Real ) coarserConstraints ) const +void Octree< Real >::SetPointConstraintsFromFiner( const SparseNodeData< PointData >& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) finerCoefficients , Pointer( Real ) coarserConstraints ) const { - const std::vector< _PointData >& points = pointInfo.points; + const std::vector< PointData >& points = pointInfo.data; // Note: We can't iterate over the finer point nodes as the point weights might be // scaled incorrectly, due to the adaptive exponent. So instead, we will iterate // over the coarser nodes and evaluate the finer solution at the associated points. @@ -1491,7 +2134,7 @@ void Octree< Real >::SetPointConstraintsFromFiner( const PointInfo& pointInfo , for( int i=sNodes.nodeCount[depth-1] ; i::SetPointConstraintsFromFiner( const PointInfo& pointInfo , } } template< class Real > -Real Octree< Real >::_WeightedFinerFunctionValue( const _PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey , const TreeOctNode* pointNode , ConstPointer( Real ) finerCoefficients ) const +Real Octree< Real >::_WeightedFinerFunctionValue( const PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey , const TreeOctNode* pointNode , ConstPointer( Real ) finerCoefficients ) const { typename TreeOctNode::Neighbors3 childNeighbors; double pointValue = 0; @@ -1565,7 +2208,7 @@ Real Octree< Real >::_WeightedFinerFunctionValue( const _PointData& pointData , return Real( pointValue * weight ); } template< class Real > -int Octree< Real >::GetSliceMatrixAndUpdateConstraints( const PointInfo& pointInfo , SparseMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine , int nStart , int nEnd ) +int Octree< Real >::GetSliceMatrixAndUpdateConstraints( const SparseNodeData< PointData >& pointInfo , SparseMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine , int nStart , int nEnd ) { size_t range = nEnd-nStart; Stencil< double , 5 > stencil , stencils[2][2][2]; @@ -1620,7 +2263,7 @@ int Octree< Real >::GetSliceMatrixAndUpdateConstraints( const PointInfo& pointIn return 1; } template< class Real > -int Octree< Real >::GetMatrixAndUpdateConstraints( const PointInfo& pointInfo , SparseSymmetricMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine ) +int Octree< Real >::GetMatrixAndUpdateConstraints( const SparseNodeData< PointData >& pointInfo , SparseSymmetricMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< 2 >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine ) { size_t start = sNodes.nodeCount[depth] , end = sNodes.nodeCount[depth+1] , range = end-start; Stencil< double , 5 > stencil , stencils[2][2][2]; @@ -1673,7 +2316,7 @@ int Octree< Real >::GetMatrixAndUpdateConstraints( const PointInfo& pointInfo , } template< class Real > -Pointer( Real ) Octree< Real >::SolveSystem( PointInfo& pointInfo , Pointer( Real ) constraints , bool showResidual , int iters , int maxSolveDepth , int cgDepth , double accuracy ) +Pointer( Real ) Octree< Real >::SolveSystem( SparseNodeData< PointData >& pointInfo , Pointer( Real ) constraints , bool showResidual , int iters , int maxSolveDepth , int cgDepth , double accuracy ) { int iter=0; typename BSplineData< 2 >::Integrator integrator; @@ -1729,10 +2372,10 @@ void Octree< Real >::_setMultiColorIndices( int start , int end , std::vector< s } } template< class Real > -int Octree< Real >::_SolveSystemGS( PointInfo& pointInfo , int depth , const typename BSplineData< 2 >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual , double* bNorm2 , double* inRNorm2 , double* outRNorm2 , bool forceSilent ) +int Octree< Real >::_SolveSystemGS( SparseNodeData< PointData >& pointInfo , int depth , const typename BSplineData< 2 >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual , double* bNorm2 , double* inRNorm2 , double* outRNorm2 , bool forceSilent ) { - Pointer( Real ) metSolution = NullPointer< Real >(); - Pointer( Real ) metConstraints = NullPointer< Real >(); + Pointer( Real ) metSolution = NullPointer( Real ); + Pointer( Real ) metConstraints = NullPointer( Real ); if( coarseToFine ) metSolution = metSolutionConstraints; // This stores the up-sampled solution up to depth-2 else metConstraints = metSolutionConstraints; // This stores the down-sampled constraints up to depth @@ -1885,10 +2528,10 @@ int Octree< Real >::_SolveSystemGS( PointInfo& pointInfo , int depth , const typ return iters; } template< class Real > -int Octree< Real >::_SolveSystemCG( PointInfo& pointInfo , int depth , const typename BSplineData< 2 >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual , double* bNorm2 , double* inRNorm2 , double* outRNorm2 , double accuracy ) +int Octree< Real >::_SolveSystemCG( SparseNodeData< PointData >& pointInfo , int depth , const typename BSplineData< 2 >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual , double* bNorm2 , double* inRNorm2 , double* outRNorm2 , double accuracy ) { - Pointer( Real ) metSolution = NullPointer< Real >(); - Pointer( Real ) metConstraints = NullPointer< Real >(); + Pointer( Real ) metSolution = NullPointer( Real ); + Pointer( Real ) metConstraints = NullPointer( Real ); if( coarseToFine ) metSolution = metSolutionConstraints; // This stores the up-sampled solution up to depth-2 else metConstraints = metSolutionConstraints; // This stores the down-sampled constraints up to depth double _maxMemoryUsage = maxMemoryUsage; @@ -1924,8 +2567,8 @@ int Octree< Real >::_SolveSystemCG( PointInfo& pointInfo , int depth , const typ systemTime = Time(); { // Get the system matrix (and adjust the right-hand-side based on the coarser solution if prolonging) - if( coarseToFine ) GetMatrixAndUpdateConstraints( pointInfo , M , constraints , integrator , depth , sNodes , metSolution , true ); - else GetMatrixAndUpdateConstraints( pointInfo , M , constraints , integrator , depth , sNodes , NullPointer< Real >() , false ); + if( coarseToFine ) GetMatrixAndUpdateConstraints( pointInfo , M , constraints , integrator , depth , sNodes , metSolution , true ); + else GetMatrixAndUpdateConstraints( pointInfo , M , constraints , integrator , depth , sNodes , NullPointer( Real ) , false ); // Set the constraint vector B.Resize( sNodes.nodeCount[depth+1]-sNodes.nodeCount[depth] ); for( int i=sNodes.nodeCount[depth] ; i::_SolveSystemCG( PointInfo& pointInfo , int depth , const typ return iter; } template< class Real > -int Octree< Real >::HasNormals( TreeOctNode* node , const NormalInfo& normalInfo ) +int Octree< Real >::HasNormals( TreeOctNode* node , const SparseNodeData< Point3D< Real > >& normalInfo ) { - int idx = normalInfo.normalIndex( node ); + int idx = normalInfo.index( node ); if( idx>=0 ) { - const Point3D< Real >& normal = normalInfo.normals[ idx ]; + const Point3D< Real >& normal = normalInfo.data[ idx ]; if( normal[0]!=0 || normal[1]!=0 || normal[2]!=0 ) return 1; } if( node->children ) for( int i=0 ; ichildren[i] , normalInfo ) ) return 1; return 0; } template< class Real > -Pointer( Real ) Octree< Real >::SetLaplacianConstraints( const NormalInfo& normalInfo ) +Pointer( Real ) Octree< Real >::SetLaplacianConstraints( const SparseNodeData< Point3D< Real > >& normalInfo ) { // To set the Laplacian constraints, we iterate over the // splatted normals and compute the dot-product of the @@ -2062,8 +2705,8 @@ Pointer( Real ) Octree< Real >::SetLaplacianConstraints( const NormalInfo& norma const TreeOctNode* _node = neighbors5.neighbors[x][y][z]; if( _node ) { - int _idx = normalInfo.normalIndex( _node ); - if( _idx>=0 ) constraints[ node->nodeData.nodeIndex ] += Point3D< Real >::Dot( stencil.values[x][y][z] , normalInfo.normals[ _idx ] ); + int _idx = normalInfo.index( _node ); + if( _idx>=0 ) constraints[ node->nodeData.nodeIndex ] += Point3D< Real >::Dot( stencil.values[x][y][z] , normalInfo.data[ _idx ] ); } } else @@ -2072,20 +2715,20 @@ Pointer( Real ) Octree< Real >::SetLaplacianConstraints( const NormalInfo& norma const TreeOctNode* _node = neighbors5.neighbors[x][y][z]; if( _node ) { - int _idx = normalInfo.normalIndex( _node ); + int _idx = normalInfo.index( _node ); if( _idx>=0 ) { int _d , _off[3]; _node->depthAndOffset( _d , _off ); - constraints[ node->nodeData.nodeIndex ] += Real( GetDivergence2( integrator , d , off , _off , false , normalInfo.normals[ _idx ] ) ); + constraints[ node->nodeData.nodeIndex ] += Real( GetDivergence2( integrator , d , off , _off , false , normalInfo.data[ _idx ] ) ); } } } UpdateCoarserSupportBounds( neighbors5.neighbors[2][2][2] , startX , endX , startY , endY , startZ , endZ ); } - int idx = normalInfo.normalIndex( node ); + int idx = normalInfo.index( node ); if( idx<0 ) continue; - const Point3D< Real >& normal = normalInfo.normals[ idx ]; + const Point3D< Real >& normal = normalInfo.data[ idx ]; if( normal[0]==0 && normal[1]==0 && normal[2]==0 ) continue; // Set the constraints for the parents @@ -2133,9 +2776,9 @@ Pointer( Real ) Octree< Real >::SetLaplacianConstraints( const NormalInfo& norma #pragma omp parallel for num_threads( threads ) for( int i=_sNodes.nodeCount[d] ; i<_sNodes.nodeCount[d+1] ; i++ ) { - int idx = normalInfo.normalIndex( _sNodes.treeNodes[i] ); + int idx = normalInfo.index( _sNodes.treeNodes[i] ); if( idx<0 ) continue; - coefficients[i] = normalInfo.normals[ idx ]; + coefficients[i] = normalInfo.data[ idx ]; } } @@ -2210,12 +2853,12 @@ template< class Real > void Octree< Real >::refineBoundary( std::vector< int >* map ){ _sNodes.set( tree , map ); } - template< class Real > -Real Octree< Real >::getCenterValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey , const TreeOctNode* node , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::template CenterEvaluator< 1 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 >& pStencil , bool isInterior ) const +template< class V > +V Octree< Real >::getCenterValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey , const TreeOctNode* node , ConstPointer( V ) solution , ConstPointer( V ) metSolution , const typename BSplineData< 2 >::template CenterEvaluator< 1 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 >& pStencil , bool isInterior ) const { if( node->children ) fprintf( stderr , "[WARNING] getCenterValue assumes leaf node\n" ); - Real value=0; + V value(0); int d , off[3]; node->depthAndOffset( d , off ); @@ -2265,9 +2908,11 @@ Real Octree< Real >::getCenterValue( const typename TreeOctNode::ConstNeighborKe return value; } template< class Real > -Real Octree< Real >::getCornerValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey , const TreeOctNode* node , int corner , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< 2 >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 > stencils[8] , bool isInterior ) const +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 { - double value = 0; + V value(0); +#pragma message ( "[WARNING] Why is this here?" ) if( _boundaryType==-1 ) value = -0.5; int d , off[3]; node->depthAndOffset( d , off ); @@ -2287,7 +2932,7 @@ Real Octree< Real >::getCornerValue( const typename TreeOctNode::ConstNeighborKe for( int x=startX ; xnodeData.nodeIndex ] * stencil.values[x][y][z]; + if( _node ) value += solution[ _node->nodeData.nodeIndex ] * Real( stencil.values[x][y][z] ); } else for( int x=startX ; x::getCornerValue( const typename TreeOctNode::ConstNeighborKe { int _d , _off[3]; _node->depthAndOffset( _d , _off ); - value += solution[ _node->nodeData.nodeIndex ] * evaluator.value( d , off[0] , cx , _off[0] , false , false ) * evaluator.value( d , off[1] , cy , _off[1] , false , false ) * evaluator.value( d , off[2] , cz , _off[2] , false , false ); + value += solution[ _node->nodeData.nodeIndex ] * Real( evaluator.value( d , off[0] , cx , _off[0] , false , false ) * evaluator.value( d , off[1] , cy , _off[1] , false , false ) * evaluator.value( d , off[2] , cz , _off[2] , false , false ) ); } } } @@ -2314,7 +2959,7 @@ Real Octree< Real >::getCornerValue( const typename TreeOctNode::ConstNeighborKe for( int x=startX ; xnodeData.nodeIndex ] * stencils[_corner].values[x][y][z]; + if( _node ) value += metSolution[ _node->nodeData.nodeIndex ] * Real( stencils[_corner].values[x][y][z] ); } else for( int x=startX ; x::getCornerValue( const typename TreeOctNode::ConstNeighborKe { int _d , _off[3]; _node->depthAndOffset( _d , _off ); - value += metSolution[ _node->nodeData.nodeIndex ] * evaluator.value( d , off[0] , cx , _off[0] , false , true ) * evaluator.value( d , off[1] , cy , _off[1] , false , true ) * evaluator.value( d , off[2] , cz , _off[2] , false , true ); + value += metSolution[ _node->nodeData.nodeIndex ] * Real( evaluator.value( d , off[0] , cx , _off[0] , false , true ) * evaluator.value( d , off[1] , cy , _off[1] , false , true ) * evaluator.value( d , off[2] , cz , _off[2] , false , true ) ); } } } @@ -2473,9 +3118,72 @@ Point3D< Real > Octree< Real >::getCornerNormal( const typename TreeOctNode::Con } template< class Real > -Real Octree< Real >::Evaluate( ConstPointer( Real ) coefficients , Point3D< Real > p , const BSplineData< 2 >* fData ) const +template< class V > +V Octree< Real >::_Evaluate( ConstPointer( V ) coefficients , Point3D< Real > p , typename TreeOctNode::ConstNeighborKey3& neighborKey3 ) const { - Real value = Real(0); + V value = V(0); + + for( int d=0 ; d<=neighborKey3.depth() ; d++ ) for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) for( int k=0 ; k<3 ; k++ ) + { + const TreeOctNode* n = neighborKey3.neighbors[d].neighbors[i][j][k]; + if( n ) + { + int dd , off[3]; + n->depthAndOffset( dd , off ); + value += + ( + coefficients[ n->nodeData.nodeIndex ] * + (Real) + ( + _fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[0] ) ]( p[0] ) * + _fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[1] ) ]( p[1] ) * + _fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[2] ) ]( p[2] ) + ) + ); + } + } + return value; +} +template< class Real > +template< class V > +V Octree< Real >::_Evaluate( const SparseNodeData< V >& coefficients , Point3D< Real > p , typename TreeOctNode::ConstNeighborKey3& neighborKey3 ) const +{ + V value = V(0); + + for( int d=0 ; d<=neighborKey3.depth() ; d++ ) for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) for( int k=0 ; k<3 ; k++ ) + { + const TreeOctNode* n = neighborKey3.neighbors[d].neighbors[i][j][k]; + if( n ) + { + int dd , off[3]; + n->depthAndOffset( dd , off ); + + int idx = coefficients.index( n ); + if( idx>=0 ) + value += + ( + coefficients.data[idx] * + (Real) + ( + _fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[0] ) ]( p[0] ) * + _fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[1] ) ]( p[1] ) * + _fData.baseFunctions[ BinaryNode::CenterIndex( dd , off[2] ) ]( p[2] ) + ) + ); + } + } + return value; +} + +template< class Real > +template< class V > +V Octree< Real >::Evaluate( ConstPointer( V ) coefficients , Point3D< Real > p , const BSplineData< 2 >* fData , int depth ) const +{ + int _depth; + if( depth<0 ) _depth = -1; + else _depth = _boundaryType==0 ? depth+1 : depth; + + V value = V(0); BSplineData< 2 > _fData; if( !fData ) _fData.set( tree.maxDepth() , _boundaryType ) , fData = &_fData; const TreeOctNode* n = tree.nextNode(); @@ -2485,27 +3193,70 @@ Real Octree< Real >::Evaluate( ConstPointer( Real ) coefficients , Point3D< Real Real w; n->centerAndWidth( c , w ); c -= p , w *= Real(1.5); - if( fabs(c[0])>w || fabs(c[1])>w || fabs(c[2])>w ) + if( fabs(c[0])>w || fabs(c[1])>w || fabs(c[2])>w || ( _depth>=0 && n->depth()>_depth ) ) { n = tree.nextBranch( n ); continue; } int d , off[3]; n->depthAndOffset( d , off ); - value += (Real) + value += ( coefficients[ n->nodeData.nodeIndex ] * - fData->baseFunctions[ BinaryNode::CenterIndex( d , off[0] ) ]( p[0] ) * - fData->baseFunctions[ BinaryNode::CenterIndex( d , off[1] ) ]( p[1] ) * - fData->baseFunctions[ BinaryNode::CenterIndex( d , off[2] ) ]( p[2] ) + (Real) + ( + fData->baseFunctions[ BinaryNode::CenterIndex( d , off[0] ) ]( p[0] ) * + fData->baseFunctions[ BinaryNode::CenterIndex( d , off[1] ) ]( p[1] ) * + fData->baseFunctions[ BinaryNode::CenterIndex( d , off[2] ) ]( p[2] ) + ) ); n = tree.nextNode( n ); } - if( _boundaryType==-1 ) value -= Real(0.5); return value; } template< class Real > -Pointer( Real ) Octree< Real >::Evaluate( ConstPointer( Real ) coefficients , int& res , Real isoValue , int depth ) +template< class V > +V Octree< Real >::Evaluate( const SparseNodeData< V >& coefficients , Point3D< Real > p , const BSplineData< 2 >* fData , int depth ) const +{ + int _depth; + if( depth<0 ) _depth = -1; + else _depth = _boundaryType==0 ? depth+1 : depth; + V value = V(0); + BSplineData< 2 > _fData; + if( !fData ) _fData.set( tree.maxDepth() , _boundaryType ) , fData = &_fData; + const TreeOctNode* n = tree.nextNode(); + while( n ) + { + Point3D< Real > c; + Real w; + n->centerAndWidth( c , w ); + c -= p , w *= Real(1.5); + if( fabs(c[0])>w || fabs(c[1])>w || fabs(c[2])>w || (_depth>=0 && n->depth()>_depth ) ) n = tree.nextBranch( n ); + else + { + int d , off[3]; + n->depthAndOffset( d , off ); + + int idx = coefficients.index( n ); + if( idx>=0 ) + value += + ( + coefficients.data[idx] * + (Real) + ( + fData->baseFunctions[ BinaryNode::CenterIndex( d , off[0] ) ]( p[0] ) * + fData->baseFunctions[ BinaryNode::CenterIndex( d , off[1] ) ]( p[1] ) * + fData->baseFunctions[ BinaryNode::CenterIndex( d , off[2] ) ]( p[2] ) + ) + ); + n = tree.nextNode( n ); + } + } + return value; +} +template< class Real > +template< class V > +Pointer( V ) Octree< Real >::Evaluate( ConstPointer( V ) coefficients , int& res , Real isoValue , int depth ) { int maxDepth = _boundaryType==0 ? tree.maxDepth()-1 : tree.maxDepth(); if( depth<=0 || depth>maxDepth ) depth = maxDepth; @@ -2539,7 +3290,7 @@ Pointer( Real ) Octree< Real >::Evaluate( ConstPointer( Real ) coefficients , in } } if( skip ) continue; - Real coefficient = coefficients[ n->nodeData.nodeIndex ]; + V coefficient = coefficients[ n->nodeData.nodeIndex ]; for( int x=start[0] ; x<=end[0] ; x+=2 ) for( int y=start[1] ; y<=end[1] ; y+=2 ) for( int z=start[2] ; z<=end[2] ; z+=2 ) @@ -2553,6 +3304,7 @@ Pointer( Real ) Octree< Real >::Evaluate( ConstPointer( Real ) coefficients , in vTables.valueTable[ idx[2] + z*vTables.functionCount ]; } } +#pragma message ( "[WARNING] Why is this here?" ) if( _boundaryType==-1 ) for( int i=0 ; i Neighbors3& setNeighbors( OctNode* root , Point3D< Real > p , int d ); @@ -275,8 +276,9 @@ public: ConstNeighborKey3( void ); ConstNeighborKey3( const ConstNeighborKey3& key3 ); ~ConstNeighborKey3( void ); + int depth( void ) const { return _depth; } - void set(int depth); + void set( int depth ); ConstNeighbors3& getNeighbors( const OctNode* node ); ConstNeighbors3& getNeighbors( const OctNode* node , int minDepth ); void getNeighbors( const OctNode* node , typename OctNode< NodeData >::ConstNeighbors5& neighbors ); diff --git a/Src/Ply.h b/Src/Ply.h index 9645e0e..5a1b8f0 100644 --- a/Src/Ply.h +++ b/Src/Ply.h @@ -234,11 +234,9 @@ typedef struct PlyFace int *vertices; int segment; } PlyFace; - -static char VertexIndicesPropName[] = "vertex_indices"; static PlyProperty face_props[] = { - { VertexIndicesPropName , PLY_INT, PLY_INT, offsetof(PlyFace, vertices), 1, PLY_UCHAR, PLY_UCHAR, offsetof(PlyFace, nr_vertices) }, + { _strdup( "vertex_indices" ) , PLY_INT , PLY_INT , offsetof( PlyFace , vertices ) , 1 , PLY_UCHAR, PLY_UCHAR , offsetof(PlyFace,nr_vertices) }, }; template< class Real > @@ -336,14 +334,17 @@ public: const static int Components=6; static PlyProperty Properties[]; - Point3D point; + Point3D< Real > point; unsigned char color[3]; - operator Point3D& () {return point;} - operator const Point3D& () const {return point;} - PlyColorVertex(void) {point.coords[0]=point.coords[1]=point.coords[2]=0,color[0]=color[1]=color[2]=0;} - PlyColorVertex(const Point3D& p) {point=p;} + operator Point3D< Real >& (){ return point; } + operator const Point3D< Real >& () const { return point; } + PlyColorVertex( void ) { point.coords[0] = point.coords[1] = point.coords[2] = 0 , color[0] = color[1] = color[2] = 0; } + PlyColorVertex( const Point3D& p ) { point=p; } + PlyColorVertex( const Point3D< Real >& p , const unsigned char c[3] ) { point = p , color[0] = c[0] , color[1] = c[1] , color[2] = c[2]; } }; +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[]= { @@ -354,6 +355,86 @@ PlyProperty PlyColorVertex< Real >::Properties[]= { _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 } }; +template< class Real > +class PlyColorAndValueVertex +{ +public: + const static int Components=7; + static PlyProperty Properties[]; + + Point3D< Real > point; + unsigned char color[3]; + Real value; + + operator Point3D< Real >& (){ return point; } + operator const Point3D< Real >& () const { return point; } + PlyColorAndValueVertex( void ) { point.coords[0] = point.coords[1] = point.coords[2] = (Real)0 , color[0] = color[1] = color[2] = 0 , value = (Real)0; } + PlyColorAndValueVertex( const Point3D< Real >& p ) { point=p; } + 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[]= +{ + { _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 } +}; + +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() ); diff --git a/Src/PointStream.h b/Src/PointStream.h index e949623..4d08d2a 100644 --- a/Src/PointStream.h +++ b/Src/PointStream.h @@ -31,54 +31,110 @@ DAMAGE. #include "Ply.h" template< class Real > -class PointStream +struct OrientedPoint3D { -public: - virtual ~PointStream( void ){} - virtual void reset( void ) = 0; - virtual bool nextPoint( Point3D< Real >& p , Point3D< Real >& n ) = 0; + Point3D< Real > p , n; + OrientedPoint3D( Point3D< Real > pp=Point3D< Real >() , Point3D< Real > nn=Point3D< Real >() ) : p(pp) , n(nn) { ; } }; template< class Real > -class MemoryPointStream : public PointStream< Real > +class OrientedPointStream { - const std::pair< Point3D< Real > , Point3D< Real > >* _points; +public: + virtual ~OrientedPointStream( void ){} + virtual void reset( void ) = 0; + virtual bool nextPoint( OrientedPoint3D< Real >& p ) = 0; +}; +template< class Real , class Data > +class OrientedPointStreamWithData : public OrientedPointStream< Real > +{ +public: + virtual ~OrientedPointStreamWithData( void ){} + virtual void reset( void ) = 0; + virtual bool nextPoint( OrientedPoint3D< Real >& p , Data& d ) = 0; + + virtual bool nextPoint( OrientedPoint3D< Real >& p ){ Data d ; return nextPoint( p , d ); } +}; + +template< class Real > +class MemoryOrientedPointStream : public OrientedPointStream< Real > +{ + const OrientedPoint3D< Real >* _points; size_t _pointCount; size_t _current; public: - MemoryPointStream( size_t pointCount , std::pair< Point3D< Real > , Point3D< Real > >* points ); - ~MemoryPointStream( void ); + MemoryOrientedPointStream( size_t pointCount , const OrientedPoint3D< Real >* points ); + ~MemoryOrientedPointStream( void ); void reset( void ); - bool nextPoint( Point3D< Real >& p , Point3D< Real >& n ); + bool nextPoint( OrientedPoint3D< Real >& p ); +}; + +template< class Real , class Data > +class MemoryOrientedPointStreamWithData : public OrientedPointStreamWithData< Real , Data > +{ + const std::pair< OrientedPoint3D< Real > , Data >* _points; + size_t _pointCount; + size_t _current; +public: + MemoryOrientedPointStreamWithData( size_t pointCount , const std::pair< OrientedPoint3D< Real > , Data >* points ); + ~MemoryOrientedPointStreamWithData( void ); + void reset( void ); + bool nextPoint( OrientedPoint3D< Real >& p , Data& d ); }; template< class Real > -class ASCIIPointStream : public PointStream< Real > +class ASCIIOrientedPointStream : public OrientedPointStream< Real > { FILE* _fp; public: - ASCIIPointStream( const char* fileName ); - ~ASCIIPointStream( void ); + ASCIIOrientedPointStream( const char* fileName ); + ~ASCIIOrientedPointStream( void ); void reset( void ); - bool nextPoint( Point3D< Real >& p , Point3D< Real >& n ); + bool nextPoint( OrientedPoint3D< Real >& p ); +}; + +template< class Real , class Data > +class ASCIIOrientedPointStreamWithData : public OrientedPointStreamWithData< Real , Data > +{ + FILE* _fp; + Data (*_readData)( FILE* ); +public: + ASCIIOrientedPointStreamWithData( const char* fileName , Data (*readData)( FILE* ) ); + ~ASCIIOrientedPointStreamWithData( void ); + void reset( void ); + bool nextPoint( OrientedPoint3D< Real >& p , Data& d ); }; template< class Real > -class BinaryPointStream : public PointStream< Real > +class BinaryOrientedPointStream : public OrientedPointStream< Real > { FILE* _fp; static const int POINT_BUFFER_SIZE=1024; - Real _pointBuffer[ POINT_BUFFER_SIZE * 2 * 3 ]; + OrientedPoint3D< Real > _pointBuffer[ POINT_BUFFER_SIZE ]; int _pointsInBuffer , _currentPointIndex; public: - BinaryPointStream( const char* filename ); - ~BinaryPointStream( void ); + BinaryOrientedPointStream( const char* filename ); + ~BinaryOrientedPointStream( void ); void reset( void ); - bool nextPoint( Point3D< Real >& p , Point3D< Real >& n ); + bool nextPoint( OrientedPoint3D< Real >& p ); +}; + +template< class Real , class Data > +class BinaryOrientedPointStreamWithData : public OrientedPointStreamWithData< Real , Data > +{ + FILE* _fp; + static const int POINT_BUFFER_SIZE=1024; + std::pair< OrientedPoint3D< Real > , Data > _pointBuffer[ POINT_BUFFER_SIZE ]; + int _pointsInBuffer , _currentPointIndex; +public: + BinaryOrientedPointStreamWithData( const char* filename ); + ~BinaryOrientedPointStreamWithData( void ); + void reset( void ); + bool nextPoint( OrientedPoint3D< Real >& p , Data& d ); }; template< class Real > -class PLYPointStream : public PointStream< Real > +class PLYOrientedPointStream : public OrientedPointStream< Real > { char* _fileName; PlyFile* _ply; @@ -88,10 +144,32 @@ class PLYPointStream : public PointStream< Real > int _pCount , _pIdx; void _free( void ); public: - PLYPointStream( const char* fileName ); - ~PLYPointStream( void ); + PLYOrientedPointStream( const char* fileName ); + ~PLYOrientedPointStream( void ); void reset( void ); - bool nextPoint( Point3D< Real >& p , Point3D< Real >& n ); + bool nextPoint( OrientedPoint3D< Real >& p ); }; + +template< class Real , class Data > +class PLYOrientedPointStreamWithData : public OrientedPointStreamWithData< Real , Data > +{ + struct _PlyOrientedVertexWithData : public PlyOrientedVertex< Real > { Data data; }; + char* _fileName; + PlyFile* _ply; + int _nr_elems; + char **_elist; + PlyProperty* _dataProperties; + int _dataPropertiesCount; + bool (*_validationFunction)( const bool* ); + + int _pCount , _pIdx; + void _free( void ); +public: + PLYOrientedPointStreamWithData( const char* fileName , const PlyProperty* dataProperties , int dataPropertiesCount , bool (*validationFunction)( const bool* )=NULL ); + ~PLYOrientedPointStreamWithData( void ); + void reset( void ); + bool nextPoint( OrientedPoint3D< Real >& p , Data& d ); +}; + #include "PointStream.inl" #endif // POINT_STREAM_INCLUDED diff --git a/Src/PointStream.inl b/Src/PointStream.inl index 0a312d5..f59685c 100644 --- a/Src/PointStream.inl +++ b/Src/PointStream.inl @@ -26,88 +26,97 @@ ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF S DAMAGE. */ + +/////////////////////////////// +// MemoryOrientedPointStream // +/////////////////////////////// template< class Real > -MemoryPointStream< Real >::MemoryPointStream( size_t pointCount , std::pair< Point3D< Real > , Point3D< Real > >* points ){ _points = points , _pointCount = pointCount , _current = 0; } +MemoryOrientedPointStream< Real >::MemoryOrientedPointStream( size_t pointCount , const OrientedPoint3D< Real >* points ){ _points = points , _pointCount = pointCount , _current = 0; } template< class Real > -MemoryPointStream< Real >::~MemoryPointStream( void ){ ; } +MemoryOrientedPointStream< Real >::~MemoryOrientedPointStream( void ){ ; } template< class Real > -void MemoryPointStream< Real >::reset( void ) { _current=0; } +void MemoryOrientedPointStream< Real >::reset( void ) { _current=0; } template< class Real > -bool MemoryPointStream< Real >::nextPoint( Point3D< Real >& p , Point3D< Real >& n ) +bool MemoryOrientedPointStream< Real >::nextPoint( OrientedPoint3D< Real >& p ) { if( _current>=_pointCount ) return false; - p = _points[_current].first , n = _points[_current].second; + p = _points[_current]; _current++; return true; } +////////////////////////////// +// ASCIIOrientedPointStream // +////////////////////////////// template< class Real > -ASCIIPointStream< Real >::ASCIIPointStream( const char* fileName ) +ASCIIOrientedPointStream< Real >::ASCIIOrientedPointStream( const char* fileName ) { _fp = fopen( fileName , "r" ); if( !_fp ) fprintf( stderr , "Failed to open file for reading: %s\n" , fileName ) , exit( 0 ); } template< class Real > -ASCIIPointStream< Real >::~ASCIIPointStream( void ) +ASCIIOrientedPointStream< Real >::~ASCIIOrientedPointStream( void ) { fclose( _fp ); _fp = NULL; } template< class Real > -void ASCIIPointStream< Real >::reset( void ) { fseek( _fp , SEEK_SET , 0 ); } +void ASCIIOrientedPointStream< Real >::reset( void ) { fseek( _fp , SEEK_SET , 0 ); } template< class Real > -bool ASCIIPointStream< Real >::nextPoint( Point3D< Real >& p , Point3D< Real >& n ) +bool ASCIIOrientedPointStream< Real >::nextPoint( OrientedPoint3D< Real >& p ) { float c[2*DIMENSION]; if( fscanf( _fp , " %f %f %f %f %f %f " , &c[0] , &c[1] , &c[2] , &c[3] , &c[4] , &c[5] )!=2*DIMENSION ) return false; - p[0] = c[0] , p[1] = c[1] , p[2] = c[2]; - n[0] = c[3] , n[1] = c[4] , n[2] = c[5]; + p.p[0] = c[0] , p.p[1] = c[1] , p.p[2] = c[2]; + p.n[0] = c[3] , p.n[1] = c[4] , p.n[2] = c[5]; return true; } + +/////////////////////////////// +// BinaryOrientedPointStream // +/////////////////////////////// template< class Real > -BinaryPointStream< Real >::BinaryPointStream( const char* fileName ) +BinaryOrientedPointStream< Real >::BinaryOrientedPointStream( const char* fileName ) { _pointsInBuffer = _currentPointIndex = 0; _fp = fopen( fileName , "rb" ); if( !_fp ) fprintf( stderr , "Failed to open file for reading: %s\n" , fileName ) , exit( 0 ); } template< class Real > -BinaryPointStream< Real >::~BinaryPointStream( void ) +BinaryOrientedPointStream< Real >::~BinaryOrientedPointStream( void ) { fclose( _fp ); _fp = NULL; } template< class Real > -void BinaryPointStream< Real >::reset( void ) +void BinaryOrientedPointStream< Real >::reset( void ) { fseek( _fp , SEEK_SET , 0 ); _pointsInBuffer = _currentPointIndex = 0; } template< class Real > -bool BinaryPointStream< Real >::nextPoint( Point3D< Real >& p , Point3D< Real >& n ) +bool BinaryOrientedPointStream< Real >::nextPoint( OrientedPoint3D< Real >& p ) { if( _currentPointIndex<_pointsInBuffer ) { - p[0] = _pointBuffer[ _currentPointIndex*6+0 ]; - p[1] = _pointBuffer[ _currentPointIndex*6+1 ]; - p[2] = _pointBuffer[ _currentPointIndex*6+2 ]; - n[0] = _pointBuffer[ _currentPointIndex*6+3 ]; - n[1] = _pointBuffer[ _currentPointIndex*6+4 ]; - n[2] = _pointBuffer[ _currentPointIndex*6+5 ]; + p = _pointBuffer[ _currentPointIndex ]; _currentPointIndex++; return true; } else { _currentPointIndex = 0; - _pointsInBuffer = int( fread( _pointBuffer , sizeof( Real ) * 6 , POINT_BUFFER_SIZE , _fp ) ); + _pointsInBuffer = int( fread( _pointBuffer , sizeof( OrientedPoint3D< Real > ) , POINT_BUFFER_SIZE , _fp ) ); if( !_pointsInBuffer ) return false; - else return nextPoint( p , n ); + else return nextPoint( p ); } } +//////////////////////////// +// PLYOrientedPointStream // +//////////////////////////// template< class Real > -PLYPointStream< Real >::PLYPointStream( const char* fileName ) +PLYOrientedPointStream< Real >::PLYOrientedPointStream( const char* fileName ) { _fileName = new char[ strlen( fileName )+1 ]; strcpy( _fileName , fileName ); @@ -115,7 +124,7 @@ PLYPointStream< Real >::PLYPointStream( const char* fileName ) reset(); } template< class Real > -void PLYPointStream< Real >::reset( void ) +void PLYOrientedPointStream< Real >::reset( void ) { int fileType; float version; @@ -166,7 +175,7 @@ void PLYPointStream< Real >::reset( void ) } } template< class Real > -void PLYPointStream< Real >::_free( void ) +void PLYOrientedPointStream< Real >::_free( void ) { if( _ply ) ply_close( _ply ) , _ply = NULL; if( _elist ) @@ -176,20 +185,222 @@ void PLYPointStream< Real >::_free( void ) } } template< class Real > -PLYPointStream< Real >::~PLYPointStream( void ) +PLYOrientedPointStream< Real >::~PLYOrientedPointStream( void ) { _free(); if( _fileName ) delete[] _fileName , _fileName = NULL; } template< class Real > -bool PLYPointStream< Real >::nextPoint( Point3D< Real >& p , Point3D< Real >& n ) +bool PLYOrientedPointStream< Real >::nextPoint( OrientedPoint3D< Real >& p ) { if( _pIdx<_pCount ) { PlyOrientedVertex< Real > op; ply_get_element( _ply, (void *)&op ); - p = op.point; - n = op.normal; + p.p = op.point; + p.n = op.normal; + _pIdx++; + return true; + } + else return false; +} + +/////////////////////////////////////// +// MemoryOrientedPointStreamWithData // +/////////////////////////////////////// +template< class Real , class Data > +MemoryOrientedPointStreamWithData< Real , Data >::MemoryOrientedPointStreamWithData( size_t pointCount , const std::pair< OrientedPoint3D< Real > , Data >* points ){ _points = points , _pointCount = pointCount , _current = 0; } +template< class Real , class Data > +MemoryOrientedPointStreamWithData< Real , Data >::~MemoryOrientedPointStreamWithData( void ){ ; } +template< class Real , class Data > +void MemoryOrientedPointStreamWithData< Real , Data >::reset( void ) { _current=0; } +template< class Real , class Data > +bool MemoryOrientedPointStreamWithData< Real , Data >::nextPoint( OrientedPoint3D< Real >& p , Data& d ) +{ + if( _current>=_pointCount ) return false; + p = _points[_current].first; + d = _points[_current].second; + _current++; + return true; +} + +////////////////////////////////////// +// ASCIIOrientedPointStreamWithData // +////////////////////////////////////// +template< class Real , class Data > +ASCIIOrientedPointStreamWithData< Real , Data >::ASCIIOrientedPointStreamWithData( const char* fileName , Data (*readData)( FILE* ) ) : _readData( readData ) +{ + _fp = fopen( fileName , "r" ); + if( !_fp ) fprintf( stderr , "Failed to open file for reading: %s\n" , fileName ) , exit( 0 ); +} +template< class Real , class Data > +ASCIIOrientedPointStreamWithData< Real , Data >::~ASCIIOrientedPointStreamWithData( void ) +{ + fclose( _fp ); + _fp = NULL; +} +template< class Real , class Data > +void ASCIIOrientedPointStreamWithData< Real , Data >::reset( void ) { fseek( _fp , SEEK_SET , 0 ); } +template< class Real , class Data > +bool ASCIIOrientedPointStreamWithData< Real , Data >::nextPoint( OrientedPoint3D< Real >& p , Data& d ) +{ + float c[2*DIMENSION]; + if( fscanf( _fp , " %f %f %f %f %f %f " , &c[0] , &c[1] , &c[2] , &c[3] , &c[4] , &c[5] )!=2*DIMENSION ) return false; + p.p[0] = c[0] , p.p[1] = c[1] , p.p[2] = c[2]; + p.n[0] = c[3] , p.n[1] = c[4] , p.n[2] = c[5]; + d = _readData( _fp ); + return true; +} + +/////////////////////////////////////// +// BinaryOrientedPointStreamWithData // +/////////////////////////////////////// +template< class Real , class Data > +BinaryOrientedPointStreamWithData< Real , Data >::BinaryOrientedPointStreamWithData( const char* fileName ) +{ + _pointsInBuffer = _currentPointIndex = 0; + _fp = fopen( fileName , "rb" ); + if( !_fp ) fprintf( stderr , "Failed to open file for reading: %s\n" , fileName ) , exit( 0 ); +} +template< class Real , class Data > +BinaryOrientedPointStreamWithData< Real , Data >::~BinaryOrientedPointStreamWithData( void ) +{ + fclose( _fp ); + _fp = NULL; +} +template< class Real , class Data > +void BinaryOrientedPointStreamWithData< Real , Data >::reset( void ) +{ + fseek( _fp , SEEK_SET , 0 ); + _pointsInBuffer = _currentPointIndex = 0; +} +template< class Real , class Data > +bool BinaryOrientedPointStreamWithData< Real , Data >::nextPoint( OrientedPoint3D< Real >& p , Data& d ) +{ + if( _currentPointIndex<_pointsInBuffer ) + { + p = _pointBuffer[ _currentPointIndex ].first; + d = _pointBuffer[ _currentPointIndex ].second; + _currentPointIndex++; + return true; + } + else + { + _currentPointIndex = 0; + _pointsInBuffer = int( fread( _pointBuffer , sizeof( std::pair< OrientedPoint3D< Real > , Data > ) , POINT_BUFFER_SIZE , _fp ) ); + if( !_pointsInBuffer ) return false; + else return nextPoint( p , d ); + } +} + +//////////////////////////////////// +// PLYOrientedPointStreamWithData // +//////////////////////////////////// +template< class Real , class Data > +PLYOrientedPointStreamWithData< Real , Data >::PLYOrientedPointStreamWithData( const char* fileName , const PlyProperty* dataProperties , int dataPropertiesCount , bool (*validationFunction)( const bool* ) ) : _dataPropertiesCount( dataPropertiesCount ) , _validationFunction( validationFunction ) +{ + _dataProperties = new PlyProperty[ _dataPropertiesCount ]; + memcpy( _dataProperties , dataProperties , sizeof(PlyProperty) * _dataPropertiesCount ); + for( int i=0 ; i<_dataPropertiesCount ; i++ ) _dataProperties[i].offset += sizeof( PlyOrientedVertex< Real > ); + _fileName = new char[ strlen( fileName )+1 ]; + strcpy( _fileName , fileName ); + _ply = NULL; + reset(); +} +template< class Real , class Data > +void PLYOrientedPointStreamWithData< Real , Data >::reset( void ) +{ + int fileType; + float version; + PlyProperty** plist; + if( _ply ) _free(); + _ply = ply_open_for_reading( _fileName, &_nr_elems, &_elist, &fileType, &version ); + if( !_ply ) + { + fprintf( stderr, "[ERROR] Failed to open ply file for reading: %s\n" , _fileName ); + exit( 0 ); + } + bool foundVertices = false; + for( int i=0 ; i<_nr_elems ; i++ ) + { + int num_elems; + int nr_props; + char* elem_name = _elist[i]; + plist = ply_get_element_description( _ply , elem_name , &num_elems , &nr_props ); + if( !plist ) + { + fprintf( stderr , "[ERROR] Failed to get element description: %s\n" , elem_name ); + exit( 0 ); + } + + if( equal_strings( "vertex" , elem_name ) ) + { + foundVertices = true; + _pCount = num_elems , _pIdx = 0; + for( int i=0 ; i::Components ; i++ ) + if( !ply_get_property( _ply , elem_name , &(PlyOrientedVertex< Real >::Properties[i]) ) ) + { + fprintf( stderr , "[ERROR] Failed to find property in ply file: %s\n" , PlyOrientedVertex< Real >::Properties[i].name ); + exit( 0 ); + } + if( _validationFunction ) + { + bool* properties = new bool[_dataPropertiesCount]; + for( int i=0 ; i<_dataPropertiesCount ; i++ ) + if( !ply_get_property( _ply , elem_name , &(_dataProperties[i]) ) ) properties[i] = false; + else properties[i] = true; + bool valid = _validationFunction( properties ); + delete[] properties; + if( !valid ) fprintf( stderr , "[ERROR] Failed to validate properties in file\n" ) , exit( 0 ); + } + else + { + for( int i=0 ; i<_dataPropertiesCount ; i++ ) + if( !ply_get_property( _ply , elem_name , &(_dataProperties[i]) ) ) + fprintf( stderr , "[WARNING] Failed to find property in ply file: %s\n" , _dataProperties[i].name ); + } + } + for( int j=0 ; jname ); + free( plist[j] ); + } + free( plist ); + if( foundVertices ) break; + } + if( !foundVertices ) + { + fprintf( stderr , "[ERROR] Could not find vertices in ply file\n" ); + exit( 0 ); + } +} +template< class Real , class Data > +void PLYOrientedPointStreamWithData< Real , Data >::_free( void ) +{ + if( _ply ) ply_close( _ply ) , _ply = NULL; + if( _elist ) + { + for( int i=0 ; i<_nr_elems ; i++ ) free( _elist[i] ); + free( _elist ); + } +} +template< class Real , class Data > +PLYOrientedPointStreamWithData< Real , Data >::~PLYOrientedPointStreamWithData( void ) +{ + _free(); + if( _fileName ) delete[] _fileName , _fileName = NULL; + if( _dataProperties ) delete[] _dataProperties , _dataProperties = NULL; +} +template< class Real , class Data > +bool PLYOrientedPointStreamWithData< Real , Data >::nextPoint( OrientedPoint3D< Real >& p , Data& d ) +{ + if( _pIdx<_pCount ) + { + _PlyOrientedVertexWithData op; + ply_get_element( _ply, (void *)&op ); + p.p = op.point; + p.n = op.normal; + d = op.data; _pIdx++; return true; } diff --git a/Src/PoissonRecon.cpp b/Src/PoissonRecon.cpp index 30eb284..923c4d8 100644 --- a/Src/PoissonRecon.cpp +++ b/Src/PoissonRecon.cpp @@ -46,8 +46,8 @@ DAMAGE. #include "omp.h" #endif // _OPENMP void DumpOutput( const char* format , ... ); -void DumpOutput2( char* str , const char* format , ... ); #include "MultiGridOctreeData.h" +void DumpOutput2( std::vector< char* >& comments , const char* format , ... ); #define DEFAULT_FULL_DEPTH 5 @@ -79,7 +79,7 @@ void DumpOutput( const char* format , ... ) va_end( args ); } } -void DumpOutput2( char* str , const char* format , ... ) +void DumpOutput2( std::vector< char* >& comments , const char* format , ... ) { if( outputFile ) { @@ -97,6 +97,8 @@ void DumpOutput2( char* str , const char* 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 ); @@ -141,7 +143,8 @@ cmdLineInt Threads( "threads" , omp_get_num_procs() ); cmdLineFloat - SamplesPerNode( "samplesPerNode" , 1.f ) , + Color( "color" , 16.f ) , + SamplesPerNode( "samplesPerNode" , 1.5f ) , Scale( "scale" , 1.1f ) , CSSolverAccuracy( "cgAccuracy" , float(1e-3) ) , PointWeight( "pointWeight" , 4.f ); @@ -159,6 +162,7 @@ cmdLineReadable* params[] = &MinDepth , &CGDepth , &Iters , &Complete , + &Color , #ifdef _WIN32 &Performance , #endif // _WIN32 @@ -200,6 +204,9 @@ void ShowUsage(char* ex) printf( "\t[--%s =%d]\n" , Iters.name , Iters.value ); printf( "\t\t This flag specifies the (maximum if CG) number of solver iterations.\n" ); + printf( "\t[--%s ]\n" , Color.name ); + printf( "\t\t This flag specifies the pull factor for color interpolation\n" ); + #ifdef _OPENMP printf( "\t[--%s =%d]\n" , Threads.name , Threads.value ); printf( "\t\t This parameter specifies the number of threads across which\n" ); @@ -260,17 +267,32 @@ void ShowUsage(char* ex) 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" ); } -template< class Real , class Vertex > + +Point3D< unsigned char > ReadASCIIColor( FILE* fp ) +{ + Point3D< unsigned char > c; + if( fscanf( fp , " %c %c %c " , &c[0] , &c[1] , &c[2] )!=3 ) fprintf( stderr , "[ERROR] Failed to read color\n" ) , exit( 0 ); + return c; +} + +PlyProperty PlyColorProperties[]= +{ + { "r" , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[0] ) ) , 0 , 0 , 0 , 0 } , + { "g" , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[1] ) ) , 0 , 0 , 0 , 0 } , + { "b" , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[2] ) ) , 0 , 0 , 0 , 0 } , + { "red" , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[0] ) ) , 0 , 0 , 0 , 0 } , + { "green" , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , coords[1] ) ) , 0 , 0 , 0 , 0 } , + { "blue" , PLY_UCHAR , PLY_UCHAR , int( offsetof( Point3D< unsigned char > , 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 , class _Vertex > int Execute( int argc , char* argv[] ) { Reset< Real >(); - int i; int paramNum = sizeof(params)/sizeof(cmdLineReadable*); - int commentNum=0; - char **comments; - - comments = new char*[paramNum+7]; - for( i=0 ; i comments; if( Verbose.set ) echoStdout=1; @@ -288,7 +310,7 @@ int Execute( int argc , char* argv[] ) for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ ) { float f; - fscanf( fp , " %f " , &f ); + if( fscanf( fp , " %f " , &f )!=1 ) fprintf( stderr , "[ERROR] Execute: Failed to read xform\n" ) , exit( 0 ); xForm(i,j) = (Real)f; } fclose( fp ); @@ -297,14 +319,14 @@ int Execute( int argc , char* argv[] ) else xForm = XForm4x4< Real >::Identity(); iXForm = xForm.inverse(); - DumpOutput2( comments[commentNum++] , "Running Screened Poisson Reconstruction (Version 6.13a)\n" ); + DumpOutput2( comments , "Running Screened Poisson Reconstruction (Version 7.0)\n" ); char str[1024]; for( int i=0 ; iset ) { params[i]->writeValue( str ); - if( strlen( str ) ) DumpOutput2( comments[commentNum++] , "\t--%s %s\n" , params[i]->name , str ); - else DumpOutput2( comments[commentNum++] , "\t--%s\n" , params[i]->name ); + if( strlen( str ) ) DumpOutput2( comments , "\t--%s %s\n" , params[i]->name , str ); + else DumpOutput2( comments , "\t--%s\n" , params[i]->name ); } double t; @@ -332,20 +354,43 @@ int Execute( int argc , char* argv[] ) double maxMemoryUsage; t=Time() , tree.maxMemoryUsage=0; - typename Octree< Real >::PointInfo* pointInfo = new typename Octree< Real >::PointInfo(); - typename Octree< Real >::NormalInfo* normalInfo = new typename Octree< Real >::NormalInfo(); + typename Octree< Real >::template SparseNodeData< typename Octree< Real >::PointData >* pointInfo = new typename Octree< Real >::template SparseNodeData< typename Octree< Real >::PointData >(); + typename Octree< Real >::template SparseNodeData< Point3D< Real > >* normalInfo = new typename Octree< Real >::template SparseNodeData< Point3D< Real > >(); std::vector< Real >* kernelDensityWeights = new std::vector< Real >(); std::vector< Real >* centerWeights = new std::vector< Real >(); - PointStream< float >* pointStream; + int pointCount; + typedef typename Octree< Real >::template ProjectiveData< Point3D< Real > > ProjectiveColor; + typename Octree< Real >::template SparseNodeData< ProjectiveColor > colorData; + char* ext = GetFileExtension( In.value ); - if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryPointStream< float >( In.value ); - else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYPointStream< float >( In.value ); - else pointStream = new ASCIIPointStream< float >( In.value ); + if( Color.set && Color.value>0 ) + { + OrientedPointStreamWithData< float , Point3D< unsigned char > >* pointStream; + if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStreamWithData< float , Point3D< unsigned char > >( In.value ); + else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStreamWithData< float , Point3D< unsigned char > >( In.value , PlyColorProperties , 6 , ValidPlyColorProperties ); + else pointStream = new ASCIIOrientedPointStreamWithData< float , Point3D< unsigned char > >( In.value , ReadASCIIColor ); + pointCount = tree.template SetTree< float >( pointStream , MinDepth.value , Depth.value , FullDepth.value , kernelDepth , Real(SamplesPerNode.value) , Scale.value , Confidence.set , NormalWeights.set , PointWeight.value , AdaptiveExponent.value , *kernelDensityWeights , *pointInfo , *normalInfo , *centerWeights , colorData , xForm , BoundaryType.value , Complete.set ); + delete pointStream; + + 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( Color.value , n->depth() ); + } + } + else + { + OrientedPointStream< float >* pointStream; + if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStream< float >( In.value ); + else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStream< float >( In.value ); + else pointStream = new ASCIIOrientedPointStream< float >( In.value ); + pointCount = tree.template SetTree< float >( pointStream , MinDepth.value , Depth.value , FullDepth.value , kernelDepth , Real(SamplesPerNode.value) , Scale.value , Confidence.set , NormalWeights.set , PointWeight.value , AdaptiveExponent.value , *kernelDensityWeights , *pointInfo , *normalInfo , *centerWeights , xForm , BoundaryType.value , Complete.set ); + delete pointStream; + } delete[] ext; - int pointCount = tree.template SetTree< float >( pointStream , MinDepth.value , Depth.value , FullDepth.value , kernelDepth , Real(SamplesPerNode.value) , Scale.value , Confidence.set , NormalWeights.set , PointWeight.value , AdaptiveExponent.value , *pointInfo , *normalInfo , *kernelDensityWeights , *centerWeights , BoundaryType.value , xForm , Complete.set ); if( !Density.set ) delete kernelDensityWeights , kernelDensityWeights = NULL; - DumpOutput2( comments[commentNum++] , "# Tree set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage ); + DumpOutput2( comments , "# Tree set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage ); DumpOutput( "Input Points: %d\n" , pointCount ); DumpOutput( "Leaves/Nodes: %d/%d\n" , tree.tree.leaves() , tree.tree.nodes() ); DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) ); @@ -354,7 +399,7 @@ int Execute( int argc , char* argv[] ) t=Time() , tree.maxMemoryUsage=0; Pointer( Real ) constraints = tree.SetLaplacianConstraints( *normalInfo ); delete normalInfo; - DumpOutput2( comments[commentNum++] , "# Constraints set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage ); + DumpOutput2( comments , "# Constraints set 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 ); @@ -362,7 +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[commentNum++] , "# Linear system solved in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage ); + 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 ); @@ -402,11 +447,11 @@ int Execute( int argc , char* argv[] ) if( Out.set ) { t = Time() , tree.maxMemoryUsage = 0; - tree.GetMCIsoSurface( kernelDensityWeights ? GetPointer( *kernelDensityWeights ) : NullPointer< Real >() , solution , isoValue , mesh , true , !NonManifold.set , PolygonMesh.set ); - if( PolygonMesh.set ) DumpOutput2( comments[commentNum++] , "# Got polygons in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage ); - else DumpOutput2( comments[commentNum++] , "# Got triangles in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage ); + tree.template GetMCIsoSurface< Vertex , _Vertex >( 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 ); - DumpOutput2( comments[commentNum++],"# Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-tt , maxMemoryUsage ); + DumpOutput2( comments , "# Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-tt , maxMemoryUsage ); if( NoComments.set ) { @@ -415,8 +460,8 @@ int Execute( int argc , char* argv[] ) } else { - if( ASCII.set ) PlyWritePolygons( Out.value , &mesh , PLY_ASCII , comments , commentNum , iXForm ); - else PlyWritePolygons( Out.value , &mesh , PLY_BINARY_NATIVE , comments , commentNum , iXForm ); + if( ASCII.set ) PlyWritePolygons( Out.value , &mesh , PLY_ASCII , &comments[0] , (int)comments.size() , iXForm ); + else PlyWritePolygons( Out.value , &mesh , PLY_BINARY_NATIVE , &comments[0] , (int)comments.size() , iXForm ); } DumpOutput( "Vertices / Polygons: %d / %d\n" , mesh.outOfCorePointCount()+mesh.inCorePoints.size() , mesh.polygonCount() ); } @@ -456,11 +501,19 @@ int main( int argc , char* argv[] ) cmdLineParse( argc-1 , &argv[1] , sizeof(params)/sizeof(cmdLineReadable*) , params , 1 ); if( Density.set ) - if( Double.set ) Execute< double , PlyValueVertex< float > >( argc , argv ); - else Execute< float , PlyValueVertex< float > >( argc , argv ); + if( Color.set ) + if( Double.set ) Execute< double , PlyColorAndValueVertex< float > , _PlyColorAndValueVertex< float > >( argc , argv ); + else Execute< float , PlyColorAndValueVertex< 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 ); else - if( Double.set ) Execute< double , PlyVertex< float > >( argc , argv ); - else Execute< float , PlyVertex< float > >( argc , argv ); + if( Color.set ) + if( Double.set ) Execute< double , PlyColorVertex< float > , _PlyColorVertex< float > >( argc , argv ); + else Execute< float , PlyColorVertex< 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 ); #ifdef _WIN32 if( Performance.set ) { diff --git a/Src/SparseMatrix.inl b/Src/SparseMatrix.inl index e8bf71e..4decf5d 100644 --- a/Src/SparseMatrix.inl +++ b/Src/SparseMatrix.inl @@ -43,8 +43,8 @@ void SparseMatrix< T >::_init( void ) _contiguous = false; _maxEntriesPerRow = 0; rows = 0; - rowSizes = NullPointer< int >( ); - m_ppElements = NullPointer< Pointer( MatrixEntry< T > ) >( ); + rowSizes = NullPointer( int ); + m_ppElements = NullPointer( Pointer( MatrixEntry< T > ) ); } template< class T > SparseMatrix< T >::SparseMatrix( void ){ _init(); } diff --git a/Src/Vector.inl b/Src/Vector.inl index db06a3a..0b042fb 100644 --- a/Src/Vector.inl +++ b/Src/Vector.inl @@ -38,13 +38,13 @@ template Vector::Vector( void ) { m_N = 0; - m_pV = NullPointer< T >(); + m_pV = NullPointer( T ); } template< class T > Vector< T >::Vector( const Vector& V ) { m_N = 0; - m_pV = NullPointer< T >(); + m_pV = NullPointer( T ); Resize( V.m_N ); memcpy( m_pV , V.m_pV , m_N*sizeof(T) ); } @@ -52,7 +52,7 @@ template Vector::Vector( size_t N ) { m_N=0; - m_pV = NullPointer< T >(); + m_pV = NullPointer( T ); Resize(N); } template diff --git a/SurfaceTrimmer.vcxproj.filters b/SurfaceTrimmer.vcxproj.filters new file mode 100644 index 0000000..b6ea680 --- /dev/null +++ b/SurfaceTrimmer.vcxproj.filters @@ -0,0 +1,47 @@ + + + + + Source Files + + + Source Files + + + Source Files + + + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + + + {237cb93d-2caa-4f6b-9bb3-5d597e5fd1dd} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {e857a52e-da0b-4eea-9e27-7b9c21e142c9} + h;hpp;hxx;hm;inl;inc;xsd + + + {61c77966-1174-4dae-9079-00d39bbe4aab} + inc;inl + + + + + Include Files + + + \ No newline at end of file