Add benchmarks

Isosurfacing_3/benchmark/Isosurfacing_3/CMakeLists.txt

@@ -0,0 +1,15 @@
+# Created by the script cgal_create_cmake_script
+# This is the CMake script for compiling a CGAL application.
+
+cmake_minimum_required(VERSION 3.1...3.14)
+project( Isosurfacing_3_benchmark )
+
+find_package(CGAL REQUIRED)
+
+create_single_source_cgal_program( "benchmark.cpp" )
+
+find_package(TBB)
+include(CGAL_TBB_support)
+if(TARGET CGAL::TBB_support)
+  target_link_libraries(benchmark PUBLIC CGAL::TBB_support)
+endif()

Isosurfacing_3/benchmark/Isosurfacing_3/Timer.h

@@ -0,0 +1,43 @@
+#ifndef CGAL_TIMER_H
+#define CGAL_TIMER_H
+
+#include <chrono>
+#include <iostream>
+
+class ScopeTimer {
+    using Clock = std::chrono::steady_clock;
+    using TimePoint = Clock::time_point;
+
+public:
+    ScopeTimer() : start(Clock::now()), msg("Duration"), running(true) {}
+
+    explicit ScopeTimer(const std::string& msg) : start(Clock::now()), msg(msg), running(true) {
+        std::cout << msg << "..." << std::endl;
+    }
+
+    void reset() {
+        start = Clock::now();
+        running = true;
+    }
+
+    int64_t stop() {
+        TimePoint end = Clock::now();
+        running = false;
+        return std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+    }
+
+    ~ScopeTimer() {
+        if (running) {
+            TimePoint end = Clock::now();
+            int64_t duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+            std::cout << msg << ": " << duration << " ms" << std::endl;
+        }
+    }
+
+private:
+    TimePoint start;
+    std::string msg;
+    bool running;
+};
+
+#endif  // CGAL_TIMER_H

Isosurfacing_3/benchmark/Isosurfacing_3/benchmark.cpp

@@ -0,0 +1,169 @@
+#include <CGAL/Cartesian.h>
+#include <CGAL/Cartesian_grid_3.h>
+#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Dual_contouring_3.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Implicit_domain.h>
+#include <CGAL/Marching_cubes_3.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/boost/graph/IO/OFF.h>
+#include <math.h>
+
+#include <iostream>
+#include <limits>
+
+#include "CLI11.hpp"
+#include "Timer.h"
+
+template <class GeomTraits>
+struct SphereValue {
+    typename GeomTraits::FT operator()(const typename GeomTraits::Point_3& point) const {
+        return CGAL::approximate_sqrt((point - CGAL::ORIGIN).squared_length());
+    }
+};
+
+template <class GeomTraits>
+struct SphereGradient {
+    typename GeomTraits::Vector_3 operator()(const typename GeomTraits::Point_3& point) const {
+        const typename GeomTraits::Vector_3 g = point - CGAL::ORIGIN;
+        return g / CGAL::approximate_sqrt(g.squared_length());
+    }
+};
+
+template <class GeomTraits>
+struct Implicit_sphere {
+
+    typedef CGAL::Isosurfacing::Implicit_domain<GeomTraits, SphereValue<GeomTraits>, SphereGradient<GeomTraits>> Domain;
+    typedef typename GeomTraits::Vector_3 Vector;
+
+    Domain domain(const std::size_t N) const {
+
+        const Vector res(2.0 / N, 2.0 / N, 2.0 / N);
+        return Domain({-1, -1, -1, 1, 1, 1}, res, SphereValue<GeomTraits>(), SphereGradient<GeomTraits>());
+    }
+
+    typename GeomTraits::FT iso() const {
+        return 0.8;
+    }
+};
+
+template <class GeomTraits>
+struct Grid_sphere {
+
+    typedef CGAL::Isosurfacing::Cartesian_grid_domain<GeomTraits> Domain;
+    typedef CGAL::Cartesian_grid_3<GeomTraits> Grid;
+    typedef typename GeomTraits::FT FT;
+    typedef typename GeomTraits::Point_3 Point;
+
+    Domain domain(const std::size_t N) const {
+
+        const FT resolution = 2.0 / N;
+        SphereValue<GeomTraits> val;
+        SphereGradient<GeomTraits> grad;
+
+        Grid grid(N, N, N, {-1, -1, -1, 1, 1, 1});
+
+        for (std::size_t x = 0; x < grid.xdim(); x++) {
+            const FT xp = x * resolution - 1.0;
+
+            for (std::size_t y = 0; y < grid.ydim(); y++) {
+                const FT yp = y * resolution - 1.0;
+
+                for (std::size_t z = 0; z < grid.zdim(); z++) {
+                    const FT zp = z * resolution - 1.0;
+
+                    grid.value(x, y, z) = val(Point(xp, yp, zp));
+                    grid.gradient(x, y, z) = grad(Point(xp, yp, zp));
+                }
+            }
+        }
+        return Domain(grid);
+    }
+
+    typename GeomTraits::FT iso() const {
+        return 0.8;
+    }
+};
+
+template <class GeomTraits>
+struct Skull_image {
+
+    typedef CGAL::Isosurfacing::Cartesian_grid_domain<GeomTraits> Domain;
+    typedef CGAL::Cartesian_grid_3<GeomTraits> Grid;
+
+    Domain domain(const std::size_t N) const {
+
+        const std::string fname = "../data/skull_2.9.inr";
+        CGAL::Image_3 image;
+        if (!image.read(fname)) {
+            std::cerr << "Error: Cannot read file " << fname << std::endl;
+            return EXIT_FAILURE;
+        }
+
+        const Grid grid(image);
+        return Domain(grid);
+    }
+
+    typename GeomTraits::FT iso() const {
+        return 2.9;
+    }
+};
+
+int main(int argc, char* argv[]) {
+    CLI::App app{"Isosurfacing benchmarks"};
+
+    std::size_t N = 2;
+    app.add_option("-N", N, "Grid size");
+
+    CLI11_PARSE(app, argc, argv);
+
+#if defined KERNEL_SIMPLE_CARTESIAN_DOUBLE
+    typedef CGAL::Simple_cartesian<double> Kernel;
+#elif defined KERNEL_SIMPLE_CARTESIAN_FLOAT
+    typedef CGAL::Simple_cartesian<float> Kernel;
+#elif defined KERNEL_CARTESIAN_DOUBLE
+    typedef CGAL::Cartesian<double> Kernel;
+#elif defined KERNEL_EPIC
+    typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+#else
+    typedef CGAL::Simple_cartesian<double> Kernel;
+#endif
+
+    typedef Kernel::Point_3 Point;
+
+    typedef std::vector<Point> Point_range;
+    typedef std::vector<std::vector<std::size_t>> Polygon_range;
+
+#if defined SCENARIO_GRID_SPHERE
+    auto scenario = Grid_sphere<Kernel>();
+#elif defined SCENARIO_IMPLICIT_SPHERE
+    auto scenario = Implicit_sphere<Kernel>();
+#elif defined SCENARIO_SKULL_IMAGE
+    auto scenario = Skull_image<Kernel>();
+#else
+    auto scenario = Implicit_sphere<Kernel>();
+#endif
+
+    Point_range points;
+    Polygon_range polygons;
+
+    ScopeTimer timer;
+
+#if defined ALGO_MARCHING_CUBES
+    CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(scenario.domain(N), scenario.iso(), points, polygons);
+#elif defined ALGO_DUAL_CONTOURING
+    CGAL::Isosurfacing::make_quad_mesh_using_dual_contouring(scenario.domain(N), scenario.iso(), points, polygons);
+#else
+    CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(scenario.domain(N), scenario.iso(), points, polygons);
+#endif
+
+    const int64_t ms = timer.stop();
+
+    if (points.size() > std::numeric_limits<std::size_t>::max() - 2) {
+        std::cout << "This should never print and only prevents optimizations" << std::endl;
+    }
+
+    std::cout << ms << std::endl;
+}

Isosurfacing_3/benchmark/Isosurfacing_3/benchmark_size.py

@@ -0,0 +1,26 @@
+from benchmark_util import *
+
+scenario = "SCENARIO_GRID_SPHERE"
+kernel = "KERNEL_SIMPLE_CARTESIAN_FLOAT"
+algorithm = "ALGO_DUAL_CONTOURING"
+threads = 1
+exponent = 2
+min_cells = 10
+max_cells = 1000
+
+build(scenario, kernel, algorithm)
+
+data = []
+
+c = min_cells
+while c < max_cells:
+	n = int(c ** (1.0 / 3.0))
+
+	time = execute(n)
+	data.append([scenario, kernel, algorithm, threads, c, time])
+
+	c *= exponent
+
+df = pd.DataFrame(data, comulmns=["scenario", "kernel", "algorithm", "threads", "cells", "time"])
+
+df.to_csv("benchmark_size.csv")

Isosurfacing_3/benchmark/Isosurfacing_3/benchmark_threads.py

@@ -0,0 +1,21 @@
+from benchmark_util import *
+
+scenario = "SCENARIO_GRID_SPHERE"
+kernel = "KERNEL_SIMPLE_CARTESIAN_FLOAT"
+algorithm = "ALGO_DUAL_CONTOURING"
+min_threads = 1
+max_threads = 8
+cells = 1000000
+
+build(scenario, kernel, algorithm)
+
+data = []
+
+for t in range(min_threads, max_threads):
+	time = execute(cells)
+
+	data.append([scenario, kernel, algorithm, t, cells, 0])
+
+df = pd.DataFrame(data, comulmns=["scenario", "kernel", "algorithm", "threads", "cells", "time"])
+
+df.to_csv("benchmark_threads.csv")

Isosurfacing_3/benchmark/Isosurfacing_3/benchmark_util.py

@@ -0,0 +1,33 @@
+import os
+import sys
+import subprocess
+import pandas as pd
+
+def print_stream(stream):
+    while True:
+        line = stream.readline()
+        if not line:
+            break
+        print(line, end="")
+
+def run(cmd, output=True):
+    process = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE, universal_newlines=True)
+    exit_code = process.wait()
+    if output:
+        print_stream(process.stdout)
+    if exit_code != 0:
+        print_stream(process.stderr)
+        sys.exit(exit_code)
+    return process
+
+def build(scenario, kernel, algorithm):
+    run(["cmake", "-E", "make_directory", "build"])
+    run(["cmake", "-B", "build", "-DCMAKE_BUILD_TYPE=Release", "-DCGAL_DIR=../../../", "-D" + scenario + "= -D" + kernel + "= -D" + algorithm + "="])
+    run(["cmake", "--build", "build"])
+
+def execute(n, times=1):
+    time = 0
+    for i in range(times):
+        run(["./build/benchmark", "-N", str(n)], False)
+        time += int(process.stdout.readline())
+    return time / times

Isosurfacing_3/benchmark/Isosurfacing_3/graphs.py

@@ -0,0 +1,63 @@
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def save_svg(file):
+    fig = plt.gcf()
+    fig.set_size_inches((10, 5), forward=False)
+    plt.savefig(file, bbox_inches="tight")
+
+
+def add_threads_graph(data, label):
+    x = data["threads"]
+    y = data["cells"] / data["time"] / 10 ** 3
+    plt.plot(x, y, label=label)
+    plt.legend()
+
+
+def add_size_graph(data, label):
+    x = data["cells"]
+    y = data["cells"] / data["time"] / 10 ** 3
+    plt.plot(x, y, label=label)
+    plt.legend()
+
+
+def plot_graph(file, name, log, ylabel, xlabel):
+    plt.title(name)
+    plt.xlabel(xlabel)
+    plt.ylabel(ylabel)
+    if log:
+        plt.xscale("log")
+    plt.gca().yaxis.grid(color='#cccccc')
+    plt.gca().xaxis.grid(color='#cccccc')
+    plt.ylim(ymin=0)
+    save_svg(file)
+    plt.show()
+
+
+latex_export = True
+if latex_export:
+    #plt.rcParams["svg.fonttype"] = "none"
+    plt.rcParams["axes.unicode_minus"] = False
+    plt.rcParams['font.size'] = "17"
+
+
+data = pd.read_csv("threads_implicit.csv", sep=";")
+add_threads_graph(data, "implicit")
+
+data = pd.read_csv("threads_grid.csv", sep=";")
+add_threads_graph(data, "grid")
+
+plt.xticks(np.arange(1, max(data["threads"]) + 0.1, 1))
+plot_graph("perf_threads.svg", "", False, "performance [10^3 cubes/s]", "cores")
+
+
+data = pd.read_csv("size_implicit.csv", sep=";")
+add_size_graph(data, "implicit")
+
+data = pd.read_csv("size_grid.csv", sep=";")
+add_size_graph(data, "grid")
+
+plot_graph("perf_size.svg", "", True, "performance [10^3 cubes/s]", "cubes")