diff --git a/CMakeLists.txt b/CMakeLists.txt index e008514..f576dcc 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.0) +cmake_minimum_required(VERSION 3.11) list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake/Modules") project(Discregrid) @@ -6,10 +6,26 @@ project(Discregrid) # Visual studio solution directories. set_property(GLOBAL PROPERTY USE_FOLDERS on) -# Require C++11 compiler -set(CMAKE_CXX_STANDARD 11) -set(CMAKE_CXX_STANDARD_REQUIRED ON) +# 2D SDF generation requires Clipper2 Library by AngusJohnson +option(USE_CLIPPER2 "Use Clipper2 library to enable 2D SDF generation" OFF) + +if (USE_CLIPPER2) + # Require C++17 compiler as Clipper2 does + set(CMAKE_CXX_STANDARD 17) + include(FetchContent) + FetchContent_Declare( + clipper2 + GIT_REPOSITORY https://github.com/AngusJohnson/Clipper2.git + GIT_TAG 81b01d2acbad7b06fb28df4fbfbd7228519b7905 + SOURCE_SUBDIR CPP + ) + FetchContent_MakeAvailable(clipper2) +else() + # Require C++11 compiler + set(CMAKE_CXX_STANDARD 11) +endif() +set(CMAKE_CXX_STANDARD_REQUIRED ON) # Enable simultaneous compilation of source files. if(MSVC) diff --git a/Dragon2D.png b/Dragon2D.png new file mode 100644 index 0000000..62ed90b Binary files /dev/null and b/Dragon2D.png differ diff --git a/README.md b/README.md index d47bd5b..24cd5bf 100755 --- a/README.md +++ b/README.md @@ -6,13 +6,16 @@ **Figure 1**: Left: Slice of a three-dimensional discrete signed distance field of the Stanford dragon. Right: Density map for SPH boundary handling of Stanford dragon. **Discregrid** is a static C++ library for the parallel discretization of (preferably smooth) functions on regular grids. -The library generates a (cubic) polynomial discretization given a box-shaped domain, a grid resolution, and a function that maps a three-dimensional position in space to a real scalar value. +The library generates a (cubic) polynomial discretization given a box-shaped domain, a grid resolution, and a function that maps a two- or three-dimensional position in space to a real scalar value. Isoparametric cubic polynomials of Serendipity type for the cell-wise discretization are employed. The coefficient vector for the discrete polynomial basis is computed using regular sampling of the input function at the higher-order grid's nodes. The algorithm to generate the discretization is moreover *fully parallelized* using OpenMP and especially well-suited for the discretization of signed distance functions. The library moreover provides the functionality to serialize and deserialize the a generated discrete grid. Discregrid ships with [TriangleMeshDistance](https://github.com/InteractiveComputerGraphics/TriangleMeshDistance) to directly provide the capability to compute and discretize signed distance fields to triangle meshes. +![](Dragon2D.png) +**Figure 2**: Left: 2D-Polygon representation of the Stanford dragon. Right: Two-dimensional discrete signed distance field of the Stanford dragon. + Besides the library, the project includes three executable programs that serve the following purposes: * *GenerateSDF*: Computes a discrete (cubic) signed distance field from a triangle mesh in OBJ format. * *DiscreteFieldToBitmap*: Generates an image in bitmap format of a two-dimensional slice of a previously computed discretization. @@ -26,9 +29,12 @@ Besides the library, the project includes three executable programs that serve t ## Build Instructions -This project is based on [CMake](https://cmake.org/). Simply generate project, Makefiles, etc. using [CMake](https://cmake.org/) and compile the project with the compiler of your choice. The code was tested with the following configurations: -- Windows 10 64-bit, CMake 3.8, Visual Studio 2017 -- Debian 9 64-bit, CMake 3.8, GCC 6.3. +This project is based on [CMake](https://cmake.org/). Simply generate project, Makefiles, etc. using [CMake](https://cmake.org/) and compile the project with the compiler of your choice. The minimum required version of CMake for this project is 3.11. The code was tested with the following configurations: +- Windows 10 64-bit, CMake 3.23.1, Visual Studio 2017 + +If the option USE_CLIPPER2 is active in CMake the [Clipper2](https://github.com/AngusJohnson/Clipper2) library is automatically downloaded and the following two executable programs can also be generated: +* *GenerateSDF2D*: Computes a discrete (cubic) 2D signed distance field from a triangle mesh in OBJ format. Clipper2 is used to convert the triangle mesh to a 2D polygon before discretization, as seen on the left in Figure 2. +* *DiscreteFieldToBitmap2D*: Generates an image in bitmap format of a previously computed 2D discretization. ## Usage In order to use the library, the main header has to be included and the static library has to be compiled and linked against the client program. diff --git a/cmd/CMakeLists.txt b/cmd/CMakeLists.txt index 1b5be71..1d78479 100755 --- a/cmd/CMakeLists.txt +++ b/cmd/CMakeLists.txt @@ -1,3 +1,8 @@ add_subdirectory(generate_sdf) add_subdirectory(discrete_field_to_bitmap) add_subdirectory(generate_density_map) + +if (USE_CLIPPER2) + add_subdirectory(generate_sdf_2d) + add_subdirectory(discrete_field_to_bitmap_2d) +endif() diff --git a/cmd/discrete_field_to_bitmap/CMakeLists.txt b/cmd/discrete_field_to_bitmap/CMakeLists.txt index a9ab1b6..ea10728 100755 --- a/cmd/discrete_field_to_bitmap/CMakeLists.txt +++ b/cmd/discrete_field_to_bitmap/CMakeLists.txt @@ -32,8 +32,6 @@ endif() add_executable(DiscreteFieldToBitmap main.cpp - bmp_file.hpp - bmp_file.cpp ) add_dependencies(DiscreteFieldToBitmap diff --git a/cmd/discrete_field_to_bitmap/main.cpp b/cmd/discrete_field_to_bitmap/main.cpp index 6ce4fe1..c6666d4 100755 --- a/cmd/discrete_field_to_bitmap/main.cpp +++ b/cmd/discrete_field_to_bitmap/main.cpp @@ -7,7 +7,7 @@ #include #include -#include "bmp_file.hpp" +#include "Discregrid/utility/bmp_file.hpp" using namespace Eigen; diff --git a/cmd/discrete_field_to_bitmap_2d/CMakeLists.txt b/cmd/discrete_field_to_bitmap_2d/CMakeLists.txt new file mode 100644 index 0000000..89d5375 --- /dev/null +++ b/cmd/discrete_field_to_bitmap_2d/CMakeLists.txt @@ -0,0 +1,45 @@ +# Eigen library. +find_package(Eigen3 REQUIRED) + +# Set include directories. +include_directories( + ../../extern + ../../discregrid/include + ${EIGEN3_INCLUDE_DIR} +) + + +if(WIN32) + add_definitions(-D_SCL_SECURE_NO_WARNINGS) + add_definitions(-D_CRT_SECURE_NO_WARNINGS) +endif(WIN32) + +if ( CMAKE_COMPILER_IS_GNUCC ) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-multichar") +endif ( CMAKE_COMPILER_IS_GNUCC ) + +# OpenMP support. +find_package(OpenMP REQUIRED) +if(OPENMP_FOUND) + if (CMAKE_VERSION VERSION_GREATER "3.8") + link_libraries(OpenMP::OpenMP_CXX) + else() + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") + endif() +endif() + +add_executable(DiscreteFieldToBitmap2D + main.cpp +) + +add_dependencies(DiscreteFieldToBitmap2D + Discregrid +) + +target_link_libraries(DiscreteFieldToBitmap2D + Discregrid +) + +set_target_properties(DiscreteFieldToBitmap2D PROPERTIES FOLDER Cmd) diff --git a/cmd/discrete_field_to_bitmap_2d/main.cpp b/cmd/discrete_field_to_bitmap_2d/main.cpp new file mode 100644 index 0000000..459d1e2 --- /dev/null +++ b/cmd/discrete_field_to_bitmap_2d/main.cpp @@ -0,0 +1,216 @@ + +#include +#include +#include + +#include +#include +#include + +#include "Discregrid/utility/bmp_file.hpp" + +using namespace Eigen; + +namespace +{ +std::array doubleToGreenBlueInverse(double v) +{ + if (v >= 0.0) + { + return {{0u, static_cast(std::min(std::max(255.0 * (1.0 - v), 0.0), 255.0)), 0u}}; + } + return {{0u, 0u, static_cast(std::min(std::max(255.0 * (1.0 + v), 0.0), 255.0))}}; +} + +std::array doubleToRedSequential(double v) +{ + return {{static_cast(std::min(std::max(255.0 * v, 0.0), 255.0)), 0u, 0u}}; +} + +std::array doubleTo5ColourHeatmap(double v, double minValue, double maxValue) +{ + const double normalised = (v - minValue) / (maxValue - minValue); + + constexpr double h1 = 0; + constexpr double h2 = 0.66667; + const double h = (1.0 - normalised) * h1 + normalised * h2; + constexpr double s = 1.0; + constexpr double l = 0.5; + + auto hue2rgb = [](double p, double q, double t){ + if(t < 0) t += 1; + if(t > 1) t -= 1; + if(t < 1/6.0) return p + (q - p) * 6 * t; + if(t < 1/2.0) return q; + if(t < 2/3.0) return p + (q - p) * (2/3.0 - t) * 6; + return p; + }; + + constexpr double q = l < 0.5 ? l * (1 + s) : l + s - l * s; + constexpr double p = 2 * l - q; + + return { + static_cast(hue2rgb(p, q, h + 1 / 3.0) * 255), + static_cast(hue2rgb(p, q, h) * 255), + static_cast(hue2rgb(p, q, h - 1 / 3.0) * 255) + }; +} + +} + +int main(int argc, char* argv[]) +{ + cxxopts::Options options(argv[0], "Transforms a slice of a discrete 2D-SDF to a bitmap image."); + options.positional_help("[input 2D-SDF file (.cdf2d)]"); + + options.add_options() + ("h,help", "Prints this help text") + ("f,field_id", "ID in which the SDF to export is stored.", cxxopts::value()->default_value("0")) + ("s,samples", "Number of samples in width direction", cxxopts::value()->default_value("1024")) + ("o,output", "Output (in bmp format)", cxxopts::value()->default_value("")) + ("c,colormap", "Color map options: redsequential (rs), green blue inverse diverging (gb) (suitable for visualisation of signed distance fields), 5 colour heatmap (hm) (suitable for visualisation of differences/errors)", cxxopts::value()->default_value("gb")) + ("input", "SDF file", cxxopts::value>()) + ; + + try + { + options.parse_positional("input"); + auto result = options.parse(argc, argv); + + if (result.count("help")) + { + std::cout << options.help() << std::endl; + std::cout << std::endl << std::endl << "Example: SDFToBitmap2D -p xz file.sdf2d" << std::endl; + exit(0); + } + if (!result.count("input")) + { + std::cout << "ERROR: No input file given." << std::endl; + std::cout << options.help() << std::endl; + std::cout << std::endl << std::endl << "Example: SDFToBitmap2D -p xz file.sdf2d" << std::endl; + exit(1); + } + + auto sdf = std::unique_ptr{}; + + auto filename = result["input"].as>().front(); + auto lastindex = filename.find_last_of("."); + auto extension = filename.substr(lastindex + 1, filename.length() - lastindex); + + std::cout << "Load SDF..."; + if (extension == "cdf2d") + { + sdf = std::make_unique(filename); + } + else + { + std::cout << "ERROR: Input file must be a '.sdf2d' file specifically." << std::endl; + std::cout << options.help() << std::endl; + std::cout << std::endl << std::endl << "Example: SDFToBitmap2D -p xz file.sdf2D" << std::endl; + exit(1); + } + std::cout << "DONE" << std::endl; + + auto const& domain = sdf->domain(); + auto diag = domain.diagonal().eval(); + + auto dir = Vector2i::Zero().eval(); + dir(1) = 1; + + auto xsamples = result["s"].as(); + auto ysamples = static_cast(std::round(diag(dir(1)) / diag(dir(0)) * static_cast(xsamples))); + + auto xwidth = diag(dir(0)) / xsamples; + auto ywidth = diag(dir(1)) / ysamples; + + auto data = std::vector{}; + data.resize(xsamples * ysamples); + + auto field_id = result["f"].as(); + + std::cout << "Sample field..."; +#pragma omp parallel for + for (int k = 0; k < static_cast(xsamples * ysamples); ++k) + { + auto i = k % xsamples; + auto j = k / xsamples; + + auto xr = static_cast(i) / static_cast(xsamples); + auto yr = static_cast(j) / static_cast(ysamples); + + auto x = domain.min()(dir(0)) + xr * diag(dir(0)) + 0.5 * xwidth; + auto y = domain.min()(dir(1)) + yr * diag(dir(1)) + 0.5 * ywidth; + + auto sample = Vector2d{}; + sample(dir(0)) = x; + sample(dir(1)) = y; + + data[k] = sdf->interpolate(field_id, sample); + if (data[k] == std::numeric_limits::max()) + { + data[k] = 0.0; + } + } + + std::cout << "DONE" << std::endl; + + auto min_v = *std::min_element(data.begin(), data.end()); + auto max_v = *std::max_element(data.begin(), data.end()); + + auto out_file = result["o"].as(); + if (out_file == "") + { + out_file = filename; + if (out_file.find(".") != std::string::npos) + { + auto lastindex = out_file.find_last_of("."); + out_file = out_file.substr(0, lastindex); + } + out_file += ".bmp"; + } + + std::cout << "Ouput file: " << out_file << std::endl; + + std::cout << "Export BMP..."; + std::transform(data.begin(), data.end(), data.begin(), [&max_v, &min_v](double v) {return v >= 0.0 ? v / std::abs(max_v) : v / std::abs(min_v); }); + + auto pixels = std::vector>(data.size()); + + auto cm = result["c"].as(); + if (cm != "gb" && cm != "rs" && cm != "hm") + { + std::cerr << "WARNING: Unknown color map option. Fallback to mode 'gb'." << std::endl; + } + + if (cm == "gb") + std::transform(data.begin(), data.end(), pixels.begin(), doubleToGreenBlueInverse); + else if (cm == "rs") + std::transform(data.begin(), data.end(), pixels.begin(), doubleToRedSequential); + else if (cm == "hm") + { + const auto min_max = std::minmax_element( + data.begin(), data.end()); + const auto min = *min_max.first; + const auto max = *min_max.second; + const auto& heatmap = [min, max](double v) { return doubleTo5ColourHeatmap(v, min, max); }; + std::transform(data.begin(), data.end(), pixels.begin(), heatmap); + } + + + BmpReaderWriter::saveFile(out_file.c_str(), xsamples, ysamples, &pixels.front()[0]); + std::cout << "DONE" << std::endl; + + std::cout << std::endl << "Statistics:" << std::endl; + std::cout << "\tdomain = " << domain.min().transpose() << ", " << domain.max().transpose() << std::endl; + std::cout << "\tmin value = " << min_v << std::endl; + std::cout << "\tmax value = " << max_v << std::endl; + std::cout << "\tbmp resolution = " << xsamples << " x " << ysamples << std::endl; + } + catch (cxxopts::OptionException const& e) + { + std::cout << "error parsing options: " << e.what() << std::endl; + exit(1); + } + + return 0; +} \ No newline at end of file diff --git a/cmd/generate_sdf_2d/CMakeLists.txt b/cmd/generate_sdf_2d/CMakeLists.txt new file mode 100644 index 0000000..3e70343 --- /dev/null +++ b/cmd/generate_sdf_2d/CMakeLists.txt @@ -0,0 +1,47 @@ +# Eigen library. +find_package(Eigen3 REQUIRED) + +# Set include directories. +include_directories( + ../../extern + ../../discregrid/include + ${EIGEN3_INCLUDE_DIR} +) + +if(WIN32) + add_definitions(-D_SCL_SECURE_NO_WARNINGS) + add_definitions(-D_USE_MATH_DEFINES) +endif(WIN32) + +# Enable PolygonDistance to use clipper2 reliant constructor. +if (USE_CLIPPER2) + add_definitions(-DIS_CLIPPER_ENABLED) +endif() + +# OpenMP support. +find_package(OpenMP REQUIRED) +if(OPENMP_FOUND) + if (CMAKE_VERSION VERSION_GREATER "3.8") + link_libraries(OpenMP::OpenMP_CXX) + else() + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") + endif() +endif() + +add_executable(GenerateSDF2D + main.cpp +) + +add_dependencies(GenerateSDF2D + Discregrid +) + +target_link_libraries(GenerateSDF2D + Discregrid + Clipper2 + Clipper2utils +) + +set_target_properties(GenerateSDF2D PROPERTIES FOLDER Cmd) diff --git a/cmd/generate_sdf_2d/main.cpp b/cmd/generate_sdf_2d/main.cpp new file mode 100644 index 0000000..7100112 --- /dev/null +++ b/cmd/generate_sdf_2d/main.cpp @@ -0,0 +1,170 @@ +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +using namespace Eigen; + +std::istream& operator>>(std::istream& is, std::array& data) +{ + is >> data[0] >> data[1]; + return is; +} + +std::istream& operator>>(std::istream& is, AlignedBox2d& data) +{ + is >> data.min()[0] >> data.min()[1] + >> data.max()[0] >> data.max()[1]; + return is; +} + +#include + +int main(int argc, char* argv[]) +{ + cxxopts::Options options(argv[0], "Generates a 2D signed distance field from a triangle mesh."); + options.positional_help("[input OBJ file]"); + + options.add_options() + ("h,help", "Prints this help text") + ("r,resolution", "Grid resolution", cxxopts::value>()->default_value("10 10")) + ("d,domain", "Domain extents (bounding box), format: \"minX minY maxX maxY\"", cxxopts::value()) + ("i,invert", "Invert SDF") + ("o,output", "Output file in cdf format", cxxopts::value()->default_value("")) + ("input", "OBJ file containing input triangle mesh", cxxopts::value>()) + ; + + try + { + options.parse_positional("input"); + auto result = options.parse(argc, argv); + + if (result.count("help")) + { + std::cout << options.help() << std::endl; + std::cout << std::endl << std::endl << "Example: GenerateSDF2D -r \"50 50\" dragon.obj" << std::endl; + exit(0); + } + if (!result.count("input")) + { + std::cout << "ERROR: No input mesh given." << std::endl; + std::cout << options.help() << std::endl; + std::cout << std::endl << std::endl << "Example: GenerateSDF2D -r \"50 50\" dragon.obj" << std::endl; + exit(1); + } + auto resolution = result["r"].as>(); + auto filename = result["input"].as>().front(); + + if (!std::ifstream(filename).good()) + { + std::cerr << "ERROR: Input file does not exist!" << std::endl; + exit(1); + } + + std::cout << "Load mesh..."; + Discregrid::TriangleMesh mesh(filename); + std::cout << "DONE" << std::endl; + + std::cout << "Weld vertices and perform backface culling..." << std::endl; + std::vector culled_vertices; + std::vector culled_triangles; + mesh.weldVerticesAndCullBackfaces2D(culled_vertices, culled_triangles); + std::cout << "N Vertices: " << culled_vertices.size() << "; N Triangles: " << culled_triangles.size() / 3 << std::endl; + std::cout << "DONE" << std::endl; + + std::cout << "Create polygon from the union of all culled triangles..." << std::endl; + Clipper2Lib::PathsD polygon(culled_triangles.size() / 3); + for (int i = 0; i < culled_triangles.size(); i += 3) + { + const auto v0 = culled_vertices[culled_triangles[i]]; + const auto v1 = culled_vertices[culled_triangles[i + 1]]; + const auto v2 = culled_vertices[culled_triangles[i + 2]]; + polygon.emplace_back(Clipper2Lib::MakePathD({v0.x(), v0.y(), v1.x(), v1.y(), v2.x(), v2.y()})); + } + polygon = Clipper2Lib::Union(polygon, Clipper2Lib::FillRule::NonZero, 8); + std::cout << "DONE" << std::endl; + + std::cout << "Set up polygon distance function..." << std::endl; + Discregrid::PolygonDistance pd(polygon); + std::cout << "DONE" << std::endl; + + // The following commented code can be used to test the signed distance function sampling and save the polygon + sample(s) to an SVG file. + // std::filesystem::path filePath(filename); + // std::string svg_name = filePath.filename().stem().string().append(".svg"); + // std::cout << "Writing polygon as SVG to " << svg_name << "..." << std::endl; + // Clipper2Lib::SvgWriter svg; + // svg.AddPaths(polygon, false, Clipper2Lib::FillRule::NonZero, 0x4080ff9C, 0xFF003300, 1.5, false); + // Eigen::Vector2d sample(0.0, 0.0); + // const Discregrid::Result2D signed_distance = pd.signed_distance(sample); + // Clipper2Lib::PathD signed_distance_path = Clipper2Lib::MakePathD({signed_distance.nearest_point.x(), signed_distance.nearest_point.y(), sample.x(), sample.y()}); + // svg.AddPath(signed_distance_path, true, Clipper2Lib::FillRule::NonZero, 0xFFFF00FF, 0xFFFF00FF, 3, false); + // svg.SaveToFile(svg_name, 1000, 1000, 10); + // std::cout << "DONE" << std::endl; + + Eigen::AlignedBox2d domain; + domain.setEmpty(); + if (result.count("d")) + { + domain = result["d"].as(); + } + if (domain.isEmpty()) + { + for (size_t i = 0; i < polygon.size(); i++) + { + + const auto &outline = polygon[i]; + for (auto const& x : outline) + { + domain.extend(Matrix{x.x, x.y}); + } + } + + domain.max() += 1.0e-3 * domain.diagonal().norm() * Vector2d::Ones(); + domain.min() -= 1.0e-3 * domain.diagonal().norm() * Vector2d::Ones(); + } + + Discregrid::CubicLagrangeDiscreteGrid2D sdf(domain, resolution); + auto func = Discregrid::DiscreteGrid2D::ContinuousFunction{}; + if (result.count("invert")) + { + func = [&pd](Vector2d const& xi) { return -1.0 * pd.signed_distance(xi).distance; }; + } + else + { + func = [&pd](Vector2d const& xi) { return pd.signed_distance(xi).distance; }; + } + + std::cout << "Generate discretization..." << std::endl; + sdf.addFunction(func, true); + std::cout << "DONE" << std::endl; + + std::cout << "Serialize discretization..."; + auto output_file = result["o"].as(); + if (output_file == "") + { + output_file = filename; + if (output_file.find(".") != std::string::npos) + { + auto lastindex = output_file.find_last_of("."); + output_file = output_file.substr(0, lastindex); + } + output_file += ".cdf2d"; + } + sdf.save(output_file); + std::cout << "DONE" << std::endl; + } + catch (cxxopts::OptionException const& e) + { + std::cout << "error parsing options: " << e.what() << std::endl; + exit(1); + } + + return 0; +} \ No newline at end of file diff --git a/discregrid/CMakeLists.txt b/discregrid/CMakeLists.txt index 29f7476..81e3b1f 100755 --- a/discregrid/CMakeLists.txt +++ b/discregrid/CMakeLists.txt @@ -2,6 +2,8 @@ set(HEADERS include/Discregrid/discrete_grid.hpp include/Discregrid/cubic_lagrange_discrete_grid.hpp + include/Discregrid/discrete_grid2D.hpp + include/Discregrid/cubic_lagrange_discrete_grid2D.hpp ) set(HEADERS_MESH @@ -13,11 +15,13 @@ set(HEADERS_MESH set(HEADERS_GEOMETRY include/Discregrid/geometry/TriangleMeshDistance.h + include/Discregrid/geometry/PolygonDistance.hpp ) set(HEADERS_UTILITY include/Discregrid/utility/serialize.hpp include/Discregrid/utility/lru_cache.hpp + include/Discregrid/utility/bmp_file.hpp src/utility/timing.hpp src/utility/spinlock.hpp @@ -26,6 +30,8 @@ set(HEADERS_UTILITY set(SOURCES src/discrete_grid.cpp src/cubic_lagrange_discrete_grid.cpp + src/discrete_grid2D.cpp + src/cubic_lagrange_discrete_grid2D.cpp ) set(SOURCES_MESH @@ -36,6 +42,7 @@ set(SOURCES_MESH set(SOURCES_UTILITY src/utility/timing.cpp + src/utility/bmp_file.cpp ) macro(SOURCEGROUP name) @@ -81,6 +88,11 @@ if(WIN32) add_definitions(-D_USE_MATH_DEFINES) endif(WIN32) +# Enable PolygonDistance to use clipper2 reliant constructor. +if (USE_CLIPPER2) + add_definitions(-DIS_CLIPPER_ENABLED) +endif() + set(DISCREGRID_SOURCE_FILES ${HEADERS} ${SOURCES} @@ -104,6 +116,7 @@ endif() # Set link libraries. target_link_libraries(Discregrid + $,Clipper2,> ) install(TARGETS Discregrid diff --git a/discregrid/include/Discregrid/All b/discregrid/include/Discregrid/All index 34431bf..23faa6e 100755 --- a/discregrid/include/Discregrid/All +++ b/discregrid/include/Discregrid/All @@ -1,3 +1,5 @@ #include "cubic_lagrange_discrete_grid.hpp" +#include "cubic_lagrange_discrete_grid2D.hpp" #include "geometry/TriangleMeshDistance.h" +#include "geometry/PolygonDistance.hpp" #include "mesh/triangle_mesh.hpp" diff --git a/discregrid/include/Discregrid/cubic_lagrange_discrete_grid2D.hpp b/discregrid/include/Discregrid/cubic_lagrange_discrete_grid2D.hpp new file mode 100644 index 0000000..f7f17b5 --- /dev/null +++ b/discregrid/include/Discregrid/cubic_lagrange_discrete_grid2D.hpp @@ -0,0 +1,71 @@ +#pragma once + +#include "discrete_grid2D.hpp" + +namespace Discregrid +{ + +class CubicLagrangeDiscreteGrid2D : public DiscreteGrid2D +{ +public: + + CubicLagrangeDiscreteGrid2D(std::string const& filename); + CubicLagrangeDiscreteGrid2D(Eigen::AlignedBox2d const& domain, + std::array const& resolution); + + void save(std::string const& filename) const override; + void load(std::string const& filename) override; + + unsigned int addFunction(ContinuousFunction const& func, bool verbose = false, + SamplePredicate const& pred = nullptr) override; + + + std::size_t nCells() const { return m_n_cells; }; + double interpolate(unsigned int field_id, Eigen::Vector2d const& xi, + Eigen::Vector2d* gradient = nullptr) const override; + + /** + * @brief Determines the shape functions for the discretization with ID field_id at point xi. + * + * @param field_id Discretization ID + * @param x Location where the shape functions should be determined + * @param cell cell of x + * @param c0 vector required for the interpolation + * @param N shape functions for the cell of x + * @param dN (Optional) derivatives of the shape functions, required to compute the gradient + * @return Success of the function. + */ + bool determineShapeFunctions(unsigned int field_id, Eigen::Vector2d const &x, + std::array &cell, Eigen::Vector2d &c0, Eigen::Matrix &N, + Eigen::Matrix *dN = nullptr) const override; + + /** + * @brief Evaluates the given discretization with ID field_id at point xi. + * + * @param field_id Discretization ID + * @param xi Location where the discrete function is evaluated + * @param cell cell of xi + * @param c0 vector required for the interpolation + * @param N shape functions for the cell of xi + * @param gradient (Optional) if a pointer to a vector is passed the gradient of the discrete function will be evaluated + * @param dN (Optional) derivatives of the shape functions, required to compute the gradient + * @return double Results of the evaluation of the discrete function at point xi + */ + double interpolate(unsigned int field_id, Eigen::Vector2d const& xi, const std::array &cell, const Eigen::Vector2d &c0, const Eigen::Matrix &N, + Eigen::Vector2d* gradient = nullptr, Eigen::Matrix *dN = nullptr) const override; + + void reduceField(unsigned int field_id, Predicate pred) override; + + void forEachCell(unsigned int field_id, + std::function const& cb) const; + +private: + + Eigen::Vector2d indexToNodePosition(unsigned int l) const; + + std::vector> m_nodes; + std::vector>> m_cells; + std::vector> m_cell_map; +}; + +} \ No newline at end of file diff --git a/discregrid/include/Discregrid/discrete_grid2D.hpp b/discregrid/include/Discregrid/discrete_grid2D.hpp new file mode 100644 index 0000000..41ab9ac --- /dev/null +++ b/discregrid/include/Discregrid/discrete_grid2D.hpp @@ -0,0 +1,100 @@ +#pragma once + +#include +#include +#include +#include + +namespace Discregrid +{ + +class DiscreteGrid2D +{ +public: + + using CoefficientVector = Eigen::Matrix; + using ContinuousFunction = std::function; + using MultiIndex = std::array; + using Predicate = std::function; + using SamplePredicate = std::function; + + DiscreteGrid2D() = default; + DiscreteGrid2D(Eigen::AlignedBox2d const& domain, std::array const& resolution) + : m_domain(domain), m_resolution(resolution), m_n_fields(0u) + { + auto n = Eigen::Matrix::Map(resolution.data()); + m_cell_size = domain.diagonal().cwiseQuotient(n.cast()); + m_inv_cell_size = m_cell_size.cwiseInverse(); + m_n_cells = n.prod(); + } + virtual ~DiscreteGrid2D() = default; + + virtual void save(std::string const& filename) const = 0; + virtual void load(std::string const& filename) = 0; + + virtual unsigned int addFunction(ContinuousFunction const& func, bool verbose = false, + SamplePredicate const& pred = nullptr) = 0; + + double interpolate(Eigen::Vector2d const& xi, Eigen::Vector2d* gradient = nullptr) const + { + return interpolate(0u, xi, gradient); + } + + virtual double interpolate(unsigned int field_id, Eigen::Vector2d const& xi, + Eigen::Vector2d* gradient = nullptr) const = 0; + + /** + * @brief Determines the shape functions for the discretization with ID field_id at point xi. + * + * @param field_id Discretization ID + * @param x Location where the shape functions should be determined + * @param cell cell of x + * @param c0 vector required for the interpolation + * @param N shape functions for the cell of x + * @param dN (Optional) derivatives of the shape functions, required to compute the gradient + * @return Success of the function. + */ + virtual bool determineShapeFunctions(unsigned int field_id, Eigen::Vector2d const &x, + std::array &cell, Eigen::Vector2d &c0, Eigen::Matrix &N, + Eigen::Matrix *dN = nullptr) const = 0; + + /** + * @brief Evaluates the given discretization with ID field_id at point xi. + * + * @param field_id Discretization ID + * @param xi Location where the discrete function is evaluated + * @param cell cell of xi + * @param c0 vector required for the interpolation + * @param N shape functions for the cell of xi + * @param gradient (Optional) if a pointer to a vector is passed the gradient of the discrete function will be evaluated + * @param dN (Optional) derivatives of the shape functions, required to compute the gradient + * @return double Results of the evaluation of the discrete function at point xi + */ + virtual double interpolate(unsigned int field_id, Eigen::Vector2d const& xi, const std::array &cell, const Eigen::Vector2d &c0, const Eigen::Matrix &N, + Eigen::Vector2d* gradient = nullptr, Eigen::Matrix *dN = nullptr) const = 0; + + virtual void reduceField(unsigned int field_id, Predicate pred) {} + + + MultiIndex singleToMultiIndex(unsigned int i) const; + unsigned int multiToSingleIndex(MultiIndex const& ij) const; + + Eigen::AlignedBox2d subdomain(MultiIndex const& ij) const; + Eigen::AlignedBox2d subdomain(unsigned int l) const; + + Eigen::AlignedBox2d const& domain() const { return m_domain; } + std::array const& resolution() const { return m_resolution; } + Eigen::Vector2d const& cellSize() const { return m_cell_size;} + Eigen::Vector2d const& invCellSize() const { return m_inv_cell_size;} + +protected: + + + Eigen::AlignedBox2d m_domain; + std::array m_resolution; + Eigen::Vector2d m_cell_size; + Eigen::Vector2d m_inv_cell_size; + std::size_t m_n_cells; + std::size_t m_n_fields; +}; +} \ No newline at end of file diff --git a/discregrid/include/Discregrid/geometry/PolygonDistance.hpp b/discregrid/include/Discregrid/geometry/PolygonDistance.hpp new file mode 100644 index 0000000..9486ea7 --- /dev/null +++ b/discregrid/include/Discregrid/geometry/PolygonDistance.hpp @@ -0,0 +1,534 @@ +#pragma once +#include +#include +#include +#include +#include "TriangleMeshDistance.h" +#ifdef IS_CLIPPER_ENABLED +#include "clipper2/clipper.core.h" +#endif + +// This file is a derivative of TriangleMeshDistance.h providing the same functionality in 2D. +// This file can be used to get the distance of a point to a closed polygon (that can contain holes). +namespace Discregrid +{ + // Point-Edge distance definitions + enum class NearestEntity2D { V0, V1, E }; + static double point_edge_sq_unsigned( + NearestEntity2D& nearest_entity, + Eigen::Vector2d& nearest_point, + const Eigen::Vector2d& point, + const Eigen::Vector2d& v0, + const Eigen::Vector2d& v1); + // ----------------------------------- + + // Struct that contains the result of a 2D distance query + struct Result2D + { + double distance = std::numeric_limits::max(); + Eigen::Vector2d nearest_point; + NearestEntity2D nearest_entity; + int edge_id = -1; + }; + // ----------------------------------- + + // A class to compute signed and unsigned distances to a triangle mesh. + class PolygonDistance + { + private: + /* Definitions */ + struct BoundingSphere + { + Eigen::Vector2d center; + double radius; + }; + + struct Node + { + BoundingSphere bv_left; + BoundingSphere bv_right; + int left = -1; // If left == -1, right is the edge_id + int right = -1; + }; + + struct Edge + { + std::array vertices; + int id = -1; + }; + + /* Fields */ + std::vector vertices; + std::vector edges; + std::vector nodes; + std::vector pseudonormals_edges; + std::vector pseudonormals_vertices; + BoundingSphere root_bv; + bool is_constructed = false; + + /* Methods */ + void _construct(); + void _build_tree(const int node_id, BoundingSphere& bounding_sphere, std::vector &edges, const int begin, const int end); + void _query(Result2D &result, const Node &node, const Eigen::Vector2d& point) const; + + public: + /* Methods */ + PolygonDistance() = default; + + /** + * @brief Constructs a new PolygonDistance object. + * + * @param vertices Pointer to the vertices coordinates array in xyxy... layout. + * @param n_vertices Number of vertices. + * @param edges Pointer to the connectivity array in ijij... layout. + * @param n_edges Number of edges. + */ + template + PolygonDistance(const FLOAT* vertices, const SIZE_T n_vertices, const INT* edges, const SIZE_T n_edges); + + /** + * @brief Constructs a new PolygonDistance object. + * + * @param vertices Vertices of the polygon. Y coordinate of the 3rd vertex should be accessed by `vertices[2][1]`. + * @param edges Edges of the polygon. Index of the 2nd vertex of the 3rd edge should be accessed by `edges[2][1]`. + */ + template + PolygonDistance(const std::vector& vertices, const std::vector& edges); + +#ifdef IS_CLIPPER_ENABLED + /** + * @brief Constructs a new PolygonDistance object. + * + * @param polygon Clipper2 PathsD (= Polygon) object. + */ + PolygonDistance(const Clipper2Lib::PathsD &polygon); +#endif + + /** + * @brief Initializes an existing PolygonDistance object (including empty ones). + * + * @param vertices Pointer to the vertices coordinates array in xyxy... layout. + * @param n_vertices Number of vertices. + * @param edges Pointer to the conectivity array in ijij... layout. + * @param n_edges Number of edges. + */ + template + void construct(const FLOAT* vertices, const SIZE_T n_vertices, const INT* edges, const SIZE_T n_edges); + + /** + * @brief Initializes an existing PolygonDistance object (including empty ones). + * + * @param vertices Vertices of the polygon. Y coordinate of the 3rd vertex should be access by `vertices[2][1]`. + * @param edges Edges of the polygon. Index of the 2nd vertex of the 3rd edge should be access by `edges[2][1]`. + */ + template + void construct(const std::vector& vertices, const std::vector& edges); + + /** + * @brief Computes the unsigned distance from a point to the polygon. Thread safe. + * + * @param point to query from. Typed to `Vector2d` but can be passed as `{x, y}`. + * + * @return Result containing distance, nearest point on the polygon, nearest entity and the nearest edge index. + */ + template + Result2D unsigned_distance(const IndexableVector2double& point) const; + Result2D unsigned_distance(const std::array& point) const; + + /** + * @brief Computes the unsigned distance from a point to the polygon. Thread safe. + * + * @param point to query from. Typed to `Vector2d` but can be passed as `{x, y}`. + * + * @return Result containing distance, nearest point on the polygon, nearest entity and the nearest edge index. + */ + template + Result2D signed_distance(const IndexableVector2double& point) const; + Result2D signed_distance(const std::array& point) const; + }; +} + +/* ========================================== DECLARATIONS ========================================== */ +template +inline Discregrid::PolygonDistance::PolygonDistance(const FLOAT* vertices, const SIZE_T n_vertices, const INT* edges, const SIZE_T n_edges) +{ + this->construct(vertices, n_vertices, edges, n_edges); +} + +template +inline Discregrid::PolygonDistance::PolygonDistance(const std::vector& vertices, const std::vector& edges) +{ + this->construct(vertices, edges); +} + +#ifdef IS_CLIPPER_ENABLED +inline Discregrid::PolygonDistance::PolygonDistance(const Clipper2Lib::PathsD &polygon) +{ + size_t size = 0; + for (size_t i = 0; i < polygon.size(); i++) + { + size += polygon[i].size(); + } + std::vector polygonVertices; + polygonVertices.reserve(size); + std::vector polygonEdges; + polygonEdges.reserve(size); + size_t counter = 0; + for (size_t i = 0; i < polygon.size(); i++) + { + auto path = polygon[i]; + size_t subCounter = counter; + for (size_t j = 0; j < path.size() - 1; j++) + { + auto point = path[j]; + polygonVertices.emplace_back(point.x, point.y); + polygonEdges.emplace_back(counter++, counter); + } + auto point = path.back(); + polygonVertices.emplace_back(point.x, point.y); + polygonEdges.emplace_back(counter++, subCounter); + } + this->construct(polygonVertices, polygonEdges); +} +#endif + +template +inline void Discregrid::PolygonDistance::construct(const FLOAT* vertices, const SIZE_T n_vertices, const INT* edges, const SIZE_T n_edges) +{ + this->vertices.resize(2 * n_vertices); + for (size_t i = 0; i < static_cast(n_vertices); i++) + { + this->vertices[i][0] = static_cast(vertices[2 * i]); + this->vertices[i][1] = static_cast(vertices[2 * i + 1]); + } + + this->edges.resize(2 * n_edges); + for (size_t i = 0; i < static_cast(n_edges); i++) + { + this->edges[i][0] = static_cast(edges[2 * i]); + this->edges[i][1] = static_cast(edges[2 * i + 1]); + } + this->_construct(); +} + +template +inline void Discregrid::PolygonDistance::construct(const std::vector& vertices, const std::vector& edges) +{ + this->vertices.resize(vertices.size()); + for (size_t i = 0; i < vertices.size(); i++) + { + this->vertices[i][0] = static_cast(vertices[i][0]); + this->vertices[i][1] = static_cast(vertices[i][1]); + } + + this->edges.resize(edges.size()); + for (size_t i = 0; i < edges.size(); i++) + { + this->edges[i][0] = static_cast(edges[i][0]); + this->edges[i][1] = static_cast(edges[i][1]); + } + this->_construct(); +} + +template +inline Discregrid::Result2D Discregrid::PolygonDistance::unsigned_distance(const IndexableVector2double& point) const +{ + return this->unsigned_distance({static_cast(point[0]), static_cast(point[1])}); +} + +inline Discregrid::Result2D Discregrid::PolygonDistance::unsigned_distance(const std::array& point) const +{ + if (!this->is_constructed) + { + std::cout << "PolygonDistance error: not constructed." << std::endl; + exit(-1); + } + + const Eigen::Vector2d p(point[0], point[1]); + Result2D result; + result.distance = std::numeric_limits::max(); + this->_query(result, this->nodes[0], p); + return result; +} + +template +inline Discregrid::Result2D Discregrid::PolygonDistance::signed_distance(const IndexableVector2double& point) const +{ + const Eigen::Vector2d p(point[0], point[1]); + Result2D result = this->unsigned_distance(p); + + const Eigen::Vector2i& edge = this->edges[result.edge_id]; + Eigen::Vector2d pseudonormal; + switch (result.nearest_entity) + { + case NearestEntity2D::V0: + pseudonormal = this->pseudonormals_vertices[edge[0]]; + break; + case NearestEntity2D::V1: + pseudonormal = this->pseudonormals_vertices[edge[1]]; + break; + case NearestEntity2D::E: + pseudonormal = this->pseudonormals_edges[result.edge_id]; + break; + default: + break; + } + + const Eigen::Vector2d u = p - result.nearest_point; + result.distance *= u.dot(pseudonormal) >= 0.0 ? 1.0 : -1.0; + + return result; +} + +inline Discregrid::Result2D Discregrid::PolygonDistance::signed_distance(const std::array& point) const +{ + return this->signed_distance({static_cast(point[0]), static_cast(point[1])}); +} + +inline void Discregrid::PolygonDistance::_construct() +{ + if (this->edges.empty()) + { + std::cout << "PolygonDistance error: Empty edge list." << std::endl; + exit(-1); + } + + // Build the tree containing the edges + std::vector edges; + + edges.resize(this->edges.size()); + for (int i = 0; i < static_cast(this->edges.size()); i++) + { + edges[i].id = i; + + const Eigen::Vector2i& edge = this->edges[i]; + edges[i].vertices[0] = this->vertices[edge[0]]; + edges[i].vertices[1] = this->vertices[edge[1]]; + } + + this->nodes.emplace_back(); + this->_build_tree(0, this->root_bv, edges, 0, static_cast(edges.size())); + + // Compute + this->pseudonormals_edges.reserve(this->edges.size()); + this->pseudonormals_vertices.resize(this->vertices.size(), { 0, 0 }); + + std::unordered_map> vertexIndexToEdgeNormalIndices; + + // (Pseudo-)Normals edges + for (int i = 0; i < static_cast(this->edges.size()); i++) { + + // Edge + const Eigen::Vector2i edge = this->edges[i]; + const int i0 = edge[0]; + const int i1 = edge[1]; + const Eigen::Vector2d a = this->vertices[i0]; + const Eigen::Vector2d b = this->vertices[i1]; + + // Edge normal calculation routine from Clipper2 library + if (a == b) + { + std::cout << "PolygonDistance::_construct error: Can not calculate edge normal because both vertices are the same!" << std::endl; + exit(-1); + } + auto dx = b.x() - a.x(); + auto dy = b.y() - a.y(); + const auto inverse_hypot = 1.0 / std::sqrt(dx * dx + dy * dy); + dx *= inverse_hypot; + dy *= inverse_hypot; + const Eigen::Vector2d normal(dy, -dx); + + for (const int index : edge) + { + auto edgeNormalIndex = static_cast(this->pseudonormals_edges.size()); + + if (vertexIndexToEdgeNormalIndices.find(index) == vertexIndexToEdgeNormalIndices.end()) + { + vertexIndexToEdgeNormalIndices[index] = std::vector(); + } + + vertexIndexToEdgeNormalIndices[index].push_back(edgeNormalIndex); + } + + this->pseudonormals_edges.push_back(normal); + } + + // Pseudonormals vertices + for (const auto& pair : vertexIndexToEdgeNormalIndices) + { + int vertexIndex = pair.first; + std::vector edgeNormalIndices = pair.second; + if (edgeNormalIndices.size() != 2) + { + std::cout << "PolygonDistance error: Polygon is not valid. At least one vertex does not have exactly two incident edges." << std::endl; + exit(-1); + } + this->pseudonormals_vertices[vertexIndex] = ((this->pseudonormals_edges[edgeNormalIndices[0]] + this->pseudonormals_edges[edgeNormalIndices[1]]) / 2.0).normalized(); + } + + this->is_constructed = true; +} + +inline void Discregrid::PolygonDistance::_build_tree(const int node_id, BoundingSphere& bounding_sphere, std::vector &edges, const int begin, const int end) +{ + const int n_edges = end - begin; + + if (n_edges <= 0) + { + std::cout << "PolygonDistance::_construct error: Empty leaf." << std::endl; + exit(-1); + } + + if (n_edges == 1) + { + // Build node leaf + this->nodes[node_id].left = -1; + this->nodes[node_id].right = edges[begin].id; + + // Bounding sphere + const Edge& edge = edges[begin]; + const Eigen::Vector2d helper = edge.vertices[0] + edge.vertices[1]; + const Eigen::Vector2d center(helper.x() / 2.0, helper.y() / 2.0); + const double radius = std::max((edge.vertices[0] - center).norm(), (edge.vertices[1] - center).norm()); + bounding_sphere.center = center; + bounding_sphere.radius = radius; + } + else + { + // Compute AABB center and largest dimension of all current edges. (AABB center is the average of all edge vertex positions) + Eigen::Vector2d top = {std::numeric_limits::lowest(), std::numeric_limits::lowest()}; + Eigen::Vector2d bottom = {std::numeric_limits::max(), std::numeric_limits::max()}; + Eigen::Vector2d center = {0, 0}; + for (int edge_i = begin; edge_i < end; edge_i++) + { + for (int vertex_i = 0; vertex_i < 2; vertex_i++) + { + const Eigen::Vector2d& p = edges[edge_i].vertices[vertex_i]; + center += p; + + for (int coord_i = 0; coord_i < 2; coord_i++) + { + top[coord_i] = std::max(top[coord_i], p[coord_i]); + bottom[coord_i] = std::min(bottom[coord_i], p[coord_i]); + } + } + } + center /= 2 * n_edges; + const Eigen::Vector2d diagonal = top - bottom; + const int split_dim = diagonal[0] < diagonal[1] ? 1 : 0; + + // Set node bounding sphere + double radius_sq = 0.0; + for (int edge_i = begin; edge_i < end; edge_i++) + { + for (int i = 0; i < 2; i++) + { + radius_sq = std::max(radius_sq, (center - edges[edge_i].vertices[i]).squaredNorm()); + } + } + bounding_sphere.center = center; + bounding_sphere.radius = std::sqrt(radius_sq); + + // Sort the triangles according to their center along the split dimension + std::sort(edges.begin() + begin, edges.begin() + end, + [split_dim](const Edge& a, const Edge& b) + { + return a.vertices[0][split_dim] < b.vertices[0][split_dim]; + } + ); + + // Children + const int mid = static_cast(0.5 * (begin + end)); + + this->nodes[node_id].left = static_cast(this->nodes.size()); + this->nodes.emplace_back(); + this->_build_tree(this->nodes[node_id].left, this->nodes[node_id].bv_left, edges, begin, mid); + + this->nodes[node_id].right = static_cast(this->nodes.size()); + this->nodes.emplace_back(); + this->_build_tree(this->nodes[node_id].right, this->nodes[node_id].bv_right, edges, mid, end); + } +} + + +inline void Discregrid::PolygonDistance::_query(Result2D &result, const Node &node, const Eigen::Vector2d& point) const +{ + // End of recursion + if (node.left == -1) + { + const int edge_id = node.right; + const Eigen::Vector2i& edge = this->edges[node.right]; // If left == -1, right is the edge + const Eigen::Vector2d& v0 = this->vertices[edge[0]]; + const Eigen::Vector2d& v1 = this->vertices[edge[1]]; + + Eigen::Vector2d nearest_point; + NearestEntity2D nearest_entity; + const double distance_sq = point_edge_sq_unsigned(nearest_entity, nearest_point, point, v0, v1); + + if (distance_sq < result.distance * result.distance) + { + result.nearest_point = nearest_point; + result.nearest_entity = nearest_entity; + result.distance = std::sqrt(distance_sq); + result.edge_id = edge_id; + } + } + // Recursion + else + { + // Find which child bounding volume is closer + const double d_left = (point - node.bv_left.center).norm() - node.bv_left.radius; + const double d_right = (point - node.bv_right.center).norm() - node.bv_right.radius; + + if (d_left < d_right) + { + // Overlap test + if (d_left < result.distance) + { + this->_query(result, this->nodes[node.left], point); + } + + if (d_right < result.distance) + { + this->_query(result, this->nodes[node.right], point); + } + } + else + { + if (d_right < result.distance) + { + this->_query(result, this->nodes[node.right], point); + } + if (d_left < result.distance) + { + this->_query(result, this->nodes[node.left], point); + } + } + } +} + +static double Discregrid::point_edge_sq_unsigned(NearestEntity2D& nearest_entity, Eigen::Vector2d& nearest_point, const Eigen::Vector2d& point, const Eigen::Vector2d& v0, const Eigen::Vector2d& v1) +{ + const Eigen::Vector2d e = v1 - v0; + const double e2 = e.dot(e); + const Eigen::Vector2d diff = point - v0; + const double t = std::clamp(diff.dot(e) / e2, 0.0, 1.0); + + // Closest point on edge + nearest_point = {v0.x() + t * e.x(), v0.y() + t * e.y()}; + + if (t <= std::numeric_limits::epsilon()) + { + nearest_entity = NearestEntity2D::V0; + } + else if (t >= 1.0 - std::numeric_limits::epsilon()) + { + nearest_entity = NearestEntity2D::V1; + } + else + { + nearest_entity = NearestEntity2D::E; + } + + return (point-nearest_point).squaredNorm(); +} \ No newline at end of file diff --git a/discregrid/include/Discregrid/mesh/triangle_mesh.hpp b/discregrid/include/Discregrid/mesh/triangle_mesh.hpp index 66d38b5..d931b0e 100755 --- a/discregrid/include/Discregrid/mesh/triangle_mesh.hpp +++ b/discregrid/include/Discregrid/mesh/triangle_mesh.hpp @@ -92,6 +92,8 @@ class TriangleMesh Eigen::Vector3d computeFaceNormal(unsigned int f) const; + void weldVerticesAndCullBackfaces2D(std::vector& culled_vertices, std::vector& culled_triangles) const; + private: void construct(); diff --git a/cmd/discrete_field_to_bitmap/bmp_file.hpp b/discregrid/include/Discregrid/utility/bmp_file.hpp old mode 100755 new mode 100644 similarity index 96% rename from cmd/discrete_field_to_bitmap/bmp_file.hpp rename to discregrid/include/Discregrid/utility/bmp_file.hpp index a3204a3..622ab00 --- a/cmd/discrete_field_to_bitmap/bmp_file.hpp +++ b/discregrid/include/Discregrid/utility/bmp_file.hpp @@ -1,42 +1,42 @@ -#pragma once - -class BmpReaderWriter -{ -public: - static bool isBigEndian(); - static unsigned short endianSwap(unsigned short nValue); - static unsigned int endianSwap(unsigned int i); - - // ------------------------------------------------------------------- - -#pragma pack(1) - - struct BMPHEADER { - unsigned short Type; - unsigned int Size; - unsigned short Reserved1; - unsigned short Reserved2; - unsigned int OffBits; - }; - - // Only Win3.0 BMPINFO (see later for OS/2) - struct BMPINFO { - unsigned int Size; - unsigned int Width; - unsigned int Height; - unsigned short Planes; - unsigned short BitCount; - unsigned int Compression; - unsigned int SizeImage; - unsigned int XPelsPerMeter; - unsigned int YPelsPerMeter; - unsigned int ClrUsed; - unsigned int ClrImportant; - }; - -#pragma pack() - - // Data is persists until the class is destructed. - static bool loadFile(const char *filename, unsigned int &width, unsigned int &height, unsigned char *&data); - static bool saveFile(const char *filename, int width, int height, unsigned char *data); -}; +#pragma once + +class BmpReaderWriter +{ +public: + static bool isBigEndian(); + static unsigned short endianSwap(unsigned short nValue); + static unsigned int endianSwap(unsigned int i); + + // ------------------------------------------------------------------- + +#pragma pack(1) + + struct BMPHEADER { + unsigned short Type; + unsigned int Size; + unsigned short Reserved1; + unsigned short Reserved2; + unsigned int OffBits; + }; + + // Only Win3.0 BMPINFO (see later for OS/2) + struct BMPINFO { + unsigned int Size; + unsigned int Width; + unsigned int Height; + unsigned short Planes; + unsigned short BitCount; + unsigned int Compression; + unsigned int SizeImage; + unsigned int XPelsPerMeter; + unsigned int YPelsPerMeter; + unsigned int ClrUsed; + unsigned int ClrImportant; + }; + +#pragma pack() + + // Data is persists until the class is destructed. + static bool loadFile(const char *filename, unsigned int &width, unsigned int &height, unsigned char *&data); + static bool saveFile(const char *filename, int width, int height, unsigned char *data); +}; diff --git a/discregrid/src/cubic_lagrange_discrete_grid2D.cpp b/discregrid/src/cubic_lagrange_discrete_grid2D.cpp new file mode 100644 index 0000000..1689b13 --- /dev/null +++ b/discregrid/src/cubic_lagrange_discrete_grid2D.cpp @@ -0,0 +1,669 @@ +#include "cubic_lagrange_discrete_grid2D.hpp" + +#include "data/z_sort_table.hpp" +#include +#include "utility/spinlock.hpp" + +#include +#include +#include +#include +#include +#include +#include + +using namespace Eigen; + +namespace Discregrid +{ + +namespace +{ + +// See Zienkiewicz The Finite Element Method its Basis and Fundamentals 2013 page 159 (2D) and 168 (3D) +Matrix +shape_function_(Vector2d const &xi, Matrix *gradient = nullptr) +{ + auto res = Matrix{}; + + const auto x = xi[0]; + const auto y = xi[1]; + + const auto x2 = x * x; + const auto y2 = y * y; + + const auto _1mx = 1.0 - x; + const auto _1my = 1.0 - y; + + const auto _1px = 1.0 + x; + const auto _1py = 1.0 + y; + + const auto _1m3x = 1.0 - 3.0 * x; + const auto _1m3y = 1.0 - 3.0 * y; + + const auto _1p3x = 1.0 + 3.0 * x; + const auto _1p3y = 1.0 + 3.0 * y; + + const auto _1mxt1my = _1mx * _1my; + const auto _1mxt1py = _1mx * _1py; + const auto _1pxt1my = _1px * _1my; + const auto _1pxt1py = _1px * _1py; + + const auto _1mx2 = 1.0 - x2; + const auto _1my2 = 1.0 - y2; + + // Corner nodes. + auto factor = (1.0 / 32.0) * (9.0 * (x2 + y2) - 10.0); + res[0] = factor * _1mxt1my; + res[1] = factor * _1pxt1my; + res[2] = factor * _1mxt1py; + res[3] = factor * _1pxt1py; + + // Edge nodes. + factor = (9.0 / 32.0) * _1mx2; + auto factort1m3x = factor * _1m3x; // Note that x is multiplied with xi (+-1/3) so 9 * xi * x becomes 1/3 * x + auto factort1p3x = factor * _1p3x; + res[4] = factort1m3x * _1my; + res[5] = factort1p3x * _1my; + res[6] = factort1m3x * _1py; + res[7] = factort1p3x * _1py; + + factor = (9.0 / 32.0) * _1my2; + auto factort1m3y = factor * _1m3y; + auto factort1p3y = factor * _1p3y; + res[8] = factort1m3y * _1mx; + res[9] = factort1p3y * _1mx; + res[10] = factort1m3y * _1px; + res[11] = factort1p3y * _1px; + + if (gradient) + { + auto &dN = *gradient; + + const auto _27tx2p9ty2m10 = 27 * x2 + 9 * y2 - 10; + const auto _27ty2p9tx2m10 = 27 * y2 + 9 * x2 - 10; + + const auto _18x = 18 * x; + const auto _18y = 18 * y; + + dN(0, 0) = _1my * (_18x - _27tx2p9ty2m10); + dN(0, 1) = _1mx * (_18y - _27ty2p9tx2m10); + dN(1, 0) = _1my * (_18x + _27tx2p9ty2m10); + dN(1, 1) = _1px * (_18y - _27ty2p9tx2m10); + dN(2, 0) = _1py * (_18x - _27tx2p9ty2m10); + dN(2, 1) = _1mx * (_18y + _27ty2p9tx2m10); + dN(3, 0) = _1py * (_18x + _27tx2p9ty2m10); + dN(3, 1) = _1px * (_18y + _27ty2p9tx2m10); + + dN.topRows(4) /= 32.0; + + const auto _9x2 = 9 * x2; + const auto _2x = 2 * x; + const auto _3x = 3 * x; + const auto _3x3 = 3 * x2 * x; + const auto _9x2m2xm3 = _9x2 - _2x - 3; + const auto _3x3mx2m3xp1 = _3x3 - x2 - _3x + 1; + const auto _9x2p2xm3 = _9x2 + _2x - 3; + const auto _3x3px2m3xm1 = _3x3 + x2 - _3x - 1; + + dN(4, 0) = _1my * _9x2m2xm3; + dN(4, 1) = -_3x3mx2m3xp1; + dN(5, 0) = _1my * -_9x2p2xm3; + dN(5, 1) = _3x3px2m3xm1; + dN(6, 0) = _1py * _9x2m2xm3; + dN(6, 1) = _3x3mx2m3xp1; + dN(7, 0) = _1py * -_9x2p2xm3; + dN(7, 1) = -_3x3px2m3xm1; + + const auto _9y2 = 9 * y2; + const auto _2y = 2 * y; + const auto _3y = 3 * y; + const auto _3y3 = 3 * y2 * y; + const auto _3y3my2m3yp1 = _3y3 - y2 - _3y + 1; + const auto _9y2m2ym3 = _9y2 - _2y - 3; + const auto _3y3py2m3ym1 = _3y3 + y2 - _3y - 1; + const auto _9y2p2ym3 = _9y2 + _2y - 3; + + dN(8, 0) = -_3y3my2m3yp1; + dN(8, 1) = _1mx * _9y2m2ym3; + dN(9, 0) = _3y3py2m3ym1; + dN(9, 1) = _1mx * -_9y2p2ym3; + dN(10, 0) = _3y3my2m3yp1; + dN(10, 1) = _1px * _9y2m2ym3; + dN(11, 0) = -_3y3py2m3ym1; + dN(11, 1) = _1px * -_9y2p2ym3; + + dN.bottomRows(12u - 4u) *= 9.0 / 32.0; + } + + return res; +} + +// Determines Morten value according to z-curve. +inline uint64_t +zValue(Vector2d const &x, double invCellSize) +{ + std::array key; + for (unsigned int i(0); i < 2; ++i) + { + if (x[i] >= 0.0) + key[i] = static_cast(invCellSize * x[i]); + else + key[i] = static_cast(invCellSize * x[i]) - 1; + } + + std::array p = { + static_cast(static_cast(key[0]) - (std::numeric_limits::lowest() + 1)), + static_cast(static_cast(key[1]) - (std::numeric_limits::lowest() + 1))}; + + return morton_lut(p); +} +} // namespace + + Vector2d CubicLagrangeDiscreteGrid2D::indexToNodePosition(unsigned int l) const +{ + auto x = Vector2d{}; + + auto n = Matrix::Map(m_resolution.data()); + + // Amount of vertices and edges + auto nv = (n[0] + 1) * (n[1] + 1); + auto ne_x = n[0] * (n[1] + 1); + + auto ij = Matrix{}; + // 0 -> nv: l is a cell's corner vertex + if (l < nv) + { + ij(1) = l / (n[0] + 1); + ij(0) = l % (n[0] + 1); + + x = m_domain.min() + m_cell_size.cwiseProduct(ij.cast()); + } + // nv -> nv + 2 * ne_x: l is a non-corner vertex in x direction + else if (l < nv + 2 * ne_x) + { + l -= nv; + // Identify cell index, then split into x and y dimensions. + auto e_ind = l / 2; + ij(1) = e_ind / n[0]; + ij(0) = e_ind % n[0]; + // Calculate cell position + x = m_domain.min() + m_cell_size.cwiseProduct(ij.cast()); + // Add sub-cell node position to cell position: Either 1/3 or 2/3 of cell size in x direction + x(0) += (1.0 + static_cast(l % 2)) / 3.0 * m_cell_size[0]; + } + // nv + 2 * ne_x -> nv + 2 * (ne_x + ne_y): l is non-corner vertex in y direction + else + { + l -= nv + 2 * ne_x; + // Identify cell index, then split into x and y dimensions. + auto e_ind = l / 2; + ij(1) = e_ind % n[1]; + ij(0) = e_ind / n[1]; + + // Calculate cell position + x = m_domain.min() + m_cell_size.cwiseProduct(ij.cast()); + // Add sub-cell node position to cell position: Either 1/3 or 2/3 of cell size in y direction + x(1) += (1.0 + static_cast(l % 2)) / 3.0 * m_cell_size[1]; + } + + return x; +} + +CubicLagrangeDiscreteGrid2D::CubicLagrangeDiscreteGrid2D(std::string const &filename) +{ + load(filename); +} + +CubicLagrangeDiscreteGrid2D::CubicLagrangeDiscreteGrid2D(AlignedBox2d const &domain, std::array const &resolution) : DiscreteGrid2D(domain, resolution) +{ +} + +void CubicLagrangeDiscreteGrid2D::save(std::string const &filename) const +{ + auto out = std::ofstream(filename, std::ios::binary); + serialize::write(*out.rdbuf(), m_domain); + serialize::write(*out.rdbuf(), m_resolution); + serialize::write(*out.rdbuf(), m_cell_size); + serialize::write(*out.rdbuf(), m_inv_cell_size); + serialize::write(*out.rdbuf(), m_n_cells); + serialize::write(*out.rdbuf(), m_n_fields); + + serialize::write(*out.rdbuf(), m_nodes.size()); + for (auto const &nodes : m_nodes) + { + serialize::write(*out.rdbuf(), nodes.size()); + for (auto const &node : nodes) + { + serialize::write(*out.rdbuf(), node); + } + } + + serialize::write(*out.rdbuf(), m_cells.size()); + for (auto const &cells : m_cells) + { + serialize::write(*out.rdbuf(), cells.size()); + for (auto const &cell : cells) + { + serialize::write(*out.rdbuf(), cell); + } + } + + serialize::write(*out.rdbuf(), m_cell_map.size()); + for (auto const &maps : m_cell_map) + { + serialize::write(*out.rdbuf(), maps.size()); + for (auto const &map : maps) + { + serialize::write(*out.rdbuf(), map); + } + } + + out.close(); +} + +void CubicLagrangeDiscreteGrid2D::load(std::string const &filename) +{ + auto in = std::ifstream(filename, std::ios::binary); + + if (!in.good()) + { + std::cerr << "ERROR: Discrete grid can not be loaded. Input file does not exist!" << std::endl; + return; + } + + serialize::read(*in.rdbuf(), m_domain); + serialize::read(*in.rdbuf(), m_resolution); + serialize::read(*in.rdbuf(), m_cell_size); + serialize::read(*in.rdbuf(), m_inv_cell_size); + serialize::read(*in.rdbuf(), m_n_cells); + serialize::read(*in.rdbuf(), m_n_fields); + + auto n_nodes = std::size_t{}; + serialize::read(*in.rdbuf(), n_nodes); + m_nodes.resize(n_nodes); + for (auto &nodes : m_nodes) + { + serialize::read(*in.rdbuf(), n_nodes); + nodes.resize(n_nodes); + for (auto &node : nodes) + { + serialize::read(*in.rdbuf(), node); + } + } + + auto n_cells = std::size_t{}; + serialize::read(*in.rdbuf(), n_cells); + m_cells.resize(n_cells); + for (auto &cells : m_cells) + { + serialize::read(*in.rdbuf(), n_cells); + cells.resize(n_cells); + for (auto &cell : cells) + { + serialize::read(*in.rdbuf(), cell); + } + } + + auto n_cell_maps = std::size_t{}; + serialize::read(*in.rdbuf(), n_cell_maps); + m_cell_map.resize(n_cell_maps); + for (auto &cell_maps : m_cell_map) + { + serialize::read(*in.rdbuf(), n_cell_maps); + cell_maps.resize(n_cell_maps); + for (auto &cell_map : cell_maps) + { + serialize::read(*in.rdbuf(), cell_map); + } + } + + in.close(); +} + +unsigned int CubicLagrangeDiscreteGrid2D::addFunction(ContinuousFunction const &func, bool verbose, SamplePredicate const &pred) +{ + using namespace std::chrono; + + const auto t0_construction = high_resolution_clock::now(); + + auto n = Matrix::Map(m_resolution.data()); + + // Amount of vertices and edges + const auto nv = (n[0] + 1) * (n[1] + 1); + const auto ne_x = n[0] * (n[1] + 1); + const auto ne_y = (n[0] + 1) * n[1]; + const auto ne = ne_x + ne_y; + const auto n_nodes = nv + 2 * ne; + + m_nodes.emplace_back(); + auto &coeffs = m_nodes.back(); + coeffs.resize(n_nodes); + + std::atomic_uint counter(0u); + SpinLock mutex; + auto t0 = high_resolution_clock::now(); + + // Evaluate function at every node +#pragma omp parallel default(shared) + { +#pragma omp for schedule(static) nowait + for (int l = 0; l < static_cast(n_nodes); ++l) + { + auto x = indexToNodePosition(l); + auto &c = coeffs[l]; + + if (!pred || pred(x)) + c = func(x); + else + c = std::numeric_limits::max(); + + if (verbose && (++counter == n_nodes || duration_cast(high_resolution_clock::now() - t0).count() > 1000u)) + { + std::async(std::launch::async, [&] + { + mutex.lock(); + t0 = high_resolution_clock::now(); + std::cout << "\r" + << "Construction " << std::setw(20) + << 100.0 * static_cast(counter) / static_cast(n_nodes) << "%"; + mutex.unlock(); + }); + } + } + } + + m_cells.emplace_back(); + auto &cells = m_cells.back(); + cells.resize(m_n_cells); + const auto nx = n[0]; + const auto ny = n[1]; + for (auto l = 0u; l < m_n_cells; ++l) + { + // Build translation table for easy querying without having to recalculate the index positions + const auto j = l / nx; + const auto i = l % nx; + + auto &cell = cells[l]; + cell[0] = (nx + 1) * j + i; + cell[1] = cell[0] + 1; + cell[2] = (nx + 1) * (j + 1) + i; + cell[3] = cell[2] + 1; + + auto offset = nv; + cell[4] = offset + 2 * (nx * j + i); + cell[5] = cell[4] + 1; + cell[6] = offset + 2 * (nx * (j + 1) + i); + cell[7] = cell[6] + 1; + + offset += 2 * ne_x; + cell[8] = offset + 2 * (ny * i + j); + cell[9] = cell[8] + 1; + cell[10] = offset + 2 * (ny * (i + 1) + j); + cell[11] = cell[10] + 1; + } + + m_cell_map.emplace_back(); + auto &cell_map = m_cell_map.back(); + cell_map.resize(m_n_cells); + std::iota(cell_map.begin(), cell_map.end(), 0u); + + if (verbose) + { + std::cout << "\rConstruction took " << std::setw(15) << static_cast(duration_cast(high_resolution_clock::now() - t0_construction).count()) / 1000.0 << "s" << std::endl; + } + + return static_cast(m_n_fields++); +} + +bool CubicLagrangeDiscreteGrid2D::determineShapeFunctions(unsigned int field_id, Eigen::Vector2d const &x, std::array &cell, Eigen::Vector2d &c0, Eigen::Matrix &N, Eigen::Matrix *dN) const +{ + if (!m_domain.contains(x)) + return false; + + auto mi = (x - m_domain.min()).cwiseProduct(m_inv_cell_size).cast().eval(); + if (mi[0] >= m_resolution[0]) + mi[0] = m_resolution[0] - 1; + if (mi[1] >= m_resolution[1]) + mi[1] = m_resolution[1] - 1; + auto i = multiToSingleIndex({ { mi(0), mi(1) } }); + auto i_ = m_cell_map[field_id][i]; + if (i_ == std::numeric_limits::max()) + return false; + + auto sd = subdomain(i); + i = i_; + + // Mapping to isoparametric space from world (actually local?) space + const auto denom = (sd.max() - sd.min()).eval(); + c0 = Vector2d::Constant(2.0).cwiseQuotient(denom).eval(); + const auto c1 = (sd.max() + sd.min()).cwiseQuotient(denom).eval(); + const auto xi = (c0.cwiseProduct(x) - c1).eval(); + + cell = m_cells[field_id][i]; + N = shape_function_(xi, dN); + return true; +} + +double CubicLagrangeDiscreteGrid2D::interpolate(unsigned int field_id, Eigen::Vector2d const& xi, const std::array &cell, const Eigen::Vector2d &c0, const Eigen::Matrix &N, Eigen::Vector2d* gradient, Eigen::Matrix *dN) const +{ + if (!gradient) + { + auto phi = 0.0; + for (auto j = 0u; j < 12u; ++j) + { + auto v = cell[j]; + auto c = m_nodes[field_id][v]; + if (c == std::numeric_limits::max()) + { + return std::numeric_limits::max(); + } + phi += c * N[j]; + } + + return phi; + } + + auto phi = 0.0; + gradient->setZero(); + for (auto j = 0u; j < 12u; ++j) + { + auto v = cell[j]; + auto c = m_nodes[field_id][v]; + if (c == std::numeric_limits::max()) + { + gradient->setZero(); + return std::numeric_limits::max(); + } + phi += c * N[j]; + (*gradient)(0) += c * (*dN)(j, 0); + (*gradient)(1) += c * (*dN)(j, 1); + } + gradient->array() *= c0.array(); + + return phi; +} + +double CubicLagrangeDiscreteGrid2D::interpolate(unsigned int field_id, Vector2d const &x, Vector2d *gradient) const +{ + if (!m_domain.contains(x)) + return std::numeric_limits::max(); + + auto mi = (x - m_domain.min()).cwiseProduct(m_inv_cell_size).cast().eval(); + if (mi[0] >= m_resolution[0]) + mi[0] = m_resolution[0] - 1; + if (mi[1] >= m_resolution[1]) + mi[1] = m_resolution[1] - 1; + auto i = multiToSingleIndex({{mi(0), mi(1)}}); + auto i_ = m_cell_map[field_id][i]; + if (i_ == std::numeric_limits::max()) + return std::numeric_limits::max(); + + auto sd = subdomain(i); + i = i_; + auto d = sd.diagonal().eval(); + + auto denom = (sd.max() - sd.min()).eval(); + auto c0 = Vector2d::Constant(2.0).cwiseQuotient(denom).eval(); + auto c1 = (sd.max() + sd.min()).cwiseQuotient(denom).eval(); + auto xi = (c0.cwiseProduct(x) - c1).eval(); + + auto const &cell = m_cells[field_id][i]; + if (!gradient) + { + auto phi = 0.0; + auto N = shape_function_(xi, nullptr); + for (auto j = 0u; j < 12u; ++j) + { + auto v = cell[j]; + auto c = m_nodes[field_id][v]; + if (c == std::numeric_limits::max()) + { + return std::numeric_limits::max(); + } + phi += c * N[j]; + } + + return phi; + } + + auto dN = Matrix{}; + auto N = shape_function_(xi, &dN); + + auto phi = 0.0; + gradient->setZero(); + for (auto j = 0u; j < 12u; ++j) + { + auto v = cell[j]; + auto c = m_nodes[field_id][v]; + if (c == std::numeric_limits::max()) + { + gradient->setZero(); + return std::numeric_limits::max(); + } + phi += c * N[j]; + (*gradient)(0) += c * dN(j, 0); + (*gradient)(1) += c * dN(j, 1); + } + gradient->array() *= c0.array(); + + return phi; +} + +void CubicLagrangeDiscreteGrid2D::reduceField(unsigned int field_id, Predicate pred) +{ + auto &coeffs = m_nodes[field_id]; + auto &cells = m_cells[field_id]; + auto keep = std::vector(coeffs.size()); + for (auto l = 0u; l < coeffs.size(); ++l) + { + auto xi = indexToNodePosition(l); + keep[l] = pred(xi, coeffs[l]) && coeffs[l] != std::numeric_limits::max(); + } + + auto &cell_map = m_cell_map[field_id]; + cell_map.resize(m_n_cells); + std::iota(cell_map.begin(), cell_map.end(), 0u); + + const auto cells_ = cells; + cells.clear(); + for (auto i = 0u; i < cells_.size(); ++i) + { + auto keep_cell = false; + auto vals = std::vector{}; + for (auto v : cells_[i]) + { + keep_cell |= keep[v]; + vals.push_back(coeffs[v]); + } + if (keep_cell) + { + cells.push_back(cells_[i]); + cell_map[i] = static_cast(cells.size() - 1); + } + else + cell_map[i] = std::numeric_limits::max(); + } + + // Reduce vertices. + auto z_values = std::vector(coeffs.size()); + for (auto l = 0u; l < coeffs.size(); ++l) + { + auto xi = indexToNodePosition(l); + z_values[l] = zValue(xi, 4.0 * m_inv_cell_size.minCoeff()); + } + + std::fill(keep.begin(), keep.end(), false); + + auto vertex_to_cell = std::vector>>(coeffs.size()); + for (auto c = 0u; c < cells.size(); ++c) + { + auto const &cell = cells[c]; + + for (auto j = 0u; j < cell.size(); ++j) + { + auto v = cell[j]; + keep[v] = true; + vertex_to_cell[v].insert({c, j}); + } + } + auto last_vertex = static_cast(coeffs.size() - 1); + for (auto i = static_cast(coeffs.size() - 1); i >= 0; --i) + { + if (!keep[i]) + { + std::swap(coeffs[i], coeffs[last_vertex]); + std::swap(z_values[i], z_values[last_vertex]); + std::swap(vertex_to_cell[i], vertex_to_cell[last_vertex]); + for (auto const &kvp : vertex_to_cell[i]) + { + cells[kvp.first][kvp.second] = i; + } + for (auto const &kvp : vertex_to_cell[last_vertex]) + { + cells[kvp.first][kvp.second] = last_vertex; + } + + last_vertex--; + } + } + coeffs.resize(last_vertex + 1); + z_values.resize(coeffs.size()); + + auto sort_pattern = std::vector(coeffs.size()); + std::iota(sort_pattern.begin(), sort_pattern.end(), 0u); + std::sort(sort_pattern.begin(), sort_pattern.end(), + [&](unsigned int i, unsigned int j) { + return z_values[i] < z_values[j]; + }); + + for (auto i = 0u; i < sort_pattern.size(); ++i) + { + auto j = sort_pattern[i]; + for (auto const &kvp : vertex_to_cell[j]) + { + assert(cells[kvp.first][kvp.second] == j); + cells[kvp.first][kvp.second] = i; + } + } + + auto coeffs_ = coeffs; + std::transform(sort_pattern.begin(), sort_pattern.end(), coeffs.begin(), + [&coeffs_](unsigned int i) { return coeffs_[i]; }); +} + +void CubicLagrangeDiscreteGrid2D::forEachCell(unsigned int field_id, std::function const &cb) const +{ + auto n = m_resolution[0] * m_resolution[1]; + for (auto i = 0u; i < n; ++i) + { + auto domain = AlignedBox2d{}; + auto mi = singleToMultiIndex(i); + domain.min() = m_domain.min() + Matrix::Map(mi.data()).cast().cwiseProduct(m_cell_size); + domain.max() = domain.min() + m_cell_size; + + cb(i, domain, 0); + } +} + +} // namespace Discregrid diff --git a/discregrid/src/data/z_sort_table.hpp b/discregrid/src/data/z_sort_table.hpp index ea22ba4..e29921c 100755 --- a/discregrid/src/data/z_sort_table.hpp +++ b/discregrid/src/data/z_sort_table.hpp @@ -1,3 +1,28 @@ +// The pre-shifted tables are from https://github.com/Forceflow/libmorton/tree/main under the MIT license: +/* +MIT License + +Copyright (c) 2016 Jeroen Baert + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + #ifndef Z_SORT_TABLE_H__ #define Z_SORT_TABLE_H__ @@ -133,4 +158,92 @@ uint64_t morton_lut(std::array const& x) return answer; } +// 2D +static const unsigned int Morton2D_encode_x_256[256] = +{ + 0, 1, 4, 5, 16, 17, 20, 21, + 64, 65, 68, 69, 80, 81, 84, 85, + 256, 257, 260, 261, 272, 273, 276, 277, + 320, 321, 324, 325, 336, 337, 340, 341, + 1024, 1025, 1028, 1029, 1040, 1041, 1044, 1045, + 1088, 1089, 1092, 1093, 1104, 1105, 1108, 1109, + 1280, 1281, 1284, 1285, 1296, 1297, 1300, 1301, + 1344, 1345, 1348, 1349, 1360, 1361, 1364, 1365, + 4096, 4097, 4100, 4101, 4112, 4113, 4116, 4117, + 4160, 4161, 4164, 4165, 4176, 4177, 4180, 4181, + 4352, 4353, 4356, 4357, 4368, 4369, 4372, 4373, + 4416, 4417, 4420, 4421, 4432, 4433, 4436, 4437, + 5120, 5121, 5124, 5125, 5136, 5137, 5140, 5141, + 5184, 5185, 5188, 5189, 5200, 5201, 5204, 5205, + 5376, 5377, 5380, 5381, 5392, 5393, 5396, 5397, + 5440, 5441, 5444, 5445, 5456, 5457, 5460, 5461, + 16384, 16385, 16388, 16389, 16400, 16401, 16404, 16405, + 16448, 16449, 16452, 16453, 16464, 16465, 16468, 16469, + 16640, 16641, 16644, 16645, 16656, 16657, 16660, 16661, + 16704, 16705, 16708, 16709, 16720, 16721, 16724, 16725, + 17408, 17409, 17412, 17413, 17424, 17425, 17428, 17429, + 17472, 17473, 17476, 17477, 17488, 17489, 17492, 17493, + 17664, 17665, 17668, 17669, 17680, 17681, 17684, 17685, + 17728, 17729, 17732, 17733, 17744, 17745, 17748, 17749, + 20480, 20481, 20484, 20485, 20496, 20497, 20500, 20501, + 20544, 20545, 20548, 20549, 20560, 20561, 20564, 20565, + 20736, 20737, 20740, 20741, 20752, 20753, 20756, 20757, + 20800, 20801, 20804, 20805, 20816, 20817, 20820, 20821, + 21504, 21505, 21508, 21509, 21520, 21521, 21524, 21525, + 21568, 21569, 21572, 21573, 21584, 21585, 21588, 21589, + 21760, 21761, 21764, 21765, 21776, 21777, 21780, 21781, + 21824, 21825, 21828, 21829, 21840, 21841, 21844, 21845 +}; + +static const unsigned int Morton2D_encode_y_256[256] = +{ + 0, 2, 8, 10, 32, 34, 40, 42, + 128, 130, 136, 138, 160, 162, 168, 170, + 512, 514, 520, 522, 544, 546, 552, 554, + 640, 642, 648, 650, 672, 674, 680, 682, + 2048, 2050, 2056, 2058, 2080, 2082, 2088, 2090, + 2176, 2178, 2184, 2186, 2208, 2210, 2216, 2218, + 2560, 2562, 2568, 2570, 2592, 2594, 2600, 2602, + 2688, 2690, 2696, 2698, 2720, 2722, 2728, 2730, + 8192, 8194, 8200, 8202, 8224, 8226, 8232, 8234, + 8320, 8322, 8328, 8330, 8352, 8354, 8360, 8362, + 8704, 8706, 8712, 8714, 8736, 8738, 8744, 8746, + 8832, 8834, 8840, 8842, 8864, 8866, 8872, 8874, + 10240, 10242, 10248, 10250, 10272, 10274, 10280, 10282, + 10368, 10370, 10376, 10378, 10400, 10402, 10408, 10410, + 10752, 10754, 10760, 10762, 10784, 10786, 10792, 10794, + 10880, 10882, 10888, 10890, 10912, 10914, 10920, 10922, + 32768, 32770, 32776, 32778, 32800, 32802, 32808, 32810, + 32896, 32898, 32904, 32906, 32928, 32930, 32936, 32938, + 33280, 33282, 33288, 33290, 33312, 33314, 33320, 33322, + 33408, 33410, 33416, 33418, 33440, 33442, 33448, 33450, + 34816, 34818, 34824, 34826, 34848, 34850, 34856, 34858, + 34944, 34946, 34952, 34954, 34976, 34978, 34984, 34986, + 35328, 35330, 35336, 35338, 35360, 35362, 35368, 35370, + 35456, 35458, 35464, 35466, 35488, 35490, 35496, 35498, + 40960, 40962, 40968, 40970, 40992, 40994, 41000, 41002, + 41088, 41090, 41096, 41098, 41120, 41122, 41128, 41130, + 41472, 41474, 41480, 41482, 41504, 41506, 41512, 41514, + 41600, 41602, 41608, 41610, 41632, 41634, 41640, 41642, + 43008, 43010, 43016, 43018, 43040, 43042, 43048, 43050, + 43136, 43138, 43144, 43146, 43168, 43170, 43176, 43178, + 43520, 43522, 43528, 43530, 43552, 43554, 43560, 43562, + 43648, 43650, 43656, 43658, 43680, 43682, 43688, 43690 +}; + +inline +uint64_t morton_lut(std::array const& coords) +{ + uint64_t answer = 0; + answer = Morton2D_encode_y_256[(coords[1] >> 16) & 0xFF] | + Morton2D_encode_x_256[(coords[0] >> 16) & 0xFF]; + answer = (answer << 16) | + Morton2D_encode_y_256[(coords[1] >> 8) & 0xFF] | + Morton2D_encode_x_256[(coords[0] >> 8) & 0xFF]; + answer = (answer << 16) | + Morton2D_encode_y_256[(coords[1]) & 0xFF] | + Morton2D_encode_x_256[(coords[0]) & 0xFF]; + return answer; +} + #endif // Z_SORT_TABLE_H__ diff --git a/discregrid/src/discrete_grid2D.cpp b/discregrid/src/discrete_grid2D.cpp new file mode 100644 index 0000000..b46c909 --- /dev/null +++ b/discregrid/src/discrete_grid2D.cpp @@ -0,0 +1,39 @@ +#include + +using namespace Eigen; + +namespace Discregrid +{ + + +DiscreteGrid2D::MultiIndex +DiscreteGrid2D::singleToMultiIndex(unsigned int l) const +{ + const auto n0 = m_resolution[0]; + const auto j = l / n0; + const auto i = l % n0; + return {{i, j}}; +} + +unsigned int +DiscreteGrid2D::multiToSingleIndex(MultiIndex const & ij) const +{ + return ij[1] * m_resolution[0] + ij[0]; +} + +AlignedBox2d +DiscreteGrid2D::subdomain(MultiIndex const& ij) const +{ + auto origin = m_domain.min() + Map const>( + ij.data()).cast().cwiseProduct(m_cell_size); + return { origin, origin + m_cell_size}; +} + +AlignedBox2d +DiscreteGrid2D::subdomain(unsigned int l) const +{ + return subdomain(singleToMultiIndex(l)); +} + + +} \ No newline at end of file diff --git a/discregrid/src/mesh/triangle_mesh.cpp b/discregrid/src/mesh/triangle_mesh.cpp index 8a59453..f95f3a5 100755 --- a/discregrid/src/mesh/triangle_mesh.cpp +++ b/discregrid/src/mesh/triangle_mesh.cpp @@ -5,6 +5,7 @@ #include #include #include +#include using namespace Eigen; @@ -214,5 +215,73 @@ TriangleMesh::computeFaceNormal(unsigned int f) const return (x1 - x0).cross(x2 - x0).normalized(); } +struct Vector2dHash +{ + std::size_t operator()(const Vector2d& v) const + { + std::size_t seed = 0; + constexpr std::hash hasher; + seed ^= hasher(v.x()) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= hasher(v.y()) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + return seed; + } +}; + +struct Vector2dEqual +{ + bool operator()(const Vector2d& v1, const Vector2d& v2) const + { + return v1.isApprox(v2); + } +}; + +void TriangleMesh::weldVerticesAndCullBackfaces2D(std::vector& culled_vertices, std::vector& culled_triangles) const +{ + const Vector3d view_vector(0, 0, 1); // XY-Plane + + std::vector culled_and_welded_vertices; + std::vector culled_and_welded_triangles; + std::unordered_map vertex_to_index; + + auto weld_to_mesh = [&](const Vector2d& vertex) + { + const auto it = vertex_to_index.find(vertex); + if (it == vertex_to_index.end()) + { + const unsigned int index = culled_and_welded_vertices.size(); + culled_and_welded_vertices.push_back(vertex); + vertex_to_index[vertex] = index; + culled_and_welded_triangles.push_back(index); + } + else + { + culled_and_welded_triangles.push_back(it->second); + } + }; + + for (size_t i = 0; i < m_faces.size(); ++i) + { + const auto& triangle = m_faces[i]; + const auto i0 = triangle[0]; + const auto i1 = triangle[1]; + const auto i2 = triangle[2]; + Vector3d v0 = m_vertices[i0]; + Vector3d v1 = m_vertices[i1]; + Vector3d v2 = m_vertices[i2]; + Vector3d e1 = v1 - v0; + Vector3d e2 = v2 - v0; + + if (e1.cross(e2).dot(view_vector) < 0.0) + { + weld_to_mesh(v0.head<2>()); + weld_to_mesh(v1.head<2>()); + weld_to_mesh(v2.head<2>()); + } + } + + culled_vertices = culled_and_welded_vertices; + culled_triangles = culled_and_welded_triangles; +} + } diff --git a/cmd/discrete_field_to_bitmap/bmp_file.cpp b/discregrid/src/utility/bmp_file.cpp old mode 100755 new mode 100644 similarity index 95% rename from cmd/discrete_field_to_bitmap/bmp_file.cpp rename to discregrid/src/utility/bmp_file.cpp index 0de1675..2d27191 --- a/cmd/discrete_field_to_bitmap/bmp_file.cpp +++ b/discregrid/src/utility/bmp_file.cpp @@ -1,147 +1,148 @@ -#include -#include -#include "bmp_file.hpp" - -// Compression Type -#define BI_RGB 0L -#define BI_RLE8 1L -#define BI_RLE4 2L - - -// ------------------------------------------------------------------- -bool BmpReaderWriter::loadFile(const char *filename, unsigned int &width, unsigned int &height, unsigned char *&data) -{ - if (data) - { - delete [] data; - data = NULL; - } - width = 0; - height = 0; - - FILE *f = fopen(filename, "rb"); - if (!f) - return false; - - size_t num; - BMPHEADER header; - num = fread(&header, sizeof(BMPHEADER), 1, f); - if(isBigEndian()) header.Type = endianSwap(header.Type); - if (num != 1) { fclose(f); return false; } - if (header.Type != 'MB') { fclose(f); return false; } - - BMPINFO info; - num = fread(&info, sizeof(BMPINFO), 1, f); - if (num != 1) { fclose(f); return false; } - if(isBigEndian()) info.Size = endianSwap(info.Size); - if(isBigEndian()) info.BitCount = endianSwap(info.BitCount); - if(isBigEndian()) info.Compression = endianSwap(info.Compression); - if(isBigEndian()) info.Width = endianSwap(info.Width); - if(isBigEndian()) info.Height = endianSwap(info.Height); - - if (info.Size != sizeof(BMPINFO)) { fclose(f); return false; } - if (info.BitCount != 24) { fclose(f); return false; } - if (info.Compression != BI_RGB) { fclose(f); return false; } - - width = info.Width; - height = info.Height; - data = new unsigned char[width * height * 3]; - - int lineLen = (((info.Width * (info.BitCount>>3)) + 3)>>2)<<2; - unsigned char *line = new unsigned char[lineLen]; - - for(int i = info.Height-1; i >= 0; i--) { - num = fread(line, lineLen, 1, f); - if (num != 1) { fclose(f); return false; } - unsigned char *src = line; - unsigned char *dest = data + i*info.Width*3; - for(unsigned int j = 0; j < info.Width; j++) { - unsigned char r,g,b; - b = *src++; g = *src++; r = *src++; - *dest++ = r; *dest++ = g; *dest++ = b; - } - } - - delete [] line; - fclose(f); - - return true; -} - -// ------------------------------------------------------------------- -bool BmpReaderWriter::saveFile(const char *filename, int width, int height, unsigned char *data) -{ - FILE *f = fopen(filename, "wb"); - if (!f) return false; - - // todo : works on pcs only, swap correctly if big endian - BMPHEADER header; - header.Type = 'MB'; - header.Size = sizeof(BMPINFO); - header.Reserved1 = 0; - header.Reserved2 = 0; - header.OffBits = sizeof(BMPHEADER) + sizeof(BMPINFO); - fwrite(&header, sizeof(BMPHEADER), 1, f); - - BMPINFO info; - info.Size = sizeof(BMPINFO); - info.Width = width; - info.Height = height; - info.Planes = 1; - info.BitCount = 24; - info.Compression = BI_RGB; - info.XPelsPerMeter = 4000; - info.YPelsPerMeter = 4000; - info.ClrUsed = 0; - info.ClrImportant = 0; - fwrite(&info, sizeof(info), 1, f); - - // padded to multiple of 4 - int lineLen = (((info.Width * (info.BitCount>>3)) + 3)>>2)<<2; - info.SizeImage = lineLen * height; - - unsigned char *line = new unsigned char[lineLen]; - - for(int i = 0; i < height; i++) - { - unsigned char *src = data + i*width*3; - unsigned char *dest = line; - for(int j = 0; j < width; j++) - { - unsigned char r,g,b; - r = *src++; g = *src++; b = *src++; - *dest++ = b; *dest++ = g; *dest++ = r; - } - for (int j = 3*width; j < lineLen; j++) - *dest++ = 0; - fwrite(line, lineLen, 1, f); - } - - delete [] line; - fclose(f); - - return true; -} - -bool BmpReaderWriter::isBigEndian() -{ - int i = 1; return *((char*)&i) == 0; -} - -unsigned short BmpReaderWriter::endianSwap(unsigned short nValue) -{ - return (((nValue >> 8)) | (nValue << 8)); -} - -unsigned int BmpReaderWriter::endianSwap(unsigned int i) -{ - unsigned char b1, b2, b3, b4; - - b1 = i & 255; - b2 = (i >> 8) & 255; - b3 = (i >> 16) & 255; - b4 = (i >> 24) & 255; - - return ((unsigned int)b1 << 24) + ((unsigned int)b2 << 16) + ((unsigned int)b3 << 8) + b4; -} - +#include "utility/bmp_file.hpp" + +#include +#include + +// Compression Type +#define BI_RGB 0L +#define BI_RLE8 1L +#define BI_RLE4 2L + + +// ------------------------------------------------------------------- +bool BmpReaderWriter::loadFile(const char *filename, unsigned int &width, unsigned int &height, unsigned char *&data) +{ + if (data) + { + delete [] data; + data = NULL; + } + width = 0; + height = 0; + + FILE *f = fopen(filename, "rb"); + if (!f) + return false; + + size_t num; + BMPHEADER header; + num = fread(&header, sizeof(BMPHEADER), 1, f); + if(isBigEndian()) header.Type = endianSwap(header.Type); + if (num != 1) { fclose(f); return false; } + if (header.Type != 'MB') { fclose(f); return false; } + + BMPINFO info; + num = fread(&info, sizeof(BMPINFO), 1, f); + if (num != 1) { fclose(f); return false; } + if(isBigEndian()) info.Size = endianSwap(info.Size); + if(isBigEndian()) info.BitCount = endianSwap(info.BitCount); + if(isBigEndian()) info.Compression = endianSwap(info.Compression); + if(isBigEndian()) info.Width = endianSwap(info.Width); + if(isBigEndian()) info.Height = endianSwap(info.Height); + + if (info.Size != sizeof(BMPINFO)) { fclose(f); return false; } + if (info.BitCount != 24) { fclose(f); return false; } + if (info.Compression != BI_RGB) { fclose(f); return false; } + + width = info.Width; + height = info.Height; + data = new unsigned char[width * height * 3]; + + int lineLen = (((info.Width * (info.BitCount>>3)) + 3)>>2)<<2; + unsigned char *line = new unsigned char[lineLen]; + + for(int i = info.Height-1; i >= 0; i--) { + num = fread(line, lineLen, 1, f); + if (num != 1) { fclose(f); return false; } + unsigned char *src = line; + unsigned char *dest = data + i*info.Width*3; + for(unsigned int j = 0; j < info.Width; j++) { + unsigned char r,g,b; + b = *src++; g = *src++; r = *src++; + *dest++ = r; *dest++ = g; *dest++ = b; + } + } + + delete [] line; + fclose(f); + + return true; +} + +// ------------------------------------------------------------------- +bool BmpReaderWriter::saveFile(const char *filename, int width, int height, unsigned char *data) +{ + FILE *f = fopen(filename, "wb"); + if (!f) return false; + + // todo : works on pcs only, swap correctly if big endian + BMPHEADER header; + header.Type = 'MB'; + header.Size = sizeof(BMPINFO); + header.Reserved1 = 0; + header.Reserved2 = 0; + header.OffBits = sizeof(BMPHEADER) + sizeof(BMPINFO); + fwrite(&header, sizeof(BMPHEADER), 1, f); + + BMPINFO info; + info.Size = sizeof(BMPINFO); + info.Width = width; + info.Height = height; + info.Planes = 1; + info.BitCount = 24; + info.Compression = BI_RGB; + info.XPelsPerMeter = 4000; + info.YPelsPerMeter = 4000; + info.ClrUsed = 0; + info.ClrImportant = 0; + fwrite(&info, sizeof(info), 1, f); + + // padded to multiple of 4 + int lineLen = (((info.Width * (info.BitCount>>3)) + 3)>>2)<<2; + info.SizeImage = lineLen * height; + + unsigned char *line = new unsigned char[lineLen]; + + for(int i = 0; i < height; i++) + { + unsigned char *src = data + i*width*3; + unsigned char *dest = line; + for(int j = 0; j < width; j++) + { + unsigned char r,g,b; + r = *src++; g = *src++; b = *src++; + *dest++ = b; *dest++ = g; *dest++ = r; + } + for (int j = 3*width; j < lineLen; j++) + *dest++ = 0; + fwrite(line, lineLen, 1, f); + } + + delete [] line; + fclose(f); + + return true; +} + +bool BmpReaderWriter::isBigEndian() +{ + int i = 1; return *((char*)&i) == 0; +} + +unsigned short BmpReaderWriter::endianSwap(unsigned short nValue) +{ + return (((nValue >> 8)) | (nValue << 8)); +} + +unsigned int BmpReaderWriter::endianSwap(unsigned int i) +{ + unsigned char b1, b2, b3, b4; + + b1 = i & 255; + b2 = (i >> 8) & 255; + b3 = (i >> 16) & 255; + b4 = (i >> 24) & 255; + + return ((unsigned int)b1 << 24) + ((unsigned int)b2 << 16) + ((unsigned int)b3 << 8) + b4; +} +