diff --git a/Alpha_wrap_3/benchmark/Alpha_wrap_3/CMakeLists.txt b/Alpha_wrap_3/benchmark/Alpha_wrap_3/CMakeLists.txt new file mode 100644 index 00000000000..a9aa0d1d63c --- /dev/null +++ b/Alpha_wrap_3/benchmark/Alpha_wrap_3/CMakeLists.txt @@ -0,0 +1,15 @@ +# Created by the script cgal_create_cmake_script +# This is the CMake script for compiling a CGAL application. + +cmake_minimum_required(VERSION 3.1...3.20) +project(Alpha_wrap_3_Benchmark) + +find_package(CGAL REQUIRED) + +include_directories (BEFORE ../../include ./Quality ./Robustness) # AW3 includes +include_directories (BEFORE ../../../CGAL-Patches/include) + +# create a target per cppfile +create_single_source_cgal_program("Performance/performance_benchmark.cpp") +create_single_source_cgal_program("Quality/quality_benchmark.cpp") +create_single_source_cgal_program("Robustness/robustness_benchmark.cpp") diff --git a/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/compute_performance_benchmark_data.py b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/compute_performance_benchmark_data.py new file mode 100644 index 00000000000..86c57d35146 --- /dev/null +++ b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/compute_performance_benchmark_data.py @@ -0,0 +1,61 @@ +# Copyright (c) 2019-2023 Google LLC (USA). +# All rights reserved. +# +# This file is part of CGAL (www.cgal.org). +# +# $URL$ +# $Id$ +# SPDX-License-Identifier: GPL-3.0-or-later +# +# +# Author(s) : Pierre Alliez +# Michael Hemmer +# Cedric Portaneri +# +#!/usr/bin/python + +import os, sys, subprocess, datetime, time, getopt + +def compute_performance_benchmark_data(execname, filename, alpha): + + output = "" + cmd = ("/usr/bin/time", "-v", + execname, "-i", + filename, "-a", alpha) + proc = subprocess.Popen( + cmd, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + start_new_session=True) + + outs, errs = proc.communicate() + output = outs.decode("utf-8") + errs.decode("utf-8") + + for output_line in output.split("\n"): + if "User time (seconds): " in output_line: + print(output_line[len("User time (seconds): "):]) + continue + if "Maximum resident set size (kbytes): " in output_line: + print(output_line[len("Maximum resident set size (kbytes): "):]) + continue + +def main(argv): + execname="" + filename="" + alpha="" + try: + opts, args = getopt.getopt(sys.argv[1:], 'e:i:a:') + except getopt.GetoptError: + sys.exit(2) + for opt, arg in opts: + if opt == "-e": + execname = arg + elif opt == "-i": + filename = arg + elif opt == "-a": + alpha = arg + + compute_performance_benchmark_data(execname, filename, alpha) + +if __name__ == "__main__": + main(sys.argv[1:]) diff --git a/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/generate_performance_benchmark_charts.py b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/generate_performance_benchmark_charts.py new file mode 100644 index 00000000000..445b6923277 --- /dev/null +++ b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/generate_performance_benchmark_charts.py @@ -0,0 +1,156 @@ +# Copyright (c) 2019-2023 Google LLC (USA). +# All rights reserved. +# +# This file is part of CGAL (www.cgal.org). +# +# $URL$ +# $Id$ +# SPDX-License-Identifier: GPL-3.0-or-later +# +# +# Author(s) : Pierre Alliez +# Michael Hemmer +# Cedric Portaneri +# +#!/usr/bin/python + +import os, sys, subprocess, datetime, time, signal, getopt +import numpy as np +import matplotlib.pyplot as plt + +def main(argv): + + inputdir="" + outputdir="" + commit_hash="" + alpha="" + do_diff=False + diffdir="" + diff_hash="" + try: + opts, args = getopt.getopt(sys.argv[1:], 'i:a:o:c:d:p:') + except getopt.GetoptError: + sys.exit(2) + for opt, arg in opts: + if opt == "-i": + inputdir = arg + elif opt == "-a": + alpha = arg + elif opt == "-o": + outputdir = arg + elif opt == "-c": + commit_hash = arg + elif opt == "-d": + diff_hash = arg + do_diff = True + elif opt == "-p": + diffdir = arg + + all_metric = { + "Time_(second)" : {}, + "Memory_Peak_(kbytes)" : {}} + num_input = 0 + for filename in os.listdir(inputdir) : + new_path = os.path.join(inputdir,filename) + new_file = open(new_path) + is_empty_new = os.path.getsize(new_path) <= 1 + if do_diff : + old_path = os.path.join(diffdir,filename) + old_file = open(old_path) + is_empty_old = os.path.getsize(old_path) <= 1 + for key in all_metric: + if is_empty_new or is_empty_old : + new_val = 0. + old_val = 0. + else : + new_val = float(new_file.readline().rstrip('\n')) + old_val = float(old_file.readline().rstrip('\n')) + mesh_id = str(filename.split('.')[0]) + all_metric[key][mesh_id] = [new_val, old_val] + else : + for key in all_metric: + if is_empty_new : + new_val = 0. + else : + new_val = float(new_file.readline().rstrip('\n')) + mesh_id = str(filename.split('.')[0]) + all_metric[key][mesh_id] = [new_val, new_val] + num_input = num_input+1 + + # update .pdf chart + date_now = datetime.datetime.now() + date_for_filename = str(date_now.year) +"_"+ str(date_now.month) +"_"+ str(date_now.day) +"_"+ str(date_now.hour) +"h"+ str(date_now.minute) +"mn" + for key in all_metric: + goal = 0 + num_el = range(len(all_metric[key])) + avg_diff_to_goal = 0. + avg = 0. + x1 = [] + x2 = [] + for value in all_metric[key].values() : + avg += value[0] + diff_to_goal = abs(value[1]-goal) - abs(value[0]-goal) + avg_diff_to_goal += diff_to_goal + x1.append(value[0]) + x2.append(value[1]) + avg_diff_to_goal /= float(len(all_metric[key])) + avg /= float(len(all_metric[key])) + + plt.figure(figsize=(8,8)) + if do_diff : + plt.hist(x2, bins=100, color='tab:green', alpha=0.5) + plt.hist(x1, bins=100, color='tab:blue', alpha=0.5) + plt.vlines(x = goal, ymin=plt.ylim()[0], ymax=plt.ylim()[1], linestyles='dashed') + + title = "" + if do_diff : + title += "Diff between " + commit_hash + " and " + diff_hash + " on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha + else : + title += "Benchmarking on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha + + avg_str = str(format(abs(avg), '.2f')) + if key == "Time_(second)" : + title += "\nIn average we spend " + avg_str + " seconds" + else : + title += "\nIn average we use up to " + avg_str + " kbytes" + + if do_diff and avg_diff_to_goal == 0. : + title += "\nNo change between the two commits" + elif do_diff : + avg_diff_str = str(format(abs(avg_diff_to_goal), '.2f')) + if key == "Time_(second)" : + if avg_diff_to_goal < 0 : + title += "\nIn average we get slower by " + else : + title += "\nIn average we get faster " + title += avg_diff_str + " seconds" + else : + if avg_diff_to_goal < 0 : + title += "\nIn average we use " + avg_diff_str + " more" + else : + title += "\nIn average we use " + avg_diff_str + " less" + title += " kbytes" + + plt.title(title, fontsize=15) + plt.xlabel(key.replace("_"," "), fontsize=14) + plt.ylabel("# of meshes", fontsize=14) + plt.tick_params(axis="x", labelsize=9) + plt.tick_params(axis="y", labelsize=9) + + chart_filename = "" + if do_diff : + chart_filename += "diff_"+commit_hash+"_"+diff_hash+"_"+key+"_"+date_for_filename+".pdf" + else : + chart_filename += "results_"+commit_hash+"_"+key+"_"+date_for_filename+".pdf" + chart_path = os.path.join(outputdir+"/charts",chart_filename) + if os.path.isfile(chart_path) : + os.remove(chart_path) + plt.savefig(chart_path, bbox_inches="tight") + plt.close() + + print("pdf updated") + + sys.exit() + +if __name__ == "__main__": + main(sys.argv[1:]) diff --git a/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/performance_benchmark.cpp b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/performance_benchmark.cpp new file mode 100644 index 00000000000..207b4704635 --- /dev/null +++ b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/performance_benchmark.cpp @@ -0,0 +1,65 @@ +#include +#include + +#include + +#include + +#include +#include +#include +#include + +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using Point_3 = K::Point_3; +using Vector_3 = K::Vector_3; + +using Mesh = CGAL::Surface_mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char** argv) +{ + const int argc_check = argc - 1; + const char* entry_name_ptr = nullptr; + double relative_alpha_ratio = 20., relative_offset_ratio = 600.; + + for(int i=1; i points; + std::vector > faces; + if(!CGAL::IO::read_polygon_soup(entry_name_ptr, points, faces) || faces.empty()) + { + std::cerr << "Error: Invalid input data." << std::endl; + return EXIT_FAILURE; + } + + CGAL::Bbox_3 bbox; + for(const Point_3& p : points) + bbox += p.bbox(); + + const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) + + CGAL::square(bbox.ymax() - bbox.ymin()) + + CGAL::square(bbox.zmax() - bbox.zmin())); + const double alpha = diag_length / relative_alpha_ratio; + const double offset = diag_length / relative_offset_ratio; + + Mesh wrap; + CGAL::alpha_wrap_3(points, faces, alpha, offset, wrap); + + return EXIT_SUCCESS; +} diff --git a/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/compute_quality_benchmark_data.py b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/compute_quality_benchmark_data.py new file mode 100644 index 00000000000..3b996cb5749 --- /dev/null +++ b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/compute_quality_benchmark_data.py @@ -0,0 +1,54 @@ +# Copyright (c) 2019-2023 Google LLC (USA). +# All rights reserved. +# +# This file is part of CGAL (www.cgal.org). +# +# $URL$ +# $Id$ +# SPDX-License-Identifier: GPL-3.0-or-later +# +# +# Author(s) : Pierre Alliez +# Michael Hemmer +# Cedric Portaneri +# +#!/usr/bin/python + +import os, sys, subprocess, datetime, time, getopt + +def compute_quality_benchmark_data(execname, filename, alpha): + + output = "" + cmd = (execname, "-i", + filename, "-a", alpha) + proc = subprocess.Popen( + cmd, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + start_new_session=True) + + outs, errs = proc.communicate() + output = outs.decode("utf-8") + errs.decode("utf-8") + + print(output) + +def main(argv): + execname="" + filename="" + alpha="" + try: + opts, args = getopt.getopt(sys.argv[1:], 'e:i:a:') + except getopt.GetoptError: + sys.exit(2) + for opt, arg in opts: + if opt == "-e": + execname = arg + elif opt == "-i": + filename = arg + elif opt == "-a": + alpha = arg + + compute_quality_benchmark_data(execname, filename, alpha) + +if __name__ == "__main__": + main(sys.argv[1:]) diff --git a/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/distance_utils.h b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/distance_utils.h new file mode 100644 index 00000000000..379573e9c90 --- /dev/null +++ b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/distance_utils.h @@ -0,0 +1,151 @@ +// Copyright (c) 2019-2022 Google LLC (USA). +// 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) : Pierre Alliez +// Michael Hemmer +// Cedric Portaneri + +#ifndef CGAL_ALPHA_WRAP_3_BENCHMARK_ALPHA_WRAP_3_QUALITY_DISTANCE_H_ +#define CGAL_ALPHA_WRAP_3_BENCHMARK_ALPHA_WRAP_3_QUALITY_DISTANCE_H_ + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Aw3i { + +enum Distance_metric { HAUSDORFF = 0, MEAN = 1, RMS = 2 }; + +template +inline double approximate_hausdorff_distance(const std::vector& sample_points, + const AABBTree& tree, + Point& hint) +{ + double hdist = 0; + for(const Point& pt : sample_points) + { + hint = tree.closest_point(pt, hint); + auto dist = CGAL::squared_distance(hint, pt); + double d = CGAL::to_double(CGAL::approximate_sqrt(dist)); + if(d > hdist) + hdist = d; + } + + return hdist; +} + +template +inline double approximate_mean_distance(const std::vector& sample_points, + const AABBTree& tree, + Point& hint) +{ + double mdist = 0; + for(const Point& pt : sample_points) + { + hint = tree.closest_point(pt, hint); + auto dist = CGAL::squared_distance(hint, pt); + double d = CGAL::to_double(CGAL::approximate_sqrt(dist)); + mdist += d; + } + + return mdist / sample_points.size(); +} + +template +inline double approximate_rms_distance(const std::vector& sample_points, + const AABBTree& tree, + Point& hint) +{ + double rmsdist = 0; + for(const Point& pt : sample_points) + { + hint = tree.closest_point(pt, hint); + auto dist = CGAL::squared_distance(hint, pt); + rmsdist += CGAL::to_double(dist); + } + + return CGAL::to_double(CGAL::approximate_sqrt(rmsdist / sample_points.size())); +} + +template +inline double approximate_distance(const TriangleMesh& tm1, + const TriangleMesh& tm2, + const Distance_metric& metric) +{ + using GT = typename CGAL::GetGeomTraits::type; + using Point_3 = typename GT::Point_3; + + using Primitive = CGAL::AABB_face_graph_triangle_primitive; + using AABB_traits = CGAL::AABB_traits; + using AABB_tree = CGAL::AABB_tree; + + using CGAL::parameters::choose_parameter; + using CGAL::parameters::get_parameter; + + std::vector original_sample_points; + CGAL::Polygon_mesh_processing::sample_triangle_mesh(tm1, std::back_inserter(original_sample_points), + CGAL::parameters::all_default()); + + std::vector sample_points(std::begin(original_sample_points), + std::end(original_sample_points)); + CGAL::spatial_sort(sample_points.begin(), sample_points.end()); + + AABB_tree tree(faces(tm2).first, faces(tm2).second, tm2); + tree.build(); + + auto vpm_2 = get(CGAL::vertex_point, tm2); + Point_3 hint = get(vpm_2, *vertices(tm2).first); + + if(metric == HAUSDORFF) + return approximate_hausdorff_distance(sample_points, tree, hint); + else if(metric == MEAN) + return approximate_mean_distance(sample_points, tree, hint); + else if(metric == RMS) + return approximate_rms_distance(sample_points, tree, hint); + else + std::cerr << "Metric unknown\n" << std::endl; + + return -1.0; +} + +template +double get_longest_diag_bbox(const TriangleMesh& tm) +{ + CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(tm); + return std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) + + CGAL::square(bbox.ymax() - bbox.ymin()) + + CGAL::square(bbox.zmax() - bbox.zmin())); +} + +template +inline double approximate_distance_relative_to_bbox(const TriangleMesh& tm1, + const TriangleMesh& tm2, + const Distance_metric& metric) +{ + double longest_diag_length = get_longest_diag_bbox(tm1); + return approximate_distance(tm1, tm2, metric) / longest_diag_length; +} + +template +inline double approximate_distance_relative_to_bbox(const TriangleMesh& tm1, + const TriangleMesh& tm2, + const Distance_metric& metric, + const FT& longest_diag_length) +{ + return approximate_distance(tm1, tm2, metric) / CGAL::to_double(longest_diag_length); +} + +} // namespace Aw3i + +#endif // CGAL_CGAL_ALPHA_WRAP_3_BENCHMARK_ALPHA_WRAP_3_QUALITY_DISTANCE_H_ diff --git a/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/generate_quality_benchmark_charts.py b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/generate_quality_benchmark_charts.py new file mode 100644 index 00000000000..0df5ed5e0f4 --- /dev/null +++ b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/generate_quality_benchmark_charts.py @@ -0,0 +1,182 @@ +# Copyright (c) 2019-2023 Google LLC (USA). +# All rights reserved. +# +# This file is part of CGAL (www.cgal.org). +# +# $URL$ +# $Id$ +# SPDX-License-Identifier: GPL-3.0-or-later +# +# +# Author(s) : Pierre Alliez +# Michael Hemmer +# Cedric Portaneri +# +#!/usr/bin/python + +import os, sys, subprocess, datetime, time, signal, getopt +import numpy as np +import matplotlib.pyplot as plt + +def main(argv): + + inputdir="" + outputdir="" + commit_hash="" + alpha="" + do_diff=False + diffdir="" + diff_hash="" + try: + opts, args = getopt.getopt(sys.argv[1:], 'i:a:o:c:d:p:') + except getopt.GetoptError: + sys.exit(2) + for opt, arg in opts: + if opt == "-i": + inputdir = arg + elif opt == "-a": + alpha = arg + elif opt == "-o": + outputdir = arg + elif opt == "-c": + commit_hash = arg + elif opt == "-d": + diff_hash = arg + do_diff = True + elif opt == "-p": + diffdir = arg + + all_metric = { + "Mean_Min_Angle_(degree)" : {}, + "Mean_Max_Angle_(degree)" : {}, + "Mean_Radius_Ratio" : {}, + "Mean_Edge_Ratio" : {}, + "Mean_Aspect_Ratio" : {}, + "Complexity_(#_of_triangle)" : {}, + "#_of_almost_degenerate_triangle" : {}, + "Hausdorff_distance_output_to_input_(%_of_bbox_diag)" : {}} + num_input = 0 + print("inputdir = ", inputdir) + for filename in os.listdir(inputdir) : + new_path = os.path.join(inputdir,filename) + new_file = open(new_path) + if do_diff : + old_path = os.path.join(diffdir,filename) + old_file = open(old_path) + is_empty_old = os.path.getsize(old_path) <= 1 + for key in all_metric : + try : + new_val = float(new_file.readline().rstrip('\n')) + old_val = float(old_file.readline().rstrip('\n')) + mesh_id = str(filename.split('.')[0]) + all_metric[key][mesh_id] = [new_val, old_val] + except ValueError: + pass + else : + for key in all_metric : + try : + new_val = float(new_file.readline().rstrip('\n')) + mesh_id = str(filename.split('.')[0]) + all_metric[key][mesh_id] = [new_val, new_val] + except ValueError: + pass + num_input = num_input+1 + + # update .pdf chart + date_now = datetime.datetime.now() + date_for_filename = str(date_now.year) +"_"+ str(date_now.month) +"_"+ str(date_now.day) +"_"+ str(date_now.hour) +"h"+ str(date_now.minute) +"mn" + for key in all_metric: + goal = 0 + if key == "Mean_Min_Angle_(degree)" or key == "Mean_Max_Angle_(degree)": + goal = 60 + elif key == "Mean_Radius_Ratio" or key == "Mean_Edge_Ratio" or key == "Mean_Aspect_Ratio" : + goal = 1 + + num_el = range(len(all_metric[key])) + avg_diff_to_goal = 0. + avg = 0. + x1 = [] + x2 = [] + for value in all_metric[key].values() : + avg += value[0] + diff_to_goal = abs(value[1]-goal) - abs(value[0]-goal) + avg_diff_to_goal += diff_to_goal + x1.append(value[0]) + x2.append(value[1]) + avg_diff_to_goal /= float(len(all_metric[key])) + avg /= float(len(all_metric[key])) + + plt.figure(figsize=(8,8)) + if do_diff : + plt.hist(x2, bins=100, color='tab:green', alpha=0.5) + plt.hist(x1, bins=100, color='tab:blue', alpha=0.5) + plt.vlines(x = goal, ymin=plt.ylim()[0], ymax=plt.ylim()[1], linestyles='dashed') + + title = "" + if do_diff : + title += "Diff between " + commit_hash + " and " + diff_hash + " on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha + else : + title += "Benchmarking on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha + + avg_str = str(format(abs(avg), '.2f')) + if key == "Mean_Min_Angle_(degree)" or key == "Mean_Max_Angle_(degree)": + title += "\nIn average we have " + avg_str + "°" + elif key == "Mean_Radius_Ratio" or key == "Mean_Edge_Ratio" or key == "Mean_Aspect_Ratio" : + title += "\nIn average we have a ratio of " + avg_str + elif key == "Hausdorff_distance_output_to_input_(%_of_bbox_diag)" : + title += "\nIn average we have a distance of " + avg_str + "% of bbox diag" + elif key == "Complexity_(#_of_triangle)" or key == "#_of_almost_degenerate_triangle" : + title += "\nIn average we have " + avg_str + " triangles" + + if do_diff and avg_diff_to_goal == 0. : + title += "\nNo change between the two commits" + elif do_diff : + avg_diff_str = str(format(abs(avg_diff_to_goal), '.2f')) + if key == "Mean_Min_Angle_(degree)" or key == "Mean_Max_Angle_(degree)": + if avg_diff_to_goal < 0 : + title += "\nIn average we loose " + else : + title += "\nIn average we gain " + title += avg_diff_str + "° toward 60°" + elif key == "Mean_Radius_Ratio" or key == "Mean_Edge_Ratio" or key == "Mean_Aspect_Ratio" : + if avg_diff_to_goal < 0 : + title += "\nIn average we loose " + else : + title += "\nIn average we gain " + title += avg_diff_str + " of ratio toward 1" + elif key == "Hausdorff_distance_output_to_input_(%_of_bbox_diag)" : + if avg_diff_to_goal < 0 : + title += "\nIn average we increase by " + else : + title += "\nIn average we reduce by " + title += avg_diff_str + " the bbox ratio" + elif key == "Complexity_(#_of_triangle)" or key == "#_of_almost_degenerate_triangle" : + if avg_diff_to_goal < 0 : + title += "\nIn average we get " + avg_diff_str + " more" + else : + title += "\nIn average we get " + avg_diff_str + " less" + title += " triangles" + + plt.title(title, fontsize=15) + plt.xlabel(key.replace("_"," "), fontsize=14) + plt.ylabel("# of meshes", fontsize=14) + plt.tick_params(axis="x", labelsize=9) + plt.tick_params(axis="y", labelsize=9) + + chart_filename = "" + if do_diff : + chart_filename += "diff_"+commit_hash+"_"+diff_hash+"_"+key+"_"+date_for_filename+".pdf" + else : + chart_filename += "results_"+commit_hash+"_"+key+"_"+date_for_filename+".pdf" + chart_path = os.path.join(outputdir+"/charts",chart_filename) + if os.path.isfile(chart_path) : + os.remove(chart_path) + plt.savefig(chart_path, bbox_inches="tight") + plt.close() + + print("pdf updated") + + sys.exit() + +if __name__ == "__main__": + main(sys.argv[1:]) diff --git a/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/quality_benchmark.cpp b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/quality_benchmark.cpp new file mode 100644 index 00000000000..f8e54c488a4 --- /dev/null +++ b/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/quality_benchmark.cpp @@ -0,0 +1,271 @@ +#include + +#include +#include + +#include + +#include + +#include +#include +#include + +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Point_3 = Kernel::Point_3; +using Vector_3 = Kernel::Vector_3; +using Triangle_3 = Kernel::Triangle_3; +using FT = Kernel::FT; + +using Mesh = CGAL::Surface_mesh; +using face_descriptor = boost::graph_traits::face_descriptor; + +using Oracle = CGAL::Alpha_wraps_3::internal::Triangle_mesh_oracle; +using Dt = CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3::Triangulation; + +namespace PMP = CGAL::Polygon_mesh_processing; + +std::array triangle_angles(const Triangle_3& tr) +{ + FT sq_a = CGAL::squared_distance(tr[0], tr[1]); + FT sq_b = CGAL::squared_distance(tr[1], tr[2]); + FT sq_c = CGAL::squared_distance(tr[2], tr[0]); + + FT two_ab = 2. * CGAL::sqrt(sq_a) * CGAL::sqrt(sq_b); + FT two_bc = 2. * CGAL::sqrt(sq_b) * CGAL::sqrt(sq_c); + FT two_ca = 2. * CGAL::sqrt(sq_c) * CGAL::sqrt(sq_a); + + FT angle_a = (sq_b + sq_c - sq_a) / two_bc; + FT angle_b = (sq_c + sq_a - sq_b) / two_ca; + FT angle_c = (sq_a + sq_b - sq_c) / two_ab; + if(angle_a < -1.) angle_a = -1.; + if(angle_b < -1.) angle_b = -1.; + if(angle_c < -1.) angle_c = -1.; + if(angle_a > 1.) angle_a = 1.; + if(angle_b > 1.) angle_b = 1.; + if(angle_c > 1.) angle_c = 1.; + angle_a = std::acos(angle_a); + angle_b = std::acos(angle_b); + angle_c = std::acos(angle_c); + + return {angle_a, angle_b, angle_c}; +} + +bool is_almost_degenerate(const Triangle_3& tr, + double threshold) +{ + FT sq_area = tr.squared_area(); + return (CGAL::sqrt(CGAL::to_double(sq_area)) < threshold); +} + +auto surface_mesh_face_to_triangle(const face_descriptor fd, + const Mesh& sm) +{ + typename boost::graph_traits::halfedge_descriptor hd = halfedge(fd,sm); + return Triangle_3(sm.point(target(hd,sm)), + sm.point(target(next(hd,sm),sm)), + sm.point(target(next(next(hd,sm),sm),sm))); +} + +double mean_min_angle(const Mesh& mesh) +{ + double mean_min_angle = 0.; + for(const face_descriptor f : faces(mesh)) + { + const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh); + std::array angles = triangle_angles(tr); + + FT min_angle = std::min({angles[0], angles[1], angles[2]}); + + min_angle = min_angle * (180.0 / CGAL_PI); + mean_min_angle += min_angle; + } + + mean_min_angle /= static_cast(mesh.number_of_faces()); + return mean_min_angle; +} + +double mean_max_angle(const Mesh& mesh) +{ + double mean_max_angle = 0.; + for(const face_descriptor f : faces(mesh)) + { + const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh); + std::array angles = triangle_angles(tr); + + FT max_angle = std::max({angles[0], angles[1], angles[2]}); + + max_angle = max_angle * (180.0 / CGAL_PI); + mean_max_angle += max_angle; + } + + mean_max_angle /= static_cast(mesh.number_of_faces()); + return mean_max_angle; +} + +double mean_radius_ratio(const Mesh& mesh, + double degenerate_threshold) +{ + double mean_radius_ratio = 0.; + size_t num_almost_degenerate_tri = 0; + for(const face_descriptor f : faces(mesh)) + { + const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh); + if(is_almost_degenerate(tr, degenerate_threshold)) + { + ++num_almost_degenerate_tri; + continue; + } + + FT circumsphere_radius = std::sqrt(CGAL::squared_radius(tr[0], tr[1], tr[2])); + + FT a = std::sqrt(CGAL::squared_distance(tr[0], tr[1])); + FT b = std::sqrt(CGAL::squared_distance(tr[1], tr[2])); + FT c = std::sqrt(CGAL::squared_distance(tr[2], tr[0])); + FT s = 0.5 * (a + b + c); + FT inscribed_radius = std::sqrt((s * (s - a) * (s - b) * (s - c)) / s); + FT radius_ratio = circumsphere_radius / inscribed_radius; + radius_ratio /= 2.; // normalized + mean_radius_ratio += radius_ratio; + } + + mean_radius_ratio /= static_cast(mesh.number_of_faces() - num_almost_degenerate_tri); + return mean_radius_ratio; +} + +double mean_edge_ratio(const Mesh& mesh, + double degenerate_threshold) +{ + double mean_edge_ratio = 0.; + size_t num_almost_degenerate_tri = 0; + + for(const face_descriptor f : faces(mesh)) + { + const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh); + if(is_almost_degenerate(tr, degenerate_threshold)) + { + ++num_almost_degenerate_tri; + continue; + } + + FT a = std::sqrt(CGAL::squared_distance(tr[0], tr[1])); + FT b = std::sqrt(CGAL::squared_distance(tr[1], tr[2])); + FT c = std::sqrt(CGAL::squared_distance(tr[2], tr[0])); + FT min_edge = std::min({a, b, c}); + FT max_edge = std::max({a, b, c}); + FT edge_ratio = max_edge / min_edge; + + mean_edge_ratio += edge_ratio; + } + + mean_edge_ratio /= static_cast(mesh.number_of_faces() - num_almost_degenerate_tri); + return mean_edge_ratio; +} + +double mean_aspect_ratio(const Mesh& mesh, + double degenerate_threshold) +{ + double mean_aspect_ratio = 0.; + size_t num_almost_degenerate_tri = 0; + for(const face_descriptor f : faces(mesh)) + { + const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh); + if(is_almost_degenerate(tr, degenerate_threshold)) + { + ++num_almost_degenerate_tri; + continue; + } + + FT a = std::sqrt(CGAL::squared_distance(tr[0], tr[1])); + FT b = std::sqrt(CGAL::squared_distance(tr[1], tr[2])); + FT c = std::sqrt(CGAL::squared_distance(tr[2], tr[0])); + FT s = 0.5 * (a + b + c); + FT inscribed_radius = std::sqrt((s * (s - a) * (s - b) * (s - c)) / s); + FT max_edge = std::max({a, b, c}); + FT aspect_ratio = max_edge / inscribed_radius; + aspect_ratio /= (2. * std::sqrt(3.)); // normalized + mean_aspect_ratio += aspect_ratio; + } + + mean_aspect_ratio /= static_cast(mesh.number_of_faces() - num_almost_degenerate_tri); + return mean_aspect_ratio; +} + +size_t num_almost_degenerate_tri(const Mesh& mesh, + double degenerate_threshold) +{ + size_t num_almost_degenerate_tri = 0; + for(const face_descriptor f : faces(mesh)) + { + const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh); + if(is_almost_degenerate(tr, degenerate_threshold)) + { + ++num_almost_degenerate_tri; + } + } + return num_almost_degenerate_tri; +} + +int main(int argc, char** argv) +{ + const int argc_check = argc - 1; + char *entry_name_ptr = nullptr; + double relative_alpha_ratio = 20.; + double relative_offset_ratio = 600.; + + for(int i=1; i +#include + +#include +#include + +#include + +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Point_3 = Kernel::Point_3; +using Mesh = CGAL::Surface_mesh; + +namespace CGAL { +namespace Alpha_wraps_3 { +namespace internal { +namespace { + +enum Robustness_benchmark_exit_code +{ + // Success + VALID_SOLID_OUTPUT = 0, + + // Failure + INPUT_IS_INVALID = 1, + OUTPUT_IS_NOT_TRIANGLE_MESH = 2, + OUTPUT_IS_COMBINATORIAL_NON_MANIFOLD = 3, + OUTPUT_HAS_BORDERS = 4, + OUTPUT_HAS_DEGENERATED_FACES = 5, + OUTPUT_HAS_GEOMETRIC_SELF_INTERSECTIONS = 6, + OUTPUT_DOES_NOT_BOUND_VOLUME = 7, + OUTPUT_DOES_NOT_CONTAIN_INPUT = 8, + OUTPUT_DISTANCE_IS_TOO_LARGE = 9, +}; + +} // namespace +} // namespace internal +} // namespace Alpha_wraps_3 +} // namespace CGAL + +namespace PMP = CGAL::Polygon_mesh_processing; +namespace AW3i = CGAL::Alpha_wraps_3::internal; + +int main(int argc, char** argv) +{ + const int argc_check = argc - 1; + char* entry_name_ptr = nullptr; + double relative_alpha_ratio = 20.; + double relative_offset_ratio = 600.; + + for(int i=1; i $2/Robustness/results/$5/$filename.log + + python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/compute_performance_benchmark_data.py \ + -e $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/build-release/performance_benchmark -i $6 -a $3 \ + > $2/Performance/results/$5/$filename.log + + python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/compute_quality_benchmark_data.py \ + -e $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/build-release/quality_benchmark -i $6 -a $3 \ + > $2/Quality/results/$5/$filename.log +} +export -f compute_benchmark_data + +# $1: directory containing the alpha wrap project +# $2: directory containing the input data folder +# $3: directory containing the output results +# $4: alpha value +# $5: timeout value for robustness benchmark in seconds +# $6: number of virtual thread used +# $7: hash of the latest commit +# $8: hash of a commit to perform the diff with latest +cd $1 + +mkdir -p $3/Robustness/results/$7 +mkdir -p $3/Performance/results/$7 +mkdir -p $3/Quality/results/$7 +mkdir -p $3/Robustness/charts_data +mkdir -p $3/Performance/charts_data +mkdir -p $3/Quality/charts_data +mkdir -p $3/Robustness/charts +mkdir -p $3/Performance/charts +mkdir -p $3/Quality/charts +mkdir -p $3/Robustness/log +mkdir -p $3/Performance/log +mkdir -p $3/Quality/log +mkdir -p $3/charts + +find $2 -mindepth 1 | parallel -j$6 compute_benchmark_data $1 $3 $4 $5 $7 ::: + +python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Robustness/generate_robustness_benchmark_charts.py -i $3/Robustness/results/$7 -o $3/Robustness -a $4 -c $7 + +if [ -z "$8" ]; then + python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/generate_performance_benchmark_charts.py -i $3/Performance/results/$7 -o $3/Performance -a $4 -c $7; +else + python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/generate_performance_benchmark_charts.py -i $3/Performance/results/$7 -o $3/Performance -a $4 -c $7 -p $3/Performance/results/$8 -d $8; +fi + +if [ -z "$8" ]; then + python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/generate_quality_benchmark_charts.py -i $3/Quality/results/$7 -o $3/Quality -a $4 -c $7; +else + python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/generate_quality_benchmark_charts.py -i $3/Quality/results/$7 -o $3/Quality -a $4 -c $7 -p $3/Quality/results/$8 -d $8; +fi + +charts_path="$(ls "$3/Robustness/charts"/* -dArt | tail -n 1) $(ls "$3/Performance/charts"/* -dArt | tail -n 2) $(ls "$3/Quality/charts"/* -dArt | tail -n 9)" + +pdfjam --nup 2x6 $charts_path --outfile $3/charts/results_$7_$8_alpha_$4_$(date '+%Y-%m-%d_%H:%M:%S').pdf diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/CMakeLists.txt b/Alpha_wrap_3/examples/Alpha_wrap_3/CMakeLists.txt index 0be54f87eed..40187ca194c 100644 --- a/Alpha_wrap_3/examples/Alpha_wrap_3/CMakeLists.txt +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/CMakeLists.txt @@ -13,3 +13,5 @@ create_single_source_cgal_program("point_set_wrap.cpp") create_single_source_cgal_program("wrap_from_cavity.cpp") create_single_source_cgal_program("mixed_inputs_wrap.cpp") create_single_source_cgal_program("volumetric_wrap.cpp") +create_single_source_cgal_program("successive_wraps.cpp") +create_single_source_cgal_program("pause_and_resume_wrapping.cpp") diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/mixed_inputs_wrap.cpp b/Alpha_wrap_3/examples/Alpha_wrap_3/mixed_inputs_wrap.cpp index ffcedc7f1cb..90b488e0407 100644 --- a/Alpha_wrap_3/examples/Alpha_wrap_3/mixed_inputs_wrap.cpp +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/mixed_inputs_wrap.cpp @@ -103,10 +103,10 @@ int main(int argc, char** argv) oracle.add_segment_soup(segments, CGAL::parameters::default_values()); oracle.add_point_set(ps_points, CGAL::parameters::default_values()); - CGAL::Alpha_wraps_3::internal::Alpha_wrap_3 aw3(oracle); + CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3 aw3(oracle); - Mesh output_mesh; - aw3(alpha, offset, output_mesh); + Mesh wrap; + aw3(alpha, offset, wrap); t.stop(); std::cout << "Took " << t.time() << std::endl; @@ -120,10 +120,11 @@ int main(int argc, char** argv) std::string ps_name = std::string(ps_filename); ps_name = ps_name.substr(ps_name.find_last_of("/") + 1, ps_name.length() - 1); ps_name = ps_name.substr(0, ps_name.find_last_of(".")); - std::string output_name = ts_name + "_" + ss_name + "_" + ps_name + "_" + std::to_string(static_cast(relative_alpha)) - + "_" + std::to_string(static_cast(relative_offset)) + ".off"; + std::string output_name = ts_name + "_" + ss_name + "_" + ps_name + "_" + + std::to_string(static_cast(relative_alpha)) + "_" + + std::to_string(static_cast(relative_offset)) + ".off"; std::cout << "Writing to " << output_name << std::endl; - CGAL::IO::write_polygon_mesh(output_name, output_mesh, CGAL::parameters::stream_precision(17)); + CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17)); return EXIT_SUCCESS; } diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/output_helper.h b/Alpha_wrap_3/examples/Alpha_wrap_3/output_helper.h new file mode 100644 index 00000000000..581ac82bff0 --- /dev/null +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/output_helper.h @@ -0,0 +1,19 @@ +#ifndef CGAL_ALPHA_WRAP_3_EXAMPLES_OUTPUT_HELPER_H +#define CGAL_ALPHA_WRAP_3_EXAMPLES_OUTPUT_HELPER_H + +#include + +std::string generate_output_name(std::string input_name, + const double alpha, + const double offset) +{ + input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1); + input_name = input_name.substr(0, input_name.find_last_of(".")); + std::string output_name = input_name + + "_" + std::to_string(static_cast(alpha)) + + "_" + std::to_string(static_cast(offset)) + ".off"; + + return output_name; +} + +#endif // CGAL_ALPHA_WRAP_3_EXAMPLES_OUTPUT_HELPER_H diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/pause_and_resume_wrapping.cpp b/Alpha_wrap_3/examples/Alpha_wrap_3/pause_and_resume_wrapping.cpp new file mode 100644 index 00000000000..cee169048c3 --- /dev/null +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/pause_and_resume_wrapping.cpp @@ -0,0 +1,164 @@ +// This example demonstrates how to interrupt the wrapping process before it has terminated, +// and how to resume afterwards. +// +// -------------------------------- !! Warning !! -------------------------------------------------- +// By default, the wrapper uses an unsorted LIFO queue of faces to refine. This means that +// the intermediate result is not very useful because the algorithm carves deep and not wide +// (somewhat like a DFS vs a BFS). +// +// The sorted queue option is enabled with the macro below to make the refinement algorithm +// more uniform. The downside is that it is slower. +// ------------------------------------------------------------------------------------------------- +#define CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + +#include "output_helper.h" + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +namespace AW3 = CGAL::Alpha_wraps_3; +namespace PMP = CGAL::Polygon_mesh_processing; + +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using Point_3 = K::Point_3; + +using Points = std::vector; +using Face = std::array; +using Faces = std::vector; + +using Mesh = CGAL::Surface_mesh; +using face_descriptor = boost::graph_traits::face_descriptor; + +struct Interrupter_visitor + : public AW3::internal::Wrapping_default_visitor +{ + using Base = AW3::internal::Wrapping_default_visitor; + + CGAL::Real_timer timer; + double max_time = -1; // in seconds + +public: + void set_max_time(double t) { max_time = t; } + +public: + template + void on_flood_fill_begin(const AlphaWrapper&) + { + std::cout << "Starting timer..." << std::endl; + timer.start(); + } + + template + bool go_further(const Wrapper&) + { + if(timer.time() > max_time) + { + timer.stop(); + std::cout << "Paused after " << timer.time() << " s." << std::endl; + return false; + } + + return true; + } +}; + +int main(int argc, char** argv) +{ + std::cout.precision(17); + std::cerr.precision(17); + + CGAL::Random rng; + std::cout << "Random seed = " << rng.get_seed() << std::endl; + + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off"); + + // = read the soup + Points points; + Faces faces; + if(!CGAL::IO::read_polygon_soup(filename, points, faces) || faces.empty()) + { + std::cerr << "Invalid soup input: " << filename << std::endl; + return EXIT_FAILURE; + } + + std::cout << "Input: " << points.size() << " points, " << faces.size() << " faces" << std::endl; + + // Compute the alpha and offset values + const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : rng.get_double(150., 200.); + const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.; + std::cout << "relative_alpha = " << relative_alpha << std::endl; + + CGAL::Bbox_3 bbox; + for(const Point_3& p : points) + bbox += p.bbox(); + + const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) + + CGAL::square(bbox.ymax() - bbox.ymin()) + + CGAL::square(bbox.zmax() - bbox.zmin())); + + const double alpha = diag_length / relative_alpha; + const double offset = diag_length / relative_offset; + + // Build the wrapper + using Oracle = CGAL::Alpha_wraps_3::internal::Triangle_soup_oracle; + Oracle oracle(alpha); + oracle.add_triangle_soup(points, faces, CGAL::parameters::default_values()); + CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3 aw3(oracle); + + // --- Launch the wrapping, and pause when the algorithm has spent 1s flooding + Interrupter_visitor interrupter; + interrupter.set_max_time(1.); + + Mesh wrap; + aw3(alpha, offset, wrap, CGAL::parameters::visitor(interrupter)); + std::cout << ">>> The current wrap has " << num_vertices(wrap) << " vertices" << std::endl; + CGAL::IO::write_polygon_mesh("stopped_1.off", wrap, CGAL::parameters::stream_precision(17)); + + // --- Restart from the previous state, and pause a bit further + interrupter.set_max_time(2.); + aw3(alpha, offset, wrap, CGAL::parameters::visitor(interrupter) + .refine_triangulation(true)); + std::cout << ">>> The current wrap has " << num_vertices(wrap) << " vertices" << std::endl; + CGAL::IO::write_polygon_mesh("stopped_2.off", wrap, CGAL::parameters::stream_precision(17)); + + // --- Restart from the previous state, and let it finish + aw3(alpha, offset, wrap, CGAL::parameters::refine_triangulation(true)); + std::cout << ">>> The final (resumed) wrap has " << num_vertices(wrap) << " vertices" << std::endl; + std::string output_name = generate_output_name(filename, relative_alpha, relative_offset); + std::cout << "Writing to " << "resumed_" + output_name << std::endl; + CGAL::IO::write_polygon_mesh("resumed_" + output_name, wrap, CGAL::parameters::stream_precision(17)); + + // --- Get the final wrap, in one go: + Mesh single_pass_wrap; + CGAL::alpha_wrap_3(points, faces, alpha, offset, single_pass_wrap); + std::cout << ">>> The final (from scratch) wrap has " << num_vertices(single_pass_wrap) << " vertices" << std::endl; + + output_name = generate_output_name(filename, relative_alpha, relative_offset); + std::cout << "Writing to " << output_name << std::endl; + CGAL::IO::write_polygon_mesh(output_name, single_pass_wrap, CGAL::parameters::stream_precision(17)); + + // --- Compare the results to ensure both approaches yield identical meshes + std::vector > common; + std::vector m1_only; + std::vector m2_only; + PMP::match_faces(wrap, single_pass_wrap, + std::back_inserter(common), + std::back_inserter(m1_only), + std::back_inserter(m2_only)); + if(!m1_only.empty() || !m2_only.empty()) + { + std::cerr << "Error: The two wraps should have been identical!" << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/point_set_wrap.cpp b/Alpha_wrap_3/examples/Alpha_wrap_3/point_set_wrap.cpp index 8742a2a2001..a602cf5c58b 100644 --- a/Alpha_wrap_3/examples/Alpha_wrap_3/point_set_wrap.cpp +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/point_set_wrap.cpp @@ -1,3 +1,5 @@ +#include "output_helper.h" + #include #include @@ -23,7 +25,7 @@ int main(int argc, char** argv) Point_container points; if(!CGAL::IO::read_points(filename, std::back_inserter(points)) || points.empty()) { - std::cerr << "Invalid input." << std::endl; + std::cerr << "Invalid input:" << filename << std::endl; return EXIT_FAILURE; } @@ -53,11 +55,7 @@ int main(int argc, char** argv) std::cout << "Took " << t.time() << " s." << std::endl; // Save the result - std::string input_name = std::string(filename); - input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1); - input_name = input_name.substr(0, input_name.find_last_of(".")); - std::string output_name = input_name + "_" + std::to_string(static_cast(relative_alpha)) - + "_" + std::to_string(static_cast(relative_offset)) + ".off"; + const std::string output_name = generate_output_name(filename, relative_alpha, relative_offset); std::cout << "Writing to " << output_name << std::endl; CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17)); diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/successive_wraps.cpp b/Alpha_wrap_3/examples/Alpha_wrap_3/successive_wraps.cpp new file mode 100644 index 00000000000..301aec46e2d --- /dev/null +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/successive_wraps.cpp @@ -0,0 +1,135 @@ +// In this example, we reuse the underlying triangulation of the previous state, and carve using +// a new (smaller) alpha value. This enables considerable speed-up: the cumulated time taken +// to run `n` successive instances of `{alpha_wrap(alpha_i)}_(i=1...n)` will be roughly equal +// to the time taken to the single instance of alpha_wrap(alpha_n) from scratch. +// +// The speed-up increases with the number of intermediate results, and on the gap between +// alpha values: if alpha_2 is close to alpha_1, practically no new computation are required, +// and the speed-up is almost 100%. +// +// -------------------------------- !! Warning !! -------------------------------------------------- +// The result of: +// > alpha_wrap(alpha_1, ...) +// > alpha_wrap(alpha_2, ..., reuse) +// is not exactly identical to calling directly: +// > alpha_wrap(alpha_2, ..., do_not_reuse) +// because the queues are sorted slightly differently and the AABB tree is rebuilt differently +// to optimize the runtime. +// ------------------------------------------------------------------------------------------------- + +#include "output_helper.h" + +#include +#include + +#include +#include +#include +#include + +#include +#include + +namespace PMP = CGAL::Polygon_mesh_processing; + +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using FT = K::FT; +using Point_3 = K::Point_3; + +using Mesh = CGAL::Surface_mesh; + +int main(int argc, char** argv) +{ + std::cout.precision(17); + std::cerr.precision(17); + + // Read the input + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cube.off"); + std::cout << "Reading " << filename << "..." << std::endl; + + Mesh mesh; + if(!PMP::IO::read_polygon_mesh(filename, mesh) || is_empty(mesh) || !is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input:" << filename << std::endl; + return EXIT_FAILURE; + } + + std::cout << "Input: " << num_vertices(mesh) << " vertices, " << num_faces(mesh) << " faces" << std::endl; + + const CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(mesh); + const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) + + CGAL::square(bbox.ymax() - bbox.ymin()) + + CGAL::square(bbox.zmax() - bbox.zmin())); + + // We want decreasing alphas, and these are relative ratios, so they need to be increasing + const std::vector relative_alphas = { 1, 50, 100, 150, 200, 250 }; + const FT relative_offset = 600; + + // =============================================================================================== + // Naive approach: + + CGAL::Real_timer t; + double total_time = 0.; + + for(std::size_t i=0; i>> [" << i << "] alpha: " << alpha << " offset: " << offset << std::endl; + + Mesh wrap; + CGAL::alpha_wrap_3(mesh, alpha, offset, wrap); + + t.stop(); + std::cout << " Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl; + std::cout << " Elapsed time: " << t.time() << " s." << std::endl; + + total_time += t.time(); + } + + std::cout << "Total elapsed time (naive): " << total_time << " s.\n" << std::endl; + + // =============================================================================================== + // Re-use approach + + total_time = 0.; + t.reset(); + + using Oracle = CGAL::Alpha_wraps_3::internal::Triangle_mesh_oracle; + using Wrapper = CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3; + Wrapper wrapper; // contains the triangulation that is being refined iteratively + + for(std::size_t i=0; i>> [" << i << "] alpha: " << alpha << " offset: " << offset << std::endl; + + // The triangle mesh oracle should be initialized with alpha to internally perform a split + // of too-big facets while building the AABB Tree. This split in fact yields a significant + // speed-up for meshes with elements that are large compared to alpha. This speed-up makes it + // faster to re-build the AABB tree for every value of alpha than to use a non-optimized tree. + Oracle oracle(alpha); + oracle.add_triangle_mesh(mesh, CGAL::parameters::default_values()); + wrapper.oracle() = oracle; + + Mesh wrap; + wrapper(alpha, offset, wrap, CGAL::parameters::refine_triangulation((i != 0))); + + t.stop(); + std::cout << " Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl; + std::cout << " Elapsed time: " << t.time() << " s." << std::endl; + + total_time += t.time(); + } + + std::cout << "Total elapsed time (successive): " << total_time << " s." << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/triangle_mesh_wrap.cpp b/Alpha_wrap_3/examples/Alpha_wrap_3/triangle_mesh_wrap.cpp index 00e2e4fd9fc..d4953490446 100644 --- a/Alpha_wrap_3/examples/Alpha_wrap_3/triangle_mesh_wrap.cpp +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/triangle_mesh_wrap.cpp @@ -1,3 +1,5 @@ +#include "output_helper.h" + #include #include @@ -25,7 +27,7 @@ int main(int argc, char** argv) Mesh mesh; if(!PMP::IO::read_polygon_mesh(filename, mesh) || is_empty(mesh) || !is_triangle_mesh(mesh)) { - std::cerr << "Invalid input." << std::endl; + std::cerr << "Invalid input:" << filename << std::endl; return EXIT_FAILURE; } @@ -56,12 +58,7 @@ int main(int argc, char** argv) std::cout << "Took " << t.time() << " s." << std::endl; // Save the result - std::string input_name = std::string(filename); - input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1); - input_name = input_name.substr(0, input_name.find_last_of(".")); - std::string output_name = input_name - + "_" + std::to_string(static_cast(relative_alpha)) - + "_" + std::to_string(static_cast(relative_offset)) + ".off"; + const std::string output_name = generate_output_name(filename, relative_alpha, relative_offset); std::cout << "Writing to " << output_name << std::endl; CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17)); diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/triangle_soup_wrap.cpp b/Alpha_wrap_3/examples/Alpha_wrap_3/triangle_soup_wrap.cpp index 626e3bdc3ba..51e04974c28 100644 --- a/Alpha_wrap_3/examples/Alpha_wrap_3/triangle_soup_wrap.cpp +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/triangle_soup_wrap.cpp @@ -1,3 +1,5 @@ +#include "output_helper.h" + #include #include @@ -30,7 +32,7 @@ int main(int argc, char** argv) std::vector > faces; if(!CGAL::IO::read_polygon_soup(filename, points, faces) || faces.empty()) { - std::cerr << "Invalid input." << std::endl; + std::cerr << "Invalid input:" << filename << std::endl; return EXIT_FAILURE; } @@ -63,12 +65,7 @@ int main(int argc, char** argv) std::cout << "Took " << t.time() << " s." << std::endl; // Save the result - std::string input_name = std::string(filename); - input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1); - input_name = input_name.substr(0, input_name.find_last_of(".")); - std::string output_name = input_name - + "_" + std::to_string(static_cast(relative_alpha)) - + "_" + std::to_string(static_cast(relative_offset)) + ".off"; + const std::string output_name = generate_output_name(filename, relative_alpha, relative_offset); std::cout << "Writing to " << output_name << std::endl; CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17)); diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/volumetric_wrap.cpp b/Alpha_wrap_3/examples/Alpha_wrap_3/volumetric_wrap.cpp index 113215b631a..ddc6f765612 100644 --- a/Alpha_wrap_3/examples/Alpha_wrap_3/volumetric_wrap.cpp +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/volumetric_wrap.cpp @@ -1,3 +1,5 @@ +#include "output_helper.h" + #include #include @@ -70,7 +72,7 @@ int main(int argc, char** argv) Faces faces; if(!CGAL::IO::read_polygon_soup(filename, points, faces) || faces.empty()) { - std::cerr << "Invalid input." << std::endl; + std::cerr << "Invalid input:" << filename << std::endl; return EXIT_FAILURE; } @@ -101,7 +103,7 @@ int main(int argc, char** argv) Oracle oracle(K{}); oracle.add_triangle_soup(points, faces, CGAL::parameters::default_values()); - CGAL::Alpha_wraps_3::internal::Alpha_wrap_3 aw3(oracle); + CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3 aw3(oracle); Mesh wrap; aw3(alpha, offset, wrap); @@ -113,12 +115,7 @@ int main(int argc, char** argv) auto dt = aw3.triangulation(); // Save the result - std::string input_name = std::string(filename); - input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1); - input_name = input_name.substr(0, input_name.find_last_of(".")); - std::string output_name = input_name - + "_" + std::to_string(static_cast(relative_alpha)) - + "_" + std::to_string(static_cast(relative_offset)) + ".off"; + const std::string output_name = generate_output_name(filename, relative_alpha, relative_offset); std::cout << "Writing to " << output_name << std::endl; CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17)); diff --git a/Alpha_wrap_3/examples/Alpha_wrap_3/wrap_from_cavity.cpp b/Alpha_wrap_3/examples/Alpha_wrap_3/wrap_from_cavity.cpp index 1c1d6c7b3d7..970bb583484 100644 --- a/Alpha_wrap_3/examples/Alpha_wrap_3/wrap_from_cavity.cpp +++ b/Alpha_wrap_3/examples/Alpha_wrap_3/wrap_from_cavity.cpp @@ -1,3 +1,5 @@ +#include "output_helper.h" + #include #include @@ -25,13 +27,13 @@ int main(int argc, char** argv) if(!PMP::IO::read_polygon_mesh(filename, input) || is_empty(input) || !is_triangle_mesh(input)) { - std::cerr << "Invalid input." << std::endl; + std::cerr << "Invalid input:" << filename << std::endl; return EXIT_FAILURE; } std::cout << "Input: " << num_vertices(input) << " vertices, " << num_faces(input) << " faces" << std::endl; - const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 30.; + const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 40.; const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.; // Compute the alpha and offset values diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h index 04b30cf7678..d6ca4fb356c 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h @@ -49,7 +49,9 @@ #include #include #include -#include +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + #include +#endif #include #include #include @@ -92,6 +94,10 @@ struct Wrapping_default_visitor template void on_flood_fill_begin(const AlphaWrapper&) { } + // Whether the flood filling process should continue + template + constexpr bool go_further(const Wrapper&) { return true; } + template void before_facet_treatment(const AlphaWrapper&, const Gate&) { } @@ -110,7 +116,7 @@ struct Wrapping_default_visitor template -class Alpha_wrap_3 +class Alpha_wrapper_3 { using Oracle = Oracle_; @@ -120,22 +126,36 @@ class Alpha_wrap_3 using Default_Vb = Alpha_wrap_triangulation_vertex_base_3; using Default_Cb = Alpha_wrap_triangulation_cell_base_3; - using Default_Cbt = Cell_base_with_timestamp; // determinism + using Default_Cbt = Cell_base_with_timestamp; // for determinism using Default_Tds = CGAL::Triangulation_data_structure_3; using Default_Triangulation = CGAL::Delaunay_triangulation_3; +public: using Triangulation = typename Default::Get::type; + // Use the geom traits from the triangulation, and trust the (advanced) user that provided it + using Geom_traits = typename Triangulation::Geom_traits; + +private: using Cell_handle = typename Triangulation::Cell_handle; using Facet = typename Triangulation::Facet; using Vertex_handle = typename Triangulation::Vertex_handle; using Locate_type = typename Triangulation::Locate_type; using Gate = internal::Gate; - using Alpha_PQ = Modifiable_priority_queue, CGAL_BOOST_PAIRING_HEAP>; - // Use the geom traits from the triangulation, and trust the (advanced) user that provided it - using Geom_traits = typename Triangulation::Geom_traits; + // A sorted queue is a priority queue sorted by circumradius, and is experimentally significantly + // slower. However, intermediate results are usable: at each point of the algorithm, the wrap + // has a somewhat uniform mesh element size. + // + // An unsorted queue is a LIFO queue, and is experimentally much faster (~35%), + // but intermediate results are not useful: a LIFO carves deep before than wide, + // and can thus for example leave scaffolding faces till almost the end of the refinement. +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + using Alpha_PQ = Modifiable_priority_queue, CGAL_BOOST_PAIRING_HEAP>; +#else + using Alpha_PQ = std::stack; +#endif using FT = typename Geom_traits::FT; using Point_3 = typename Geom_traits::Point_3; @@ -149,33 +169,48 @@ class Alpha_wrap_3 using SC_Iso_cuboid_3 = SC::Iso_cuboid_3; using SC2GT = Cartesian_converter; + using Seeds = std::vector; protected: - const Oracle m_oracle; + Oracle m_oracle; SC_Iso_cuboid_3 m_bbox; - FT m_alpha, m_sq_alpha; - FT m_offset, m_sq_offset; + FT m_alpha = FT(-1), m_sq_alpha = FT(-1); + FT m_offset = FT(-1), m_sq_offset = FT(-1); + + Seeds m_seeds; Triangulation m_tr; + Alpha_PQ m_queue; public: - // Main constructor - Alpha_wrap_3(const Oracle& oracle) - : - m_oracle(oracle), - m_tr(Geom_traits(oracle.geom_traits())), - // used to set up the initial MPQ, use some arbitrary not-too-small value - m_queue(4096) + Alpha_wrapper_3() +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + // '4096' is an arbitrary, not-too-small value for the largest ID in queue initialization + : m_queue(4096) +#endif { // Due to the Steiner point computation being a dichotomy, the algorithm is inherently inexact // and passing exact kernels is explicitly disabled to ensure no misunderstanding. static_assert(std::is_floating_point::value); } + Alpha_wrapper_3(const Oracle& oracle) + : + m_oracle(oracle), + m_tr(Geom_traits(oracle.geom_traits())) +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + , m_queue(4096) +#endif + { + static_assert(std::is_floating_point::value); + } + public: const Geom_traits& geom_traits() const { return m_tr.geom_traits(); } + Oracle& oracle() { return m_oracle; } + const Oracle& oracle() const { return m_oracle; } Triangulation& triangulation() { return m_tr; } const Triangulation& triangulation() const { return m_tr; } const Alpha_PQ& queue() const { return m_queue; } @@ -222,26 +257,53 @@ public: using parameters::get_parameter_reference; using parameters::choose_parameter; + // using OVPM = typename CGAL::GetVertexPointMap::type; OVPM ovpm = choose_parameter(get_parameter(out_np, internal_np::vertex_point), get_property_map(vertex_point, output_mesh)); - typedef typename internal_np::Lookup_named_param_def < - internal_np::visitor_t, - InputNamedParameters, - Wrapping_default_visitor // default - >::reference Visitor; + // + using Visitor = typename internal_np::Lookup_named_param_def< + internal_np::visitor_t, + InputNamedParameters, + Wrapping_default_visitor // default + >::reference; Wrapping_default_visitor default_visitor; Visitor visitor = choose_parameter(get_parameter_reference(in_np, internal_np::visitor), default_visitor); - std::vector no_seeds; - using Seeds = typename internal_np::Lookup_named_param_def< - internal_np::seed_points_t, InputNamedParameters, std::vector >::reference; - Seeds seeds = choose_parameter(get_parameter_reference(in_np, internal_np::seed_points), no_seeds); + // Points used to create initial cavities + m_seeds = choose_parameter(get_parameter_reference(in_np, internal_np::seed_points), Seeds()); + // Whether or not some cells should be reflagged as "inside" after the refinement+carving loop + // as ended, as to ensure that the result is not only combinatorially manifold, but also + // geometrically manifold. + // + // -- Warning -- + // Manifoldness postprocessing will be performed even if the wrapping is interrupted (and + // this option is enabled). const bool do_enforce_manifoldness = choose_parameter(get_parameter(in_np, internal_np::do_enforce_manifoldness), true); + // Whether to keep pockets of "outside" cells that are not connected to the exterior (or to the + // initial cavities, if used). + // + // -- Warning -- + // Pockets of "outside" cells will be purged even if the wrapping is interrupted (and + // this option is enabled). + const bool keep_inner_ccs = choose_parameter(get_parameter(in_np, internal_np::keep_inner_connected_components), false); + + // This parameter enables avoiding recomputing the triangulation from scratch when wrapping + // the same input for multiple values of alpha (and typically the same offset values). + // + // -- Warning -- + // If this is enabled, the 3D triangulation will NOT be re-initialized at launch. + // This means that the triangulation is NOT cleared, even if: + // - you use an alpha value that is greater than what was used in a previous run; you will + // obtain the same result as the last run. + // - you use a different offset value between runs, you might then get points that are not + // on the offset surface corresponding to that corresponding to the latter offset value. + const bool refining = choose_parameter(get_parameter(in_np, internal_np::refine_triangulation), false); + #ifdef CGAL_AW3_TIMER CGAL::Real_timer t; t.start(); @@ -249,76 +311,59 @@ public: visitor.on_alpha_wrapping_begin(*this); - if(!initialize(alpha, offset, seeds)) + if(!initialize(alpha, offset, refining)) return; -#ifdef CGAL_AW3_DEBUG_DUMP_EVERY_STEP - extract_surface(output_mesh, ovpm, true /*tolerate non manifoldness*/); - CGAL::IO::write_polygon_mesh("initial_cavities.off", output_mesh, - CGAL::parameters::vertex_point_map(ovpm).stream_precision(17)); +#ifdef CGAL_AW3_TIMER + t.stop(); + std::cout << "Initialization took: " << t.time() << " s." << std::endl; + t.reset(); + t.start(); +#endif + +#ifdef CGAL_AW3_DEBUG_DUMP_INTERMEDIATE_WRAPS + dump_triangulation_faces("starting_wrap.off", true /*only_boundary_faces*/); #endif alpha_flood_fill(visitor); +#ifdef CGAL_AW3_DEBUG_DUMP_INTERMEDIATE_WRAPS + dump_triangulation_faces("flood_filled_wrap.off", true /*only_boundary_faces*/); +#endif + #ifdef CGAL_AW3_TIMER t.stop(); std::cout << "Flood filling took: " << t.time() << " s." << std::endl; + t.reset(); + t.start(); #endif if(do_enforce_manifoldness) { -#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS - std::cout << "> Make manifold..." << std::endl; - - extract_surface(output_mesh, ovpm, true /*tolerate non manifoldness*/); - -#ifdef CGAL_AW3_DEBUG_DUMP_EVERY_STEP - dump_triangulation_faces("intermediate_dt3.off", false /*only_boundary_faces*/); - IO::write_polygon_mesh("intermediate.off", output_mesh, - CGAL::parameters::vertex_point_map(ovpm).stream_precision(17)); -#endif - - FT base_vol = 0; - if(is_closed(output_mesh)) // might not be due to manifoldness - base_vol = PMP::volume(output_mesh, CGAL::parameters::vertex_point_map(ovpm)); - else - std::cerr << "Warning: couldn't compute volume before manifoldness fixes (mesh is not closed)" << std::endl; -#endif - -#ifdef CGAL_AW3_TIMER - t.reset(); - t.start(); -#endif - make_manifold(); +#ifdef CGAL_AW3_DEBUG_DUMP_INTERMEDIATE_WRAPS + dump_triangulation_faces("manifold_wrap.off", true /*only_boundary_faces*/); +#endif + #ifdef CGAL_AW3_TIMER t.stop(); - std::cout << "\nManifoldness post-processing took: " << t.time() << " s." << std::endl; + std::cout << "Manifoldness post-processing took: " << t.time() << " s." << std::endl; + t.reset(); + t.start(); #endif + } -#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS - if(!is_zero(base_vol)) - { - extract_surface(output_mesh, ovpm, false /*do not tolerate non-manifoldness*/); + if(!keep_inner_ccs) + { + // We could purge *before* manifold enforcement, but making the mesh manifold is + // very cheap in most cases, so it is better to keep the code simpler. + purge_inner_connected_components(); - const FT manifold_vol = PMP::volume(output_mesh, CGAL::parameters::vertex_point_map(ovpm)); - const FT ratio = manifold_vol / base_vol; - - std::cout << "Volumes post-manifoldness fix:\n" - << "before: " << base_vol << "\n" - << "after: " << manifold_vol << "\n" - << "ratio: " << ratio << std::endl; - if(ratio > 1.1) // more than 10% extra volume - std::cerr << "Warning: large increase of volume after manifoldness resolution" << std::endl; - } -#endif - } // do_enforce_manifoldness - -#ifdef CGAL_AW3_TIMER - t.reset(); - t.start(); +#ifdef CGAL_AW3_DEBUG_DUMP_INTERMEDIATE_WRAPS + dump_triangulation_faces("purged_wrap.off", true /*only_boundary_faces*/); #endif + } extract_surface(output_mesh, ovpm, !do_enforce_manifoldness); @@ -331,9 +376,9 @@ public: std::cout << "Alpha wrap vertices: " << vertices(output_mesh).size() << std::endl; std::cout << "Alpha wrap faces: " << faces(output_mesh).size() << std::endl; - #ifdef CGAL_AW3_DEBUG_DUMP_EVERY_STEP - IO::write_polygon_mesh("final.off", output_mesh, CGAL::parameters::stream_precision(17)); - dump_triangulation_faces("final_dt3.off", false /*only_boundary_faces*/); + #ifdef CGAL_AW3_DEBUG_DUMP_INTERMEDIATE_WRAPS + IO::write_polygon_mesh("final_wrap.off", output_mesh, CGAL::parameters::stream_precision(17)); + dump_triangulation_faces("final_tr.off", false /*only_boundary_faces*/); #endif #endif @@ -382,6 +427,20 @@ public: return bbox; } +private: + // The distinction between inside boundary and outside boundary is the presence of cells + // being flagged for manifoldness: inside boundary considers those outside, and outside + // boundary considers them inside. + bool is_on_inside_boundary(Cell_handle ch, Cell_handle nh) const + { + return (ch->is_inside() != nh->is_inside()); // one is "inside", the other is not + } + + bool is_on_outside_boundary(Cell_handle ch, Cell_handle nh) const + { + return (ch->is_outside() != nh->is_outside()); // one is "outside", the other is not + } + private: // Adjust the bbox & insert its corners to construct the starting triangulation void insert_bbox_corners() @@ -409,20 +468,20 @@ private: // - Cells whose circumcenter is in the offset volume are inside: this is because // we need to have outside cell circumcenters outside of the volume to ensure // that the refinement point is separated from the existing point set. - bool cavity_cell_outside_tag(const Cell_handle ch) + Cell_label cavity_cell_label(const Cell_handle ch) { CGAL_precondition(!m_tr.is_infinite(ch)); const Tetrahedron_with_outside_info tet(ch, geom_traits()); if(m_oracle.do_intersect(tet)) - return false; + return Cell_label::INSIDE; const Point_3& ch_cc = circumcenter(ch); typename Geom_traits::Construct_ball_3 ball = geom_traits().construct_ball_3_object(); const Ball_3 ch_cc_offset_ball = ball(ch_cc, m_sq_offset); const bool is_cc_in_offset = m_oracle.do_intersect(ch_cc_offset_ball); - return !is_cc_in_offset; + return is_cc_in_offset ? Cell_label::INSIDE : Cell_label::OUTSIDE; } // Create a cavity using seeds rather than starting from the infinity. @@ -461,15 +520,14 @@ private: // // Another way is to simply make faces incident to the seed always traversable, and then // we only have to ensure faces opposite of the seed are traversable (i.e., radius ~= 1.65 * alpha) - template - bool initialize_with_cavities(const SeedRange& seeds) + bool initialize_with_cavities() { #ifdef CGAL_AW3_DEBUG_INITIALIZATION std::cout << "> Dig cavities" << std::endl; - std::cout << seeds.size() << " seed(s)" << std::endl; + std::cout << m_seeds.size() << " seed(s)" << std::endl; #endif - CGAL_precondition(!seeds.empty()); + CGAL_precondition(!m_seeds.empty()); // Get a double value approximating the scaling factors // std::cout << sqrt(3) * sin(2pi / 5) << std::endl; @@ -478,7 +536,7 @@ private: Iso_cuboid_3 bbox = SC2GT()(m_bbox); std::vector seed_vs; - for(const Point_3& seed_p : seeds) + for(const Point_3& seed_p : m_seeds) { #ifdef CGAL_AW3_DEBUG_INITIALIZATION std::cout << "Initialize from seed " << seed_p << std::endl; @@ -504,7 +562,7 @@ private: continue; } - // Mark the seeds and icosahedron vertices as "artificial vertices" such that the facets + // Mark the seeds and icosahedron vertices as "scaffolding" vertices such that the facets // incident to these vertices are always traversable regardless of their circumcenter. // This is done because otherwise some cavities can appear on the mesh: non-traversable facets // with two vertices on the offset, and the third being a deeper inside seed / ico_seed. @@ -573,17 +631,19 @@ private: std::vector inc_cells; inc_cells.reserve(64); m_tr.incident_cells(seed_v, std::back_inserter(inc_cells)); + for(Cell_handle ch : inc_cells) - ch->is_outside() = cavity_cell_outside_tag(ch); + ch->set_label(cavity_cell_label(ch)); } - // Might as well go through the full triangulation since only seeds should have been inserted + // Should be cheap enough to go through the full triangulation as only seeds have been inserted for(Cell_handle ch : m_tr.all_cell_handles()) { - if(!ch->is_outside()) + if(ch->is_inside()) continue; - // When the algorithm starts from a manually dug hole, infinite cells are tagged "inside" + // When the algorithm starts from a manually dug hole, infinite cells are initialized + // as "inside" such that they do not appear on the boundary CGAL_assertion(!m_tr.is_infinite(ch)); for(int i=0; i<4; ++i) @@ -601,7 +661,7 @@ private: return true; } - // tag all infinite cells OUTSIDE and all finite cells INSIDE + // tag all infinite cells "outside" and all finite cells "inside" // init queue with all convex hull facets bool initialize_from_infinity() { @@ -609,20 +669,52 @@ private: { if(m_tr.is_infinite(ch)) { - ch->is_outside() = true; + ch->set_label(Cell_label::OUTSIDE); const int inf_index = ch->index(m_tr.infinite_vertex()); push_facet(std::make_pair(ch, inf_index)); } else { - ch->is_outside() = false; + ch->set_label(Cell_label::INSIDE); } } return true; } -public: + void reset_manifold_labels() + { + // No erase counter increment, or it will mess up with a possibly non-empty queue. + for(Cell_handle ch : m_tr.all_cell_handles()) + if(ch->label() == Cell_label::MANIFOLD) + ch->set_label(Cell_label::OUTSIDE); + } + + // This function is used in the case of resumption of a previous run: m_tr is not cleared, + // and we fill the queue with the new parameters. + bool initialize_from_existing_triangulation() + { +#ifdef CGAL_AW3_DEBUG_INITIALIZATION + std::cout << "Restart from a DT of " << m_tr.number_of_cells() << " cells" << std::endl; +#endif + + for(Cell_handle ch : m_tr.all_cell_handles()) + { + if(ch->is_inside()) + continue; + + for(int i=0; i<4; ++i) + { + if(ch->neighbor(i)->is_inside()) + push_facet(std::make_pair(ch, i)); + + } + } + + return true; + } + +private: // Manifoldness is tolerated while debugging and extracting at intermediate states // Not the preferred way because it uses 3*nv storage template @@ -662,7 +754,6 @@ public: if(cell->tds_data().processed()) continue; - cell->tds_data().mark_processed(); for(int fid=0; fid<4; ++fid) @@ -689,8 +780,6 @@ public: } } - PMP::duplicate_non_manifold_edges_in_polygon_soup(points, faces); - CGAL_assertion(PMP::is_polygon_soup_a_polygon_mesh(faces)); PMP::polygon_soup_to_polygon_mesh(points, faces, output_mesh, CGAL::parameters::default_values(), @@ -708,6 +797,8 @@ public: CGAL_postcondition(is_closed(output_mesh)); PMP::orient_to_bound_a_volume(output_mesh, CGAL::parameters::vertex_point_map(ovpm)); + + collect_garbage(output_mesh); } template @@ -717,7 +808,7 @@ public: namespace PMP = Polygon_mesh_processing; #ifdef CGAL_AW3_DEBUG - std::cout << "> Extract wrap... ()" << std::endl; + std::cout << "> Extract manifold wrap... ()" << std::endl; #endif CGAL_assertion_code(for(Vertex_handle v : m_tr.finite_vertex_handles())) @@ -725,35 +816,30 @@ public: clear(output_mesh); - // boundary faces to polygon soup std::vector points; std::vector > faces; std::unordered_map vertex_to_id; - std::size_t nv = 0; - for(auto fit=m_tr.finite_facets_begin(), fend=m_tr.finite_facets_end(); fit!=fend; ++fit) + for(auto fit=m_tr.all_facets_begin(), fend=m_tr.all_facets_end(); fit!=fend; ++fit) { Facet f = *fit; if(!f.first->is_outside()) f = m_tr.mirror_facet(f); - const Cell_handle c = f.first; + const Cell_handle ch = f.first; const int s = f.second; - const Cell_handle nh = c->neighbor(s); - if(c->is_outside() == nh->is_outside()) + const Cell_handle nh = ch->neighbor(s); + if(!is_on_outside_boundary(ch, nh)) continue; std::array ids; for(int pos=0; pos<3; ++pos) { - Vertex_handle vh = c->vertex(Triangulation::vertex_triple_index(s, pos)); - auto insertion_res = vertex_to_id.emplace(vh, nv); + Vertex_handle vh = ch->vertex(Triangulation::vertex_triple_index(s, pos)); + auto insertion_res = vertex_to_id.emplace(vh, vertex_to_id.size()); if(insertion_res.second) // successful insertion, never-seen-before vertex - { points.push_back(m_tr.point(vh)); - ++nv; - } ids[pos] = insertion_res.first->second; } @@ -761,12 +847,22 @@ public: faces.emplace_back(std::array{ids[0], ids[1], ids[2]}); } +#ifdef CGAL_AW3_DEBUG + std::cout << "\t" << points.size() << " points" << std::endl; + std::cout << "\t" << faces.size() << " faces" << std::endl; +#endif + if(faces.empty()) + { +#ifdef CGAL_AW3_DEBUG + std::cerr << "Warning: empty wrap?..." << std::endl; +#endif return; + } if(!PMP::is_polygon_soup_a_polygon_mesh(faces)) { - CGAL_warning_msg(false, "Could NOT extract mesh..."); + CGAL_warning_msg(false, "Failed to extract a manifold boundary!"); return; } @@ -777,10 +873,12 @@ public: CGAL_postcondition(!is_empty(output_mesh)); CGAL_postcondition(is_valid_polygon_mesh(output_mesh)); CGAL_postcondition(is_closed(output_mesh)); + CGAL_postcondition(PMP::does_bound_a_volume(output_mesh, CGAL::parameters::vertex_point_map(ovpm))); PMP::orient_to_bound_a_volume(output_mesh, CGAL::parameters::vertex_point_map(ovpm)); } +public: template void extract_surface(OutputMesh& output_mesh, OVPM ovpm, @@ -886,43 +984,76 @@ private: } private: - enum Facet_queue_status + // A permissive gate is a gate that we can traverse without checking its circumradius + enum class Facet_status { IRRELEVANT = 0, - ARTIFICIAL_FACET, + HAS_INFINITE_NEIGHBOR, // the cell incident to the mirrored facet is infinite (permissive) + SCAFFOLDING, // incident to a SEED or BBOX vertex (permissive) TRAVERSABLE }; + inline const char* get_status_message(const Facet_status status) + { + constexpr std::size_t status_count = 4; + + // Messages corresponding to Error_code list above. Must be kept in sync! + static const char* message[status_count] = { + "Irrelevant facet", + "Facet incident to infinite neighbor", + "Facet with a bbox/seed vertex", + "Traversable facet" + }; + + const std::size_t status_id = static_cast(status); + if(status_id > status_count || status_id < 0) + return "Unknown status"; + else + return message[status_id]; + } + +public: // @speed some decent time may be spent constructing Facet (pairs) for no reason as it's always // just to grab the .first and .second as soon as it's constructed, and not due to API requirements - // e.g. from DT3 - Facet_queue_status facet_status(const Facet& f) const + Facet_status facet_status(const Facet& f) const { CGAL_precondition(!m_tr.is_infinite(f)); #ifdef CGAL_AW3_DEBUG_FACET_STATUS std::cout << "facet status: " - << f.first->vertex((f.second + 1)&3)->point() << " " - << f.first->vertex((f.second + 2)&3)->point() << " " - << f.first->vertex((f.second + 3)&3)->point() << std::endl; + << m_tr.point(f.first, Triangulation::vertex_triple_index(f.second, 0)) << " " + << m_tr.point(f.first, Triangulation::vertex_triple_index(f.second, 1)) << " " + << m_tr.point(f.first, Triangulation::vertex_triple_index(f.second, 2)) << std::endl; #endif - // skip if neighbor is OUTSIDE or infinite + // skip if neighbor is "outside" or infinite const Cell_handle ch = f.first; const int id = f.second; + CGAL_precondition(ch->label() == Cell_label::INSIDE || ch->label() == Cell_label::OUTSIDE); + + if(!ch->is_outside()) // "inside" or "manifold" + { +#ifdef CGAL_AW3_DEBUG_FACET_STATUS + std::cout << "Facet is inside" << std::endl; +#endif + return Facet_status::IRRELEVANT; + } + const Cell_handle nh = ch->neighbor(id); + CGAL_precondition(ch->label() == Cell_label::INSIDE || ch->label() == Cell_label::OUTSIDE); + if(m_tr.is_infinite(nh)) - return TRAVERSABLE; + return Facet_status::HAS_INFINITE_NEIGHBOR; if(nh->is_outside()) { #ifdef CGAL_AW3_DEBUG_FACET_STATUS std::cout << "Neighbor already outside" << std::endl; #endif - return IRRELEVANT; + return Facet_status::IRRELEVANT; } - // push if facet is connected to artificial vertices + // push if facet is connected to scaffolding vertices for(int i=0; i<3; ++i) { const Vertex_handle vh = ch->vertex(Triangulation::vertex_triple_index(id, i)); @@ -930,9 +1061,9 @@ private: vh->type() == AW3i::Vertex_type:: SEED_VERTEX) { #ifdef CGAL_AW3_DEBUG_FACET_STATUS - std::cout << "artificial facet due to artificial vertex #" << i << std::endl; + std::cout << "Scaffolding facet due to vertex #" << i << std::endl; #endif - return ARTIFICIAL_FACET; + return Facet_status::SCAFFOLDING; } } @@ -942,78 +1073,143 @@ private: #ifdef CGAL_AW3_DEBUG_FACET_STATUS std::cout << "traversable" << std::endl; #endif - return TRAVERSABLE; + return Facet_status::TRAVERSABLE; } #ifdef CGAL_AW3_DEBUG_FACET_STATUS std::cout << "not traversable" << std::endl; #endif - return IRRELEVANT; + return Facet_status::IRRELEVANT; } +private: bool push_facet(const Facet& f) { CGAL_precondition(f.first->is_outside()); +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE // skip if f is already in queue if(m_queue.contains_with_bounds_check(Gate(f))) return false; +#endif - const Facet_queue_status s = facet_status(f); - if(s == IRRELEVANT) + // @todo could avoid useless facet_status() calls by doing it after the zombie check + // for the unsorted priority queue, but AFAIR, it doesn't save noticeable time (and that + // increases the queue size). + const Facet_status status = facet_status(f); + if(status == Facet_status::IRRELEVANT) return false; - const Cell_handle ch = f.first; - const int id = f.second; - const Point_3& p0 = m_tr.point(ch, (id+1)&3); - const Point_3& p1 = m_tr.point(ch, (id+2)&3); - const Point_3& p2 = m_tr.point(ch, (id+3)&3); +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + const FT sqr = smallest_squared_radius_3(f, m_tr); + const bool is_permissive = (status == Facet_status::HAS_INFINITE_NEIGHBOR || + status == Facet_status::SCAFFOLDING); + m_queue.resize_and_push(Gate(f, sqr, is_permissive)); +#else + m_queue.push(Gate(f, m_tr)); +#endif - // @todo should prob be the real value we compare to alpha instead of squared_radius - const FT sqr = geom_traits().compute_squared_radius_3_object()(p0, p1, p2); - m_queue.resize_and_push(Gate(f, sqr, (s == ARTIFICIAL_FACET))); +#ifdef CGAL_AW3_DEBUG_QUEUE + const Cell_handle ch = f.first; + const int s = f.second; + const Point_3& p0 = m_tr.point(ch, Triangulation::vertex_triple_index(s, 0)); + const Point_3& p1 = m_tr.point(ch, Triangulation::vertex_triple_index(s, 1)); + const Point_3& p2 = m_tr.point(ch, Triangulation::vertex_triple_index(s, 2)); + + static int gid = 0; + std::cout << "Queue insertion #" << gid++ << "\n" + << " ch = " << &*ch << " (" << m_tr.is_infinite(ch) << ") " << "\n" + << "\t" << p0 << "\n\t" << p1 << "\n\t" << p2 << std::endl; + std::cout << " Status: " << get_status_message(status) << std::endl; + #ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + std::cout << " SQR: " << sqr << std::endl; + std::cout << " Permissiveness: " << is_permissive << std::endl; + + CGAL_assertion(is_permissive || sqr >= m_sq_alpha); + #endif // CGAL_AW3_USE_SORTED_PRIORITY_QUEUE +#endif // CGAL_AW3_DEBUG_QUEUE return true; } private: - template bool initialize(const double alpha, const double offset, - const SeedRange& seeds) + const bool refining) { #ifdef CGAL_AW3_DEBUG std::cout << "> Initialize..." << std::endl; - std::cout << "Alpha: " << alpha << std::endl; - std::cout << "Offset: " << offset << std::endl; +#endif + + const bool resuming = refining && (alpha == m_alpha) && (offset == m_offset); + +#ifdef CGAL_AW3_DEBUG + std::cout << "\tAlpha: " << alpha << std::endl; + std::cout << "\tOffset: " << offset << std::endl; + std::cout << "\tRefining? " << refining << std::endl; + std::cout << "\tResuming? " << resuming << std::endl; #endif if(!is_positive(alpha) || !is_positive(offset)) { #ifdef CGAL_AW3_DEBUG - std::cout << "Error: invalid input parameters" << std::endl; + std::cerr << "Error: invalid input parameters: " << alpha << " and" << offset << std::endl; #endif return false; } + if(refining && alpha > m_alpha) + std::cerr << "Warning: refining with an alpha greater than the last iteration's!" << std::endl; + if(refining && offset != m_offset) + std::cerr << "Warning: refining with a different offset value!" << std::endl; + m_alpha = FT(alpha); m_sq_alpha = square(m_alpha); m_offset = FT(offset); m_sq_offset = square(m_offset); - m_tr.clear(); + // Resuming means that we do not need to re-initialize the queue: either we have finished + // and there is nothing to do, or the interruption was due to a user callback in the visitor, + // and we can resume with the current queue + if(resuming) + { +#ifdef CGAL_AW3_DEBUG + std::cout << "Resuming with a queue of size: " << m_queue.size() << std::endl; +#endif + + reset_manifold_labels(); + return true; + } + +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE m_queue.clear(); +#else + m_queue = { }; +#endif - insert_bbox_corners(); + if(refining) + { + // If we are re-using the triangulation, change the label of the extra elements + // that we have added to ensure a manifold result back to external ("manifold" -> "outside") + reset_manifold_labels(); - if(seeds.empty()) - return initialize_from_infinity(); + return initialize_from_existing_triangulation(); + } else - return initialize_with_cavities(seeds); + { + m_tr.clear(); + + insert_bbox_corners(); + + if(m_seeds.empty()) + return initialize_from_infinity(); + else + return initialize_with_cavities(); + } } template - void alpha_flood_fill(Visitor& visitor) + bool alpha_flood_fill(Visitor& visitor) { #ifdef CGAL_AW3_DEBUG std::cout << "> Flood fill..." << std::endl; @@ -1024,33 +1220,47 @@ private: // Explore all finite cells that are reachable from one of the initial outside cells. while(!m_queue.empty()) { + if(!visitor.go_further(*this)) + return false; + #ifdef CGAL_AW3_DEBUG_QUEUE_PP check_queue_sanity(); #endif // const& to something that will be popped, but safe as `ch` && `id` are extracted before the pop const Gate& gate = m_queue.top(); + +#ifndef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + if(gate.is_zombie()) + { + m_queue.pop(); + continue; + } +#endif + const Facet& f = gate.facet(); CGAL_precondition(!m_tr.is_infinite(f)); const Cell_handle ch = f.first; - const int id = f.second; - const Cell_handle neighbor = ch->neighbor(id); + const int s = f.second; + CGAL_precondition(ch->is_outside()); + + const Cell_handle nh = ch->neighbor(s); + CGAL_precondition(nh->label() == Cell_label::INSIDE || nh->label() == Cell_label::OUTSIDE); #ifdef CGAL_AW3_DEBUG_QUEUE static int fid = 0; std::cout << m_tr.number_of_vertices() << " DT vertices" << std::endl; std::cout << m_queue.size() << " facets in the queue" << std::endl; std::cout << "Face " << fid++ << "\n" - << "c = " << &*ch << " (" << m_tr.is_infinite(ch) << "), n = " << &*neighbor << " (" << m_tr.is_infinite(neighbor) << ")" << "\n" - << m_tr.point(ch, (id+1)&3) << "\n" << m_tr.point(ch, (id+2)&3) << "\n" << m_tr.point(ch, (id+3)&3) << std::endl; - std::cout << "Priority: " << gate.priority() << std::endl; + << "c = " << &*ch << " (" << m_tr.is_infinite(ch) << "), n = " << &*nh << " (" << m_tr.is_infinite(nh) << ")" << "\n" + << m_tr.point(ch, Triangulation::vertex_triple_index(s, 0)) << "\n" + << m_tr.point(ch, Triangulation::vertex_triple_index(s, 1)) << "\n" + << m_tr.point(ch, Triangulation::vertex_triple_index(s, 2)) << std::endl; + std::cout << "Priority: " << gate.priority() << " (sq alpha: " << m_sq_alpha << ")" << std::endl; + std::cout << "Permissiveness: " << gate.is_permissive_facet() << std::endl; #endif - visitor.before_facet_treatment(*this, gate); - - m_queue.pop(); - #ifdef CGAL_AW3_DEBUG_DUMP_EVERY_STEP static int i = 0; std::string step_name = "results/steps/step_" + std::to_string(static_cast(i)) + ".off"; @@ -1059,18 +1269,27 @@ private: std::string face_name = "results/steps/face_" + std::to_string(static_cast(i++)) + ".xyz"; std::ofstream face_out(face_name); face_out.precision(17); - face_out << "3\n" << m_tr.point(ch, (id+1)&3) << "\n" << m_tr.point(ch, (id+2)&3) << "\n" << m_tr.point(ch, (id+3)&3) << std::endl; + face_out << "3\n" << m_tr.point(ch, Triangulation::vertex_triple_index(s, 0)) << "\n" + << m_tr.point(ch, Triangulation::vertex_triple_index(s, 1)) << "\n" + << m_tr.point(ch, Triangulation::vertex_triple_index(s, 2)) << std::endl; face_out.close(); #endif - if(m_tr.is_infinite(neighbor)) + visitor.before_facet_treatment(*this, gate); + + m_queue.pop(); + + if(m_tr.is_infinite(nh)) { - neighbor->is_outside() = true; + nh->set_label(Cell_label::OUTSIDE); +#ifndef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + nh->increment_erase_counter(); +#endif continue; } Point_3 steiner_point; - if(compute_steiner_point(ch, neighbor, steiner_point)) + if(compute_steiner_point(ch, nh, steiner_point)) { // std::cout << CGAL::abs(CGAL::approximate_sqrt(m_oracle.squared_distance(steiner_point)) - m_offset) // << " vs " << 1e-2 * m_offset << std::endl; @@ -1079,7 +1298,7 @@ private: // locate cells that are going to be destroyed and remove their facet from the queue int li, lj = 0; Locate_type lt; - const Cell_handle conflict_cell = m_tr.locate(steiner_point, lt, li, lj, neighbor); + const Cell_handle conflict_cell = m_tr.locate(steiner_point, lt, li, lj, nh); CGAL_assertion(lt != Triangulation::VERTEX); // Using small vectors like in Triangulation_3 does not bring any runtime improvement @@ -1092,6 +1311,7 @@ private: std::back_inserter(boundary_facets), std::back_inserter(conflict_zone)); +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE // Purge the queue of facets that will be deleted/modified by the Steiner point insertion, // and which might have been gates for(const Cell_handle& cch : conflict_zone) @@ -1110,6 +1330,7 @@ private: if(m_queue.contains_with_bounds_check(Gate(mf))) m_queue.erase(Gate(mf)); } +#endif visitor.before_Steiner_point_insertion(*this, steiner_point); @@ -1124,10 +1345,10 @@ private: std::vector new_cells; new_cells.reserve(32); m_tr.incident_cells(vh, std::back_inserter(new_cells)); - for(const Cell_handle& ch : new_cells) + for(const Cell_handle& new_ch : new_cells) { - // std::cout << "new cell has time stamp " << ch->time_stamp() << std::endl; - ch->is_outside() = m_tr.is_infinite(ch); + // std::cout << "new cell has time stamp " << new_ch->time_stamp() << std::endl; + new_ch->set_label(m_tr.is_infinite(new_ch) ? Cell_label::OUTSIDE : Cell_label::INSIDE); } // Push all new boundary facets to the queue. @@ -1135,34 +1356,37 @@ private: // because we need to handle internal facets, infinite facets, and also more subtle changes // such as a new cell being marked inside which now creates a boundary // with its incident "outside" flagged cell. - for(Cell_handle ch : new_cells) + for(Cell_handle new_ch : new_cells) { for(int i=0; i<4; ++i) { - if(m_tr.is_infinite(ch, i)) + if(m_tr.is_infinite(new_ch, i)) continue; - const Cell_handle nh = ch->neighbor(i); - if(nh->is_outside() == ch->is_outside()) // not on the boundary + const Cell_handle new_nh = new_ch->neighbor(i); + if(new_nh->label() == new_ch->label()) // not on a boundary continue; - const Facet boundary_f = std::make_pair(ch, i); - if(ch->is_outside()) + const Facet boundary_f = std::make_pair(new_ch, i); + if(new_ch->is_outside()) push_facet(boundary_f); else push_facet(m_tr.mirror_facet(boundary_f)); } } } - else + else // no need for a Steiner point, carve through and continue { - // tag neighbor as OUTSIDE - neighbor->is_outside() = true; + nh->set_label(Cell_label::OUTSIDE); +#ifndef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + nh->increment_erase_counter(); +#endif // for each finite facet of neighbor, push it to the queue - for(int i=0; i<4; ++i) + const int mi = m_tr.mirror_index(ch, s); + for(int i=1; i<4; ++i) { - const Facet neighbor_f = std::make_pair(neighbor, i); + const Facet neighbor_f = std::make_pair(nh, (mi+i)&3); push_facet(neighbor_f); } } @@ -1172,11 +1396,103 @@ private: // Check that no useful facet has been ignored CGAL_postcondition_code(for(auto fit=m_tr.finite_facets_begin(), fend=m_tr.finite_facets_end(); fit!=fend; ++fit) {) - CGAL_postcondition_code( if(fit->first->is_outside() == fit->first->neighbor(fit->second)->is_outside()) continue;) + CGAL_postcondition_code( Cell_handle ch = fit->first; Cell_handle nh = fit->first->neighbor(fit->second); ) + CGAL_postcondition_code( if(ch->label() == nh->label()) continue;) CGAL_postcondition_code( Facet f = *fit;) - CGAL_postcondition_code( if(!fit->first->is_outside()) f = m_tr.mirror_facet(f);) - CGAL_postcondition( facet_status(f) == IRRELEVANT); + CGAL_postcondition_code( if(ch->is_inside()) f = m_tr.mirror_facet(f);) + CGAL_postcondition( facet_status(f) == Facet_status::IRRELEVANT); CGAL_postcondition_code(}) + + return true; + } + + // Any outside cell that isn't reachable from infinity is a cavity that can be discarded. + std::size_t purge_inner_connected_components() + { +#ifdef CGAL_AW3_DEBUG + std::cout << "> Purge inner islands..." << std::endl; + + std::size_t pre_counter = 0; + for(Cell_handle ch : m_tr.all_cell_handles()) + if(ch->is_outside()) + ++pre_counter; + std::cout << pre_counter << " / " << m_tr.all_cell_handles().size() << " (pre purge)" << std::endl; +#endif + + std::size_t label_change_counter = 0; + + std::stack cells_to_visit; + + if(!m_seeds.empty()) + { + for(const Point_3& seed : m_seeds) + { + Locate_type lt; + int li, lj; + Cell_handle ch = m_tr.locate(seed, lt, li, lj); + + if(!ch->is_outside()) + { + std::cerr << "Warning: cell containing seed is not outside?!" << std::endl; + continue; + } + + cells_to_visit.push(ch); + } + } + else // typical flooding from outside + { + cells_to_visit.push(m_tr.infinite_vertex()->cell()); + } + + while(!cells_to_visit.empty()) + { + Cell_handle curr_c = cells_to_visit.top(); + cells_to_visit.pop(); + + CGAL_assertion(curr_c->is_outside()); + + if(curr_c->tds_data().processed()) + continue; + curr_c->tds_data().mark_processed(); + + for(int j=0; j<4; ++j) + { + Cell_handle neigh_c = curr_c->neighbor(j); + if(neigh_c->tds_data().processed() || !neigh_c->is_outside()) + continue; + + cells_to_visit.push(neigh_c); + } + } + + for(Cell_handle ch : m_tr.all_cell_handles()) + { + if(ch->tds_data().is_clear() && ch->is_outside()) + { + ch->set_label(Cell_label::INSIDE); +#ifndef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + ch->increment_erase_counter(); +#endif + ++label_change_counter; + } + } + + // reset the conflict flags + for(Cell_handle ch : m_tr.all_cell_handles()) + ch->tds_data().clear(); + +#ifdef CGAL_AW3_DEBUG + std::size_t post_counter = 0; + for(Cell_handle ch : m_tr.all_cell_handles()) + if(ch->is_outside()) + ++post_counter; + std::cout << post_counter << " / " << m_tr.all_cell_handles().size() << " (pre purge)" << std::endl; + + std::cout << label_change_counter << " label changes" << std::endl; +#endif + + return label_change_counter; } private: @@ -1190,7 +1506,7 @@ private: inc_cells.reserve(64); m_tr.incident_cells(v, std::back_inserter(inc_cells)); - // Flood one inside and outside CC. + // Flood one inside and outside CC within the cell umbrella of the vertex. // Process both an inside and an outside CC to also detect edge pinching. // If there are still unprocessed afterwards, there is a non-manifoldness issue. // @@ -1204,7 +1520,7 @@ private: if(ic->is_outside()) outside_start = ic; else if(inside_start == Cell_handle()) - inside_start = ic; + inside_start = ic; // can be "inside" or "manifold" } // fully inside / outside @@ -1236,8 +1552,10 @@ private: CGAL_assertion(neigh_c->has_vertex(v)); if(neigh_c->tds_data().processed() || - neigh_c->is_outside() != curr_c->is_outside()) // do not cross the boundary + is_on_outside_boundary(neigh_c, curr_c)) // do not cross the boundary + { continue; + } cells_to_visit.push(neigh_c); } @@ -1330,22 +1648,47 @@ public: // Not the best complexity, but it's very cheap compared to the rest of the algorithm. void make_manifold() { - namespace PMP = Polygon_mesh_processing; +#ifdef CGAL_AW3_DEBUG + std::cout << "> Make manifold..." << std::endl; - // This seems more harmful than useful after the priority queue has been introduced since - // it adds a lot of flat cells into the triangulation, which then get added to the mesh - // during manifoldness fixing. + auto wrap_volume = [&]() + { + FT vol = 0; + for(Cell_handle ch : m_tr.finite_cell_handles()) + if(!ch->is_outside()) + vol += volume(m_tr.point(ch, 0), m_tr.point(ch, 1), m_tr.point(ch, 2), m_tr.point(ch, 3)); + + return vol; + }; + + #ifdef CGAL_AW3_DEBUG_DUMP_INTERMEDIATE_WRAPS + dump_triangulation_faces("carved_tr.off", true /*only_boundary_faces*/); + #endif + + FT base_vol = wrap_volume(); + if(!is_positive(base_vol)) + std::cerr << "Warning: wrap with non-positive volume?" << std::endl; +#endif + + // This ends up more harmful than useful after the priority queue has been introduced since + // it usually results in a lot of flat cells into the triangulation, which then get added + // to the mesh during manifoldness fixing. // remove_bbox_vertices(); std::stack non_manifold_vertices; // @todo sort somehow? for(Vertex_handle v : m_tr.finite_vertex_handles()) { if(is_non_manifold(v)) + { +#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS_PP + std::cout << v->point() << " is non-manifold" << std::endl; +#endif non_manifold_vertices.push(v); + } } // Some lambdas for the comparer - auto has_artificial_vertex = [](Cell_handle c) -> bool + auto has_scaffolding_vertex = [](Cell_handle c) -> bool { for(int i=0; i<4; ++i) { @@ -1359,79 +1702,77 @@ public: return false; }; - auto is_on_boundary = [](Cell_handle c, int i) -> bool - { - return (c->is_outside() != c->neighbor(i)->is_outside()); - }; + // This originally seemed like a good idea, but in the end it can have strong cascading issues, + // whereas some cells with much smaller volume could have solved the non-manifoldness. +// auto is_on_boundary = [](Cell_handle c, int i) -> bool +// { +// return is_on_outside_boundary(c, c->neighbor(i)); +// }; +// +// auto count_boundary_facets = [&](Cell_handle c, Vertex_handle v) -> int +// { +// const int v_index_in_c = c->index(v); +// int boundary_facets = 0; +// for(int i=0; i<3; ++i) // also take into account the opposite facet? +// { +// if(i == v_index_in_c) +// continue; +// +// boundary_facets += is_on_boundary(c, i); +// } +// +// return boundary_facets; +// }; - auto count_boundary_facets = [&](Cell_handle c, Vertex_handle v) -> int - { - const int v_index_in_c = c->index(v); - int boundary_facets = 0; - for(int i=0; i<3; ++i) // also take into account the opposite facet? - { - if(i == v_index_in_c) - continue; - - boundary_facets += is_on_boundary(c, i); - } - - return boundary_facets; - }; - - // longest edge works better + // Experimentally, longest edge works better // auto sq_circumradius = [&](Cell_handle c) -> FT // { // const Point_3& cc = circumcenter(c); // return geom_traits().compute_squared_distance_3_object()(m_tr.point(c, 0), cc); // }; + // The reasoning behind using longest edge rather than volume is that we want to avoid + // spikes (which would have a small volume), and can often happen since we do not spend + // any care on the quality of tetrahedra. auto sq_longest_edge = [&](Cell_handle c) -> FT { return (std::max)({ squared_distance(m_tr.point(c, 0), m_tr.point(c, 1)), squared_distance(m_tr.point(c, 0), m_tr.point(c, 2)), squared_distance(m_tr.point(c, 0), m_tr.point(c, 3)), squared_distance(m_tr.point(c, 1), m_tr.point(c, 2)), - squared_distance(m_tr.point(c, 3), m_tr.point(c, 3)), + squared_distance(m_tr.point(c, 1), m_tr.point(c, 3)), squared_distance(m_tr.point(c, 2), m_tr.point(c, 3)) }); }; -#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS +#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS_PP std::cout << non_manifold_vertices.size() << " initial NMV" << std::endl; #endif while(!non_manifold_vertices.empty()) { -#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS +#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS_PP std::cout << non_manifold_vertices.size() << " NMV in queue" << std::endl; #endif Vertex_handle v = non_manifold_vertices.top(); non_manifold_vertices.pop(); -#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS - std::cout << "·"; -#endif - if(!is_non_manifold(v)) continue; // Prioritize: // - cells without bbox vertices - // - cells that already have a large number of boundary facets // - small cells when equal number of boundary facets - // @todo give topmost priority to cells with > 1 non-manifold vertex? + // + // Note that these are properties that do not depend on the cell labels, and so we only need + // to sort once. However, if a criterion such as the number of incident inside cells were added, + // we would need to sort after each modification of "inside"/"outside" labels. auto comparer = [&](Cell_handle l, Cell_handle r) -> bool { - if(has_artificial_vertex(l)) - return false; - if(has_artificial_vertex(r)) - return true; + CGAL_precondition(!m_tr.is_infinite(l) && !m_tr.is_infinite(r)); - const int l_bf_count = count_boundary_facets(l, v); - const int r_bf_count = count_boundary_facets(r, v); - if(l_bf_count != r_bf_count) - return l_bf_count > r_bf_count; + if(has_scaffolding_vertex(l) != has_scaffolding_vertex(r)) + return has_scaffolding_vertex(r); return sq_longest_edge(l) < sq_longest_edge(r); }; @@ -1440,22 +1781,22 @@ public: inc_cells.reserve(64); m_tr.finite_incident_cells(v, std::back_inserter(inc_cells)); -#define CGAL_AW3_USE_BRUTE_FORCE_MUTABLE_PRIORITY_QUEUE -#ifndef CGAL_AW3_USE_BRUTE_FORCE_MUTABLE_PRIORITY_QUEUE - std::sort(inc_cells.begin(), inc_cells.end(), comparer); // sort once -#endif + std::vector finite_outside_inc_cells; + finite_outside_inc_cells.reserve(64); + std::copy_if(inc_cells.begin(), inc_cells.end(), std::back_inserter(finite_outside_inc_cells), + [&](Cell_handle c) -> bool { return !m_tr.is_infinite(c) && c->is_outside(); }); - for(auto cit=inc_cells.begin(), cend=inc_cells.end(); cit!=cend; ++cit) + // 'std::stable_sort' to have determinism without having to write something like: + // if(longest_edge(l) == longest_edge(r)) return ... + // in the comparer. It's almost always a small range, so the extra cost does not matter. + std::stable_sort(finite_outside_inc_cells.begin(), finite_outside_inc_cells.end(), comparer); + + for(Cell_handle ic : finite_outside_inc_cells) { -#ifdef CGAL_AW3_USE_BRUTE_FORCE_MUTABLE_PRIORITY_QUEUE - // sort at every iteration since the number of boundary facets evolves - std::sort(cit, cend, comparer); -#endif - Cell_handle ic = *cit; - CGAL_assertion(!m_tr.is_infinite(ic)); + CGAL_assertion(!m_tr.is_infinite(ic) && ic->is_outside()); // This is where new material is added - ic->is_outside() = false; + ic->set_label(Cell_label::MANIFOLD); #ifdef CGAL_AW3_DEBUG_DUMP_EVERY_STEP static int i = 0; @@ -1470,6 +1811,8 @@ public: CGAL_assertion(!is_non_manifold(v)); + // Check if the new material has not created a non-manifold configuration. + // @speed this could be done on only the vertices of cells whose labels have changed. std::vector adj_vertices; adj_vertices.reserve(64); m_tr.finite_adjacent_vertices(v, std::back_inserter(adj_vertices)); @@ -1481,12 +1824,34 @@ public: CGAL_assertion_code(for(Vertex_handle v : m_tr.finite_vertex_handles())) CGAL_assertion(!is_non_manifold(v)); + +#ifdef CGAL_AW3_DEBUG + std::size_t nm_cells_counter = 0; + for(Cell_handle ch : m_tr.all_cell_handles()) + if(ch->label() == Cell_label::MANIFOLD) + ++nm_cells_counter; + std::cout << "Number of added cells: " << nm_cells_counter << std::endl; + + if(!is_zero(base_vol)) + { + const FT manifold_vol = wrap_volume(); + const FT ratio = manifold_vol / base_vol; + + std::cout << "Volumes post-manifoldness fix:\n" + << "before: " << base_vol << "\n" + << "after: " << manifold_vol << "\n" + << "ratio: " << ratio << std::endl; + if(ratio > 1.1) // more than 10% extra volume + std::cerr << "Warning: large increase of volume after manifoldness resolution" << std::endl; + } +#endif } private: void check_queue_sanity() { - std::cout << "Check queue sanity..." << std::endl; + std::cout << "\t~~~ Check queue sanity ~~~" << std::endl; + std::vector queue_gates; Gate previous_top_gate = m_queue.top(); while(!m_queue.empty()) @@ -1496,24 +1861,31 @@ private: const Facet& current_f = current_gate.facet(); const Cell_handle ch = current_f.first; const int id = current_f.second; - const Point_3& p0 = m_tr.point(ch, (id+1)&3); - const Point_3& p1 = m_tr.point(ch, (id+2)&3); - const Point_3& p2 = m_tr.point(ch, (id+3)&3); - const FT sqr = geom_traits().compute_squared_radius_3_object()(p0, p1, p2); + const Point_3& p0 = m_tr.point(ch, Triangulation::vertex_triple_index(id, 0)); + const Point_3& p1 = m_tr.point(ch, Triangulation::vertex_triple_index(id, 1)); + const Point_3& p2 = m_tr.point(ch, Triangulation::vertex_triple_index(id, 2)); - std::cout << "At Facet with VID " << get(Gate_ID_PM(), current_gate) << std::endl; + std::cout << "Facet with VID " << get(Gate_ID_PM(), current_gate) << "\n"; + std::cout << "\t" << p0 << "\n\t" << p1 << "\n\t" << p2 << "\n"; - if(current_gate.priority() != sqr) - std::cerr << "Error: facet in queue has wrong priority" << std::endl; +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + std::cout << " Permissiveness: " << current_gate.is_permissive_facet() << "\n"; + std::cout << " SQR: " << geom_traits().compute_squared_radius_3_object()(p0, p1, p2) << "\n"; + std::cout << " Priority " << current_gate.priority() << std::endl; if(Less_gate()(current_gate, previous_top_gate)) std::cerr << "Error: current gate has higher priority than the previous top" << std::endl; previous_top_gate = current_gate; +#else + if(current_gate.is_zombie()) + std::cout << "Gate is a zombie!" << std::endl; +#endif m_queue.pop(); } - std::cout << "End sanity check" << std::endl; + + std::cout << "\t~~~ End queue sanity check ~~~" << std::endl; // Rebuild CGAL_assertion(m_queue.empty()); @@ -1536,17 +1908,17 @@ private: for(auto fit=m_tr.finite_facets_begin(), fend=m_tr.finite_facets_end(); fit!=fend; ++fit) { - Cell_handle c = fit->first; + Cell_handle ch = fit->first; int s = fit->second; - Cell_handle nc = c->neighbor(s); - if(only_boundary_faces && (c->is_outside() == nc->is_outside())) + Cell_handle nh = ch->neighbor(s); + if(only_boundary_faces && ch->label() == nh->label()) continue; std::array ids; for(std::size_t pos=0; pos<3; ++pos) { - Vertex_handle v = c->vertex((s+pos+1)&3); + Vertex_handle v = ch->vertex((s+pos+1)&3); auto insertion_res = vertex_to_id.emplace(v, nv); if(insertion_res.second) { diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_triangulation_cell_base_3.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_triangulation_cell_base_3.h index db1df61c8f6..efaeb82d330 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_triangulation_cell_base_3.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_triangulation_cell_base_3.h @@ -20,18 +20,26 @@ namespace CGAL { namespace Alpha_wraps_3 { namespace internal { +enum class Cell_label +{ + // Cells that have been carved + OUTSIDE, + // Cells that have not yet been carved + INSIDE, + // OUTSIDE cells that have been labeled "inside" again as to make the result manifold + MANIFOLD +}; + template < typename GT, typename Cb = CGAL::Delaunay_triangulation_cell_base_with_circumcenter_3 > class Alpha_wrap_triangulation_cell_base_3 : public Cb { -private: - bool outside = false; - public: typedef typename Cb::Vertex_handle Vertex_handle; typedef typename Cb::Cell_handle Cell_handle; +public: template < typename TDS2 > struct Rebind_TDS { @@ -39,6 +47,14 @@ public: using Other = Alpha_wrap_triangulation_cell_base_3; }; +private: + Cell_label m_label = Cell_label::INSIDE; + +#ifndef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + unsigned int m_erase_counter; +#endif + +public: Alpha_wrap_triangulation_cell_base_3() : Cb() {} @@ -55,8 +71,26 @@ public: : Cb(v0, v1, v2, v3, n0, n1, n2, n3) {} - bool is_outside() const { return outside; } - bool& is_outside() { return outside; } +public: + Cell_label label() const { return m_label; } + void set_label(const Cell_label label) { m_label = label; } + bool is_inside() const { return m_label == Cell_label::INSIDE; } + bool is_outside() const { return m_label == Cell_label::OUTSIDE; } + +#ifndef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + unsigned int erase_counter() const + { + return m_erase_counter; + } + void set_erase_counter(unsigned int c) + { + m_erase_counter = c; + } + void increment_erase_counter() + { + ++m_erase_counter; + } +#endif }; template diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Oracle_base.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Oracle_base.h index 85bd11c00b6..7f3534a484c 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Oracle_base.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Oracle_base.h @@ -321,7 +321,6 @@ public: typename AABB_tree::Bounding_box bbox() const { CGAL_precondition(!empty()); - return tree().bbox(); } diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Point_set_oracle.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Point_set_oracle.h index 7bad2ff313d..33f8a52b178 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Point_set_oracle.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Point_set_oracle.h @@ -28,6 +28,7 @@ #include #include #include +#include #include namespace CGAL { @@ -41,6 +42,7 @@ struct PS_oracle_traits using Geom_traits = Alpha_wrap_AABB_geom_traits; // Wrap the kernel to add Ball_3 + custom Do_intersect_3 using Points = std::vector; + using Points_ptr = std::shared_ptr; using PR_iterator = typename Points::const_iterator; using Primitive = AABB_primitive; private: - Points m_points; + Points_ptr m_points_ptr; public: // Constructors - Point_set_oracle() - : Oracle_base(BaseOracle(), Base_GT()) - { } - Point_set_oracle(const BaseOracle& base_oracle, const Base_GT& gt = Base_GT()) : Oracle_base(base_oracle, gt) - { } + { + m_points_ptr = std::make_shared(); + } Point_set_oracle(const Base_GT& gt, const BaseOracle& base_oracle = BaseOracle()) - : Oracle_base(base_oracle, gt) + : Point_set_oracle(base_oracle, gt) + { } + + Point_set_oracle() + : Point_set_oracle(BaseOracle(), Base_GT()) { } public: @@ -101,21 +106,27 @@ public: if(points.empty()) { #ifdef CGAL_AW3_DEBUG - std::cout << "Warning: Input is empty " << std::endl; + std::cout << "Warning: Input is empty (PS)" << std::endl; #endif return; } - const std::size_t old_size = m_points.size(); - m_points.insert(std::cend(m_points), std::cbegin(points), std::cend(points)); + const std::size_t old_size = m_points_ptr->size(); + m_points_ptr->insert(std::cend(*m_points_ptr), std::cbegin(points), std::cend(points)); #ifdef CGAL_AW3_DEBUG std::cout << "Insert into AABB tree (points)..." << std::endl; #endif - this->tree().insert(std::next(std::cbegin(m_points), old_size), std::cend(m_points)); + this->tree().insert(std::next(std::cbegin(*m_points_ptr), old_size), std::cend(*m_points_ptr)); - CGAL_postcondition(this->tree().size() == m_points.size()); + // Manually constructing it here purely for profiling reasons: if we keep the lazy approach, + // it will be done at the first treatment of a facet that needs a Steiner point. + // So if one wanted to bench the flood fill runtime, it would be skewed by the time it takes + // to accelerate the tree. + this->tree().accelerate_distance_queries(); + + CGAL_postcondition(this->tree().size() == m_points_ptr->size()); } }; diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Segment_soup_oracle.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Segment_soup_oracle.h index d02a9f9faaf..349cd012ef7 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Segment_soup_oracle.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Segment_soup_oracle.h @@ -28,6 +28,7 @@ #include #include #include +#include #include namespace CGAL { @@ -40,7 +41,9 @@ struct SS_oracle_traits { using Geom_traits = Alpha_wrap_AABB_geom_traits; // Wrap the kernel to add Ball_3 + custom Do_intersect_3 - using Segments = std::vector; + using Segment = typename GT_::Segment_3; + using Segments = std::vector; + using Segments_ptr = std::shared_ptr; using SR_iterator = typename Segments::const_iterator; using Primitive = AABB_primitive; private: - Segments m_segments; + Segments_ptr m_segments_ptr; public: // Constructors - Segment_soup_oracle() - : Oracle_base(BaseOracle(), Base_GT()) - { } - Segment_soup_oracle(const BaseOracle& base_oracle, const Base_GT& gt = Base_GT()) : Oracle_base(base_oracle, gt) - { } + { + m_segments_ptr = std::make_shared(); + } Segment_soup_oracle(const Base_GT& gt, const BaseOracle& base_oracle = BaseOracle()) - : Oracle_base(base_oracle, gt) + : Segment_soup_oracle(base_oracle, gt) + { } + + Segment_soup_oracle() + : Segment_soup_oracle(BaseOracle(), Base_GT()) { } public: @@ -100,20 +107,40 @@ public: if(segments.empty()) { #ifdef CGAL_AW3_DEBUG - std::cout << "Warning: Input is empty " << std::endl; + std::cout << "Warning: Input is empty (SS)" << std::endl; #endif return; } - const std::size_t old_size = m_segments.size(); - m_segments.insert(std::cend(m_segments), std::cbegin(segments), std::cend(segments)); + typename Geom_traits::Is_degenerate_3 is_degenerate = this->geom_traits().is_degenerate_3_object(); + + const std::size_t old_size = m_segments_ptr->size(); + + for(const Segment& s : segments) + { + if(is_degenerate(s)) + { +#ifdef CGAL_AW3_DEBUG + std::cerr << "Warning: ignoring degenerate segment " << s << std::endl; +#endif + continue; + } + + m_segments_ptr->push_back(s); + } #ifdef CGAL_AW3_DEBUG std::cout << "Insert into AABB tree (segments)..." << std::endl; #endif - this->tree().insert(std::next(std::cbegin(m_segments), old_size), std::cend(m_segments)); + this->tree().insert(std::next(std::cbegin(*m_segments_ptr), old_size), std::cend(*m_segments_ptr)); - CGAL_postcondition(this->tree().size() == m_segments.size()); + // Manually constructing it here purely for profiling reasons: if we keep the lazy approach, + // it will be done at the first treatment of a facet that needs a Steiner point. + // So if one wanted to bench the flood fill runtime, it would be skewed by the time it takes + // to accelerate the tree. + this->tree().accelerate_distance_queries(); + + CGAL_postcondition(this->tree().size() == m_segments_ptr->size()); } }; diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Triangle_mesh_oracle.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Triangle_mesh_oracle.h index c87f82ac75f..d7f325e1428 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Triangle_mesh_oracle.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Triangle_mesh_oracle.h @@ -135,7 +135,7 @@ public: if(is_empty(tmesh)) { #ifdef CGAL_AW3_DEBUG - std::cout << "Warning: Input is empty " << std::endl; + std::cout << "Warning: Input is empty (TM)" << std::endl; #endif return; } @@ -153,7 +153,12 @@ public: for(face_descriptor f : faces(tmesh)) { if(Polygon_mesh_processing::is_degenerate_triangle_face(f, tmesh, np)) + { +#ifdef CGAL_AW3_DEBUG + std::cerr << "Warning: ignoring degenerate face " << f << std::endl; +#endif continue; + } const Point_ref p0 = get(vpm, source(halfedge(f, tmesh), tmesh)); const Point_ref p1 = get(vpm, target(halfedge(f, tmesh), tmesh)); @@ -164,6 +169,12 @@ public: Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits()); } + // Manually constructing it here purely for profiling reasons: if we keep the lazy approach, + // it will be done at the first treatment of a facet that needs a Steiner point. + // So if one wanted to bench the flood fill runtime, it would be skewed by the time it takes + // to accelerate the tree. + this->tree().accelerate_distance_queries(); + #ifdef CGAL_AW3_DEBUG std::cout << "Tree: " << this->tree().size() << " primitives (" << num_faces(tmesh) << " faces in input)" << std::endl; #endif diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Triangle_soup_oracle.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Triangle_soup_oracle.h index 0a8f589fc2d..dadfb5d8be1 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Triangle_soup_oracle.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Triangle_soup_oracle.h @@ -133,7 +133,7 @@ public: if(points.empty() || faces.empty()) { #ifdef CGAL_AW3_DEBUG - std::cout << "Warning: Input is empty " << std::endl; + std::cout << "Warning: Input is empty (TS)" << std::endl; #endif return; } @@ -164,11 +164,22 @@ public: const Triangle_3 tr = triangle(p0, p1, p2); if(is_degenerate(tr)) + { +#ifdef CGAL_AW3_DEBUG + std::cerr << "Warning: ignoring degenerate face " << tr << std::endl; +#endif continue; + } Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits()); } + // Manually constructing it here purely for profiling reasons: if we keep the lazy approach, + // it will be done at the first treatment of a facet that needs a Steiner point. + // So if one wanted to bench the flood fill runtime, it would be skewed by the time it takes + // to accelerate the tree. + this->tree().accelerate_distance_queries(); + #ifdef CGAL_AW3_DEBUG std::cout << "Tree: " << this->tree().size() << " primitives (" << faces.size() << " faces in input)" << std::endl; #endif @@ -179,12 +190,31 @@ public: void add_triangle_soup(const TriangleRange& triangles, const CGAL_NP_CLASS& /*np*/ = CGAL::parameters::default_values()) { + if(triangles.empty()) + { +#ifdef CGAL_AW3_DEBUG + std::cout << "Warning: Input is empty (TS)" << std::endl; +#endif + return; + } + +#ifdef CGAL_AW3_DEBUG + std::cout << "Insert into AABB Tree (triangles)..." << std::endl; +#endif + typename Geom_traits::Is_degenerate_3 is_degenerate = this->geom_traits().is_degenerate_3_object(); + Splitter_base::reserve(triangles.size()); + for(const Triangle_3& tr : triangles) { if(is_degenerate(tr)) + { +#ifdef CGAL_AW3_DEBUG + std::cerr << "Warning: ignoring degenerate triangle " << tr << std::endl; +#endif continue; + } Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits()); } diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/gate_priority_queue.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/gate_priority_queue.h index 8d63e34c9e3..ec48b4f458b 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/gate_priority_queue.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/gate_priority_queue.h @@ -27,27 +27,29 @@ namespace CGAL { namespace Alpha_wraps_3 { namespace internal { +#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + // Represents an alpha-traversable facet in the mutable priority queue -template +template class Gate { - using Facet = typename DT3::Facet; - using FT = typename DT3::Geom_traits::FT; + using Facet = typename Tr::Facet; + using FT = typename Tr::Geom_traits::FT; private: Facet m_facet; FT m_priority; // circumsphere sq_radius - bool m_is_artificial_facet; + bool m_is_permissive_facet; public: // Constructors Gate(const Facet& facet, const FT& priority, - const bool is_artificial_facet) + const bool is_permissive_facet) : m_facet(facet), m_priority(priority), - m_is_artificial_facet(is_artificial_facet) + m_is_permissive_facet(is_permissive_facet) { CGAL_assertion(priority >= 0); } @@ -60,34 +62,85 @@ public: public: const Facet& facet() const { return m_facet; } const FT& priority() const { return m_priority; } - bool is_artificial_facet() const { return m_is_artificial_facet; } + bool is_permissive_facet() const { return m_is_permissive_facet; } }; struct Less_gate { - template - bool operator()(const Gate& a, const Gate& b) const + template + bool operator()(const Gate& a, const Gate& b) const { - // @fixme? make it a total order by comparing addresses if both gates are bbox facets - if(a.is_artificial_facet()) - return true; - else if(b.is_artificial_facet()) - return false; + // If one is permissive and the other is not, give priority to the permissive facet. + // + // The permissive facet are given highest priority because they need to be treated + // regardless of their circumradius. Treating them first allow the part that depends + // on alpha to be treated uniformly in a way: whatever the alpha, all permissive faces + // will first be treated. + if(a.is_permissive_facet() != b.is_permissive_facet()) + return a.is_permissive_facet(); + + if(a.priority() == b.priority()) + { + // arbitrary, the sole purpose is to make it a total order for determinism + if(a.facet().first->time_stamp() == b.facet().first->time_stamp()) + return a.facet().second < b.facet().second; + + return a.facet().first->time_stamp() < b.facet().first->time_stamp(); + } + return a.priority() > b.priority(); } }; -template +#else // CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + +// Represents an alpha-traversable facet in the mutable priority queue +template +class Gate +{ + using Facet = typename Tr::Facet; + using FT = typename Tr::Geom_traits::FT; + +private: + Facet m_facet, m_mirror_facet; + const unsigned int m_erase_counter_mem; + const unsigned int m_mirror_erase_counter_mem; + +public: + // Constructors + Gate(const Facet& facet, + const Tr& tr) + : + m_facet(facet), + m_mirror_facet(tr.mirror_facet(facet)), + m_erase_counter_mem(m_facet.first->erase_counter()), + m_mirror_erase_counter_mem(m_mirror_facet.first->erase_counter()) + { + } + +public: + const Facet& facet() const { return m_facet; } + + bool is_zombie() const + { + return (m_facet.first->erase_counter() != m_erase_counter_mem) || + (m_mirror_facet.first->erase_counter() != m_mirror_erase_counter_mem); + } +}; + +#endif // CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + +template struct Gate_ID_PM { - using key_type = Gate; + using key_type = Gate; using value_type = std::size_t; using reference = std::size_t; using category = boost::readable_property_map_tag; inline friend value_type get(Gate_ID_PM, const key_type& k) { - using Facet = typename DT3::Facet; + using Facet = typename Tr::Facet; const Facet& f = k.facet(); return (4 * f.first->time_stamp() + f.second); diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/geometry_utils.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/geometry_utils.h index d3814d0f3b2..7f5f60de701 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/geometry_utils.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/geometry_utils.h @@ -40,16 +40,16 @@ struct Orientation_of_circumcenter } }; -template +template bool -less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha, - const typename Dt::Facet& fh, - const Dt& dt) +less_squared_radius_of_min_empty_sphere(typename Tr::Geom_traits::FT sq_alpha, + const typename Tr::Facet& fh, + const Tr& tr) { - using Cell_handle = typename Dt::Cell_handle; - using Point = typename Dt::Point; + using Cell_handle = typename Tr::Cell_handle; + using Point = typename Tr::Point; - using CK = typename Dt::Geom_traits; + using CK = typename Tr::Geom_traits; using Exact_kernel = typename Exact_kernel_selector::Exact_kernel; using Approximate_kernel = Simple_cartesian; using C2A = Cartesian_converter; @@ -61,21 +61,30 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha, Orientation_of_circumcenter orientation_of_circumcenter; +#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY + std::cout << "Checking for traversability of facet" << std::endl; +#endif + const Cell_handle c = fh.first; const int ic = fh.second; const Cell_handle n = c->neighbor(ic); - const Point& p1 = dt.point(c, Dt::vertex_triple_index(ic,0)); - const Point& p2 = dt.point(c, Dt::vertex_triple_index(ic,1)); - const Point& p3 = dt.point(c, Dt::vertex_triple_index(ic,2)); + const Point& p1 = tr.point(c, Tr::vertex_triple_index(ic,0)); + const Point& p2 = tr.point(c, Tr::vertex_triple_index(ic,1)); + const Point& p3 = tr.point(c, Tr::vertex_triple_index(ic,2)); // This is not actually possible in the context of alpha wrapping, but keeping it for genericity // and because it does not cost anything. - if(dt.is_infinite(n)) + if(tr.is_infinite(n)) { +#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY + std::cerr << "Warning: computing less_squared_radius_of_min_empty_sphere() with an infinite neighbor?" << std::endl; +#endif + CGAL_assertion(!tr.is_infinite(c)); + Orientation ori = orientation_of_circumcenter(p1, p2, p3, - dt.point(c, 0), dt.point(c, 1), - dt.point(c, 2), dt.point(c, 3)); + tr.point(c, 0), tr.point(c, 1), + tr.point(c, 2), tr.point(c, 3)); if(ori == POSITIVE) { @@ -84,18 +93,22 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha, } else { - Comparison_result cr = compare_squared_radius(dt.point(c, 0), dt.point(c, 1), - dt.point(c, 2), dt.point(c, 3), + Comparison_result cr = compare_squared_radius(tr.point(c, 0), tr.point(c, 1), + tr.point(c, 2), tr.point(c, 3), sq_alpha); return cr == LARGER; } } - if(dt.is_infinite(c)) + if(tr.is_infinite(c)) { Orientation ori = orientation_of_circumcenter(p1, p2, p3, - dt.point(n, 0), dt.point(n, 1), - dt.point(n, 2), dt.point(n, 3)); + tr.point(n, 0), tr.point(n, 1), + tr.point(n, 2), tr.point(n, 3)); + +#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY + std::cout << "Cell 'c' is infinite; Orientation: " << ori << std::endl; +#endif if(ori == NEGATIVE) { @@ -104,8 +117,8 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha, } else { - Comparison_result cr = compare_squared_radius(dt.point(n, 0), dt.point(n, 1), - dt.point(n, 2), dt.point(n, 3), + Comparison_result cr = compare_squared_radius(tr.point(n, 0), tr.point(n, 1), + tr.point(n, 2), tr.point(n, 3), sq_alpha); return cr == LARGER; } @@ -113,40 +126,40 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha, // both c and n are finite if(orientation_of_circumcenter(p1, p2, p3, - dt.point(c, 0), dt.point(c, 1), dt.point(c, 2), dt.point(c, 3)) != + tr.point(c, 0), tr.point(c, 1), tr.point(c, 2), tr.point(c, 3)) != orientation_of_circumcenter(p1, p2, p3, - dt.point(n, 0), dt.point(n, 1), dt.point(n, 2), dt.point(n, 3))) + tr.point(n, 0), tr.point(n, 1), tr.point(n, 2), tr.point(n, 3))) { Comparison_result cr = compare_squared_radius(p1, p2, p3, sq_alpha); #ifdef CGAL_AW3_DEBUG_TRAVERSABILITY std::cout << "dual crosses the face; CR: " - << typename Dt::Geom_traits().compute_squared_radius_3_object()(p1, p2, p3) + << typename Tr::Geom_traits().compute_squared_radius_3_object()(p1, p2, p3) << " sq alpha " << sq_alpha << std::endl; #endif return cr == LARGER; } else { - Comparison_result cr = compare_squared_radius(dt.point(c, 0), dt.point(c, 1), - dt.point(c, 2), dt.point(c, 3), + Comparison_result cr = compare_squared_radius(tr.point(c, 0), tr.point(c, 1), + tr.point(c, 2), tr.point(c, 3), sq_alpha); #ifdef CGAL_AW3_DEBUG_TRAVERSABILITY std::cout << "dual does not cross the face; CR(c): " - << typename Dt::Geom_traits().compute_squared_radius_3_object()(dt.point(c, 0), dt.point(c, 1), - dt.point(c, 2), dt.point(c, 3)) + << typename Tr::Geom_traits().compute_squared_radius_3_object()(tr.point(c, 0), tr.point(c, 1), + tr.point(c, 2), tr.point(c, 3)) << " sq alpha " << sq_alpha << std::endl; #endif if(cr != LARGER) return false; - cr = compare_squared_radius(dt.point(n, 0), dt.point(n, 1), - dt.point(n, 2), dt.point(n, 3), + cr = compare_squared_radius(tr.point(n, 0), tr.point(n, 1), + tr.point(n, 2), tr.point(n, 3), sq_alpha); #ifdef CGAL_AW3_DEBUG_TRAVERSABILITY std::cout << "dual does not cross the face; CR(n): " - << typename Dt::Geom_traits().compute_squared_radius_3_object()(dt.point(n, 0), dt.point(n, 1), - dt.point(n, 2), dt.point(n, 3)) + << typename Tr::Geom_traits().compute_squared_radius_3_object()(tr.point(n, 0), tr.point(n, 1), + tr.point(n, 2), tr.point(n, 3)) << " sq alpha " << sq_alpha << std::endl; #endif @@ -154,6 +167,100 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha, } } +template +typename Tr::Geom_traits::FT +smallest_squared_radius_3(const typename Tr::Facet& fh, + const Tr& tr) +{ + using Cell_handle = typename Tr::Cell_handle; + using Point = typename Tr::Point; + using FT = typename Tr::Geom_traits::FT; + + using CK = typename Tr::Geom_traits; + using Exact_kernel = typename Exact_kernel_selector::Exact_kernel; + using Approximate_kernel = Simple_cartesian; + using C2A = Cartesian_converter; + using C2E = typename Exact_kernel_selector::C2E; + + using Orientation_of_circumcenter = Filtered_predicate, + Orientation_of_circumcenter, + C2E, C2A>; + + Orientation_of_circumcenter orientation_of_circumcenter; + + auto squared_radius = tr.geom_traits().compute_squared_radius_3_object(); + +#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY + std::cout << "Computing circumradius of facet" << std::endl; +#endif + + CGAL_precondition(!tr.is_infinite(fh)); + + const Cell_handle c = fh.first; + const int ic = fh.second; + const Cell_handle n = c->neighbor(ic); + + const Point& p1 = tr.point(c, Tr::vertex_triple_index(ic,0)); + const Point& p2 = tr.point(c, Tr::vertex_triple_index(ic,1)); + const Point& p3 = tr.point(c, Tr::vertex_triple_index(ic,2)); + + // This is not actually possible in the context of alpha wrapping, but keeping it for genericity + // and because it does not cost anything. + if(tr.is_infinite(n)) + { + CGAL_assertion(!tr.is_infinite(c)); + + Orientation ori = orientation_of_circumcenter(p1, p2, p3, + tr.point(c, 0), tr.point(c, 1), + tr.point(c, 2), tr.point(c, 3)); + if(ori == POSITIVE) + return squared_radius(p1, p2, p3); + else + return squared_radius(tr.point(c, 0), tr.point(c, 1), tr.point(c, 2), tr.point(c, 3)); + } + + if(tr.is_infinite(c)) + { + Orientation ori = orientation_of_circumcenter(p1, p2, p3, + tr.point(n, 0), tr.point(n, 1), + tr.point(n, 2), tr.point(n, 3)); + +#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY + std::cout << "Cell 'c' is infinite; Orientation: " << ori << std::endl; +#endif + + if(ori == NEGATIVE) + return squared_radius(p1, p2, p3); + else + return squared_radius(tr.point(n, 0), tr.point(n, 1), tr.point(n, 2), tr.point(n, 3)); + } + + // both c and n are finite + if(orientation_of_circumcenter(p1, p2, p3, + tr.point(c, 0), tr.point(c, 1), tr.point(c, 2), tr.point(c, 3)) != + orientation_of_circumcenter(p1, p2, p3, + tr.point(n, 0), tr.point(n, 1), tr.point(n, 2), tr.point(n, 3))) + { +#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY + std::cout << "dual crosses the face; CR: " << squared_radius(p1, p2, p3) << std::endl; +#endif + + return squared_radius(p1, p2, p3); + } + else + { + const FT cr = squared_radius(tr.point(c, 0), tr.point(c, 1), tr.point(c, 2), tr.point(c, 3)); + const FT cnr = squared_radius(tr.point(n, 0), tr.point(n, 1), tr.point(n, 2), tr.point(n, 3)); + +#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY + std::cout << "dual does not cross the face; CR(c): " << cr << " CRn: " << cnr << std::endl; +#endif + + return (CGAL::min)(cr, cnr); + } +} + + } // namespace internal } // namespace Alpha_wraps_3 } // namespace CGAL diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/alpha_wrap_validation.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/validation.h similarity index 97% rename from Alpha_wrap_3/test/Alpha_wrap_3/alpha_wrap_validation.h rename to Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/validation.h index c9445c9e9c5..6d9a9191d5f 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/alpha_wrap_validation.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/validation.h @@ -117,7 +117,7 @@ bool has_expected_Hausdorff_distance(const TriangleMesh& wrap, template bool is_valid_wrap(const TriangleMesh& wrap, - const bool check_manifoldness = true, + const bool check_manifoldness, const NamedParameters& np = parameters::default_values()) { namespace PMP = CGAL::Polygon_mesh_processing; @@ -203,6 +203,13 @@ bool is_valid_wrap(const TriangleMesh& wrap, return true; } +template +bool is_valid_wrap(const TriangleMesh& wrap, + const NamedParameters& np = parameters::default_values()) +{ + return is_valid_wrap(wrap, true /*consider manifoldness*/, np); +} + template diff --git a/Alpha_wrap_3/include/CGAL/alpha_wrap_3.h b/Alpha_wrap_3/include/CGAL/alpha_wrap_3.h index 638c5fb9fa1..8be78675422 100644 --- a/Alpha_wrap_3/include/CGAL/alpha_wrap_3.h +++ b/Alpha_wrap_3/include/CGAL/alpha_wrap_3.h @@ -105,7 +105,7 @@ void alpha_wrap_3(const PointRange& points, using NP_helper = Point_set_processing_3_np_helper; using Geom_traits = typename NP_helper::Geom_traits; using Oracle = Alpha_wraps_3::internal::Triangle_soup_oracle; - using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3; + using AW3 = Alpha_wraps_3::internal::Alpha_wrapper_3; Geom_traits gt = choose_parameter(get_parameter(in_np, internal_np::geom_traits)); @@ -254,7 +254,7 @@ void alpha_wrap_3(const TriangleMesh& tmesh, using Geom_traits = typename GetGeomTraits::type; using Oracle = Alpha_wraps_3::internal::Triangle_mesh_oracle; - using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3; + using AW3 = Alpha_wraps_3::internal::Alpha_wrapper_3; Geom_traits gt = choose_parameter(get_parameter(in_np, internal_np::geom_traits)); @@ -350,7 +350,7 @@ void alpha_wrap_3(const PointRange& points, using NP_helper = Point_set_processing_3_np_helper; using Geom_traits = typename NP_helper::Geom_traits; using Oracle = Alpha_wraps_3::internal::Point_set_oracle; - using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3; + using AW3 = Alpha_wraps_3::internal::Alpha_wrapper_3; Geom_traits gt = choose_parameter(get_parameter(in_np, internal_np::geom_traits)); diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp index d7567bab205..af94eecd25a 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp +++ b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp @@ -1,12 +1,13 @@ #define CGAL_AW3_TIMER #define CGAL_AW3_DEBUG +#define CGAL_AW3_DEBUG_MANIFOLDNESS // #define CGAL_AW3_DEBUG_INITIALIZATION #include #include #include -#include "alpha_wrap_validation.h" +#include #include @@ -27,7 +28,7 @@ void generate_random_seeds(const Oracle& oracle, Seeds& seeds, CGAL::Random& r) { - const auto bbox = CGAL::Alpha_wraps_3::internal::Alpha_wrap_3(oracle).construct_bbox(offset); + const auto bbox = CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3(oracle).construct_bbox(offset); const double sq_offset = CGAL::square(offset); while(seeds.size() < 3) @@ -69,7 +70,7 @@ void alpha_wrap_triangle_mesh(Mesh& input_mesh, Oracle oracle; oracle.add_triangle_mesh(input_mesh); - AW3::internal::Alpha_wrap_3 aw3(oracle); + AW3::internal::Alpha_wrapper_3 aw3(oracle); if(seeds.empty()) generate_random_seeds(oracle, offset, seeds, r); diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp index 638431e3056..954fe6163df 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp +++ b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp @@ -1,11 +1,12 @@ #define CGAL_AW3_TIMER #define CGAL_AW3_DEBUG +#define CGAL_AW3_DEBUG_MANIFOLDNESS //#define CGAL_AW3_DEBUG_STEINER_COMPUTATION //#define CGAL_AW3_DEBUG_INITIALIZATION //#define CGAL_AW3_DEBUG_QUEUE #include -#include "alpha_wrap_validation.h" +#include #include #include diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp index e2abc6f1f12..35a954f17b2 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp +++ b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp @@ -1,10 +1,11 @@ #define CGAL_AW3_TIMER #define CGAL_AW3_DEBUG +#define CGAL_AW3_DEBUG_MANIFOLDNESS #include #include -#include "alpha_wrap_validation.h" +#include #include #include @@ -52,7 +53,7 @@ void alpha_wrap_triangle_soup(Points& pr, // AW3 Oracle oracle; oracle.add_triangle_soup(pr, fr); - AW3::internal::Alpha_wrap_3 aw3(oracle); + AW3::internal::Alpha_wrapper_3 aw3(oracle); Mesh wrap; aw3(alpha, offset, wrap, CGAL::parameters::do_enforce_manifoldness(false)); diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp b/Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp index 470d48e40d4..6209019a7a4 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp +++ b/Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp @@ -1,12 +1,12 @@ #define CGAL_AW3_TIMER //#define CGAL_AW3_DEBUG -//#define CGAL_AW3_DEBUG_MANIFOLDNESS +#define CGAL_AW3_DEBUG_MANIFOLDNESS //#define CGAL_AW3_DEBUG_STEINER_COMPUTATION //#define CGAL_AW3_DEBUG_INITIALIZATION //#define CGAL_AW3_DEBUG_QUEUE #include -#include "alpha_wrap_validation.h" +#include #include #include diff --git a/BGL/include/CGAL/boost/graph/helpers.h b/BGL/include/CGAL/boost/graph/helpers.h index 76135fd82a7..43016f9889f 100644 --- a/BGL/include/CGAL/boost/graph/helpers.h +++ b/BGL/include/CGAL/boost/graph/helpers.h @@ -953,6 +953,11 @@ void swap_edges(const typename boost::graph_traits::halfedge_descript if (fo2 != nf && halfedge(fo2, g)==oh2) set_halfedge(fo2, oh1, g); } +template +void collect_garbage(Graph&) +{ + // nothing by default +} } //end of internal namespace diff --git a/Data/data/meshes/regular_tetrahedron.off b/Data/data/meshes/regular_tetrahedron.off new file mode 100644 index 00000000000..08a25ad43b3 --- /dev/null +++ b/Data/data/meshes/regular_tetrahedron.off @@ -0,0 +1,11 @@ +OFF +4 4 0 +-1 0 -0.707107 +1 0 -0.707107 +0 -1 0.707107 +0 1 0.707107 +3 0 1 2 +3 0 3 1 +3 0 2 3 +3 1 3 2 + diff --git a/Documentation/doc/biblio/cgal_manual.bib b/Documentation/doc/biblio/cgal_manual.bib index e16202c69dd..51dfb109004 100644 --- a/Documentation/doc/biblio/cgal_manual.bib +++ b/Documentation/doc/biblio/cgal_manual.bib @@ -778,7 +778,7 @@ Teillaud" , volume = 44 , year = 2011 , pages = "160--168" -, url = "https://theses.hal.science/inria-00560388/" +, url = "https://hal.science/inria-00560388/" , doi = "10.1016/j.comgeo.2010.09.010" } diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 6c6c5f84401..09bde261d24 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -44,6 +44,7 @@ Release date: October 2023 the mean and Gaussian curvatures, as well as the principal curvature and directions. - Added the option to use a variable sizing field for `CGAL::Polygon_mesh_processing::isotropic_remeshing()`, and a sizing function based on a measure of local curvature for adaptive remeshing. +- Added the function `CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel()` that refines a polygon mesh along an isocurve. - Added the function `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` that refines a soup of triangles so that no pair of triangles intersects in their interiors. Also added, the function `autorefine()` operating directly on a triangle mesh and updating it diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 96f7bb9dfbc..084a6578913 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -978,11 +978,42 @@ namespace CommonKernelFunctors { const Vector_3 ad = vector(a,d); const Vector_3 abad = cross_product(ab,ad); - const double x = CGAL::to_double(scalar_product(cross_product(ab,ac), abad)); - const double l_ab = CGAL::sqrt(CGAL::to_double(sq_distance(a,b))); - const double y = l_ab * CGAL::to_double(scalar_product(ac,abad)); + const Vector_3 abac = cross_product(ab,ac); - return FT(std::atan2(y, x) * 180 / CGAL_PI ); + // The dihedral angle we are interested in is the angle around the oriented + // edge ab which is the same (in absolute value) as the angle between the + // vectors ab^ac and ab^ad (cross-products). + // (abac points inside the tetra abcd if its orientation is positive and outside otherwise) + // + // We consider the vector abad in the basis defined by the three vectors + // (, , ) + // where denote the normalized vector u/|u|. + // + // In this orthonormal basis, the vector adab has the coordinates + // x = * abad + // y = * abad + // z = * abad + // We have x == 0, because abad and ab are orthogonal, and thus abad is in + // the plane (yz) of the new basis. + // + // In that basis, the dihedral angle is the angle between the y axis and abad + // which is the arctan of y/z, or atan2(z, y). + // + // (Note that ab^abac is in the plane abc, pointing outside the tetra if + // its orientation is positive and inside otherwise). + // + // For the normalization, abad appears in both scalar products + // in the quotient so we can ignore its norm. For the second + // terms of the scalar products, we are left with ab^abac and abac. + // Since ab and abac are orthogonal, the sinus of the angle between the + // two vectors is 1. + // So the norms are |ab|.|abac| vs |abac|, which is why we have a + // multiplication by |ab| in y below. + const double l_ab = CGAL::sqrt(CGAL::to_double(sq_distance(a,b))); + const double y = l_ab * CGAL::to_double(scalar_product(abac, abad)); + const double z = CGAL::to_double(scalar_product(cross_product(ab,abac),abad)); + + return FT(std::atan2(z, y) * 180 / CGAL_PI ); } }; diff --git a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp index b39ce5c01e7..dc3ccdaf989 100644 --- a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp +++ b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp @@ -10,10 +10,51 @@ struct query { double expected_angle; }; +void sign_test() +{ + K::Point_3 a(0,0,0), b(1,0,0), c(0,1, 0), d(0,0,1); + + assert( CGAL::approximate_dihedral_angle(a, b, c, d) > 0); + assert( CGAL::approximate_dihedral_angle(c, a, b, d) > 0); + assert( CGAL::approximate_dihedral_angle(a, d, b, c) > 0); + assert( CGAL::approximate_dihedral_angle(c, b, d, a) > 0); + assert( CGAL::approximate_dihedral_angle(d, b, a, c) > 0); + assert( CGAL::approximate_dihedral_angle(d, c, b, a) > 0); + + assert( CGAL::approximate_dihedral_angle(a, b, d, c) < 0); + assert( CGAL::approximate_dihedral_angle(c, a, d, b) < 0); + assert( CGAL::approximate_dihedral_angle(a, d, c, b) < 0); + assert( CGAL::approximate_dihedral_angle(c, b, a, d) < 0); + assert( CGAL::approximate_dihedral_angle(d, b, c, a) < 0); + assert( CGAL::approximate_dihedral_angle(d, c, a, b) < 0); +} + +auto almost_equal_angle(double a, double b) { + return (std::min)(std::abs(a - b), std::abs(a + 360 - b)) < 0.1; +} + +void test_regular_tetrahedron() +{ + auto half_root_of_2 = std::sqrt(2) / 2; + + // Regular tetrahedron + Point_3 a{ -1, 0, -half_root_of_2}; + Point_3 b{ 1, 0, -half_root_of_2}; + Point_3 c{ 0, 1, half_root_of_2}; + Point_3 d{ 0, -1, half_root_of_2}; + assert(orientation(a, b, c, d) == CGAL::POSITIVE); + assert(almost_equal_angle(CGAL::approximate_dihedral_angle(a, b, c, d), 70.5288)); +} + int main() { + std::cout.precision(17); + sign_test(); + test_regular_tetrahedron(); + Point_3 a = {0, 0, 0}; - Point_3 b = {0, 1, 0}; - Point_3 c = {1, 0, 0}; + Point_3 b = {0, -1, 0}; // ab is oriented so that it sees the plan xz positively. + [[maybe_unused]] Point_3 c = {1, 0, 0}; + // c can be any point in the half-plane xy, with x>0 const query queries[] = { { { 1, 0, 0}, 0.}, @@ -26,11 +67,27 @@ int main() { { { 1, 0, -1}, -45.}, }; - for(auto query: queries) { - const auto& expected = query.expected_angle; - const auto& p = query.p; - auto approx = CGAL::approximate_dihedral_angle(a, b, c, p); - std::cout << approx << " -- " << expected << '\n'; - assert( std::abs(approx - expected) < 0.1 ); + auto cnt = 0u; + for(double yc = -10; yc < 10; yc += 0.1) { + Point_3 c{1, yc, 0}; + // std::cout << "c = " << c << '\n'; + for(const auto& query : queries) { + for(double yp = -10; yp < 10; yp += 0.3) { + const auto& expected = query.expected_angle; + const Point_3 p{query.p.x(), yp, query.p.z()}; + // std::cout << "p = " << p << '\n'; + auto approx = CGAL::approximate_dihedral_angle(a, b, c, p); + // std::cout << approx << " -- " << expected << '\n'; + if(!almost_equal_angle(approx, expected)) { + std::cout << "ERROR:\n"; + std::cout << "CGAL::approximate_dihedral_angle(" << a << ", " << b << ", " << c << ", " << p << ") = " << approx << '\n'; + std::cout << "expected: " << expected << '\n'; + return 1; + } + ++cnt; + } + } } + std::cout << "OK (" << cnt << " tests)\n"; + assert(cnt > 10000); } diff --git a/Mesh_3/include/CGAL/Mesh_criteria_3.h b/Mesh_3/include/CGAL/Mesh_criteria_3.h index 96511da5520..c2d5d3d9b13 100644 --- a/Mesh_3/include/CGAL/Mesh_criteria_3.h +++ b/Mesh_3/include/CGAL/Mesh_criteria_3.h @@ -70,22 +70,15 @@ public: // is not a problem template Mesh_criteria_3_impl(const CGAL_NP_CLASS& np) - :edge_criteria_(parameters::choose_parameter(parameters::get_parameter(np, internal_np::edge_size_param), - parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::edge_sizing_field_param), - parameters::choose_parameter(parameters::get_parameter(np, internal_np::sizing_field_param), FT(DBL_MAX)))), + :edge_criteria_(parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::edge_size_param), FT(DBL_MAX)), parameters::choose_parameter(parameters::get_parameter(np, internal_np::edge_min_size_param), FT(0))), facet_criteria_(parameters::choose_parameter(parameters::get_parameter(np, internal_np::facet_angle_param), FT(0)), - parameters::choose_parameter(parameters::get_parameter(np, internal_np::facet_size_param), - parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::facet_sizing_field_param), - parameters::choose_parameter(parameters::get_parameter(np, internal_np::sizing_field_param), FT(0)))), + parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::facet_size_param), FT(0)), parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::facet_distance_param), FT(0)), parameters::choose_parameter(parameters::get_parameter(np, internal_np::facet_topology_param), CGAL::FACET_VERTICES_ON_SURFACE), parameters::choose_parameter(parameters::get_parameter(np, internal_np::facet_min_size_param), FT(0))), - cell_criteria_(parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_radius_edge_ratio_param), - parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_radius_edge_param), FT(0))), - parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_size_param), - parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::cell_sizing_field_param), - parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::sizing_field_param), FT(0)))), + cell_criteria_(parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_radius_edge_ratio_param), FT(0)), + parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::cell_size_param), FT(0)), parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_min_size_param), FT(0))) { } diff --git a/Mesh_3/test/Mesh_3/test_mesh_capsule_var_distance_bound.cpp b/Mesh_3/test/Mesh_3/test_mesh_capsule_var_distance_bound.cpp index 14390abda0d..9fa7eaad0c3 100644 --- a/Mesh_3/test/Mesh_3/test_mesh_capsule_var_distance_bound.cpp +++ b/Mesh_3/test/Mesh_3/test_mesh_capsule_var_distance_bound.cpp @@ -58,7 +58,9 @@ int main() facet_distance=field); // Mesh generation - C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, + perturb(time_limit = 0, sliver_bound = 0), + exude(time_limit = 0, sliver_bound = 0)); // // Output // std::ofstream medit_file("out.mesh"); diff --git a/Mesh_3/test/Mesh_3/test_mesh_criteria_creation.cpp b/Mesh_3/test/Mesh_3/test_mesh_criteria_creation.cpp index b81befd7a56..8436b51c64d 100644 --- a/Mesh_3/test/Mesh_3/test_mesh_criteria_creation.cpp +++ b/Mesh_3/test/Mesh_3/test_mesh_criteria_creation.cpp @@ -75,28 +75,12 @@ int main() Mc ec1(edge_size = 1); assert( ec1.edge_criteria_object().sizing_field(bp1,1,index) == 1 ); - Mc ec2(edge_sizing_field = Esf(2)); + Mc ec2(edge_size = Esf(2)); assert( ec2.edge_criteria_object().sizing_field(bp1,1,index) == 2 ); - Mc ec3(edge_sizing_field = 3.); + Mc ec3(edge_size = 3.); assert( ec3.edge_criteria_object().sizing_field(bp1,1,index) == 3 ); - Mc ec4(edge_size = 4., - edge_sizing_field = Esf(5.)); - assert( ec4.edge_criteria_object().sizing_field(bp1,1,index) == 4. ); - - Mc ec5(sizing_field = 5.); - assert( ec5.edge_criteria_object().sizing_field(bp1,1,index) == 5 ); - - Mc ec6(sizing_field = 6., - edge_sizing_field = 7.); - assert( ec6.edge_criteria_object().sizing_field(bp1,1,index) == 7. ); - - Mc ec7(sizing_field = 7., - edge_size = 8.); - assert( ec7.edge_criteria_object().sizing_field(bp1,1,index) == 8. ); - - // ----------------------------------- // Test facet criteria // ----------------------------------- @@ -109,48 +93,29 @@ int main() cp(tr.point(ch, k+3)))); FT facet_size_ok = radius_facet*FT(10); - FT facet_size_nok = radius_facet/FT(10); Mc fc1(facet_size = facet_size_ok); assert( ! fc1.facet_criteria_object()(tr, f) ); - Mc fc2(facet_sizing_field = facet_size_ok); - assert( ! fc2.facet_criteria_object()(tr, f) ); - - Mc fc3(facet_sizing_field = Fsf(facet_size_ok)); + Mc fc3(facet_size = Fsf(facet_size_ok)); assert( ! fc3.facet_criteria_object()(tr, f) ); - Mc fc4(facet_sizing_field = facet_size_nok, - facet_size = facet_size_ok); - assert( ! fc4.facet_criteria_object()(tr, f) ); - - Mc fc5(sizing_field = facet_size_ok); - assert( ! fc5.facet_criteria_object()(tr, f) ); - - Mc fc6(facet_size = facet_size_ok, - facet_sizing_field = facet_size_nok, - sizing_field = facet_size_nok); - assert( ! fc6.facet_criteria_object()(tr, f) ); - - Mc fc7(facet_sizing_field = Fsf(facet_size_ok), - sizing_field = facet_size_nok); - assert( ! fc7.facet_criteria_object()(tr, f) ); - Mc fc8(facet_distance = 8.); Mc fc9(facet_angle = 9.); Mc fc10(facet_angle = 10.1, facet_distance = 10.2, facet_size = 10.3, - facet_min_size = 0.2, - facet_sizing_field = Fsf(10.4), - sizing_field = 10.5); + facet_min_size = 0.2); + Mc fc1Ob(facet_angle = 10.11, + facet_distance = 10.12, + facet_min_size = 0.3, + facet_size = Fsf(10.14)); // Test construction from int Mc fc11(facet_size = 11); Mc fc11b(facet_size = 11, facet_min_size = 1); - Mc fc12(facet_sizing_field = 12); - Mc fc12b(facet_sizing_field = 12, facet_min_size = 2); - Mc fc13(sizing_field = 13); + Mc fc12(facet_size = Fsf(12)); + Mc fc12b(facet_size = Fsf(12), facet_min_size = 1); // Test topological criterion creation Mc fc14(facet_topology = CGAL::FACET_VERTICES_ON_SURFACE); @@ -168,46 +133,27 @@ int main() cp(tr.point(ch, 3)))); FT cell_size_ok = radius_cell*FT(10); - FT cell_size_nok = radius_cell/FT(10); Mc cc1(cell_size = cell_size_ok); assert( ! cc1.cell_criteria_object()(tr, ch) ); - Mc cc2(cell_sizing_field = cell_size_ok); - assert( ! cc2.cell_criteria_object()(tr, ch) ); - - Mc cc3(cell_sizing_field = Fsf(cell_size_ok)); + Mc cc3(cell_size = Fsf(cell_size_ok)); assert( ! cc3.cell_criteria_object()(tr, ch) ); - Mc cc4(cell_sizing_field = cell_size_nok, - cell_size = cell_size_ok); - assert( ! cc4.cell_criteria_object()(tr, ch) ); - - Mc cc5(sizing_field = cell_size_ok); - assert( ! cc5.cell_criteria_object()(tr, ch) ); - - Mc cc6(cell_size = cell_size_ok, - cell_sizing_field = cell_size_nok, - sizing_field = cell_size_nok); - assert( ! cc6.cell_criteria_object()(tr, ch) ); - - Mc cc7(cell_sizing_field = Csf(cell_size_ok), - sizing_field = cell_size_nok); + Mc cc7(cell_size = Csf(cell_size_ok)); assert( ! cc7.cell_criteria_object()(tr, ch) ); Mc cc8(cell_radius_edge_ratio = 8.); Mc cc9(cell_radius_edge_ratio = 9.1, - sizing_field = Csf(9.2) ); + cell_size = Csf(9.2) ); Mc cc10(cell_radius_edge_ratio = 10.1, cell_size = 10.2, + cell_min_size = 0.1); + Mc cc10b(cell_radius_edge_ratio = 10.1, cell_min_size = 0.1, - cell_sizing_field = Csf(10.3), - sizing_field = 10.4); + cell_size = Csf(10.3)); // Test construction from int Mc cc11(cell_size = 11); Mc cc11b(cell_size = 11, cell_min_size = 1); - Mc cc12(cell_sizing_field = 12); - Mc cc12b(cell_sizing_field = 12, cell_min_size = 2); - Mc cc13(sizing_field = 13); } diff --git a/Nef_S2/include/CGAL/Nef_S2/Sphere_segment.h b/Nef_S2/include/CGAL/Nef_S2/Sphere_segment.h index e6ff36dd0ce..a526ec5085d 100644 --- a/Nef_S2/include/CGAL/Nef_S2/Sphere_segment.h +++ b/Nef_S2/include/CGAL/Nef_S2/Sphere_segment.h @@ -33,7 +33,10 @@ template class Sphere_segment_rep typedef Sphere_segment_rep Rep; friend class Sphere_segment; public: -Sphere_segment_rep() { ps_ = pt_ = Point(); c_ = Circle(); } + +Sphere_segment_rep() : + ps_(), pt_(), c_() +{} Sphere_segment_rep(const Point& p1, const Point& p2, bool shorter_arc=true) : diff --git a/Periodic_3_mesh_3/test/Periodic_3_mesh_3/test_implicit_shapes_bunch.cpp b/Periodic_3_mesh_3/test/Periodic_3_mesh_3/test_implicit_shapes_bunch.cpp index f09fffd0d4f..87669ceb91d 100644 --- a/Periodic_3_mesh_3/test/Periodic_3_mesh_3/test_implicit_shapes_bunch.cpp +++ b/Periodic_3_mesh_3/test/Periodic_3_mesh_3/test_implicit_shapes_bunch.cpp @@ -230,7 +230,7 @@ int main() Periodic_mesh_criteria criteria(facet_angle = 30, facet_size = 0.05, facet_distance = 0.025, - cell_radius_edge = 2, + cell_radius_edge_ratio = 2, cell_size = 0.05); // Mesh generation diff --git a/Poisson_surface_reconstruction_3/include/CGAL/Reconstruction_triangulation_3.h b/Poisson_surface_reconstruction_3/include/CGAL/Reconstruction_triangulation_3.h index 9f605978f7c..5b625f8042e 100644 --- a/Poisson_surface_reconstruction_3/include/CGAL/Reconstruction_triangulation_3.h +++ b/Poisson_surface_reconstruction_3/include/CGAL/Reconstruction_triangulation_3.h @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -148,6 +149,12 @@ struct Reconstruction_triangulation_default_geom_traits_3 : public BaseGt }; +template < class BaseGt > +struct Triangulation_structural_filtering_traits > { + typedef typename Triangulation_structural_filtering_traits::Use_structural_filtering_tag Use_structural_filtering_tag; +}; + + /// \internal /// The Reconstruction_triangulation_3 class /// provides the interface requested by the Poisson_reconstruction_function class: diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSizingField.h b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSizingField.h index 8641c30c6ca..b1216f3281e 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSizingField.h +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSizingField.h @@ -16,6 +16,7 @@ class PMPSizingField{ public: /// @name Types +/// These types are used for the documentation of the functions of the concept and not needed implementation wise. /// @{ /// Vertex descriptor type @@ -38,27 +39,30 @@ typedef unspecified_type FT; /// @name Functions /// @{ -/// a function that returns the sizing value at `v`. -FT at(const vertex_descriptor v) const; +/// returns the sizing value at `v` (used during tangential relaxation). +FT at(const vertex_descriptor v, const PolygonMesh& pmesh) const; -/// a function controlling edge split and edge collapse, -/// returning the ratio of the current edge length and the local target edge length between -/// the points of `va` and `vb` in case the current edge is too long, and `std::nullopt` otherwise. +/// returns the ratio of the current edge squared length and the local target edge squared length between +/// the points of `va` and `vb` in case the current edge is too long, and `std::nullopt` otherwise +/// (used for triggering edge splits and preventing some edge collapses). std::optional is_too_long(const vertex_descriptor va, - const vertex_descriptor vb) const; + const vertex_descriptor vb, + const PolygonMesh& pmesh) const; -/// a function controlling edge collapse by returning the ratio of the squared length of `h` and the -/// local target edge length if it is too short, and `std::nullopt` otherwise. +/// returns the ratio of the squared length of `h` and the +/// local target edge squared length if it is too short, and `std::nullopt` otherwise +/// (used for triggering edge collapses). std::optional is_too_short(const halfedge_descriptor h, const PolygonMesh& pmesh) const; -/// a function returning the location of the split point of the edge of `h`. +/// returns the position of the new vertex created when splitting the edge of `h`. Point_3 split_placement(const halfedge_descriptor h, const PolygonMesh& pmesh) const; -/// a function that updates the sizing field value at the vertex `v`. -void update(const vertex_descriptor v, - const PolygonMesh& pmesh); +/// function called after the addition of the split vertex `v` in `pmesh`. +/// This function can be used for example to update a pre-computed sizing field. +void register_split_vertex(const vertex_descriptor v, + const PolygonMesh& pmesh); /// @} }; diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index 5acb6148714..ad143bf835d 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -270,6 +270,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - `CGAL::Polygon_mesh_processing::triangle()` - `CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces()` - `CGAL::Polygon_mesh_processing::detect_corners_of_regions()` +- `CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel()` \cgalCRPSection{I/O Functions} - \link PMP_IO_grp `CGAL::Polygon_mesh_processing::IO::read_polygon_mesh()`\endlink diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 24151e1533d..f7062cd3a54 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -51,6 +51,7 @@ create_single_source_cgal_program("match_faces.cpp") create_single_source_cgal_program("cc_compatible_orientations.cpp") create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp") create_single_source_cgal_program("hausdorff_bounded_error_distance_example.cpp") +create_single_source_cgal_program("isotropic_remeshing_with_custom_sizing_example.cpp") create_single_source_cgal_program("triangle_mesh_autorefinement.cpp") create_single_source_cgal_program("soup_autorefinement.cpp") @@ -77,6 +78,8 @@ if(TARGET CGAL::Eigen3_support) target_link_libraries(delaunay_remeshing_example PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("remesh_almost_planar_patches.cpp") target_link_libraries(remesh_almost_planar_patches PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("geodesic_isolevel_refinement.cpp") + target_link_libraries(geodesic_isolevel_refinement PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("interpolated_corrected_curvatures_SM.cpp") target_link_libraries(interpolated_corrected_curvatures_SM PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("interpolated_corrected_curvatures_PH.cpp") diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp index c7e1fbdb21c..949fcc8dd64 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp @@ -7,14 +7,14 @@ #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef K::Point_3 Point; -typedef K::Vector_3 Vector; +typedef K::Point_3 Point; +typedef K::Vector_3 Vector; -typedef CGAL::Surface_mesh Surface_mesh; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; -typedef boost::graph_traits::face_descriptor face_descriptor; +typedef CGAL::Surface_mesh Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::face_descriptor face_descriptor; namespace PMP = CGAL::Polygon_mesh_processing; @@ -22,7 +22,7 @@ int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off"); - Surface_mesh mesh; + Mesh mesh; if(!PMP::IO::read_polygon_mesh(filename, mesh)) { std::cerr << "Invalid input." << std::endl; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp index d21d675eb2b..4418c6096d1 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp @@ -16,9 +16,9 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef K::Point_3 Point; typedef K::Vector_3 Vector; -typedef CGAL::Polyhedron_3 Polyhedron; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; -typedef boost::graph_traits::face_descriptor face_descriptor; +typedef CGAL::Polyhedron_3 Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::face_descriptor face_descriptor; namespace PMP = CGAL::Polygon_mesh_processing; @@ -26,7 +26,7 @@ int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off"); - Polyhedron mesh; + Mesh mesh; if(!PMP::IO::read_polygon_mesh(filename, mesh)) { std::cerr << "Invalid input." << std::endl; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_LCC.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_LCC.cpp index 232164a1b95..d4d774d1105 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_LCC.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_LCC.cpp @@ -12,7 +12,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::Point_3 Point; typedef CGAL::Linear_cell_complex_traits<3, Kernel> MyTraits; -typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper<2, 3, MyTraits>::type LCC; +typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper<2, 3, MyTraits>::type Mesh; namespace PMP = CGAL::Polygon_mesh_processing; @@ -21,7 +21,7 @@ int main(int argc, char* argv[]) const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off"); - LCC mesh1, mesh2; + Mesh mesh1, mesh2; if(!PMP::IO::read_polygon_mesh(filename1, mesh1) || !PMP::IO::read_polygon_mesh(filename2, mesh2)) { std::cerr << "Invalid input." << std::endl; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/extrude.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/extrude.cpp index aa7c2bd40ea..baaa0612452 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/extrude.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/extrude.cpp @@ -8,13 +8,13 @@ #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -typedef CGAL::Surface_mesh SM; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef CGAL::Surface_mesh Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; -typedef Kernel::Vector_3 Vector; -typedef boost::property_map::type VPMap; -typedef SM::template Property_map VNMap; +typedef Kernel::Vector_3 Vector; +typedef boost::property_map::type VPMap; +typedef Mesh::template Property_map VNMap; struct Bottom { @@ -52,7 +52,7 @@ struct Top int main(int argc, char* argv[]) { - SM in, out; + Mesh in, out; std::string filename = (argc > 1) ? std::string(argv[1]) : CGAL::data_file_path("meshes/cube-ouvert.off"); double vlen = (argc > 2) ? std::stod(argv[2]) : 0.1; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/geodesic_isolevel_refinement.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/geodesic_isolevel_refinement.cpp new file mode 100644 index 00000000000..29f43453e87 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/geodesic_isolevel_refinement.cpp @@ -0,0 +1,72 @@ +#include +#include + +#include + +#include +#include + +#include +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef Kernel::Point_3 Point_3; +typedef CGAL::Surface_mesh Triangle_mesh; + +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::edge_descriptor edge_descriptor; + +typedef Triangle_mesh::Property_map Vertex_distance_map; +typedef CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3 Heat_method; + +int main(int argc, char* argv[]) +{ + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off"); + + Triangle_mesh tm; + if(!CGAL::IO::read_polygon_mesh(filename, tm) || + CGAL::is_empty(tm) || !CGAL::is_triangle_mesh(tm)) + { + std::cerr << "Invalid input file." << std::endl; + return EXIT_FAILURE; + } + + // default isovalues for cutting the mesh + std::vector isovalues = {0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 10}; + + if (argc>2) + { + isovalues.clear(); + for (int i=2; i("v:distance", 0).first; + + Heat_method hm(tm); + + //use heat method to compute approximated geodesic distances to the source vertex `s` + vertex_descriptor s = *(vertices(tm).first); + hm.add_source(s); + hm.estimate_geodesic_distances(vertex_distance); + + // property map to flag new cut edge added in the mesh + auto ecm = tm.add_property_map("e:is_constrained", 0).first; + + // refine the mesh along isovalues + for (double isovalue : isovalues) + CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(tm, vertex_distance, isovalue, CGAL::parameters::edge_is_constrained_map(ecm)); + + // split the mesh in connected components bounded by the isocurves + std::vector edges_split; + CGAL::Polygon_mesh_processing::split_connected_components(tm, edges_split, CGAL::parameters::edge_is_constrained_map(ecm)); + + assert(argc!=1 || edges_split.size() == 22); + + // export each submesh in a file + for(std::size_t i=0; i Polyhedron; +typedef CGAL::Polyhedron_3 Mesh; -typedef Polyhedron::Vertex_handle Vertex_handle; -typedef Polyhedron::Halfedge_handle Halfedge_handle; -typedef Polyhedron::Facet_handle Facet_handle; +typedef Mesh::Vertex_handle Vertex_handle; +typedef Mesh::Halfedge_handle Halfedge_handle; +typedef Mesh::Facet_handle Facet_handle; namespace PMP = CGAL::Polygon_mesh_processing; @@ -24,7 +24,7 @@ int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mech-holes-shark.off"); - Polyhedron poly; + Mesh poly; if(!PMP::IO::read_polygon_mesh(filename, poly)) { std::cerr << "Invalid input." << std::endl; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_LCC.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_LCC.cpp index be7b6166e4e..89825600510 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_LCC.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_LCC.cpp @@ -14,11 +14,11 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::Point_3 Point; typedef CGAL::Linear_cell_complex_traits<3, Kernel> MyTraits; -typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper<2, 3, MyTraits>::type LCC; +typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper<2, 3, MyTraits>::type Mesh; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; -typedef boost::graph_traits::halfedge_descriptor halfedge_descriptor; -typedef boost::graph_traits::face_descriptor face_descriptor; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::halfedge_descriptor halfedge_descriptor; +typedef boost::graph_traits::face_descriptor face_descriptor; namespace PMP = CGAL::Polygon_mesh_processing; @@ -26,7 +26,7 @@ int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mech-holes-shark.off"); - LCC mesh; + Mesh mesh; if(!PMP::IO::read_polygon_mesh(filename, mesh)) { std::cerr << "Invalid input." << std::endl; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp index 09bc213068d..678597ad255 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp @@ -7,17 +7,17 @@ #include #include -#include +#include namespace PMP = CGAL::Polygon_mesh_processing; typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; -typedef CGAL::Polyhedron_3 Polyhedron; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef CGAL::Polyhedron_3 Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; int main(int argc, char* argv[]) { - Polyhedron polyhedron; + Mesh polyhedron; const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/sphere.off"); @@ -29,11 +29,11 @@ int main(int argc, char* argv[]) } // define property map to store curvature value and directions - boost::property_map>::type + boost::property_map>::type mean_curvature_map = get(CGAL::dynamic_vertex_property_t(), polyhedron), Gaussian_curvature_map = get(CGAL::dynamic_vertex_property_t(), polyhedron); - boost::property_map>>::type + boost::property_map>>::type principal_curvatures_and_directions_map = get(CGAL::dynamic_vertex_property_t>(), polyhedron); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp index 74a678ea062..ac8378781fe 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp @@ -11,12 +11,12 @@ namespace PMP = CGAL::Polygon_mesh_processing; typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; -typedef CGAL::Surface_mesh Surface_Mesh; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef CGAL::Surface_mesh Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; int main(int argc, char* argv[]) { - Surface_Mesh smesh; + Mesh smesh; const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/sphere.off"); @@ -29,7 +29,7 @@ int main(int argc, char* argv[]) // creating and tying surface mesh property maps for curvatures (with defaults = 0) bool created = false; - Surface_Mesh::Property_map + Mesh::Property_map mean_curvature_map, Gaussian_curvature_map; boost::tie(mean_curvature_map, created) = @@ -41,7 +41,7 @@ int main(int argc, char* argv[]) assert(created); // we use a tuple of 2 scalar values and 2 vectors for principal curvatures and directions - Surface_Mesh::Property_map> + Mesh::Property_map> principal_curvatures_and_directions_map; boost::tie(principal_curvatures_and_directions_map, created) = @@ -67,4 +67,5 @@ int main(int argc, char* argv[]) << ", GC = " << Gaussian_curvature_map[v] << "\n" << ", PC = [ " << PC.min_curvature << " , " << PC.max_curvature << " ]\n"; } + return 0; } diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp index bdcdba6b7fd..dc78152059b 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp @@ -6,17 +6,18 @@ #include #include +#include namespace PMP = CGAL::Polygon_mesh_processing; typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; -typedef CGAL::Surface_mesh Surface_Mesh; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef CGAL::Surface_mesh Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; int main(int argc, char* argv[]) { // instantiating and reading mesh - Surface_Mesh smesh; + Mesh smesh; const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/sphere.off"); @@ -47,4 +48,5 @@ int main(int argc, char* argv[]) << ", GC = " << g << "\n" << ", PC = [ " << p.min_curvature << " , " << p.max_curvature << " ]\n"; } + return 0; } diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_with_custom_sizing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_with_custom_sizing_example.cpp new file mode 100644 index 00000000000..53d70af603f --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_with_custom_sizing_example.cpp @@ -0,0 +1,101 @@ +#include +#include +#include +#include +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +// a sizing fied that is increasing the size of edge along the y-axis +// starting at a minimum size at y-max and ending at a maximum size at +// y-min, with a linear interpolation of sizes in between the two extreme +// sizing values +struct My_sizing_field +{ + double min_size, max_size; + double ymin, ymax; + + My_sizing_field(double min_size, double max_size, double ymin, double ymax) + : min_size(min_size) + , max_size(max_size) + , ymin(ymin) + , ymax(ymax) + {} + + double at(K::Point_3 p) const + { + double y=p.y(); + return CGAL::square( (y-ymin)/(ymax-ymin) * (min_size - max_size) + max_size ); + } + double at(const Mesh::Vertex_index v, const Mesh& mesh) const { return at(mesh.point(v)); } + + std::optional is_too_long(const Mesh::Vertex_index va, + const Mesh::Vertex_index vb, + const Mesh& mesh) const + { + // TODO: no mesh as parameters? + K::Point_3 mp = CGAL::midpoint(mesh.point(va), mesh.point(vb)); + double sql_at = at(mp); + double sql = CGAL::squared_distance(mesh.point(va), mesh.point(vb)); + if (sql > sql_at) + return sql / sql_at; + return std::nullopt; + } + + std::optional is_too_short(const Mesh::Halfedge_index h, + const Mesh& mesh) const + { + K::Point_3 mp = CGAL::midpoint(mesh.point(source(h, mesh)), mesh.point(target(h, mesh))); + double sql_at = at(mp); + double sql = CGAL::squared_distance(mesh.point(source(h, mesh)), mesh.point(target(h, mesh))); + if (sql < sql_at) + return sql / sql_at; + return std::nullopt; + } + + K::Point_3 split_placement(const Mesh::Halfedge_index h, + const Mesh& mesh) const + { + return CGAL::midpoint(mesh.point(source(h, mesh)), mesh.point(target(h, mesh))); + } + + void register_split_vertex(const Mesh::Vertex_index, const Mesh&) {} +}; + + +int main(int argc, char* argv[]) +{ + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elk.off"); + + Mesh mesh; + if (!PMP::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) { + std::cerr << "Not a valid input file." << std::endl; + return 1; + } + + std::cout << "Start remeshing of " << filename + << " (" << num_faces(mesh) << " faces)..." << std::endl; + + CGAL::Bbox_3 bb = PMP::bbox(mesh); + My_sizing_field sizing_field(0.1, 30, bb.ymin(), bb.ymax()); + unsigned int nb_iter = 5; + + PMP::isotropic_remeshing( + faces(mesh), + sizing_field, + mesh, + CGAL::parameters::number_of_iterations(nb_iter) + .number_of_relaxation_steps(3) + ); + + CGAL::IO::write_polygon_mesh("custom_remesh_out.off", mesh, CGAL::parameters::stream_precision(17)); + + std::cout << "Remeshing done." << std::endl; + + return 0; +} diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/orient_polygon_soup_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/orient_polygon_soup_example.cpp index f603925bb6b..072ec4c6360 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/orient_polygon_soup_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/orient_polygon_soup_example.cpp @@ -15,7 +15,7 @@ #include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef CGAL::Polyhedron_3 Polyhedron; +typedef CGAL::Polyhedron_3 Mesh; // Optional visitor for orientating a polygon soup to demonstrate usage for some functions. // inherits from the default class as some functions are not overloaded @@ -59,12 +59,12 @@ int main(int argc, char* argv[]) Visitor visitor; CGAL::Polygon_mesh_processing::orient_polygon_soup(points, polygons, CGAL::parameters::visitor(visitor)); - Polyhedron mesh; + Mesh mesh; CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh); // Number the faces because 'orient_to_bound_a_volume' needs a face <--> index map int index = 0; - for(Polyhedron::Face_iterator fb=mesh.facets_begin(), fe=mesh.facets_end(); fb!=fe; ++fb) + for(Mesh::Face_iterator fb=mesh.facets_begin(), fe=mesh.facets_end(); fb!=fe; ++fb) fb->id() = index++; if(CGAL::is_closed(mesh)) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/point_inside_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/point_inside_example.cpp index 93b158ce192..1a8b207cd6b 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/point_inside_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/point_inside_example.cpp @@ -12,14 +12,14 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef K::Point_3 Point; -typedef CGAL::Polyhedron_3 Polyhedron; +typedef CGAL::Polyhedron_3 Mesh; namespace PMP = CGAL::Polygon_mesh_processing; -double max_coordinate(const Polyhedron& poly) +double max_coordinate(const Mesh& poly) { double max_coord = -std::numeric_limits::infinity(); - for(Polyhedron::Vertex_handle v : vertices(poly)) + for(Mesh::Vertex_handle v : vertices(poly)) { Point p = v->point(); max_coord = (std::max)(max_coord, CGAL::to_double(p.x())); @@ -33,14 +33,14 @@ int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off"); - Polyhedron poly; + Mesh poly; if(!PMP::IO::read_polygon_mesh(filename, poly) || CGAL::is_empty(poly) || !CGAL::is_triangle_mesh(poly)) { std::cerr << "Invalid input." << std::endl; return 1; } - CGAL::Side_of_triangle_mesh inside(poly); + CGAL::Side_of_triangle_mesh inside(poly); double size = max_coordinate(poly); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope.cpp index afbb938afb5..5cc9379992f 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope.cpp @@ -9,13 +9,13 @@ int main(int argc, char* argv[]) { typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::Point_3 Point_3; - typedef CGAL::Surface_mesh Surface_mesh; - typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef CGAL::Surface_mesh Mesh; + typedef boost::graph_traits::vertex_descriptor vertex_descriptor; typedef CGAL::Polyhedral_envelope Envelope; std::ifstream in((argc>1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off")); - Surface_mesh tmesh; + Mesh tmesh; in >> tmesh; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope_mesh_containment.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope_mesh_containment.cpp index 4e814c7dccd..ede5edf61e8 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope_mesh_containment.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope_mesh_containment.cpp @@ -13,19 +13,19 @@ int main(int argc, char* argv[]) { typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::Point_3 Point_3; - typedef CGAL::Surface_mesh Surface_mesh; + typedef CGAL::Surface_mesh Mesh; typedef CGAL::Polyhedral_envelope Envelope; std::ifstream in((argc>1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off")); - Surface_mesh tmesh; + Mesh tmesh; in >> tmesh; // remesh the input using the longest edge size as target edge length - Surface_mesh query = tmesh; - Surface_mesh::Edge_iterator longest_edge_it = + Mesh query = tmesh; + Mesh::Edge_iterator longest_edge_it = std::max_element(edges(query).begin(), edges(query).end(), - [&query](Surface_mesh::Edge_index e1, Surface_mesh::Edge_index e2) + [&query](Mesh::Edge_index e1, Mesh::Edge_index e2) { return PMP::edge_length(halfedge(e1, query), query) < PMP::edge_length(halfedge(e2, query), query); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp index 742685e7485..e043f6a9a77 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp @@ -10,9 +10,9 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef K::Point_3 Point; -typedef CGAL::Surface_mesh Surface_mesh; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; -typedef boost::graph_traits::face_descriptor face_descriptor; +typedef CGAL::Surface_mesh Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::face_descriptor face_descriptor; namespace PMP = CGAL::Polygon_mesh_processing; @@ -20,7 +20,7 @@ int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off"); - Surface_mesh mesh; + Mesh mesh; if(!PMP::IO::read_polygon_mesh(filename, mesh)) { std::cerr << "Invalid input." << std::endl; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp index 23053a6766c..30f1ed16a14 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp @@ -13,8 +13,8 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -typedef CGAL::Polyhedron_3 Polyhedron; -typedef Polyhedron::Vertex_handle Vertex_handle; +typedef CGAL::Polyhedron_3 Mesh; +typedef Mesh::Vertex_handle Vertex_handle; namespace PMP = CGAL::Polygon_mesh_processing; @@ -34,7 +34,7 @@ void extract_k_ring(Vertex_handle v, { v = qv[current_index++]; - Polyhedron::Halfedge_around_vertex_circulator e(v->vertex_begin()), e_end(e); + Mesh::Halfedge_around_vertex_circulator e(v->vertex_begin()), e_end(e); do { Vertex_handle new_v = e->opposite()->vertex(); if (D.insert(std::make_pair(new_v, dist_v + 1)).second) @@ -47,14 +47,14 @@ int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); - Polyhedron poly; + Mesh poly; if(!PMP::IO::read_polygon_mesh(filename, poly) || !CGAL::is_triangle_mesh(poly)) { std::cerr << "Invalid input." << std::endl; return 1; } - std::vector new_facets; + std::vector new_facets; std::vector new_vertices; PMP::refine(poly, faces(poly), @@ -68,7 +68,7 @@ int main(int argc, char* argv[]) refined_off.close(); std::cout << "Refinement added " << new_vertices.size() << " vertices." << std::endl; - Polyhedron::Vertex_iterator v = poly.vertices_begin(); + Mesh::Vertex_iterator v = poly.vertices_begin(); std::advance(v, 82/*e.g.*/); std::vector region; extract_k_ring(v, 12/*e.g.*/, region); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_almost_planar_patches.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_almost_planar_patches.cpp index e3924c3e095..381183d0c99 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_almost_planar_patches.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_almost_planar_patches.cpp @@ -14,13 +14,13 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::Point_3 Point_3; -typedef CGAL::Surface_mesh Surface_mesh; +typedef CGAL::Surface_mesh Mesh; namespace PMP = CGAL::Polygon_mesh_processing; int main() { - Surface_mesh sm; + Mesh sm; CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/fandisk.off"), sm); //apply a perturbation to input vertices so that points are no longer coplanar @@ -51,7 +51,7 @@ int main() edge_is_constrained_map(CGAL::make_random_access_property_map(ecm))); // run the remeshing algorithm using filled properties - Surface_mesh out; + Mesh out; PMP::remesh_almost_planar_patches(sm, out, nb_regions, nb_corners, diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_planar_patches.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_planar_patches.cpp index 9f591e6a269..92aa3b15c6f 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_planar_patches.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_planar_patches.cpp @@ -12,12 +12,12 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::Point_3 Point_3; -typedef CGAL::Surface_mesh Surface_mesh; +typedef CGAL::Surface_mesh Mesh; namespace PMP = CGAL::Polygon_mesh_processing; int main() { - Surface_mesh sm; + Mesh sm; CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/cube_quad.off"), sm); // triangulate faces; @@ -25,8 +25,8 @@ int main() std::cout << "Input mesh has " << faces(sm).size() << " faces" << std::endl; assert(faces(sm).size()==12); - Surface_mesh::Property_map ecm = - sm.add_property_map("ecm",false).first; + Mesh::Property_map ecm = + sm.add_property_map("ecm",false).first; // detect sharp edges of the cube PMP::detect_sharp_edges(sm, 60, ecm); @@ -37,7 +37,7 @@ int main() assert(faces(sm).size()>100); // decimate the mesh - Surface_mesh out; + Mesh out; PMP::remesh_planar_patches(sm, out); CGAL::IO::write_polygon_mesh("cube_decimated.off", out, CGAL::parameters::stream_precision(17)); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp index a7afc963ed7..28d770626e2 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp @@ -9,7 +9,7 @@ #include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef CGAL::Polyhedron_3 Polyhedron; +typedef CGAL::Polyhedron_3 Mesh; namespace PMP = CGAL::Polygon_mesh_processing; @@ -17,7 +17,7 @@ int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/quads_to_stitch.off"); - Polyhedron mesh; + Mesh mesh; if(!PMP::IO::read_polygon_mesh(filename, mesh)) { std::cerr << "Invalid input." << std::endl; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example_OM.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example_OM.cpp index b731e8c0220..cf97918fd3b 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example_OM.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example_OM.cpp @@ -10,8 +10,8 @@ #include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; - typedef OpenMesh::PolyMesh_ArrayKernelT< > Mesh; + int main(int argc, char* argv[]) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/quads_to_stitch.off"); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp index 7587d012363..9c7ef677bc0 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp @@ -11,7 +11,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::Point_3 Point; -typedef CGAL::Surface_mesh Surface_mesh; +typedef CGAL::Surface_mesh Mesh; namespace PMP = CGAL::Polygon_mesh_processing; @@ -20,7 +20,7 @@ int main(int argc, char* argv[]) const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/P.off"); const char* outfilename = (argc > 2) ? argv[2] : "P_tri.off"; - Surface_mesh mesh; + Mesh mesh; if(!PMP::IO::read_polygon_mesh(filename, mesh)) { std::cerr << "Error: Invalid input." << std::endl; @@ -41,7 +41,7 @@ int main(int argc, char* argv[]) PMP::triangulate_faces(mesh); // Confirm that all faces are triangles. - for(boost::graph_traits::face_descriptor f : faces(mesh)) + for(boost::graph_traits::face_descriptor f : faces(mesh)) { if(!CGAL::is_triangle(halfedge(f, mesh), mesh)) std::cerr << "Error: non-triangular face left in mesh." << std::endl; diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp index 9b246f7678d..84695b32e32 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp @@ -6,15 +6,14 @@ #include #include -#include #include #include #include typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Surface_mesh Surface_mesh; -typedef boost::graph_traits::face_descriptor face_descriptor; +typedef Kernel::Point_3 Point; +typedef CGAL::Surface_mesh Mesh; +typedef boost::graph_traits::face_descriptor face_descriptor; class Insert_iterator { @@ -41,7 +40,7 @@ public: }; -struct Visitor : public CGAL::Polygon_mesh_processing::Triangulate_faces::Default_visitor +struct Visitor : public CGAL::Polygon_mesh_processing::Triangulate_faces::Default_visitor { typedef std::unordered_map Container; @@ -78,7 +77,7 @@ int main(int argc, char* argv[]) const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/P.off"); std::ifstream input(filename); - Surface_mesh mesh; + Mesh mesh; if (!input || !(input >> mesh) || mesh.is_empty()) { std::cerr << "Not a valid off file." << std::endl; @@ -87,7 +86,7 @@ int main(int argc, char* argv[]) std::unordered_map t2q; - Surface_mesh copy; + Mesh copy; CGAL::copy_face_graph(mesh, copy, CGAL::parameters::face_to_face_output_iterator(Insert_iterator(t2q))); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/volume_connected_components.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/volume_connected_components.cpp index bd4181c3ff4..e07c803b83e 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/volume_connected_components.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/volume_connected_components.cpp @@ -12,7 +12,7 @@ #include typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -typedef CGAL::Surface_mesh Surface_mesh; +typedef CGAL::Surface_mesh Mesh; namespace PMP = CGAL::Polygon_mesh_processing; namespace params = CGAL::parameters; @@ -21,7 +21,7 @@ int main(int argc, char** argv) { const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); - Surface_mesh mesh; + Mesh mesh; if(!PMP::IO::read_polygon_mesh(filename, mesh)) { std::cerr << "Invalid input." << std::endl; @@ -29,8 +29,8 @@ int main(int argc, char** argv) } // property map to assign a volume id for each face - Surface_mesh::Property_map vol_id_map = - mesh.add_property_map().first; + Mesh::Property_map vol_id_map = + mesh.add_property_map().first; // fill the volume id map std::vector err_codes; @@ -39,7 +39,7 @@ int main(int argc, char** argv) std::cout << "Found " << nb_vol << " volumes\n"; // write each volume in an OFF file - typedef CGAL::Face_filtered_graph Filtered_graph; + typedef CGAL::Face_filtered_graph Filtered_graph; Filtered_graph vol_mesh(mesh, 0, vol_id_map); for(std::size_t id = 0; id < nb_vol; ++id) { @@ -49,11 +49,12 @@ int main(int argc, char** argv) if(id > 0) vol_mesh.set_selected_faces(id, vol_id_map); - Surface_mesh out; + Mesh out; CGAL::copy_face_graph(vol_mesh, out); std::ostringstream oss; oss << "vol_" << id <<".off"; std::ofstream os(oss.str().data()); os << out; } + return 0; } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Adaptive_sizing_field.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Adaptive_sizing_field.h index a61edc699d0..e673ac9dfbb 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Adaptive_sizing_field.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Adaptive_sizing_field.h @@ -38,6 +38,8 @@ namespace Polygon_mesh_processing * Edges too long with respect to the local target edge length are split in two, while * edges that are too short are collapsed. * +* This class depends on the \eigen library. +* * \cgalModels{PMPSizingField} * * \sa `isotropic_remeshing()` @@ -213,13 +215,13 @@ private: } public: - FT at(const vertex_descriptor v) const + FT at(const vertex_descriptor v, const PolygonMesh& /* pmesh */) const { CGAL_assertion(get(m_vertex_sizing_map, v)); return get(m_vertex_sizing_map, v); } - std::optional is_too_long(const vertex_descriptor va, const vertex_descriptor vb) const + std::optional is_too_long(const vertex_descriptor va, const vertex_descriptor vb, const PolygonMesh& /* pmesh */) const { const FT sqlen = sqlength(va, vb); FT sqtarg_len = CGAL::square(4./3. * (CGAL::min)(get(m_vertex_sizing_map, va), @@ -251,7 +253,7 @@ public: get(m_vpmap, source(h, pmesh))); } - void update(const vertex_descriptor v, const PolygonMesh& pmesh) + void register_split_vertex(const vertex_descriptor v, const PolygonMesh& pmesh) { // calculating it as the average of two vertices on other ends // of halfedges as updating is done during an edge split diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Uniform_sizing_field.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Uniform_sizing_field.h index a4144ceee65..d14f4be9666 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Uniform_sizing_field.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Uniform_sizing_field.h @@ -102,12 +102,12 @@ private: } public: - FT at(const vertex_descriptor /* v */) const + FT at(const vertex_descriptor /* v */, const PolygonMesh& /* pmesh */) const { return m_size; } - std::optional is_too_long(const vertex_descriptor va, const vertex_descriptor vb) const + std::optional is_too_long(const vertex_descriptor va, const vertex_descriptor vb, const PolygonMesh& /* pmesh */) const { const FT sqlen = sqlength(va, vb); if (sqlen > m_sq_long) @@ -133,7 +133,7 @@ public: get(m_vpmap, source(h, pmesh))); } - void update(const vertex_descriptor /* v */, const PolygonMesh& /* pmesh */) + void register_split_vertex(const vertex_descriptor /* v */, const PolygonMesh& /* pmesh */) {} private: diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index c1696428d75..30d201d1ef5 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -246,7 +246,7 @@ namespace internal { get(ecmap, e) || get(fpm, face(h,pmesh))!=get(fpm, face(opposite(h,pmesh),pmesh)) ) { - if (sizing.is_too_long(source(h, pmesh), target(h, pmesh))) + if (sizing.is_too_long(source(h, pmesh), target(h, pmesh), pmesh)) { return false; } @@ -400,7 +400,7 @@ namespace internal { for(edge_descriptor e : edge_range) { const halfedge_descriptor he = halfedge(e, mesh_); - std::optional sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_)); + std::optional sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_), mesh_); if(sqlen != std::nullopt) long_edges.emplace(he, sqlen.value()); } @@ -433,16 +433,16 @@ namespace internal { std::cout << " refinement point : " << refinement_point << std::endl; #endif //update sizing field with the new point - sizing.update(vnew, mesh_); + sizing.register_split_vertex(vnew, mesh_); //check sub-edges //if it was more than twice the "long" threshold, insert them - std::optional sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_)); + std::optional sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_), mesh_); if(sqlen_new != std::nullopt) long_edges.emplace(hnew, sqlen_new.value()); const halfedge_descriptor hnext = next(hnew, mesh_); - sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_)); + sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_), mesh_); if (sqlen_new != std::nullopt) long_edges.emplace(hnext, sqlen_new.value()); @@ -500,7 +500,7 @@ namespace internal { if (!is_split_allowed(e)) continue; const halfedge_descriptor he = halfedge(e, mesh_); - std::optional sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_)); + std::optional sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_), mesh_); if(sqlen != std::nullopt) long_edges.emplace(halfedge(e, mesh_), sqlen.value()); } @@ -550,16 +550,16 @@ namespace internal { halfedge_added(hnew_opp, status(opposite(he, mesh_))); //update sizing field with the new point - sizing.update(vnew, mesh_); + sizing.register_split_vertex(vnew, mesh_); //check sub-edges //if it was more than twice the "long" threshold, insert them - std::optional sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_)); + std::optional sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_), mesh_); if(sqlen_new != std::nullopt) long_edges.emplace(hnew, sqlen_new.value()); const halfedge_descriptor hnext = next(hnew, mesh_); - sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_)); + sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_), mesh_); if (sqlen_new != std::nullopt) long_edges.emplace(hnext, sqlen_new.value()); @@ -580,7 +580,7 @@ namespace internal { if (snew == PATCH) { - std::optional sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_)); + std::optional sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_), mesh_); if(sql != std::nullopt) long_edges.emplace(hnew2, sql.value()); } @@ -603,7 +603,7 @@ namespace internal { if (snew == PATCH) { - std::optional sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_)); + std::optional sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_), mesh_); if (sql != std::nullopt) long_edges.emplace(hnew2, sql.value()); } @@ -747,7 +747,7 @@ namespace internal { for(halfedge_descriptor ha : halfedges_around_target(va, mesh_)) { vertex_descriptor va_i = source(ha, mesh_); - std::optional sqha = sizing.is_too_long(vb, va_i); + std::optional sqha = sizing.is_too_long(vb, va_i, mesh_); if (sqha != std::nullopt) { collapse_ok = false; diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Sizing_field_base.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Sizing_field_base.h index f13135cd873..c06f0568e68 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Sizing_field_base.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Sizing_field_base.h @@ -60,13 +60,14 @@ public: typedef typename K::FT FT; public: - virtual FT at(const vertex_descriptor v) const = 0; + virtual FT at(const vertex_descriptor v, const PolygonMesh&) const = 0; virtual std::optional is_too_long(const vertex_descriptor va, - const vertex_descriptor vb) const = 0; + const vertex_descriptor vb, + const PolygonMesh&) const = 0; virtual std::optional is_too_short(const halfedge_descriptor h, const PolygonMesh& pmesh) const = 0; virtual Point_3 split_placement(const halfedge_descriptor h, const PolygonMesh& pmesh) const = 0; - virtual void update(const vertex_descriptor v, const PolygonMesh& pmesh) = 0; + virtual void register_split_vertex(const vertex_descriptor v, const PolygonMesh& pmesh) = 0; }; diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h new file mode 100644 index 00000000000..d86b470a260 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h @@ -0,0 +1,174 @@ +// Copyright (c) 2021-2023 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) : Sebastien Loriot + +#ifndef CGAL_POLYGON_MESH_PROCESSING_REFINE_MESH_AT_ISOLEVEL_H +#define CGAL_POLYGON_MESH_PROCESSING_REFINE_MESH_AT_ISOLEVEL_H + +#include + +#include +#include + +#include +#include + +namespace CGAL { +namespace Polygon_mesh_processing { + +/*! + * \ingroup PkgPolygonMeshProcessingRef + * + * refines `pm` by adding new vertices on edges having their incident vertices associated with + * values respectively larger and smaller than `isovalue` in `value_map`. + * The placement of new vertices on edges will be done by linear interpolation + * using the aforementioned values. + * New vertices will be associated `isovalue` in `value_map` when created. + * Additionally, new edges will be added by connecting new vertices created sharing + * a common incident face. Note that in case more than two new vertices are added + * on a face boundary, no edges will be created in that face. + * + * @tparam PolygonMesh a model of the concepts `EdgeListGraph` and `FaceListGraph` + * @tparam ValueMap a model of the concept `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` + * as key type and with its value type being the type of the coordinates of points associated with vertices + * in the vertex map provided to the `vertex_point_map()` named parameter. + * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" for `pm` + * + * @param pm the polygon mesh to be refined. + * @param value_map the property map containing a value at each vertex for a given function defined over the mesh. + * @param isovalue the value used to refine + * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamNBegin{edge_is_constrained_map} + * \cgalParamDescription{an output property map associating `true` to all edges connecting vertices on the isolevel, + * and `false` for all other edges.} + * \cgalParamType{a class model of `WritablePropertyMap` with `boost::graph_traits::%edge_descriptor` + * as key type and `bool` as value type} + * \cgalParamDefault{No marks on edges will be put} + * \cgalParamNEnd + * + * \cgalParamNBegin{vertex_point_map} + * \cgalParamDescription{a property map associating points to the vertices of `pm`} + * \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` + * as key type and `%Point_3` as value type} + * \cgalParamDefault{`boost::get(CGAL::vertex_point, pm)`} + * \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` + * must be available in `PolygonMesh`.} + * \cgalParamNEnd + * \cgalNamedParamsEnd + * + */ +template +void refine_mesh_at_isolevel(PolygonMesh& pm, + ValueMap value_map, + typename boost::property_traits::value_type isovalue, + const NamedParameters& np = parameters::default_values()) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::property_traits::value_type FT; + + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::is_default_parameter; + + typedef Static_boolean_property_map Default_ECM; + typedef typename internal_np::Lookup_named_param_def::type ECM; + typedef typename GetVertexPointMap < PolygonMesh, NamedParameters>::type VPM; + + VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, pm)); + + ECM ecm = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), Default_ECM()); + + std::unordered_map > faces_to_split; + std::vector to_split; + std::unordered_set vertices_on_isoline; + for (edge_descriptor e : edges(pm)) + { + vertex_descriptor src = source(e, pm), tgt = target(e, pm); + + if (get(value_map, src)==isovalue) + { + vertices_on_isoline.insert(src); + if (get(value_map, tgt)==isovalue) + { + put(ecm, e, true); // special case for faces entirely on an isovalue + continue; + } + continue; + } + if (get(value_map, tgt)==isovalue) + { + vertices_on_isoline.insert(tgt); + continue; + } + if ( (get(value_map, tgt) < isovalue) != (get(value_map, src) < isovalue) ) + { + to_split.push_back(e); + } + } + + for (edge_descriptor e : to_split) + { + vertex_descriptor src = source(e, pm), tgt = target(e, pm); + FT ds = get(value_map, src); + FT dt = get(value_map, tgt); + FT alpha = (isovalue - dt) / (ds - dt); + halfedge_descriptor hnew = CGAL::Euler::split_edge(halfedge(e, pm), pm); + put(vpm, target(hnew, pm), barycenter(get(vpm,src), alpha, get(vpm, tgt), 1-alpha)); + put(value_map, target(hnew, pm) , isovalue); + face_descriptor f = face(hnew, pm); + if (f!=boost::graph_traits::null_face()) + faces_to_split[f].push_back(hnew); + hnew=pm.prev(opposite(hnew, pm)); + f = face(hnew, pm); + if (f!=boost::graph_traits::null_face()) + faces_to_split[f].push_back(hnew); + } + + for (vertex_descriptor vh : vertices_on_isoline) + { + for (halfedge_descriptor h : halfedges_around_target(vh, pm)) + { + face_descriptor f = face(h, pm); + if (f!=boost::graph_traits::null_face()) + faces_to_split[f].push_back(h); + } + } + + for (const auto& p : faces_to_split) + { + if(p.second.size()!=2) continue; + + std::pair res = edge(target(p.second[0],pm), + target(p.second[1],pm), pm); + if (res.second) + { + // no split as the edge already exists (the two vertices are on the isolevel) + put(ecm, res.first, true); + continue; + } + + halfedge_descriptor hnew = CGAL::Euler::split_face(p.second[0], p.second[1], pm); + put(ecm, edge(hnew, pm), true); + } +} + +} } // end of CGAL::Polygon_mesh_processing + + +#endif diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh_planar_patches.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh_planar_patches.h index 03171af70bc..33514d6f2de 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh_planar_patches.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh_planar_patches.h @@ -299,7 +299,7 @@ template void mark_constrained_edges( - TriangleMesh& tm, + const TriangleMesh& tm, EdgeIsConstrainedMap edge_is_constrained, double coplanar_cos_threshold, const VertexPointMap& vpm) @@ -319,7 +319,7 @@ template std::size_t mark_corner_vertices( - TriangleMesh& tm, + const TriangleMesh& tm, EdgeIsConstrainedMap& edge_is_constrained, VertexCornerIdMap& vertex_corner_id, double coplanar_cos_threshold, @@ -546,7 +546,7 @@ template std::pair -tag_corners_and_constrained_edges(TriangleMesh& tm, +tag_corners_and_constrained_edges(const TriangleMesh& tm, double coplanar_cos_threshold, VertexCornerIdMap& vertex_corner_id, EdgeIsConstrainedMap& edge_is_constrained, diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h index 6693524324a..145d5a24135 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h @@ -288,9 +288,9 @@ void tangential_relaxation(const VertexRange& vertices, const double tri_area = gt_area(get(vpm, v), get(vpm, v1), get(vpm, v2)); const double face_weight = tri_area - / (1. / 3. * (sizing.at(v) - + sizing.at(v1) - + sizing.at(v2))); + / (1. / 3. * (sizing.at(v, tm) + + sizing.at(v1, tm) + + sizing.at(v2, tm))); weight += face_weight; const Point_3 centroid = gt_centroid(get(vpm, v), get(vpm, v1), get(vpm, v2)); diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 248831b5051..44e32386664 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -65,6 +65,7 @@ create_single_source_cgal_program("triangulate_hole_with_cdt_2_test.cpp") create_single_source_cgal_program("test_pmp_polyhedral_envelope.cpp") create_single_source_cgal_program("test_pmp_np_function.cpp") create_single_source_cgal_program("test_degenerate_pmp_clip_split_corefine.cpp") +create_single_source_cgal_program("test_isolevel_refinement.cpp") # create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp") find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_isolevel_refinement.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_isolevel_refinement.cpp new file mode 100644 index 00000000000..9e07c831435 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_isolevel_refinement.cpp @@ -0,0 +1,68 @@ +#include +#include + +#include +#include + +#include +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef Kernel::Point_3 Point_3; +typedef CGAL::Surface_mesh Triangle_mesh; + +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::edge_descriptor edge_descriptor; + +typedef Triangle_mesh::Property_map Vertex_distance_map; + +int main() +{ + const std::string filename = CGAL::data_file_path("meshes/elephant.off"); + + Triangle_mesh tm; + if(!CGAL::IO::read_polygon_mesh(filename, tm) || + CGAL::is_empty(tm) || !CGAL::is_triangle_mesh(tm)) + { + std::cerr << "Invalid input file." << std::endl; + return EXIT_FAILURE; + } + //property map for the distance values to the source set + Vertex_distance_map vertex_distance = tm.add_property_map("v:distance", 5).first; + std::vector zero_vids={ + /*a closed polyline*/ 2144,145,2690,1752,339,215,1395,338,77,2145,2052,2054,343,1936,22,1751,214,1499,142,358,2694,1750,301,65,59,2650,2060,205,2651,2061,2490,1939,898,13,298, + /*two adjacent triangles*/ 532, 185, 534, 2735, + /*another patch with missing crossed edges*/134,73,1883,2533,72,532,185,131,534 + }; + + std::vector minus_vids = {132, 364}; + + for (int i : zero_vids) + put(vertex_distance, Triangle_mesh::Vertex_index(i), 0); + for (int i : minus_vids) + put(vertex_distance, Triangle_mesh::Vertex_index(i), -5); + + // property map to flag new cut edge added in the mesh + auto ecm = tm.add_property_map("e:is_constrained", 0).first; + + CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(tm, vertex_distance, 0, CGAL::parameters::edge_is_constrained_map(ecm)); + + std::ofstream debug("edges.polylines.txt"); + for (Triangle_mesh::Edge_index e : edges(tm)) + if (get(ecm, e)) + debug << "2 " << tm.point(source(e, tm)) << " " << tm.point(target(e, tm)) << "\n"; + debug.close(); + + // split the mesh in connected components bounded by the isocurves + std::vector edges_split; + CGAL::Polygon_mesh_processing::split_connected_components(tm, edges_split, CGAL::parameters::edge_is_constrained_map(ecm)); + + // export each submesh in a file + for(std::size_t i=0; i patch_vertices; // use all param - read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); + read_poly_with_borders("elephant_quad_hole_no_DT3.off", poly, border_reps); CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole (poly, border_reps[0], CGAL::parameters:: @@ -345,7 +345,7 @@ void test_triangulate_refine_and_fair_hole_compile() { sparse_linear_solver(Default_solver())); // default solver - read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); + read_poly_with_borders("elephant_quad_hole_no_DT3.off", poly, border_reps); CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole (poly, border_reps[0], CGAL::parameters:: @@ -354,7 +354,7 @@ void test_triangulate_refine_and_fair_hole_compile() { weight_calculator(CGAL::Weights::Uniform_weight())); // default solver and weight - read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); + read_poly_with_borders("elephant_quad_hole_no_DT3.off", poly, border_reps); CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole (poly, border_reps[0], CGAL::parameters:: @@ -376,11 +376,11 @@ void generate_elephant_with_hole() { Halfedge_handle nh=opposite(halfedge(fd,poly), poly); CGAL::Euler::remove_face(halfedge(fd, poly), poly); - std::ofstream output("elephant_triangle_hole.off"); + std::ofstream output("elephant_triangle_hole_no_DT3.off"); output << poly; output.close(); CGAL::Euler::remove_face(nh, poly); - output.open("elephant_quad_hole.off"); + output.open("elephant_quad_hole_no_DT3.off"); output << poly; return; } @@ -396,8 +396,8 @@ typedef CGAL::Surface_mesh Polyhedron; generate_elephant_with_hole(); std::vector input_files; - input_files.push_back("elephant_triangle_hole.off"); - input_files.push_back("elephant_quad_hole.off"); + input_files.push_back("elephant_triangle_hole_no_DT3.off"); + input_files.push_back("elephant_quad_hole_no_DT3.off"); input_files.push_back(CGAL::data_file_path("meshes/mech-holes-shark.off")); // std::cerr.precision(15); for(std::vector::iterator it = input_files.begin(); it != input_files.end(); ++it) { diff --git a/Polyhedron/demo/Polyhedron/CGAL_double_edit.cpp b/Polyhedron/demo/Polyhedron/CGAL_double_edit.cpp index a50aef4180b..9f6f4c9b5a2 100644 --- a/Polyhedron/demo/Polyhedron/CGAL_double_edit.cpp +++ b/Polyhedron/demo/Polyhedron/CGAL_double_edit.cpp @@ -19,7 +19,6 @@ public: } QValidator::State validate ( QString & input, int & pos ) const { - fixup(input); return QDoubleValidator::validate(input, pos); } }; @@ -59,7 +58,7 @@ public: void DoubleEdit::setRange(double rmin, double rmax) { - this->validator->setRange(rmin, rmax, this->validator->decimals()); + this->validator->setRange(rmin, rmax, -1); } double DoubleEdit::getValue() diff --git a/Polyhedron/demo/Polyhedron/Plugins/Alpha_wrap_3/Alpha_wrap_3_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Alpha_wrap_3/Alpha_wrap_3_plugin.cpp index 563885a376c..91390c759dd 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Alpha_wrap_3/Alpha_wrap_3_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Alpha_wrap_3/Alpha_wrap_3_plugin.cpp @@ -9,6 +9,10 @@ #include "Scene_polylines_item.h" #include "Scene_points_with_normal_item.h" +// Since we want to do visualization and interruption, it's better to use the sorted priority queue, +// even if it is slower +#define CGAL_AW3_USE_SORTED_PRIORITY_QUEUE + #include #include @@ -35,13 +39,13 @@ using TS_Oracle = CGAL::Alpha_wraps_3::internal::Triangle_soup_oracle; using SS_Oracle = CGAL::Alpha_wraps_3::internal::Segment_soup_oracle; using Oracle = CGAL::Alpha_wraps_3::internal::Point_set_oracle; -using Wrapper = CGAL::Alpha_wraps_3::internal::Alpha_wrap_3; +using Wrapper = CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3; // Here is the pipeline for the interruption box: // - The main window is connected to a wrapping thread, which performs the wrapping. // - The wrapping has a visitor, AW3_interrupter_visitor, which has a shared_ptr to a Boolean -// - When the user clicks the box, the Boolean is switched to *false*, and the visitor throws -// - The wrapping thread catches the exception, and creates the wip mesh +// - When the user clicks the box, the Boolean is switched to *false* +// - The wrapping thread creates the wip mesh // Here is the pipeline for the iterative visualization: // - The main window is connected to a wrapping thread, which performs the wrapping. @@ -113,10 +117,10 @@ public: if(!points || !faces || !fcolors || !vcolors) return; - // If the next top of the queue has vertices on the bbox, don't draw (as to avoid producing + // If the next top of the queue has vertices on the bbox, don't draw (try to avoid producing // spikes in the visualization) // const auto& gate = wrapper.queue().top(); -// if(wrapper.triangulation().number_of_vertices() > 500 && gate.is_artificial_facet()) +// if(wrapper.triangulation().number_of_vertices() > 500 && gate.is_permissive_facet()) // return; // Skip some... @@ -201,9 +205,6 @@ public: } }; -// Use a throw to get out of the AW3 refinement loop -class Out_of_patience_exception : public std::exception { }; - template struct AW3_interrupter_visitor : BaseVisitor @@ -215,15 +216,10 @@ struct AW3_interrupter_visitor : BaseVisitor(base) { } - // Only overload this one because it gives a better state of the wrap (for other visitor calls, - // we often get tetrahedral spikes because there are artificial gates in the queue) - template - void before_Steiner_point_insertion(const Wrapper& wrapper, const Point& p) + template + constexpr bool go_further(const Wrapper&) { - if(*should_stop) - throw Out_of_patience_exception(); - - return BaseVisitor::before_Steiner_point_insertion(wrapper, p); + return !(*should_stop); } }; @@ -273,25 +269,14 @@ public: QElapsedTimer elapsed_timer; elapsed_timer.start(); - // try-catch because the stop visitor currently uses a throw - try - { - wrapper(alpha, offset, wrap, - CGAL::parameters::do_enforce_manifoldness(enforce_manifoldness) - .visitor(visitor)); + wrapper(alpha, offset, wrap, + CGAL::parameters::do_enforce_manifoldness(enforce_manifoldness) + .visitor(visitor)); + if(wrapper.queue().empty()) Q_EMIT done(this); - } - catch(const Out_of_patience_exception&) - { - if(enforce_manifoldness) - wrapper.make_manifold(); - - // extract the wrap in its current state - wrapper.extract_surface(wrap, CGAL::get(CGAL::vertex_point, wrap), !enforce_manifoldness); - + else Q_EMIT interrupted(this); - } std::cout << "Wrapping took " << elapsed_timer.elapsed() / 1000. << "s" << std::endl; } @@ -562,8 +547,6 @@ public Q_SLOTS: triangles.emplace_back(get(vpm, target(h, *pMesh)), get(vpm, target(next(h, *pMesh), *pMesh)), get(vpm, source(h, *pMesh))); - - m_wrap_bbox += triangles.back().bbox(); } continue; @@ -586,8 +569,6 @@ public Q_SLOTS: triangles.emplace_back(soup_item->points()[p[0]], soup_item->points()[p[1]], soup_item->points()[p[2]]); - - m_wrap_bbox += triangles.back().bbox(); } continue; @@ -614,8 +595,6 @@ public Q_SLOTS: triangles.emplace_back(get(vpm, target(h, *pMesh)), get(vpm, target(next(h, *pMesh), *pMesh)), get(vpm, source(h, *pMesh))); - - m_wrap_bbox += triangles.back().bbox(); } segments.reserve(segments.size() + selection_item->selected_edges.size()); @@ -623,16 +602,12 @@ public Q_SLOTS: { segments.emplace_back(get(vpm, target(halfedge(e, *pMesh), *pMesh)), get(vpm, target(opposite(halfedge(e, *pMesh), *pMesh), *pMesh))); - - m_wrap_bbox += segments.back().bbox(); } points.reserve(points.size() + selection_item->selected_vertices.size()); for(const auto& v : selection_item->selected_vertices) { points.push_back(get(vpm, v)); - - m_wrap_bbox += points.back().bbox(); } continue; @@ -671,6 +646,22 @@ public Q_SLOTS: std::cout << segments.size() << " edges" << std::endl; std::cout << points.size() << " points" << std::endl; + if(!triangles.empty()) + m_wrap_bbox = triangles.front().bbox(); + else if(!segments.empty()) + m_wrap_bbox = segments.front().bbox(); + else if(!points.empty()) + m_wrap_bbox = points.front().bbox(); + + for(const Kernel::Triangle_3& tr : triangles) + m_wrap_bbox += tr.bbox(); + for(const Kernel::Segment_3& sg : segments) + m_wrap_bbox += sg.bbox(); + for(const Kernel::Point_3& pt : points) + m_wrap_bbox += pt.bbox(); + + std::cout << "Bbox:\n" << m_wrap_bbox << std::endl; + // The relative value uses the bbox of the full scene and not that of selected items to wrap // This is intentional, both because it's tedious to make it otherwise, and because it seems // to be simpler to compare between "all wrapped" / "some wrapped" @@ -703,8 +694,8 @@ public Q_SLOTS: const bool use_message_box = ui.enableMessageBox->isChecked(); - std::cout << "Wrapping edges? " << std::boolalpha << wrap_segments << std::endl; std::cout << "Wrapping faces? " << std::boolalpha << wrap_triangles << std::endl; + std::cout << "Wrapping edges? " << std::boolalpha << wrap_segments << std::endl; if(!wrap_triangles) { @@ -756,15 +747,11 @@ public Q_SLOTS: if(alpha <= 0. || offset <= 0.) { - print_message("Warning: alpha/offset must be strictly positive - nothing to wrap"); + print_message("Warning: alpha/offset must be strictly positive"); QApplication::restoreOverrideCursor(); return; } - // Switch from 'wait' to 'busy' - QApplication::restoreOverrideCursor(); - QApplication::setOverrideCursor(Qt::BusyCursor); - for(int index : this->scene->selectionIndices()) { Scene_surface_mesh_item* sm_item = qobject_cast(this->scene->item(index)); @@ -824,6 +811,10 @@ public Q_SLOTS: // Create message box with stop button if(use_message_box) { + // Switch from 'wait' to 'busy' + QApplication::restoreOverrideCursor(); + QApplication::setOverrideCursor(Qt::BusyCursor); + m_message_box = new QMessageBox(QMessageBox::NoIcon, "Wrapping", "Wrapping in progress...", @@ -841,6 +832,8 @@ public Q_SLOTS: } // Actual start + QApplication::setOverrideCursor(Qt::WaitCursor); + wrapper_thread->start(); CGAL::Three::Three::getMutex()->lock(); diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp index 521576aa657..47f220f6e58 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp @@ -624,6 +624,8 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type, ui.approx->setRange(diag * 10e-7, // min diag); // max ui.approx->setValue(approx); + ui.approx->setToolTip(tr("Approximation error: in [%1; %2]") + .arg(diag * 10e-7).arg(diag)); ui.protect->setEnabled(features_protection_available); ui.protect->setChecked(features_protection_available); diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp index abc82f15033..d48aef5bb1f 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp @@ -134,6 +134,8 @@ Meshing_thread* cgal_code_mesh_3(QList pMeshes, param.manifold = manifold; param.protect_features = protect_features || protect_borders; param.use_sizing_field_with_aabb_tree = polylines.empty() && protect_features; + param.image_3_ptr = nullptr; + param.weights_ptr = nullptr; typedef ::Mesh_function Mesh_function; @@ -237,6 +239,8 @@ Meshing_thread* cgal_code_mesh_3(const QList pMeshes, param.manifold = manifold; param.protect_features = protect_features || protect_borders; param.use_sizing_field_with_aabb_tree = protect_features; + param.image_3_ptr = nullptr; + param.weights_ptr = nullptr; typedef ::Mesh_function Mesh_function; @@ -292,7 +296,8 @@ Meshing_thread* cgal_code_mesh_3(const Implicit_function_interface* pfunction, param.manifold = manifold; param.detect_connected_components = false; // to avoid random values // in the debug displays - + param.image_3_ptr = nullptr; + param.weights_ptr = nullptr; typedef ::Mesh_function Mesh_function; diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h index dfd31aadd84..8bc5a6c3ff6 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h @@ -138,20 +138,40 @@ QStringList Mesh_parameters:: log() const { - return QStringList() - << QString("edge max size: %1").arg(edge_sizing) - << QString("edge min size: %1").arg(edge_min_sizing) - << QString("facet min angle: %1").arg(facet_angle) - << QString("facet max size: %1").arg(facet_sizing) - << QString("facet min size: %1").arg(facet_min_sizing) - << QString("facet approx error: %1").arg(facet_approx) - << QString("tet shape (radius-edge): %1").arg(tet_shape) - << QString("tet max size: %1").arg(tet_sizing) - << QString("tet min size: %1").arg(tet_min_sizing) - << QString("detect connected components: %1") - .arg(detect_connected_components) - << QString("use weights: %1").arg(weights_ptr != nullptr) - << QString("protect features: %1").arg(protect_features); + QStringList res("Mesh criteria"); + + // doubles + if(edge_sizing > 0) + res << QString("edge max size: %1").arg(edge_sizing); + if(edge_min_sizing > 0) + res << QString("edge min size: %1").arg(edge_min_sizing); + if(facet_angle > 0) + res << QString("facet min angle: %1").arg(facet_angle); + if(facet_sizing > 0) + res << QString("facet max size: %1").arg(facet_sizing); + if(facet_min_sizing > 0) + res << QString("facet min size: %1").arg(facet_min_sizing); + if(facet_approx > 0) + res << QString("facet approx error: %1").arg(facet_approx); + if(tet_shape > 0) + res << QString("tet shape (radius-edge): %1").arg(tet_shape); + if(tet_sizing > 0) + res << QString("tet max size: %1").arg(tet_sizing); + if(tet_min_sizing > 0) + res << QString("tet min size: %1").arg(tet_min_sizing); + + // booleans + res << QString("protect features: %1").arg(protect_features); + if(image_3_ptr != nullptr) + { + res << QString("detect connected components: %1") + .arg(detect_connected_components); + res << QString("use weights: %1").arg(weights_ptr != nullptr); + } + res << QString("use aabb tree: %1").arg(use_sizing_field_with_aabb_tree); + res << QString("manifold: %1").arg(manifold); + + return res; } diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/mesh_parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/mesh_parameters_interface.h index ece342c321b..8fc2baca0c3 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/mesh_parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/mesh_parameters_interface.h @@ -55,7 +55,7 @@ perturb(const CGAL_NP_CLASS& np = parameters::default_values()) } template -Named_function_parameters<::CGAL::parameters::internal::Exude_options, ::CGAL::internal_np::exude_options_param_t, CGAL_NP_BASE> +Named_function_parameters<::CGAL::parameters::internal::Perturb_options, ::CGAL::internal_np::perturb_options_param_t, CGAL_NP_BASE> perturb(const CGAL_NP_CLASS_1& np1, const CGAL_NP_CLASS_2& np2, const NP& ... nps) { return perturb(::CGAL::internal_np::combine_named_parameters(np1, np2, nps...)); @@ -91,7 +91,7 @@ exude(const CGAL_NP_CLASS& np = parameters::default_values()) using ::CGAL::parameters::choose_parameter; using ::CGAL::parameters::get_parameter; double time_limit = choose_parameter(get_parameter(np,::CGAL::internal_np::maximum_running_time),::CGAL::parameters::internal::undef_parameter); - double sliver_bound = choose_parameter(get_parameter(np,::CGAL::internal_np::lower_sliver_bound),::CGAL::parameters::default_values_for_mesh_3::perturb_sliver_bound); + double sliver_bound = choose_parameter(get_parameter(np,::CGAL::internal_np::lower_sliver_bound),::CGAL::parameters::default_values_for_mesh_3::exude_sliver_bound); ::CGAL::parameters::internal::Exude_options options(true); @@ -304,7 +304,7 @@ mesh_3_options(const CGAL_NP_CLASS& np = parameters::default_values()) options.dump_after_refine_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_refine_prefix_param), ""); options.dump_after_glob_opt_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_glob_opt_prefix_param), ""); options.dump_after_perturb_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_perturb_prefix_param), ""); - options.dump_after_exude_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_refine_surface_prefix_param), ""); + options.dump_after_exude_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_exude_prefix_param), ""); options.number_of_initial_points=choose_parameter(get_parameter(np, ::CGAL::internal_np::number_of_initial_points_param), -1); options.nonlinear_growth_of_balls = choose_parameter(get_parameter(np, ::CGAL::internal_np::nonlinear_growth_of_balls_param), false); options.maximal_number_of_vertices=choose_parameter(get_parameter(np, ::CGAL::internal_np::maximal_number_of_vertices_param), 0); 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 e2850c63713..9b172f5b4c1 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -248,6 +248,8 @@ CGAL_add_named_parameter(smooth_constrained_edges_t, smooth_constrained_edges, s // List of named parameters used in Alpha_wrap_3 CGAL_add_named_parameter(do_enforce_manifoldness_t, do_enforce_manifoldness, do_enforce_manifoldness) CGAL_add_named_parameter(seed_points_t, seed_points, seed_points) +CGAL_add_named_parameter(refine_triangulation_t, refine_triangulation, refine_triangulation) +CGAL_add_named_parameter(keep_inner_connected_components_t, keep_inner_connected_components, keep_inner_connected_components) // SMDS_3 parameters CGAL_add_named_parameter(surface_facets_t, surface_facets, surface_facets) @@ -348,19 +350,14 @@ CGAL_add_named_parameter_with_compatibility(features_detector_param_t, features_ CGAL_add_named_parameter_with_compatibility(input_features_param_t, input_features_param, input_features) CGAL_add_named_parameter_with_compatibility(edge_size_param_t, edge_size_param, edge_size) -CGAL_add_named_parameter_with_compatibility_cref_only(edge_sizing_field_param_t, edge_sizing_field_param, edge_sizing_field) CGAL_add_named_parameter_with_compatibility(edge_min_size_param_t, edge_min_size_param, edge_min_size) CGAL_add_named_parameter_with_compatibility(facet_angle_param_t, facet_angle_param, facet_angle) CGAL_add_named_parameter_with_compatibility(facet_size_param_t, facet_size_param, facet_size) -CGAL_add_named_parameter_with_compatibility_cref_only(facet_sizing_field_param_t, facet_sizing_field_param, facet_sizing_field) CGAL_add_named_parameter_with_compatibility_cref_only(facet_distance_param_t, facet_distance_param, facet_distance) CGAL_add_named_parameter_with_compatibility(facet_min_size_param_t, facet_min_size_param, facet_min_size) CGAL_add_named_parameter_with_compatibility(facet_topology_param_t, facet_topology_param, facet_topology) -CGAL_add_named_parameter_with_compatibility(cell_radius_edge_param_t, cell_radius_edge_param, cell_radius_edge) CGAL_add_named_parameter_with_compatibility(cell_radius_edge_ratio_param_t, cell_radius_edge_ratio_param, cell_radius_edge_ratio) CGAL_add_named_parameter_with_compatibility(cell_size_param_t, cell_size_param, cell_size) -CGAL_add_named_parameter_with_compatibility_cref_only(cell_sizing_field_param_t, cell_sizing_field_param, cell_sizing_field) -CGAL_add_named_parameter_with_compatibility_cref_only(sizing_field_param_t, sizing_field_param, sizing_field) CGAL_add_named_parameter_with_compatibility(cell_min_size_param_t, cell_min_size_param, cell_min_size) CGAL_add_named_parameter_with_compatibility(function_param_t, function_param, function) diff --git a/Scripts/developer_scripts/add_toc_to_github_wiki_page.py b/Scripts/developer_scripts/add_toc_to_github_wiki_page.py index ce61f49db7e..928703c502e 100644 --- a/Scripts/developer_scripts/add_toc_to_github_wiki_page.py +++ b/Scripts/developer_scripts/add_toc_to_github_wiki_page.py @@ -1,125 +1,134 @@ -from sys import argv -from sys import exit import codecs import re import argparse +import sys +from urllib.parse import quote parser = argparse.ArgumentParser() -parser.add_argument("filename", - help="the Mardown file to process") -parser.add_argument("--codebase", - help="for a Markdown file of Codebase instead of Github", - action="store_true") -parser.add_argument("--h1", - help="support level one sections (h1)", - action="store_true") +parser.add_argument("filename", help="the Markdown file to process") +parser.add_argument( + "--codebase", + help="for a Markdown file of Codebase instead of Github", + action="store_true", +) +parser.add_argument("--h1", help="support level one sections (h1)", action="store_true") +parser.add_argument("--max-level", help="maximum level of sections", type=int, default = 5) args = parser.parse_args() + # a probably incomplete version to generate an anchor from a section name def get_anchor(s): - s = s.replace("`","") - s = s.replace("(","") - s = s.replace(")","") - s = s.replace(".","") - s = s.replace("#","") - s = s.replace(":","") - s = s.replace(",","") - s = s.replace(";","") - if args.codebase: - s = s.replace("/","-") - else: - s = s.replace("/","") - s = s.replace("<","") - s = s.replace(">","") - s = s.replace("+","") - s = s.replace("=","") - s = s.replace("?","") - s = s.replace("@","") - s = s.lstrip(" ") - s = s.rstrip("\n") - s = s.rstrip(" ") - s = re.sub(r'\s+','-',s) - if not args.codebase: - s = s.lower() - if args.codebase: - s = s.replace("'","-and-39-") - return "#"+s + s = s.replace("`", "") + s = s.replace("(", "") + s = s.replace(")", "") + s = s.replace(".", "") + s = s.replace("#", "") + s = s.replace(":", "") + s = s.replace(",", "") + s = s.replace(";", "") + if args.codebase: + s = s.replace("/", "-") + else: + s = s.replace("/", "") + s = s.replace("<", "") + s = s.replace(">", "") + s = s.replace("+", "") + s = s.replace("=", "") + s = s.replace("?", "") + s = s.replace("@", "") + s = s.lstrip(" ") + s = s.rstrip("\n") + s = s.rstrip(" ") + s = re.sub(r"\s+", "-", s) + if not args.codebase: + s = s.lower() + if args.codebase: + s = s.replace("'", "-and-39-") + return "#" + quote(s) + # indices the nesting level (first level allowed is ##) def get_level(s): - m = re.search('^(#+)\s', s) - if m: - return len(m.group(1)) - else: - return 0 + m = re.search(r"^(#+)\s", s) + if m: + return len(m.group(1)) + else: + return 0 + def get_name(s): - m = re.search('^#+\s+(.*)\s*$', s) - if m: - return m.group(1) - else: - return "ERROR: Section name extraction" - -#generate the entry for one section -def get_toc_entry(s): - name = get_name(s) - if args.h1: - level = get_level(s)-1 - else: - level = get_level(s)-2 - anchor = get_anchor(s) - - if level<0: - return "ERROR: h1 sections are not allowed" - - res="* ["+name+"]("+anchor+")" - for i in range(0,level): - res=" "+res - return res - -#now the main -input = args.filename - -f = codecs.open(input, 'r', encoding='utf-8') - -if not f: - print("Cannot open "+input+"\n") - exit() - -#look for the begin of the file -line=f.readline() -if line.find("")==-1: - exit() - -#skip current TOC -line=f.readline() -while line and line.find("")==-1: - line=f.readline() - -if not line: - exit() - -buffer="" -TOC="\n\n# Table of Contents\n" - -verbatim_mode=False # to ignore verbatim mode while looking for sections -TOC_empty=True -for line in f.readlines(): - buffer+=line - if verbatim_mode: - if line[:3]=="```": - verbatim_mode=False - else: - if line[:3]=="```": - verbatim_mode=True + m = re.search(r"^#+\s+(.*)\s*$", s) + if m: + return m.group(1) else: - if line[0]=="#": - TOC+=(get_toc_entry(line)+"\n") - TOC_empty=False -TOC+="\n\n" + return "ERROR: Section name extraction" -if not TOC_empty: - f.close() - f = codecs.open(input, 'w', encoding='utf-8') - f.write(TOC) - f.write(buffer) + +# generate the entry for one section +def get_toc_entry(s): + name = get_name(s) + if args.h1: + level = get_level(s) - 1 + else: + level = get_level(s) - 2 + anchor = get_anchor(s) + + if level < 0: + return "ERROR: h1 sections are not allowed" + + res = "* [" + name + "](" + anchor + ")" + for _ in range(0, level): + res = " " + res + return res + + +# now the main +def main(): + filename = args.filename + + f = codecs.open(filename, "r", encoding="utf-8") + + if not f: + print("Cannot open " + input + "\n") + sys.exit() + + # look for the begin of the file + line = f.readline() + if line.find("") == -1: + sys.exit() + + # skip current TOC + line = f.readline() + while line and line.find("") == -1: + line = f.readline() + + if not line: + sys.exit() + + buffer = "" + toc = "\n\n# Table of Contents\n" + + verbatim_mode = False # to ignore verbatim mode while looking for sections + toc_empty = True + for line in f.readlines(): + buffer += line + if verbatim_mode: + if line[:3] == "```": + verbatim_mode = False + else: + if line[:3] == "```": + verbatim_mode = True + else: + if line[0] == "#" and get_level(line) <= args.max_level: + toc += get_toc_entry(line) + "\n" + toc_empty = False + toc += "\n\n" + + if not toc_empty: + f.close() + f = codecs.open(filename, "w", encoding="utf-8") + f.write(toc) + f.write(buffer) + +if __name__ == "__main__": + main() diff --git a/Shape_regularization/include/CGAL/Shape_regularization/regularize_segments.h b/Shape_regularization/include/CGAL/Shape_regularization/regularize_segments.h index 2042cff97b5..ad290801019 100644 --- a/Shape_regularization/include/CGAL/Shape_regularization/regularize_segments.h +++ b/Shape_regularization/include/CGAL/Shape_regularization/regularize_segments.h @@ -70,6 +70,9 @@ namespace Segments { that should be addressed. The function is based on the class `QP_regularization`. Please refer to that class and these concepts for more information. + This class requires a `QPSolver` model which defaults to the \ref thirdpartyOSQP "OSQP" + library, which must be available on the system. + \tparam InputRange a model of `ConstRange` whose iterator type is `RandomAccessIterator` diff --git a/Spatial_searching/doc/Spatial_searching/examples.txt b/Spatial_searching/doc/Spatial_searching/examples.txt index defc30a2079..8c3c7920ff7 100644 --- a/Spatial_searching/doc/Spatial_searching/examples.txt +++ b/Spatial_searching/doc/Spatial_searching/examples.txt @@ -11,6 +11,7 @@ \example Spatial_searching/searching_with_circular_query.cpp \example Spatial_searching/searching_surface_mesh_vertices.cpp \example Spatial_searching/searching_polyhedron_vertices.cpp +\example Spatial_searching/searching_triangulation_vertices.cpp \example Spatial_searching/searching_with_point_with_info.cpp \example Spatial_searching/searching_with_point_with_info_inplace.cpp \example Spatial_searching/searching_with_point_with_info_pmap.cpp diff --git a/Spatial_searching/examples/Spatial_searching/CMakeLists.txt b/Spatial_searching/examples/Spatial_searching/CMakeLists.txt index f9546507bd4..b273ac3477b 100644 --- a/Spatial_searching/examples/Spatial_searching/CMakeLists.txt +++ b/Spatial_searching/examples/Spatial_searching/CMakeLists.txt @@ -23,6 +23,7 @@ create_single_source_cgal_program("searching_with_point_with_info_inplace.cpp") create_single_source_cgal_program("searching_with_point_with_info_pmap.cpp") create_single_source_cgal_program("searching_surface_mesh_vertices.cpp") create_single_source_cgal_program("searching_polyhedron_vertices.cpp") +create_single_source_cgal_program("searching_triangulation_vertices.cpp") create_single_source_cgal_program("searching_polyhedron_vertices_with_fuzzy_sphere.cpp") create_single_source_cgal_program("user_defined_point_and_distance.cpp") create_single_source_cgal_program("using_fair_splitting_rule.cpp") diff --git a/Spatial_searching/examples/Spatial_searching/searching_triangulation_vertices.cpp b/Spatial_searching/examples/Spatial_searching/searching_triangulation_vertices.cpp new file mode 100644 index 00000000000..0dfd38e71e1 --- /dev/null +++ b/Spatial_searching/examples/Spatial_searching/searching_triangulation_vertices.cpp @@ -0,0 +1,65 @@ +#include +#include +#include +#include +#include +#include + +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::Point_3 Point_3; +typedef CGAL::Delaunay_triangulation_3 Delaunay; +typedef Delaunay::Vertex_handle Vertex_handle; + + +struct Project { + typedef Vertex_handle argument_type; + typedef const Point_3& Point; + typedef Point result_type; + Point operator()( Vertex_handle v) const { return v->point(); } +}; + + +typedef boost::function_property_map Vertex_point_pmap; +typedef CGAL::Search_traits_3 Traits_base; +typedef CGAL::Search_traits_adapter Traits; +typedef CGAL::Orthogonal_k_neighbor_search K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef Tree::Splitter Splitter; +typedef K_neighbor_search::Distance Distance; + +int main(int argc, char* argv[]) +{ + std::ifstream in((argc>1)?argv[1]:CGAL::data_file_path("points_3/kitten.xyz")); + std::vector points; + Point_3 p; + while(in >> p){ + points.push_back(p); + } + + Delaunay dt(points.begin(), points.end()); + + const unsigned int K = 5; + + Project project; + Vertex_point_pmap vppmap(project); + + // Insert number_of_data_points in the tree + Tree tree(dt.finite_vertex_handles().begin(), dt.finite_vertex_handles().end(), Splitter(), Traits(vppmap)); + + // search K nearest neighbors + Point_3 query(0.0, 0.0, 0.0); + Distance tr_dist(vppmap); + + K_neighbor_search search(tree, query, K,0,true,tr_dist); + std::cout <<"The "<< K << " nearest vertices to the query point at (0,0,0) are:" << std::endl; + for(K_neighbor_search::iterator it = search.begin(); it != search.end(); it++){ + std::cout << "vertex " << &*(it->first) << " : " << vppmap[it->first] << " at distance " + << tr_dist.inverse_of_transformed_distance(it->second) << std::endl; + } + + return 0; +} diff --git a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h index d0a77fa1a78..a2a2228f6f5 100644 --- a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h @@ -324,8 +324,8 @@ private: CGAL_precondition( (! (face(edge,mesh)==boost::graph_traits::null_face())) && (! (face(opposite(edge,mesh),mesh)==boost::graph_traits::null_face())) ); - const Point a = get(vertex_point_pmap,target(edge,mesh)); - const Point b = get(vertex_point_pmap,target(prev(edge,mesh),mesh)); + const Point a = get(vertex_point_pmap,source(edge,mesh)); + const Point b = get(vertex_point_pmap,target(edge,mesh)); const Point c = get(vertex_point_pmap,target(next(edge,mesh),mesh)); const Point d = get(vertex_point_pmap,target(next(opposite(edge,mesh),mesh),mesh)); // As far as I check: if, say, dihedral angle is 5, this returns 175, diff --git a/TDS_3/include/CGAL/Triangulation_utils_3.h b/TDS_3/include/CGAL/Triangulation_utils_3.h index d563fda6954..80c8c507977 100644 --- a/TDS_3/include/CGAL/Triangulation_utils_3.h +++ b/TDS_3/include/CGAL/Triangulation_utils_3.h @@ -80,6 +80,7 @@ struct Triangulation_utils_3 CGAL_precondition( ( i >= 0 && i < 4 ) && ( j >= 0 && j < 4 ) && ( i != j ) ); + CGAL_assume(i!=j); return tab_next_around_edge[i][j]; } diff --git a/Triangulation_3/include/CGAL/Robust_weighted_circumcenter_filtered_traits_3.h b/Triangulation_3/include/CGAL/Robust_weighted_circumcenter_filtered_traits_3.h index 0ab43248e21..beddf4c0b3d 100644 --- a/Triangulation_3/include/CGAL/Robust_weighted_circumcenter_filtered_traits_3.h +++ b/Triangulation_3/include/CGAL/Robust_weighted_circumcenter_filtered_traits_3.h @@ -18,6 +18,7 @@ #include #include #include +#include #include namespace CGAL { @@ -531,6 +532,13 @@ public: Robust_circumcenter_filtered_traits_3(const Kernel& k = Kernel()) : Kernel(k) { } }; + +template < class BaseGt > +struct Triangulation_structural_filtering_traits > { + typedef typename Triangulation_structural_filtering_traits::Use_structural_filtering_tag Use_structural_filtering_tag; +}; + + template class Robust_weighted_circumcenter_filtered_traits_3 : public Robust_circumcenter_filtered_traits_3