From 55f0cca8309a04927b57b31025a4ef0e54ea4de8 Mon Sep 17 00:00:00 2001 From: Daniel Girardeau-Montaut Date: Fri, 6 Dec 2019 00:09:22 +0100 Subject: [PATCH] Updated version of the PoissonRecon wrapper --- CMakeLists.txt | 5 +- Src_CC_wrap/PoissonReconLib.cpp | 1065 ++++++++++++++++++++++--------- Src_CC_wrap/PoissonReconLib.h | 67 +- 3 files changed, 776 insertions(+), 361 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a753edc..38f16ea 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,17 +13,14 @@ set( SAMPLE_PROJECT_VERSION_MAJOR 12 ) set( SAMPLE_PROJECT_VERSION_MINOR 0 ) 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/MarchingCubes.h Src/MAT.h Src/MemoryUsage.h Src/MultiGridOctreeData.h Src/MyTime.h Src/Octree.h Src/Ply.h Src/PointStream.h Src/Polynomial.h Src/PPolynomial.h Src/SparseMatrix.h ) -#list( APPEND inline_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 ) file( GLOB source_list Src_CC_wrap/*.cpp ) -#add_library( ${PROJECT_NAME} STATIC ${header_list} ${inline_list} ${source_list} ) add_library( ${PROJECT_NAME} STATIC ${header_list} ${source_list} ) # Add preprocessor 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 ) +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) diff --git a/Src_CC_wrap/PoissonReconLib.cpp b/Src_CC_wrap/PoissonReconLib.cpp index 2509aa6..3c6709a 100644 --- a/Src_CC_wrap/PoissonReconLib.cpp +++ b/Src_CC_wrap/PoissonReconLib.cpp @@ -17,38 +17,21 @@ #include "PoissonReconLib.h" -#ifdef _WIN32 -#include -#include -#endif // _WIN32 +//PoissonRecon +#include "../Src/FEMTree.h" -#ifdef WITH_OPENMP -#include -#endif #include -#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 DIMENSION 3 // The dimension of the system - -//PoissonRecon -#include "../Src/MyMiscellany.h" -#include "../Src/PPolynomial.h" -#include "../Src/FEMTree.h" -#include "../Src/Ply.h" -#include "../Src/Array.h" -#include "../Src/PointStreamData.h" -#include "../Src/Image.h" - -//#define DumpOutput(...) ((void)0) -//#include "../Src/MultiGridOctreeData.h" //only after DumpOutput has been defined! - -#define BSPLINE_DEGREE 2 +// The order of the B-Spline used to splat in data for color interpolation +static const int DATA_DEGREE = 0; +// The order of the B-Spline used to splat in the weights for density estimation +static const int WEIGHT_DEGREE = 2; +// The order of the B-Spline used to splat in the normals for constructing the Laplacian constraints +static const int NORMAL_DEGREE = 2; +// The default finite-element degree +static const int DEFAULT_FEM_DEGREE = 1; +// The dimension of the system +static const int DIMENSION = 3; PoissonReconLib::Parameters::Parameters() : depth(8) //8 @@ -75,320 +58,764 @@ PoissonReconLib::Parameters::Parameters() #endif } -static bool Reconstruct(Parameters params, const std::vector< std::pair< Pointf, TotalPointSampleDataf > >& inCorePoints, const TotalPointSampleDataf::Transform& xForm); - - -template< class Real > -XForm4x4< Real > GetPointXForm(OrientedPointStream< Real >& stream, Real scaleFactor) -{ - Point3D< Real > min, max; - stream.boundingBox(min, max); - - Real scale = std::max< Real >(max[0] - min[0], std::max< Real >(max[1] - min[1], max[2] - min[2])); - scale *= scaleFactor; - - XForm4x4< Real > tXForm = XForm4x4< Real >::Identity(); - XForm4x4< Real > sXForm = XForm4x4< Real >::Identity(); - Point3D< Real > center = (max + min) / 2; - for (int i = 0; i < 3; i++) +template +class PointData { +public: + PointData() : normal{ 0, 0, 0 }, color{ 0, 0, 0 } {} + PointData(const Real _normal[3], const Real _color[3], Real scale = 1.0) { - sXForm(i, i) = static_cast(1.0 / scale); - tXForm(3, i) = -center[i] + scale / 2; + normal[0] = scale * _normal[0]; + normal[1] = scale * _normal[1]; + normal[2] = scale * _normal[2]; + color[0] = scale * _color[0]; + color[1] = scale * _color[1]; + color[2] = scale * _color[2]; + } + + PointData operator * (Real s) const + { + return PointData(normal, color, s); + } + + PointData operator / (Real s) const + { + return PointData(normal, color, 1 / s); + } + + PointData& operator += (const PointData& d) + { + normal[0] += d.normal[0]; + normal[1] += d.normal[1]; + normal[2] += d.normal[2]; + color[0] += d.color[0]; + color[1] += d.color[1]; + color[2] += d.color[2]; + return *this; + } + PointData& operator *= (Real s) + { + normal[0] *= s; + normal[1] *= s; + normal[2] *= s; + color[0] *= s; + color[1] *= s; + color[2] *= s; + return *this; + } + +public: + Real normal[3]; + Real color[3]; +}; + +template +class Vertex : public PointData<_Real> +{ +public: + + typedef _Real Real; + + Vertex(const Point& point) + : PointData() + , point(point) + , w(0) + {} + + Vertex(const Point& point, const PointData& data, double _w = 0.0) + : PointData(data.normal, data.color) + , point(point) + , w(_w) + {} + + Vertex() : Vertex(Point(0, 0, 0)) + {} + + Vertex& operator *= (Real s) + { + PointData::operator *= (s); + point *= s; + w *= s; + + return *this; + } + + Vertex& operator /= (Real s) + { + PointData::operator *= (1 / s); + point /= s; + w /= s; + return *this; + } + + Vertex& operator+=(const Vertex& p) + { + PointData::operator += (p); + point += p.point; + w += p.w; + + return *this; + } + +public: + Point point; + double w; +}; + +template +class PointStream : public InputPointStreamWithData > +{ +public: + PointStream(const PoissonReconLib::ICloud& _cloud) + : cloud(_cloud), xform(nullptr), currentIndex(0) + {} + + void reset(void) override + { + currentIndex = 0; + } + + bool nextPoint(Point& p, PointData& d) override + { + if (currentIndex >= cloud.size()) + { + return false; + } + cloud.getPoint(currentIndex, p.coords); + + if (xform != nullptr) + { + p = (*xform) * p; + } + + if (cloud.hasNormals()) + { + cloud.getNormal(currentIndex, d.normal); + } + else + { + d.normal[0] = d.normal[1] = d.normal[2]; + } + + if (cloud.hasColors()) + { + cloud.getColor(currentIndex, d.color); + } + else + { + d.color[0] = d.color[1] = d.color[2]; + } + + currentIndex++; + return true; + } + +public: + const PoissonReconLib::ICloud& cloud; + XForm* xform; + size_t currentIndex; +}; + +template +struct FEMTreeProfiler { + FEMTree& tree; + double t; + + FEMTreeProfiler(FEMTree& t) : tree(t) {} + void start(void) { + t = Time(), FEMTree::ResetLocalMemoryUsage(); + } + void dumpOutput(const char* header) const { + FEMTree::MemoryUsage(); + //if (header) { + // utility::LogDebug("{} {} (s), {} (MB) / {} (MB) / {} (MB)", header, + // Time() - t, + // FEMTree::LocalMemoryUsage(), + // FEMTree::MaxMemoryUsage(), + // MemoryInfo::PeakMemoryUsageMB()); + //} + //else { + // utility::LogDebug("{} (s), {} (MB) / {} (MB) / {} (MB)", Time() - t, + // FEMTree::LocalMemoryUsage(), + // FEMTree::MaxMemoryUsage(), + // MemoryInfo::PeakMemoryUsageMB()); + //} + } +}; + +template +XForm GetBoundingBoxXForm(Point min, + Point max, + Real scaleFactor) { + Point center = (max + min) / 2; + Real scale = max[0] - min[0]; + for (unsigned int d = 1; d < Dim; d++) { + scale = std::max(scale, max[d] - min[d]); + } + scale *= scaleFactor; + for (unsigned int i = 0; i < Dim; i++) { + center[i] -= scale / 2; + } + XForm tXForm = XForm::Identity(), + sXForm = XForm::Identity(); + for (unsigned int i = 0; i < Dim; i++) { + sXForm(i, i) = (Real)(1. / scale), tXForm(Dim, i) = -center[i]; } return sXForm * tXForm; } -template< class Real, int Degree, BoundaryType BType, class Vertex > -bool Execute( PoissonReconLib::Parameters params, - OrientedPointStream< Real >* pointStream, - bool withColors, - CoredVectorMeshData< Vertex >& mesh, - XForm4x4< Real >& iXForm) +template +XForm GetBoundingBoxXForm(Point min, + Point max, + Real width, + Real scaleFactor, + int& depth) { + // Get the target resolution (along the largest dimension) + Real resolution = (max[0] - min[0]) / width; + for (unsigned int d = 1; d < Dim; d++) { + resolution = std::max(resolution, (max[d] - min[d]) / width); + } + resolution *= scaleFactor; + depth = 0; + while ((1 << depth) < resolution) { + depth++; + } + + Point center = (max + min) / 2; + Real scale = (1 << depth) * width; + + for (unsigned int i = 0; i < Dim; i++) { + center[i] -= scale / 2; + } + XForm tXForm = XForm::Identity(), + sXForm = XForm::Identity(); + for (unsigned int i = 0; i < Dim; i++) { + sXForm(i, i) = (Real)(1. / scale), tXForm(Dim, i) = -center[i]; + } + return sXForm * tXForm; +} + +template +XForm GetPointXForm(InputPointStream& stream, + Real width, + Real scaleFactor, + int& depth) { + Point min, max; + stream.boundingBox(min, max); + return GetBoundingBoxXForm(min, max, width, scaleFactor, depth); +} + +template +XForm GetPointXForm(InputPointStream& stream, + Real scaleFactor) { + Point min, max; + stream.boundingBox(min, max); + return GetBoundingBoxXForm(min, max, scaleFactor); +} + +template +struct ConstraintDual { + Real target, weight; + ConstraintDual(Real t, Real w) : target(t), weight(w) {} + CumulativeDerivativeValues operator()( + const Point& p) const { + return CumulativeDerivativeValues(target * weight); + }; +}; + +template +struct SystemDual { + Real weight; + SystemDual(Real w) : weight(w) {} + CumulativeDerivativeValues operator()( + const Point& p, + const CumulativeDerivativeValues& dValues) const { + return dValues * weight; + }; + CumulativeDerivativeValues operator()( + const Point& p, + const CumulativeDerivativeValues& dValues) const { + return dValues * weight; + }; +}; + +template +struct SystemDual { + typedef double Real; + Real weight; + SystemDual(Real w) : weight(w) {} + CumulativeDerivativeValues operator()( + const Point& p, + const CumulativeDerivativeValues& dValues) const { + return dValues * weight; + }; +}; + +template + void ExtractMesh( + float datax, + bool linear_fit, + UIntPack, + std::tuple, + FEMTree& tree, + const DenseNodeData>& solution, + Real isoValue, + const std::vector::PointSample>* samples, + std::vector< PointData >* sampleData, + const typename FEMTree::template DensityEstimator* + density, + const SetVertexFunction& SetVertex, + XForm iXForm, + PoissonReconLib::IMesh& out_mesh) { - typedef typename Octree< Real >::template DensityEstimator< WEIGHT_DEGREE > DensityEstimator; - typedef typename Octree< Real >::template InterpolationInfo< false > InterpolationInfo; - typedef OrientedPointStreamWithData< Real, Point3D< Real > > PointStreamWithData; - typedef TransformedOrientedPointStream< Real > XPointStream; - typedef TransformedOrientedPointStreamWithData< Real, Point3D< Real > > XPointStreamWithData; - Reset< Real >(); + static const int Dim = sizeof...(FEMSigs); + typedef UIntPack Sigs; + static const unsigned int DataSig = + FEMDegreeAndBType::Signature; + typedef typename FEMTree::template DensityEstimator + DensityEstimator; - //DGM: do this begore initializing the octree!! - OctNode< TreeNodeData >::SetAllocator(MEMORY_ALLOCATOR_BLOCK_SIZE); + FEMTreeProfiler profiler(tree); - Octree< Real > tree; - tree.threads = params.threads; + CoredMeshData* mesh; + mesh = new CoredVectorMeshData(); - if (params.maxSolveDepth == 0) - params.maxSolveDepth = params.depth; + bool non_manifold = true; + bool polygon_mesh = false; - 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); - - try - { - XForm4x4< Real > xForm = XForm4x4< Real >::Identity(); - { - xForm = GetPointXForm(*pointStream, static_cast(params.scale)); + profiler.start(); + typename IsoSurfaceExtractor::IsoStats isoStats; + if (sampleData) { + SparseNodeData, Real>, + IsotropicUIntPack> + _sampleData = + tree.template setMultiDepthDataField( + *samples, *sampleData, (DensityEstimator*)NULL); + for (const RegularTreeNode* + n = tree.tree().nextNode(); + n; n = tree.tree().nextNode(n)) { + ProjectiveData, Real>* clr = _sampleData(n); + if (clr) (*clr) *= (Real)pow(datax, tree.depth(n)); } - iXForm = xForm.inverse(); + isoStats = IsoSurfaceExtractor::template Extract< + PointData >(Sigs(), UIntPack(), + UIntPack(), tree, density, &_sampleData, + solution, isoValue, *mesh, SetVertex, !linear_fit, + !non_manifold, polygon_mesh, false); + } + else { + isoStats = IsoSurfaceExtractor::template Extract< + PointData >(Sigs(), UIntPack(), + UIntPack(), tree, density, NULL, solution, + isoValue, *mesh, SetVertex, !linear_fit, + !non_manifold, polygon_mesh, false); + } - std::vector< typename Octree< Real >::PointSample > samples; - std::vector< ProjectiveData< Point3D< Real >, Real > > sampleData; + mesh->resetIterator(); + for (size_t vidx = 0; vidx < mesh->outOfCorePointCount(); ++vidx) { + Vertex v; + mesh->nextOutOfCorePoint(v); + v.point = iXForm * v.point; - if (withColors) + out_mesh.addVertex(v.point.coords); + out_mesh.addNormal(v.normal); + out_mesh.addColor(v.color); + out_mesh.addDensity(v.w); + } + for (size_t tidx = 0; tidx < mesh->polygonCount(); ++tidx) { + std::vector> triangle; + mesh->nextPolygon(triangle); + if (triangle.size() == 3) { - XPointStreamWithData _pointStream(xForm, *((PointStreamWithData*)pointStream)); - int pointCount = tree.template init< Point3D< Real > >( - _pointStream, - params.depth, - params.confidence, - samples, - &sampleData); + out_mesh.addTriangle(triangle[0].idx, triangle[1].idx, triangle[2].idx); } else { - XPointStream _pointStream(xForm, *pointStream); - int pointCount = tree.template init< Point3D< Real > >( - _pointStream, - params.depth, - params.confidence, - samples, - 0); - } - - DenseNodeData< Real, Degree > solution; - DensityEstimator* density = NULL; - { - int solveDepth = params.depth; - - tree.resetNodeIndices(); - - // Get the kernel density estimator [If discarding, compute anew. Otherwise, compute once.] - density = tree.template setDensityEstimator< WEIGHT_DEGREE >(samples, kernelDepth, params.samplesPerNode); - - // Transform the Hermite samples into a vector field - Real pointWeightSum = 0; - SparseNodeData< Point3D< Real >, NORMAL_DEGREE > normalInfo = tree.template setNormalField< NORMAL_DEGREE >(samples, *density, pointWeightSum, BType == BOUNDARY_NEUMANN); - - if (!params.density) - { - delete density; - density = 0; - } - - // Trim the tree and prepare for multigrid - { - std::vector< int > indexMap; - - tree.template inalizeForBroodedMultigrid< NORMAL_DEGREE, Degree, BType >(params.fullDepth, typename Octree< Real >::template HasNormalDataFunctor< NORMAL_DEGREE >(normalInfo), &indexMap); - - normalInfo.remapIndices(indexMap); - if (params.density) - { - density->remapIndices(indexMap); - } - } - - // Add the FEM constraints - DenseNodeData< Real, Degree > constraints; - { - constraints = tree.template initDenseNodeData< Degree >(); - tree.template addFEMConstraints< Degree, BType, NORMAL_DEGREE, BType >(FEMVFConstraintFunctor< NORMAL_DEGREE, BType, Degree, BType >(1., 0.), normalInfo, constraints, solveDepth); - } - - // Free up the normal info [If we don't need it for subseequent iterations.] - normalInfo.clear(); - - // Add the interpolation constraints - InterpolationInfo* iInfo = NULL; - if (params.pointWeight > 0) - { - Real targetValue = static_cast(0.5); - iInfo = new InterpolationInfo(tree, samples, targetValue, params.adaptiveExp, static_cast(params.pointWeight) * pointWeightSum, (Real)0); - tree.template addInterpolationConstraints< Degree, BType >(*iInfo, constraints, solveDepth); - } - - //DumpOutput("Leaf Nodes / Active Nodes / Ghost Nodes: %d / %d / %d\n", (int)tree.leaves(), (int)tree.nodes(), (int)tree.ghostNodes()); - //DumpOutput("Memory Usage: %.3f MB\n", float(MemoryInfo::Usage()) / (1 << 20)); - - // Solve the linear system - double lowResIterMultiplier = 1.0; - { - typename Octree< Real >::SolverInfo solverInfo; - solverInfo.cgDepth = params.cgDepth; - solverInfo.iters = params.iters; - solverInfo.cgAccuracy = params.cgAccuracy; - solverInfo.verbose = false; - solverInfo.showResidual = params.showResidual; - solverInfo.lowResIterMultiplier = std::max< double >(1.0, lowResIterMultiplier); - solution = tree.template solveSystem< Degree, BType >(FEMSystemFunctor< Degree, BType >(0, 1., 0), iInfo, constraints, solveDepth, solverInfo); - if (iInfo) - { - delete iInfo; - iInfo = NULL; - } - } - } - - Real isoValue = 0; - { - double valueSum = 0, weightSum = 0; - typename Octree< Real >::template MultiThreadedEvaluator< Degree, BType > evaluator(&tree, solution, params.threads); - -#pragma omp parallel for num_threads( params.threads ) reduction( + : valueSum , weightSum ) - for (int j = 0; j < samples.size(); j++) - { - const ProjectiveData< OrientedPoint3D< Real >, Real >& sample = samples[j].sample; - if (sample.weight > 0) - { - weightSum += sample.weight; - valueSum += evaluator.value(sample.data.p / sample.weight, omp_get_thread_num(), samples[j].node) * sample.weight; - } - } - isoValue = static_cast(valueSum / weightSum); - //DumpOutput("Iso-Value: %e\n", isoValue); - } - - SparseNodeData< ProjectiveData< Point3D< Real >, Real >, DATA_DEGREE >* colorData = NULL; - if (withColors) - { - colorData = new SparseNodeData< ProjectiveData< Point3D< Real >, Real >, DATA_DEGREE >(); - *colorData = tree.template setDataField< DATA_DEGREE, false >(samples, sampleData, (DensityEstimator*)NULL); - for (const OctNode< TreeNodeData >* n = tree.tree().nextNode(); n; n = tree.tree().nextNode(n)) - { - ProjectiveData< Point3D< Real >, Real >* clr = (*colorData)(n); - if (clr) (*clr) *= static_cast(pow(params.colorInterp, tree.depth(n))); - } - } - - bool linearFit = false; - bool polygonMesh = false; - tree.template getMCIsoSurface< Degree, BType, WEIGHT_DEGREE, DATA_DEGREE >( - density, - colorData, - solution, - isoValue, - mesh, - !linearFit, - !params.nonManifold, - polygonMesh); - - if (density) - { - delete density; - density = NULL; - } - - //DumpOutput("Vertices / Polygons: %d / %d\n", mesh.outOfCorePointCount() + mesh.inCorePoints.size(), mesh.polygonCount()); - if (colorData) - { - delete colorData; - colorData = NULL; + assert(false); } } - catch (const std::bad_alloc&) + + delete mesh; +} + +template +static Real ComputeNorm(const Real vec[3]) +{ + return sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); +} + +template +static void Execute(PointStream& pointStream, + PoissonReconLib::IMesh& out_mesh, + int depth, + Real width, + float scale, + bool linear_fit, + UIntPack) { + static const int Dim = sizeof...(FEMSigs); + typedef UIntPack Sigs; + typedef UIntPack::Degree...> Degrees; + typedef UIntPack::BType, + 1>::BType>::Signature...> + NormalSigs; + typedef typename FEMTree::template DensityEstimator + DensityEstimator; + typedef typename FEMTree::template InterpolationInfo + InterpolationInfo; + + XForm xForm, iXForm; + xForm = XForm::Identity(); + + float datax = 32.f; + int base_depth = 0; + int base_v_cycles = 1; + float confidence = 0.f; + float point_weight = 2.f * DEFAULT_FEM_DEGREE; + float confidence_bias = 0.f; + float samples_per_node = 1.5f; + float cg_solver_accuracy = 1e-3f; + int full_depth = 5; + int iters = 8; + bool exact_interpolation = false; + + double startTime = Time(); + Real isoValue = 0; + + FEMTree tree(MEMORY_ALLOCATOR_BLOCK_SIZE); + FEMTreeProfiler profiler(tree); + + size_t pointCount; + + Real pointWeightSum; + std::vector::PointSample> samples; + std::vector< PointData > sampleData; + DensityEstimator* density = NULL; + SparseNodeData, NormalSigs>* normalInfo = NULL; + Real targetValue = (Real)0.5; + + // Read in the samples (and color data) { - //not enough memory + if (width > 0) { + xForm = GetPointXForm(pointStream, width, + static_cast(scale > 0 ? scale : 1.0), + depth) * + xForm; + } + else { + xForm = scale > 0 ? GetPointXForm(pointStream, + (Real)scale) * + xForm + : xForm; + } + + pointStream.xform = &xForm; + + { + auto ProcessDataWithConfidence = [&](const Point& p, + PointData& d) { + Real l = ComputeNorm(d.normal); + if (!l || l != l) return (Real)-1.; + return (Real)pow(l, confidence); + }; + auto ProcessData = [](const Point& p, PointData& d) { + Real l = ComputeNorm(d.normal); + if (!l || l != l) return (Real)-1.; + d.normal[0] /= l; + d.normal[1] /= l; + d.normal[2] /= l; + return (Real)1.; + }; + if (confidence > 0) { + pointCount = FEMTreeInitializer::template Initialize< + PointData>(tree.spaceRoot(), pointStream, depth, + samples, sampleData, true, + tree.nodeAllocators[0], tree.initializer(), + ProcessDataWithConfidence); + } + else { + pointCount = FEMTreeInitializer::template Initialize< + PointData>(tree.spaceRoot(), pointStream, depth, + samples, sampleData, true, + tree.nodeAllocators[0], tree.initializer(), + ProcessData); + } + } + iXForm = xForm.inverse(); + + //utility::LogDebug("Input Points / Samples: {} / {}", pointCount, + // samples.size()); + } + + int kernelDepth = depth - 2; + if (kernelDepth < 0) { + //utility::LogError( + // "[CreateFromPointCloudPoisson] depth (={}) has to be >= 2", + // depth); + } + + DenseNodeData solution; + { + DenseNodeData constraints; + InterpolationInfo* iInfo = NULL; + int solveDepth = depth; + + tree.resetNodeIndices(); + + // Get the kernel density estimator + { + profiler.start(); + density = tree.template setDensityEstimator( + samples, kernelDepth, samples_per_node, 1); + profiler.dumpOutput("# Got kernel density:"); + } + + // Transform the Hermite samples into a vector field + { + profiler.start(); + normalInfo = new SparseNodeData, NormalSigs>(); + std::function, Point&)> + ConversionFunction = + [](PointData in, Point& out) { + // Point n = in.template data<0>(); + Point n(in.normal[0], in.normal[1], in.normal[2]); + 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, Point&, Real&)> + ConversionAndBiasFunction = [&](PointData in, + Point& out, + Real& bias) { + // Point n = in.template data<0>(); + Point n(in.normal[0], in.normal[1], in.normal[2]); + 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) * confidence_bias / + log(1 << (Dim - 1))); + return true; + }; + if (confidence_bias > 0) { + *normalInfo = tree.setDataField( + NormalSigs(), samples, sampleData, density, + pointWeightSum, ConversionAndBiasFunction); + } + else { + *normalInfo = tree.setDataField( + NormalSigs(), samples, sampleData, density, + pointWeightSum, ConversionFunction); + } + ThreadPool::Parallel_for(0, normalInfo->size(), + [&](unsigned int, size_t i) { + (*normalInfo)[i] *= (Real)-1.; + }); + profiler.dumpOutput("# Got normal field:"); + //utility::LogDebug("Point weight / Estimated Area: {:e} / {:e}", + // pointWeightSum, pointCount * pointWeightSum); + } + + // Trim the tree and prepare for multigrid + { + profiler.start(); + constexpr int MAX_DEGREE = NORMAL_DEGREE > Degrees::Max() + ? NORMAL_DEGREE + : Degrees::Max(); + tree.template finalizeForMultigrid( + full_depth, + typename FEMTree::template HasNormalDataFunctor< + NormalSigs>(*normalInfo), + normalInfo, density); + profiler.dumpOutput("# Finalized tree:"); + } + + // Add the FEM constraints + { + profiler.start(); + constraints = tree.initDenseNodeData(Sigs()); + typename FEMIntegrator::template Constraint< + Sigs, IsotropicUIntPack, NormalSigs, + IsotropicUIntPack, Dim> + F; + unsigned int derivatives2[Dim]; + for (unsigned int d = 0; d < Dim; d++) derivatives2[d] = 0; + typedef IsotropicUIntPack Derivatives1; + typedef IsotropicUIntPack Derivatives2; + for (unsigned int d = 0; d < Dim; d++) { + unsigned int derivatives1[Dim]; + for (unsigned int dd = 0; dd < Dim; dd++) + derivatives1[dd] = dd == d ? 1 : 0; + F.weights[d] + [TensorDerivatives::Index(derivatives1)] + [TensorDerivatives::Index( + derivatives2)] = 1; + } + tree.addFEMConstraints(F, *normalInfo, constraints, solveDepth); + profiler.dumpOutput("# Set FEM constraints:"); + } + + // Free up the normal info + delete normalInfo, normalInfo = NULL; + + // Add the interpolation constraints + if (point_weight > 0) { + profiler.start(); + if (exact_interpolation) { + iInfo = FEMTree:: + template InitializeExactPointInterpolationInfo( + tree, samples, + ConstraintDual( + targetValue, + (Real)point_weight * pointWeightSum), + SystemDual((Real)point_weight * + pointWeightSum), + true, false); + } + else { + iInfo = FEMTree:: + template InitializeApproximatePointInterpolationInfo< + Real, 0>( + tree, samples, + ConstraintDual( + targetValue, + (Real)point_weight * pointWeightSum), + SystemDual((Real)point_weight * + pointWeightSum), + true, 1); + } + tree.addInterpolationConstraints(constraints, solveDepth, *iInfo); + profiler.dumpOutput("#Set point constraints:"); + } + + //utility::LogDebug( + // "Leaf Nodes / Active Nodes / Ghost Nodes: {} / {} / {}", + // tree.leaves(), tree.nodes(), tree.ghostNodes()); + //utility::LogDebug("Memory Usage: {:.3f} MB", + // float(MemoryInfo::Usage()) / (1 << 20)); + + // Solve the linear system + { + profiler.start(); + typename FEMTree::SolverInfo sInfo; + sInfo.cgDepth = 0, sInfo.cascadic = true, sInfo.vCycles = 1, + sInfo.iters = iters, sInfo.cgAccuracy = cg_solver_accuracy, + sInfo.verbose = false/* utility::Logger::i().verbosity_level_ == + utility::VerbosityLevel::Debug */, + sInfo.showResidual = false/*utility::Logger::i().verbosity_level_ == + utility::VerbosityLevel::Debug*/, + sInfo.showGlobalResidual = SHOW_GLOBAL_RESIDUAL_NONE, + sInfo.sliceBlockSize = 1; + sInfo.baseDepth = base_depth, sInfo.baseVCycles = base_v_cycles; + typename FEMIntegrator::template System> + F({ 0., 1. }); + solution = tree.solveSystem(Sigs(), F, constraints, solveDepth, + sInfo, iInfo); + profiler.dumpOutput("# Linear system solved:"); + if (iInfo) delete iInfo, iInfo = NULL; + } + } + + { + profiler.start(); + double valueSum = 0, weightSum = 0; + typename FEMTree::template MultiThreadedEvaluator + evaluator(&tree, solution); + std::vector valueSums(ThreadPool::NumThreads(), 0), + weightSums(ThreadPool::NumThreads(), 0); + ThreadPool::Parallel_for( + 0, samples.size(), [&](unsigned int thread, size_t j) { + ProjectiveData, 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 < valueSums.size(); t++) + valueSum += valueSums[t], weightSum += weightSums[t]; + isoValue = (Real)(valueSum / weightSum); + profiler.dumpOutput("Got average:"); + //utility::LogDebug("Iso-Value: {:e} = {:e} / {:e}", isoValue, valueSum, + // weightSum); + } + + auto SetVertex = [](Vertex& v, Point p, double w, + PointData d) { + v = Vertex(p, d, w); + }; + ExtractMesh, Real>( + datax, linear_fit, UIntPack(), + std::tuple(), tree, solution, isoValue, &samples, + &sampleData, density, SetVertex, iXForm, out_mesh); + + if (density) + { + delete density; + density = nullptr; + } + //utility::LogDebug("# Total Solve: {:9.1f} (s), {:9.1f} (MB)", + // Time() - startTime, FEMTree::MaxMemoryUsage()); +} + + +bool PoissonReconLib::Reconstruct( const Parameters& params, + const ICloud& inCloud, + IMesh& outMesh ) +{ + if (!inCloud.hasNormals()) + { + //we need normals return false; } - catch (std::exception e) + +#ifdef WITH_OPENMP + ThreadPool::Init((ThreadPool::ParallelType)(int)ThreadPool::OPEN_MP, + std::thread::hardware_concurrency()); +#else + ThreadPool::Init((ThreadPool::ParallelType)(int)ThreadPool::THREAD_POOL, + std::thread::hardware_concurrency()); +#endif + + PointStream pointStream(inCloud); + + switch (params.boundary) { - //not enough memory - return false; - } - catch (...) - { - //not enough memory - return false; + case Parameters::FREE: + typedef IsotropicUIntPack::Signature> FEMSigsFree; + Execute(pointStream, outMesh, params.depth, params.width, params.scale, params.linear_fit, FEMSigsFree()); + break; + case Parameters::DIRICHLET: + typedef IsotropicUIntPack::Signature> FEMSigsDirichlet; + Execute(pointStream, outMesh, params.depth, params.width, params.scale, params.linear_fit, FEMSigsDirichlet()); + break; + case Parameters::NEUMANN: + typedef IsotropicUIntPack::Signature> FEMSigsNeumann; + Execute(pointStream, outMesh, params.depth, params.width, params.scale, params.linear_fit, FEMSigsNeumann()); + break; + default: + assert(false); + break; } + ThreadPool::Terminate(); + return true; } - -bool PoissonReconLib::Reconstruct( Parameters params, - OrientedPointStreamWithData< float , Point3D< float > >* pointStream, - CoredVectorMeshData< PlyColorAndValueVertex< float > >& mesh, - XForm4x4< float >& iXForm) -{ - switch (params.boundary) - { - case Parameters::FREE: - return Execute< float, BSPLINE_DEGREE, BOUNDARY_FREE, PlyColorAndValueVertex< float > >(params, pointStream, true, mesh, iXForm); - case Parameters::DIRICHLET: - return Execute< float, BSPLINE_DEGREE, BOUNDARY_DIRICHLET, PlyColorAndValueVertex< float > >(params, pointStream, true, mesh, iXForm); - case Parameters::NEUMANN: - return Execute< float, BSPLINE_DEGREE, BOUNDARY_NEUMANN, PlyColorAndValueVertex< float > >(params, pointStream, true, mesh, iXForm); - default: - assert(false); - break; - } - - return false; -} - -bool PoissonReconLib::Reconstruct( Parameters params, - OrientedPointStream< float >* pointStream, - CoredVectorMeshData< PlyValueVertex< float > >& mesh, - XForm4x4< float >& iXForm) -{ - switch (params.boundary) - { - case Parameters::FREE: - return Execute< float, BSPLINE_DEGREE, BOUNDARY_FREE, PlyValueVertex< float > >(params, pointStream, false, mesh, iXForm); - case Parameters::DIRICHLET: - return Execute< float, BSPLINE_DEGREE, BOUNDARY_DIRICHLET, PlyValueVertex< float > >(params, pointStream, false, mesh, iXForm); - case Parameters::NEUMANN: - return Execute< float, BSPLINE_DEGREE, BOUNDARY_NEUMANN, PlyValueVertex< float > >(params, pointStream, false, mesh, iXForm); - default: - assert(false); - break; - } - - return false; -} - -bool PoissonReconLib::Reconstruct( Parameters params, - OrientedPointStreamWithData< double , Point3D< double > >* pointStream, - CoredVectorMeshData< PlyColorAndValueVertex< double > >& mesh, - XForm4x4< double >& iXForm) -{ - switch (params.boundary) - { - case Parameters::FREE: - return Execute< double, BSPLINE_DEGREE, BOUNDARY_FREE, PlyColorAndValueVertex< double > >(params, pointStream, true, mesh, iXForm); - case Parameters::DIRICHLET: - return Execute< double, BSPLINE_DEGREE, BOUNDARY_DIRICHLET, PlyColorAndValueVertex< double > >(params, pointStream, true, mesh, iXForm); - case Parameters::NEUMANN: - return Execute< double, BSPLINE_DEGREE, BOUNDARY_NEUMANN, PlyColorAndValueVertex< double > >(params, pointStream, true, mesh, iXForm); - default: - assert(false); - break; - } - - return false; -} - -bool PoissonReconLib::Reconstruct( Parameters params, - OrientedPointStream< double >* pointStream, - CoredVectorMeshData< PlyValueVertex< double > >& mesh, - XForm4x4< double >& iXForm) -{ - switch (params.boundary) - { - case Parameters::FREE: - return Execute< double, BSPLINE_DEGREE, BOUNDARY_FREE, PlyValueVertex< double > >(params, pointStream, false, mesh, iXForm); - case Parameters::DIRICHLET: - return Execute< double, BSPLINE_DEGREE, BOUNDARY_DIRICHLET, PlyValueVertex< double > >(params, pointStream, false, mesh, iXForm); - case Parameters::NEUMANN: - return Execute< double, BSPLINE_DEGREE, BOUNDARY_NEUMANN, PlyValueVertex< double > >(params, pointStream, false, mesh, iXForm); - default: - assert(false); - break; - } - - return false; -} diff --git a/Src_CC_wrap/PoissonReconLib.h b/Src_CC_wrap/PoissonReconLib.h index f5f9722..0d76e8c 100644 --- a/Src_CC_wrap/PoissonReconLib.h +++ b/Src_CC_wrap/PoissonReconLib.h @@ -18,15 +18,6 @@ #ifndef CC_POISSON_RECON_LIB_WRAPPER #define CC_POISSON_RECON_LIB_WRAPPER -#include "../Src/PointStream.h" - -//#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 DIMENSION 3 // The dimension of the system - //! Wrapper to use PoissonRecon (Kazhdan et. al) as a library class PoissonReconLib { @@ -86,7 +77,7 @@ public: //! Pull factor for color interpolation float colorInterp; - //DGM: the above parameters are not documented in PoissonRecon + //DGM: the parameters below are not documented in PoissonRecon bool showResidual; int kernelDepth; @@ -95,7 +86,7 @@ public: enum BoundaryType { FREE, DIRICHLET, NEUMANN }; BoundaryType boundary; - //DGM: the above parameters are hidden in PoissonRecon + //DGM: the below 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; @@ -103,38 +94,38 @@ public: float cgAccuracy; //! This flag specifies the exponent scale for the adaptive weighting int adaptiveExp; + + float width = 0; + bool linear_fit = false; }; - static const int Degree = DEFAULT_FEM_DEGREE; - static const BoundaryType BType = DEFAULT_FEM_BOUNDARY; - typedef IsotropicUIntPack< DIMENSION, FEMDegreeAndBType< Degree, BType >::Signature > FEMSigs; - static const int Dim = sizeof(FEMSigs); + template + class ICloud + { + public: + virtual size_t size() const = 0; + virtual bool hasNormals() const = 0; + virtual bool hasColors() const = 0; + virtual void getPoint(size_t index, Real* coords) const = 0; + virtual void getNormal(size_t index, Real* coords) const = 0; + virtual void getColor(size_t index, Real* rgb) const = 0; + }; - typedef Point< float, Dim > Pointf; - typedef Point< double, Dim > Pointd; + template + class IMesh + { + public: + virtual void addVertex(const Real* coords) = 0; + virtual void addNormal(const Real* coords) = 0; + virtual void addColor(const Real* rgb) = 0; + virtual void addDensity(double d) = 0; + virtual void addTriangle(size_t i1, size_t i2, size_t i3) = 0; + }; - typedef PointStreamNormal< float, Dim > NormalPointSampleDataf; - typedef PointStreamColor< float > SampleDataf; - typedef MultiPointStreamData< float, SampleDataf > AdditionalPointSampleDataf; - typedef MultiPointStreamData< float, NormalPointSampleDataf, AdditionalPointSampleDataf > TotalPointSampleDataf; - - static bool Reconstruct(Parameters params, const std::vector< std::pair< Pointf, TotalPointSampleDataf > >& inCorePoints, const TotalPointSampleDataf::Transform& xForm); - - //typedef PointStreamColor< double > SampleDatad; - //BinaryInputPointStreamWithData - - //typedef PointStreamNormal< float, 3 > NormalPointSampleData; - - ////! Main entry point (shortcut to Execute) - //static bool Reconstruct(Parameters params, OrientedPointStream< float >* pointStream, CoredVectorMeshData< PlyValueVertex< float > >& mesh, XForm4x4< float >& iXForm); - ////! Main entry point (shortcut to Execute) for colored clouds - //static bool Reconstruct(Parameters params, OrientedPointStreamWithData< float, Point3D< float > >* pointStream, CoredVectorMeshData< PlyColorAndValueVertex< float > >& mesh, XForm4x4< float >& iXForm); - - ////! Main entry point (shortcut to Execute) - //static bool Reconstruct(Parameters params, OrientedPointStream< double >* pointStream, CoredVectorMeshData< PlyValueVertex< double > >& mesh, XForm4x4< double >& iXForm); - ////! Main entry point (shortcut to Execute) for colored clouds - //static bool Reconstruct(Parameters params, OrientedPointStreamWithData< double, Point3D< double > >* pointStream, CoredVectorMeshData< PlyColorAndValueVertex< double > >& mesh, XForm4x4< double >& iXForm); + static bool Reconstruct(const Parameters& params, + const PoissonReconLib::ICloud& inCloud, + PoissonReconLib::IMesh& ouMesh); }; #endif // CC_POISSON_RECON_LIB_6_11_WRAPPER