diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..26e1ac5 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,34 @@ +cmake_minimum_required(VERSION 2.8) + +include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ) + +project( QPOISSON_RECON_LIB ) + +#to remove very annoying "NO ARGUMENTS THAT DEPEND ON A TEMPLATE PARAMETER" errors with G++! +if( UNIX OR MINGW ) + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fpermissive -w") +endif() + +set( SAMPLE_PROJECT_VERSION_MAJOR 6 ) +set( SAMPLE_PROJECT_VERSION_MINOR 13a ) + +file( GLOB header_list Src_CC_wrap/*.h ) +list (APPEND header_list Src/Allocator.h Src/Array.h Src/BinaryNode.h Src/BSplineData.h Src/CmdLineParser.h Src/Factor.h Src/FunctionData.h Src/Geometry.h Src/Hash.h Src/MarchingCubes.h Src/MAT.h Src/MemoryUsage.h Src/MultiGridOctreeData.h Src/Octree.h Src/Ply.h Src/PointStream.h Src/Polynomial.h Src/PPolynomial.h Src/SparseMatrix.h Src/Vector.h ) +list (APPEND header_list Src/Array.inl Src/BSplineData.inl Src/CmdLineParser.inl Src/FunctionData.inl Src/Geometry.inl Src/MAT.inl Src/MultiGridOctreeData.inl Src/MultiGridOctreeData.IsoSurface.inl Src/MultiGridOctreeData.SortedTreeNodes.inl Src/Octree.inl Src/PointStream.inl Src/Polynomial.inl Src/PPolynomial.inl Src/SparseMatrix.inl Src/Vector.inl ) +file( GLOB source_list Src_CC_wrap/*.cpp ) +list (APPEND source_list Src/CmdLineParser.cpp Src/Factor.cpp Src/Geometry.cpp Src/MarchingCubes.cpp Src/PlyFile.cpp ) + +add_library( ${PROJECT_NAME} STATIC ${header_list} ${inline_list} ${source_list} ) + +# Add prepocessor definitions +set_property( TARGET ${PROJECT_NAME} APPEND PROPERTY COMPILE_DEFINITIONS _CRT_SECURE_NO_DEPRECATE _CRT_SECURE_NO_WARNINGS NOMINMAX ) + +option( POISSON_RECON_WITH_OPEN_MP "Check to compile PoissonRecon plugin with OpenMP support" OFF ) +if ( POISSON_RECON_WITH_OPEN_MP ) + find_package(OpenMP) + if (OPENMP_FOUND) + set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + set_property( TARGET ${PROJECT_NAME} APPEND PROPERTY COMPILE_DEFINITIONS WITH_OPENMP ) + endif() +endif () diff --git a/Src/Geometry.h b/Src/Geometry.h index 0d73eda..beb361a 100644 --- a/Src/Geometry.h +++ b/Src/Geometry.h @@ -322,6 +322,7 @@ class CoredVectorMeshData : public CoredMeshData< Vertex > std::vector< std::vector< int > > polygons; int polygonIndex; int oocPointIndex; + void* writelock; public: CoredVectorMeshData(void); diff --git a/Src/Geometry.inl b/Src/Geometry.inl index 97be936..6b27cf2 100644 --- a/Src/Geometry.inl +++ b/Src/Geometry.inl @@ -27,6 +27,7 @@ DAMAGE. */ #include +#include template Real Random(void){return Real(rand())/RAND_MAX;} @@ -426,7 +427,11 @@ int Triangulation::flipMinimize(int eIndex){ // CoredVectorMeshData // ///////////////////////// template< class Vertex > -CoredVectorMeshData< Vertex >::CoredVectorMeshData( void ) { oocPointIndex = polygonIndex = 0; } +CoredVectorMeshData< Vertex >::CoredVectorMeshData( void ) : writelock(0) +{ + oocPointIndex = polygonIndex = 0; + omp_init_lock(&writelock); +} template< class Vertex > void CoredVectorMeshData< Vertex >::resetIterator ( void ) { oocPointIndex = polygonIndex = 0; } template< class Vertex > @@ -450,11 +455,13 @@ template< class Vertex > int CoredVectorMeshData< Vertex >::addPolygon_s( const std::vector< int >& polygon ) { size_t sz; -#pragma omp critical (CoredVectorMeshData_addPolygon_s) +//#pragma omp critical (CoredVectorMeshData_addPolygon_s) + omp_set_lock(reinterpret_cast(&writelock)); { sz = polygon.size(); polygons.push_back( polygon ); } + omp_unset_lock(reinterpret_cast(&writelock)); return (int)sz; } template< class Vertex > diff --git a/Src_CC_wrap/PoissonReconLib.cpp b/Src_CC_wrap/PoissonReconLib.cpp new file mode 100644 index 0000000..72104f5 --- /dev/null +++ b/Src_CC_wrap/PoissonReconLib.cpp @@ -0,0 +1,311 @@ +//########################################################################## +//# # +//# CLOUDCOMPARE WRAPPER: PoissonReconLib # +//# # +//# This program is free software; you can redistribute it and/or modify # +//# it under the terms of the GNU General Public License as published by # +//# the Free Software Foundation; version 2 of the License. # +//# # +//# This program is distributed in the hope that it will be useful, # +//# but WITHOUT ANY WARRANTY; without even the implied warranty of # +//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +//# GNU General Public License for more details. # +//# # +//# COPYRIGHT: Daniel Girardeau-Montaut # +//# # +//########################################################################## + +#include "PoissonReconLib.h" + +#ifdef _WIN32 +#include +#include +#endif // _WIN32 + +#ifdef WITH_OPENMP +#include +#endif + +//PoissonRecon +#include "../Src/Ply.h" +#include "../Src/Array.h" +#include "../Src/Octree.h" +#include "../Src/SparseMatrix.h" + +#define DumpOutput(...) ((void)0) +#include "../Src/MultiGridOctreeData.h" //only after DumpOutput has been defined! + +PoissonReconLib::Parameters::Parameters() + : depth(8) //8 + , cgDepth(0) //0 + , kernelDepth(0) //? + , adaptiveExp(1) //AdaptiveExponent (1) + , iters(8) //8 + , fullDepth(5) //5 + , minDepth(0) //0 + , maxSolveDepth(0) //? + , boundary(1) //1 + , threads(1) //ideally omp_get_num_procs() + , samplesPerNode(1.0f) //1.0f + , scale(1.1f) //1.1f + , cgAccuracy(1.0e-3f) //1.0e-3f + , pointWeight(4.0f) //4.0f + , complete(false) + , showResidual(false) + , confidence(false) + , normalWeights(false) + , nonManifold(false) + , density(false) + , colorInterp(16.0f) +{ +#ifdef WITH_OPENMP + threads = omp_get_num_procs(); +#endif +} + +template< class PointCoordinateType, class Real, class Vertex > +bool Execute(PoissonReconLib::Parameters params, OrientedPointStream< PointCoordinateType >* pointStream, CoredVectorMeshData< Vertex >& mesh) +{ + XForm4x4< Real > xForm = XForm4x4< Real >::Identity(); + XForm4x4< Real > iXForm = xForm.inverse(); + + //DGM: reset static parameters!!! + TreeNodeData::NodeCount = 0; + + Octree< Real > tree; + tree.threads = params.threads; + + if (params.maxSolveDepth == 0) + params.maxSolveDepth = params.depth; + + OctNode< TreeNodeData >::SetAllocator( MEMORY_ALLOCATOR_BLOCK_SIZE ); + + if (params.maxSolveDepth < 2) + return false; + int kernelDepth = params.kernelDepth != 0 ? params.kernelDepth : params.maxSolveDepth-2; + if( kernelDepth > params.depth ) + return false; + params.fullDepth = std::min(params.fullDepth, params.depth); + + tree.maxMemoryUsage = 0; + typename Octree< Real >::template SparseNodeData< typename Octree< Real >::PointData >* pointInfo = new typename Octree< Real >::template SparseNodeData< typename Octree< Real >::PointData >(); + typename Octree< Real >::template SparseNodeData< Point3D< Real > >* normalInfo = new typename Octree< Real >::template SparseNodeData< Point3D< Real > >(); + std::vector< Real >* kernelDensityWeights = new std::vector< Real >(); + std::vector< Real >* centerWeights = new std::vector< Real >(); + typedef typename Octree< Real >::template ProjectiveData< Point3D< Real > > ProjectiveColor; + + int pointCount = tree.template SetTree< PointCoordinateType >( + pointStream, + params.minDepth, + params.depth, + params.fullDepth, + kernelDepth, + static_cast(params.samplesPerNode), + params.scale, + params.confidence, + params.normalWeights, + params.pointWeight, + params.adaptiveExp, + *kernelDensityWeights, + *pointInfo, + *normalInfo, + *centerWeights, + xForm, + params.boundary, + params.complete ); + + if( !params.density ) + { + delete kernelDensityWeights; + kernelDensityWeights = NULL; + } + + //DumpOutput( "Input Points: %d\n" , pointCount ); + //DumpOutput( "Leaves/Nodes: %d/%d\n" , tree.tree.leaves() , tree.tree.nodes() ); + + double maxMemoryUsage = tree.maxMemoryUsage; + tree.maxMemoryUsage = 0; + + Pointer( Real ) constraints = tree.SetLaplacianConstraints( *normalInfo ); + delete normalInfo; + normalInfo = 0; + + maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage ); + tree.maxMemoryUsage = 0; + + Pointer( Real ) solution = tree.SolveSystem( *pointInfo , constraints , params.showResidual , params.iters, params.maxSolveDepth, params.cgDepth, params.cgAccuracy ); + + delete pointInfo; + pointInfo = 0; + FreePointer( constraints ); + + maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage ); + + Real isoValue = tree.GetIsoValue( solution , *centerWeights ); + delete centerWeights; + centerWeights = 0; + + //DumpOutput( "Iso-Value: %e\n" , isoValue ); + + //output + tree.maxMemoryUsage = 0; + + tree.GetMCIsoSurface( kernelDensityWeights ? GetPointer( *kernelDensityWeights ) : NullPointer( Real ), + NULL, + solution, + isoValue, + mesh, + true, + !params.nonManifold, + false ); + + maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage ); + + //DumpOutput( "Vertices / Polygons: %d / %d\n" , mesh.outOfCorePointCount()+mesh.inCorePoints.size() , mesh.polygonCount() ); + + FreePointer( solution ); + + return true; +} + +template< class PointCoordinateType, class Real, class Vertex > +bool Execute(PoissonReconLib::Parameters params, OrientedPointStreamWithData< PointCoordinateType , Point3D< unsigned char > >* pointStream, CoredVectorMeshData< Vertex >& mesh) +{ + XForm4x4< Real > xForm = XForm4x4< Real >::Identity(); + XForm4x4< Real > iXForm = xForm.inverse(); + + //DGM: reset static parameters!!! + TreeNodeData::NodeCount = 0; + + Octree< Real > tree; + tree.threads = params.threads; + + if (params.maxSolveDepth == 0) + params.maxSolveDepth = params.depth; + + OctNode< TreeNodeData >::SetAllocator( MEMORY_ALLOCATOR_BLOCK_SIZE ); + + if (params.maxSolveDepth < 2) + return false; + int kernelDepth = params.kernelDepth != 0 ? params.kernelDepth : params.maxSolveDepth-2; + if( kernelDepth > params.depth ) + return false; + params.fullDepth = std::min(params.fullDepth, params.depth); + + tree.maxMemoryUsage = 0; + typename Octree< Real >::template SparseNodeData< typename Octree< Real >::PointData >* pointInfo = new typename Octree< Real >::template SparseNodeData< typename Octree< Real >::PointData >(); + typename Octree< Real >::template SparseNodeData< Point3D< Real > >* normalInfo = new typename Octree< Real >::template SparseNodeData< Point3D< Real > >(); + std::vector< Real >* kernelDensityWeights = new std::vector< Real >(); + std::vector< Real >* centerWeights = new std::vector< Real >(); + typedef typename Octree< Real >::template ProjectiveData< Point3D< Real > > ProjectiveColor; + typename Octree< Real >::template SparseNodeData< ProjectiveColor > colorData; + + int pointCount = tree.template SetTree< PointCoordinateType >( + pointStream, + params.minDepth, + params.depth, + params.fullDepth, + kernelDepth, + static_cast(params.samplesPerNode), + params.scale, + params.confidence, + params.normalWeights, + params.pointWeight, + params.adaptiveExp, + *kernelDensityWeights, + *pointInfo, + *normalInfo, + *centerWeights, + colorData, + xForm, + params.boundary, + params.complete ); + + for (const OctNode< TreeNodeData >* n = tree.tree.nextNode(); n != NULL; n = tree.tree.nextNode(n)) + { + int idx = colorData.index(n); + if (idx >= 0) + colorData.data[idx] *= static_cast(pow(params.colorInterp, n->depth())); + } + + if( !params.density ) + { + delete kernelDensityWeights; + kernelDensityWeights = NULL; + } + + //DumpOutput( "Input Points: %d\n" , pointCount ); + //DumpOutput( "Leaves/Nodes: %d/%d\n" , tree.tree.leaves() , tree.tree.nodes() ); + + double maxMemoryUsage = tree.maxMemoryUsage; + tree.maxMemoryUsage = 0; + + Pointer( Real ) constraints = tree.SetLaplacianConstraints( *normalInfo ); + delete normalInfo; + normalInfo = 0; + + maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage ); + tree.maxMemoryUsage = 0; + + Pointer( Real ) solution = tree.SolveSystem( *pointInfo , constraints , params.showResidual , params.iters, params.maxSolveDepth, params.cgDepth, params.cgAccuracy ); + + delete pointInfo; + pointInfo = 0; + FreePointer( constraints ); + + maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage ); + + Real isoValue = tree.GetIsoValue( solution , *centerWeights ); + delete centerWeights; + centerWeights = 0; + + //DumpOutput( "Iso-Value: %e\n" , isoValue ); + + //output + tree.maxMemoryUsage = 0; + + tree.GetMCIsoSurface( kernelDensityWeights ? GetPointer( *kernelDensityWeights ) : NullPointer( Real ), + &colorData, + solution, + isoValue, + mesh, + true, + !params.nonManifold, + false ); + + maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage ); + + //DumpOutput( "Vertices / Polygons: %d / %d\n" , mesh.outOfCorePointCount()+mesh.inCorePoints.size() , mesh.polygonCount() ); + + FreePointer( solution ); + + return true; +} + +bool PoissonReconLib::Reconstruct(Parameters params, OrientedPointStreamWithData< float , Point3D< unsigned char > >* pointStream, CoredVectorMeshData< PlyColorAndValueVertex< float > >& mesh) +{ + return Execute< float, + float, + PlyColorAndValueVertex< float > > (params, pointStream, mesh); +} + +bool PoissonReconLib::Reconstruct(Parameters params, OrientedPointStream< float >* pointStream, CoredVectorMeshData< PlyValueVertex< float > >& mesh) +{ + return Execute< float, + float, + PlyValueVertex< float > > (params, pointStream, mesh); +} + +bool PoissonReconLib::Reconstruct(Parameters params, OrientedPointStreamWithData< double , Point3D< unsigned char > >* pointStream, CoredVectorMeshData< PlyColorAndValueVertex< double > >& mesh) +{ + return Execute< double, + double, + PlyColorAndValueVertex< double > > (params, pointStream, mesh); +} + +bool PoissonReconLib::Reconstruct(Parameters params, OrientedPointStream< double >* pointStream, CoredVectorMeshData< PlyValueVertex< double > >& mesh) +{ + return Execute< double, + double, + PlyValueVertex< double > > (params, pointStream, mesh); +} diff --git a/Src_CC_wrap/PoissonReconLib.h b/Src_CC_wrap/PoissonReconLib.h new file mode 100644 index 0000000..854455c --- /dev/null +++ b/Src_CC_wrap/PoissonReconLib.h @@ -0,0 +1,118 @@ +//########################################################################## +//# # +//# CLOUDCOMPARE WRAPPER: PoissonReconLib # +//# # +//# This program is free software; you can redistribute it and/or modify # +//# it under the terms of the GNU General Public License as published by # +//# the Free Software Foundation; version 2 of the License. # +//# # +//# This program is distributed in the hope that it will be useful, # +//# but WITHOUT ANY WARRANTY; without even the implied warranty of # +//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +//# GNU General Public License for more details. # +//# # +//# COPYRIGHT: Daniel Girardeau-Montaut # +//# # +//########################################################################## + +#ifndef CC_POISSON_RECON_LIB_6_11_WRAPPER +#define CC_POISSON_RECON_LIB_6_11_WRAPPER + +#include "../Src/Ply.h" +#include "../Src/Geometry.h" +#include "../Src/PointStream.h" + +//! Wrapper to use PoissonRecon (Kazhdan et. al) as a library +class PoissonReconLib +{ +public: + + //! Algorithm parameters + struct Parameters + { + //! Default initializer + Parameters(); + + //! The maximum depth of the tree that will be used for surface reconstruction + /** Running at depth d corresponds to solving on a 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. + **/ + int depth; + + //! The depth beyond which the octree will be adapted. + /** At coarser depths, the octree will be complete, containing all 2^d x 2^d x 2^d nodes. + **/ + int fullDepth; + + //! 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. + **/ + int cgDepth; + + //! The ratio between the diameter of the cube used for reconstruction and the diameter of the samples' bounding cube. + /** Specifies the factor of the bounding cube that the input samples should fit into. + **/ + float scale; + + //! The minimum number of sample points that should fall within an octree node as the octree construction is adapted to sampling density. + /** This parameter specifies the minimum number of points that should fall within an octree node. + 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. + **/ + float samplesPerNode; + + //! The importance that interpolation of the point samples is given in the formulation of the screened Poisson equation. + /** The results of the original (unscreened) Poisson Reconstruction can be obtained by setting this value to 0. + **/ + float pointWeight; + + //! The (maximum if CG) number of solver iterations + int iters; + + //! This parameter specifies the number of threads across which the solver should be parallelized + int threads; + + //! If this flag is enabled, the size of a sample's normals is used as a confidence value, affecting the sample's constribution to the reconstruction process + bool confidence; + //! If this flag is enabled, the size of a sample's normals is used as to modulate the interpolation weight + bool normalWeights; + //! If this flag is enabled, the sampling density is written out with the vertices + bool density; + + //! Pull factor for color interpolation + float colorInterp; + + //DGM: the above parameters are not documented in PoissonRecon + + bool complete; + bool showResidual; + int kernelDepth; + int maxSolveDepth; + int boundary; + + //DGM: the above parameters are hidden in PoissonRecon + + //! If this flag is enabled, the isosurface extraction does not add a planar polygon's barycenter in order to ensure that the output mesh is manifold + bool nonManifold; + //! This flag specifies the coarsest depth at which the system is to be solved + int minDepth; + //! This flag specifies the accuracy cut-off to be used for CG + float cgAccuracy; + //! This flag specifies the exponent scale for the adaptive weighting + int adaptiveExp; + }; + + //! Main entry point (shortcut to Execute) + static bool Reconstruct(Parameters params, OrientedPointStream< float >* pointStream, CoredVectorMeshData< PlyValueVertex< float > >& mesh); + //! Main entry point (shortcut to Execute) for colored clouds + static bool Reconstruct(Parameters params, OrientedPointStreamWithData< float , Point3D< unsigned char > >* pointStream, CoredVectorMeshData< PlyColorAndValueVertex< float > >& mesh); + + //! Main entry point (shortcut to Execute) + static bool Reconstruct(Parameters params, OrientedPointStream< double >* pointStream, CoredVectorMeshData< PlyValueVertex< double > >& mesh); + //! Main entry point (shortcut to Execute) for colored clouds + static bool Reconstruct(Parameters params, OrientedPointStreamWithData< double , Point3D< unsigned char > >* pointStream, CoredVectorMeshData< PlyColorAndValueVertex< double > >& mesh); + +}; + +#endif // CC_POISSON_RECON_LIB_6_11_WRAPPER