From 0952160ee6ce84e76b4fec36ae5ce113f56da0ac Mon Sep 17 00:00:00 2001 From: Harsh Bhatia Date: Wed, 21 Jul 2021 17:05:50 -0700 Subject: [PATCH 1/5] added a cmake file for PoissonRecon target --- CMakeLists.txt | 54 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 00000000..f5f27e20 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,54 @@ +## ----------------------------------------------------------------------------- +cmake_minimum_required(VERSION 3.20) +project(AdaptiveMultigridSolvers) + + +set(CMAKE_VERBOSE_MAKEFILE ON) + +## ----------------------------------------------------------------------------- +set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -std=c++14 -pthread -Wno-invalid-offsetof -fPIC -Wno-deprecated") +set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${CMAKE_CXX_FLAGS} -g") +set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${CMAKE_CXX_FLAGS} -Ofast -DRELEASE -funroll-loops -ffast-math") + + +## ----------------------------------------------------------------------------- +find_package(OpenMP) + + +## ----------------------------------------------------------------------------- +# all include paths +include_directories(${PROJECT_SOURCE_DIR}) +include_directories(${PROJECT_SOURCE_DIR}/Src) +include_directories(${PROJECT_SOURCE_DIR}/PNG) +include_directories(${PROJECT_SOURCE_DIR}/JPEG) +include_directories(${PROJECT_SOURCE_DIR}/ZLIB) + +## ----------------------------------------------------------------------------- +## build dependencies as libraries +## ----------------------------------------------------------------------------- +set(TRG_PNG pr_png) +set(TRG_JPEG pr_jpeg) +set(TRG_ZLIB pr_zlib) + +file(GLOB_RECURSE SRCS_PNG "PNG/*.c") +file(GLOB_RECURSE SRCS_JPEG "JPEG/*.cpp") +file(GLOB_RECURSE SRCS_ZLIB "ZLIB/*.c") + +add_library(${TRG_PNG} STATIC ${SRCS_PNG}) +add_library(${TRG_JPEG} SHARED ${SRCS_JPEG}) +add_library(${TRG_ZLIB} SHARED ${SRCS_ZLIB}) + + +## ----------------------------------------------------------------------------- +## build poisson reconstruction +## ----------------------------------------------------------------------------- +add_executable(poisson_recon + ${PROJECT_SOURCE_DIR}/Src/PoissonRecon.cpp) + +target_link_libraries(poisson_recon + ${TRG_PNG} ${TRG_JPEG} ${TRG_ZLIB}) + +## ----------------------------------------------------------------------------- From 1a8d5717b53db2a4a1bb7799ffbfb24c5386e5bf Mon Sep 17 00:00:00 2001 From: Harsh Bhatia Date: Wed, 21 Jul 2021 17:09:00 -0700 Subject: [PATCH 2/5] split PoissonRecon target into a shlib and exe --- CMakeLists.txt | 11 +- Src/PoissonRecon.cpp | 1139 +------------------------------------- Src/PoissonReconLib.cpp | 1166 +++++++++++++++++++++++++++++++++++++++ Src/PoissonReconLib.h | 31 ++ 4 files changed, 1209 insertions(+), 1138 deletions(-) create mode 100644 Src/PoissonReconLib.cpp create mode 100644 Src/PoissonReconLib.h diff --git a/CMakeLists.txt b/CMakeLists.txt index f5f27e20..fa694ae4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -43,12 +43,19 @@ add_library(${TRG_ZLIB} SHARED ${SRCS_ZLIB}) ## ----------------------------------------------------------------------------- -## build poisson reconstruction +## build poisson reconstruction library and executable ## ----------------------------------------------------------------------------- +add_library(poisson_recon_lib SHARED + ${PROJECT_SOURCE_DIR}/Src/PoissonReconLib.h + ${PROJECT_SOURCE_DIR}/Src/PoissonReconLib.cpp) + +target_link_libraries(poisson_recon_lib + ${TRG_PNG} ${TRG_JPEG} ${TRG_ZLIB}) + add_executable(poisson_recon ${PROJECT_SOURCE_DIR}/Src/PoissonRecon.cpp) target_link_libraries(poisson_recon - ${TRG_PNG} ${TRG_JPEG} ${TRG_ZLIB}) + poisson_recon_lib ${TRG_PNG} ${TRG_JPEG} ${TRG_ZLIB}) ## ----------------------------------------------------------------------------- diff --git a/Src/PoissonRecon.cpp b/Src/PoissonRecon.cpp index 5659b083..0eb541a6 100644 --- a/Src/PoissonRecon.cpp +++ b/Src/PoissonRecon.cpp @@ -26,1140 +26,7 @@ ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF S DAMAGE. */ -#include "PreProcessor.h" - -#undef USE_DOUBLE // If enabled, double-precesion is used - -#define DATA_DEGREE 0 // The order of the B-Spline used to splat in data for color interpolation -#define WEIGHT_DEGREE 2 // The order of the B-Spline used to splat in the weights for density estimation -#define NORMAL_DEGREE 2 // The order of the B-Spline used to splat in the normals for constructing the Laplacian constraints -#define DEFAULT_FEM_DEGREE 1 // The default finite-element degree -#define DEFAULT_FEM_BOUNDARY BOUNDARY_NEUMANN // The default finite-element boundary type -#define DEFAULT_DIMENSION 3 // The dimension of the system - -#include -#include -#include -#include -#include "MyMiscellany.h" -#include "CmdLineParser.h" -#include "PPolynomial.h" -#include "FEMTree.h" -#include "Ply.h" -#include "VertexFactory.h" -#include "Image.h" -#include "RegularGrid.h" - -MessageWriter messageWriter; - -const float DefaultPointWeightMultiplier = 2.f; - -enum NormalType -{ - NORMALS_NONE , - NORMALS_SAMPLES , - NORMALS_GRADIENTS , - NORMALS_COUNT -}; -const char* NormalsNames[] = { "none" , "samples" , "gradients" }; - -cmdLineParameter< char* > - In( "in" ) , - Out( "out" ) , - TempDir( "tempDir" ) , - Grid( "grid" ) , - Tree( "tree" ) , - Envelope( "envelope" ) , - EnvelopeGrid( "envelopeGrid" ), - Transform( "xForm" ); - -cmdLineReadable - Performance( "performance" ) , - ShowResidual( "showResidual" ) , - NoComments( "noComments" ) , - PolygonMesh( "polygonMesh" ) , - NonManifold( "nonManifold" ) , - ASCII( "ascii" ) , - Density( "density" ) , - LinearFit( "linearFit" ) , - PrimalGrid( "primalGrid" ) , - ExactInterpolation( "exact" ) , - Colors( "colors" ) , - InCore( "inCore" ) , - NoDirichletErode( "noErode" ) , - Verbose( "verbose" ); - -cmdLineParameter< int > -#ifndef FAST_COMPILE - Degree( "degree" , DEFAULT_FEM_DEGREE ) , -#endif // !FAST_COMPILE - Depth( "depth" , 8 ) , - KernelDepth( "kernelDepth" ) , - SolveDepth( "solveDepth" ) , - EnvelopeDepth( "envelopeDepth" ) , - Iters( "iters" , 8 ) , - FullDepth( "fullDepth" , 5 ) , - BaseDepth( "baseDepth" ) , - BaseVCycles( "baseVCycles" , 1 ) , -#ifndef FAST_COMPILE - BType( "bType" , DEFAULT_FEM_BOUNDARY+1 ) , -#endif // !FAST_COMPILE - Normals( "normals" , NORMALS_NONE ) , - MaxMemoryGB( "maxMemory" , 0 ) , -#ifdef _OPENMP - ParallelType( "parallel" , (int)ThreadPool::OPEN_MP ) , -#else // !_OPENMP - ParallelType( "parallel" , (int)ThreadPool::THREAD_POOL ) , -#endif // _OPENMP - ScheduleType( "schedule" , (int)ThreadPool::DefaultSchedule ) , - ThreadChunkSize( "chunkSize" , (int)ThreadPool::DefaultChunkSize ) , - Threads( "threads" , (int)std::thread::hardware_concurrency() ); - -cmdLineParameter< float > - DataX( "data" , 32.f ) , - SamplesPerNode( "samplesPerNode" , 1.5f ) , - Scale( "scale" , 1.1f ) , - Width( "width" , 0.f ) , - Confidence( "confidence" , 0.f ) , - ConfidenceBias( "confidenceBias" , 0.f ) , - CGSolverAccuracy( "cgAccuracy" , 1e-3f ) , - LowDepthCutOff( "lowDepthCutOff" , 0.f ) , - PointWeight( "pointWeight" ); - -cmdLineReadable* params[] = -{ -#ifndef FAST_COMPILE - &Degree , &BType , -#endif // !FAST_COMPILE - &In , &Depth , &Out , &Transform , - &SolveDepth , - &Envelope , - &EnvelopeGrid , - &Width , - &Scale , &Verbose , &CGSolverAccuracy , &NoComments , - &KernelDepth , &SamplesPerNode , &Confidence , &NonManifold , &PolygonMesh , &ASCII , &ShowResidual , - &EnvelopeDepth , - &NoDirichletErode , - &ConfidenceBias , - &BaseDepth , &BaseVCycles , - &PointWeight , - &Grid , &Threads , - &Tree , - &Density , - &FullDepth , - &Iters , - &DataX , - &Colors , - &Normals , - &LinearFit , - &PrimalGrid , - &TempDir , - &ExactInterpolation , - &Performance , - &MaxMemoryGB , - &InCore , - &ParallelType , - &ScheduleType , - &ThreadChunkSize , - &LowDepthCutOff , - NULL -}; - -void ShowUsage(char* ex) -{ - printf( "Usage: %s\n" , ex ); - printf( "\t --%s \n" , In.name ); - printf( "\t[--%s \n" , Envelope.name ); - printf( "\t[--%s ]\n" , Out.name ); - printf( "\t[--%s ]\n" , Grid.name ); - printf( "\t[--%s \n" , EnvelopeGrid.name ); - printf( "\t[--%s ]\n" , Tree.name ); -#ifndef FAST_COMPILE - printf( "\t[--%s =%d]\n" , Degree.name , Degree.value ); - printf( "\t[--%s =%d]\n" , BType.name , BType.value ); - for( int i=0 ; i=%d]\n" , Depth.name , Depth.value ); - printf( "\t[--%s =%d]\n" , SolveDepth.name , SolveDepth.value ); - printf( "\t[--%s ]\n" , Width.name ); - printf( "\t[--%s =%d]\n" , FullDepth.name , FullDepth.value ); - printf( "\t[--%s =%d]\n" , EnvelopeDepth.name , EnvelopeDepth.value ); - printf( "\t[--%s ]\n" , BaseDepth.name ); - printf( "\t[--%s =%d]\n" , BaseVCycles.name , BaseVCycles.value ); - printf( "\t[--%s =%f]\n" , Scale.name , Scale.value ); - printf( "\t[--%s =%f]\n" , SamplesPerNode.name, SamplesPerNode.value ); - printf( "\t[--%s =%.3e * ]\n" , PointWeight.name , DefaultPointWeightMultiplier * DEFAULT_FEM_DEGREE ); - printf( "\t[--%s =%d]\n" , Iters.name , Iters.value ); - printf( "\t[--%s]\n" , ExactInterpolation.name ); - printf( "\t[--%s =%f]\n" , DataX.name , DataX.value ); - printf( "\t[--%s]\n" , Colors.name ); - printf( "\t[--%s =%d]\n" , Normals.name , Normals.value ); - for( int i=0 ; i=%d]\n" , Threads.name , Threads.value ); - printf( "\t[--%s =%d]\n" , ParallelType.name , ParallelType.value ); - for( size_t i=0 ; i=%d]\n" , ScheduleType.name , ScheduleType.value ); - for( size_t i=0 ; i=%d]\n" , ThreadChunkSize.name , ThreadChunkSize.value ); - printf( "\t[--%s =%f]\n" , LowDepthCutOff.name , LowDepthCutOff.value ); - printf( "\t[--%s =%f]\n" , Confidence.name , Confidence.value ); - printf( "\t[--%s =%f]\n" , ConfidenceBias.name , ConfidenceBias.value ); - printf( "\t[--%s]\n" , NonManifold.name ); - printf( "\t[--%s]\n" , PolygonMesh.name ); - printf( "\t[--%s =%g]\n" , CGSolverAccuracy.name , CGSolverAccuracy.value ); - printf( "\t[--%s =%d]\n" , MaxMemoryGB.name , MaxMemoryGB.value ); - printf( "\t[--%s]\n" , NoDirichletErode.name ); - printf( "\t[--%s]\n" , Performance.name ); - printf( "\t[--%s]\n" , Density.name ); - printf( "\t[--%s]\n" , LinearFit.name ); - printf( "\t[--%s]\n" , PrimalGrid.name ); - printf( "\t[--%s]\n" , ASCII.name ); - printf( "\t[--%s]\n" , NoComments.name ); - printf( "\t[--%s]\n" , TempDir.name ); - printf( "\t[--%s]\n" , InCore.name ); - printf( "\t[--%s]\n" , Verbose.name ); -} - -double Weight( double v , double start , double end ) -{ - v = ( v - start ) / ( end - start ); - if ( v<0 ) return 1.; - else if( v>1 ) return 0.; - else - { - // P(x) = a x^3 + b x^2 + c x + d - // P (0) = 1 , P (1) = 0 , P'(0) = 0 , P'(1) = 0 - // => d = 1 , a + b + c + d = 0 , c = 0 , 3a + 2b + c = 0 - // => c = 0 , d = 1 , a + b = -1 , 3a + 2b = 0 - // => a = 2 , b = -3 , c = 0 , d = 1 - // => P(x) = 2 x^3 - 3 x^2 + 1 - return 2. * v * v * v - 3. * v * v + 1.; - } -} - -template< unsigned int Dim , class Real > -struct FEMTreeProfiler -{ - double t; - - void start( void ){ t = Time() , FEMTree< Dim , Real >::ResetLocalMemoryUsage(); } - void print( const char* header ) const - { - FEMTree< Dim , Real >::MemoryUsage(); - if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); - else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); - } - void dumpOutput( const char* header ) const - { - FEMTree< Dim , Real >::MemoryUsage(); - if( header ) messageWriter( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); - else messageWriter( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); - } - void dumpOutput2( std::vector< std::string >& comments , const char* header ) const - { - FEMTree< Dim , Real >::MemoryUsage(); - if( header ) messageWriter( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); - else messageWriter( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); - } -}; - -template< class Real , unsigned int Dim > -XForm< Real , Dim+1 > GetBoundingBoxXForm( Point< Real , Dim > min , Point< Real , Dim > max , Real scaleFactor ) -{ - Point< Real , Dim > center = ( max + min ) / 2; - Real scale = max[0] - min[0]; - for( int d=1 ; d( scale , max[d]-min[d] ); - scale *= scaleFactor; - for( int i=0 ; i tXForm = XForm< Real , Dim+1 >::Identity() , sXForm = XForm< Real , Dim+1 >::Identity(); - for( int i=0 ; i -XForm< Real , Dim+1 > GetBoundingBoxXForm( Point< Real , Dim > min , Point< Real , Dim > max , Real width , Real scaleFactor , int& depth ) -{ - // Get the target resolution (along the largest dimension) - Real resolution = ( max[0]-min[0] ) / width; - for( int d=1 ; d( resolution , ( max[d]-min[d] ) / width ); - resolution *= scaleFactor; - depth = 0; - while( (1< center = ( max + min ) / 2; - Real scale = (1< tXForm = XForm< Real , Dim+1 >::Identity() , sXForm = XForm< Real , Dim+1 >::Identity(); - for( int i=0 ; i -using InputOrientedPointStreamInfo = typename FEMTreeInitializer< Dim , Real >::template InputPointStream< VectorTypeUnion< Real , typename VertexFactory::NormalFactory< Real , Dim >::VertexType , AuxData > >; - -template< typename Real , unsigned int Dim , typename AuxData > -using InputOrientedPointStream = typename InputOrientedPointStreamInfo< Real , Dim , AuxData >::StreamType; - -template< class Real , unsigned int Dim , typename AuxData > -XForm< Real , Dim+1 > GetPointXForm( InputOrientedPointStream< Real , Dim , AuxData > &stream , Real width , Real scaleFactor , int& depth ) -{ - Point< Real , Dim > min , max; - InputOrientedPointStreamInfo< Real , Dim , AuxData >::BoundingBox( stream , min , max ); - return GetBoundingBoxXForm( min , max , width , scaleFactor , depth ); -} -template< class Real , unsigned int Dim , typename AuxData > -XForm< Real , Dim+1 > GetPointXForm( InputOrientedPointStream< Real , Dim , AuxData > &stream , Real scaleFactor ) -{ - Point< Real , Dim > min , max; - InputOrientedPointStreamInfo< Real , Dim , AuxData >::BoundingBox( stream , min , max ); - return GetBoundingBoxXForm( min , max , scaleFactor ); -} - -template< unsigned int Dim , typename Real > -struct ConstraintDual -{ - Real target , weight; - ConstraintDual( Real t , Real w ) : target(t) , weight(w){ } - CumulativeDerivativeValues< Real , Dim , 0 > operator()( const Point< Real , Dim >& p ) const { return CumulativeDerivativeValues< Real , Dim , 0 >( target*weight ); }; -}; -template< unsigned int Dim , typename Real > -struct SystemDual -{ - Real weight; - SystemDual( Real w ) : weight(w){ } - CumulativeDerivativeValues< Real , Dim , 0 > operator()( const Point< Real , Dim >& p , const CumulativeDerivativeValues< Real , Dim , 0 >& dValues ) const { return dValues * weight; }; - CumulativeDerivativeValues< double , Dim , 0 > operator()( const Point< Real , Dim >& p , const CumulativeDerivativeValues< double , Dim , 0 >& dValues ) const { return dValues * weight; }; -}; -template< unsigned int Dim > -struct SystemDual< Dim , double > -{ - typedef double Real; - Real weight; - SystemDual( Real w ) : weight(w){ } - CumulativeDerivativeValues< Real , Dim , 0 > operator()( const Point< Real , Dim >& p , const CumulativeDerivativeValues< Real , Dim , 0 >& dValues ) const { return dValues * weight; }; -}; - -template< typename Real , typename SetVertexFunction , typename InputSampleDataType , typename VertexFactory , unsigned int ... FEMSigs > -void ExtractMesh -( - UIntPack< FEMSigs ... > , - FEMTree< sizeof ... ( FEMSigs ) , Real >& tree , - const DenseNodeData< Real , UIntPack< FEMSigs ... > >& solution , - Real isoValue , - const std::vector< typename FEMTree< sizeof ... ( FEMSigs ) , Real >::PointSample > *samples , - std::vector< InputSampleDataType > *sampleData , - const typename FEMTree< sizeof ... ( FEMSigs ) , Real >::template DensityEstimator< WEIGHT_DEGREE > *density , - const VertexFactory &vertexFactory , - const InputSampleDataType &zeroInputSampleDataType , - SetVertexFunction SetVertex , - std::vector< std::string > &comments , - XForm< Real , sizeof...(FEMSigs)+1 > unitCubeToModel -) -{ - static const int Dim = sizeof ... ( FEMSigs ); - typedef UIntPack< FEMSigs ... > Sigs; - typedef typename VertexFactory::VertexType Vertex; - - static const unsigned int DataSig = FEMDegreeAndBType< DATA_DEGREE , BOUNDARY_FREE >::Signature; - typedef typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE > DensityEstimator; - - FEMTreeProfiler< Dim , Real > profiler; - - char tempHeader[1024]; - { - char tempPath[1024]; - tempPath[0] = 0; - if( TempDir.set ) strcpy( tempPath , TempDir.value ); - else SetTempDirectory( tempPath , sizeof(tempPath) ); - if( strlen(tempPath)==0 ) sprintf( tempPath , ".%c" , FileSeparator ); - if( tempPath[ strlen( tempPath )-1 ]==FileSeparator ) sprintf( tempHeader , "%sPR_" , tempPath ); - else sprintf( tempHeader , "%s%cPR_" , tempPath , FileSeparator ); - } - CoredMeshData< Vertex , node_index_type > *mesh; - if( InCore.set ) mesh = new CoredVectorMeshData< Vertex , node_index_type >(); - else mesh = new CoredFileMeshData< node_index_type , VertexFactory >( vertexFactory , tempHeader ); - - profiler.start(); - typename IsoSurfaceExtractor< Dim , Real , Vertex >::IsoStats isoStats; - if( sampleData ) - { - SparseNodeData< ProjectiveData< InputSampleDataType , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setExtrapolatedDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL ); - for( const RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) ) - { - ProjectiveData< InputSampleDataType , Real >* clr = _sampleData( n ); - if( clr ) (*clr) *= (Real)pow( DataX.value , tree.depth( n ) ); - } - isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< InputSampleDataType >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , density , &_sampleData , solution , isoValue , *mesh , zeroInputSampleDataType , SetVertex , !LinearFit.set , Normals.value==NORMALS_GRADIENTS , !NonManifold.set , PolygonMesh.set , false ); - } -#if defined( __GNUC__ ) && __GNUC__ < 5 -#ifdef SHOW_WARNINGS -#warning "you've got me gcc version<5" -#endif // SHOW_WARNINGS - else isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< InputSampleDataType >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , density , (SparseNodeData< ProjectiveData< InputSampleDataType , Real > , IsotropicUIntPack< Dim , DataSig > > *)NULL , solution , isoValue , *mesh , zeroInputSampleDataType , SetVertex , !LinearFit.set , Normals.value==NORMALS_GRADIENTS , !NonManifold.set , PolygonMesh.set , false ); -#else // !__GNUC__ || __GNUC__ >=5 - else isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< InputSampleDataType >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , density , NULL , solution , isoValue , *mesh , zeroInputSampleDataType , SetVertex , !LinearFit.set , Normals.value==NORMALS_GRADIENTS , !NonManifold.set , PolygonMesh.set , false ); -#endif // __GNUC__ || __GNUC__ < 4 - messageWriter( "Vertices / Polygons: %llu / %llu\n" , (unsigned long long)( mesh->outOfCoreVertexNum()+mesh->inCoreVertices.size() ) , (unsigned long long)mesh->polygonNum() ); - std::string isoStatsString = isoStats.toString() + std::string( "\n" ); - messageWriter( isoStatsString.c_str() ); - if( PolygonMesh.set ) profiler.dumpOutput2( comments , "# Got polygons:" ); - else profiler.dumpOutput2( comments , "# Got triangles:" ); - - std::vector< std::string > noComments; - typename VertexFactory::Transform unitCubeToModelTransform( unitCubeToModel ); - PLY::WritePolygons< VertexFactory , node_index_type , Real , Dim >( Out.value , vertexFactory , mesh , ASCII.set ? PLY_ASCII : PLY_BINARY_NATIVE , NoComments.set ? noComments : comments , unitCubeToModelTransform ); - delete mesh; -} - -template< typename Real , unsigned int Dim > -void WriteGrid( const char *fileName , ConstPointer( Real ) values , unsigned int res , XForm< Real , Dim+1 > voxelToModel , bool verbose ) -{ - char *ext = GetFileExtension( fileName ); - - if( Dim==2 && ImageWriter::ValidExtension( ext ) ) - { - unsigned int totalResolution = 1; - for( int d=0 ; d avgs( ThreadPool::NumThreads() , 0 ); - ThreadPool::Parallel_for( 0 , totalResolution , [&]( unsigned int thread , size_t i ){ avgs[thread] += values[i]; } ); - for( unsigned int t=0 ; t stds( ThreadPool::NumThreads() , 0 ); - ThreadPool::Parallel_for( 0 , totalResolution , [&]( unsigned int thread , size_t i ){ stds[thread] += ( values[i] - avg ) * ( values[i] - avg ); } ); - for( unsigned int t=0 ; t [0,255]\n" , avg - 2*std , avg + 2*std ); - printf( "Transform:\n" ); - for( int i=0 ; i( (Real)1. , std::max< Real >( (Real)-1. , ( values[i] - avg ) / (2*std ) ) ); - v = (Real)( ( v + 1. ) / 2. * 256. ); - unsigned char color = (unsigned char )std::min< Real >( (Real)255. , std::max< Real >( (Real)0. , v ) ); - for( int c=0 ; c<3 ; c++ ) pixels[i*3+c ] = color; - } - ); - ImageWriter::Write( fileName , pixels , res , res , 3 ); - delete[] pixels; - } - else if( !strcasecmp( ext , "iso" ) ) - { - FILE *fp = fopen( fileName , "wb" ); - if( !fp ) ERROR_OUT( "Failed to open file for writing: " , fileName ); - int r = (int)res; - fwrite( &r , sizeof(int) , 1 , fp ); - size_t count = 1; - for( unsigned int d=0 ; d::Write( fileName , _res , values , voxelToModel ); - } - delete[] ext; -} - -template< class Real , typename AuxDataFactory , unsigned int ... FEMSigs > -void Execute( UIntPack< FEMSigs ... > , const AuxDataFactory &auxDataFactory ) -{ - static const int Dim = sizeof ... ( FEMSigs ); - typedef UIntPack< FEMSigs ... > Sigs; - typedef UIntPack< FEMSignature< FEMSigs >::Degree ... > Degrees; - typedef UIntPack< FEMDegreeAndBType< NORMAL_DEGREE , DerivativeBoundary< FEMSignature< FEMSigs >::BType , 1 >::BType >::Signature ... > NormalSigs; - static const unsigned int DataSig = FEMDegreeAndBType< DATA_DEGREE , BOUNDARY_FREE >::Signature; - typedef typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE > DensityEstimator; - typedef typename FEMTree< Dim , Real >::template InterpolationInfo< Real , 0 > InterpolationInfo; - using namespace VertexFactory; - - // The factory for constructing an input sample - typedef Factory< Real , PositionFactory< Real , Dim > , Factory< Real , NormalFactory< Real , Dim > , AuxDataFactory > > InputSampleFactory; - - // The factory for constructing an input sample's data - typedef Factory< Real , NormalFactory< Real , Dim > , AuxDataFactory > InputSampleDataFactory; - - // The input point stream information: First piece of data is the normal; the remainder is the auxiliary data - typedef InputOrientedPointStreamInfo< Real , Dim , typename AuxDataFactory::VertexType > InputPointStreamInfo; - - // The type of the input sample - typedef typename InputPointStreamInfo::PointAndDataType InputSampleType; - - // The type of the input sample's data - typedef typename InputPointStreamInfo::DataType InputSampleDataType; - - typedef InputDataStream< InputSampleType > InputPointStream; - typedef TransformedInputDataStream< InputSampleType > XInputPointStream; - - InputSampleFactory inputSampleFactory( PositionFactory< Real , Dim >() , InputSampleDataFactory( NormalFactory< Real , Dim >() , auxDataFactory ) ); - InputSampleDataFactory inputSampleDataFactory( NormalFactory< Real , Dim >() , auxDataFactory ); - - typedef RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type > FEMTreeNode; - typedef typename FEMTreeInitializer< Dim , Real >::GeometryNodeType GeometryNodeType; - std::vector< std::string > comments; - messageWriter( comments , "*************************************************************\n" ); - messageWriter( comments , "*************************************************************\n" ); - messageWriter( comments , "** Running Screened Poisson Reconstruction (Version %s) **\n" , VERSION ); - messageWriter( comments , "*************************************************************\n" ); - messageWriter( comments , "*************************************************************\n" ); - if( !Threads.set ) messageWriter( comments , "Running with %d threads\n" , Threads.value ); - - bool needNormalData = DataX.value>0 && Normals.value; - bool needAuxData = DataX.value>0 && auxDataFactory.bufferSize(); - - XForm< Real , Dim+1 > modelToUnitCube , unitCubeToModel; - if( Transform.set ) - { - FILE* fp = fopen( Transform.value , "r" ); - if( !fp ) - { - WARN( "Could not read x-form from: " , Transform.value ); - modelToUnitCube = XForm< Real , Dim+1 >::Identity(); - } - else - { - for( int i=0 ; i::Identity(); - - char str[1024]; - for( int i=0 ; params[i] ; i++ ) - if( params[i]->set ) - { - params[i]->writeValue( str ); - if( strlen( str ) ) messageWriter( comments , "\t--%s %s\n" , params[i]->name , str ); - else messageWriter( comments , "\t--%s\n" , params[i]->name ); - } - - double startTime = Time(); - Real isoValue = 0; - - FEMTree< Dim , Real > tree( MEMORY_ALLOCATOR_BLOCK_SIZE ); - FEMTreeProfiler< Dim , Real > profiler; - - if( Depth.set && Width.value>0 ) - { - WARN( "Both --" , Depth.name , " and --" , Width.name , " set, ignoring --" , Width.name ); - Width.value = 0; - } - - size_t pointCount; - - Real pointWeightSum; - std::vector< typename FEMTree< Dim , Real >::PointSample >* samples = new std::vector< typename FEMTree< Dim , Real >::PointSample >(); - DenseNodeData< GeometryNodeType , IsotropicUIntPack< Dim , FEMTrivialSignature > > geometryNodeDesignators; - std::vector< InputSampleDataType >* sampleData = NULL; - DensityEstimator* density = NULL; - SparseNodeData< Point< Real , Dim > , NormalSigs >* normalInfo = NULL; - Real targetValue = (Real)0.5; - - // Read in the samples (and color data) - { - profiler.start(); - InputPointStream* pointStream; - char* ext = GetFileExtension( In.value ); - sampleData = new std::vector< InputSampleDataType >(); - std::vector< InputSampleType > inCorePoints; - if( InCore.set ) - { - InputPointStream *_pointStream; - if ( !strcasecmp( ext , "bnpts" ) ) _pointStream = new BinaryInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); - else if( !strcasecmp( ext , "ply" ) ) _pointStream = new PLYInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); - else _pointStream = new ASCIIInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); - InputSampleType p; - while( _pointStream->next( p ) ) inCorePoints.push_back( p ); - delete _pointStream; - - pointStream = new MemoryInputDataStream< InputSampleType >( inCorePoints.size() , &inCorePoints[0] ); - } - else - { - if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); - else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); - else pointStream = new ASCIIInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); - } - delete[] ext; - - typename InputSampleFactory::Transform _modelToUnitCube( modelToUnitCube ); - auto XFormFunctor = [&]( InputSampleType &p ){ p = _modelToUnitCube( p ); }; - XInputPointStream _pointStream( XFormFunctor , *pointStream ); - if( Width.value>0 ) modelToUnitCube = GetPointXForm< Real , Dim , typename AuxDataFactory::VertexType >( _pointStream , Width.value , (Real)( Scale.value>0 ? Scale.value : 1. ) , Depth.value ) * modelToUnitCube; - else modelToUnitCube = Scale.value>0 ? GetPointXForm< Real , Dim , typename AuxDataFactory::VertexType >( _pointStream , (Real)Scale.value ) * modelToUnitCube : modelToUnitCube; - - if( !SolveDepth.set ) SolveDepth.value = Depth.value; - if( SolveDepth.value>Depth.value ) - { - WARN( "Solution depth cannot exceed system depth: " , SolveDepth.value , " <= " , Depth.value ); - SolveDepth.value = Depth.value; - } - - { - typename InputSampleFactory::Transform _modelToUnitCube( modelToUnitCube ); - auto XFormFunctor = [&]( InputSampleType &p ){ p = _modelToUnitCube( p ); }; - XInputPointStream _pointStream( XFormFunctor , *pointStream ); - auto ProcessDataWithConfidence = [&]( const Point< Real , Dim > &p , typename InputPointStreamInfo::DataType &d ) - { - Real l = (Real)Length( d.template get<0>() ); - if( !l || !std::isfinite( l ) ) return (Real)-1.; - return (Real)pow( l , Confidence.value ); - }; - auto ProcessData = []( const Point< Real , Dim > &p , typename InputPointStreamInfo::DataType &d ) - { - Real l = (Real)Length( d.template get<0>() ); - if( !l || !std::isfinite( l ) ) return (Real)-1.; - d.template get<0>() /= l; - return (Real)1.; - }; - if( Confidence.value>0 ) pointCount = FEMTreeInitializer< Dim , Real >::template Initialize< InputSampleDataType >( tree.spaceRoot() , _pointStream , Depth.value , *samples , *sampleData , true , tree.nodeAllocators[0] , tree.initializer() , ProcessDataWithConfidence ); - else pointCount = FEMTreeInitializer< Dim , Real >::template Initialize< InputSampleDataType >( tree.spaceRoot() , _pointStream , Depth.value , *samples , *sampleData , true , tree.nodeAllocators[0] , tree.initializer() , ProcessData ); - } - - unitCubeToModel = modelToUnitCube.inverse(); - delete pointStream; - - messageWriter( "Input Points / Samples: %llu / %llu\n" , (unsigned long long)pointCount , (unsigned long long)samples->size() ); - profiler.dumpOutput2( comments , "# Read input into tree:" ); - } - - DenseNodeData< Real , Sigs > solution; - { - DenseNodeData< Real , Sigs > constraints; - InterpolationInfo* iInfo = NULL; - int solveDepth = Depth.value; - - tree.resetNodeIndices( 0 , std::make_tuple() ); - - // Get the kernel density estimator - { - profiler.start(); - density = tree.template setDensityEstimator< 1 , WEIGHT_DEGREE >( *samples , KernelDepth.value , SamplesPerNode.value ); - profiler.dumpOutput2( comments , "# Got kernel density:" ); - } - - // Transform the Hermite samples into a vector field - { - profiler.start(); - normalInfo = new SparseNodeData< Point< Real , Dim > , NormalSigs >(); - std::function< bool ( InputSampleDataType , Point< Real , Dim >& ) > ConversionFunction = []( InputSampleDataType in , Point< Real , Dim > &out ) - { - Point< Real , Dim > n = in.template get<0>(); - Real l = (Real)Length( n ); - // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals... - if( !l ) return false; - out = n / l; - return true; - }; - std::function< bool ( InputSampleDataType , Point< Real , Dim >& , Real & ) > ConversionAndBiasFunction = []( InputSampleDataType in , Point< Real , Dim > &out , Real &bias ) - { - Point< Real , Dim > n = in.template get<0>(); - Real l = (Real)Length( n ); - // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals... - if( !l ) return false; - out = n / l; - bias = (Real)( log( l ) * ConfidenceBias.value / log( 1<<(Dim-1) ) ); - return true; - }; -#if 1 - if( ConfidenceBias.value>0 ) *normalInfo = tree.setInterpolatedDataField( NormalSigs() , *samples , *sampleData , density , BaseDepth.value , Depth.value , (Real)LowDepthCutOff.value , pointWeightSum , ConversionAndBiasFunction ); - else *normalInfo = tree.setInterpolatedDataField( NormalSigs() , *samples , *sampleData , density , BaseDepth.value , Depth.value , (Real)LowDepthCutOff.value , pointWeightSum , ConversionFunction ); -#else - if( ConfidenceBias.value>0 ) *normalInfo = tree.setInterpolatedDataField( NormalSigs() , *samples , *sampleData , density , 0 , Depth.value , (Real)LowDepthCutOff.value , pointWeightSum , ConversionAndBiasFunction ); - else *normalInfo = tree.setInterpolatedDataField( NormalSigs() , *samples , *sampleData , density , 0 , Depth.value , (Real)LowDepthCutOff.value , pointWeightSum , ConversionFunction ); -#endif - ThreadPool::Parallel_for( 0 , normalInfo->size() , [&]( unsigned int , size_t i ){ (*normalInfo)[i] *= (Real)-1.; } ); - profiler.dumpOutput2( comments , "# Got normal field:" ); - messageWriter( "Point weight / Estimated Measure: %g / %g\n" , pointWeightSum , pointCount*pointWeightSum ); - } - - // Get the geometry designators indicating if the space node are interior to, exterior to, or contain the envelope boundary - if( Envelope.set ) - { - profiler.start(); - { - // Make the octree complete up to the base depth - std::function< void ( FEMTreeNode * , unsigned int ) > MakeComplete = [&]( FEMTreeNode *node , unsigned int depth ) - { - if( node->depth()<(int)depth ) - { - if( !node->children ) node->template initChildren< false >( tree.nodeAllocators.size() ? tree.nodeAllocators[0] : NULL , tree.initializer() ); - for( int c=0 ; c<(1<children+c , depth ); - } - }; - MakeComplete( &tree.spaceRoot() , BaseDepth.value ); - - // Read in the envelope geometry - std::vector< Point< Real , Dim > > vertices; - std::vector< SimplexIndex< Dim-1 , node_index_type > > simplices; - { - std::vector< typename PositionFactory< Real , Dim >::VertexType > _vertices; - std::vector< std::vector< int > > polygons; - std::vector< std::string > comments; - int file_type; - PLY::ReadPolygons( Envelope.value , PositionFactory< Real , Dim >() , _vertices , polygons , file_type , comments ); - vertices.resize( _vertices.size() ); - for( int i=0 ; i::template GetGeometryNodeDesignators( &tree.spaceRoot() , vertices , simplices , BaseDepth.value , EnvelopeDepth.value , tree.nodeAllocators , tree.initializer() ); - - // Make nodes in the support of the vector field @{ExactDepth} interior - if( !NoDirichletErode.set ) - { - // What to do if we find a node in the support of the vector field - auto SetScratchFlag = [&]( FEMTreeNode *node ) - { - if( node ) - { - while( node->depth()>BaseDepth.value ) node = node->parent; - node->nodeData.setScratchFlag( true ); - } - }; - - std::function< void ( FEMTreeNode * ) > PropagateToLeaves = [&]( const FEMTreeNode *node ) - { - geometryNodeDesignators[ node ] = GeometryNodeType::INTERIOR; - if( node->children ) for( int c=0 ; c<(1<children+c ); - }; - - // Flags indicating if a node contains a non-zero vector field coefficient - std::vector< bool > isVectorFieldElement( tree.nodeCount() , false ); - - // Get the set of base nodes - std::vector< FEMTreeNode * > baseNodes; - auto TerminationLambda = []( const FEMTreeNode *node ){ return node->depth()==BaseDepth.value; }; - for( FEMTreeNode *node=tree.spaceRoot().nextNode( TerminationLambda , NULL ) ; node ; node=tree.spaceRoot().nextNode( TerminationLambda , node ) ) if( node->depth()==BaseDepth.value ) baseNodes.push_back( node ); - - std::vector< node_index_type > vectorFieldElementCounts( baseNodes.size() ); - for( int i=0 ; inextNode() ; node ; node=baseNodes[i]->nextNode(node) ) - { - Point< Real , Dim > *n = (*normalInfo)( node ); - if( n && Point< Real , Dim >::SquareNorm( *n ) ) isVectorFieldElement[ node->nodeData.nodeIndex ] = true , vectorFieldElementCounts[i]++; - } - } ); - size_t vectorFieldElementCount = 0; - for( int i=0 ; i vectorFieldElements; - vectorFieldElements.reserve( vectorFieldElementCount ); - { - std::vector< std::vector< FEMTreeNode * > > _vectorFieldElements( baseNodes.size() ); - for( int i=0 ; i<_vectorFieldElements.size() ; i++ ) _vectorFieldElements[i].reserve( vectorFieldElementCounts[i] ); - ThreadPool::Parallel_for( 0 , baseNodes.size() , [&]( unsigned int t , size_t i ) - { - for( FEMTreeNode *node=baseNodes[i]->nextNode() ; node ; node=baseNodes[i]->nextNode(node) ) - { - if( isVectorFieldElement[ node->nodeData.nodeIndex ] ) _vectorFieldElements[i].push_back( node ); - node->nodeData.setScratchFlag( false ); - } - } ); - for( int i=0 ; i<_vectorFieldElements.size() ; i++ ) vectorFieldElements.insert( vectorFieldElements.end() , _vectorFieldElements[i].begin() , _vectorFieldElements[i].end() ); - } - - // Set the scratch flag for the base nodes on which the vector field is supported -#ifdef SHOW_WARNINGS -#pragma message( "[WARNING] In principal, we should unlock finite elements whose support overlaps the vector field" ) -#endif // SHOW_WARNINGS - tree.template processNeighboringLeaves< -BSplineSupportSizes< NORMAL_DEGREE >::SupportStart , BSplineSupportSizes< NORMAL_DEGREE >::SupportEnd >( &vectorFieldElements[0] , vectorFieldElements.size() , SetScratchFlag , false ); - - // Set sub-trees rooted at interior nodes @ ExactDepth to interior - ThreadPool::Parallel_for( 0 , baseNodes.size() , [&]( unsigned int , size_t i ){ if( baseNodes[i]->nodeData.getScratchFlag() ) PropagateToLeaves( baseNodes[i] ); } ); - - // Adjust the coarser node designators in case exterior nodes have become boundary. - ThreadPool::Parallel_for( 0 , baseNodes.size() , [&]( unsigned int , size_t i ){ FEMTreeInitializer< Dim , Real >::PullGeometryNodeDesignatorsFromFiner( baseNodes[i] , geometryNodeDesignators ); } ); - FEMTreeInitializer< Dim , Real >::PullGeometryNodeDesignatorsFromFiner( &tree.spaceRoot() , geometryNodeDesignators , BaseDepth.value ); - } - } - profiler.dumpOutput2( comments , "# Initialized envelope constraints:" ); - } - - if( !Density.set ) delete density , density = NULL; - if( !needNormalData && !needAuxData ) delete sampleData , sampleData = NULL; - - // Add the interpolation constraints - if( PointWeight.value>0 ) - { - profiler.start(); - if( ExactInterpolation.set ) iInfo = FEMTree< Dim , Real >::template InitializeExactPointInterpolationInfo< Real , 0 > ( tree , *samples , ConstraintDual< Dim , Real >( targetValue , (Real)PointWeight.value * pointWeightSum ) , SystemDual< Dim , Real >( (Real)PointWeight.value * pointWeightSum ) , true , false ); - else iInfo = FEMTree< Dim , Real >::template InitializeApproximatePointInterpolationInfo< Real , 0 > ( tree , *samples , ConstraintDual< Dim , Real >( targetValue , (Real)PointWeight.value * pointWeightSum ) , SystemDual< Dim , Real >( (Real)PointWeight.value * pointWeightSum ) , true , 1 ); - profiler.dumpOutput2( comments , "#Initialized point interpolation constraints:" ); - } - - // Trim the tree and prepare for multigrid - { - profiler.start(); - constexpr int MAX_DEGREE = NORMAL_DEGREE > Degrees::Max() ? NORMAL_DEGREE : Degrees::Max(); - typename FEMTree< Dim , Real >::template HasNormalDataFunctor< NormalSigs > hasNormalDataFunctor( *normalInfo ); - auto hasDataFunctor = [&]( const FEMTreeNode *node ){ return hasNormalDataFunctor( node ); }; - if( geometryNodeDesignators.size() ) tree.template finalizeForMultigrid< MAX_DEGREE , Degrees::Max() >( BaseDepth.value , FullDepth.value , hasDataFunctor , [&]( const FEMTreeNode *node ){ return node->nodeData.nodeIndex<(node_index_type)geometryNodeDesignators.size() && geometryNodeDesignators[node]==GeometryNodeType::EXTERIOR; } , std::make_tuple( iInfo ) , std::make_tuple( normalInfo , density , &geometryNodeDesignators ) ); - else tree.template finalizeForMultigrid< MAX_DEGREE , Degrees::Max() >( BaseDepth.value , FullDepth.value , hasDataFunctor , []( const FEMTreeNode * ){ return false; } , std::make_tuple( iInfo ) , std::make_tuple( normalInfo , density ) ); - - if( geometryNodeDesignators.size() && EnvelopeGrid.set ) - { - FEMTreeInitializer< Dim , Real >::PushGeometryNodeDesignatorsToFiner( &tree.spaceRoot() , geometryNodeDesignators ); - FEMTreeInitializer< Dim , Real >::PullGeometryNodeDesignatorsFromFiner( &tree.spaceRoot() , geometryNodeDesignators ); - - auto WriteEnvelopeGrid = [&]( bool showFinest ) - { - int res = 0; - DenseNodeData< Real , IsotropicUIntPack< Dim , FEMTrivialSignature > > coefficients = tree.initDenseNodeData( IsotropicUIntPack< Dim , FEMTrivialSignature >() ); - for( const FEMTreeNode* n = tree.spaceRoot().nextNode() ; n ; n=tree.spaceRoot().nextNode( n ) ) - if( n->nodeData.nodeIndex!=-1 && ( ( showFinest && !n->children ) || ( !showFinest && geometryNodeDesignators[n->parent]==GeometryNodeType::BOUNDARY ) ) ) - { -#if 0 // Randomize the colors - auto Value = []( double v , double eps ){ return (Real)( v + Random< double >() * 2. * eps - eps ); }; - // Show the octree structure - if ( geometryNodeDesignators[n]==GeometryNodeType::INTERIOR ) coefficients[n] = Value( 0.75 , 0.25 ); - else if( geometryNodeDesignators[n]==GeometryNodeType::EXTERIOR ) coefficients[n] = Value( -0.75 , 0.25 ); - -#else - if ( geometryNodeDesignators[n]==GeometryNodeType::INTERIOR ) coefficients[n] = (Real) 1.; - else if( geometryNodeDesignators[n]==GeometryNodeType::EXTERIOR ) coefficients[n] = (Real)-1.; -#endif - } - Pointer( Real ) values = tree.template regularGridEvaluate< true >( coefficients , res , -1 , false ); - XForm< Real , Dim+1 > voxelToUnitCube = XForm< Real , Dim+1 >::Identity(); - for( int d=0 ; d( EnvelopeGrid.value , values , res , unitCubeToModel * voxelToUnitCube , Verbose.set ); - DeletePointer( values ); - }; - - WriteEnvelopeGrid( true ); - } - - profiler.dumpOutput2( comments , "# Finalized tree:" ); - } - // Add the FEM constraints - { - profiler.start(); - constraints = tree.initDenseNodeData( Sigs() ); - typename FEMIntegrator::template Constraint< Sigs , IsotropicUIntPack< Dim , 1 > , NormalSigs , IsotropicUIntPack< Dim , 0 > , Dim > F; - unsigned int derivatives2[Dim]; - for( int d=0 ; d Derivatives1; - typedef IsotropicUIntPack< Dim , 0 > Derivatives2; - for( int d=0 ; d::Index( derivatives1 ) ][ TensorDerivatives< Derivatives2 >::Index( derivatives2 ) ] = 1; - } - tree.addFEMConstraints( F , *normalInfo , constraints , solveDepth ); - profiler.dumpOutput2( comments , "# Set FEM constraints:" ); - } - - // Free up the normal info - delete normalInfo , normalInfo = NULL; - - // Add the interpolation constraints - if( PointWeight.value>0 ) - { - profiler.start(); - tree.addInterpolationConstraints( constraints , solveDepth , std::make_tuple( iInfo ) ); - profiler.dumpOutput2( comments , "#Set point constraints:" ); - } - - messageWriter( "Leaf Nodes / Active Nodes / Ghost Nodes / Dirichlet Supported Nodes: %llu / %llu / %llu / %llu\n" , (unsigned long long)tree.leaves() , (unsigned long long)tree.nodes() , (unsigned long long)tree.ghostNodes() , (unsigned long long)tree.dirichletElements() ); - messageWriter( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) ); - - // Solve the linear system - { - profiler.start(); - typename FEMTree< Dim , Real >::SolverInfo sInfo; - sInfo.cgDepth = 0 , sInfo.cascadic = true , sInfo.vCycles = 1 , sInfo.iters = Iters.value , sInfo.cgAccuracy = CGSolverAccuracy.value , sInfo.verbose = Verbose.set , sInfo.showResidual = ShowResidual.set , sInfo.showGlobalResidual = SHOW_GLOBAL_RESIDUAL_NONE , sInfo.sliceBlockSize = 1; - sInfo.baseVCycles = BaseVCycles.value; - typename FEMIntegrator::template System< Sigs , IsotropicUIntPack< Dim , 1 > > F( { 0. , 1. } ); - solution = tree.solveSystem( Sigs() , F , constraints , SolveDepth.value , sInfo , std::make_tuple( iInfo ) ); - profiler.dumpOutput2( comments , "# Linear system solved:" ); - if( iInfo ) delete iInfo , iInfo = NULL; - } - } - - { - profiler.start(); - double valueSum = 0 , weightSum = 0; - typename FEMTree< Dim , Real >::template MultiThreadedEvaluator< Sigs , 0 > evaluator( &tree , solution ); - std::vector< double > valueSums( ThreadPool::NumThreads() , 0 ) , weightSums( ThreadPool::NumThreads() , 0 ); - ThreadPool::Parallel_for( 0 , samples->size() , [&]( unsigned int thread , size_t j ) - { - ProjectiveData< Point< Real , Dim > , Real >& sample = (*samples)[j].sample; - Real w = sample.weight; - if( w>0 ) weightSums[thread] += w , valueSums[thread] += evaluator.values( sample.data / sample.weight , thread , (*samples)[j].node )[0] * w; - } - ); - for( size_t t=0 ; t::WriteParameter( fp ); - DenseNodeData< Real , Sigs >::WriteSignatures( fp ); - tree.write( fp , modelToUnitCube ); - solution.write( fp ); - if( sampleData ) - { - SparseNodeData< ProjectiveData< InputSampleDataType , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setExtrapolatedDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL ); - for( const RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) ) - { - ProjectiveData< InputSampleDataType , Real >* clr = _sampleData( n ); - if( clr ) (*clr) *= (Real)pow( DataX.value , tree.depth( n ) ); - } - _sampleData.write( fp ); - } - if( density ) density->write( fp ); - fclose( fp ); - } - - if( Grid.set ) - { - int res = 0; - profiler.start(); - Pointer( Real ) values = tree.template regularGridEvaluate< true >( solution , res , -1 , PrimalGrid.set ); - profiler.dumpOutput( "Got grid:" ); - XForm< Real , Dim+1 > voxelToUnitCube = XForm< Real , Dim+1 >::Identity(); - if( PrimalGrid.set ) for( int d=0 ; d( Grid.value , values , res , unitCubeToModel * voxelToUnitCube , Verbose.set ); - DeletePointer( values ); - } - - if( Out.set ) - { - if( Normals.value ) - { - if( Density.set ) - { - typedef Factory< Real , PositionFactory< Real , Dim > , NormalFactory< Real , Dim > , ValueFactory< Real > , AuxDataFactory > VertexFactory; - VertexFactory vertexFactory( PositionFactory< Real , Dim >() , NormalFactory< Real , Dim >() , ValueFactory< Real >() , auxDataFactory ); - if( Normals.value==NORMALS_SAMPLES ) - { - auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = d.template get<0>() , v.template get<2>() = w , v.template get<3>() = d.template get<1>(); }; - ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); - } - else if( Normals.value==NORMALS_GRADIENTS ) - { - auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = -g/(1<() = w , v.template get<3>() = d.template get<1>(); }; - ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); - } - } - else - { - typedef Factory< Real , PositionFactory< Real , Dim > , NormalFactory< Real , Dim > , AuxDataFactory > VertexFactory; - VertexFactory vertexFactory( PositionFactory< Real , Dim >() , NormalFactory< Real , Dim >() , auxDataFactory ); - if( Normals.value==NORMALS_SAMPLES ) - { - auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = d.template get<0>() , v.template get<2>() = d.template get<1>(); }; - ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); - } - else if( Normals.value==NORMALS_GRADIENTS ) - { - auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = -g/(1<() = d.template get<1>(); }; - ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); - } - } - } - else - { - if( Density.set ) - { - typedef Factory< Real , PositionFactory< Real , Dim > , ValueFactory< Real > , AuxDataFactory > VertexFactory; - VertexFactory vertexFactory( PositionFactory< Real , Dim >() , ValueFactory< Real >() , auxDataFactory ); - auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = w , v.template get<2>() = d.template get<1>(); }; - ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); - } - else - { - typedef Factory< Real , PositionFactory< Real , Dim > , AuxDataFactory > VertexFactory; - VertexFactory vertexFactory( PositionFactory< Real , Dim >() , auxDataFactory ); - auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = d.template get<1>(); }; - ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); - } - } - if( sampleData ){ delete sampleData ; sampleData = NULL; } - } - if( density ) delete density , density = NULL; - messageWriter( comments , "# Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-startTime , FEMTree< Dim , Real >::MaxMemoryUsage() ); -} - -#ifndef FAST_COMPILE -template< unsigned int Dim , class Real , BoundaryType BType , typename AuxDataFactory > -void Execute( const AuxDataFactory &auxDataFactory ) -{ - switch( Degree.value ) - { - case 1: return Execute< Real >( IsotropicUIntPack< Dim , FEMDegreeAndBType< 1 , BType >::Signature >() , auxDataFactory ); - case 2: return Execute< Real >( IsotropicUIntPack< Dim , FEMDegreeAndBType< 2 , BType >::Signature >() , auxDataFactory ); -// case 3: return Execute< Real >( IsotropicUIntPack< Dim , FEMDegreeAndBType< 3 , BType >::Signature >() , auxDataFactory ); -// case 4: return Execute< Real >( IsotropicUIntPack< Dim , FEMDegreeAndBType< 4 , BType >::Signature >() , auxDataFactory ); - default: ERROR_OUT( "Only B-Splines of degree 1 - 2 are supported" ); - } -} - -template< unsigned int Dim , class Real , typename AuxDataFactory > -void Execute( const AuxDataFactory &auxDataFactory ) -{ - switch( BType.value ) - { - case BOUNDARY_FREE+1: return Execute< Dim , Real , BOUNDARY_FREE >( auxDataFactory ); - case BOUNDARY_NEUMANN+1: return Execute< Dim , Real , BOUNDARY_NEUMANN >( auxDataFactory ); - case BOUNDARY_DIRICHLET+1: return Execute< Dim , Real , BOUNDARY_DIRICHLET >( auxDataFactory ); - default: ERROR_OUT( "Not a valid boundary type: " , BType.value ); - } -} -#endif // !FAST_COMPILE - -int main( int argc , char* argv[] ) -{ - Timer timer; -#ifdef USE_SEG_FAULT_HANDLER - WARN( "using seg-fault handler" ); - StackTracer::exec = argv[0]; - signal( SIGSEGV , SignalHandler ); -#endif // USE_SEG_FAULT_HANDLER -#ifdef ARRAY_DEBUG - WARN( "Array debugging enabled" ); -#endif // ARRAY_DEBUG - cmdLineParse( argc-1 , &argv[1] , params ); - if( MaxMemoryGB.value>0 ) SetPeakMemoryMB( MaxMemoryGB.value<<10 ); - ThreadPool::DefaultChunkSize = ThreadChunkSize.value; - ThreadPool::DefaultSchedule = (ThreadPool::ScheduleType)ScheduleType.value; - ThreadPool::Init( (ThreadPool::ParallelType)ParallelType.value , Threads.value ); - - messageWriter.echoSTDOUT = Verbose.set; - if( !In.set ) - { - ShowUsage( argv[0] ); - return 0; - } - - if( !BaseDepth.set ) BaseDepth.value = FullDepth.value; - if( BaseDepth.value>FullDepth.value ) - { - if( BaseDepth.set ) WARN( "Base depth must be smaller than full depth: " , BaseDepth.value , " <= " , FullDepth.value ); - BaseDepth.value = FullDepth.value; - } - if( !KernelDepth.set ) KernelDepth.value = Depth.value-2; - if( KernelDepth.value>Depth.value ) - { - WARN( "Kernel depth should not exceed depth: " , KernelDepth.name , " <= " , KernelDepth.value ); - KernelDepth.value = Depth.value; - } - - if( !EnvelopeDepth.set ) EnvelopeDepth.value = BaseDepth.value; - if( EnvelopeDepth.value>Depth.value ) - { - WARN( EnvelopeDepth.name , " can't be greater than " , Depth.name , ": " , EnvelopeDepth.value , " <= " , Depth.value ); - EnvelopeDepth.value = Depth.value; - } - if( EnvelopeDepth.value= " , BaseDepth.value ); - EnvelopeDepth.value = BaseDepth.value; - } - -#ifdef USE_DOUBLE - typedef double Real; -#else // !USE_DOUBLE - typedef float Real; -#endif // USE_DOUBLE - -#ifdef FAST_COMPILE - static const int Degree = DEFAULT_FEM_DEGREE; - static const BoundaryType BType = DEFAULT_FEM_BOUNDARY; - typedef IsotropicUIntPack< DEFAULT_DIMENSION , FEMDegreeAndBType< Degree , BType >::Signature > FEMSigs; - WARN( "Compiled for degree-" , Degree , ", boundary-" , BoundaryNames[ BType ] , ", " , sizeof(Real)==4 ? "single" : "double" , "-precision _only_" ); - if( !PointWeight.set ) PointWeight.value = DefaultPointWeightMultiplier*Degree; - char *ext = GetFileExtension( In.value ); - if( !strcasecmp( ext , "ply" ) ) - { - typedef VertexFactory::Factory< Real , typename VertexFactory::PositionFactory< Real , DEFAULT_DIMENSION > , typename VertexFactory::NormalFactory< Real , DEFAULT_DIMENSION > > Factory; - Factory factory; - bool *readFlags = new bool[ factory.plyReadNum() ]; - std::vector< PlyProperty > unprocessedProperties; - PLY::ReadVertexHeader( In.value , factory , readFlags , unprocessedProperties ); - if( !factory.plyValidReadProperties<0>( readFlags ) ) ERROR_OUT( "Ply file does not contain positions" ); - if( !factory.plyValidReadProperties<1>( readFlags ) ) ERROR_OUT( "Ply file does not contain normals" ); - delete[] readFlags; - - if( unprocessedProperties.size() ) Execute< Real >( FEMSigs() , VertexFactory::DynamicFactory< Real >( unprocessedProperties ) ); - else Execute< Real >( FEMSigs() , VertexFactory::EmptyFactory< Real >() ); - } - else - { - if( Colors.set ) Execute< Real >( FEMSigs() , VertexFactory::RGBColorFactory< Real >() ); - else Execute< Real >( FEMSigs() , VertexFactory::EmptyFactory< Real >() ); - } - delete[] ext; -#else // !FAST_COMPILE - if( !PointWeight.set ) PointWeight.value = DefaultPointWeightMultiplier*Degree.value; - char *ext = GetFileExtension( In.value ); - if( !strcasecmp( ext , "ply" ) ) - { - typedef VertexFactory::Factory< Real , typename VertexFactory::PositionFactory< Real , DEFAULT_DIMENSION > , typename VertexFactory::NormalFactory< Real , DEFAULT_DIMENSION > > Factory; - Factory factory; - bool *readFlags = new bool[ factory.plyReadNum() ]; - std::vector< PlyProperty > unprocessedProperties; - PLY::ReadVertexHeader( In.value , factory , readFlags , unprocessedProperties ); - if( !factory.plyValidReadProperties<0>( readFlags ) ) ERROR_OUT( "Ply file does not contain positions" ); - if( !factory.plyValidReadProperties<1>( readFlags ) ) ERROR_OUT( "Ply file does not contain normals" ); - delete[] readFlags; - - if( unprocessedProperties.size() ) Execute< DEFAULT_DIMENSION , Real >( VertexFactory::DynamicFactory< Real >( unprocessedProperties ) ); - else Execute< DEFAULT_DIMENSION , Real >( VertexFactory::EmptyFactory< Real >() ); - } - else - { - if( Colors.set ) Execute< DEFAULT_DIMENSION , Real >( VertexFactory::RGBColorFactory< Real >() ); - else Execute< DEFAULT_DIMENSION , Real >( VertexFactory::EmptyFactory< Real >() ); - } - delete[] ext; -#endif // FAST_COMPILE - if( Performance.set ) - { - printf( "Time (Wall/CPU): %.2f / %.2f\n" , timer.wallTime() , timer.cpuTime() ); - printf( "Peak Memory (MB): %d\n" , MemoryInfo::PeakMemoryUsageMB() ); - } - - ThreadPool::Terminate(); - return EXIT_SUCCESS; +#include "PoissonReconLib.h" +int main(int argc, char **argv) { + return PoissonReconLib(argc, argv); } diff --git a/Src/PoissonReconLib.cpp b/Src/PoissonReconLib.cpp new file mode 100644 index 00000000..8bacda1e --- /dev/null +++ b/Src/PoissonReconLib.cpp @@ -0,0 +1,1166 @@ +/* +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 "PreProcessor.h" + +#undef USE_DOUBLE // If enabled, double-precesion is used + +#define DATA_DEGREE 0 // The order of the B-Spline used to splat in data for color interpolation +#define WEIGHT_DEGREE 2 // The order of the B-Spline used to splat in the weights for density estimation +#define NORMAL_DEGREE 2 // The order of the B-Spline used to splat in the normals for constructing the Laplacian constraints +#define DEFAULT_FEM_DEGREE 1 // The default finite-element degree +#define DEFAULT_FEM_BOUNDARY BOUNDARY_NEUMANN // The default finite-element boundary type +#define DEFAULT_DIMENSION 3 // The dimension of the system + +#include +#include +#include +#include +#include "MyMiscellany.h" +#include "CmdLineParser.h" +#include "PPolynomial.h" +#include "FEMTree.h" +#include "Ply.h" +#include "VertexFactory.h" +#include "Image.h" +#include "RegularGrid.h" + +MessageWriter messageWriter; + +const float DefaultPointWeightMultiplier = 2.f; + +enum NormalType +{ + NORMALS_NONE , + NORMALS_SAMPLES , + NORMALS_GRADIENTS , + NORMALS_COUNT +}; +const char* NormalsNames[] = { "none" , "samples" , "gradients" }; + +cmdLineParameter< char* > + In( "in" ) , + Out( "out" ) , + TempDir( "tempDir" ) , + Grid( "grid" ) , + Tree( "tree" ) , + Envelope( "envelope" ) , + EnvelopeGrid( "envelopeGrid" ), + Transform( "xForm" ); + +cmdLineReadable + Performance( "performance" ) , + ShowResidual( "showResidual" ) , + NoComments( "noComments" ) , + PolygonMesh( "polygonMesh" ) , + NonManifold( "nonManifold" ) , + ASCII( "ascii" ) , + Density( "density" ) , + LinearFit( "linearFit" ) , + PrimalGrid( "primalGrid" ) , + ExactInterpolation( "exact" ) , + Colors( "colors" ) , + InCore( "inCore" ) , + NoDirichletErode( "noErode" ) , + Verbose( "verbose" ); + +cmdLineParameter< int > +#ifndef FAST_COMPILE + Degree( "degree" , DEFAULT_FEM_DEGREE ) , +#endif // !FAST_COMPILE + Depth( "depth" , 8 ) , + KernelDepth( "kernelDepth" ) , + SolveDepth( "solveDepth" ) , + EnvelopeDepth( "envelopeDepth" ) , + Iters( "iters" , 8 ) , + FullDepth( "fullDepth" , 5 ) , + BaseDepth( "baseDepth" ) , + BaseVCycles( "baseVCycles" , 1 ) , +#ifndef FAST_COMPILE + BType( "bType" , DEFAULT_FEM_BOUNDARY+1 ) , +#endif // !FAST_COMPILE + Normals( "normals" , NORMALS_NONE ) , + MaxMemoryGB( "maxMemory" , 0 ) , +#ifdef _OPENMP + ParallelType( "parallel" , (int)ThreadPool::OPEN_MP ) , +#else // !_OPENMP + ParallelType( "parallel" , (int)ThreadPool::THREAD_POOL ) , +#endif // _OPENMP + ScheduleType( "schedule" , (int)ThreadPool::DefaultSchedule ) , + ThreadChunkSize( "chunkSize" , (int)ThreadPool::DefaultChunkSize ) , + Threads( "threads" , (int)std::thread::hardware_concurrency() ); + +cmdLineParameter< float > + DataX( "data" , 32.f ) , + SamplesPerNode( "samplesPerNode" , 1.5f ) , + Scale( "scale" , 1.1f ) , + Width( "width" , 0.f ) , + Confidence( "confidence" , 0.f ) , + ConfidenceBias( "confidenceBias" , 0.f ) , + CGSolverAccuracy( "cgAccuracy" , 1e-3f ) , + LowDepthCutOff( "lowDepthCutOff" , 0.f ) , + PointWeight( "pointWeight" ); + +cmdLineReadable* params[] = +{ +#ifndef FAST_COMPILE + &Degree , &BType , +#endif // !FAST_COMPILE + &In , &Depth , &Out , &Transform , + &SolveDepth , + &Envelope , + &EnvelopeGrid , + &Width , + &Scale , &Verbose , &CGSolverAccuracy , &NoComments , + &KernelDepth , &SamplesPerNode , &Confidence , &NonManifold , &PolygonMesh , &ASCII , &ShowResidual , + &EnvelopeDepth , + &NoDirichletErode , + &ConfidenceBias , + &BaseDepth , &BaseVCycles , + &PointWeight , + &Grid , &Threads , + &Tree , + &Density , + &FullDepth , + &Iters , + &DataX , + &Colors , + &Normals , + &LinearFit , + &PrimalGrid , + &TempDir , + &ExactInterpolation , + &Performance , + &MaxMemoryGB , + &InCore , + &ParallelType , + &ScheduleType , + &ThreadChunkSize , + &LowDepthCutOff , + NULL +}; + +void ShowUsage(char* ex) +{ + printf( "Usage: %s\n" , ex ); + printf( "\t --%s \n" , In.name ); + printf( "\t[--%s \n" , Envelope.name ); + printf( "\t[--%s ]\n" , Out.name ); + printf( "\t[--%s ]\n" , Grid.name ); + printf( "\t[--%s \n" , EnvelopeGrid.name ); + printf( "\t[--%s ]\n" , Tree.name ); +#ifndef FAST_COMPILE + printf( "\t[--%s =%d]\n" , Degree.name , Degree.value ); + printf( "\t[--%s =%d]\n" , BType.name , BType.value ); + for( int i=0 ; i=%d]\n" , Depth.name , Depth.value ); + printf( "\t[--%s =%d]\n" , SolveDepth.name , SolveDepth.value ); + printf( "\t[--%s ]\n" , Width.name ); + printf( "\t[--%s =%d]\n" , FullDepth.name , FullDepth.value ); + printf( "\t[--%s =%d]\n" , EnvelopeDepth.name , EnvelopeDepth.value ); + printf( "\t[--%s ]\n" , BaseDepth.name ); + printf( "\t[--%s =%d]\n" , BaseVCycles.name , BaseVCycles.value ); + printf( "\t[--%s =%f]\n" , Scale.name , Scale.value ); + printf( "\t[--%s =%f]\n" , SamplesPerNode.name, SamplesPerNode.value ); + printf( "\t[--%s =%.3e * ]\n" , PointWeight.name , DefaultPointWeightMultiplier * DEFAULT_FEM_DEGREE ); + printf( "\t[--%s =%d]\n" , Iters.name , Iters.value ); + printf( "\t[--%s]\n" , ExactInterpolation.name ); + printf( "\t[--%s =%f]\n" , DataX.name , DataX.value ); + printf( "\t[--%s]\n" , Colors.name ); + printf( "\t[--%s =%d]\n" , Normals.name , Normals.value ); + for( int i=0 ; i=%d]\n" , Threads.name , Threads.value ); + printf( "\t[--%s =%d]\n" , ParallelType.name , ParallelType.value ); + for( size_t i=0 ; i=%d]\n" , ScheduleType.name , ScheduleType.value ); + for( size_t i=0 ; i=%d]\n" , ThreadChunkSize.name , ThreadChunkSize.value ); + printf( "\t[--%s =%f]\n" , LowDepthCutOff.name , LowDepthCutOff.value ); + printf( "\t[--%s =%f]\n" , Confidence.name , Confidence.value ); + printf( "\t[--%s =%f]\n" , ConfidenceBias.name , ConfidenceBias.value ); + printf( "\t[--%s]\n" , NonManifold.name ); + printf( "\t[--%s]\n" , PolygonMesh.name ); + printf( "\t[--%s =%g]\n" , CGSolverAccuracy.name , CGSolverAccuracy.value ); + printf( "\t[--%s =%d]\n" , MaxMemoryGB.name , MaxMemoryGB.value ); + printf( "\t[--%s]\n" , NoDirichletErode.name ); + printf( "\t[--%s]\n" , Performance.name ); + printf( "\t[--%s]\n" , Density.name ); + printf( "\t[--%s]\n" , LinearFit.name ); + printf( "\t[--%s]\n" , PrimalGrid.name ); + printf( "\t[--%s]\n" , ASCII.name ); + printf( "\t[--%s]\n" , NoComments.name ); + printf( "\t[--%s]\n" , TempDir.name ); + printf( "\t[--%s]\n" , InCore.name ); + printf( "\t[--%s]\n" , Verbose.name ); +} + +double Weight( double v , double start , double end ) +{ + v = ( v - start ) / ( end - start ); + if ( v<0 ) return 1.; + else if( v>1 ) return 0.; + else + { + // P(x) = a x^3 + b x^2 + c x + d + // P (0) = 1 , P (1) = 0 , P'(0) = 0 , P'(1) = 0 + // => d = 1 , a + b + c + d = 0 , c = 0 , 3a + 2b + c = 0 + // => c = 0 , d = 1 , a + b = -1 , 3a + 2b = 0 + // => a = 2 , b = -3 , c = 0 , d = 1 + // => P(x) = 2 x^3 - 3 x^2 + 1 + return 2. * v * v * v - 3. * v * v + 1.; + } +} + +template< unsigned int Dim , class Real > +struct FEMTreeProfiler +{ + double t; + + void start( void ){ t = Time() , FEMTree< Dim , Real >::ResetLocalMemoryUsage(); } + void print( const char* header ) const + { + FEMTree< Dim , Real >::MemoryUsage(); + if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); + else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); + } + void dumpOutput( const char* header ) const + { + FEMTree< Dim , Real >::MemoryUsage(); + if( header ) messageWriter( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); + else messageWriter( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); + } + void dumpOutput2( std::vector< std::string >& comments , const char* header ) const + { + FEMTree< Dim , Real >::MemoryUsage(); + if( header ) messageWriter( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); + else messageWriter( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %d (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() ); + } +}; + +template< class Real , unsigned int Dim > +XForm< Real , Dim+1 > GetBoundingBoxXForm( Point< Real , Dim > min , Point< Real , Dim > max , Real scaleFactor ) +{ + Point< Real , Dim > center = ( max + min ) / 2; + Real scale = max[0] - min[0]; + for( int d=1 ; d( scale , max[d]-min[d] ); + scale *= scaleFactor; + for( int i=0 ; i tXForm = XForm< Real , Dim+1 >::Identity() , sXForm = XForm< Real , Dim+1 >::Identity(); + for( int i=0 ; i +XForm< Real , Dim+1 > GetBoundingBoxXForm( Point< Real , Dim > min , Point< Real , Dim > max , Real width , Real scaleFactor , int& depth ) +{ + // Get the target resolution (along the largest dimension) + Real resolution = ( max[0]-min[0] ) / width; + for( int d=1 ; d( resolution , ( max[d]-min[d] ) / width ); + resolution *= scaleFactor; + depth = 0; + while( (1< center = ( max + min ) / 2; + Real scale = (1< tXForm = XForm< Real , Dim+1 >::Identity() , sXForm = XForm< Real , Dim+1 >::Identity(); + for( int i=0 ; i +using InputOrientedPointStreamInfo = typename FEMTreeInitializer< Dim , Real >::template InputPointStream< VectorTypeUnion< Real , typename VertexFactory::NormalFactory< Real , Dim >::VertexType , AuxData > >; + +template< typename Real , unsigned int Dim , typename AuxData > +using InputOrientedPointStream = typename InputOrientedPointStreamInfo< Real , Dim , AuxData >::StreamType; + +template< class Real , unsigned int Dim , typename AuxData > +XForm< Real , Dim+1 > GetPointXForm( InputOrientedPointStream< Real , Dim , AuxData > &stream , Real width , Real scaleFactor , int& depth ) +{ + Point< Real , Dim > min , max; + InputOrientedPointStreamInfo< Real , Dim , AuxData >::BoundingBox( stream , min , max ); + return GetBoundingBoxXForm( min , max , width , scaleFactor , depth ); +} +template< class Real , unsigned int Dim , typename AuxData > +XForm< Real , Dim+1 > GetPointXForm( InputOrientedPointStream< Real , Dim , AuxData > &stream , Real scaleFactor ) +{ + Point< Real , Dim > min , max; + InputOrientedPointStreamInfo< Real , Dim , AuxData >::BoundingBox( stream , min , max ); + return GetBoundingBoxXForm( min , max , scaleFactor ); +} + +template< unsigned int Dim , typename Real > +struct ConstraintDual +{ + Real target , weight; + ConstraintDual( Real t , Real w ) : target(t) , weight(w){ } + CumulativeDerivativeValues< Real , Dim , 0 > operator()( const Point< Real , Dim >& p ) const { return CumulativeDerivativeValues< Real , Dim , 0 >( target*weight ); }; +}; +template< unsigned int Dim , typename Real > +struct SystemDual +{ + Real weight; + SystemDual( Real w ) : weight(w){ } + CumulativeDerivativeValues< Real , Dim , 0 > operator()( const Point< Real , Dim >& p , const CumulativeDerivativeValues< Real , Dim , 0 >& dValues ) const { return dValues * weight; }; + CumulativeDerivativeValues< double , Dim , 0 > operator()( const Point< Real , Dim >& p , const CumulativeDerivativeValues< double , Dim , 0 >& dValues ) const { return dValues * weight; }; +}; +template< unsigned int Dim > +struct SystemDual< Dim , double > +{ + typedef double Real; + Real weight; + SystemDual( Real w ) : weight(w){ } + CumulativeDerivativeValues< Real , Dim , 0 > operator()( const Point< Real , Dim >& p , const CumulativeDerivativeValues< Real , Dim , 0 >& dValues ) const { return dValues * weight; }; +}; + +template< typename Real , typename SetVertexFunction , typename InputSampleDataType , typename VertexFactory , unsigned int ... FEMSigs > +void ExtractMesh +( + UIntPack< FEMSigs ... > , + FEMTree< sizeof ... ( FEMSigs ) , Real >& tree , + const DenseNodeData< Real , UIntPack< FEMSigs ... > >& solution , + Real isoValue , + const std::vector< typename FEMTree< sizeof ... ( FEMSigs ) , Real >::PointSample > *samples , + std::vector< InputSampleDataType > *sampleData , + const typename FEMTree< sizeof ... ( FEMSigs ) , Real >::template DensityEstimator< WEIGHT_DEGREE > *density , + const VertexFactory &vertexFactory , + const InputSampleDataType &zeroInputSampleDataType , + SetVertexFunction SetVertex , + std::vector< std::string > &comments , + XForm< Real , sizeof...(FEMSigs)+1 > unitCubeToModel +) +{ + static const int Dim = sizeof ... ( FEMSigs ); + typedef UIntPack< FEMSigs ... > Sigs; + typedef typename VertexFactory::VertexType Vertex; + + static const unsigned int DataSig = FEMDegreeAndBType< DATA_DEGREE , BOUNDARY_FREE >::Signature; + typedef typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE > DensityEstimator; + + FEMTreeProfiler< Dim , Real > profiler; + + char tempHeader[1024]; + { + char tempPath[1024]; + tempPath[0] = 0; + if( TempDir.set ) strcpy( tempPath , TempDir.value ); + else SetTempDirectory( tempPath , sizeof(tempPath) ); + if( strlen(tempPath)==0 ) sprintf( tempPath , ".%c" , FileSeparator ); + if( tempPath[ strlen( tempPath )-1 ]==FileSeparator ) sprintf( tempHeader , "%sPR_" , tempPath ); + else sprintf( tempHeader , "%s%cPR_" , tempPath , FileSeparator ); + } + CoredMeshData< Vertex , node_index_type > *mesh; + if( InCore.set ) mesh = new CoredVectorMeshData< Vertex , node_index_type >(); + else mesh = new CoredFileMeshData< node_index_type , VertexFactory >( vertexFactory , tempHeader ); + + profiler.start(); + typename IsoSurfaceExtractor< Dim , Real , Vertex >::IsoStats isoStats; + if( sampleData ) + { + SparseNodeData< ProjectiveData< InputSampleDataType , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setExtrapolatedDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL ); + for( const RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) ) + { + ProjectiveData< InputSampleDataType , Real >* clr = _sampleData( n ); + if( clr ) (*clr) *= (Real)pow( DataX.value , tree.depth( n ) ); + } + isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< InputSampleDataType >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , density , &_sampleData , solution , isoValue , *mesh , zeroInputSampleDataType , SetVertex , !LinearFit.set , Normals.value==NORMALS_GRADIENTS , !NonManifold.set , PolygonMesh.set , false ); + } +#if defined( __GNUC__ ) && __GNUC__ < 5 +#ifdef SHOW_WARNINGS +#warning "you've got me gcc version<5" +#endif // SHOW_WARNINGS + else isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< InputSampleDataType >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , density , (SparseNodeData< ProjectiveData< InputSampleDataType , Real > , IsotropicUIntPack< Dim , DataSig > > *)NULL , solution , isoValue , *mesh , zeroInputSampleDataType , SetVertex , !LinearFit.set , Normals.value==NORMALS_GRADIENTS , !NonManifold.set , PolygonMesh.set , false ); +#else // !__GNUC__ || __GNUC__ >=5 + else isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< InputSampleDataType >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , density , NULL , solution , isoValue , *mesh , zeroInputSampleDataType , SetVertex , !LinearFit.set , Normals.value==NORMALS_GRADIENTS , !NonManifold.set , PolygonMesh.set , false ); +#endif // __GNUC__ || __GNUC__ < 4 + messageWriter( "Vertices / Polygons: %llu / %llu\n" , (unsigned long long)( mesh->outOfCoreVertexNum()+mesh->inCoreVertices.size() ) , (unsigned long long)mesh->polygonNum() ); + std::string isoStatsString = isoStats.toString() + std::string( "\n" ); + messageWriter( isoStatsString.c_str() ); + if( PolygonMesh.set ) profiler.dumpOutput2( comments , "# Got polygons:" ); + else profiler.dumpOutput2( comments , "# Got triangles:" ); + + std::vector< std::string > noComments; + typename VertexFactory::Transform unitCubeToModelTransform( unitCubeToModel ); + PLY::WritePolygons< VertexFactory , node_index_type , Real , Dim >( Out.value , vertexFactory , mesh , ASCII.set ? PLY_ASCII : PLY_BINARY_NATIVE , NoComments.set ? noComments : comments , unitCubeToModelTransform ); + delete mesh; +} + +template< typename Real , unsigned int Dim > +void WriteGrid( const char *fileName , ConstPointer( Real ) values , unsigned int res , XForm< Real , Dim+1 > voxelToModel , bool verbose ) +{ + char *ext = GetFileExtension( fileName ); + + if( Dim==2 && ImageWriter::ValidExtension( ext ) ) + { + unsigned int totalResolution = 1; + for( int d=0 ; d avgs( ThreadPool::NumThreads() , 0 ); + ThreadPool::Parallel_for( 0 , totalResolution , [&]( unsigned int thread , size_t i ){ avgs[thread] += values[i]; } ); + for( unsigned int t=0 ; t stds( ThreadPool::NumThreads() , 0 ); + ThreadPool::Parallel_for( 0 , totalResolution , [&]( unsigned int thread , size_t i ){ stds[thread] += ( values[i] - avg ) * ( values[i] - avg ); } ); + for( unsigned int t=0 ; t [0,255]\n" , avg - 2*std , avg + 2*std ); + printf( "Transform:\n" ); + for( int i=0 ; i( (Real)1. , std::max< Real >( (Real)-1. , ( values[i] - avg ) / (2*std ) ) ); + v = (Real)( ( v + 1. ) / 2. * 256. ); + unsigned char color = (unsigned char )std::min< Real >( (Real)255. , std::max< Real >( (Real)0. , v ) ); + for( int c=0 ; c<3 ; c++ ) pixels[i*3+c ] = color; + } + ); + ImageWriter::Write( fileName , pixels , res , res , 3 ); + delete[] pixels; + } + else if( !strcasecmp( ext , "iso" ) ) + { + FILE *fp = fopen( fileName , "wb" ); + if( !fp ) ERROR_OUT( "Failed to open file for writing: " , fileName ); + int r = (int)res; + fwrite( &r , sizeof(int) , 1 , fp ); + size_t count = 1; + for( unsigned int d=0 ; d::Write( fileName , _res , values , voxelToModel ); + } + delete[] ext; +} + +template< class Real , typename AuxDataFactory , unsigned int ... FEMSigs > +void Execute( UIntPack< FEMSigs ... > , const AuxDataFactory &auxDataFactory ) +{ + static const int Dim = sizeof ... ( FEMSigs ); + typedef UIntPack< FEMSigs ... > Sigs; + typedef UIntPack< FEMSignature< FEMSigs >::Degree ... > Degrees; + typedef UIntPack< FEMDegreeAndBType< NORMAL_DEGREE , DerivativeBoundary< FEMSignature< FEMSigs >::BType , 1 >::BType >::Signature ... > NormalSigs; + static const unsigned int DataSig = FEMDegreeAndBType< DATA_DEGREE , BOUNDARY_FREE >::Signature; + typedef typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE > DensityEstimator; + typedef typename FEMTree< Dim , Real >::template InterpolationInfo< Real , 0 > InterpolationInfo; + using namespace VertexFactory; + + // The factory for constructing an input sample + typedef Factory< Real , PositionFactory< Real , Dim > , Factory< Real , NormalFactory< Real , Dim > , AuxDataFactory > > InputSampleFactory; + + // The factory for constructing an input sample's data + typedef Factory< Real , NormalFactory< Real , Dim > , AuxDataFactory > InputSampleDataFactory; + + // The input point stream information: First piece of data is the normal; the remainder is the auxiliary data + typedef InputOrientedPointStreamInfo< Real , Dim , typename AuxDataFactory::VertexType > InputPointStreamInfo; + + // The type of the input sample + typedef typename InputPointStreamInfo::PointAndDataType InputSampleType; + + // The type of the input sample's data + typedef typename InputPointStreamInfo::DataType InputSampleDataType; + + typedef InputDataStream< InputSampleType > InputPointStream; + typedef TransformedInputDataStream< InputSampleType > XInputPointStream; + + InputSampleFactory inputSampleFactory( PositionFactory< Real , Dim >() , InputSampleDataFactory( NormalFactory< Real , Dim >() , auxDataFactory ) ); + InputSampleDataFactory inputSampleDataFactory( NormalFactory< Real , Dim >() , auxDataFactory ); + + typedef RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type > FEMTreeNode; + typedef typename FEMTreeInitializer< Dim , Real >::GeometryNodeType GeometryNodeType; + std::vector< std::string > comments; + messageWriter( comments , "*************************************************************\n" ); + messageWriter( comments , "*************************************************************\n" ); + messageWriter( comments , "** Running Screened Poisson Reconstruction (Version %s) **\n" , VERSION ); + messageWriter( comments , "*************************************************************\n" ); + messageWriter( comments , "*************************************************************\n" ); + if( !Threads.set ) messageWriter( comments , "Running with %d threads\n" , Threads.value ); + + bool needNormalData = DataX.value>0 && Normals.value; + bool needAuxData = DataX.value>0 && auxDataFactory.bufferSize(); + + XForm< Real , Dim+1 > modelToUnitCube , unitCubeToModel; + if( Transform.set ) + { + FILE* fp = fopen( Transform.value , "r" ); + if( !fp ) + { + WARN( "Could not read x-form from: " , Transform.value ); + modelToUnitCube = XForm< Real , Dim+1 >::Identity(); + } + else + { + for( int i=0 ; i::Identity(); + + char str[1024]; + for( int i=0 ; params[i] ; i++ ) + if( params[i]->set ) + { + params[i]->writeValue( str ); + if( strlen( str ) ) messageWriter( comments , "\t--%s %s\n" , params[i]->name , str ); + else messageWriter( comments , "\t--%s\n" , params[i]->name ); + } + + double startTime = Time(); + Real isoValue = 0; + + FEMTree< Dim , Real > tree( MEMORY_ALLOCATOR_BLOCK_SIZE ); + FEMTreeProfiler< Dim , Real > profiler; + + if( Depth.set && Width.value>0 ) + { + WARN( "Both --" , Depth.name , " and --" , Width.name , " set, ignoring --" , Width.name ); + Width.value = 0; + } + + size_t pointCount; + + Real pointWeightSum; + std::vector< typename FEMTree< Dim , Real >::PointSample >* samples = new std::vector< typename FEMTree< Dim , Real >::PointSample >(); + DenseNodeData< GeometryNodeType , IsotropicUIntPack< Dim , FEMTrivialSignature > > geometryNodeDesignators; + std::vector< InputSampleDataType >* sampleData = NULL; + DensityEstimator* density = NULL; + SparseNodeData< Point< Real , Dim > , NormalSigs >* normalInfo = NULL; + Real targetValue = (Real)0.5; + + // Read in the samples (and color data) + { + profiler.start(); + InputPointStream* pointStream; + char* ext = GetFileExtension( In.value ); + sampleData = new std::vector< InputSampleDataType >(); + std::vector< InputSampleType > inCorePoints; + if( InCore.set ) + { + InputPointStream *_pointStream; + if ( !strcasecmp( ext , "bnpts" ) ) _pointStream = new BinaryInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); + else if( !strcasecmp( ext , "ply" ) ) _pointStream = new PLYInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); + else _pointStream = new ASCIIInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); + InputSampleType p; + while( _pointStream->next( p ) ) inCorePoints.push_back( p ); + delete _pointStream; + + pointStream = new MemoryInputDataStream< InputSampleType >( inCorePoints.size() , &inCorePoints[0] ); + } + else + { + if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); + else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); + else pointStream = new ASCIIInputDataStream< InputSampleFactory >( In.value , inputSampleFactory ); + } + delete[] ext; + + typename InputSampleFactory::Transform _modelToUnitCube( modelToUnitCube ); + auto XFormFunctor = [&]( InputSampleType &p ){ p = _modelToUnitCube( p ); }; + XInputPointStream _pointStream( XFormFunctor , *pointStream ); + if( Width.value>0 ) modelToUnitCube = GetPointXForm< Real , Dim , typename AuxDataFactory::VertexType >( _pointStream , Width.value , (Real)( Scale.value>0 ? Scale.value : 1. ) , Depth.value ) * modelToUnitCube; + else modelToUnitCube = Scale.value>0 ? GetPointXForm< Real , Dim , typename AuxDataFactory::VertexType >( _pointStream , (Real)Scale.value ) * modelToUnitCube : modelToUnitCube; + + if( !SolveDepth.set ) SolveDepth.value = Depth.value; + if( SolveDepth.value>Depth.value ) + { + WARN( "Solution depth cannot exceed system depth: " , SolveDepth.value , " <= " , Depth.value ); + SolveDepth.value = Depth.value; + } + + { + typename InputSampleFactory::Transform _modelToUnitCube( modelToUnitCube ); + auto XFormFunctor = [&]( InputSampleType &p ){ p = _modelToUnitCube( p ); }; + XInputPointStream _pointStream( XFormFunctor , *pointStream ); + auto ProcessDataWithConfidence = [&]( const Point< Real , Dim > &p , typename InputPointStreamInfo::DataType &d ) + { + Real l = (Real)Length( d.template get<0>() ); + if( !l || !std::isfinite( l ) ) return (Real)-1.; + return (Real)pow( l , Confidence.value ); + }; + auto ProcessData = []( const Point< Real , Dim > &p , typename InputPointStreamInfo::DataType &d ) + { + Real l = (Real)Length( d.template get<0>() ); + if( !l || !std::isfinite( l ) ) return (Real)-1.; + d.template get<0>() /= l; + return (Real)1.; + }; + if( Confidence.value>0 ) pointCount = FEMTreeInitializer< Dim , Real >::template Initialize< InputSampleDataType >( tree.spaceRoot() , _pointStream , Depth.value , *samples , *sampleData , true , tree.nodeAllocators[0] , tree.initializer() , ProcessDataWithConfidence ); + else pointCount = FEMTreeInitializer< Dim , Real >::template Initialize< InputSampleDataType >( tree.spaceRoot() , _pointStream , Depth.value , *samples , *sampleData , true , tree.nodeAllocators[0] , tree.initializer() , ProcessData ); + } + + unitCubeToModel = modelToUnitCube.inverse(); + delete pointStream; + + messageWriter( "Input Points / Samples: %llu / %llu\n" , (unsigned long long)pointCount , (unsigned long long)samples->size() ); + profiler.dumpOutput2( comments , "# Read input into tree:" ); + } + + DenseNodeData< Real , Sigs > solution; + { + DenseNodeData< Real , Sigs > constraints; + InterpolationInfo* iInfo = NULL; + int solveDepth = Depth.value; + + tree.resetNodeIndices( 0 , std::make_tuple() ); + + // Get the kernel density estimator + { + profiler.start(); + density = tree.template setDensityEstimator< 1 , WEIGHT_DEGREE >( *samples , KernelDepth.value , SamplesPerNode.value ); + profiler.dumpOutput2( comments , "# Got kernel density:" ); + } + + // Transform the Hermite samples into a vector field + { + profiler.start(); + normalInfo = new SparseNodeData< Point< Real , Dim > , NormalSigs >(); + std::function< bool ( InputSampleDataType , Point< Real , Dim >& ) > ConversionFunction = []( InputSampleDataType in , Point< Real , Dim > &out ) + { + Point< Real , Dim > n = in.template get<0>(); + Real l = (Real)Length( n ); + // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals... + if( !l ) return false; + out = n / l; + return true; + }; + std::function< bool ( InputSampleDataType , Point< Real , Dim >& , Real & ) > ConversionAndBiasFunction = []( InputSampleDataType in , Point< Real , Dim > &out , Real &bias ) + { + Point< Real , Dim > n = in.template get<0>(); + Real l = (Real)Length( n ); + // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals... + if( !l ) return false; + out = n / l; + bias = (Real)( log( l ) * ConfidenceBias.value / log( 1<<(Dim-1) ) ); + return true; + }; +#if 1 + if( ConfidenceBias.value>0 ) *normalInfo = tree.setInterpolatedDataField( NormalSigs() , *samples , *sampleData , density , BaseDepth.value , Depth.value , (Real)LowDepthCutOff.value , pointWeightSum , ConversionAndBiasFunction ); + else *normalInfo = tree.setInterpolatedDataField( NormalSigs() , *samples , *sampleData , density , BaseDepth.value , Depth.value , (Real)LowDepthCutOff.value , pointWeightSum , ConversionFunction ); +#else + if( ConfidenceBias.value>0 ) *normalInfo = tree.setInterpolatedDataField( NormalSigs() , *samples , *sampleData , density , 0 , Depth.value , (Real)LowDepthCutOff.value , pointWeightSum , ConversionAndBiasFunction ); + else *normalInfo = tree.setInterpolatedDataField( NormalSigs() , *samples , *sampleData , density , 0 , Depth.value , (Real)LowDepthCutOff.value , pointWeightSum , ConversionFunction ); +#endif + ThreadPool::Parallel_for( 0 , normalInfo->size() , [&]( unsigned int , size_t i ){ (*normalInfo)[i] *= (Real)-1.; } ); + profiler.dumpOutput2( comments , "# Got normal field:" ); + messageWriter( "Point weight / Estimated Measure: %g / %g\n" , pointWeightSum , pointCount*pointWeightSum ); + } + + // Get the geometry designators indicating if the space node are interior to, exterior to, or contain the envelope boundary + if( Envelope.set ) + { + profiler.start(); + { + // Make the octree complete up to the base depth + std::function< void ( FEMTreeNode * , unsigned int ) > MakeComplete = [&]( FEMTreeNode *node , unsigned int depth ) + { + if( node->depth()<(int)depth ) + { + if( !node->children ) node->template initChildren< false >( tree.nodeAllocators.size() ? tree.nodeAllocators[0] : NULL , tree.initializer() ); + for( int c=0 ; c<(1<children+c , depth ); + } + }; + MakeComplete( &tree.spaceRoot() , BaseDepth.value ); + + // Read in the envelope geometry + std::vector< Point< Real , Dim > > vertices; + std::vector< SimplexIndex< Dim-1 , node_index_type > > simplices; + { + std::vector< typename PositionFactory< Real , Dim >::VertexType > _vertices; + std::vector< std::vector< int > > polygons; + std::vector< std::string > comments; + int file_type; + PLY::ReadPolygons( Envelope.value , PositionFactory< Real , Dim >() , _vertices , polygons , file_type , comments ); + vertices.resize( _vertices.size() ); + for( int i=0 ; i::template GetGeometryNodeDesignators( &tree.spaceRoot() , vertices , simplices , BaseDepth.value , EnvelopeDepth.value , tree.nodeAllocators , tree.initializer() ); + + // Make nodes in the support of the vector field @{ExactDepth} interior + if( !NoDirichletErode.set ) + { + // What to do if we find a node in the support of the vector field + auto SetScratchFlag = [&]( FEMTreeNode *node ) + { + if( node ) + { + while( node->depth()>BaseDepth.value ) node = node->parent; + node->nodeData.setScratchFlag( true ); + } + }; + + std::function< void ( FEMTreeNode * ) > PropagateToLeaves = [&]( const FEMTreeNode *node ) + { + geometryNodeDesignators[ node ] = GeometryNodeType::INTERIOR; + if( node->children ) for( int c=0 ; c<(1<children+c ); + }; + + // Flags indicating if a node contains a non-zero vector field coefficient + std::vector< bool > isVectorFieldElement( tree.nodeCount() , false ); + + // Get the set of base nodes + std::vector< FEMTreeNode * > baseNodes; + auto TerminationLambda = []( const FEMTreeNode *node ){ return node->depth()==BaseDepth.value; }; + for( FEMTreeNode *node=tree.spaceRoot().nextNode( TerminationLambda , NULL ) ; node ; node=tree.spaceRoot().nextNode( TerminationLambda , node ) ) if( node->depth()==BaseDepth.value ) baseNodes.push_back( node ); + + std::vector< node_index_type > vectorFieldElementCounts( baseNodes.size() ); + for( int i=0 ; inextNode() ; node ; node=baseNodes[i]->nextNode(node) ) + { + Point< Real , Dim > *n = (*normalInfo)( node ); + if( n && Point< Real , Dim >::SquareNorm( *n ) ) isVectorFieldElement[ node->nodeData.nodeIndex ] = true , vectorFieldElementCounts[i]++; + } + } ); + size_t vectorFieldElementCount = 0; + for( int i=0 ; i vectorFieldElements; + vectorFieldElements.reserve( vectorFieldElementCount ); + { + std::vector< std::vector< FEMTreeNode * > > _vectorFieldElements( baseNodes.size() ); + for( int i=0 ; i<_vectorFieldElements.size() ; i++ ) _vectorFieldElements[i].reserve( vectorFieldElementCounts[i] ); + ThreadPool::Parallel_for( 0 , baseNodes.size() , [&]( unsigned int t , size_t i ) + { + for( FEMTreeNode *node=baseNodes[i]->nextNode() ; node ; node=baseNodes[i]->nextNode(node) ) + { + if( isVectorFieldElement[ node->nodeData.nodeIndex ] ) _vectorFieldElements[i].push_back( node ); + node->nodeData.setScratchFlag( false ); + } + } ); + for( int i=0 ; i<_vectorFieldElements.size() ; i++ ) vectorFieldElements.insert( vectorFieldElements.end() , _vectorFieldElements[i].begin() , _vectorFieldElements[i].end() ); + } + + // Set the scratch flag for the base nodes on which the vector field is supported +#ifdef SHOW_WARNINGS +#pragma message( "[WARNING] In principal, we should unlock finite elements whose support overlaps the vector field" ) +#endif // SHOW_WARNINGS + tree.template processNeighboringLeaves< -BSplineSupportSizes< NORMAL_DEGREE >::SupportStart , BSplineSupportSizes< NORMAL_DEGREE >::SupportEnd >( &vectorFieldElements[0] , vectorFieldElements.size() , SetScratchFlag , false ); + + // Set sub-trees rooted at interior nodes @ ExactDepth to interior + ThreadPool::Parallel_for( 0 , baseNodes.size() , [&]( unsigned int , size_t i ){ if( baseNodes[i]->nodeData.getScratchFlag() ) PropagateToLeaves( baseNodes[i] ); } ); + + // Adjust the coarser node designators in case exterior nodes have become boundary. + ThreadPool::Parallel_for( 0 , baseNodes.size() , [&]( unsigned int , size_t i ){ FEMTreeInitializer< Dim , Real >::PullGeometryNodeDesignatorsFromFiner( baseNodes[i] , geometryNodeDesignators ); } ); + FEMTreeInitializer< Dim , Real >::PullGeometryNodeDesignatorsFromFiner( &tree.spaceRoot() , geometryNodeDesignators , BaseDepth.value ); + } + } + profiler.dumpOutput2( comments , "# Initialized envelope constraints:" ); + } + + if( !Density.set ) delete density , density = NULL; + if( !needNormalData && !needAuxData ) delete sampleData , sampleData = NULL; + + // Add the interpolation constraints + if( PointWeight.value>0 ) + { + profiler.start(); + if( ExactInterpolation.set ) iInfo = FEMTree< Dim , Real >::template InitializeExactPointInterpolationInfo< Real , 0 > ( tree , *samples , ConstraintDual< Dim , Real >( targetValue , (Real)PointWeight.value * pointWeightSum ) , SystemDual< Dim , Real >( (Real)PointWeight.value * pointWeightSum ) , true , false ); + else iInfo = FEMTree< Dim , Real >::template InitializeApproximatePointInterpolationInfo< Real , 0 > ( tree , *samples , ConstraintDual< Dim , Real >( targetValue , (Real)PointWeight.value * pointWeightSum ) , SystemDual< Dim , Real >( (Real)PointWeight.value * pointWeightSum ) , true , 1 ); + profiler.dumpOutput2( comments , "#Initialized point interpolation constraints:" ); + } + + // Trim the tree and prepare for multigrid + { + profiler.start(); + constexpr int MAX_DEGREE = NORMAL_DEGREE > Degrees::Max() ? NORMAL_DEGREE : Degrees::Max(); + typename FEMTree< Dim , Real >::template HasNormalDataFunctor< NormalSigs > hasNormalDataFunctor( *normalInfo ); + auto hasDataFunctor = [&]( const FEMTreeNode *node ){ return hasNormalDataFunctor( node ); }; + if( geometryNodeDesignators.size() ) tree.template finalizeForMultigrid< MAX_DEGREE , Degrees::Max() >( BaseDepth.value , FullDepth.value , hasDataFunctor , [&]( const FEMTreeNode *node ){ return node->nodeData.nodeIndex<(node_index_type)geometryNodeDesignators.size() && geometryNodeDesignators[node]==GeometryNodeType::EXTERIOR; } , std::make_tuple( iInfo ) , std::make_tuple( normalInfo , density , &geometryNodeDesignators ) ); + else tree.template finalizeForMultigrid< MAX_DEGREE , Degrees::Max() >( BaseDepth.value , FullDepth.value , hasDataFunctor , []( const FEMTreeNode * ){ return false; } , std::make_tuple( iInfo ) , std::make_tuple( normalInfo , density ) ); + + if( geometryNodeDesignators.size() && EnvelopeGrid.set ) + { + FEMTreeInitializer< Dim , Real >::PushGeometryNodeDesignatorsToFiner( &tree.spaceRoot() , geometryNodeDesignators ); + FEMTreeInitializer< Dim , Real >::PullGeometryNodeDesignatorsFromFiner( &tree.spaceRoot() , geometryNodeDesignators ); + + auto WriteEnvelopeGrid = [&]( bool showFinest ) + { + int res = 0; + DenseNodeData< Real , IsotropicUIntPack< Dim , FEMTrivialSignature > > coefficients = tree.initDenseNodeData( IsotropicUIntPack< Dim , FEMTrivialSignature >() ); + for( const FEMTreeNode* n = tree.spaceRoot().nextNode() ; n ; n=tree.spaceRoot().nextNode( n ) ) + if( n->nodeData.nodeIndex!=-1 && ( ( showFinest && !n->children ) || ( !showFinest && geometryNodeDesignators[n->parent]==GeometryNodeType::BOUNDARY ) ) ) + { +#if 0 // Randomize the colors + auto Value = []( double v , double eps ){ return (Real)( v + Random< double >() * 2. * eps - eps ); }; + // Show the octree structure + if ( geometryNodeDesignators[n]==GeometryNodeType::INTERIOR ) coefficients[n] = Value( 0.75 , 0.25 ); + else if( geometryNodeDesignators[n]==GeometryNodeType::EXTERIOR ) coefficients[n] = Value( -0.75 , 0.25 ); + +#else + if ( geometryNodeDesignators[n]==GeometryNodeType::INTERIOR ) coefficients[n] = (Real) 1.; + else if( geometryNodeDesignators[n]==GeometryNodeType::EXTERIOR ) coefficients[n] = (Real)-1.; +#endif + } + Pointer( Real ) values = tree.template regularGridEvaluate< true >( coefficients , res , -1 , false ); + XForm< Real , Dim+1 > voxelToUnitCube = XForm< Real , Dim+1 >::Identity(); + for( int d=0 ; d( EnvelopeGrid.value , values , res , unitCubeToModel * voxelToUnitCube , Verbose.set ); + DeletePointer( values ); + }; + + WriteEnvelopeGrid( true ); + } + + profiler.dumpOutput2( comments , "# Finalized tree:" ); + } + // Add the FEM constraints + { + profiler.start(); + constraints = tree.initDenseNodeData( Sigs() ); + typename FEMIntegrator::template Constraint< Sigs , IsotropicUIntPack< Dim , 1 > , NormalSigs , IsotropicUIntPack< Dim , 0 > , Dim > F; + unsigned int derivatives2[Dim]; + for( int d=0 ; d Derivatives1; + typedef IsotropicUIntPack< Dim , 0 > Derivatives2; + for( int d=0 ; d::Index( derivatives1 ) ][ TensorDerivatives< Derivatives2 >::Index( derivatives2 ) ] = 1; + } + tree.addFEMConstraints( F , *normalInfo , constraints , solveDepth ); + profiler.dumpOutput2( comments , "# Set FEM constraints:" ); + } + + // Free up the normal info + delete normalInfo , normalInfo = NULL; + + // Add the interpolation constraints + if( PointWeight.value>0 ) + { + profiler.start(); + tree.addInterpolationConstraints( constraints , solveDepth , std::make_tuple( iInfo ) ); + profiler.dumpOutput2( comments , "#Set point constraints:" ); + } + + messageWriter( "Leaf Nodes / Active Nodes / Ghost Nodes / Dirichlet Supported Nodes: %llu / %llu / %llu / %llu\n" , (unsigned long long)tree.leaves() , (unsigned long long)tree.nodes() , (unsigned long long)tree.ghostNodes() , (unsigned long long)tree.dirichletElements() ); + messageWriter( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) ); + + // Solve the linear system + { + profiler.start(); + typename FEMTree< Dim , Real >::SolverInfo sInfo; + sInfo.cgDepth = 0 , sInfo.cascadic = true , sInfo.vCycles = 1 , sInfo.iters = Iters.value , sInfo.cgAccuracy = CGSolverAccuracy.value , sInfo.verbose = Verbose.set , sInfo.showResidual = ShowResidual.set , sInfo.showGlobalResidual = SHOW_GLOBAL_RESIDUAL_NONE , sInfo.sliceBlockSize = 1; + sInfo.baseVCycles = BaseVCycles.value; + typename FEMIntegrator::template System< Sigs , IsotropicUIntPack< Dim , 1 > > F( { 0. , 1. } ); + solution = tree.solveSystem( Sigs() , F , constraints , SolveDepth.value , sInfo , std::make_tuple( iInfo ) ); + profiler.dumpOutput2( comments , "# Linear system solved:" ); + if( iInfo ) delete iInfo , iInfo = NULL; + } + } + + { + profiler.start(); + double valueSum = 0 , weightSum = 0; + typename FEMTree< Dim , Real >::template MultiThreadedEvaluator< Sigs , 0 > evaluator( &tree , solution ); + std::vector< double > valueSums( ThreadPool::NumThreads() , 0 ) , weightSums( ThreadPool::NumThreads() , 0 ); + ThreadPool::Parallel_for( 0 , samples->size() , [&]( unsigned int thread , size_t j ) + { + ProjectiveData< Point< Real , Dim > , Real >& sample = (*samples)[j].sample; + Real w = sample.weight; + if( w>0 ) weightSums[thread] += w , valueSums[thread] += evaluator.values( sample.data / sample.weight , thread , (*samples)[j].node )[0] * w; + } + ); + for( size_t t=0 ; t::WriteParameter( fp ); + DenseNodeData< Real , Sigs >::WriteSignatures( fp ); + tree.write( fp , modelToUnitCube ); + solution.write( fp ); + if( sampleData ) + { + SparseNodeData< ProjectiveData< InputSampleDataType , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setExtrapolatedDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL ); + for( const RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) ) + { + ProjectiveData< InputSampleDataType , Real >* clr = _sampleData( n ); + if( clr ) (*clr) *= (Real)pow( DataX.value , tree.depth( n ) ); + } + _sampleData.write( fp ); + } + if( density ) density->write( fp ); + fclose( fp ); + } + + if( Grid.set ) + { + int res = 0; + profiler.start(); + Pointer( Real ) values = tree.template regularGridEvaluate< true >( solution , res , -1 , PrimalGrid.set ); + profiler.dumpOutput( "Got grid:" ); + XForm< Real , Dim+1 > voxelToUnitCube = XForm< Real , Dim+1 >::Identity(); + if( PrimalGrid.set ) for( int d=0 ; d( Grid.value , values , res , unitCubeToModel * voxelToUnitCube , Verbose.set ); + DeletePointer( values ); + } + + if( Out.set ) + { + if( Normals.value ) + { + if( Density.set ) + { + typedef Factory< Real , PositionFactory< Real , Dim > , NormalFactory< Real , Dim > , ValueFactory< Real > , AuxDataFactory > VertexFactory; + VertexFactory vertexFactory( PositionFactory< Real , Dim >() , NormalFactory< Real , Dim >() , ValueFactory< Real >() , auxDataFactory ); + if( Normals.value==NORMALS_SAMPLES ) + { + auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = d.template get<0>() , v.template get<2>() = w , v.template get<3>() = d.template get<1>(); }; + ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); + } + else if( Normals.value==NORMALS_GRADIENTS ) + { + auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = -g/(1<() = w , v.template get<3>() = d.template get<1>(); }; + ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); + } + } + else + { + typedef Factory< Real , PositionFactory< Real , Dim > , NormalFactory< Real , Dim > , AuxDataFactory > VertexFactory; + VertexFactory vertexFactory( PositionFactory< Real , Dim >() , NormalFactory< Real , Dim >() , auxDataFactory ); + if( Normals.value==NORMALS_SAMPLES ) + { + auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = d.template get<0>() , v.template get<2>() = d.template get<1>(); }; + ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); + } + else if( Normals.value==NORMALS_GRADIENTS ) + { + auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = -g/(1<() = d.template get<1>(); }; + ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); + } + } + } + else + { + if( Density.set ) + { + typedef Factory< Real , PositionFactory< Real , Dim > , ValueFactory< Real > , AuxDataFactory > VertexFactory; + VertexFactory vertexFactory( PositionFactory< Real , Dim >() , ValueFactory< Real >() , auxDataFactory ); + auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = w , v.template get<2>() = d.template get<1>(); }; + ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); + } + else + { + typedef Factory< Real , PositionFactory< Real , Dim > , AuxDataFactory > VertexFactory; + VertexFactory vertexFactory( PositionFactory< Real , Dim >() , auxDataFactory ); + auto SetVertex = []( typename VertexFactory::VertexType &v , Point< Real , Dim > p , Point< Real , Dim > g , Real w , InputSampleDataType d ){ v.template get<0>() = p , v.template get<1>() = d.template get<1>(); }; + ExtractMesh( UIntPack< FEMSigs ... >() , tree , solution , isoValue , samples , sampleData , density , vertexFactory , inputSampleDataFactory() , SetVertex , comments , unitCubeToModel ); + } + } + if( sampleData ){ delete sampleData ; sampleData = NULL; } + } + if( density ) delete density , density = NULL; + messageWriter( comments , "# Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-startTime , FEMTree< Dim , Real >::MaxMemoryUsage() ); +} + +#ifndef FAST_COMPILE +template< unsigned int Dim , class Real , BoundaryType BType , typename AuxDataFactory > +void Execute( const AuxDataFactory &auxDataFactory ) +{ + switch( Degree.value ) + { + case 1: return Execute< Real >( IsotropicUIntPack< Dim , FEMDegreeAndBType< 1 , BType >::Signature >() , auxDataFactory ); + case 2: return Execute< Real >( IsotropicUIntPack< Dim , FEMDegreeAndBType< 2 , BType >::Signature >() , auxDataFactory ); +// case 3: return Execute< Real >( IsotropicUIntPack< Dim , FEMDegreeAndBType< 3 , BType >::Signature >() , auxDataFactory ); +// case 4: return Execute< Real >( IsotropicUIntPack< Dim , FEMDegreeAndBType< 4 , BType >::Signature >() , auxDataFactory ); + default: ERROR_OUT( "Only B-Splines of degree 1 - 2 are supported" ); + } +} + +template< unsigned int Dim , class Real , typename AuxDataFactory > +void Execute( const AuxDataFactory &auxDataFactory ) +{ + switch( BType.value ) + { + case BOUNDARY_FREE+1: return Execute< Dim , Real , BOUNDARY_FREE >( auxDataFactory ); + case BOUNDARY_NEUMANN+1: return Execute< Dim , Real , BOUNDARY_NEUMANN >( auxDataFactory ); + case BOUNDARY_DIRICHLET+1: return Execute< Dim , Real , BOUNDARY_DIRICHLET >( auxDataFactory ); + default: ERROR_OUT( "Not a valid boundary type: " , BType.value ); + } +} +#endif // !FAST_COMPILE + + +int PoissonReconLib( int argc , char* argv[] ) +{ + Timer timer; +#ifdef USE_SEG_FAULT_HANDLER + WARN( "using seg-fault handler" ); + StackTracer::exec = argv[0]; + signal( SIGSEGV , SignalHandler ); +#endif // USE_SEG_FAULT_HANDLER +#ifdef ARRAY_DEBUG + WARN( "Array debugging enabled" ); +#endif // ARRAY_DEBUG + cmdLineParse( argc-1 , &argv[1] , params ); + if( MaxMemoryGB.value>0 ) SetPeakMemoryMB( MaxMemoryGB.value<<10 ); + ThreadPool::DefaultChunkSize = ThreadChunkSize.value; + ThreadPool::DefaultSchedule = (ThreadPool::ScheduleType)ScheduleType.value; + ThreadPool::Init( (ThreadPool::ParallelType)ParallelType.value , Threads.value ); + + messageWriter.echoSTDOUT = Verbose.set; + if( !In.set ) + { + ShowUsage( argv[0] ); + return 0; + } + + if( !BaseDepth.set ) BaseDepth.value = FullDepth.value; + if( BaseDepth.value>FullDepth.value ) + { + if( BaseDepth.set ) WARN( "Base depth must be smaller than full depth: " , BaseDepth.value , " <= " , FullDepth.value ); + BaseDepth.value = FullDepth.value; + } + if( !KernelDepth.set ) KernelDepth.value = Depth.value-2; + if( KernelDepth.value>Depth.value ) + { + WARN( "Kernel depth should not exceed depth: " , KernelDepth.name , " <= " , KernelDepth.value ); + KernelDepth.value = Depth.value; + } + + if( !EnvelopeDepth.set ) EnvelopeDepth.value = BaseDepth.value; + if( EnvelopeDepth.value>Depth.value ) + { + WARN( EnvelopeDepth.name , " can't be greater than " , Depth.name , ": " , EnvelopeDepth.value , " <= " , Depth.value ); + EnvelopeDepth.value = Depth.value; + } + if( EnvelopeDepth.value= " , BaseDepth.value ); + EnvelopeDepth.value = BaseDepth.value; + } + +#ifdef USE_DOUBLE + typedef double Real; +#else // !USE_DOUBLE + typedef float Real; +#endif // USE_DOUBLE + +#ifdef FAST_COMPILE + static const int Degree = DEFAULT_FEM_DEGREE; + static const BoundaryType BType = DEFAULT_FEM_BOUNDARY; + typedef IsotropicUIntPack< DEFAULT_DIMENSION , FEMDegreeAndBType< Degree , BType >::Signature > FEMSigs; + WARN( "Compiled for degree-" , Degree , ", boundary-" , BoundaryNames[ BType ] , ", " , sizeof(Real)==4 ? "single" : "double" , "-precision _only_" ); + if( !PointWeight.set ) PointWeight.value = DefaultPointWeightMultiplier*Degree; + char *ext = GetFileExtension( In.value ); + if( !strcasecmp( ext , "ply" ) ) + { + typedef VertexFactory::Factory< Real , typename VertexFactory::PositionFactory< Real , DEFAULT_DIMENSION > , typename VertexFactory::NormalFactory< Real , DEFAULT_DIMENSION > > Factory; + Factory factory; + bool *readFlags = new bool[ factory.plyReadNum() ]; + std::vector< PlyProperty > unprocessedProperties; + PLY::ReadVertexHeader( In.value , factory , readFlags , unprocessedProperties ); + if( !factory.plyValidReadProperties<0>( readFlags ) ) ERROR_OUT( "Ply file does not contain positions" ); + if( !factory.plyValidReadProperties<1>( readFlags ) ) ERROR_OUT( "Ply file does not contain normals" ); + delete[] readFlags; + + if( unprocessedProperties.size() ) Execute< Real >( FEMSigs() , VertexFactory::DynamicFactory< Real >( unprocessedProperties ) ); + else Execute< Real >( FEMSigs() , VertexFactory::EmptyFactory< Real >() ); + } + else + { + if( Colors.set ) Execute< Real >( FEMSigs() , VertexFactory::RGBColorFactory< Real >() ); + else Execute< Real >( FEMSigs() , VertexFactory::EmptyFactory< Real >() ); + } + delete[] ext; +#else // !FAST_COMPILE + if( !PointWeight.set ) PointWeight.value = DefaultPointWeightMultiplier*Degree.value; + char *ext = GetFileExtension( In.value ); + if( !strcasecmp( ext , "ply" ) ) + { + typedef VertexFactory::Factory< Real , typename VertexFactory::PositionFactory< Real , DEFAULT_DIMENSION > , typename VertexFactory::NormalFactory< Real , DEFAULT_DIMENSION > > Factory; + Factory factory; + bool *readFlags = new bool[ factory.plyReadNum() ]; + std::vector< PlyProperty > unprocessedProperties; + PLY::ReadVertexHeader( In.value , factory , readFlags , unprocessedProperties ); + if( !factory.plyValidReadProperties<0>( readFlags ) ) ERROR_OUT( "Ply file does not contain positions" ); + if( !factory.plyValidReadProperties<1>( readFlags ) ) ERROR_OUT( "Ply file does not contain normals" ); + delete[] readFlags; + + if( unprocessedProperties.size() ) Execute< DEFAULT_DIMENSION , Real >( VertexFactory::DynamicFactory< Real >( unprocessedProperties ) ); + else Execute< DEFAULT_DIMENSION , Real >( VertexFactory::EmptyFactory< Real >() ); + } + else + { + if( Colors.set ) Execute< DEFAULT_DIMENSION , Real >( VertexFactory::RGBColorFactory< Real >() ); + else Execute< DEFAULT_DIMENSION , Real >( VertexFactory::EmptyFactory< Real >() ); + } + delete[] ext; +#endif // FAST_COMPILE + if( Performance.set ) + { + printf( "Time (Wall/CPU): %.2f / %.2f\n" , timer.wallTime() , timer.cpuTime() ); + printf( "Peak Memory (MB): %d\n" , MemoryInfo::PeakMemoryUsageMB() ); + } + + ThreadPool::Terminate(); + return EXIT_SUCCESS; +} diff --git a/Src/PoissonReconLib.h b/Src/PoissonReconLib.h new file mode 100644 index 00000000..30c563c2 --- /dev/null +++ b/Src/PoissonReconLib.h @@ -0,0 +1,31 @@ +/* +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. +*/ + +#pragma once + +int PoissonReconLib(int argc, char **argv); From f4e41eee33fefd2adcc00b443215b55dfa2ce520 Mon Sep 17 00:00:00 2001 From: Harsh Bhatia Date: Wed, 21 Jul 2021 18:25:14 -0700 Subject: [PATCH 3/5] add cxx flags and use existing libs if possible --- CMakeLists.txt | 87 ++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 63 insertions(+), 24 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index fa694ae4..ffbcace7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,44 +2,81 @@ cmake_minimum_required(VERSION 3.20) project(AdaptiveMultigridSolvers) - set(CMAKE_VERBOSE_MAKEFILE ON) ## ----------------------------------------------------------------------------- set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -std=c++14 -pthread -Wno-invalid-offsetof -fPIC -Wno-deprecated") -set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${CMAKE_CXX_FLAGS} -g") -set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${CMAKE_CXX_FLAGS} -Ofast -DRELEASE -funroll-loops -ffast-math") +# set custom compiler flags +if(CMAKE_COMPILER_IS_GNUCXX) + set(CFLAGS "-std=c++14 -pthread -Wno-deprecated -Wno-invalid-offsetof") + set(LFLAGS "-lstdc++ -lpthread") +else() + set(CFLAGS, "-std=c++14 -pthread -Wno-deprecated -Wno-invalid-offsetof -Wno-dangling-else") + set(LFLAGS, "-lstdc++") +endif() -## ----------------------------------------------------------------------------- -find_package(OpenMP) +set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${CFLAGS} -g3") +set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${CFLAGS} -Ofast -DRELEASE -funroll-loops -ffast-math") +set(CMAKE_EXE_LINKER_FLAGS_DEBUG "${CMAKE_EXE_LINKER_FLAGS_DEBUG} ${LFLAGS}") +set(CMAKE_EXE_LINKER_FLAGS_RELEASE "${CMAKE_EXE_LINKER_FLAGS_RELEASE} ${LFLAGS} -Ofast") +## ----------------------------------------------------------------------------- +find_package(OpenMP REQUIRED) ## ----------------------------------------------------------------------------- -# all include paths -include_directories(${PROJECT_SOURCE_DIR}) -include_directories(${PROJECT_SOURCE_DIR}/Src) -include_directories(${PROJECT_SOURCE_DIR}/PNG) -include_directories(${PROJECT_SOURCE_DIR}/JPEG) -include_directories(${PROJECT_SOURCE_DIR}/ZLIB) +find_package(PNG) # PNG depends on zlib and adds zlib to path already +if(PNG_FOUND) + #message("-- png lib found: include_dir=(${PNG_INCLUDE_DIRS}), libs=(${PNG_LIBRARIES})") +else() + + set(PNG_PATH "${PROJECT_SOURCE_DIR}/PNG") + message("-- png lib not found: will build from (${PNG_PATH})") + + set(PNG_INCLUDE_DIRS ${PNG_PATH}) + set(PNG_LIBRARIES "png") # will result in libpng + + file(GLOB_RECURSE SRCS_PNG "${PNG_PATH}/*.c") + add_library(${PNG_LIBRARIES} SHARED ${SRCS_PNG}) +endif() ## ----------------------------------------------------------------------------- -## build dependencies as libraries +find_package(ZLIB) # dont really need to check for zlib if png was success +if(ZLIB_FOUND) + #message("-- zlib found: include_dir=(${ZLIB_INCLUDE_DIRS}), libs=(${ZLIB_LIBRARIES})") +else() + set(ZLIB_PATH "${PROJECT_SOURCE_DIR}/ZLIB") + message("-- zlib not found: will build from (${ZLIB_PATH})") + + set(ZLIB_INCLUDE_DIRS ${ZLIB_PATH}) + set(ZLIB_LIBRARIES "z") # will result in libz + + file(GLOB_RECURSE SRCS_ZLIB "${ZLIB_PATH}/*.c") + add_library(${ZLIB_LIBRARIES} SHARED ${SRCS_ZLIB}) +endif() + ## ----------------------------------------------------------------------------- -set(TRG_PNG pr_png) -set(TRG_JPEG pr_jpeg) -set(TRG_ZLIB pr_zlib) +find_package(JPEG) +if(JPEG_FOUND) + #message("-- jpeg lib found: include_dir=(${JPEG_INCLUDE_DIRS}), libs=(${JPEG_LIBRARIES})") +else() + set(JPEG_PATH "${PROJECT_SOURCE_DIR}/JPEG") + message("-- jpeg lib not found: will build from (${JPEG_PATH})") -file(GLOB_RECURSE SRCS_PNG "PNG/*.c") -file(GLOB_RECURSE SRCS_JPEG "JPEG/*.cpp") -file(GLOB_RECURSE SRCS_ZLIB "ZLIB/*.c") + set(JPEG_INCLUDE_DIRS ${JPEG}) + set(JPEG_LIBRARIES "jpeg") # will result in libjpeg -add_library(${TRG_PNG} STATIC ${SRCS_PNG}) -add_library(${TRG_JPEG} SHARED ${SRCS_JPEG}) -add_library(${TRG_ZLIB} SHARED ${SRCS_ZLIB}) + file(GLOB_RECURSE SRCS_JPEG "${JPEG_PATH}/*.c") + add_library(${JPEG_LIBRARIES} SHARED ${SRCS_JPEG}) +endif() + + +## ----------------------------------------------------------------------------- +# all include paths +include_directories(${PROJECT_SOURCE_DIR} ${PROJECT_SOURCE_DIR}/Src) +include_directories(${PNG_INCLUDE_DIRS} ${ZLIB_INCLUDE_DIRS} ${JPEG_INCLUDE_DIRS}) ## ----------------------------------------------------------------------------- @@ -50,12 +87,14 @@ add_library(poisson_recon_lib SHARED ${PROJECT_SOURCE_DIR}/Src/PoissonReconLib.cpp) target_link_libraries(poisson_recon_lib - ${TRG_PNG} ${TRG_JPEG} ${TRG_ZLIB}) + ${PNG_LIBRARIES} ${ZLIB_LIBRARIES} ${JPEG_LIBRARIES}) + add_executable(poisson_recon ${PROJECT_SOURCE_DIR}/Src/PoissonRecon.cpp) target_link_libraries(poisson_recon - poisson_recon_lib ${TRG_PNG} ${TRG_JPEG} ${TRG_ZLIB}) + poisson_recon_lib + ${PNG_LIBRARIES} ${ZLIB_LIBRARIES} ${JPEG_LIBRARIES}) ## ----------------------------------------------------------------------------- From 407ffcacaf0af79550307046a2a8aba35d00b86c Mon Sep 17 00:00:00 2001 From: Harsh Bhatia Date: Wed, 21 Jul 2021 19:50:31 -0700 Subject: [PATCH 4/5] Added python module update gitignore for python temp files --- .gitignore | 10 +- Python/ReadMe.md | 13 +++ Python/pypoissonrecon.pyx | 193 ++++++++++++++++++++++++++++++++++++++ Python/setup.py | 68 ++++++++++++++ Src/PoissonReconLib.cpp | 9 ++ Src/PoissonReconLib.h | 6 ++ 6 files changed, 298 insertions(+), 1 deletion(-) create mode 100644 Python/ReadMe.md create mode 100644 Python/pypoissonrecon.pyx create mode 100644 Python/setup.py diff --git a/.gitignore b/.gitignore index 73ae1738..62c8871f 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,12 @@ *.opensdf *.bat *.zip -.vs/* \ No newline at end of file +.vs/* + + +*.DS_Store +CMakeLists.txt.user +Python/build/ +Python/dist/ +Python/pypoissonrecon.cpp +Python/pypoissonrecon.egg-info/ diff --git a/Python/ReadMe.md b/Python/ReadMe.md new file mode 100644 index 00000000..dd0b0f33 --- /dev/null +++ b/Python/ReadMe.md @@ -0,0 +1,13 @@ +### Python Bindings for Poisson Reconstruction Library + +##### Installation +``` +$ python setup.py build_ext +$ python setup.py install --prefix= +``` + +##### Contributed by + +* The original [Matlab bindings](https://github.com/daeyun/poisson-surface-reconstruction) were developed by [Daeyun Shin](https://github.com/daeyun). +* Adapted to Python by [Miguel Molero](https://github.com/mmolero/pypoisson) as a separate project [pypoisson](https://github.com/mmolero/pypoisson), which ultimately got outdated. +* Modernized to version 13.72 and contributed directly to the main repository by [Harsh Bhatia](https://github.com/bhatiaharsh) for improved versioning and management. diff --git a/Python/pypoissonrecon.pyx b/Python/pypoissonrecon.pyx new file mode 100644 index 00000000..93efbc5d --- /dev/null +++ b/Python/pypoissonrecon.pyx @@ -0,0 +1,193 @@ +#cython: language_level=3str + +cimport numpy as np +import numpy as np +from libcpp.vector cimport vector +from libcpp.string cimport string +from libc.stdlib cimport malloc, free +from libcpp cimport bool + +np.import_array() + +cdef extern from "../Src/PoissonReconLib.h": + cdef int PoissonReconLib(int argc, char* argv[]) + cdef vector[double] double_data + cdef vector[int] int_data + cdef vector[double] mem_data + +def poisson_reconstruction(points, normals, + depth=8, full_depth=5, scale=1.1, + samples_per_node=1.0, cg_depth=0.0, + enable_polygon_mesh=False, enable_density=False, + nthreads=0, verbose=False): + + """ + Python Binding of Poisson Surface Reconstruction + Usage: + + faces, vertices = poisson_reconstruction(points, normals, depth=10) + + + Parameters + ---------- + + points: array-like + + list of oriented vertices with the x-, y-, and z-coordinates of the positions + + normals: array-like + + list of the x-, y-, and z-coordinates of the normals + + depth: Integer + + This integer is the maximum depth of the tree that will be used for surface reconstruction. + Running at depth d corresponds to solving on a voxel grid whose resolution is no larger than 2^d x 2^d x 2^d. + Note that since the reconstructor adapts the octree to the sampling density, the specified reconstruction depth is only an upper bound. + The default value for this parameter is 8. + + full_depth: Integer + + This integer specifies the depth beyond depth the octree will be adapted. + At coarser depths, the octree will be complete, containing all 2^d x 2^d x 2^d nodes. + The default value for this parameter is 5. + + scale: float + + This floating point value specifies the ratio between the diameter of the cube used for reconstruction and the diameter of the samples' bounding cube. + The default value is 1.1. + + samples_per_node: float + + This floating point value specifies the minimum number of sample points that should fall within an octree node as the octree construction is adapted to sampling density. + For noise-free samples, small values in the range [1.0 - 5.0] can be used. + For more noisy samples, larger values in the range [15.0 - 20.0] may be needed to provide a smoother, noise-reduced, reconstruction. + The default value is 1.0. + + cg_depth: Integer + + This integer is the depth up to which a conjugate-gradients solver will be used to solve the linear system. + Beyond this depth Gauss-Seidel relaxation will be used. + The default value for this parameter is 0. + + enable_polygon_mesh: Bool + Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes). + The default value for this parameter is False. + + enable_density: Bool + Enabling this flag tells the reconstructor to output the estimated depth values of the iso-surface vertices + The default value for this parameter is False. + + nthreads: int + Number of omp threads to use. Default = 0 = all available threads. + + verbose: Bool + Enable verbose mode. + + + + Returns + ------- + + faces: array-like + faces of the reconstructed mesh + + vertices: array-like + + vertices of the reconstructed mesh + + + + + """ + + return _poisson_reconstruction(np.ascontiguousarray(np.float64(points)), + np.ascontiguousarray(np.float64(normals)), + depth, full_depth, scale, samples_per_node, + cg_depth, enable_polygon_mesh, enable_density, + nthreads, verbose) + + + +cdef _poisson_reconstruction(np.float64_t[:, ::1] points, + np.float64_t[:, ::1] normals, + int depth=8, + int full_depth=5, + double scale=1.10, + double samples_per_node=1.0, + double cg_depth=0.0, + bool enable_polygon_mesh=False, + bool enable_density=False, + int nthreads=0, + bool verbose=False): + + cdef: + char **c_argv + string arg_depth = str(depth).encode() + string arg_full_depth = str(full_depth).encode() + string arg_scale = str(scale).encode() + string arg_samples_per_node = str(samples_per_node).encode() + string arg_cg_depth = str(cg_depth).encode() + string arg_nthreads = str(nthreads).encode() + + int_data.clear() + double_data.clear() + mem_data.clear() + + point_nrows, point_ncols = np.shape(points) + normal_nrows, normal_ncols = np.shape(normals) + + mem_data.resize(point_ncols * point_nrows + normal_ncols * normal_nrows) + + for i in range(point_nrows): + for j in range(point_ncols): + mem_data[j + i*(point_ncols + normal_ncols)] = points[i,j] + mem_data[j + point_ncols + i *(point_ncols + normal_ncols)] = normals[i,j] + + + args = [b"PoissonRecon", b"--in", b"none", b"--out", b"none", b"--depth", arg_depth.c_str(), + b"--fullDepth", arg_full_depth.c_str(), b"--scale", arg_scale.c_str(), + b"--samplesPerNode", arg_samples_per_node.c_str(), + b"--cgDepth", arg_cg_depth.c_str()] + + if verbose == True: + args += [b"--verbose"] + + if nthreads > 0: + args += [b"--threads", arg_nthreads.c_str()] + + if enable_polygon_mesh: + args += [b"--polygonMesh"] + if enable_density: + args += [b"--density"] + + c_argv = malloc(sizeof(char*) * len(args)) + for idx, s in enumerate(args): + c_argv[idx] = s + + try: + PoissonReconLib(len(args), c_argv) + finally: + free(c_argv) + + + face_cols, vertex_cols = 3, 3 + face_rows = int_data.size() // face_cols + vertex_rows = double_data.size() // vertex_cols + + cdef int *ptr_faces = &int_data[0] + cdef double *ptr_vertices = &double_data[0] + + faces = np.zeros((face_rows*face_cols,), dtype=np.int32 ) + vertices = np.zeros((vertex_rows*vertex_cols,), dtype=np.float64) + + for i in range(face_rows*face_cols): + faces[i] = ptr_faces[i] + + for i in range(vertex_rows*vertex_cols): + vertices[i] = ptr_vertices[i] + + int_data.clear() + double_data.clear() + + return faces.reshape(face_rows,face_cols), vertices.reshape(vertex_rows,vertex_cols) diff --git a/Python/setup.py b/Python/setup.py new file mode 100644 index 00000000..e2c3e7a5 --- /dev/null +++ b/Python/setup.py @@ -0,0 +1,68 @@ +import os +import numpy + +from setuptools import setup, Extension +from setuptools.command.build_ext import build_ext + +# ------------------------------------------------------------------------------ +# paths to these external libraries +libs = ['png', 'z', 'jpeg'] + +# TODO: need to properly pass in the paths to these external libs +ext_inc_dir = ['/opt/local/include'] +ext_lib_dir = ['/opt/local/lib'] + +# ------------------------------------------------------------------------------ +# path to the project directory +proj_dir = os.path.dirname(os.path.realpath('.')) + +# code to be build for this extension module +sources = [os.path.join(proj_dir, 'Src', 'PoissonReconLib.cpp'), + 'pypoissonrecon.pyx'] + +include_dirs = ext_inc_dir + [proj_dir] +lib_dirs = ext_lib_dir + +# ------------------------------------------------------------------------------ +# flags taken from original Makefile +CXX_FLAGS = "-fopenmp -std=c++14 -pthread -fPIC -Wno-deprecated -Wno-invalid-offsetof " +CXX_LINKER_FLAGS = "-fopenmp -lstdc++ -lpthread " + + # let's just assume we are building release +if True: + CXX_FLAGS += "-Ofast -DRELEASE -funroll-loops -ffast-math " + CXX_LINKER_FLAGS += "-Ofast " + +# if we wanted to build debug +else: + CXX_FLAGS += "-g3 " + CXX_LINKER_FLAGS += "" + +# flags to suppress further warnings +if True: + CXX_FLAGS += "-Wno-unused-variable -Wno-unused-but-set-variable -Wno-unused-local-typedefs " + CXX_FLAGS += "-Wno-delete-non-virtual-dtor -Wno-class-memaccess -Wno-strict-aliasing " + CXX_FLAGS += "-Wno-sign-compare -Wno-reorder -Wno-dangling-else " + CXX_FLAGS += "-Wno-maybe-uninitialized -Wno-format-overflow " + +# ------------------------------------------------------------------------------ +extn_mod = Extension('pypoissonrecon', language='c++', + sources = sources, + include_dirs = include_dirs + [numpy.get_include()], + extra_compile_args = CXX_FLAGS.split(), + extra_link_args = CXX_LINKER_FLAGS.split(), + libraries=libs, + library_dirs=lib_dirs) + +# ------------------------------------------------------------------------------ +setup( + name='pypoissonrecon', + version='13.72', + description='Poisson Surface Reconstruction (Adaptive Multigrid Solvers)', + author='', + author_email='', + url='https://github.com/mkazhdan/PoissonRecon', + cmdclass = {'build_ext': build_ext}, + ext_modules = [extn_mod] +) +# ------------------------------------------------------------------------------ diff --git a/Src/PoissonReconLib.cpp b/Src/PoissonReconLib.cpp index 8bacda1e..ecf87456 100644 --- a/Src/PoissonReconLib.cpp +++ b/Src/PoissonReconLib.cpp @@ -26,6 +26,15 @@ ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF S DAMAGE. */ +#include +std::vector double_data; +std::vector int_data; +std::vector mem_data; +bool* is_verbose; +bool is_writing_to_memory = true; +bool is_reading_from_memory = true; +bool is_writing_int; + #include "PreProcessor.h" #undef USE_DOUBLE // If enabled, double-precesion is used diff --git a/Src/PoissonReconLib.h b/Src/PoissonReconLib.h index 30c563c2..4ba1d697 100644 --- a/Src/PoissonReconLib.h +++ b/Src/PoissonReconLib.h @@ -28,4 +28,10 @@ DAMAGE. #pragma once +#include +extern std::vector double_data; +extern std::vector int_data; +extern std::vector mem_data; +extern bool* is_verbose; + int PoissonReconLib(int argc, char **argv); From 6a928050e41a428cbf5fa087783920489eed9654 Mon Sep 17 00:00:00 2001 From: Harsh Bhatia Date: Wed, 21 Jul 2021 22:22:07 -0700 Subject: [PATCH 5/5] removed cg_depth parameter from pyx interface --- Python/pypoissonrecon.pyx | 35 +++++++++++++++-------------------- 1 file changed, 15 insertions(+), 20 deletions(-) diff --git a/Python/pypoissonrecon.pyx b/Python/pypoissonrecon.pyx index 93efbc5d..09c38161 100644 --- a/Python/pypoissonrecon.pyx +++ b/Python/pypoissonrecon.pyx @@ -1,5 +1,8 @@ #cython: language_level=3str +## 2021.07.21: code adapted from https://github.com/mmolero/pypoisson +## https://github.com/mmolero/pypoisson/blob/4f6a651d20ff2bab5ec0d0f8d0a14cc6cef3ef51/src/pypoisson.pyx + cimport numpy as np import numpy as np from libcpp.vector cimport vector @@ -17,9 +20,9 @@ cdef extern from "../Src/PoissonReconLib.h": def poisson_reconstruction(points, normals, depth=8, full_depth=5, scale=1.1, - samples_per_node=1.0, cg_depth=0.0, + samples_per_node=1.0, enable_polygon_mesh=False, enable_density=False, - nthreads=0, verbose=False): + nthreads=8, verbose=False): """ Python Binding of Poisson Surface Reconstruction @@ -64,12 +67,6 @@ def poisson_reconstruction(points, normals, For more noisy samples, larger values in the range [15.0 - 20.0] may be needed to provide a smoother, noise-reduced, reconstruction. The default value is 1.0. - cg_depth: Integer - - This integer is the depth up to which a conjugate-gradients solver will be used to solve the linear system. - Beyond this depth Gauss-Seidel relaxation will be used. - The default value for this parameter is 0. - enable_polygon_mesh: Bool Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes). The default value for this parameter is False. @@ -79,7 +76,8 @@ def poisson_reconstruction(points, normals, The default value for this parameter is False. nthreads: int - Number of omp threads to use. Default = 0 = all available threads. + This parameter specifies the number of threads across which the solver should be parallelized. + The default value for this parameter is 8. verbose: Bool Enable verbose mode. @@ -104,7 +102,7 @@ def poisson_reconstruction(points, normals, return _poisson_reconstruction(np.ascontiguousarray(np.float64(points)), np.ascontiguousarray(np.float64(normals)), depth, full_depth, scale, samples_per_node, - cg_depth, enable_polygon_mesh, enable_density, + enable_polygon_mesh, enable_density, nthreads, verbose) @@ -115,7 +113,6 @@ cdef _poisson_reconstruction(np.float64_t[:, ::1] points, int full_depth=5, double scale=1.10, double samples_per_node=1.0, - double cg_depth=0.0, bool enable_polygon_mesh=False, bool enable_density=False, int nthreads=0, @@ -127,7 +124,6 @@ cdef _poisson_reconstruction(np.float64_t[:, ::1] points, string arg_full_depth = str(full_depth).encode() string arg_scale = str(scale).encode() string arg_samples_per_node = str(samples_per_node).encode() - string arg_cg_depth = str(cg_depth).encode() string arg_nthreads = str(nthreads).encode() int_data.clear() @@ -145,17 +141,16 @@ cdef _poisson_reconstruction(np.float64_t[:, ::1] points, mem_data[j + point_ncols + i *(point_ncols + normal_ncols)] = normals[i,j] - args = [b"PoissonRecon", b"--in", b"none", b"--out", b"none", b"--depth", arg_depth.c_str(), - b"--fullDepth", arg_full_depth.c_str(), b"--scale", arg_scale.c_str(), - b"--samplesPerNode", arg_samples_per_node.c_str(), - b"--cgDepth", arg_cg_depth.c_str()] + args = [b"PoissonRecon", b"--in", b"none", b"--out", b"none", + b"--depth", arg_depth.c_str(), + b"--fullDepth", arg_full_depth.c_str(), + b"--scale", arg_scale.c_str(), + b"--samplesPerNode", arg_samples_per_node.c_str(), + b"--threads", arg_nthreads.c_str() + ] if verbose == True: args += [b"--verbose"] - - if nthreads > 0: - args += [b"--threads", arg_nthreads.c_str()] - if enable_polygon_mesh: args += [b"--polygonMesh"] if enable_density: