Add tests

Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt

@@ -7,6 +7,7 @@ project( Isosurfacing_3_tests )
 find_package(CGAL REQUIRED)
 
 create_single_source_cgal_program( "test.cpp" )
+create_single_source_cgal_program( "test_marching_cubes.cpp" )
 
 find_package(OpenMP)
 
@@ -15,4 +16,6 @@ include(CGAL_TBB_support)
 if(TARGET CGAL::TBB_support)
   target_link_libraries(test PUBLIC CGAL::TBB_support)
   target_link_libraries(test PRIVATE OpenMP::OpenMP_CXX)
+
+  target_link_libraries(test_marching_cubes PUBLIC CGAL::TBB_support)
 endif()

Isosurfacing_3/test/Isosurfacing_3/test.cpp

@@ -11,7 +11,6 @@
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/TC_marching_cubes_3.h>
 #include <CGAL/boost/graph/IO/OFF.h>
-
 #include <tbb/concurrent_vector.h>
 
 #include "Timer.h"
@@ -27,21 +26,23 @@ typedef tbb::concurrent_vector<Point> Point_range;
 typedef tbb::concurrent_vector<std::vector<std::size_t>> Polygon_range;
 
 int main() {
-    const Vector spacing(0.002f, 0.002f, 0.02f);
+    const Vector spacing(0.002, 0.002, 0.02);
     const CGAL::Bbox_3 bbox = {-1, -1, -1, 1, 1, 1};
 
     auto sphere_function = [](const Point& point) {
         return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z());
     };
-    CGAL::Isosurfacing::Implicit_domain<
-        Kernel, decltype(sphere_function),
-        decltype(CGAL::Isosurfacing::Default_gradient<Kernel, decltype(sphere_function)>(sphere_function))>
-        implicit_domain({-1, -1, -1, 1, 1, 1}, spacing, sphere_function,
-                        CGAL::Isosurfacing::Default_gradient<Kernel, decltype(sphere_function)>(
-                            sphere_function));  // TODO: this is ugly
 
+    typedef CGAL::Isosurfacing::Finite_difference_gradient<Kernel, decltype(sphere_function)> Gradient;
 
-    Grid grid(2.f / spacing.x(), 2.f / spacing.y(), 2.f / spacing.z(), bbox);
+    CGAL::Isosurfacing::Implicit_domain<Kernel, decltype(sphere_function), Gradient> implicit_domain(
+        {-1, -1, -1, 1, 1, 1}, spacing, sphere_function, Gradient(sphere_function, 0.0001));  // TODO: this is ugly
+
+    const std::size_t nx = static_cast<std::size_t>(2.0 / spacing.x());
+    const std::size_t ny = static_cast<std::size_t>(2.0 / spacing.y());
+    const std::size_t nz = static_cast<std::size_t>(2.0 / spacing.z());
+
+    Grid grid(nx, ny, nz, bbox);
 
     for (std::size_t x = 0; x < grid.xdim(); x++) {
         for (std::size_t y = 0; y < grid.ydim(); y++) {
@@ -71,15 +72,15 @@ int main() {
 
     {
         ScopeTimer timer;
-        CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes<CGAL::Parallel_tag>(implicit_domain, 0.8f, points,
-                                                                                        polygons);
+        CGAL::Isosurfacing::make_quad_mesh_using_dual_contouring<CGAL::Parallel_tag>(implicit_domain, 0.8f, points,
+                                                                                     polygons);
     }
 
     // TODO: compare results with mesh_3
 
-    //Mesh mesh;
-    //CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh);
+    // Mesh mesh;
+    // CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh);
 
-    //CGAL::IO::write_OFF("result.off", mesh);
+    // CGAL::IO::write_OFF("result.off", mesh);
     CGAL::IO::write_OFF("result.off", points, polygons);
 }

Isosurfacing_3/test/Isosurfacing_3/test_marching_cubes.cpp

@@ -0,0 +1,96 @@
+
+#include <CGAL/Marching_cubes_3.h>
+
+#include "test_util.h"
+
+#define WRITE_OFF
+
+struct Sphere_function {
+    FT operator()(const Point& point) const {
+        return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z());
+    }
+};
+
+template <class Domain_>
+void run(const Domain_& domain, const FT iso_value, Point_range& points, Polygon_range& polygons) {
+    CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes<CGAL::Parallel_tag>(domain, iso_value, points,
+                                                                                    polygons);
+}
+
+void test_implicit_sphere() {
+    const std::string test_name = "test_implicit_sphere()";
+
+    const Vector spacing(0.2, 0.2, 0.2);
+    const CGAL::Bbox_3 bbox = {-1, -1, -1, 1, 1, 1};
+
+    CGAL::Isosurfacing::Implicit_domain<Kernel, Sphere_function> domain(bbox, spacing,
+                                                                        Sphere_function());  // TODO: this is ugly
+
+    Point_range points;
+    Polygon_range polygons;
+    run(domain, 0.8, points, polygons);
+
+#ifdef WRITE_OFF
+    CGAL::IO::write_OFF(test_name + ".off", points, polygons);
+#endif
+
+    assert(is_polygon_mesh(polygons));
+    Mesh m = to_mesh(points, polygons);
+
+    assert(is_manifold(m));
+    assert(!has_degenerate_faces(m));
+
+    std::cout << "Test passed: " << test_name << std::endl;
+}
+
+void test_grid_sphere(const std::size_t n) {
+    const std::string test_name = "test_grid_sphere(" + std::to_string(n) + ")";
+
+    const CGAL::Bbox_3 bbox = {-1, -1, -1, 1, 1, 1};
+    const Vector spacing(2.0 / (n - 1), 2.0 / (n - 1), 2.0 / (n - 1));
+
+    Sphere_function sphere_function;
+
+    Grid grid(n, n, n, bbox);
+
+    for (std::size_t x = 0; x < grid.xdim(); x++) {
+        for (std::size_t y = 0; y < grid.ydim(); y++) {
+            for (std::size_t z = 0; z < grid.zdim(); z++) {
+
+                const Point pos(x * spacing.x() + bbox.xmin(), y * spacing.y() + bbox.ymin(),
+                                z * spacing.z() + bbox.zmin());
+
+                grid.value(x, y, z) = sphere_function(pos);
+            }
+        }
+    }
+
+    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
+
+    Point_range points;
+    Polygon_range polygons;
+    run(domain, 0.777, points, polygons);
+
+#ifdef WRITE_OFF
+    CGAL::IO::write_OFF(test_name + ".off", points, polygons);
+#endif
+
+    assert(is_polygon_mesh(polygons));
+    Mesh m = to_mesh(points, polygons);
+
+    assert(is_manifold(m));
+    assert(!has_degenerate_faces(m));
+
+    std::cout << "Test passed: " << test_name << std::endl;
+}
+
+int main() {
+    test_implicit_sphere();
+    test_grid_sphere(2);
+    test_grid_sphere(3);
+    test_grid_sphere(10);
+    test_grid_sphere(11);
+    test_grid_sphere(100);
+
+    std::cout << "All tests passed" << std::endl;
+}

Isosurfacing_3/test/Isosurfacing_3/test_util.h

@@ -0,0 +1,70 @@
+#ifndef CGAL_ISOSURFACING_TEST_UTIL_H
+#define CGAL_ISOSURFACING_TEST_UTIL_H
+
+#include <CGAL/Cartesian_grid_3.h>
+#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Implicit_domain.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+#include <CGAL/Polygon_mesh_processing/manifoldness.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/repair.h>
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Surface_mesh.h>
+#include <tbb/concurrent_vector.h>
+
+#include <cassert>
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef typename Kernel::FT FT;
+typedef typename Kernel::Vector_3 Vector;
+typedef typename Kernel::Point_3 Point;
+
+typedef CGAL::Surface_mesh<Point> Mesh;
+typedef CGAL::Cartesian_grid_3<Kernel> Grid;
+
+typedef std::vector<Point> Point_range;
+typedef std::vector<std::vector<std::size_t>> Polygon_range;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+bool has_duplicate_points(Point_range points, Polygon_range polygons) {
+    return PMP::merge_duplicate_points_in_polygon_soup(points, polygons) != 0;
+}
+
+bool has_duplicate_polygons(Point_range points, Polygon_range polygons) {
+    return PMP::merge_duplicate_polygons_in_polygon_soup(points, polygons) != 0;
+}
+
+bool has_isolated_vertices(Point_range points, Polygon_range polygons) {
+    return PMP::remove_isolated_points_in_polygon_soup(points, polygons) != 0;
+}
+
+bool is_polygon_mesh(const Polygon_range& polygons) {
+    return PMP::is_polygon_soup_a_polygon_mesh(polygons);
+}
+
+Mesh to_mesh(const Point_range& points, const Polygon_range& polygons) {
+    Mesh m;
+    PMP::polygon_soup_to_polygon_mesh(points, polygons, m);
+    return m;
+}
+
+bool is_manifold(Mesh& m) {
+    return PMP::duplicate_non_manifold_vertices(m, CGAL::parameters::dry_run(true)) == 0;
+}
+
+bool has_degenerate_faces(Mesh& m) {
+    return PMP::remove_connected_components_of_negligible_size(
+               m, CGAL::parameters::dry_run(true).area_threshold(std::numeric_limits<FT>::epsilon())) != 0;
+}
+
+bool check_mesh_distance(const Mesh& m0, const Mesh& m1) {
+    auto dist = PMP::approximate_Hausdorff_distance<CGAL::Sequential_tag>(
+        m0, m1, CGAL::parameters::number_of_points_per_area_unit(4000));
+    std::cout << dist << std::endl;
+    return true;
+}
+
+#endif  // CGAL_ISOSURFACING_TEST_UTIL_H