adding triangulation and tetrahedralization functions

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.10)
 
 project(Bindings
 	DESCRIPTION
-		"Python libigl bindings"
+		"Python bindings"
 )
 
 set(CMAKE_CXX_STANDARD 17)
@@ -23,10 +23,20 @@ if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
 endif()
 
 
+# Include libigl
 option(LIBIGL_COPYLEFT_CGAL "Use CGAL"          ON)
+option(LIBIGL_COPYLEFT_TETGEN "Use Tetgen"      ON)
 option(LIBIGL_EMBREE "Build target igl::embree" ON)
 include(libigl)
 
+# Add triangulation library.
+include(CDT)
+add_subdirectory("${CMAKE_BINARY_DIR}/_deps/cdt-src/CDT/" CDT)
+
+# List of all libraries to link
+set(LIBRARIES_TO_LINK igl::core igl::embree CDT)
+set(COPYLEFT_LIBRARIES_TO_LINK igl::core igl_copyleft::cgal igl_copyleft::tetgen)
+
 add_library(cpytoolbox
 	STATIC
 	# SHARED
@@ -81,9 +91,10 @@ add_library(cpytoolbox
 	src/cpp/swept_volume/swept_volume.h
 	)
 
-# target_link_libraries(cpytoolbox igl::core igl_copyleft::cgal igl::embree)
-target_link_libraries(cpytoolbox igl::core igl::embree )
-target_link_libraries(cpytoolbox_copyleft igl::core igl_copyleft::cgal)
+target_link_libraries(cpytoolbox ${LIBRARIES_TO_LINK})
+set(LIBRARIES_TO_LINK cpytoolbox ${LIBRARIES_TO_LINK})
+target_link_libraries(cpytoolbox_copyleft ${COPYLEFT_LIBRARIES_TO_LINK})
+set(COPYLEFT_LIBRARIES_TO_LINK cpytoolbox_copyleft ${COPYLEFT_LIBRARIES_TO_LINK})
 
 # Otman's python bindings using pybind11
 add_subdirectory(./ext/pybind11/)
@@ -106,16 +117,18 @@ pybind11_add_module(gpytoolbox_bindings
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_read_ply.cpp"
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_write_ply.cpp"
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_curved_hessian_intrinsic.cpp"
+	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_triangulate.cpp"
 )
 
 pybind11_add_module(gpytoolbox_bindings_copyleft
-    "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/gpytoolbox_bindings_copyleft_core.cpp"
-    "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_swept_volume.cpp"
+  "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/gpytoolbox_bindings_copyleft_core.cpp"
+  "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_swept_volume.cpp"
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_booleans.cpp"
+	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_tetrahedralize.cpp"
 )
 
-target_link_libraries(gpytoolbox_bindings PUBLIC cpytoolbox igl::core)
-target_link_libraries(gpytoolbox_bindings_copyleft PUBLIC cpytoolbox_copyleft igl::core igl_copyleft::cgal)
+target_link_libraries(gpytoolbox_bindings PUBLIC ${LIBRARIES_TO_LINK})
+target_link_libraries(gpytoolbox_bindings_copyleft PUBLIC ${COPYLEFT_LIBRARIES_TO_LINK})
 
 target_include_directories(gpytoolbox_bindings PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/")
 if(${CMAKE_SYSTEM_NAME} MATCHES "Windows")

LICENSE.GPL

@@ -1,4 +1,4 @@
-Copyright (c) 2022 Silvia Sellán and Oded Stein
+Copyright (c) 2024 Silvia Sellán and Oded Stein
 
                    GNU GENERAL PUBLIC LICENSE
                        Version 3, 29 June 2007
@@ -673,4 +673,4 @@ into proprietary programs.  If your program is a subroutine library, you
 may consider it more useful to permit linking proprietary applications with
 the library.  If this is what you want to do, use the GNU Lesser General
 Public License instead of this License.  But first, please read
-<http://www.gnu.org/philosophy/why-not-lgpl.html>.
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.

LICENSE.MIT

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2022 Silvia Sellán and Oded Stein
+Copyright (c) 2024 Silvia Sellán and Oded Stein
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
@@ -28,3 +28,12 @@ This software contains the ColorBrewer color maps by Cynthia Brewer
 (Pennsylvania State University), released under the Apache 2.0 license
 (http://www.personal.psu.edu/cab38/ColorBrewer/ColorBrewer_updates.html),
 and modifications by AxisMaps (https://colorbrewer2.org)
+
+========================
+
+This software includes code from libigl, released under MPL2 (https://github.com/libigl/libigl/blob/main/LICENSE.MPL2)
+
+========================
+
+This software includes code from CDT, released under MPL2 (https://github.com/artem-ogre/CDT/blob/master/LICENSE), and an active pull request by Islam0mar (https://github.com/Islam0mar/CDT/tree/ruppert-delaunay-refinement)
+

README.md

@@ -105,8 +105,8 @@ not you individually).
 ## License
 
 Gpytoolbox's is released under an MIT license ([see details](/LICENSE.MIT)),
-except for files in the `gpytoolbox.copyleft` module, which are under a GPL one
-([see details](/LICENSE.GPL)). Functions in the copyleft module must be imported
+except for files in the `gpytoolbox.copyleft` module, which are under a GPL/AGPL one
+(see [here](/LICENSE.GPL) and [here](/LICENSE.AGPL)). Functions in the copyleft module must be imported
 explicitly; this way, if you import only the main Gpytoolbox module
 ```python
 import gpytoolbox
@@ -120,7 +120,7 @@ import any functionality from `gpytoolbox.copyleft`; e.g.,
 ```python
 from gpytoolbox.copyleft import mesh_boolean
 ```
-you will be bound by the more restrictive GPL license.
+you will be bound by the more restrictive GPL/AGPL license.
 
 ## Attribution
 

cmake/CDT.cmake

@@ -0,0 +1,11 @@
+if(TARGET CDT)
+    return()
+endif()
+
+include(FetchContent)
+FetchContent_Declare(
+    CDT
+    GIT_REPOSITORY https://github.com/artem-ogre/CDT.git
+    GIT_TAG 24ff7d7969dded0243c9e1e992a5398a6d0293dd
+)
+FetchContent_MakeAvailable(CDT)

src/cpp/binding_tetrahedralize.cpp

@@ -0,0 +1,43 @@
+#include <pybind11/stl.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/eigen.h>
+#include <pybind11/functional.h>
+#include <string>
+#include <igl/copyleft/tetgen/tetrahedralize.h>
+
+using namespace Eigen;
+namespace py = pybind11;
+using EigenDStride = Stride<Eigen::Dynamic, Eigen::Dynamic>;
+template <typename MatrixType>
+using EigenDRef = Ref<MatrixType, 0, EigenDStride>; //allows passing column/row order matrices easily
+
+void binding_tetrahedralize(py::module& m) {
+    m.def("_tetrahedralize_cpp_impl",[](EigenDRef<MatrixXd> _V,
+                         EigenDRef<MatrixXi> _F,
+                         EigenDRef<MatrixXd> _H,
+                         double max_volume,
+                         double min_rad_edge_ratio)
+        {
+
+            std::ostringstream params_stream;
+            params_stream << "Q";
+            if(max_volume>0) {
+                params_stream << "a" << max_volume;
+            }
+            if(min_rad_edge_ratio>0) {
+                params_stream << "q" << min_rad_edge_ratio;
+            }
+            Eigen::MatrixXd V(_V), H(_H);
+            Eigen::MatrixXi F(_F);
+            Eigen::MatrixXd R,W;
+            Eigen::MatrixXi TR,PT;
+            Eigen::MatrixXi T,TF,FT,TN;
+            size_t numRegions;
+            int status = igl::copyleft::tetgen::tetrahedralize(V,F,H,
+                R,
+                params_stream.str(),
+                W,T,TF,TR,TN,PT,FT,numRegions);
+            return std::make_tuple(status, W, T, TF);
+        });
+
+}

src/cpp/binding_triangulate.cpp

@@ -0,0 +1,55 @@
+#include <CDT.h>
+#include <pybind11/stl.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/eigen.h>
+#include <pybind11/functional.h>
+#include <string>
+
+using namespace Eigen;
+namespace py = pybind11;
+using EigenDStride = Stride<Eigen::Dynamic, Eigen::Dynamic>;
+template <typename MatrixType>
+using EigenDRef = Ref<MatrixType, 0, EigenDStride>; //allows passing column/row order matrices easily
+
+void binding_triangulate(py::module& m) {
+    m.def("_triangulate_cpp_impl",[](EigenDRef<MatrixXd> V,
+                         EigenDRef<MatrixXi> E,
+                         double max_area,
+                         double min_angle,
+                         int max_steiner_points)
+        {
+
+            std::vector<CDT::V2d<double> > vertices;
+            for(int i=0; i<V.rows(); ++i) {
+                vertices.push_back(CDT::V2d<double>::make(V(i,0), V(i,1)));
+            }
+            std::vector<CDT::Edge> edges;
+            if(E.size()>0) {
+                for(int i=0; i<E.rows(); ++i) {
+                    edges.emplace_back(E(i,0), E(i,1));
+                }
+            }
+            CDT::RemoveDuplicatesAndRemapEdges(vertices,edges);
+            CDT::Triangulation<double> cdt;
+            cdt.insertVertices(vertices);
+            if(E.size()>0) {
+                cdt.conformToEdges(edges);
+                cdt.eraseOuterTrianglesAndHoles();
+            } else {
+                cdt.eraseSuperTriangle();
+            }
+            Eigen::MatrixXd W(cdt.vertices.size(), 2);
+            for(int i=0; i<cdt.vertices.size(); ++i) {
+                W.row(i) << cdt.vertices[i].x, cdt.vertices[i].y;
+            }
+            Eigen::MatrixXi F(cdt.triangles.size(), 3);
+            for(int i=0; i<cdt.triangles.size(); ++i) {
+                for(int j=0; j<3; ++j) {
+                    F(i,j) = cdt.triangles[i].vertices[j];
+                }
+            }
+
+            return std::make_tuple(W, F);
+        });
+
+}

src/cpp/gpytoolbox_bindings_copyleft_core.cpp

@@ -12,13 +12,15 @@ namespace py = pybind11;
 //forward declare all bindings
 void binding_swept_volume(py::module& m);
 void binding_booleans(py::module& m);
+void binding_tetrahedralize(py::module& m);
 
 
 PYBIND11_MODULE(gpytoolbox_bindings_copyleft, m) {
 
     /// call all bindings declared above  
     binding_swept_volume(m);
     binding_booleans(m);
+    binding_tetrahedralize(m);
 
     m.def("help", [&]() {printf("hi"); });
 }

src/cpp/gpytoolbox_bindings_core.cpp

@@ -27,6 +27,7 @@ void binding_upper_envelope(py::module& m);
 void binding_read_ply(py::module& m);
 void binding_write_ply(py::module& m);
 void binding_curved_hessian_intrinsic(py::module& m);
+void binding_triangulate(py::module& m);
 
 PYBIND11_MODULE(gpytoolbox_bindings, m) {
 
@@ -48,6 +49,7 @@ PYBIND11_MODULE(gpytoolbox_bindings, m) {
     binding_read_ply(m);
     binding_write_ply(m);
     binding_curved_hessian_intrinsic(m);
+    binding_triangulate(m);
 
     m.def("help", [&]() {printf("hi"); });
 }

src/gpytoolbox/__init__.py

@@ -117,4 +117,5 @@
 from .adjacency_matrix import adjacency_matrix
 from .non_manifold_edges import non_manifold_edges
 from .connected_components import connected_components
-from .rotation_matrix_from_vectors import rotation_matrix_from_vectors
+from .rotation_matrix_from_vectors import rotation_matrix_from_vectors
+from .triangulate import triangulate

src/gpytoolbox/copyleft/__init__.py

@@ -1,4 +1,5 @@
 from .mesh_boolean import mesh_boolean
 from .lazy_cage import lazy_cage
 from .do_meshes_intersect import do_meshes_intersect
-from .swept_volume import swept_volume
+from .swept_volume import swept_volume
+from .tetrahedralize import tetrahedralize

src/gpytoolbox/copyleft/tetrahedralize.py

@@ -0,0 +1,80 @@
+import numpy as np
+
+def tetrahedralize(V,F=None,
+    H=None,
+    max_volume=None,
+    min_rad_edge_ratio=None):
+    """Create a tetrahedralization of the input mesh in 3D.
+
+    This is a wrapper for libigl's Tetgen implementation.
+
+    Parameters
+    ----------
+    V : (n,3) numpy array
+        vertex list of points to mesh
+    F : (m,3) numpy int array, optional (default None)
+        face index list of a triangle.
+        If this is None, will simply mesh the convex hull of V.
+    H : (h,3) numpy array
+        list of seed points inside holes
+    max_volume: float, optional (default None)
+        If this is not None, the method will refine until this is the maximal
+        volume of any tet.
+    min_rad_edge_ratio : float, optional (default None)
+        If this is not None, the method will refine until this is the minimal
+        radius-edge ratio (in radians) of any tet. 
+        If this value is too aggressive, the method might silently fail to respect
+        the threshold, or hang forever.
+
+    Returns
+    -------
+    W : numpy double array
+        Matrix of mesh vertices
+    T : numpy int array
+        Matrix of tet indices
+    TF : numpy int array
+        Matrix of tet face indices
+
+    Notes
+    -----
+    In the future, we will hopefully use
+    "CDT - Constrained Delaunay Tetrahedrization made robust and practical" by
+    Diazzi et al. 2023.
+
+
+    Examples
+    --------
+    ```python
+    import gpytoolbox as gpy
+    # Mesh in V,F
+    W,T,TF = gpy.copyleft.tetrahedralize(V,F)
+    # Tet mesh in W,T,TF
+    ```
+    """
+
+    # Try to import C++ binding
+    try:
+        from gpytoolbox_bindings_copyleft import _tetrahedralize_cpp_impl
+    except:
+        raise ImportError("Gpytoolbox cannot import its C++ tetrahedralize binding.")
+
+    assert len(V.shape)==2 and V.shape[1]==3, "V must be a list of points in 3D."
+    assert max_volume is None or max_volume>0, "max_volume must be either None or a positive number."
+    if max_volume is None:
+        max_volume = -1.
+    assert min_rad_edge_ratio is None or min_rad_edge_ratio>=0, "max_angle must be either None or a nonnegative number."
+    if min_rad_edge_ratio is None:
+        min_rad_edge_ratio = -1.
+
+    if F is None:
+        F = np.array([], dtype=int)
+    if H is None:
+        H = np.array([], dtype=int)
+
+    status,W,T,TF = _tetrahedralize_cpp_impl(V.astype(np.float64),F.astype(np.int32),
+        H.astype(np.float64),
+        max_volume, min_rad_edge_ratio)
+    if status != 0:
+        raise RuntimeError("Tetgen failed.")
+
+    return W,T,TF