.gitlab-ci.yml

@@ -11,7 +11,7 @@ tests:
   stage: test 
   tags:
     - ci.inria.fr
-    - small
+    - medium
   script:
     - pip install pytest coverage
     - pip install -e .
@@ -33,7 +33,7 @@ build:
   stage: test 
   tags:
     - ci.inria.fr
-    - small
+    - medium
   script:
     - pip install --upgrade build
     - python -m build

CMakeLists.txt

@@ -9,7 +9,7 @@ project(
   VERSION ${SKBUILD_PROJECT_VERSION}
   LANGUAGES CXX)
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 # Find the module development requirements (requires FindPython from 3.17 or

cpp/CMakeLists.txt

@@ -2,8 +2,10 @@ project(_cpp)
 
 add_library(_cpp SHARED
     CgalLabelledImage.h
+    CgalImageMeshing.h
     LabelledTetrahedralMesh.h
     CgalLabelledImage.cpp
+    CgalImageMeshing.cpp
     LabelledTetrahedralMesh.cpp
 )
 target_link_libraries(_cpp PUBLIC Eigen3::Eigen)

cpp/CgalImageMeshing.cpp

+#include <iostream>
+#include <map>
+#include <chrono>
+
+#include "CgalLabelledImage.h"
+#include "LabelledTetrahedralMesh.h"
+
+#include "CgalImageMeshing.h"
+
+
+C3t3 meshCgalImage(CGAL::Image_3 image, double facet_angle, double facet_size, double facet_distance, double cell_radius_edge_ratio, double cell_size)
+{
+    auto start_time = std::chrono::system_clock::now();
+    std::cout<<"  --> Computing tetrahedrization"<<std::endl;
+    Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image);
+
+    // Mesh criteria
+    Mesh_criteria::Facet_criteria facet_criteria(facet_angle, facet_size, facet_distance); // angle, size, approximation
+    Mesh_criteria::Cell_criteria cell_criteria(cell_radius_edge_ratio, cell_size); // radius-edge ratio, size
+    Mesh_criteria criteria(facet_criteria, cell_criteria);
+
+    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+
+    auto end_time = std::chrono::system_clock::now();
+    std::cout<<"  <-- Computing tetrahedrization ("<<c3t3.number_of_cells()<<" cells, "<<c3t3.triangulation().number_of_vertices()<<" vertices)  ["<<(end_time - start_time).count()/1000.<<"s]"<<std::endl;
+
+    return c3t3;
+}
+
+LabelledTetrahedralMesh tetraMeshFromCgal(C3t3 c3t3)
+{
+    LabelledTetrahedralMesh mesh;
+
+    Tr tri = c3t3.triangulation();
+
+    auto start_time = std::chrono::system_clock::now();
+    std::cout<<"  --> Adding points"<<std::endl;
+    std::vector<Tr::Point> points;
+    std::unordered_map<double, std::unordered_map<double, std::unordered_map<double, int> > > point_vertices;
+    for (auto vertex_it = tri.all_vertices_begin(); vertex_it != tri.all_vertices_end(); ++vertex_it) {
+        auto point = vertex_it->point();
+        if (!point_vertices.contains(point.x())) {
+            point_vertices[point.x()] = std::unordered_map<double, std::unordered_map<double, int> >();
+        }
+        if (!point_vertices[point.x()].contains(point.y())) {
+            point_vertices[point.x()][point.y()] = std::unordered_map<double, int>();
+        }
+        point_vertices[point.x()][point.y()][point.z()] = points.size();
+        points.push_back(point);
+        mesh.addPoint({point.x(), point.y(), point.z()});
+    }
+    auto end_time = std::chrono::system_clock::now();
+    std::cout<<"  <-- Adding points ["<<(end_time - start_time).count()/1000.<<"s]"<<std::endl;
+
+    std::vector<std::vector<int> > triangle_vertices;
+    std::unordered_map<int, std::unordered_map<int, std::unordered_map<int, int> > > triangle_indices;
+    for (auto facet_it = c3t3.facets_in_complex_begin(); facet_it != c3t3.facets_in_complex_end(); ++facet_it) {
+        std::pair<Tr::Cell_handle, int> facet = *facet_it;
+        auto cell = facet.first;
+        auto f_id = facet.second;
+        auto triangle = tri.triangle(cell, f_id);
+        //std::cout<<triangle[0]<<" "<<triangle[1]<<" "<<triangle[2]<<std::endl;
+        std::vector<int> vertices;
+        for (int v=0; v<3; v++) {
+            auto point = triangle[v];
+            int vertex = point_vertices[point.x()][point.y()][point.z()];
+            vertices.push_back(vertex);
+        }
+        triangle_indices[vertices[0]][vertices[1]][vertices[2]] = triangle_vertices.size();
+        triangle_vertices.push_back(vertices);
+    }
+    // std::cout<<triangle_vertices.size()<<" triangles"<<std::endl;
+
+    std::vector<std::vector<int> > tetra_vertices;
+    std::vector<int> tetra_labels;
+    start_time = std::chrono::system_clock::now();
+    std::cout<<"  --> Adding tetras"<<std::endl;
+    for (auto cell_it = c3t3.cells_in_complex_begin(); cell_it != c3t3.cells_in_complex_end(); ++cell_it) {
+        auto tetra = tri.tetrahedron(cell_it);
+        std::vector<int> vertices;
+        for (int v=0; v<4; v++) {
+            auto point = tetra[v];
+            int vertex = point_vertices[point.x()][point.y()][point.z()];
+            vertices.push_back(vertex);
+        }
+        tetra_vertices.push_back(vertices);
+        int label = cell_it->subdomain_index();
+        tetra_labels.push_back(label);
+
+        mesh.addTetra({vertices[0], vertices[1], vertices[2], vertices[3]}, label);
+
+        // std::cout<<"  Tetra "<<vertices[0]<<" "<<vertices[1]<<" "<<vertices[2]<<" "<<vertices[3]<<std::endl;
+        for (int t=0; t<4; t++) {
+            auto triangle = tri.triangle(cell_it, t);
+            std::vector<int> t_vertices;
+            for (int v=0; v<3; v++) {
+                auto point = triangle[v];
+                int vertex = point_vertices[point.x()][point.y()][point.z()];
+                t_vertices.push_back(vertex);
+            }
+            // std::cout<<"  Triangle "<<triangle_indices[t_vertices[0]][t_vertices[1]][t_vertices[2]]<<std::endl;
+        }
+    }
+    // std::cout<<tetra_vertices.size()<<" tetras"<<std::endl;
+
+    end_time = std::chrono::system_clock::now();
+    std::cout<<"  <-- Adding tetras ["<<(end_time - start_time).count()/1000.<<"s]"<<std::endl;
+
+    return mesh;
+}
+
+LabelledTetrahedralMesh cgalImageMeshing(Eigen::Array<unsigned short, Eigen::Dynamic, 1> image_array, Eigen::Array<int, 3, 1> shape, Eigen::Array<double, 3, 1> voxelsize, double facet_angle, double facet_size, double facet_distance, double cell_radius_edge_ratio, double cell_size)
+{
+    CGAL::Image_3 image = buildCgalImage(image_array, shape, voxelsize);
+    C3t3 c3t3 = meshCgalImage(image, facet_angle, facet_size, facet_distance, cell_radius_edge_ratio, cell_size);
+    LabelledTetrahedralMesh mesh = tetraMeshFromCgal(c3t3);
+
+    return mesh;
+}

cpp/CgalImageMeshing.h

+#pragma once
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Labeled_image_mesh_domain_3.h>
+#include <CGAL/make_mesh_3.h>
+
+#include <CGAL/Image_3.h>
+
+#include <Eigen/Dense>
+
+class LabelledTetrahedralMesh;
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, CGAL::Sequential_tag>::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+
+C3t3 meshCgalImage(CGAL::Image_3 image, double facet_angle=30, double facet_size=6, double facet_distance=4, double cell_radius_edge_ratio=3, double cell_size=8);
+
+LabelledTetrahedralMesh tetraMeshFromCgal(C3t3 c3t3);
+
+LabelledTetrahedralMesh cgalImageMeshing(
+    Eigen::Array<unsigned short, Eigen::Dynamic, 1> image_array,
+    Eigen::Array<int, 3, 1> shape,
+    Eigen::Array<double, 3, 1> voxelsize={1.0, 1.0, 1.0},
+    double facet_angle=30,
+    double facet_size=6,
+    double facet_distance=4,
+    double cell_radius_edge_ratio=3,
+    double cell_size=8
+);

cpp/LabelledTetrahedralMesh.cpp

 #include "LabelledTetrahedralMesh.h"
 
+#include <iostream>
+
+std::vector< std::vector<int> > LabelledTetrahedralMesh::_tetra_triangle_indices = {{0, 1, 2}, {1, 3, 2}, {2, 3, 0}, {3, 1, 0}};
+
 
 Eigen::ArrayX3d LabelledTetrahedralMesh::vertexPoints(void) const
 {
@@ -38,9 +42,53 @@ Eigen::Array<unsigned short, Eigen::Dynamic, 1> LabelledTetrahedralMesh::tetraLa
     return tetra_labels;
 }
 
+Eigen::Array<unsigned short, Eigen::Dynamic, 3> LabelledTetrahedralMesh::triangleVertices(void) const
+{
+    Eigen::Array<unsigned short, Eigen::Dynamic, 3> triangle_vertices(this->_triangle_vertices.size(), 3);
+    for (int i=0; i<this->_triangle_vertices.size(); i++) {
+        for (int v=0; v<3; v++) {
+            triangle_vertices(i, v) = this->_triangle_vertices[i][v];
+        }
+    }
+    return triangle_vertices;
+}
+
+Eigen::Array<int, Eigen::Dynamic, 2> LabelledTetrahedralMesh::triangleLabels(void) const
+{
+    Eigen::Array<int, Eigen::Dynamic, 2> triangle_labels(this->_triangle_tetras.size(), 2);
+    for (int i=0; i<this->_triangle_tetras.size(); i++) {
+        triangle_labels(i, 0) = this->_tetra_labels[this->_triangle_tetras[i][0]];
+        if (this->_triangle_tetras[i].size() > 1) {
+            triangle_labels(i, 1) = this->_tetra_labels[this->_triangle_tetras[i][1]];
+        } else {
+            triangle_labels(i, 1) = -1;
+        }
+    }
+    return triangle_labels;
+}
+
 void LabelledTetrahedralMesh::addTetra(Eigen::Array<unsigned short, 4, 1>  tetra, unsigned short label)
 {
     std::vector<unsigned short> _tetra = std::vector<unsigned short> {tetra(0), tetra(1), tetra(2), tetra(3)};
+    unsigned short tetra_index = static_cast<unsigned short>(this->_tetra_vertices.size());
+
+    for (auto triangle_indices : LabelledTetrahedralMesh::_tetra_triangle_indices) {
+        std::vector<unsigned short> _triangle;
+        for (auto i_t : triangle_indices) {
+            _triangle.push_back(tetra(i_t));
+        }
+        std::sort(_triangle.begin(), _triangle.end());
+
+        auto triangle_match = std::find(this->_triangle_vertices.begin(), this->_triangle_vertices.end(), _triangle);
+        if (triangle_match == this->_triangle_vertices.end()) {
+            this->_triangle_vertices.push_back(_triangle);
+            this->_triangle_tetras.push_back({ tetra_index });
+        } else {
+            int triangle_index = triangle_match-_triangle_vertices.begin();
+            this->_triangle_tetras[triangle_index].push_back(tetra_index);
+        }
+    }
+
     this->_tetra_vertices.push_back(_tetra);
     this->_tetra_labels.push_back(label);
 }

cpp/LabelledTetrahedralMesh.h

@@ -17,10 +17,21 @@ public:
     Eigen::Array<unsigned short, Eigen::Dynamic, 1> tetraLabels(void) const;
     void addTetra(Eigen::Array<unsigned short, 4, 1>  tetra, unsigned short label);
 
+public:
+    Eigen::Array<unsigned short, Eigen::Dynamic, 3> triangleVertices(void) const;
+    Eigen::Array<int, Eigen::Dynamic, 2> triangleLabels(void) const;
+
 private:
     std::vector< std::vector<double> > _vertex_points;
     std::vector< std::vector<unsigned short> > _tetra_vertices;
     std::vector< unsigned short > _tetra_labels;
+
+    std::vector< std::vector<unsigned short> > _triangle_vertices;
+    std::vector< std::vector<unsigned short> > _triangle_tetras;
+
+
+private:
+    static std::vector< std::vector<int> > _tetra_triangle_indices;
 };
 
 LabelledTetrahedralMesh makeTetrahedron(void);
\ No newline at end of file

src/pybind_example/labelled_tetrahedral_mesh.py

-from typing import Iterable, Optional
+from typing import Iterable
 
 import numpy as np
 
 from ._wrp import LabelledTetrahedralMesh, make_tetrahedron
 from ._wrp import build_cgal_image as _build_cgal_image
+from ._wrp import cgal_image_meshing as _cgal_image_meshing
 
 
-def build_cgal_image(img: np.ndarray, voxelsize: Optional[Iterable[float]]=(1., 1., 1.)) -> bool:
+def build_cgal_image(img: np.ndarray, voxelsize: Iterable[float]=(1., 1., 1.)) -> bool:
     return _build_cgal_image(img.ravel().astype(np.uint16), img.shape, voxelsize)
+
+
+def cgal_image_meshing(
+    img: np.ndarray, voxelsize: Iterable[float]=(1., 1., 1.),
+    facet_angle: float=30., facet_size: float=6., facet_distance: float=4.,
+    cell_radius_edge_ratio: float=3., cell_size: float=8.
+) -> LabelledTetrahedralMesh:
+    return _cgal_image_meshing(
+        img.ravel().astype(np.uint16), img.shape, voxelsize,
+        facet_angle, facet_size, facet_distance, cell_radius_edge_ratio, cell_size
+    )

tests/test_cgal_meshing.py

+import unittest
+
+import numpy as np
+
+from pybind_example.labelled_tetrahedral_mesh import cgal_image_meshing
+
+
+class TestCgalImage(unittest.TestCase):
+
+    def setUp(self):
+        self.image = np.zeros((10, 8, 6), dtype=np.uint16)
+        self.image[:, :, 4:] = 1
+        self.image[:, :4, :4] = 2
+        self.image[:, 4:, :4] = 3
+
+        self.voxelsize = [0.5, 0.5, 1]
+
+    def tearDown(self):
+        pass
+
+    def test_cgal_meshing(self):
+        mesh = cgal_image_meshing(self.image, voxelsize=self.voxelsize, cell_size=1.)
+        tetra_labels = mesh.tetra_labels()
+        assert len(tetra_labels) > 0
+        assert all([l in tetra_labels for l in np.unique(self.image)])
+
+        n_tetras = len(tetra_labels)
+
+        refined_mesh = cgal_image_meshing(self.image, voxelsize=self.voxelsize, cell_size=0.5)
+        assert len(refined_mesh.tetra_labels()) > n_tetras

wrp/LabelledTetrahedralMeshWrapper.cpp

@@ -3,6 +3,7 @@
 #include <pybind11/eigen.h>
 
 #include "CgalLabelledImage.h"
+#include "CgalImageMeshing.h"
 #include "LabelledTetrahedralMesh.h"
 
 namespace py = pybind11;
@@ -24,9 +25,13 @@ PYBIND11_MODULE(_wrp, m) {
         .def("add_point", &LabelledTetrahedralMesh::addPoint, py::arg("point"))
         .def("tetra_vertices", &LabelledTetrahedralMesh::tetraVertices)
         .def("tetra_labels", &LabelledTetrahedralMesh::tetraLabels)
+        .def("triangle_vertices", &LabelledTetrahedralMesh::triangleVertices)
+        .def("triangle_labels", &LabelledTetrahedralMesh::triangleLabels)
         .def("add_tetra", &LabelledTetrahedralMesh::addTetra, py::arg("tetra"), py::arg("label"));
 
     m.def("make_tetrahedron", &makeTetrahedron);
 
     m.def("build_cgal_image", &cgalLabelledImage, py::arg("image_array"), py::arg("shape"), py::arg("voxelsize"));
+
+    m.def("cgal_image_meshing", &cgalImageMeshing, py::arg("image_array"), py::arg("shape"), py::arg("voxelsize"), py::arg("facet_angle"), py::arg("facet_size"), py::arg("facet_distance"), py::arg("cell_radius_edge_ratio"), py::arg("cell_size"));
 }
\ No newline at end of file