Poisson/Delaunay meshing (#93)

* Added poisson/delaunay meshing

* fixed poisson meshing

* Updated README.

---------

Co-authored-by: Mihai Dusmanu <mihai.dusmanu@gmail.com>
Co-authored-by: Paul-Edouard Sarlin <Skydes@users.noreply.github.com>

CMakeLists.txt

@@ -22,6 +22,14 @@ if (${COLMAP_CUDA_ENABLED})
 else()
     message("Compiling bindings without CUDA support.")
 endif()
+
+if (${COLMAP_CGAL_ENABLED})
+    message("Compiling bindings with CGAL support.")
+    add_definitions(-DCGAL_ENABLED)
+else()
+    message("Compiling bindings without CGAL support.")
+endif()
+
 pybind11_add_module(pycolmap main.cc)
 target_link_libraries(pycolmap PRIVATE ${COLMAP_LIBRARIES})
 if(MSVC)

README.md

@@ -70,7 +70,10 @@ pycolmap.patch_match_stereo(mvs_path)  # requires compilation with CUDA
 pycolmap.stereo_fusion(mvs_path / "dense.ply", mvs_path)
 ```
 
-PyCOLMAP can leverage the GPU for feature extraction, matching, and multi-view stereo if COLMAP was compiled with CUDA support. This requires to build the package from source and is not available with the PyPI wheels.
+PyCOLMAP can leverage the GPU for feature extraction, matching, and multi-view stereo if COLMAP was compiled with CUDA support.
+
+Similarly, PyCOLMAP can run Delauney Triangulation if COLMAP was compiled with CGAL support.
+This requires to build the package from source and is not available with the PyPI wheels.
 
 All of the above steps are easily configurable with python dicts which are recursively merged into
 their respective defaults, e.g.
@@ -312,6 +315,7 @@ camera_dict = {
 }
 ```
 
+
 ## SIFT feature extraction
 
 ```python

main.cc

@@ -21,6 +21,7 @@ namespace py = pybind11;
 
 #include "pipeline/sfm.cc"
 #include "pipeline/mvs.cc"
+#include "pipeline/meshing.cc"
 
 #include "reconstruction/reconstruction.cc"
 #include "reconstruction/correspondence_graph.cc"
@@ -98,6 +99,7 @@ PYBIND11_MODULE(pycolmap, m) {
     // Main reconstruction steps
     init_sfm(m);
     init_mvs(m);
+    init_meshing(m);
 
     // SIFT feature detector and descriptor
     init_sift(m);

pipeline/meshing.cc

@@ -0,0 +1,139 @@
+// Author: Philipp Lindenberger (Phil26AT)
+#include "colmap/mvs/meshing.h"
+
+using namespace colmap;
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+namespace py = pybind11;
+using namespace pybind11::literals;
+
+#include <string>
+
+#include "helpers.h"
+#include "log_exceptions.h"
+
+void init_meshing(py::module& m) {
+    using PoissonMOpts = mvs::PoissonMeshingOptions;
+    auto PyPoissonMeshingOptions =
+        py::class_<PoissonMOpts>(m, "PoissonMeshingOptions")
+            .def(py::init<>())
+            .def_readwrite("point_weight", &PoissonMOpts::point_weight,
+                "This floating point value specifies 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.")
+            .def_readwrite("depth", &PoissonMOpts::depth,
+                "This integer is the maximum depth of the tree that will be used for surface"
+                "reconstruction. Running at depth d corresponds to solving on a voxel grid"
+                "whose resolution is no larger than 2^d x 2^d x 2^d. Note that since the"
+                "reconstructor adapts the octree to the sampling density, the specified"
+                "reconstruction depth is only an upper bound.")
+            .def_readwrite("color", &PoissonMOpts::color,
+                "If specified, the reconstruction code assumes that the input is equipped"
+                "with colors and will extrapolate the color values to the vertices of the"
+                "reconstructed mesh. The floating point value specifies the relative"
+                "importance of finer color estimates over lower ones.")
+            .def_readwrite("trim", &PoissonMOpts::trim,
+                "This floating point values specifies the value for mesh trimming. The"
+                "subset of the mesh with signal value less than the trim value is discarded.")
+            .def_readwrite("num_threads", &PoissonMOpts::num_threads,
+                "The number of threads used for the Poisson reconstruction.");
+    make_dataclass(PyPoissonMeshingOptions);
+    auto poisson_options = PyPoissonMeshingOptions().cast<PoissonMOpts>();
+
+    using DMOpts = mvs::DelaunayMeshingOptions;
+    auto PyDelaunayMeshingOptions = 
+        py::class_<DMOpts>(m, "DelaunayMeshingOptions")
+            .def(py::init<>())
+            .def_readwrite("max_proj_dist", &DMOpts::max_proj_dist,
+                "Unify input points into one cell in the Delaunay triangulation that fall"
+                "within a reprojected radius of the given pixels.")
+            .def_readwrite("max_depth_dist", &DMOpts::max_depth_dist,
+                "Maximum relative depth difference between input point and a vertex of an"
+                "existing cell in the Delaunay triangulation, otherwise a new vertex is"
+                "created in the triangulation.")
+            .def_readwrite("visibility_sigma", &DMOpts::visibility_sigma,
+                "The standard deviation of wrt. the number of images seen by each point."
+                "Increasing this value decreases the influence of points seen in few images.")
+            .def_readwrite("distance_sigma_factor",
+                           &DMOpts::distance_sigma_factor,
+                "The factor that is applied to the computed distance sigma, which is"
+                "automatically computed as the 25th percentile of edge lengths. A higher"
+                "value will increase the smoothness of the surface.")
+            .def_readwrite("quality_regularization",
+                           &DMOpts::quality_regularization,
+                "A higher quality regularization leads to a smoother surface.")
+            .def_readwrite("max_side_length_factor",
+                           &DMOpts::max_side_length_factor,
+                "Filtering thresholds for outlier surface mesh faces. If the longest side of"
+                "a mesh face (longest out of 3) exceeds the side lengths of all faces at a"
+                "certain percentile by the given factor, then it is considered an outlier"
+                "mesh face and discarded.")
+            .def_readwrite("max_side_length_percentile",
+                           &DMOpts::max_side_length_percentile,
+                "Filtering thresholds for outlier surface mesh faces. If the longest side of"
+                "a mesh face (longest out of 3) exceeds the side lengths of all faces at a"
+                "certain percentile by the given factor, then it is considered an outlier"
+                "mesh face and discarded.")
+            .def_readwrite("num_threads", &DMOpts::num_threads,
+                "The number of threads to use for reconstruction. Default is all threads.");
+    make_dataclass(PyDelaunayMeshingOptions);
+    auto delaunay_options = PyDelaunayMeshingOptions().cast<DMOpts>();
+
+    m.def("poisson_meshing",
+          [](py::object input_path_,
+             py::object output_path_,
+             PoissonMOpts options) -> void {
+            std::string input_path = py::str(input_path_).cast<std::string>();
+            THROW_CHECK_HAS_FILE_EXTENSION(input_path, ".ply")
+            THROW_CHECK_FILE_EXISTS(input_path);
+
+            std::string output_path = py::str(output_path_).cast<std::string>();
+            THROW_CHECK_HAS_FILE_EXTENSION(output_path, ".ply")
+            THROW_CHECK_FILE_OPEN(output_path);
+            PoissonMeshing(options, input_path, output_path);
+          },
+          py::arg("input_path"),
+          py::arg("output_path"),
+          py::arg("options") = poisson_options,
+          "Perform Poisson surface reconstruction and return true if successful.");
+
+#ifdef CGAL_ENABLED
+    m.def("sparse_delaunay_meshing",
+          [](py::object input_path_,
+             py::object output_path_,
+             DMOpts options) -> void {
+            std::string input_path = py::str(input_path_).cast<std::string>();
+            THROW_CHECK_DIR_EXISTS(input_path);
+
+            std::string output_path = py::str(output_path_).cast<std::string>();
+            THROW_CHECK_HAS_FILE_EXTENSION(output_path, ".ply")
+            THROW_CHECK_FILE_OPEN(output_path);
+            mvs::SparseDelaunayMeshing(options, input_path, output_path);
+          },
+          py::arg("input_path"),
+          py::arg("output_path"),
+          py::arg("options") = delaunay_options,
+          "Delaunay meshing of sparse COLMAP reconstructions.");
+
+    m.def("dense_delaunay_meshing",
+          [](py::object input_path_,
+             py::object output_path_,
+             DMOpts options) -> void {
+            std::string input_path = py::str(input_path_).cast<std::string>();
+            THROW_CHECK_DIR_EXISTS(input_path);
+
+            std::string output_path = py::str(output_path_).cast<std::string>();
+            THROW_CHECK_HAS_FILE_EXTENSION(output_path, ".bin")
+            THROW_CHECK_FILE_OPEN(output_path);
+            mvs::DenseDelaunayMeshing(options, input_path, output_path);
+          },
+          py::arg("input_path"),
+          py::arg("output_path"),
+          py::arg("options") = delaunay_options,
+          "Delaunay meshing of dense COLMAP reconstructions.");
+#endif
+};