diff --git a/Documentation/doc/biblio/geom.bib b/Documentation/doc/biblio/geom.bib index b1e2985919c..bfb054de162 100644 --- a/Documentation/doc/biblio/geom.bib +++ b/Documentation/doc/biblio/geom.bib @@ -152087,6 +152087,18 @@ keywords = {polygonal surface mesh, Surface reconstruction, kinetic framework, s publisher={Elsevier} } +@unpublished{lazard:hal-04907149, + TITLE = {{Removing self-intersections in 3D meshes while preserving floating-point coordinates}}, + AUTHOR = {Lazard, Sylvain and Valque, Leo}, + URL = {https://inria.hal.science/hal-04907149}, + NOTE = {working paper or preprint}, + YEAR = {2025}, + MONTH = Jan, + KEYWORDS = {Snap rounding ; mesh intersection ; robustness}, + PDF = {https://inria.hal.science/hal-04907149v1/file/Snap-HAL.pdf}, + HAL_ID = {hal-04907149}, + HAL_VERSION = {v1}, +} @inproceedings{si2005meshing, title={Meshing piecewise linear complexes by constrained {Delaunay} tetrahedralizations}, diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 8669729c808..ed6e884697b 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -39,6 +39,7 @@ - New implementation of `CGAL::Polygon_mesh_processing::clip()` with a plane as clipper that is much faster and is now able to handle non-triangulated surface meshes. - New implementation of `CGAL::Polygon_mesh_processing::split()` with a plane as clipper that is much faster and is now able to handle non-triangulated surface meshes. - Added the function `CGAL::Polygon_mesh_processing::refine_with_plane()`, which enables users to refine a mesh with their intersection with a plane. +- Added the parameter `apply_iterative_snap_rounding` to the function `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()`. When set to `true`, the coordinates are rounded to fit in double and may perform additional subdivisions to ensure the output is free of self-intersections. ### [Point Set Processing](https://doc.cgal.org/6.1/Manual/packages.html#PkgPointSetProcessing3) - Added `poisson_eliminate()` to downsample a point cloud to a target size while providing Poisson disk property, i.e., a larger minimal distance between points. diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt index c24a90d7904..ab1a6b06092 100644 --- a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt @@ -37,7 +37,21 @@ else() endif() create_single_source_cgal_program("fast.cpp") +create_single_source_cgal_program("rotated_cubes_autorefinement.cpp") +create_single_source_cgal_program("coplanar_cubes_autorefinement.cpp") + +create_single_source_cgal_program("Performance/performance_snap_polygon_soup.cpp") +create_single_source_cgal_program("Robustness/robustness_snap_polygon_soup.cpp") +create_single_source_cgal_program("Quality/quality_snap_polygon_soup.cpp") create_single_source_cgal_program("polygon_mesh_slicer.cpp") target_link_libraries(polygon_mesh_slicer PRIVATE CGAL::Eigen3_support) +find_package(TBB QUIET) +include(CGAL_TBB_support) +if(TARGET CGAL::TBB_support) + target_link_libraries(rotated_cubes_autorefinement PRIVATE CGAL::TBB_support) + target_link_libraries(coplanar_cubes_autorefinement PRIVATE CGAL::TBB_support) +else() + message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.") +endif() diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Performance/performance_snap_polygon_soup.cpp b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Performance/performance_snap_polygon_soup.cpp new file mode 100644 index 00000000000..d64cb27da9b --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Performance/performance_snap_polygon_soup.cpp @@ -0,0 +1,50 @@ +#include +#include +#include +#include +#include + +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef typename Kernel::Point_3 Point_3; +namespace PMP = CGAL::Polygon_mesh_processing; + +enum EXIT_CODES { VALID_OUTPUT=0, + INVALID_INPUT=1, + ROUNDING_FAILED=2, + SIGSEGV=10, + SIGSABRT=11, + SIGFPE=12, + TIMEOUT=13 + }; + +int main(int argc, char** argv) +{ + if(argc<4){ + std::cout << "Invalid argument" << std::endl; + return 1; + } + + const std::string filename = std::string(argv[1]); + const int grid_size = std::stoi(std::string(argv[2])); + const bool erase_duplicate = std::stoi(argv[3])==1; + + std::vector points; + std::vector> triangles; + + if (!CGAL::IO::read_polygon_soup(filename, points, triangles)) + { + std::cerr << "Cannot read " << filename << "\n"; + return 1; + } + + PMP::repair_polygon_soup(points, triangles); + PMP::triangulate_polygons(points, triangles); + + PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(erase_duplicate).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15)); + + return 0; +} diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Performance/run_performance.sh b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Performance/run_performance.sh new file mode 100755 index 00000000000..aa6ea1a28ba --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Performance/run_performance.sh @@ -0,0 +1,28 @@ +#!/bin/bash + +set -e + +if [ "$#" -lt 4 ]; then + echo "Usage: $0 [component_params...]" + exit 1 +fi + +INPUT_FILE=$1 +TIMEOUT=$2 +GRID_SIZE=$3 +ERASE_ALL_DUPLICATE=$4 + +# Use /usr/bin/time for memory usage (maximum resident set size in KB) +TMP_LOG=$(mktemp) + +# Run the benchmarked command +/usr/bin/time -f "TIME:%e\nMEM:%M" timeout "$TIMEOUT"s performance_snap_polygon_soup "$INPUT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" 2> "$TMP_LOG" + +# Parse time and memory +SECONDS=$(grep "TIME" "$TMP_LOG" | cut -d':' -f2) +MEMORY=$(grep "MEM" "$TMP_LOG" | cut -d':' -f2) + +rm -f "$TMP_LOG" + +# Output JSON +echo "{\"seconds\": \"$SECONDS\", \"memory_peaks\": \"$MEMORY\"}" diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Quality/quality_snap_polygon_soup.cpp b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Quality/quality_snap_polygon_soup.cpp new file mode 100644 index 00000000000..8e4bb99c2ba --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Quality/quality_snap_polygon_soup.cpp @@ -0,0 +1,60 @@ +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +typedef typename Kernel::Point_3 Point_3; +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char** argv) +{ + if(argc<4){ + std::cout << "Invalid argument" << std::endl; + return 1; + } + const std::string filename = std::string(argv[1]); + const int grid_size = std::stoi(std::string(argv[2])); + const bool erase_duplicate = std::stoi(argv[3])==1; + + std::vector points; + std::vector> triangles; + + CGAL::Bbox_3 bb = CGAL::bbox_3(points.begin(), points.end()); + double diag_length=std::sqrt((bb.xmax()-bb.xmin())*(bb.xmax()-bb.xmin()) + (bb.ymax()-bb.ymin())*(bb.ymax()-bb.ymin()) + (bb.zmax()-bb.zmin())*(bb.zmax()-bb.zmin())); + if (!CGAL::IO::read_polygon_soup(filename, points, triangles)) + { + std::cerr << "Cannot read " << filename << "\n"; + return 1; + } + + std::vector input_points(points.begin(), points.end()); + + PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(erase_duplicate).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15)); + + + std::cout << "{" << + "\"Nb_output_points\": \"" << points.size() << "\",\n" << + "\"Nb_output_triangles\": \"" << triangles.size() << "\",\n" << + "\"Is_2_manifold\": \"" << (PMP::orient_polygon_soup(points, triangles)?"True":"False") << "\",\n"; + CGAL::Surface_mesh sm; + PMP::polygon_soup_to_polygon_mesh(points, triangles, sm); + + std::cout << std::setprecision(17) << + "\"Hausdorff_distance_output_to_input_(divide_by_bbox_diag)\": \"" << PMP::max_distance_to_triangle_mesh(input_points, sm) / diag_length << "\",\n" << + "\"Closed_output\": \"" << (CGAL::is_closed(sm)?"True":"False") << "\",\n" << + "\"Ouput_bound_a_volume\": \"" << (PMP::does_bound_a_volume(sm)?"True":"False") << "\"\n}" + << std::endl; + + return 0; +} \ No newline at end of file diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Quality/run_quality.sh b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Quality/run_quality.sh new file mode 100755 index 00000000000..c3a194483eb --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Quality/run_quality.sh @@ -0,0 +1,19 @@ +#!/bin/bash + +set -e + +if [ "$#" -lt 4 ]; then + echo "Usage: $0 [component_params...]" + exit 1 +fi + +INPUT_FILE=$1 +TIMEOUT=$2 +GRID_SIZE=$3 +ERASE_ALL_DUPLICATE=$4 + +TMP_LOG=$(mktemp) +timeout "$TIMEOUT"s quality_snap_polygon_soup "$INPUT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" > "$TMP_LOG" + +cat $TMP_LOG +rm -f "$TMP_LOG" \ No newline at end of file diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Robustness/robustness_snap_polygon_soup.cpp b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Robustness/robustness_snap_polygon_soup.cpp new file mode 100644 index 00000000000..da424eb6422 --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Robustness/robustness_snap_polygon_soup.cpp @@ -0,0 +1,53 @@ +#include +#include +#include +#include +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef typename Kernel::Point_3 Point_3; +namespace PMP = CGAL::Polygon_mesh_processing; + +enum EXIT_CODES { VALID_OUTPUT=0, + INVALID_INPUT=1, + ROUNDING_FAILED=2, + SELF_INTERSECTING_OUTPUT=3, + SIGSEGV=10, + SIGSABRT=11, + SIGFPE=12, + TIMEOUT=13 + }; + +int main(int argc, char** argv) +{ + if(argc<4){ + std::cout << "Invalid argument" << std::endl; + return 1; + } + + const std::string filename = std::string(argv[1]); + const int grid_size = std::stoi(std::string(argv[2])); + const bool erase_duplicate = std::stoi(argv[3])==1; + + std::vector points; + std::vector> triangles; + + if (!CGAL::IO::read_polygon_soup(filename, points, triangles) || points.size()==0 || triangles.size()==0) + { + return INVALID_INPUT; + } + + PMP::repair_polygon_soup(points, triangles); + PMP::triangulate_polygons(points, triangles); + + bool success=PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(erase_duplicate).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15)); + + if(!success) + return ROUNDING_FAILED; + if( PMP::does_triangle_soup_self_intersect(points, triangles) ) + return SELF_INTERSECTING_OUTPUT; + + return VALID_OUTPUT; +} diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Robustness/run_robustness.sh b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Robustness/run_robustness.sh new file mode 100755 index 00000000000..4e4e28dcdfd --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/Robustness/run_robustness.sh @@ -0,0 +1,38 @@ +#!/bin/bash + +if [ "$#" -lt 4 ]; then + echo "Usage: $0 [component_params...]" + exit 1 +fi + +timeout_bis() { + timeout 5 sleep 10 +} + +INPUT_FILE=$1 +TIMEOUT=$2 +GRID_SIZE=$3 +ERASE_ALL_DUPLICATE=$4 + +# Run with timeout, capture exit code +timeout "--foreground" "$TIMEOUT"s robustness_snap_polygon_soup "$INPUT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" +EXIT_CODE=$? + +# Interpret exit codes +declare -A TAGS +TAGS[0]="VALID_OUTPUT" +TAGS[1]="INPUT_IS_INVALID" +TAGS[2]="ROUNDING_FAILED" +TAGS[3]="SELF_INTERSECTING_OUTPUT" +TAGS[139]="SIGSEGV" +TAGS[11]="SIGSEGV" +TAGS[6]="SIGABRT" +TAGS[8]="SIGFPE" +TAGS[132]="SIGILL" +TAGS[124]="TIMEOUT" + +TAG_NAME=${TAGS[$EXIT_CODE]:-UNKNOWN} +TAG_DESC=$([[ "$EXIT_CODE" -eq 0 ]] && echo "OK" || echo "Error") + +# Output JSON +echo "{\"TAG_NAME\": \"$TAG_NAME\", \"TAG\": \"$TAG_DESC\"}" \ No newline at end of file diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/benchmarking_snap_polygon_soup.sh b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/benchmarking_snap_polygon_soup.sh new file mode 100755 index 00000000000..da3f363566a --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/benchmarking_snap_polygon_soup.sh @@ -0,0 +1,116 @@ +#!/bin/bash + +# Temp directory for individual result JSONs +TMP_RESULT_DIR=$(mktemp -d) + +# Job control +JOBS=0 +MAX_JOBS=$NUM_THREADS + +# Function to process a single file +process_file() { + INPUT_PATH="$1" + INPUT_ID=$(basename "$INPUT_PATH" | cut -d. -f1) + COMPONENT_NAME="$2" + PROJECT_DIR="$3" + TIMEOUT="$4" + OUTPUT_DIR="$5" + TMP_RESULT_FILE="$6" + GRID_SIZE="$7" + ERASE_ALL_DUPLICATE="$8" + { + echo " \"$INPUT_ID\": {" + echo " \"path\": \"$INPUT_PATH\"," + + PERF_OUTPUT=$(bash "$PROJECT_DIR/Performance/run_performance.sh" "$INPUT_PATH" "$TIMEOUT" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" 2>> "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Performance/$INPUT_ID.log") + echo " \"Performance\": $PERF_OUTPUT," + + QUALITY_OUTPUT=$(bash "$PROJECT_DIR/Quality/run_quality.sh" "$INPUT_PATH" "$TIMEOUT" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" 2>> "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Quality/$INPUT_ID.log") + echo " \"Quality\": $QUALITY_OUTPUT," + + ROBUST_OUTPUT=$(bash "$PROJECT_DIR/Robustness/run_robustness.sh" "$INPUT_PATH" "$TIMEOUT" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" 2>> "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Robustness/$INPUT_ID.log") + echo " \"Robustness\": $ROBUST_OUTPUT" + + echo " }" + } > "$TMP_RESULT_FILE" +} +export -f process_file + +# Usage function +usage() { + echo "Usage: $0 [component_params...]" + exit 1 +} + +# Check parameters +if [ "$#" -lt 5 ]; then + usage +fi + +# Arguments +PROJECT_DIR=$1 +INPUT_DIR=$2 +OUTPUT_DIR=$3 +TIMEOUT=$4 +NUM_THREADS=$5 +GRID_SIZE=$6 +ERASE_ALL_DUPLICATE=$7 + +# Get component name from the project directory name +COMPONENT_NAME=$(basename "$PROJECT_DIR") +DATE_TAG=$(date +"%Y-%m-%d") +TIMESTAMP=$(date +"%Y-%m-%d %H:%M:%S") +RESULT_JSON="$OUTPUT_DIR/${COMPONENT_NAME}_results_${DATE_TAG}.json" + +# Compile +# Do not forget to define CGAL_DIR +cmake "$PROJECT_DIR" "-DCMAKE_BUILD_TYPE=Release" "-DCMAKE_CXX_FLAGS=-O3" +make -j $NUM_THREADS + +# Prepare log directories +mkdir -p "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Performance" +mkdir -p "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Quality" +mkdir -p "$OUTPUT_DIR/Logs/$COMPONENT_NAME/Robustness" + +# Initialize JSON +echo "{" > "$RESULT_JSON" +echo " \"$COMPONENT_NAME\": {" >> "$RESULT_JSON" +echo " \"Thingi10K\": {" >> "$RESULT_JSON" + +#process_file "$INPUT_DIR/100036.stl" "$COMPONENT_NAME" "$PROJECT_DIR" "$TIMEOUT" "$OUTPUT_DIR" "$TMP_RESULT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" +# Loop input files and spawn parallel jobs +for INPUT_FILE in "$INPUT_DIR"/*; do + INPUT_ID=$(basename "$INPUT_FILE" | cut -d. -f1) + TMP_RESULT_FILE="$TMP_RESULT_DIR/$INPUT_ID.json" + + process_file "$INPUT_FILE" "$COMPONENT_NAME" "$PROJECT_DIR" "$TIMEOUT" "$OUTPUT_DIR" "$TMP_RESULT_FILE" "$GRID_SIZE" "$ERASE_ALL_DUPLICATE" + + ((JOBS+=1)) + if [ "$JOBS" -ge "$NUM_THREADS" ]; then + wait + JOBS=0 + fi +done + +wait + +# Merge all partial JSONs +echo "{" > "$RESULT_JSON" +echo " \"$COMPONENT_NAME\": {" >> "$RESULT_JSON" +echo " \"Thingi10K\": {" >> "$RESULT_JSON" + +FIRST_ENTRY=true +for FILE in "$TMP_RESULT_DIR"/*.json; do + if [ "$FIRST_ENTRY" = true ]; then + FIRST_ENTRY=false + else + echo "," >> "$RESULT_JSON" + fi + cat "$FILE" >> "$RESULT_JSON" +done + +echo "" >> "$RESULT_JSON" +echo " }," >> "$RESULT_JSON" +echo " \"finished_at\": \"$TIMESTAMP\"" >> "$RESULT_JSON" +echo " }" >> "$RESULT_JSON" +echo "}" >> "$RESULT_JSON" \ No newline at end of file diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/coplanar_cubes_autorefinement.cpp b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/coplanar_cubes_autorefinement.cpp new file mode 100644 index 00000000000..b97ed8f0365 --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/coplanar_cubes_autorefinement.cpp @@ -0,0 +1,135 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef K::Point_3 Point_3; + +typedef CGAL::Simple_cartesian Cartesian; +typedef Cartesian::Point_3 Double_Point_3; + +namespace PMP = CGAL::Polygon_mesh_processing; +namespace params = CGAL::parameters; + +struct Sphere_function { + double radius; + Sphere_function(double r) : radius(r) {} + Kernel::FT operator()(const Kernel::Point_3& p) const { + return p.x()*p.x() + p.y()*p.y() + p.z()*p.z() - radius*radius; + } +}; + +//Thanks Roberto! +K::Aff_transformation_3 +random_rotation(CGAL::Random &gen) +{ + double a=gen.get_double(0,2*CGAL_PI); + + double ca = cos(a); + double sa = sin(a); + + K::Aff_transformation_3 aff(1, 0, 0, + 0, ca,-sa, + 0, sa, ca); + std::cout << "Rotation by " << a << std::endl; + return aff; +} + +int main(int argc, char** argv) +{ + Mesh cube; + std::vector points; + std::vector> faces; + + CGAL::make_hexahedron( + Point_3(-1,-1,-1), + Point_3(1,-1,-1), + Point_3(1,1,-1), + Point_3(-1,1,-1), + Point_3(-1,1,1), + Point_3(-1,-1,1), + Point_3(1,-1,1), + Point_3(1,1,1), + cube, + CGAL::parameters::do_not_triangulate_faces(false) + ); + + std::cout << "Iterative intersection of rotative cubes with snapping" << std::endl; + + int i=0; + CGAL::Random random_gen = argc == 1 ? CGAL::get_default_random() : CGAL::Random(std::stoi(argv[1])); + + Mesh inter=cube; + while(true) + { + std::cout << "Iteration " << i << std::endl; + + CGAL::Real_timer t; + t.start(); + + std::cout << "Add a randomly rotated cube to the scene" << std::endl; + Mesh rotated_cube=cube; + PMP::transform(random_rotation(random_gen), rotated_cube); + + std::cout << "compute_intersection" << std::endl; + bool OK = PMP::corefine_and_compute_intersection(inter, rotated_cube, inter); + + if(!OK){ + std::cout << "No manifold, stop experiment" << std::endl; + exit(0); + } + + points.clear(); + faces.clear(); + PMP::polygon_mesh_to_polygon_soup(inter, points, faces); + + std::cout << "Snapped the points on double" << std::endl; + bool success=PMP::autorefine_triangle_soup(points, faces, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag())); + t.stop(); + if(!success){ + std::cout << "Rounding failed" << std::endl; + exit(0); + } + + //dump model every 100 iterations + if(i%100==0){ + std::cout << "Dump model" << std::endl; + std::vector double_points; + for(auto &p: points) + double_points.emplace_back(CGAL::to_double(p.x()),CGAL::to_double(p.y()),CGAL::to_double(p.z())); + std::ofstream outfile("cubes_"+std::to_string(i)+".off"); + outfile.precision(17); + outfile << "OFF\n" << points.size() << " " << faces.size() << " 0\n"; + for(auto p : points) + outfile << p.x() << " " << p.y() << " " << p.z() << std::endl; + + for(auto &t : faces) + outfile << "3" << " " << t[0] << " " << t[1] << " " << t[2] << std::endl; + outfile.close();// + } + + std::cout << "#points = " << points.size() << " and #triangles = " << faces.size() << " in " << t.time() << " sec.\n\n" << std::endl; + + inter.clear(); + PMP::polygon_soup_to_polygon_mesh(points, faces, inter); + ++i; + } + + return 0; +} + diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/rotated_cubes_autorefinement.cpp b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/rotated_cubes_autorefinement.cpp new file mode 100644 index 00000000000..69ae0b8c83e --- /dev/null +++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/rotated_cubes_autorefinement.cpp @@ -0,0 +1,137 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel K; +typedef CGAL::Surface_mesh Mesh; +typedef K::Point_3 Point_3; + +typedef CGAL::Simple_cartesian Cartesian; +typedef Cartesian::Point_3 Double_Point_3; + +namespace PMP = CGAL::Polygon_mesh_processing; +namespace params = CGAL::parameters; + +struct Sphere_function { + double radius; + Sphere_function(double r) : radius(r) {} + Kernel::FT operator()(const Kernel::Point_3& p) const { + return p.x()*p.x() + p.y()*p.y() + p.z()*p.z() - radius*radius; + } +}; + +//Thanks Roberto! +K::Aff_transformation_3 +random_rotation(CGAL::Random &gen) +{ + double a=gen.get_double(0,2*CGAL_PI); + double b=gen.get_double(0,2*CGAL_PI); + double c=gen.get_double(0,2*CGAL_PI); + + double ca = cos(a), cb = cos(b), cc = cos(c); + double sa = sin(a), sb = sin(b), sc = sin(c); + + K::Aff_transformation_3 aff(cb * cc, cc* sa* sb - ca * sc, ca* cc* sb + sa * sc, + cb* sc, ca* cc + sa * sb * sc, ca* sb* sc - cc * sa, + -sb, cb* sa, ca* cb); + std::cout << "Rotation by " << a << " " << b << " " << c << std::endl; + return aff; +} + +int main(int argc, char** argv) +{ + Mesh cube; + std::vector points; + std::vector> faces; + + CGAL::make_hexahedron( + Point_3(-1,-1,-1), + Point_3(1,-1,-1), + Point_3(1,1,-1), + Point_3(-1,1,-1), + Point_3(-1,1,1), + Point_3(-1,-1,1), + Point_3(1,-1,1), + Point_3(1,1,1), + cube, + CGAL::parameters::do_not_triangulate_faces(false) + ); + + std::cout << "Iterative intersection of rotative cubes with snapping" << std::endl; + + int i=0; + CGAL::Random random_gen = argc == 1 ? CGAL::get_default_random() : CGAL::Random(std::stoi(argv[1])); + + Mesh inter=cube; + while(true) + { + std::cout << "Iteration " << i << std::endl; + + CGAL::Real_timer t; + t.start(); + + std::cout << "Add a randomly rotated cube to the scene" << std::endl; + Mesh rotated_cube=cube; + PMP::transform(random_rotation(random_gen), rotated_cube); + + std::cout << "compute_intersection" << std::endl; + bool OK = PMP::corefine_and_compute_intersection(inter, rotated_cube, inter); + + if(!OK){ + std::cout << "No manifold, stop experiment" << std::endl; + exit(0); + } + + points.clear(); + faces.clear(); + PMP::polygon_mesh_to_polygon_soup(inter, points, faces); + + std::cout << "Snapped the points on double" << std::endl; + bool success=PMP::autorefine_triangle_soup(points, faces, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(true).concurrency_tag(CGAL::Parallel_if_available_tag())); + t.stop(); + if(!success){ + std::cout << "Rounding failed" << std::endl; + exit(0); + } + + //dump model every 100 iterations + if(i%100==0){ + std::cout << "Dump model" << std::endl; + std::vector double_points; + for(auto &p: points) + double_points.emplace_back(CGAL::to_double(p.x()),CGAL::to_double(p.y()),CGAL::to_double(p.z())); + std::ofstream outfile("cubes_"+std::to_string(i)+".off"); + outfile.precision(17); + outfile << "OFF\n" << points.size() << " " << faces.size() << " 0\n"; + for(auto p : points) + outfile << p.x() << " " << p.y() << " " << p.z() << std::endl; + + for(auto &t : faces) + outfile << "3" << " " << t[0] << " " << t[1] << " " << t[2] << std::endl; + outfile.close();// + } + + std::cout << "#points = " << points.size() << " and #triangles = " << faces.size() << " in " << t.time() << " sec.\n\n" << std::endl; + + inter.clear(); + PMP::polygon_soup_to_polygon_mesh(points, faces, inter); + ++i; + } + + return 0; +} + diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h index ef2f124329b..7cff8ea4bc0 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h @@ -23,5 +23,7 @@ public: /// called for each subtriangle created from a triangle with intersection, `tgt_id` is the position in the triangle container after calling /// `autorefine_triangle_soup()` of the subtriangle, while `src_id` was the position of the original support triangle before calling the function. void new_subtriangle(std::size_t tgt_id, std::size_t src_id); + /// called for each input triangle absent in the output because it was degenerate. `src_id` is the position of that triangle in the input range. Additionally, if `apply_iterative_snap_rounding()` is set to `true`, some extra triangles might become degenerate and will be absent in the output. Those triangles are also reported by this function. + void delete_triangle(std::size_t src_id); /// @} }; diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index ef3f3c544b8..883a40175d1 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt @@ -4,7 +4,7 @@ namespace CGAL { \anchor Chapter_PolygonMeshProcessing \cgalAutoToc -\authors David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katrioplas, Sébastien Loriot, Ivan Pađen, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, Sébastien Valette, and Ilker %O. Yaz +\authors David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katrioplas, Sébastien Loriot, Ivan Pađen, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, Sébastien Valette, Léo Valque, and Ilker %O. Yaz \image html neptun_head.jpg \image latex neptun_head.jpg @@ -911,8 +911,10 @@ would then also includes overlaps of duplicated points. The function `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` provides a way to refine a triangle soup using the intersections of the triangles from the soup. In particular, if some points are duplicated they will be merged. Note that if a kernel with exact predicates but inexact constructions is used, some new self-intersections -might be introduced due to rounding issues of points coordinates. -To guarantee that the triangle soup is free from self-intersections, a kernel with exact constructions must be used. +might be introduced due to the rounding of the coordinates of intersection points. The `apply_iterative_snap_rounding` option can be used to resolve this issue. +When set to `true`, it ensures that the coordinates are rounded to fit in `double` with potential additional subdivisions, +preventing any self-intersections from occurring. + \subsection PMPRemoveCapsNeedles Removal of Almost Degenerate Triangle Faces Triangle faces of a mesh made up of almost collinear points are badly shaped elements that @@ -1477,6 +1479,7 @@ supervision of Sébastien Valette and Sébastien Loriot, who later finalized the The implementation is based on \cgalCite{cgal:vcp-grtmmdvd-08}. and preceding work. ACVD's implementation was also used as a reference during the project. +The `apply_iterative_snap_rounding` option for autorefinement was implemented in 2025, by Léo Valque, based on his work with Sylvain Lazard \cgalCite{lazard:hal-04907149}. */ } /* namespace CGAL */ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt index f2005fa8b4f..4d657bcef5e 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt @@ -53,4 +53,5 @@ \example Polygon_mesh_processing/acvd_remeshing_example.cpp \example Polygon_mesh_processing/sample_example.cpp \example Polygon_mesh_processing/soup_autorefinement.cpp +\example Polygon_mesh_processing/snap_polygon_soup.cpp */ diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index f5035cd158c..717ff715394 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -58,6 +58,7 @@ create_single_source_cgal_program("isotropic_remeshing_with_custom_sizing_exampl create_single_source_cgal_program("isotropic_remeshing_with_allow_move.cpp") create_single_source_cgal_program("triangle_mesh_autorefinement.cpp") create_single_source_cgal_program("soup_autorefinement.cpp") +create_single_source_cgal_program("snap_polygon_soup.cpp") find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater) include(CGAL_Eigen3_support) @@ -140,6 +141,7 @@ if(TARGET CGAL::TBB_support) target_link_libraries(hausdorff_distance_remeshing_example PRIVATE CGAL::TBB_support) target_link_libraries(hausdorff_bounded_error_distance_example PRIVATE CGAL::TBB_support) target_link_libraries(soup_autorefinement PRIVATE CGAL::TBB_support) + target_link_libraries(snap_polygon_soup PRIVATE CGAL::TBB_support) create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp") target_link_libraries(corefinement_parallel_union_meshes PRIVATE CGAL::TBB_support) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/snap_polygon_soup.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/snap_polygon_soup.cpp new file mode 100644 index 00000000000..61a69b83ee7 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/snap_polygon_soup.cpp @@ -0,0 +1,62 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char** argv) +{ + const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off") + : std::string(argv[1]); + + const int grid_size = argc <= 2 ? 23 + : std::stoi(std::string(argv[2])); + + const std::string out_file = "rounded_soup.off"; + + std::vector points; + std::vector> triangles; + + std::cout << "Snap rounding on " << filename << "\n"; + if (!CGAL::IO::read_polygon_soup(filename, points, triangles)) + { + std::cerr << "Cannot read " << filename << "\n"; + return 1; + } + + PMP::repair_polygon_soup(points, triangles); + PMP::triangulate_polygons(points, triangles); + + std::cout << "#points = " << points.size() << " and #triangles = " << triangles.size() << std::endl; + std::cout << "Is 2-manifold: " << PMP::is_polygon_soup_a_polygon_mesh(triangles) << std::endl; + + std::vector> pairs_of_intersecting_triangles; + PMP::triangle_soup_self_intersections(points, triangles, std::back_inserter(pairs_of_intersecting_triangles)); + std::cout << "Nb of pairs of intersecting triangles: " << pairs_of_intersecting_triangles.size() << std::endl; + + CGAL::Real_timer t; + t.start(); + bool success=PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::apply_iterative_snap_rounding(true).erase_all_duplicates(false).concurrency_tag(CGAL::Parallel_if_available_tag()).snap_grid_size(grid_size).number_of_iterations(15)); + t.stop(); + + std::cout << "\nOutput:" << std::endl; + std::cout << "#points = " << points.size() << " and #triangles = " << triangles.size() << " in " << t.time() << " sec." << std::endl; + if(success) + std::cout << "Does self-intersect: " << PMP::does_triangle_soup_self_intersect(points, triangles) << std::endl; + else + std::cout << "ROUNDING FAILED" << std::endl; + + CGAL::IO::write_polygon_soup(out_file, points, triangles, CGAL::parameters::stream_precision(17)); + std::cout << "Is 2-manifold: " << PMP::orient_polygon_soup(points, triangles) << "\n\n" << std::endl; + return 0; +} diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index a2dafa3dd36..00a48f919e7 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -33,7 +33,6 @@ #include #include - #ifdef CGAL_PMP_AUTOREFINE_USE_DEFAULT_VERBOSE #define CGAL_PMP_AUTOREFINE_VERBOSE(X) std::cout << X << "\n"; #endif @@ -90,6 +89,7 @@ struct Default_visitor inline void number_of_output_triangles(std::size_t /*nbt*/) {} inline void verbatim_triangle_copy(std::size_t /*tgt_id*/, std::size_t /*src_id*/) {} inline void new_subtriangle(std::size_t /*tgt_id*/, std::size_t /*src_id*/) {} + inline void delete_triangle(std::size_t /*src_id*/) {} }; } // end of Autorefinement visitor @@ -975,59 +975,19 @@ void generate_subtriangles(std::size_t ti, } // end of autorefine_impl #endif -/** -* \ingroup PMP_corefinement_grp -* -* refines a soup of triangles so that no pair of triangles intersects. -* Output triangles may share a common edge or a common vertex (but with the same indexed position in `points`). -* Note that points in `soup_points` can only be added (intersection points) at the end of the container, with the initial order preserved. -* Note that if `soup_points` contains two or more identical points then only the first copy (following the order in the `soup_points`) -* will be used in `soup_triangles`. -* `soup_triangles` will be updated to contain both the input triangles and the new subdivided triangles. Degenerate triangles will be removed. -* Also triangles in `soup_triangles` will be triangles without intersection first, followed by triangles coming from a subdivision induced -* by an intersection. The named parameter `visitor()` can be used to track -* -* @tparam PointRange a model of the concept `RandomAccessContainer` -* whose value type is the point type -* @tparam TriangleRange a model of the concepts `RandomAccessContainer`, `BackInsertionSequence` and `Swappable`, whose -* value type is a model of the concept `RandomAccessContainer` whose value type is convertible to `std::size_t` and that -* is constructible from an `std::initializer_list` of size 3. -* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" -* -* @param soup_points points of the soup of polygons -* @param soup_triangles each element in the range describes a triangle using the indexed position of the points in `soup_points` -* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below -* -* \cgalNamedParamsBegin -* \cgalParamNBegin{concurrency_tag} -* \cgalParamDescription{a tag indicating if the task should be done using one or several threads.} -* \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`} -* \cgalParamDefault{`CGAL::Sequential_tag`} -* \cgalParamNEnd -* \cgalParamNBegin{point_map} -* \cgalParamDescription{a property map associating points to the elements of the range `soup_points`} -* \cgalParamType{a model of `ReadWritePropertyMap` whose value type is a point type} -* \cgalParamDefault{`CGAL::Identity_property_map`} -* \cgalParamNEnd -* \cgalParamNBegin{geom_traits} -* \cgalParamDescription{an instance of a geometric traits class} -* \cgalParamType{a class model of `Kernel`} -* \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} -* \cgalParamExtra{The geometric traits class must be compatible with the point type.} -* \cgalParamNEnd -* \cgalParamNBegin{visitor} -* \cgalParamDescription{a visitor used to track the creation of new faces} -* \cgalParamType{a class model of `PMPAutorefinementVisitor`} -* \cgalParamDefault{`Autorefinement::Default_visitor`} -* \cgalParamExtra{The visitor will be copied.} -* \cgalParamNEnd -* \cgalNamedParamsEnd -* -*/ +namespace autorefine_impl{ +// Forward declaration +struct Wrap_snap_visitor; + +template +bool polygon_soup_snap_rounding(PointRange &points, + PolygonRange &triangles, + const NamedParameters& np = parameters::default_values()); + template -void autorefine_triangle_soup(PointRange& soup_points, - TriangleRange& soup_triangles, - const NamedParameters& np = parameters::default_values()) +bool autorefine_triangle_soup(PointRange& soup_points, + TriangleRange& soup_triangles, + const NamedParameters& np = parameters::default_values()) { using parameters::choose_parameter; using parameters::get_parameter; @@ -1050,7 +1010,6 @@ void autorefine_triangle_soup(PointRange& soup_points, > ::type Visitor; Visitor visitor(choose_parameter(get_parameter(np, internal_np::visitor))); - constexpr bool parallel_execution = std::is_same_v; #ifndef CGAL_LINKED_WITH_TBB @@ -1077,7 +1036,7 @@ void autorefine_triangle_soup(PointRange& soup_points, for(std::size_t i=0; i tri_inter_ids_inverse(triangles.size()); for (Input_TID f=0; f& r) { for (size_t ti = r.begin(); ti != r.end(); ++ti) { const std::array& t = new_triangles[ti].first; - visitor.new_subtriangle(offset+ti, tri_inter_ids_inverse[new_triangles[ti].second]); triangle_buffer[ti] = CGAL::make_array(concurrent_get_point_id(t[0]), concurrent_get_point_id(t[1]), concurrent_get_point_id(t[2])); } } ); - tbb::parallel_for(tbb::blocked_range(0, new_triangles.size()), - [&](const tbb::blocked_range& r) { - for (size_t ti = r.begin(); ti != r.end(); ++ti) - { - soup_triangles_out[offset + ti] = - { triangle_buffer[ti][0]->second, - triangle_buffer[ti][1]->second, - triangle_buffer[ti][2]->second }; + + // The constexpr was initially inside the lambda, but that did not compile with VC 2017 + if constexpr(std::is_same_v){ + tbb::parallel_for(tbb::blocked_range(0, new_triangles.size()), + [&](const tbb::blocked_range& r) { + for (size_t ti = r.begin(); ti != r.end(); ++ti) + { + soup_triangles_out[offset + ti] = + { triangle_buffer[ti][0]->second, + triangle_buffer[ti][1]->second, + triangle_buffer[ti][2]->second }; + visitor.new_subdivision(soup_triangles_out[offset + ti], soup_triangles[tri_inter_ids_inverse[new_triangles[ti].second]]); + } } - } - ); + ); + } + else + { + tbb::parallel_for(tbb::blocked_range(0, new_triangles.size()), + [&](const tbb::blocked_range& r) { + for (size_t ti = r.begin(); ti != r.end(); ++ti) + { + soup_triangles_out[offset + ti] = + { triangle_buffer[ti][0]->second, + triangle_buffer[ti][1]->second, + triangle_buffer[ti][2]->second }; + visitor.new_subtriangle(offset+ti, tri_inter_ids_inverse[new_triangles[ti].second]); + } + } + ); + } #else //option 2 (without mutex) /// Lambda concurrent_get_point_id() @@ -1591,17 +1571,33 @@ void autorefine_triangle_soup(PointRange& soup_points, } ); - tbb::parallel_for(tbb::blocked_range(0, new_triangles.size()), - [&](const tbb::blocked_range& r) { - for (size_t ti = r.begin(); ti != r.end(); ++ti) - { - soup_triangles_out[offset + ti] = - { triangle_buffer[ti][0]->second, - triangle_buffer[ti][1]->second, - triangle_buffer[ti][2]->second }; + // The constexpr was initially inside the lambda, but that did not compile with VC 2017 + if constexpr(std::is_same_v){ + tbb::parallel_for(tbb::blocked_range(0, new_triangles.size()), + [&](const tbb::blocked_range& r) { + for (size_t ti = r.begin(); ti != r.end(); ++ti) + { + soup_triangles_out[offset + ti] = + { triangle_buffer[ti][0]->second, + triangle_buffer[ti][1]->second, + triangle_buffer[ti][2]->second }; + visitor.new_subdivision(soup_triangles_out[offset + ti], soup_triangles[tri_inter_ids_inverse[new_triangles[ti].second]]); + } } - } - ); + ); + }else{ + tbb::parallel_for(tbb::blocked_range(0, new_triangles.size()), + [&](const tbb::blocked_range& r) { + for (size_t ti = r.begin(); ti != r.end(); ++ti) + { + soup_triangles_out[offset + ti] = + { triangle_buffer[ti][0]->second, + triangle_buffer[ti][1]->second, + triangle_buffer[ti][2]->second }; + } + } + ); + } #endif } else @@ -1613,10 +1609,13 @@ void autorefine_triangle_soup(PointRange& soup_points, soup_triangles_out.reserve(offset + new_triangles.size()); for (const std::pair, std::size_t>& t_and_id : new_triangles) { - visitor.new_subtriangle(soup_triangles_out.size(), tri_inter_ids_inverse[t_and_id.second]); soup_triangles_out.push_back({ get_point_id(t_and_id.first[0]), get_point_id(t_and_id.first[1]), get_point_id(t_and_id.first[2]) }); + if constexpr(std::is_same_v) + visitor.new_subdivision(soup_triangles_out[soup_triangles_out.size()-1], soup_triangles[tri_inter_ids_inverse[t_and_id.second]]); + else + visitor.new_subtriangle(soup_triangles_out.size()-1, tri_inter_ids_inverse[t_and_id.second]); } } @@ -1642,8 +1641,108 @@ void autorefine_triangle_soup(PointRange& soup_points, swap(soup_triangles, soup_triangles_out); CGAL_PMP_AUTOREFINE_VERBOSE("done"); + return true; } +} // end of autorefine_impl + +/** +* \ingroup PMP_corefinement_grp +* +* refines a soup of triangles so that no pair of triangles intersects. +* Output triangles may share a common edge or a common vertex (but with the same indexed position in `points`). +* Note that if `apply_iterative_snap_rounding` option is set to `false`, points in `soup_points` can only be added (intersection points) at the end of the container, with the initial order preserved. +* Note that if `soup_points` contains two or more identical points then only the first copy (following the order in the `soup_points`) +* will be used in `soup_triangles`. if `apply_iterative_snap_rounding` is set to `true`, all duplicates points are removed. +* `soup_triangles` will be updated to contain both the input triangles and the new subdivided triangles. Degenerate triangles will be removed. +* Also if `apply_iterative_snap_rounding` option is set to `false`, triangles in `soup_triangles` will be triangles without intersection first, followed by triangles coming from a subdivision induced +* by an intersection. The named parameter `visitor()` can be used to track the creation and removal of triangles independantly of +* the `apply_iterative_snap_rounding` option. +* If the `apply_iterative_snap_rounding` parameter is set to `true`, the coordinates of the vertices are rounded to fit within the precision of a double-precision floating point, +* while trying to make the triangle soup free of intersections. The `snap_grid_size()` parameter limits the drift of the snapped vertices. +* A smaller value is more likely to output an intersection free output and perform more vertex collapses, but it may increase the Hausdorff distance from the input. +* +* @tparam PointRange a model of the concept `RandomAccessContainer` +* whose value type is the point type +* @tparam TriangleRange a model of the concepts `RandomAccessContainer`, `BackInsertionSequence` and `Swappable`, whose +* value type is a model of the concept `RandomAccessContainer` whose value type is convertible to `std::size_t` and that +* is constructible from an `std::initializer_list` of size 3. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" +* +* @param soup_points points of the soup of polygons +* @param soup_triangles each element in the range describes a triangle using the indexed position of the points in `soup_points` +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* +* \cgalNamedParamsBegin +* \cgalParamNBegin{concurrency_tag} +* \cgalParamDescription{a tag indicating if the task should be done using one or several threads.} +* \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`} +* \cgalParamDefault{`CGAL::Sequential_tag`} +* \cgalParamNEnd +* \cgalParamNBegin{point_map} +* \cgalParamDescription{a property map associating points to the elements of the range `soup_points`} +* \cgalParamType{a model of `ReadWritePropertyMap` whose value type is a point type} +* \cgalParamDefault{`CGAL::Identity_property_map`} +* \cgalParamNEnd +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class} +* \cgalParamType{a class model of `Kernel`} +* \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} +* \cgalParamExtra{The geometric traits class must be compatible with the point type.} +* \cgalParamNEnd +* \cgalParamNBegin{visitor} +* \cgalParamDescription{a visitor used to track the creation of new faces} +* \cgalParamType{a class model of `PMPAutorefinementVisitor`} +* \cgalParamDefault{`Autorefinement::Default_visitor`} +* \cgalParamExtra{The visitor will be copied.} +* \cgalParamNEnd +* \cgalParamNBegin{apply_iterative_snap_rounding} +* \cgalParamDescription{Enable the rounding of the coordinates so that they fit in doubles.} +* \cgalParamType{Boolean} +* \cgalParamDefault{false} +* \cgalParamNEnd +* \cgalParamNBegin{snap_grid_size} +* \cgalParamDescription{A value `gs` used to scale the points to `[-2^gs, 2^gs]` before rounding them on integers. Used only if `apply_iterative_snap_rounding()` is set to `true`} +* \cgalParamType{`unsigned int`} +* \cgalParamDefault{23} +* \cgalParamExtra{Must be lower than 52.} +* \cgalParamNEnd +* \cgalParamNBegin{number_of_iterations} +* \cgalParamDescription{Maximum number of iterations performed by the snap rounding algorithm. Used only if `apply_iterative_snap_rounding` is true.} +* \cgalParamType{unsigned int} +* \cgalParamDefault{5} +* \cgalParamNEnd +* \cgalNamedParamsEnd +* +* \return `true` if `apply_iterative_snap_rounding` is set to `false`. +* Otherwise, returns `true` if the modified triangle soup is free of self-intersections, +* or `false` if the algorithm was unable to produce such a result within the allowed number of iterations. +* In the latter case, the output triangle soup represents a partial result from the final iteration, +* with no guarantee of its validity. +* +*/ +template +bool autorefine_triangle_soup(PointRange& soup_points, + TriangleRange& soup_triangles, + const NamedParameters& np = parameters::default_values()) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + //Dispatch the execution according the apply_iterative_snap_rounding parameter + const bool do_snap = choose_parameter(get_parameter(np, internal_np::apply_iterative_snap_rounding), false); + if(do_snap) + { + CGAL_PMP_AUTOREFINE_VERBOSE("Snap polygon soup"); + return autorefine_impl::polygon_soup_snap_rounding(soup_points, soup_triangles, np); + } + else + { + return autorefine_impl::autorefine_triangle_soup(soup_points, soup_triangles, np); + } +} + + /** * \ingroup PMP_corefinement_grp * refines a triangle mesh so that no triangles intersects in their interior. @@ -1717,4 +1816,6 @@ autorefine( TriangleMesh& tm, #endif #endif +#include + #endif // CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h new file mode 100644 index 00000000000..2d5272364b5 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h @@ -0,0 +1,585 @@ +// Copyright (c) 2025 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Léo Valque + +#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_POLYGON_SOUP_SNAP_ROUNDING_H +#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_POLYGON_SOUP_SNAP_ROUNDING_H + +#include + +#include +#include + +#include +#include + +#include +#include + +namespace CGAL +{ + +namespace Polygon_mesh_processing +{ + +namespace autorefine_impl +{ + +// Certified ceil function for exact number types +template double double_ceil(const Lazy_exact_nt< NT > &x); +template double double_ceil(const NT &x); + +template +double double_ceil(const Lazy_exact_nt< NT > &x){ + // If both sides are in the same ceil, return this ceil + double ceil_left=std::ceil(to_interval(x).first); + if(ceil_left==std::ceil(to_interval(x).second)) + return ceil_left; + // If not refine the interval by contracting the DAG and try again + x.exact(); + ceil_left=std::ceil(to_interval(x).first); + if(ceil_left==std::ceil(to_interval(x).second)) + return ceil_left; + // If not return the ceil of the exact value + return double_ceil( x.exact()); +}; + +template +double double_ceil(const NT &x){ + using FT = Fraction_traits; + if constexpr(FT::Is_fraction::value){ + // If NT is a fraction, the ceil value is the result of the euclidian division of the numerator and the denominator. + typename FT::Numerator_type num, r, e; + typename FT::Denominator_type denom; + typename FT::Decompose()(x,num,denom); + div_mod(num, denom, r, e); + if((r>=0) && e!=0) //If the result is positive, the ceil value is one above + return to_double(r+1); + return to_double(r); + } else { + // Return the ceil of the approximation + return std::ceil(to_double(x)); + } +}; + +// Provides an index to the range, which is used by the visitor to track the correspondence between the input and the output +template +class Indexes_range : public Range{ + typedef std::remove_cv_t::value_type> Value_type; +public: + typedef typename Range::const_iterator const_iterator; + typedef typename Range::iterator iterator; + + Indexes_range() = default; + Indexes_range(const std::initializer_list &l): Range(l), m_id(0), modified(true){} + Indexes_range(Range &p): Range(p), modified(true){} + Indexes_range(Range &&p): Range(std::move(p)), modified(true){} + Indexes_range(Range &p, std::size_t id): Range(p), m_id(id),modified(false){} + Indexes_range(Range &&p, std::size_t id): Range(std::move(p)), m_id(id),modified(false){} + + inline std::size_t id() const { return m_id; } + inline void set_id(std::size_t id){ m_id=id; } + inline bool was_modified() const { return modified; } + +private: + std::size_t m_id; + bool modified; +}; + +// Repair_polygon_soup without remove_pinched_polygons since our polygon are triangles +template ::Polygon_3> +struct Triangle_soup_fixer +{ + template + void operator()(PointRange& points, + PolygonRange& polygons, + const NamedParameters& np) const + { + using parameters::get_parameter; + using parameters::choose_parameter; + + using Traits = typename GetPolygonGeomTraits::type; + Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits)); + + CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup(points, polygons, np); + CGAL::Polygon_mesh_processing::internal::simplify_polygons_in_polygon_soup(points, polygons, traits); + CGAL::Polygon_mesh_processing::internal::remove_invalid_polygons_in_polygon_soup(points, polygons); + CGAL::Polygon_mesh_processing::merge_duplicate_polygons_in_polygon_soup(points, polygons, np); + CGAL::Polygon_mesh_processing::remove_isolated_points_in_polygon_soup(points, polygons); + } +}; + +// Specialization for array and Indexes_range of array +template +struct Triangle_soup_fixer > +{ + template + void operator()(PointRange& points, + PolygonRange& polygons, + const NamedParameters& np) const + { + repair_polygon_soup(points, polygons, np); + } +}; + +template +struct Triangle_soup_fixer > > +{ + template + void operator()(PointRange& points, + PolygonRange& polygons, + const NamedParameters& np) const + { + using parameters::get_parameter; + using parameters::choose_parameter; + + using Traits = typename GetPolygonGeomTraits::type; + Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits)); + + CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup(points, polygons, np); + CGAL::Polygon_mesh_processing::internal::remove_invalid_polygons_in_array_polygon_soup(points, polygons, traits); + CGAL::Polygon_mesh_processing::merge_duplicate_polygons_in_polygon_soup(points, polygons, np); + CGAL::Polygon_mesh_processing::remove_isolated_points_in_polygon_soup(points, polygons); + } +}; + +template +void repair_triangle_soup(PointRange& points, + PolygonRange& polygons, + const NamedParameters& np) +{ + Triangle_soup_fixer fixer; + fixer(points, polygons, np); +} + +// A visitor of Autorefinement to track the correspondance between input and output triangles +struct Wrap_snap_visitor : public Autorefinement::Default_visitor +{ + template< typename Triangle> + void new_subdivision(Triangle& new_t, const Triangle& old_t) { + new_t.set_id(old_t.id()); + } +}; + +/** +* +* rounds the coordinates of the points so that they fit in doubles while making and keeping the model intersection free by potentially subdividing the triangles. +* The input can be any triangle soup and the output is an intersection-free triangle soup with Hausdorff distance +* between the input and the output bounded by `M*2^-gs*k` where `M`is the maximum absolute coordinate in the model, `gs` the `snap_grid_size` (see below) and `k` the number of iterations +* performed by the algorithm. +* +* @tparam PointRange a model of the concept `RandomAccessContainer` +* whose value type is the point type +* @tparam TriangleRange a model of the concepts `RandomAccessContainer`, `BackInsertionSequence` and `Swappable`, whose +* value type is a model of the concept `RandomAccessContainer` whose value type is convertible to `std::size_t` and that +* is constructible from an `std::initializer_list` of size 3. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" +* +* @param soup_points points of the soup of polygons +* @param soup_triangles each element in the range describes a triangle using the indexed position of the points in `soup_points` +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* +* \cgalNamedParamsBegin +* \cgalParamNBegin{concurrency_tag} +* \cgalParamDescription{a tag indicating if the task should be done using one or several threads.} +* \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`} +* \cgalParamDefault{`CGAL::Sequential_tag`} +* \cgalParamNEnd +* \cgalParamNBegin{point_map} +* \cgalParamDescription{a property map associating points to the elements of the range `soup_points`} +* \cgalParamType{a model of `ReadWritePropertyMap` whose value type is a point type} +* \cgalParamDefault{`CGAL::Identity_property_map`} +* \cgalParamNEnd +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class} +* \cgalParamType{a class model of `Kernel`} +* \cgalParamDefault{a \cgal kernel deduced from the point type, using `CGAL::Kernel_traits`} +* \cgalParamExtra{The geometric traits class must be compatible with the point type.} +* \cgalParamNEnd +* \cgalParamNBegin{snap_grid_size} +* \cgalParamDescription{Scale the points to [-2^gs, 2^gs] where gs is the snap_grid_size before to round them on integer.} +* \cgalParamType{unsigned int} +* \cgalParamDefault{23} +* \cgalParamExtra{Must be lower than 52.} +* \cgalParamNEnd +* \cgalParamNBegin{number_of_iterations} +* \cgalParamDescription{Maximum number of iteration performed by the algorithm.} +* \cgalParamType{unsigned int} +* \cgalParamDefault{5} +* \cgalParamNEnd +* \cgalNamedParamsEnd +* +* \return `true` if the modified triangle soup is free from self-intersection, and `false` if the algorithm was not +* able to provide such a triangle soup within the number of iterations. +*/ +template +bool polygon_soup_snap_rounding_impl(PointRange &points, + PolygonRange &triangles, + const NamedParameters& np) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + // typedef typename GetPolygonSoupGeomTraits::type GT; + typedef typename GetPointMap::const_type Point_map; + Point_map pm = choose_parameter(get_parameter(np, internal_np::point_map)); + + typedef typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NamedParameters, + Sequential_tag + > ::type Concurrency_tag; + + // visitor + typedef typename internal_np::Lookup_named_param_def < + internal_np::visitor_t, + NamedParameters, + Autorefinement::Default_visitor + >::type Visitor; + Visitor visitor(choose_parameter(get_parameter(np, internal_np::visitor))); + + constexpr bool has_visitor = !std::is_same_v; + constexpr bool parallel_execution = std::is_same_v; + const size_t number_of_input_triangles = triangles.size(); + +#ifndef CGAL_LINKED_WITH_TBB + static_assert (!parallel_execution, + "Parallel_tag is enabled but TBB is unavailable."); +#endif + + using Point_3 = std::remove_cv_t::value_type>; + using Kernel = typename Kernel_traits::Kernel; + + // Get the grid size from the named parameter, the grid size could not be greater than 52 + const unsigned int grid_size = (std::min)(52,choose_parameter(get_parameter(np, internal_np::snap_grid_size), 23)); + const unsigned int max_nb_of_iteration = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 5); + +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "Compute the scaling of the coordinates" << std::endl; + std::cout << grid_size << std::endl; +#endif + + auto exp = [](const double v) + { + int n; + frexp(v, &n); + return n; + }; + auto pow_2 = [](const int k) + { + return (k>=0)?std::pow(2,k):1./std::pow(2,-k); + }; + + // Get max absolute value of each dimension + Bbox_3 bb = bbox_3(points.begin(), points.end()); + std::array max_abs{(std::max)(-bb.xmin(), bb.xmax()), + (std::max)(-bb.ymin(), bb.ymax()), + (std::max)(-bb.zmin(), bb.zmax())}; + // Compute scale so that the exponent of max absolute value are 52-1. + std::array scale{pow_2(grid_size - exp(max_abs[0]) - 1), + pow_2(grid_size - exp(max_abs[1]) - 1), + pow_2(grid_size - exp(max_abs[2]) - 1)}; + +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "Scaling: " << scale[0] << " " << scale[1] << " " << scale[2] << std::endl; +#endif + + auto snap = [](typename Kernel::FT x, double scale) + { + // Scale the coordinate, round to nearest integer and scale back + return double_ceil((x * scale) - 0.5) / scale; + }; + auto snap_p = [scale, snap](const Point_3 &p) + { + return Point_3( snap(p.x(),scale[0]), + snap(p.y(),scale[1]), + snap(p.z(),scale[2]) ); + }; + + for (std::size_t i = 0; i < max_nb_of_iteration; ++i) + { +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "Start Iteration " << i << std::endl; + std::cout << "Model size: " << points.size() << " " << triangles.size() << std::endl; +#endif + +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "Round all coordinates on doubles" << std::endl; +#endif + +#ifdef CGAL_LINKED_WITH_TBB + if constexpr (parallel_execution) + { + tbb::parallel_for(tbb::blocked_range(0, points.size()), + [&](const tbb::blocked_range& r){ + for(size_t pi = r.begin(); pi != r.end(); ++pi) + points[pi] = Point_3(to_double(points[pi].x()), to_double(points[pi].y()), to_double(points[pi].z())); + } + ); + } else +#endif + for (Point_3 &p : points) + p = Point_3(to_double(p.x()), to_double(p.y()), to_double(p.z())); + repair_triangle_soup(points, triangles, np); + + // Get all intersecting triangles + std::vector> pairs_of_intersecting_triangles; + triangle_soup_self_intersections(points, triangles, std::back_inserter(pairs_of_intersecting_triangles), np); + + if (pairs_of_intersecting_triangles.empty()) + { +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "End of the snapping" << std::endl; +#endif + if constexpr(has_visitor) + { + std::vector > map_io(number_of_input_triangles); + size_t id=0; + for(auto &t: triangles) + map_io[t.id()].push_back(id++); + + visitor.number_of_output_triangles(triangles.size()); + for(size_t src_id=0; src_id!=map_io.size(); ++src_id){ + if(map_io[src_id].size()==0) + visitor.delete_triangle(src_id); + else if(map_io[src_id].size()==1 && !triangles[map_io[src_id][0]].was_modified()) + visitor.verbatim_triangle_copy(map_io[src_id][0],src_id); + else + for(size_t new_id: map_io[src_id]) + visitor.new_subtriangle(new_id,src_id); + } + } + return true; + } + +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "Number of pairs of intersecting triangles: " << pairs_of_intersecting_triangles.size() << std::endl; +#endif + +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "Snap the coordinates of the vertices of the intersecting triangles" << std::endl; +#endif + +// Variants of the rounding process depending of the macros defined +#if defined(PMP_ROUNDING_VERTICES_NAIVE_SNAP) + // Nothing + +#elif defined(PMP_ROUNDING_VERTICES_CLOSEST_POINT_SNAP) + // Version where points in a voxel are rounded to the closest point + + // Group the points of the vertices of the intersecting triangles by their voxel + std::map snap_points; + std::size_t index=0; + for (auto &pair : pairs_of_intersecting_triangles) + { + for (size_t pi : triangles[pair.first]) + { + auto res=snap_points.emplace(snap_p(points[pi]), index); + if(res.second) + ++index; + } + for (size_t pi : triangles[pair.second]) + { + auto res=snap_points.emplace(snap_p(points[pi]), index); + if(res.second) + ++index; + } + } + + std::vector> identical_points(index); + for (size_t i=0; i!=points.size(); ++i) + { + Point_3 p_snap = snap_p(points[i]); + auto it=snap_points.find(p_snap); + if (it!=snap_points.end()){ + identical_points[it->second].push_back(i); + } + } + + // Replace all points in a voxel by the closest point + for(const auto &v: identical_points){ + if(v.size()>1){ + Point_3 center = snap_p(points[v[0]]); + size_t argmin(0); + for(size_t i=1; i!=v.size(); ++i){ + if(Kernel().compare_distance_3_object()(center, points[v[i]], points[v[argmin]])==SMALLER) + { + argmin=i; + } + } + for(auto i: v){ + points[i]=points[v[argmin]]; + } + } + } + +#elif defined(PMP_ROUNDING_VERTICES_BARYCENTER_SNAP) + // Version where points in a voxel are rounded to their barycenter. + + // Group the points of the vertices of the intersecting triangles by their voxel + std::map snap_points; + std::size_t index=0; + for (const auto &pair : pairs_of_intersecting_triangles) + { + for (size_t pi : triangles[pair.first]) + { + auto res=snap_points.emplace(snap_p(points[pi]), index); + if(res.second) + ++index; + } + for (size_t pi : triangles[pair.second]) + { + auto res=snap_points.emplace(snap_p(points[pi]), index); + if(res.second) + ++index; + } + } + + std::vector> identical_points(index); + for (size_t i=0; i!=points.size(); ++i) + { + Point_3 p_snap = snap_p(points[i]); + auto it=snap_points.find(p_snap); + if (it!=snap_points.end()){ + identical_points[it->second].push_back(i); + } + } + + // Replace all points in a voxel by their barycenter + for(const auto &v: identical_points){ + if(v.size()>1){ + std::array a{0,0,0}; + for(auto i: v){ + a[0]+=to_double(points[i].x()); + a[1]+=to_double(points[i].y()); + a[2]+=to_double(points[i].z()); + } + a[0]/=v.size(); + a[1]/=v.size(); + a[2]/=v.size(); + for(auto i: v){ + points[i]=Point_3(a[0],a[1],a[2]); + } + } + } +#elif defined(PMP_ROUNDING_VERTICES_FLOAT_SNAP) + // Version where points are rounded to the closest simple precision float. + + for (const auto &pair : pairs_of_intersecting_triangles) + { + for (size_t pi : triangles[pair.first]) + points[pi]=Point_3( (float) to_double(points[pi].x()), (float) to_double(points[pi].y()), (float) to_double(points[pi].z())); + for (size_t pi : triangles[pair.second]) + points[pi]=Point_3( (float) to_double(points[pi].x()), (float) to_double(points[pi].y()), (float) to_double(points[pi].z())); + } +#else // Default Version + // Version where points are rounded to the center of their voxel. + + // Get all the snap version of the points of the vertices of the intersecting triangles + // Note: points will not be modified here, they will be modified in the next for loop + + std::vector snap_points; + snap_points.reserve(pairs_of_intersecting_triangles.size() * 3); + + for (const auto &pair : pairs_of_intersecting_triangles) + { + for (size_t pi : triangles[pair.first]) + snap_points.emplace_back(snap_p(points[pi])); + for (size_t pi : triangles[pair.second]) + snap_points.emplace_back(snap_p(points[pi])); + } + +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "Snap the coordinates of the vertices close-by the previous ones" << std::endl; +#endif + + std::sort(snap_points.begin(), snap_points.end()); + snap_points.erase(std::unique(snap_points.begin(), snap_points.end()), snap_points.end()); + + // If the snapped version of a point correspond to one of the previous point, we snap it +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(parallel_execution) + { + tbb::parallel_for(tbb::blocked_range(0, points.size()), + [&](const tbb::blocked_range& r){ + for(size_t pi = r.begin(); pi != r.end(); ++pi){ + Point_3 p_snap=snap_p(points[pi]); + if (std::binary_search(snap_points.begin(), snap_points.end(), p_snap)) + points[pi] = p_snap; + } + } + ); + } else +#endif + for (Point_3 &p : points) + { + Point_3 p_snap = snap_p(p); + if (std::binary_search(snap_points.begin(), snap_points.end(), p_snap)) + p = p_snap; + } +#endif + + repair_triangle_soup(points, triangles, np); +#ifdef PMP_ROUNDING_VERTICES_IN_POLYGON_SOUP_VERBOSE + std::cout << "Model size: " << points.size() << " " << triangles.size() << std::endl; + std::cout << "Autorefine the soup" << std::endl; +#endif + if constexpr(has_visitor) + { + Wrap_snap_visitor visitor; + autorefine_triangle_soup(points, triangles, parameters::point_map(pm).concurrency_tag(Concurrency_tag()).visitor(visitor)); + } + else + { + autorefine_triangle_soup(points, triangles, parameters::point_map(pm).concurrency_tag(Concurrency_tag())); + } + } + return false; +} + +template +bool polygon_soup_snap_rounding(PointRange &soup_points, + PolygonRange &soup_triangles, + const NamedParameters& np) +{ + typedef typename internal_np::Lookup_named_param_def < + internal_np::visitor_t, + NamedParameters, + Autorefinement::Default_visitor//default + > ::type Visitor; + + constexpr bool has_visitor = !std::is_same_v; + if constexpr(has_visitor) + { + //If a visitor is provided, we color the triangles with the index of their source to correctly track the modification + using Triangle = std::remove_cv_t::value_type>; + std::vector > indexes_soup_triangles; + size_t id=0; + for(typename PolygonRange::iterator it=soup_triangles.begin(); it!=soup_triangles.end(); ++it) + indexes_soup_triangles.emplace_back((*it), id++); + + bool res=polygon_soup_snap_rounding_impl(soup_points, indexes_soup_triangles, np); + + soup_triangles.clear(); + soup_triangles.reserve(indexes_soup_triangles.size()); + for(const Indexes_range &t: indexes_soup_triangles) + soup_triangles.push_back({t[0],t[1],t[2]}); + return res; + } + else + { + return polygon_soup_snap_rounding_impl(soup_points, soup_triangles, np); + } +} + +} } } //end of CGAL::Polygon_mesh_processing::autorefine_impl namespace + +#endif //CGAL_POLYGON_MESH_PROCESSING_POLYGON_SOUP_SNAP_ROUNDING_H diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 0ed0ae8b2c7..51e2e1491a2 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -43,6 +43,7 @@ create_single_source_cgal_program("test_coref_epic_points_identity.cpp") create_single_source_cgal_program("test_does_bound_a_volume.cpp") create_single_source_cgal_program("test_pmp_clip.cpp") create_single_source_cgal_program("test_autorefinement.cpp") +create_single_source_cgal_program("test_snap_rounding.cpp") create_single_source_cgal_program("autorefinement_sm.cpp") create_single_source_cgal_program( "corefine_non_manifold.cpp" ) create_single_source_cgal_program("triangulate_hole_polyline_test.cpp") @@ -115,6 +116,7 @@ if(TARGET CGAL::TBB_support) target_link_libraries(orient_polygon_soup_test PRIVATE CGAL::TBB_support) target_link_libraries(self_intersection_surface_mesh_test PRIVATE CGAL::TBB_support) target_link_libraries(test_autorefinement PRIVATE CGAL::TBB_support) + target_link_libraries(test_snap_rounding PRIVATE CGAL::TBB_support) else() message(STATUS "NOTICE: Intel TBB was not found. Tests will use sequential code.") endif() diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_1.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_1.off new file mode 100644 index 00000000000..b0df88993a7 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_1.off @@ -0,0 +1,16 @@ +OFF +8 4 0 + +0 0 0 +0 1 0 +1 1 0 +1 0 0 +0 0 1 +0 1 1 +1 1 -1e-10 +1 0 -1e-10 +3 0 1 2 +3 0 2 3 +3 4 5 7 +3 5 6 7 + diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_2.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_2.off new file mode 100644 index 00000000000..212333059d3 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_2.off @@ -0,0 +1,19 @@ +OFF +12 4 0 + +-1 -1 -1 +-1 -1 1 +1 -1e-15 0 +-1 1 -1 +-1 1 1 +1 1e-15 0 +0 -1e-15 -1 +0 -1e-15 1 +0 -2 0 +0 1e-15 -1 +0 1e-15 1 +0 2 0 +3 0 1 2 +3 3 4 5 +3 6 7 8 +3 9 10 11 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_3.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_3.off new file mode 100644 index 00000000000..9e70e0d3303 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/collapse_3.off @@ -0,0 +1,19 @@ +OFF +12 4 0 + +-1 -1 -1 +-1 -1 1 +1 -1e-15 0 +1 1 -1 +1 1 1 +1 1e-15 0 +0 -1e-15 -1 +0 -1e-15 1 +0 -2 0 +0 1e-15 -1 +0 1e-15 1 +0 2 0 +3 0 1 2 +3 3 4 5 +3 6 7 8 +3 9 10 11 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/coplanar_cubes.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/coplanar_cubes.off new file mode 100644 index 00000000000..dfc01c0eea7 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/coplanar_cubes.off @@ -0,0 +1,202 @@ +OFF +80 120 0 +-1 -5.627423701743102e-06 1.4142135623618981 +1 -5.627423701743102e-06 1.4142135623618981 +1 -1.4142135623618981 -5.627423701743102e-06 +-1 -1.4142135623618981 -5.627423701743102e-06 +-1 5.627423701743102e-06 -1.4142135623618981 +-1 1.4142135623618981 5.627423701743102e-06 +1 1.4142135623618981 5.627423701743102e-06 +1 5.627423701743102e-06 -1.4142135623618981 +-1 -4.9339086289293424e-06 1.4142135623644878 +1 -4.9339086289293424e-06 1.4142135623644878 +1 -1.4142135623644878 -4.9339086289293424e-06 +-1 -1.4142135623644878 -4.9339086289293424e-06 +-1 4.9339086289293424e-06 -1.4142135623644878 +-1 1.4142135623644878 4.9339086289293424e-06 +1 1.4142135623644878 4.9339086289293424e-06 +1 4.9339086289293424e-06 -1.4142135623644878 +-1 -4.3072127542901049e-06 1.4142135623665355 +1 -4.3072127542901049e-06 1.4142135623665355 +1 -1.4142135623665355 -4.3072127542901049e-06 +-1 -1.4142135623665355 -4.3072127542901049e-06 +-1 4.3072127542901049e-06 -1.4142135623665355 +-1 1.4142135623665355 4.3072127542901049e-06 +1 1.4142135623665355 4.3072127542901049e-06 +1 4.3072127542901049e-06 -1.4142135623665355 +-1 -3.7949909043547426e-06 1.4142135623680028 +1 -3.7949909043547426e-06 1.4142135623680028 +1 -1.4142135623680028 -3.7949909043547426e-06 +-1 -1.4142135623680028 -3.7949909043547426e-06 +-1 3.7949909043547426e-06 -1.4142135623680028 +-1 1.4142135623680028 3.7949909043547426e-06 +1 1.4142135623680028 3.7949909043547426e-06 +1 3.7949909043547426e-06 -1.4142135623680028 +-1 -3.6070768958898969e-06 1.4142135623684944 +1 -3.6070768958898969e-06 1.4142135623684944 +1 -1.4142135623684944 -3.6070768958898969e-06 +-1 -1.4142135623684944 -3.6070768958898969e-06 +-1 3.6070768958898969e-06 -1.4142135623684944 +-1 1.4142135623684944 3.6070768958898969e-06 +1 1.4142135623684944 3.6070768958898969e-06 +1 3.6070768958898969e-06 -1.4142135623684944 +-1 -2.4904507319126939e-06 1.4142135623709018 +1 -2.4904507319126939e-06 1.4142135623709018 +1 -1.4142135623709018 -2.4904507319126939e-06 +-1 -1.4142135623709018 -2.4904507319126939e-06 +-1 2.4904507319126939e-06 -1.4142135623709018 +-1 1.4142135623709018 2.4904507319126939e-06 +1 1.4142135623709018 2.4904507319126939e-06 +1 2.4904507319126939e-06 -1.4142135623709018 +-1 0.99999891029865429 1.0000010897001574 +1 0.99999891029865429 1.0000010897001574 +1 -1.0000010897001574 0.99999891029865429 +-1 -1.0000010897001574 0.99999891029865429 +-1 -0.99999891029865429 -1.0000010897001574 +-1 1.0000010897001574 -0.99999891029865429 +1 1.0000010897001574 -0.99999891029865429 +1 -0.99999891029865429 -1.0000010897001574 +-1 1.4142135623724554 1.3448878106118537e-06 +1 1.4142135623724554 1.3448878106118537e-06 +1 -1.3448878106118537e-06 1.4142135623724554 +-1 -1.3448878106118537e-06 1.4142135623724554 +-1 -1.4142135623724554 -1.3448878106118537e-06 +-1 1.3448878106118537e-06 -1.4142135623724554 +1 1.3448878106118537e-06 -1.4142135623724554 +1 -1.4142135623724554 -1.3448878106118537e-06 +-1 1.0000000341206472 -0.99999996587935136 +1 1.0000000341206472 -0.99999996587935136 +1 0.99999996587935136 1.0000000341206472 +-1 0.99999996587935136 1.0000000341206472 +-1 -1.0000000341206472 0.99999996587935136 +-1 -0.99999996587935136 -1.0000000341206472 +1 -0.99999996587935136 -1.0000000341206472 +1 -1.0000000341206472 0.99999996587935136 +-1 1.5942993642314858e-08 -1.4142135623730947 +1 1.5942993642314858e-08 -1.4142135623730947 +1 1.4142135623730947 1.5942993642314858e-08 +-1 1.4142135623730947 1.5942993642314858e-08 +-1 -1.5942993642314858e-08 1.4142135623730947 +-1 -1.4142135623730947 -1.5942993642314858e-08 +1 -1.4142135623730947 -1.5942993642314858e-08 +1 -1.5942993642314858e-08 1.4142135623730947 +3 6 7 4 +3 3 2 1 +3 6 1 2 +3 5 0 1 +3 4 3 0 +3 7 2 3 +3 1 0 3 +3 4 5 6 +3 2 7 6 +3 1 6 5 +3 0 5 4 +3 3 4 7 +3 14 15 12 +3 11 10 9 +3 14 9 10 +3 13 8 9 +3 12 11 8 +3 15 10 11 +3 9 8 11 +3 12 13 14 +3 10 15 14 +3 9 14 13 +3 8 13 12 +3 11 12 15 +3 22 23 20 +3 19 18 17 +3 22 17 18 +3 21 16 17 +3 20 19 16 +3 23 18 19 +3 17 16 19 +3 20 21 22 +3 18 23 22 +3 17 22 21 +3 16 21 20 +3 19 20 23 +3 30 31 28 +3 27 26 25 +3 30 25 26 +3 29 24 25 +3 28 27 24 +3 31 26 27 +3 25 24 27 +3 28 29 30 +3 26 31 30 +3 25 30 29 +3 24 29 28 +3 27 28 31 +3 38 39 36 +3 35 34 33 +3 38 33 34 +3 37 32 33 +3 36 35 32 +3 39 34 35 +3 33 32 35 +3 36 37 38 +3 34 39 38 +3 33 38 37 +3 32 37 36 +3 35 36 39 +3 46 47 44 +3 43 42 41 +3 46 41 42 +3 45 40 41 +3 44 43 40 +3 47 42 43 +3 41 40 43 +3 44 45 46 +3 42 47 46 +3 41 46 45 +3 40 45 44 +3 43 44 47 +3 54 55 52 +3 51 50 49 +3 54 49 50 +3 53 48 49 +3 52 51 48 +3 55 50 51 +3 49 48 51 +3 52 53 54 +3 50 55 54 +3 49 54 53 +3 48 53 52 +3 51 52 55 +3 62 63 60 +3 59 58 57 +3 62 57 58 +3 61 56 57 +3 60 59 56 +3 63 58 59 +3 57 56 59 +3 60 61 62 +3 58 63 62 +3 57 62 61 +3 56 61 60 +3 59 60 63 +3 70 71 68 +3 67 66 65 +3 70 65 66 +3 69 64 65 +3 68 67 64 +3 71 66 67 +3 65 64 67 +3 68 69 70 +3 66 71 70 +3 65 70 69 +3 64 69 68 +3 67 68 71 +3 78 79 76 +3 75 74 73 +3 78 73 74 +3 77 72 73 +3 76 75 72 +3 79 74 75 +3 73 72 75 +3 76 77 78 +3 74 79 78 +3 73 78 77 +3 72 77 76 +3 75 76 79 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/cubes.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/cubes.off new file mode 100644 index 00000000000..493092d1a6a --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/cubes.off @@ -0,0 +1,204 @@ +OFF +80 120 0 + +-0.61594615543034026 -1.0701741879286757 -1.2146347356723264 +1.0452131429007676 -1.3737111650653515 -0.14299342950741956 +1.5883443594490732 0.52704573278436739 -0.44652546552251943 +-0.072814938882034802 0.83058270992104299 -1.5181667716874265 +-1.0452131429007676 1.3737111650653515 0.14299342950741956 +-1.5883443594490732 -0.52704573278436739 0.44652546552251943 +0.072814938882034802 -0.83058270992104299 1.5181667716874265 +0.61594615543034026 1.0701741879286757 1.2146347356723264 +-0.61594639934604423 -1.0701743934138028 -1.2146344309355146 +1.0452122820086462 -1.3737119347288513 -0.14299232819208521 +1.5883448836888012 0.52704444686667884 -0.44652511853937799 +-0.072813797665889982 0.83058198818172735 -1.518167221282807 +-1.0452122820086462 1.3737119347288513 0.14299232819208521 +-1.5883448836888012 -0.52704444686667884 0.44652511853937799 +0.072813797665889982 -0.83058198818172735 1.518167221282807 +0.61594639934604423 1.0701743934138028 1.2146344309355146 +-0.6159463762173395 -1.07017459735539 -1.2146342629779363 +1.0452117929893034 -1.3737123992443807 -0.14299144015755649 +1.5883452549553072 0.52704365666978259 -0.44652473058622827 +-0.072812914251335584 0.83058145855877352 -1.518167553406609 +-1.0452117929893034 1.3737123992443807 0.14299144015755649 +-1.5883452549553072 -0.52704365666978259 0.44652473058622827 +0.072812914251335584 -0.83058145855877352 1.518167553406609 +0.6159463762173395 1.07017459735539 1.2146342629779363 +-0.61594668912789197 -1.0701744000985252 -1.2146342780961472 +1.0452104492672072 -1.3737135467380566 -0.14299023830665158 +1.5883461040242048 0.52704182836364122 -0.44652386832759344 +-0.072811034370893465 0.83058097500317263 -1.5181679081170896 +-1.0452104492672072 1.3737135467380566 0.14299023830665158 +-1.5883461040242048 -0.52704182836364122 0.44652386832759344 +0.072811034370893465 -0.83058097500317263 1.5181679081170896 +0.61594668912789197 1.0701744000985252 1.2146342780961472 +-0.6159468868925585 -1.070174149353629 -1.2146343987318162 +1.0452090174071582 -1.3737147790282049 -0.14298886603597585 +1.5883470704382798 0.52703987181399814 -0.44652274001266495 +-0.072808833861436925 0.83058050148857432 -1.5181682727085062 +-1.0452090174071582 1.3737147790282049 0.14298886603597585 +-1.5883470704382798 -0.52703987181399814 0.44652274001266495 +0.072808833861436925 -0.83058050148857432 1.5181682727085062 +0.6159468868925585 1.070174149353629 1.2146343987318162 +-0.61594620872112704 -1.0701743005704609 -1.2146346094034179 +1.0452091168725113 -1.3737147792450597 -0.14298813688653067 +1.5883473052656241 0.52703984477504562 -0.44652193661109801 +-0.072808020328014339 0.83058032344964505 -1.5181684091279861 +-1.0452091168725113 1.3737147792450597 0.14298813688653067 +-1.5883473052656241 -0.52703984477504562 0.44652193661109801 +0.072808020328014339 -0.83058032344964505 1.5181684091279861 +0.61594620872112704 1.0701743005704609 1.2146346094034179 +0.015812704102963069 -1.4305657319567158 -0.97633582590218304 +-0.063147971057831695 -1.3970994075162435 1.0218246323460018 +1.3734324159237437 -0.0060147457399582988 1.0552948506097966 +1.4523930910845384 -0.03948107018043074 -0.94286560763838845 +0.063147971057831695 1.3970994075162435 -1.0218246323460018 +-1.3734324159237437 0.0060147457399582988 -1.0552948506097966 +-1.4523930910845384 0.03948107018043074 0.94286560763838845 +-0.015812704102963069 1.4305657319567158 0.97633582590218304 +-0.017005997503037738 -1.7115480160220284 -0.26516822000379986 +-1.4526640017493855 -0.672335677974221 0.66161320584190753 +-0.062496589720092359 0.32138030895690106 1.7008259385628288 +1.3731614145262554 -0.71783202909090671 0.77404451271712071 +1.4526640017493855 0.672335677974221 -0.66161320584190753 +0.062496589720092359 -0.32138030895690106 -1.7008259385628288 +-1.3731614145262554 0.71783202909090671 -0.77404451271712071 +0.017005997503037738 1.7115480160220284 0.26516822000379986 +-0.67677662809399264 -1.5909903432533568 0.10355251489119292 +-1.5303314706074931 -0.030328946305495721 -0.81065760039877899 +-1.0732219475484153 1.1338825933415291 0.75000327720449067 +-0.21966710503491477 -0.42677880360633264 1.6642133924944624 +1.5303314706074931 0.030328946305495721 0.81065760039877899 +1.0732219475484153 -1.1338825933415291 -0.75000327720449067 +0.21966710503491477 0.42677880360633264 -1.6642133924944624 +0.67677662809399264 1.5909903432533568 -0.10355251489119292 +-1.207106772440474 -1.2071067303442111 -0.29289346439627145 +-0.20710789493838763 -0.20710532012834212 -1.7071068233208859 +-0.50000073181889593 1.5000005345367655 -0.70710512978621964 +-1.4999996093209822 0.4999991243208966 0.70710822913839488 +0.20710789493838763 0.20710532012834212 1.7071068233208859 +0.50000073181889593 -1.5000005345367655 0.70710512978621964 +1.4999996093209822 -0.4999991243208966 -0.70710822913839488 +1.207106772440474 1.2071067303442111 0.29289346439627145 +3 6 7 4 +3 3 2 1 +3 6 1 2 +3 5 0 1 +3 4 3 0 +3 7 2 3 +3 1 0 3 +3 4 5 6 +3 2 7 6 +3 1 6 5 +3 0 5 4 +3 3 4 7 +3 14 15 12 +3 11 10 9 +3 14 9 10 +3 13 8 9 +3 12 11 8 +3 15 10 11 +3 9 8 11 +3 12 13 14 +3 10 15 14 +3 9 14 13 +3 8 13 12 +3 11 12 15 +3 22 23 20 +3 19 18 17 +3 22 17 18 +3 21 16 17 +3 20 19 16 +3 23 18 19 +3 17 16 19 +3 20 21 22 +3 18 23 22 +3 17 22 21 +3 16 21 20 +3 19 20 23 +3 30 31 28 +3 27 26 25 +3 30 25 26 +3 29 24 25 +3 28 27 24 +3 31 26 27 +3 25 24 27 +3 28 29 30 +3 26 31 30 +3 25 30 29 +3 24 29 28 +3 27 28 31 +3 38 39 36 +3 35 34 33 +3 38 33 34 +3 37 32 33 +3 36 35 32 +3 39 34 35 +3 33 32 35 +3 36 37 38 +3 34 39 38 +3 33 38 37 +3 32 37 36 +3 35 36 39 +3 46 47 44 +3 43 42 41 +3 46 41 42 +3 45 40 41 +3 44 43 40 +3 47 42 43 +3 41 40 43 +3 44 45 46 +3 42 47 46 +3 41 46 45 +3 40 45 44 +3 43 44 47 +3 54 55 52 +3 51 50 49 +3 54 49 50 +3 53 48 49 +3 52 51 48 +3 55 50 51 +3 49 48 51 +3 52 53 54 +3 50 55 54 +3 49 54 53 +3 48 53 52 +3 51 52 55 +3 62 63 60 +3 59 58 57 +3 62 57 58 +3 61 56 57 +3 60 59 56 +3 63 58 59 +3 57 56 59 +3 60 61 62 +3 58 63 62 +3 57 62 61 +3 56 61 60 +3 59 60 63 +3 70 71 68 +3 67 66 65 +3 70 65 66 +3 69 64 65 +3 68 67 64 +3 71 66 67 +3 65 64 67 +3 68 69 70 +3 66 71 70 +3 65 70 69 +3 64 69 68 +3 67 68 71 +3 78 79 76 +3 75 74 73 +3 78 73 74 +3 77 72 73 +3 76 75 72 +3 79 74 75 +3 73 72 75 +3 76 77 78 +3 74 79 78 +3 73 78 77 +3 72 77 76 +3 75 76 79 + diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_1.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_1.off new file mode 100644 index 00000000000..bd804a21a3a --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_1.off @@ -0,0 +1,14 @@ +OFF +7 3 0 + +0 0 0 +1 2 0 +1 0 0 +0 1 0 +1e-30 0 -8 +1 2 16 +2 2 2 +3 0 1 3 +3 0 2 3 +3 4 5 6 + diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_2.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_2.off new file mode 100644 index 00000000000..622155e02cb --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_2.off @@ -0,0 +1,15 @@ +OFF +9 3 0 + +0 1 0 +0 0 0 +2 1 0 +1 1 0 +2 0 0 +1.00000000000002 0 0 +1.00000000000001 1e-30 0 +2 -1 0 +1.00000000000002 -1 0 +3 0 1 2 +3 3 4 5 +3 6 7 8 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_3.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_3.off new file mode 100644 index 00000000000..9cd7437d4e0 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/intersection_when_rounded_3.off @@ -0,0 +1,19 @@ +OFF +12 4 0 + +0 0 0 +0 8.000000000001 0 +8.000000000001 0 0 +4.000000000002 4.000000000002 0 +4 8.000000000001 0 +8.000000000001 4 0 +6.000000000002 6.000000000002 0 +6 8.000000000001 0 +8.000000000001 6 0 +6.5 6.5 0 +6.5 8.000000000001 0 +8.000000000001 6.5 0 +3 0 1 2 +3 3 4 5 +3 6 7 8 +3 9 10 11 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/no_collapse_1.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/no_collapse_1.off new file mode 100644 index 00000000000..62b413a8497 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/no_collapse_1.off @@ -0,0 +1,16 @@ +OFF +8 4 0 + +0 0 0 +0 1 0 +1 1 0 +1 0 0 +0 0 1 +0 1 1 +1 1 1e-10 +1 0 1e-10 +3 0 1 2 +3 0 2 3 +3 4 5 7 +3 5 6 7 + diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/no_collapse_2.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/no_collapse_2.off new file mode 100644 index 00000000000..754d0a50ba0 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data-snap/no_collapse_2.off @@ -0,0 +1,19 @@ +OFF +12 4 0 + +1 -1 -1 +1 -1 1 +3 -1e-15 0 +1 1 -1 +1 1 1 +3 1e-15 0 +0 -1e-15 -1 +0 -1e-15 1 +0 -2 0 +0 1e-15 -1 +0 1e-15 1 +0 2 0 +3 0 1 2 +3 3 4 5 +3 6 7 8 +3 9 10 11 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp index e04e07ba556..e67dc15a4d8 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp @@ -30,7 +30,7 @@ struct My_exp_visitor : std::shared_ptr i; }; -struct My_visitor +struct My_visitor : public PMP::Autorefinement::Default_visitor { My_visitor(std::size_t nb_input, std::size_t expected_nb_output) : nb_input(nb_input) @@ -66,6 +66,11 @@ struct My_visitor tgt_check[tgt_id]=1; } + void delete_triangle(std::size_t src_id) + { + assert(src_id tgt_check; diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding.cmd b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding.cmd new file mode 100644 index 00000000000..2d9a58e2f6e --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding.cmd @@ -0,0 +1,37 @@ +data-autoref/test_01.off 4 2 +data-autoref/test_02.off 10 17 +data-autoref/test_03.off 13 20 +data-autoref/test_04.off 12 20 +data-autoref/test_05.off 12 20 +data-autoref/test_06.off 89 248 +data-autoref/test_07.off 8 12 +data-autoref/test_08.off 57 148 +data-autoref/test_09.off 4 3 +data-autoref/test_10.off 10 13 +data-autoref/test_11.off 9 12 +data-autoref/test_12.off 9 6 +data-autoref/test_13.off 12 22 +data-autoref/test_14.off 12 22 +data-autoref/test_15.off 12 22 +data-autoref/test_16.off 4 2 +data-autoref/test_17.off 4 2 +data-autoref/triple_inter_exception/triple.off 15 18 +data-autoref/triple_inter_exception/cubes_cpln_1.off 66 206 +data-autoref/triple_inter_exception/cubes_cpln_2.off 54 158 +data-autoref/triple_inter_exception/cubes_cpln_3.off 61 188 +data-autoref/open_01.off 1313 2622 +data-autoref/open_02.off 562 1056 +data-autoref/cpln_01.off 30 78 +data-autoref/cpln_02.off 24 56 +data-autoref/cpln_03.off 27 68 +data-autoref/four_cubes.off 106 352 +data-autoref/spiral.off 19 31 +data-snap/intersection_when_rounded_1.off 9 6 +data-snap/intersection_when_rounded_2.off 10 7 +data-snap/intersection_when_rounded_3.off 16 14 +data-snap/collapse_1.off 6 4 +data-snap/collapse_2.off 13 16 +data-snap/collapse_3.off 11 10 +data-snap/no_collapse_1.off 8 4 +data-snap/no_collapse_2.off 12 4 +data-snap/coplanar_cubes.off 1514 5544 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding.cpp new file mode 100644 index 00000000000..b8d8689a752 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding.cpp @@ -0,0 +1,100 @@ + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +namespace PMP = CGAL::Polygon_mesh_processing; + +struct My_visitor : public PMP::Autorefinement::Default_visitor +{ + My_visitor(std::size_t nb_input, std::size_t expected_nb_output) + : nb_input(nb_input) + , expected_nb_output(expected_nb_output) + {} + + ~My_visitor() + { + for(std::size_t i=0; i tgt_check; +}; + +template +void test(const char* fname, std::size_t nb_vertices_after_autorefine, std::size_t expected_nb_output, Tag tag) +{ + std::cout << "Running tests on " << fname; + if (std::is_same_v) + std::cout << " (Sequential)\n"; + else + std::cout << " (Parallel)\n"; + + std::vector points; + std::vector< std::vector > triangles; + if (!CGAL::IO::read_polygon_soup(fname, points, triangles)) + { + std::cerr << " Input mesh is not a valid file." << std::endl; + exit(EXIT_FAILURE); + } + +// Testing autorefine() + My_visitor visitor(triangles.size(), expected_nb_output); + PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::visitor(visitor).concurrency_tag(tag).apply_iterative_snap_rounding(true)); + assert( nb_vertices_after_autorefine==points.size()); + assert( expected_nb_output==triangles.size()); + assert( !PMP::does_triangle_soup_self_intersect(points, triangles) ); +// CGAL::IO::write_polygon_soup("/tmp/debug.off", points, triangles); +} + +int main(int argc, const char** argv) +{ + // file expected_nb_of_vertices expected_nb_of_faces (after repair) + for (int i=0;i<(argc-1)/3; ++i) + { + test(argv[1+3*i], atoi(argv[1+3*i+1]), atoi(argv[1+3*i+2]), CGAL::Sequential_tag()); +#ifdef CGAL_LINKED_WITH_TBB + test(argv[1+3*i], atoi(argv[1+3*i+1]), atoi(argv[1+3*i+2]), CGAL::Parallel_tag()); +#endif + } +} diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding_full.cmd b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding_full.cmd new file mode 100644 index 00000000000..49eb20e9a00 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_snap_rounding_full.cmd @@ -0,0 +1,38 @@ +data-autoref/test_01.off 4 2 +data-autoref/test_02.off 10 17 +data-autoref/test_03.off 13 20 +data-autoref/test_04.off 12 20 +data-autoref/test_05.off 12 20 +data-autoref/test_06.off 89 248 +data-autoref/test_07.off 8 12 +data-autoref/test_08.off 57 148 +data-autoref/test_09.off 4 3 +data-autoref/test_10.off 10 13 +data-autoref/test_11.off 9 12 +data-autoref/test_12.off 9 6 +data-autoref/test_13.off 12 22 +data-autoref/test_14.off 12 22 +data-autoref/test_15.off 12 22 +data-autoref/test_16.off 4 2 +data-autoref/test_17.off 4 2 +data-autoref/triple_inter_exception/triple.off 15 18 +data-autoref/triple_inter_exception/cubes_cpln_1.off 66 206 +data-autoref/triple_inter_exception/cubes_cpln_2.off 54 158 +data-autoref/triple_inter_exception/cubes_cpln_3.off 61 188 +data-autoref/open_01.off 1313 2622 +data-autoref/open_02.off 562 1056 +data-autoref/cpln_01.off 30 78 +data-autoref/cpln_02.off 24 56 +data-autoref/cpln_03.off 27 68 +data-autoref/four_cubes.off 106 352 +data-autoref/spiral.off 19 31 +data-snap/intersection_when_rounded_1.off 9 6 +data-snap/intersection_when_rounded_2.off 10 7 +data-snap/intersection_when_rounded_3.off 16 14 +data-snap/collapse_1.off 6 4 +data-snap/collapse_2.off 13 16 +data-snap/collapse_3.off 11 10 +data-snap/no_collapse_1.off 8 4 +data-snap/no_collapse_2.off 12 4 +data-snap/cubes.off 5042 26036 +data-snap/coplanar_cubes.off 1514 5544 diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index 8adf0d25f8b..ffe386d086c 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -177,7 +177,8 @@ CGAL_add_named_parameter(weight_limiting_t, weight_limiting, weight_limiting) CGAL_add_named_parameter(progressive_t, progressive, progressive) CGAL_add_named_parameter(tiling_t, tiling, tiling) CGAL_add_named_parameter(dimension_t, dimension, dimension) - +CGAL_add_named_parameter(apply_iterative_snap_rounding_t, apply_iterative_snap_rounding, apply_iterative_snap_rounding) +CGAL_add_named_parameter(snap_grid_size_t, snap_grid_size, snap_grid_size) // List of named parameters that we use in the package 'Surface Mesh Simplification' CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)