From a7080f101f4317aa047096bbae1bfeded9d72139 Mon Sep 17 00:00:00 2001 From: Gernot Date: Tue, 26 Nov 2019 19:32:41 +0100 Subject: [PATCH] Poisson Surface Reconstruction (#1317) * Poisson Surface Reconstruction * addressing CI warnings * disable -Werror * merge upstream * remove poisson files to setup submodule * PoissonRecon as submodule * Merge remote-tracking branch 'upstream/master' into poisson * fix and enable warning, and misc * use /W3 without /WX for windows * addressing review comments * added method to remove low density vertices --- .gitignore | 1 + .gitmodules | 3 + 3rdparty/CMakeLists.txt | 5 + 3rdparty/PoissonRecon | 1 + 3rdparty/README.txt | 4 + .../surface_reconstruction_poisson.py | 47 ++ examples/Python/Misc/meshes.py | 10 + src/CMakeLists.txt | 8 +- .../Geometry/SurfaceReconstructionPoisson.cpp | 763 ++++++++++++++++++ src/Open3D/Geometry/TriangleMesh.cpp | 61 ++ src/Open3D/Geometry/TriangleMesh.h | 44 + .../open3d_pybind/geometry/trianglemesh.cpp | 53 ++ src/UnitTest/Geometry/TriangleMesh.cpp | 245 +++++- src/UnitTest/TestUtility/UnitTest.cpp | 22 +- src/UnitTest/TestUtility/UnitTest.h | 16 +- 15 files changed, 1262 insertions(+), 21 deletions(-) create mode 160000 3rdparty/PoissonRecon create mode 100644 examples/Python/Advanced/surface_reconstruction_poisson.py create mode 100644 src/Open3D/Geometry/SurfaceReconstructionPoisson.cpp diff --git a/.gitignore b/.gitignore index 1983ea0bcaf..aacac320fa9 100644 --- a/.gitignore +++ b/.gitignore @@ -41,6 +41,7 @@ compile_commands.json examples/Python/ReconstructionSystem/dataset/ examples/TestData/Armadillo.ply examples/TestData/Bunny.ply +examples/TestData/eagle.ply examples/TestData/voxelized.ply examples/Python/Basic/voxel_grid_test.ply examples/Python/Benchmark/testdata/ diff --git a/.gitmodules b/.gitmodules index 302ac25f00a..05ddb7184c5 100644 --- a/.gitmodules +++ b/.gitmodules @@ -28,3 +28,6 @@ [submodule "3rdparty/libjpeg-turbo/libjpeg-turbo"] path = 3rdparty/libjpeg-turbo/libjpeg-turbo url = https://github.com/libjpeg-turbo/libjpeg-turbo.git +[submodule "3rdparty/PoissonRecon/Open3D-PoissonRecon"] + path = 3rdparty/PoissonRecon + url = https://github.com/intel-isl/Open3D-PoissonRecon.git diff --git a/3rdparty/CMakeLists.txt b/3rdparty/CMakeLists.txt index a884d69979c..1f3d71b34c3 100644 --- a/3rdparty/CMakeLists.txt +++ b/3rdparty/CMakeLists.txt @@ -313,6 +313,10 @@ endif() set(fmt_INCLUDE_DIRS ${CMAKE_CURRENT_SOURCE_DIR}/fmt/include) INSTALL_HEADERS(fmt) +# PoissonRecon +set(poisson_INCLUDE_DIRS ${CMAKE_CURRENT_SOURCE_DIR}/) + + list(APPEND 3RDPARTY_INCLUDE_DIRS # ${dirent_INCLUDE_DIRS} # fails on Linux, seems to require Windows headers ${EIGEN3_INCLUDE_DIRS} @@ -333,6 +337,7 @@ list(APPEND 3RDPARTY_INCLUDE_DIRS ${googletest_INCLUDE_DIRS} ${fmt_INCLUDE_DIRS} ${k4a_INCLUDE_DIRS} + ${poisson_INCLUDE_DIRS} ) # set 3RDPARTY_INCLUDE_DIRS_AT_INSTALL diff --git a/3rdparty/PoissonRecon b/3rdparty/PoissonRecon new file mode 160000 index 00000000000..fd273ea8c77 --- /dev/null +++ b/3rdparty/PoissonRecon @@ -0,0 +1 @@ +Subproject commit fd273ea8c77a36973d6565a495c9969ccfb12d3b diff --git a/3rdparty/README.txt b/3rdparty/README.txt index f3022cc5d10..ce319b8ecec 100644 --- a/3rdparty/README.txt +++ b/3rdparty/README.txt @@ -72,3 +72,7 @@ https://github.com/syoyo/tinyobjloader pybind11 2.2 BSD license Python binding for C++11 https://github.com/pybind/pybind11 +-------------------------------------------------------------------------------- +PoissonReco 12.0 BSD license +Poisson Surface Reconstruction +https://github.com/mkazhdan/PoissonRecon diff --git a/examples/Python/Advanced/surface_reconstruction_poisson.py b/examples/Python/Advanced/surface_reconstruction_poisson.py new file mode 100644 index 00000000000..0042a4c1c3c --- /dev/null +++ b/examples/Python/Advanced/surface_reconstruction_poisson.py @@ -0,0 +1,47 @@ +# Open3D: www.open3d.org +# The MIT License (MIT) +# See license file or visit www.open3d.org for details + +# examples/Python/Advanced/surface_reconstruction_poisson.py + +import open3d as o3d +import numpy as np +import matplotlib.pyplot as plt +import os + +import sys + +sys.path.append( + os.path.join(os.path.dirname(os.path.realpath(__file__)), "../Misc")) +import meshes + +if __name__ == "__main__": + o3d.utility.set_verbosity_level(o3d.utility.VerbosityLevel.Debug) + + pcd = meshes.eagle() + print(pcd) + o3d.visualization.draw_geometries([pcd]) + + print('run Poisson surface reconstruction') + mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( + pcd, depth=8) + print(mesh) + o3d.visualization.draw_geometries([mesh]) + + print('visualize densities') + densities = np.asarray(densities) + density_colors = plt.get_cmap('plasma')( + (densities - densities.min()) / (densities.max() - densities.min())) + density_colors = density_colors[:, :3] + density_mesh = o3d.geometry.TriangleMesh() + density_mesh.vertices = mesh.vertices + density_mesh.triangles = mesh.triangles + density_mesh.triangle_normals = mesh.triangle_normals + density_mesh.vertex_colors = o3d.utility.Vector3dVector(density_colors) + o3d.visualization.draw_geometries([density_mesh]) + + print('remove low density vertices') + vertices_to_remove = densities < np.quantile(densities, 0.1) + mesh.remove_vertices_by_mask(vertices_to_remove) + print(mesh) + o3d.visualization.draw_geometries([mesh]) diff --git a/examples/Python/Misc/meshes.py b/examples/Python/Misc/meshes.py index 55964440a5d..1bff416eca7 100644 --- a/examples/Python/Misc/meshes.py +++ b/examples/Python/Misc/meshes.py @@ -176,6 +176,16 @@ def bunny(): return mesh +def eagle(): + path = _relative_path("../../TestData/eagle.ply") + if not os.path.exists(path): + print("downloading eagle pcl") + url = "http://www.cs.jhu.edu/~misha/Code/PoissonRecon/eagle.points.ply" + urllib.request.urlretrieve(url, path) + pcd = o3d.io.read_point_cloud(path) + return pcd + + def center_and_scale(mesh): vertices = np.asarray(mesh.vertices) vertices = vertices / max(vertices.max(axis=0) - vertices.min(axis=0)) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1be08a2cd2e..761a6219973 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,10 +1,10 @@ function(ShowAndAbortOnWarning trgt) if(MSVC) - # target_compile_options(${trgt} PRIVATE /W4 /WX) - target_compile_options(${trgt} PRIVATE /W3) + # target_compile_options(${trgt} PRIVATE /W4 /WX) + target_compile_options(${trgt} PRIVATE /W3) else() - # target_compile_options(${trgt} PRIVATE -Wall -Wextra -Werror) - target_compile_options(${trgt} PRIVATE -Wall -Werror) + # target_compile_options(${trgt} PRIVATE -Wall -Wextra -Werror) + target_compile_options(${trgt} PRIVATE -Wall -Werror) endif() endfunction() diff --git a/src/Open3D/Geometry/SurfaceReconstructionPoisson.cpp b/src/Open3D/Geometry/SurfaceReconstructionPoisson.cpp new file mode 100644 index 00000000000..ac8263c012f --- /dev/null +++ b/src/Open3D/Geometry/SurfaceReconstructionPoisson.cpp @@ -0,0 +1,763 @@ +// ---------------------------------------------------------------------------- +// - Open3D: www.open3d.org - +// ---------------------------------------------------------------------------- +// The MIT License (MIT) +// +// Copyright (c) 2018 www.open3d.org +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. +// ---------------------------------------------------------------------------- + +#include "Open3D/Geometry/PointCloud.h" +#include "Open3D/Geometry/TriangleMesh.h" +#include "Open3D/Utility/Console.h" + +#include +#include +#include +#include +#include +#include +#include + +// clang-format off +#include "PoissonRecon/Src/PreProcessor.h" +#include "PoissonRecon/Src/MyMiscellany.h" +#include "PoissonRecon/Src/CmdLineParser.h" +#include "PoissonRecon/Src/FEMTree.h" +#include "PoissonRecon/Src/PPolynomial.h" +#include "PoissonRecon/Src/PointStreamData.h" +// clang-format on + +namespace open3d { +namespace geometry { +namespace poisson { + +// 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 default finite-element boundary type +static const BoundaryType DEFAULT_FEM_BOUNDARY = BOUNDARY_NEUMANN; +// The dimension of the system +static const int DIMENSION = 3; + +class Open3DData { +public: + Open3DData() : normal_(0, 0, 0), color_(0, 0, 0) {} + Open3DData(const Eigen::Vector3d& normal, const Eigen::Vector3d& color) + : normal_(normal), color_(color) {} + + Open3DData operator*(double s) const { + return Open3DData(s * normal_, s * color_); + } + Open3DData operator/(double s) const { + return Open3DData(normal_ / s, (1 / s) * color_); + } + Open3DData& operator+=(const Open3DData& d) { + normal_ += d.normal_; + color_ += d.color_; + return *this; + } + Open3DData& operator*=(double s) { + normal_ *= s; + color_ *= s; + return *this; + } + +public: + Eigen::Vector3d normal_; + Eigen::Vector3d color_; +}; + +template +class Open3DPointStream + : public InputPointStreamWithData { +public: + Open3DPointStream(const open3d::geometry::PointCloud* pcd) + : pcd_(pcd), xform_(nullptr), current_(0) {} + void reset(void) { current_ = 0; } + bool nextPoint(Point& p, Open3DData& d) { + if (current_ >= pcd_->points_.size()) { + return false; + } + p.coords[0] = pcd_->points_[current_](0); + p.coords[1] = pcd_->points_[current_](1); + p.coords[2] = pcd_->points_[current_](2); + + if (xform_ != nullptr) { + p = (*xform_) * p; + } + + if (pcd_->HasNormals()) { + d.normal_ = pcd_->normals_[current_]; + } else { + d.normal_ = Eigen::Vector3d(0, 0, 0); + } + + if (pcd_->HasColors()) { + d.color_ = pcd_->colors_[current_]; + } else { + d.color_ = Eigen::Vector3d(0, 0, 0); + } + + current_++; + return true; + } + +public: + const open3d::geometry::PointCloud* pcd_; + XForm* xform_; + size_t current_; +}; + +template +class Open3DVertex { +public: + typedef _Real Real; + + Open3DVertex() : Open3DVertex(Point(0, 0, 0)) {} + Open3DVertex(Point point) + : point(point), normal_(0, 0, 0), color_(0, 0, 0), w_(0) {} + + Open3DVertex& operator*=(Real s) { + point *= s; + normal_ *= s; + color_ *= s; + w_ *= s; + return *this; + } + + Open3DVertex& operator+=(const Open3DVertex& p) { + point += p.point; + normal_ += p.normal_; + color_ += p.color_; + w_ += p.w_; + return *this; + } + + Open3DVertex& operator/=(Real s) { + point /= s; + normal_ /= s; + color_ /= s; + w_ /= s; + return *this; + } + +public: + // point can not have trailing _, because template methods assume that it is + // named this way + Point point; + Eigen::Vector3d normal_; + Eigen::Vector3d color_; + double w_; +}; + +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 +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* sampleData, + const typename FEMTree::template DensityEstimator* + density, + const SetVertexFunction& SetVertex, + XForm iXForm, + std::shared_ptr& out_mesh, + std::vector& out_densities) { + static const int Dim = sizeof...(FEMSigs); + typedef UIntPack Sigs; + static const unsigned int DataSig = + FEMDegreeAndBType::Signature; + typedef typename FEMTree::template DensityEstimator + DensityEstimator; + + FEMTreeProfiler profiler(tree); + + CoredMeshData* mesh; + mesh = new CoredVectorMeshData(); + + bool non_manifold = true; + bool polygon_mesh = false; + + profiler.start(); + typename IsoSurfaceExtractor::IsoStats isoStats; + if (sampleData) { + SparseNodeData, + IsotropicUIntPack> + _sampleData = + tree.template setMultiDepthDataField( + *samples, *sampleData, (DensityEstimator*)NULL); + for (const RegularTreeNode* + n = tree.tree().nextNode(); + n; n = tree.tree().nextNode(n)) { + ProjectiveData* clr = _sampleData(n); + if (clr) (*clr) *= (Real)pow(datax, tree.depth(n)); + } + isoStats = IsoSurfaceExtractor::template Extract< + Open3DData>(Sigs(), UIntPack(), + UIntPack(), tree, density, &_sampleData, + solution, isoValue, *mesh, SetVertex, !linear_fit, + !non_manifold, polygon_mesh, false); + } else { + isoStats = IsoSurfaceExtractor::template Extract< + Open3DData>(Sigs(), UIntPack(), + UIntPack(), tree, density, NULL, solution, + isoValue, *mesh, SetVertex, !linear_fit, + !non_manifold, polygon_mesh, false); + } + + mesh->resetIterator(); + out_densities.clear(); + for (size_t vidx = 0; vidx < mesh->outOfCorePointCount(); ++vidx) { + Vertex v; + mesh->nextOutOfCorePoint(v); + v.point = iXForm * v.point; + out_mesh->vertices_.push_back( + Eigen::Vector3d(v.point[0], v.point[1], v.point[2])); + out_mesh->vertex_normals_.push_back(v.normal_); + out_mesh->vertex_colors_.push_back(v.color_); + out_densities.push_back(v.w_); + } + for (size_t tidx = 0; tidx < mesh->polygonCount(); ++tidx) { + std::vector> triangle; + mesh->nextPolygon(triangle); + if (triangle.size() != 3) { + open3d::utility::LogError("got polygon"); + } else { + out_mesh->triangles_.push_back(Eigen::Vector3i( + triangle[0].idx, triangle[1].idx, triangle[2].idx)); + } + } + + delete mesh; +} + +template +void Execute(const open3d::geometry::PointCloud& pcd, + std::shared_ptr& out_mesh, + std::vector& out_densities, + int depth, + size_t 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.f; + 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 sampleData; + DensityEstimator* density = NULL; + SparseNodeData, NormalSigs>* normalInfo = NULL; + Real targetValue = (Real)0.5; + + // Read in the samples (and color data) + { + Open3DPointStream pointStream(&pcd); + + if (width > 0) { + xForm = GetPointXForm(pointStream, width, + (Real)(scale > 0 ? scale : 1.), + depth) * + xForm; + } else { + xForm = scale > 0 ? GetPointXForm(pointStream, + (Real)scale) * + xForm + : xForm; + } + + pointStream.xform_ = &xForm; + + { + auto ProcessDataWithConfidence = [&](const Point& p, + Open3DData& d) { + Real l = (Real)d.normal_.norm(); + if (!l || l != l) return (Real)-1.; + return (Real)pow(l, confidence); + }; + auto ProcessData = [](const Point& p, Open3DData& d) { + Real l = (Real)d.normal_.norm(); + if (!l || l != l) return (Real)-1.; + d.normal_ /= l; + return (Real)1.; + }; + if (confidence > 0) { + pointCount = FEMTreeInitializer::template Initialize< + Open3DData>(tree.spaceRoot(), pointStream, depth, + samples, sampleData, true, + tree.nodeAllocators[0], tree.initializer(), + ProcessDataWithConfidence); + } else { + pointCount = FEMTreeInitializer::template Initialize< + Open3DData>(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&)> + ConversionFunction = + [](Open3DData 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&, Real&)> + ConversionAndBiasFunction = [&](Open3DData 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 = utility::Logger::i().verbosity_level_ == + utility::VerbosityLevel::Debug, + sInfo.showResidual = 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 = [](Open3DVertex& v, Point p, Real w, + Open3DData d) { + v.point = p; + v.normal_ = d.normal_; + v.color_ = d.color_; + v.w_ = w; + }; + ExtractMesh, Real>( + datax, linear_fit, UIntPack(), + std::tuple(), tree, solution, isoValue, &samples, + &sampleData, density, SetVertex, iXForm, out_mesh, out_densities); + + if (density) delete density, density = NULL; + utility::LogDebug("# Total Solve: {:9.1f} (s), {:9.1f} (MB)", + Time() - startTime, FEMTree::MaxMemoryUsage()); +} + +} // namespace poisson + +std::tuple, std::vector> +TriangleMesh::CreateFromPointCloudPoisson(const PointCloud& pcd, + size_t depth, + size_t width, + float scale, + bool linear_fit) { + static const BoundaryType BType = poisson::DEFAULT_FEM_BOUNDARY; + typedef IsotropicUIntPack< + poisson::DIMENSION, + FEMDegreeAndBType::Signature> + FEMSigs; + + if (!pcd.HasNormals()) { + utility::LogError("[CreateFromPointCloudPoisson] pcd has no normals"); + } + +#ifdef _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 + + auto mesh = std::make_shared(); + std::vector densities; + poisson::Execute(pcd, mesh, densities, depth, width, scale, + linear_fit, FEMSigs()); + + ThreadPool::Terminate(); + + return std::make_tuple(mesh, densities); +} + +} // namespace geometry +} // namespace open3d diff --git a/src/Open3D/Geometry/TriangleMesh.cpp b/src/Open3D/Geometry/TriangleMesh.cpp index 45f7db2485b..8fa8cbc7bca 100644 --- a/src/Open3D/Geometry/TriangleMesh.cpp +++ b/src/Open3D/Geometry/TriangleMesh.cpp @@ -1465,6 +1465,67 @@ void TriangleMesh::RemoveTrianglesByMask( } } +void TriangleMesh::RemoveVerticesByIndex( + const std::vector &vertex_indices) { + std::vector vertex_mask(vertices_.size(), false); + for (auto vidx : vertex_indices) { + if (vidx >= 0 && vidx < vertices_.size()) { + vertex_mask[vidx] = true; + } else { + utility::LogWarning( + "[RemoveVerticessByIndex] contains vertex index {} that is " + "not within the bounds", + vidx); + } + } + + RemoveVerticesByMask(vertex_mask); +} + +void TriangleMesh::RemoveVerticesByMask(const std::vector &vertex_mask) { + if (vertex_mask.size() != vertices_.size()) { + utility::LogError("vertex_mask has a different size than vertices_"); + } + + bool has_normal = HasVertexNormals(); + bool has_color = HasVertexColors(); + int to_vidx = 0; + std::unordered_map vertex_map; + for (size_t from_vidx = 0; from_vidx < vertices_.size(); ++from_vidx) { + if (!vertex_mask[from_vidx]) { + vertex_map[from_vidx] = to_vidx; + vertices_[to_vidx] = vertices_[from_vidx]; + if (has_normal) { + vertex_normals_[to_vidx] = vertex_normals_[from_vidx]; + } + if (has_color) { + vertex_colors_[to_vidx] = vertex_colors_[from_vidx]; + } + to_vidx++; + } + } + vertices_.resize(to_vidx); + if (has_normal) { + vertex_normals_.resize(to_vidx); + } + if (has_color) { + vertex_colors_.resize(to_vidx); + } + + std::vector triangle_mask(triangles_.size()); + for (size_t tidx = 0; tidx < triangles_.size(); ++tidx) { + auto &tria = triangles_[tidx]; + triangle_mask[tidx] = vertex_mask[tria(0)] || vertex_mask[tria(1)] || + vertex_mask[tria(2)]; + if (!triangle_mask[tidx]) { + tria(0) = vertex_map[tria(0)]; + tria(1) = vertex_map[tria(1)]; + tria(2) = vertex_map[tria(2)]; + } + } + RemoveTrianglesByMask(triangle_mask); +} + std::unordered_map> diff --git a/src/Open3D/Geometry/TriangleMesh.h b/src/Open3D/Geometry/TriangleMesh.h index 8364b50e1d7..5edce638f87 100644 --- a/src/Open3D/Geometry/TriangleMesh.h +++ b/src/Open3D/Geometry/TriangleMesh.h @@ -385,6 +385,22 @@ class TriangleMesh : public MeshBase { /// Should have same size as \ref triangles_. void RemoveTrianglesByMask(const std::vector &triangle_mask); + /// \brief This function removes the vertices with index in + /// \p vertex_indices. Note that also all triangles associated with the + /// vertices are removeds. + /// + /// \param triangle_indices Indices of the triangles that should be + /// removed. + void RemoveVerticesByIndex(const std::vector &vertex_indices); + + /// \brief This function removes the vertices that are masked in + /// \p vertex_mask. Note that also all triangles associated with the + /// vertices are removed.. + /// + /// \param vertex_mask Mask of vertices that should be removed. + /// Should have same size as \ref vertices_. + void RemoveVerticesByMask(const std::vector &vertex_mask); + /// \brief This function deforms the mesh using the method by /// Sorkine and Alexa, "As-Rigid-As-Possible Surface Modeling", 2007. /// @@ -429,6 +445,34 @@ class TriangleMesh : public MeshBase { static std::shared_ptr CreateFromPointCloudBallPivoting( const PointCloud &pcd, const std::vector &radii); + /// \brief Function that computes a triangle mesh from a oriented PointCloud + /// pcd. This implements the Screened Poisson Reconstruction proposed in + /// Kazhdan and Hoppe, "Screened Poisson Surface Reconstruction", 2013. + /// This function uses the original implementation by Kazhdan. See + /// https://github.com/mkazhdan/PoissonRecon + /// + /// \param pcd PointCloud with normals and optionally colors. + /// \param depth Maximum depth of the tree that will be used for surface + /// reconstruction. Running at depth d corresponds to solving on a 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. + /// \param width Specifies the + /// target width of the finest level octree cells. This parameter is ignored + /// if depth is specified. + /// \param scale Specifies the ratio between the + /// diameter of the cube used for reconstruction and the diameter of the + /// samples' bounding cube. \param linear_fit If true, the reconstructor use + /// linear interpolation to estimate the positions of iso-vertices. + /// \return The estimated TriangleMesh, and per vertex densitie values that + /// can be used to to trim the mesh. + static std::tuple, std::vector> + CreateFromPointCloudPoisson(const PointCloud &pcd, + size_t depth = 8, + size_t width = 0, + float scale = 1.1f, + bool linear_fit = false); + /// Factory function to create a tetrahedron mesh (trianglemeshfactory.cpp). /// the mesh centroid will be at (0,0,0) and \param radius defines the /// distance from the center to the mesh vertices. diff --git a/src/Python/open3d_pybind/geometry/trianglemesh.cpp b/src/Python/open3d_pybind/geometry/trianglemesh.cpp index 777e36233ce..bb14d9963f2 100644 --- a/src/Python/open3d_pybind/geometry/trianglemesh.cpp +++ b/src/Python/open3d_pybind/geometry/trianglemesh.cpp @@ -293,6 +293,18 @@ void pybind_trianglemesh(py::module &m) { "set to true. Call remove_unreferenced_vertices to clean up " "vertices afterwards.", "triangle_mask"_a) + .def("remove_vertices_by_index", + &geometry::TriangleMesh::RemoveVerticesByIndex, + "This function removes the vertices with index in " + "vertex_indices. Note that also all triangles associated with " + "the vertices are removed.", + "vertex_indices"_a) + .def("remove_vertices_by_mask", + &geometry::TriangleMesh::RemoveVerticesByMask, + "This function removes the vertices that are masked in " + "vertex_mask. Note that also all triangles associated with " + "the vertices are removed.", + "vertex_mask"_a) .def("deform_as_rigid_as_possible", &geometry::TriangleMesh::DeformAsRigidAsPossible, "This function deforms the mesh using the method by Sorkine " @@ -321,6 +333,16 @@ void pybind_trianglemesh(py::module &m) { "radius over the point cloud, whenever the ball touches " "three points a triangle is created.", "pcd"_a, "radii"_a) + .def_static("create_from_point_cloud_poisson", + &geometry::TriangleMesh::CreateFromPointCloudPoisson, + "Function that computes a triangle mesh from a " + "oriented PointCloud pcd. This implements the Screened " + "Poisson Reconstruction proposed in Kazhdan and Hoppe, " + "\"Screened Poisson Surface Reconstruction\", 2013. " + "This function uses the original implementation by " + "Kazhdan. See https://github.com/mkazhdan/PoissonRecon", + "pcd"_a, "depth"_a = 8, "width"_a = 0, "scale"_a = 1.1, + "linear_fit"_a = false) .def_static("create_box", &geometry::TriangleMesh::CreateBox, "Factory function to create a box. The left bottom " "corner on the " @@ -568,6 +590,16 @@ void pybind_trianglemesh(py::module &m) { {{"triangle_mask", "1D bool array, True values indicate " "triangles that should be removed."}}); + docstring::ClassMethodDocInject( + m, "TriangleMesh", "remove_vertices_by_index", + {{"vertex_indices", + "1D array of vertex indices that should be removed from the " + "TriangleMesh."}}); + docstring::ClassMethodDocInject(m, "TriangleMesh", + "remove_vertices_by_mask", + {{"vertex_mask", + "1D bool array, True values indicate " + "vertices that should be removed."}}); docstring::ClassMethodDocInject( m, "TriangleMesh", "deform_as_rigid_as_possible", {{"constraint_vertex_indices", @@ -599,6 +631,27 @@ void pybind_trianglemesh(py::module &m) { {"radii", "The radii of the ball that are used for the surface " "reconstruction."}}); + docstring::ClassMethodDocInject( + m, "TriangleMesh", "create_from_point_cloud_poisson", + {{"pcd", + "PointCloud from which the TriangleMesh surface is " + "reconstructed. Has to contain normals."}, + {"depth", + "Maximum depth of the tree that will be used for surface " + "reconstruction. Running at depth d corresponds to solving on a " + "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."}, + {"width", + "Specifies the target width of the finest level octree cells. " + "This parameter is ignored if depth is specified"}, + {"scale", + "Specifies the ratio between the diameter of the cube used for " + "reconstruction and the diameter of the samples' bounding cube."}, + {"linear_fit", + "If true, the reconstructor use linear interpolation to estimate " + "the positions of iso-vertices."}}); docstring::ClassMethodDocInject(m, "TriangleMesh", "create_box", {{"width", "x-directional length."}, {"height", "y-directional length."}, diff --git a/src/UnitTest/Geometry/TriangleMesh.cpp b/src/UnitTest/Geometry/TriangleMesh.cpp index 3164a005f91..681c2a2c130 100644 --- a/src/UnitTest/Geometry/TriangleMesh.cpp +++ b/src/UnitTest/Geometry/TriangleMesh.cpp @@ -35,12 +35,13 @@ using namespace std; using namespace unit_test; void ExpectEQ(const open3d::geometry::TriangleMesh& mesh0, - const open3d::geometry::TriangleMesh& mesh1) { - ExpectEQ(mesh0.vertices_, mesh1.vertices_); - ExpectEQ(mesh0.vertex_normals_, mesh1.vertex_normals_); - ExpectEQ(mesh0.vertex_colors_, mesh1.vertex_colors_); + const open3d::geometry::TriangleMesh& mesh1, + double threshold = 1e-6) { + ExpectEQ(mesh0.vertices_, mesh1.vertices_, threshold); + ExpectEQ(mesh0.vertex_normals_, mesh1.vertex_normals_, threshold); + ExpectEQ(mesh0.vertex_colors_, mesh1.vertex_colors_, threshold); ExpectEQ(mesh0.triangles_, mesh1.triangles_); - ExpectEQ(mesh0.triangle_normals_, mesh1.triangle_normals_); + ExpectEQ(mesh0.triangle_normals_, mesh1.triangle_normals_, threshold); } TEST(TriangleMesh, Constructor) { @@ -1578,6 +1579,240 @@ TEST(TriangleMesh, CropTriangleMesh) { ExpectEQ(ref_triangle_normals, output_tm->triangle_normals_); } +TEST(TriangleMesh, CreateFromPointCloudPoisson) { + geometry::PointCloud pcd; + pcd.points_ = { + {-0.215279, 0.121252, 0.965784}, {0.079266, 0.643799, 0.755848}, + {0.001534, -0.691225, 0.720568}, {0.663793, 0.567055, 0.478929}, + {-0.929397, -0.262081, 0.238929}, {0.741441, -0.628480, 0.209423}, + {-0.844085, 0.510828, 0.131410}, {0.001478, -0.999346, 0.006834}, + {-0.075616, 0.987936, -0.077669}, {-0.599516, -0.738522, -0.297292}, + {0.952159, -0.061734, -0.284187}, {0.616751, -0.543117, -0.562931}, + {-0.066283, 0.653032, -0.748836}, {0.408743, 0.133839, -0.900872}, + {-0.165481, -0.187750, -0.964982}}; + pcd.normals_ = { + {-0.227789, 0.134744, 0.961608}, {0.084291, 0.646929, 0.752716}, + {-0.002898, -0.694359, 0.717585}, {0.666542, 0.565546, 0.476952}, + {-0.930073, -0.260750, 0.237944}, {0.740792, -0.629608, 0.208559}, + {-0.843751, 0.511569, 0.130869}, {0.001436, -0.999352, 0.006806}, + {-0.076095, 0.987987, -0.077349}, {-0.598343, -0.739977, -0.296067}, + {0.952641, -0.060018, -0.283016}, {0.620242, -0.541600, -0.560605}, + {-0.071145, 0.656184, -0.745733}, {0.414599, 0.141599, -0.897109}, + {-0.172774, -0.204293, -0.960815}}; + pcd.colors_ = { + {0.380208, 0.759545, 0.171340}, {0.288485, 0.063605, 0.643738}, + {0.997344, 0.888908, 0.802631}, {0.249312, 0.621767, 0.718174}, + {0.953116, 0.032469, 0.828249}, {0.763606, 0.914161, 0.458050}, + {0.735258, 0.974643, 0.458129}, {0.985589, 0.649781, 0.064284}, + {0.714224, 0.413067, 0.800399}, {0.433070, 0.962528, 0.138826}, + {0.066043, 0.867413, 0.276809}, {0.321559, 0.662692, 0.011849}, + {0.353338, 0.784374, 0.899556}, {0.248052, 0.431480, 0.339511}, + {0.880225, 0.614243, 0.379607}}; + geometry::TriangleMesh mesh_gt; + mesh_gt.vertices_ = {{-0.535122, -0.678988, -0.546102}, + {-0.856971, -0.552207, -0.546102}, + {0.011381, -0.852840, -0.546102}, + {-1.081624, -0.478258, -0.546102}, + {0.557883, -0.571511, -0.546102}, + {0.576765, -0.552207, -0.546102}, + {0.739454, -0.005705, -0.546102}, + {-1.081624, 0.339532, -0.546102}, + {-0.720499, 0.540798, -0.546102}, + {-0.535122, 0.637316, -0.546102}, + {0.011381, 0.788968, -0.546102}, + {0.559266, 0.540798, -0.546102}, + {0.557883, 0.542438, -0.546102}, + {-0.535122, -0.552207, -0.761037}, + {0.011381, -0.552207, -0.838612}, + {-1.081624, -0.005705, -1.024482}, + {-0.535122, -0.005705, -1.038080}, + {0.011381, -0.005705, -1.001911}, + {0.557883, -0.552207, -0.565591}, + {0.557883, -0.005705, -0.755925}, + {-0.535122, 0.540798, -0.678056}, + {0.011381, 0.540798, -0.807635}, + {0.557883, 0.540798, -0.547422}, + {-0.535122, -0.868303, 0.000401}, + {-0.957409, -0.552207, 0.000401}, + {0.011381, -0.981478, 0.000401}, + {-1.081624, -0.413615, 0.000401}, + {0.557883, -0.735421, 0.000401}, + {0.756038, -0.552207, 0.000401}, + {1.045917, -0.005705, 0.000401}, + {-1.081624, 0.279246, 0.000401}, + {-0.856915, 0.540798, 0.000401}, + {-0.535122, 0.806682, 0.000401}, + {0.011381, 0.995774, 0.000401}, + {0.876659, 0.540798, 0.000401}, + {0.557883, 0.851917, 0.000401}, + {0.011381, -0.775497, 0.546903}, + {-0.494233, -0.552207, 0.546903}, + {-0.535122, -0.513868, 0.546903}, + {-0.768247, -0.005705, 0.546903}, + {0.557883, -0.601114, 0.546903}, + {0.621286, -0.552207, 0.546903}, + {0.891579, -0.005705, 0.546903}, + {-0.535122, 0.523186, 0.546903}, + {-0.516131, 0.540798, 0.546903}, + {0.011381, 0.784610, 0.546903}, + {0.613233, 0.540798, 0.546903}, + {0.557883, 0.585415, 0.546903}, + {-0.535122, -0.552207, 0.520166}, + {-1.081624, -0.005705, 0.145186}, + {-0.535122, 0.540798, 0.529053}, + {0.011381, -0.552207, 0.785752}, + {-0.535122, -0.005705, 0.773541}, + {0.011381, -0.005705, 1.045411}, + {0.557883, -0.552207, 0.612610}, + {0.557883, -0.005705, 0.880973}, + {0.011381, 0.540798, 0.802148}, + {0.557883, 0.540798, 0.593895}}; + mesh_gt.vertex_normals_ = {{-0.363328, -0.449218, -0.596538}, + {-0.363328, -0.449218, -0.596538}, + {-0.363328, -0.449218, -0.596538}, + {-0.363328, -0.449218, -0.596538}, + {0.748097, -0.287048, -0.401288}, + {0.748097, -0.287048, -0.401288}, + {0.624682, -0.539519, -0.560146}, + {-0.067839, 0.778945, -0.391656}, + {-0.067839, 0.778945, -0.391656}, + {-0.067839, 0.778945, -0.391656}, + {-0.071381, 0.660859, -0.743014}, + {0.375628, 0.123591, -0.806610}, + {0.375628, 0.123591, -0.806610}, + {-0.363328, -0.449218, -0.596538}, + {-0.363328, -0.449218, -0.596538}, + {-0.363328, -0.449218, -0.596538}, + {-0.363328, -0.449218, -0.596538}, + {-0.176322, -0.208808, -0.957378}, + {0.748097, -0.287048, -0.401288}, + {0.748097, -0.287048, -0.401288}, + {-0.067839, 0.778945, -0.391656}, + {-0.067839, 0.778945, -0.391656}, + {0.414932, 0.141672, -0.897776}, + {-0.596520, -0.737720, -0.301946}, + {-0.596520, -0.737720, -0.301946}, + {-0.363328, -0.449218, -0.596538}, + {-0.363328, -0.449218, -0.596538}, + {0.748097, -0.287048, -0.401288}, + {0.748097, -0.287048, -0.401288}, + {0.953452, -0.063973, -0.286170}, + {-0.067839, 0.778945, -0.391656}, + {-0.067839, 0.778945, -0.391656}, + {-0.067839, 0.778945, -0.391656}, + {-0.076419, 0.990524, -0.082928}, + {0.375628, 0.123591, -0.806610}, + {0.375628, 0.123591, -0.806610}, + {-0.005705, -0.991160, 0.013986}, + {-0.005705, -0.991160, 0.013986}, + {-0.919808, -0.270837, 0.240837}, + {-0.919808, -0.270837, 0.240837}, + {0.670787, -0.570145, 0.185455}, + {0.744380, -0.632658, 0.209540}, + {0.670787, -0.570145, 0.185455}, + {-0.842126, 0.510549, 0.137740}, + {-0.506191, 0.304790, 0.516282}, + {-0.506191, 0.304790, 0.516282}, + {0.358188, 0.573235, 0.581938}, + {0.358188, 0.573235, 0.581938}, + {-0.298994, -0.629163, 0.308635}, + {-0.919808, -0.270837, 0.240837}, + {-0.842126, 0.510549, 0.137740}, + {-0.009941, -0.693804, 0.708891}, + {-0.298994, -0.629163, 0.308635}, + {-0.298994, -0.629163, 0.308635}, + {0.670787, -0.570145, 0.185455}, + {0.670787, -0.570145, 0.185455}, + {-0.232912, 0.137861, 0.956863}, + {0.358188, 0.573235, 0.581938}}; + mesh_gt.vertex_colors_ = { + {0.651185, 0.780322, 0.270666}, {0.651185, 0.780322, 0.270666}, + {0.651185, 0.780322, 0.270666}, {0.651185, 0.780322, 0.270666}, + {0.213958, 0.758281, 0.162138}, {0.213958, 0.758281, 0.162138}, + {0.319808, 0.664247, 0.014295}, {0.535119, 0.601156, 0.828728}, + {0.535119, 0.601156, 0.828728}, {0.535119, 0.601156, 0.828728}, + {0.356296, 0.781393, 0.898403}, {0.280560, 0.453637, 0.352787}, + {0.280560, 0.453637, 0.352787}, {0.651185, 0.780322, 0.270666}, + {0.651185, 0.780322, 0.270666}, {0.651185, 0.780322, 0.270666}, + {0.651185, 0.780322, 0.270666}, {0.876498, 0.616945, 0.377835}, + {0.213958, 0.758281, 0.162138}, {0.213958, 0.758281, 0.162138}, + {0.535119, 0.601156, 0.828728}, {0.535119, 0.601156, 0.828728}, + {0.248333, 0.431672, 0.339626}, {0.436619, 0.959563, 0.140971}, + {0.436619, 0.959563, 0.140971}, {0.651185, 0.780322, 0.270666}, + {0.651185, 0.780322, 0.270666}, {0.213958, 0.758281, 0.162138}, + {0.213958, 0.758281, 0.162138}, {0.068450, 0.865637, 0.274943}, + {0.535119, 0.601156, 0.828728}, {0.535119, 0.601156, 0.828728}, + {0.535119, 0.601156, 0.828728}, {0.711309, 0.416128, 0.800860}, + {0.280560, 0.453637, 0.352787}, {0.280560, 0.453637, 0.352787}, + {0.985049, 0.646890, 0.076101}, {0.985049, 0.646890, 0.076101}, + {0.953348, 0.044255, 0.821904}, {0.953348, 0.044255, 0.821904}, + {0.742035, 0.885688, 0.458892}, {0.763420, 0.913915, 0.458058}, + {0.742035, 0.885688, 0.458892}, {0.732370, 0.972691, 0.455932}, + {0.557746, 0.854674, 0.323112}, {0.557746, 0.854674, 0.323112}, + {0.284898, 0.359292, 0.669062}, {0.284898, 0.359292, 0.669062}, + {0.962866, 0.528193, 0.561334}, {0.953348, 0.044255, 0.821904}, + {0.732370, 0.972691, 0.455932}, {0.996524, 0.880332, 0.796895}, + {0.962866, 0.528193, 0.561334}, {0.962866, 0.528193, 0.561334}, + {0.742035, 0.885688, 0.458892}, {0.742035, 0.885688, 0.458892}, + {0.383097, 0.761093, 0.173810}, {0.284898, 0.359292, 0.669062}}; + mesh_gt.triangles_ = { + {1, 13, 0}, {0, 14, 2}, {13, 14, 0}, {1, 15, 13}, + {15, 16, 13}, {1, 3, 15}, {14, 13, 16}, {16, 17, 14}, + {2, 18, 4}, {14, 18, 2}, {4, 18, 5}, {18, 14, 17}, + {17, 19, 18}, {18, 19, 6}, {6, 5, 18}, {7, 16, 15}, + {7, 8, 16}, {8, 20, 16}, {21, 16, 20}, {17, 16, 21}, + {9, 20, 8}, {10, 20, 9}, {21, 20, 10}, {22, 17, 21}, + {19, 17, 22}, {19, 22, 11}, {11, 6, 19}, {22, 21, 10}, + {10, 12, 22}, {11, 22, 12}, {24, 0, 23}, {1, 0, 24}, + {0, 2, 25}, {25, 23, 0}, {3, 1, 24}, {24, 26, 3}, + {2, 4, 27}, {27, 25, 2}, {27, 5, 28}, {4, 5, 27}, + {28, 6, 29}, {5, 6, 28}, {8, 7, 30}, {30, 31, 8}, + {32, 8, 31}, {9, 8, 32}, {33, 9, 32}, {10, 9, 33}, + {6, 11, 34}, {34, 29, 6}, {12, 10, 33}, {33, 35, 12}, + {34, 12, 35}, {11, 12, 34}, {24, 23, 48}, {36, 23, 25}, + {36, 37, 23}, {37, 48, 23}, {38, 24, 48}, {38, 39, 24}, + {39, 26, 24}, {39, 49, 26}, {38, 48, 37}, {25, 27, 40}, + {40, 36, 25}, {27, 28, 41}, {41, 40, 27}, {28, 29, 42}, + {42, 41, 28}, {39, 30, 49}, {39, 31, 30}, {39, 50, 31}, + {39, 43, 50}, {44, 50, 43}, {32, 31, 50}, {44, 32, 50}, + {44, 45, 32}, {45, 33, 32}, {42, 34, 46}, {29, 34, 42}, + {47, 33, 45}, {35, 33, 47}, {34, 35, 47}, {47, 46, 34}, + {37, 36, 51}, {39, 38, 52}, {53, 52, 51}, {52, 37, 51}, + {52, 38, 37}, {36, 40, 54}, {54, 51, 36}, {54, 40, 41}, + {51, 54, 55}, {55, 53, 51}, {55, 41, 42}, {54, 41, 55}, + {43, 39, 52}, {56, 52, 53}, {56, 44, 52}, {44, 43, 52}, + {45, 44, 56}, {56, 55, 57}, {53, 55, 56}, {55, 42, 46}, + {46, 57, 55}, {45, 56, 57}, {57, 47, 45}, {57, 46, 47}}; + std::vector densities_gt = { + 0.39865168929100037, 0.32580316066741943, 0.39778709411621094, + 0.2200755625963211, 0.428702175617218, 0.4288075268268585, + 0.44350555539131165, 0.24208465218544006, 0.3794262409210205, + 0.407772421836853, 0.41914406418800354, 0.4330996870994568, + 0.4330378770828247, 0.3689049780368805, 0.4012255072593689, + 0.07905912399291992, 0.30074167251586914, 0.3844510018825531, + 0.4286557137966156, 0.43353089690208435, 0.3970388174057007, + 0.41389259696006775, 0.4331146478652954, 0.39110293984413147, + 0.3399316370487213, 0.3989846110343933, 0.294955313205719, + 0.44293972849845886, 0.4377133548259735, 0.36585745215415955, + 0.31271663308143616, 0.3885934352874756, 0.4121553897857666, + 0.3849177360534668, 0.3899257779121399, 0.3932742774486542, + 0.42859214544296265, 0.44281646609306335, 0.4422561824321747, + 0.4249078631401062, 0.4210644066333771, 0.41715115308761597, + 0.38326939940452576, 0.43757864832878113, 0.4378965198993683, + 0.42082998156547546, 0.41933944821357727, 0.42220038175582886, + 0.439664751291275, 0.3250156342983246, 0.437633752822876, + 0.4227336645126343, 0.42734524607658386, 0.35958021879196167, + 0.41651853919029236, 0.3831002116203308, 0.41637951135635376, + 0.42157289385795593}; + + std::shared_ptr mesh_es; + std::vector densities_es; + std::tie(mesh_es, densities_es) = + geometry::TriangleMesh::CreateFromPointCloudPoisson(pcd, 2); + ExpectEQ(*mesh_es, mesh_gt, 1e-4); + ExpectEQ(densities_es, densities_gt, 1e-4); +} + TEST(TriangleMesh, CreateFromPointCloudAlphaShape) { geometry::PointCloud pcd; pcd.points_ = { diff --git a/src/UnitTest/TestUtility/UnitTest.cpp b/src/UnitTest/TestUtility/UnitTest.cpp index 4c2736ae572..89d8752159a 100644 --- a/src/UnitTest/TestUtility/UnitTest.cpp +++ b/src/UnitTest/TestUtility/UnitTest.cpp @@ -85,16 +85,19 @@ void unit_test::ExpectEQ(const vector& v0, const vector& v1) { // ---------------------------------------------------------------------------- void unit_test::ExpectEQ(const float* const v0, const float* const v1, - const size_t& size) { - for (size_t i = 0; i < size; i++) EXPECT_NEAR(v0[i], v1[i], THRESHOLD_1E_6); + const size_t& size, + float threshold) { + for (size_t i = 0; i < size; i++) EXPECT_NEAR(v0[i], v1[i], threshold); } // ---------------------------------------------------------------------------- // Test equality of two vectors of float. // ---------------------------------------------------------------------------- -void unit_test::ExpectEQ(const vector& v0, const vector& v1) { +void unit_test::ExpectEQ(const vector& v0, + const vector& v1, + float threshold) { EXPECT_EQ(v0.size(), v1.size()); - ExpectEQ(v0.data(), v1.data(), v0.size()); + ExpectEQ(v0.data(), v1.data(), v0.size(), threshold); } // ---------------------------------------------------------------------------- @@ -102,14 +105,17 @@ void unit_test::ExpectEQ(const vector& v0, const vector& v1) { // ---------------------------------------------------------------------------- void unit_test::ExpectEQ(const double* const v0, const double* const v1, - const size_t& size) { - for (size_t i = 0; i < size; i++) EXPECT_NEAR(v0[i], v1[i], THRESHOLD_1E_6); + const size_t& size, + double threshold) { + for (size_t i = 0; i < size; i++) EXPECT_NEAR(v0[i], v1[i], threshold); } // ---------------------------------------------------------------------------- // Test equality of two vectors of double. // ---------------------------------------------------------------------------- -void unit_test::ExpectEQ(const vector& v0, const vector& v1) { +void unit_test::ExpectEQ(const vector& v0, + const vector& v1, + double threshold) { EXPECT_EQ(v0.size(), v1.size()); - ExpectEQ(v0.data(), v1.data(), v0.size()); + ExpectEQ(v0.data(), v1.data(), v0.size(), threshold); } diff --git a/src/UnitTest/TestUtility/UnitTest.h b/src/UnitTest/TestUtility/UnitTest.h index f0b3f73d624..35084106ce0 100644 --- a/src/UnitTest/TestUtility/UnitTest.h +++ b/src/UnitTest/TestUtility/UnitTest.h @@ -135,18 +135,26 @@ void ExpectEQ(const int* const v0, const int* const v1, const size_t& size); void ExpectEQ(const std::vector& v0, const std::vector& v1); // Test equality of two arrays of float. -void ExpectEQ(const float* const v0, const float* const v1, const size_t& size); +void ExpectEQ(const float* const v0, + const float* const v1, + const size_t& size, + float threshold = THRESHOLD_1E_6); // Test equality of two vectors of float. -void ExpectEQ(const std::vector& v0, const std::vector& v1); +void ExpectEQ(const std::vector& v0, + const std::vector& v1, + float threshold = THRESHOLD_1E_6); // Test equality of two arrays of double. void ExpectEQ(const double* const v0, const double* const v1, - const size_t& size); + const size_t& size, + double threshold = THRESHOLD_1E_6); // Test equality of two vectors of double. -void ExpectEQ(const std::vector& v0, const std::vector& v1); +void ExpectEQ(const std::vector& v0, + const std::vector& v1, + double threshold = THRESHOLD_1E_6); // Reinterpret cast from uint8_t* to float*. template