This commit is contained in:
Sebastien Loriot 2025-12-03 17:09:27 +00:00 committed by GitHub
commit 26d56aa70e
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
9 changed files with 1359 additions and 42 deletions

View File

@ -28,49 +28,59 @@ namespace internal {
// triple plane intersection
template <class K>
std::optional<typename K::Point_3>
intersection_point(const typename K::Plane_3& plane1,
const typename K::Plane_3& plane2,
const typename K::Plane_3& plane3,
const K&)
intersection_point_RT(const typename K::Plane_3& plane1,
const typename K::Plane_3& plane2,
const typename K::Plane_3& plane3,
const K&)
{
typedef typename K::FT FT;
typedef typename K::RT RT;
const FT& m00 = plane1.a();
const FT& m01 = plane1.b();
const FT& m02 = plane1.c();
const FT& b0 = - plane1.d();
const FT& m10 = plane2.a();
const FT& m11 = plane2.b();
const FT& m12 = plane2.c();
const FT& b1 = - plane2.d();
const FT& m20 = plane3.a();
const FT& m21 = plane3.b();
const FT& m22 = plane3.c();
const FT& b2 = - plane3.d();
const RT& m00 = plane1.a();
const RT& m01 = plane1.b();
const RT& m02 = plane1.c();
const RT& b0 = - plane1.d();
const RT& m10 = plane2.a();
const RT& m11 = plane2.b();
const RT& m12 = plane2.c();
const RT& b1 = - plane2.d();
const RT& m20 = plane3.a();
const RT& m21 = plane3.b();
const RT& m22 = plane3.c();
const RT& b2 = - plane3.d();
// Minors common to two determinants
const FT minor_0 = m00*m11 - m10*m01;
const FT minor_1 = m00*m21 - m20*m01;
const FT minor_2 = m10*m21 - m20*m11;
const RT minor_0 = m00*m11 - m10*m01;
const RT minor_1 = m00*m21 - m20*m01;
const RT minor_2 = m10*m21 - m20*m11;
const FT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M
const RT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M
if(is_zero(den)){
return std::nullopt;
}
const FT num3 = minor_0*b2 - minor_1*b1 + minor_2*b0; // determinant of M with M[x:2] swapped with [b0,b1,b2]
const RT num3 = minor_0*b2 - minor_1*b1 + minor_2*b0; // determinant of M with M[x:2] swapped with [b0,b1,b2]
// Minors common to two determinants
const FT minor_3 = b0*m12 - b1*m02;
const FT minor_4 = b0*m22 - b2*m02;
const FT minor_5 = b1*m22 - b2*m12;
const RT minor_3 = b0*m12 - b1*m02;
const RT minor_4 = b0*m22 - b2*m02;
const RT minor_5 = b1*m22 - b2*m12;
// num1 has opposite signs because b0 and M[:1] have been swapped
const FT num1 = - minor_3*m21 + minor_4*m11 - minor_5*m01; // determinant of M with M[x:0] swapped with [b0,b1,b2]
const FT num2 = minor_3*m20 - minor_4*m10 + minor_5*m00; // determinant of M with M[x:1] swapped with [b0,b1,b2]
const RT num1 = - minor_3*m21 + minor_4*m11 - minor_5*m01; // determinant of M with M[x:0] swapped with [b0,b1,b2]
const RT num2 = minor_3*m20 - minor_4*m10 + minor_5*m00; // determinant of M with M[x:1] swapped with [b0,b1,b2]
return std::make_optional(typename K::Point_3(num1/den, num2/den, num3/den));
return typename K::Point_3(num1, num2, num3, den);
}
template <class K>
std::optional<typename K::Point_3>
intersection_point(const typename K::Plane_3& plane1,
const typename K::Plane_3& plane2,
const typename K::Plane_3& plane3,
const K& k)
{
return intersection_point_RT(plane1, plane2, plane3, k);
}
template <class K>

View File

@ -2097,18 +2097,24 @@ namespace CommonKernelFunctors {
template <typename K>
class Construct_plane_line_intersection_point_3
{
typedef typename K::Plane_3 Plane;
typedef typename K::FT FT;
typedef typename K::Line_3 Line;
typedef typename K::Point_3 Point;
typename K::Construct_plane_3 construct_plane;
typename K::Construct_line_3 construct_line;
typedef typename K::Plane_3 Plane;
typedef typename K::Vector_3 Vector;
typename K::Construct_cross_product_vector_3 construct_cross_product_vector_3;
typename K::Compute_scalar_product_3 scalar_product;
typename K::Construct_line_3 construct_line;
typename K::Construct_vector_3 vector;
typename K::Construct_plane_3 construct_plane;
public:
Point
operator()(const Point& p1, const Point& p2, const Point& p3,
const Point& l1, const Point& l2) const
{
Plane plane = construct_plane(p1, p2, p3);
Plane plane = construct_plane( p1, p2, p3 );
Line line = construct_line( l1, l2 );
const auto res = typename K::Intersect_3()(plane,line);
@ -2116,6 +2122,13 @@ namespace CommonKernelFunctors {
const Point* e_pt = std::get_if<Point>(&(*res));
CGAL_assertion(e_pt!=nullptr);
return *e_pt;
// Vector n = construct_cross_product_vector_3( vector(p1, p2), vector(p1, p3));
// CGAL_assertion(n!=CGAL::NULL_VECTOR);
// Vector l = vector(l1, l2);
// CGAL_assertion(l!=CGAL::NULL_VECTOR);
// FT t = scalar_product(n, vector(p1,l1)) / scalar_product(n, l);
// return l1 + t * l;
}
Point

View File

@ -265,6 +265,24 @@ bool close(PolygonMesh& pm, VertexPointMap vpm, typename Traits::Vector_3 plane_
#endif
template <class PolygonMesh,
class Visitor>
void naive_close(PolygonMesh& pm, Visitor& visitor)
{
using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
std::vector< halfedge_descriptor > border_cycles;
extract_boundary_cycles(pm, std::back_inserter(border_cycles));
std::vector< Bbox_3 > bboxes;
for (halfedge_descriptor h : border_cycles)
{
visitor.before_face_copy(boost::graph_traits<PolygonMesh>::null_face(), pm, pm);
Euler::fill_hole(h, pm);
visitor.after_face_copy(boost::graph_traits<PolygonMesh>::null_face(), pm, face(h, pm), pm);
}
}
template <class Plane_3,
class TriangleMesh,
@ -994,7 +1012,7 @@ bool clip(PolygonMesh& pm,
{
using parameters::choose_parameter;
using parameters::get_parameter;
using parameters::get_parameter_reference ;
using parameters::get_parameter_reference;
using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
@ -1026,11 +1044,12 @@ bool clip(PolygonMesh& pm,
const bool allow_self_intersections =
choose_parameter(get_parameter(np, internal_np::allow_self_intersections), false);
bool triangulate = !choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), false);
constexpr bool traits_supports_cdt2 = !internal::Has_member_Does_not_support_CDT2<GT>::value;
const bool used_for_kernel = choose_parameter(get_parameter(np, internal_np::used_for_kernel), false);
auto vos = get(dynamic_vertex_property_t<Oriented_side>(), pm);
auto ecm = get(dynamic_edge_property_t<bool>(), pm, false);
if (triangulate && !is_triangle_mesh(pm))
if (traits_supports_cdt2 && triangulate && !is_triangle_mesh(pm))
triangulate = false;
refine_with_plane(pm, plane, parameters::vertex_oriented_side_map(vos)
@ -1082,15 +1101,26 @@ bool clip(PolygonMesh& pm,
remove_connected_components(pm, ccs_to_remove, fcc);
if (clip_volume)
if constexpr (traits_supports_cdt2)
{
//TODO: add in the traits construct_orthogonal_vector
if (triangulate)
internal::close_and_triangulate<GT>(pm, vpm, plane.orthogonal_vector(), visitor);
else
if (!internal::close<GT>(pm, vpm, plane.orthogonal_vector(), visitor))
internal::close_and_triangulate<GT>(pm, vpm, plane.orthogonal_vector(), visitor);
if (clip_volume)
{
if (!used_for_kernel)
{
//TODO: add in the traits construct_orthogonal_vector
if (triangulate)
internal::close_and_triangulate<GT>(pm, vpm, plane.orthogonal_vector(), visitor);
else
if (!internal::close<GT>(pm, vpm, plane.orthogonal_vector(), visitor))
internal::close_and_triangulate<GT>(pm, vpm, plane.orthogonal_vector(), visitor);
}
else
internal::naive_close(pm, visitor);
}
}
else
if (clip_volume)
internal::naive_close(pm, visitor);
return true;
}

View File

@ -0,0 +1,302 @@
// Copyright (c) 2025 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Léo Valque
#ifndef CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H
#define CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H
#include <CGAL/license/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/boost/graph/named_params_helper.h>
namespace CGAL {
namespace Polygon_mesh_processing {
template <class PolygonMesh, class NamedParameters = parameters::Default_named_parameters>
void clip_convex(PolygonMesh& pm,
const typename GetGeomTraits<PolygonMesh, NamedParameters>::type::Plane_3& plane,
const NamedParameters& np = parameters::default_values())
{
using parameters::choose_parameter;
using parameters::get_parameter;
using parameters::get_parameter_reference;
using parameters::is_default_parameter;
// graph typedefs
using BGT = boost::graph_traits<PolygonMesh>;
using face_descriptor = typename BGT::face_descriptor;
using edge_descriptor = typename BGT::edge_descriptor;
using halfedge_descriptor = typename BGT::halfedge_descriptor;
using vertex_descriptor = typename BGT::vertex_descriptor;
// np typedefs
// using Default_ecm = Static_boolean_property_map<edge_descriptor, false>;
using Default_visitor = Default_cut_visitor<PolygonMesh>;
using Visitor_ref = typename internal_np::Lookup_named_param_def<internal_np::visitor_t, NamedParameters, Default_visitor>::reference;
using GT = typename GetGeomTraits<PolygonMesh, NamedParameters>::type;
GT traits = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
using FT = typename GT::FT;
using Point_3 = typename GT::Point_3;
Default_visitor default_visitor;
Visitor_ref visitor = choose_parameter(get_parameter_reference(np, internal_np::visitor), default_visitor);
constexpr bool has_visitor = !std::is_same_v<Default_visitor, std::remove_cv_t<std::remove_reference_t<Visitor_ref>>>;
auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(vertex_point, pm));
bool triangulate = !choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), true);
if (triangulate && !is_triangle_mesh(pm))
triangulate = false;
bool throw_on_self_intersection = choose_parameter(get_parameter(np, internal_np::use_compact_clipper), false);
if (throw_on_self_intersection && !is_triangle_mesh(pm))
throw_on_self_intersection = false;
// typedef typename internal_np::Lookup_named_param_def <
// internal_np::concurrency_tag_t,
// NamedParameters,
// Sequential_tag
// > ::type Concurrency_tag;
// constexpr bool parallel_execution = std::is_same_v<Parallel_tag, Concurrency_tag>;
auto oriented_side = traits.oriented_side_3_object();
auto intersection_point = traits.construct_plane_line_intersection_point_3_object();
auto sq = traits.compute_squared_distance_3_object();
// auto csq = traits.compare_squared_distance_3_object();
// auto vector_3 = traits.construct_vector_3_object();
// TODO: the default is not thread-safe for example for Polyhedron
using V_os_tag = dynamic_vertex_property_t<Oriented_side>;
static constexpr bool use_default_vosm =
is_default_parameter<NamedParameters, internal_np::vertex_oriented_side_map_t>::value;
using Vertex_oriented_side_map =
std::conditional_t<use_default_vosm,
typename boost::property_map<PolygonMesh, V_os_tag>::type,
typename internal_np::Get_param<typename NamedParameters::base,
internal_np::vertex_oriented_side_map_t>::type>;
Vertex_oriented_side_map vertex_os;
// ____________________ Find a crossing edge _____________________
vertex_descriptor src=*vertices(pm).begin();
FT sp_src = sq(plane, get(vpm, src)); // Not normalized distance
Sign direction_to_zero = sign(sp_src);
vertex_descriptor trg;
FT sp_trg;
bool is_crossing_edge=false;
if(direction_to_zero!=EQUAL){
do{
bool is_local_max=true;
for(auto v: vertices_around_target(src ,pm)){
sp_trg = sq(plane, get(vpm, v));
CGAL_assertion(sq(plane, get(vpm, v)) == sp_trg);
// TODO with EPICK, use compare_distance(plane, src, plane, trg) (But no possibility to memorize some computations for the next)
// Check if v in the direction to the plane
if(compare(sp_src, sp_trg)==direction_to_zero){
if(sign(sp_trg)!=direction_to_zero){
// Fund a crossing edge
trg = v;
is_crossing_edge=true;
} else {
// Continue from v
sp_src = sp_trg;
src = v;
}
is_local_max=false; // repeat with the new vertex
break;
}
}
// No intersection with the plane, kernel is either empty or full
if(is_local_max){
if(direction_to_zero==POSITIVE)
clear(pm); // The result is empty
return;
}
} while(!is_crossing_edge);
} else {
// src on the plane
trg=src;
sp_trg = sp_src;
}
if(sign(sp_trg)==EQUAL && direction_to_zero!=POSITIVE){
// Search a vertex around trg coming from positive side
bool no_positive_side = true;
for(auto v: vertices_around_target(trg ,pm)){
Oriented_side side_v = oriented_side(plane, get(vpm, v));
if(side_v==ON_POSITIVE_SIDE){
src = v;
sp_src = sq(plane, get(vpm, src));
no_positive_side = false;
break;
}
}
// Nothing to clip
if(no_positive_side)
return;
} else if(direction_to_zero==NEGATIVE){
// Orient the edge from negative to positive
std::swap(src, trg);
}
CGAL_assertion(oriented_side(plane, get(vpm, src)) == ON_POSITIVE_SIDE);
CGAL_assertion(oriented_side(plane, get(vpm, trg)) != ON_POSITIVE_SIDE);
// Cut the convex along the plane by marching along crossing edges starting from the previous edge
std::vector<halfedge_descriptor> boundaries;
std::vector<vertex_descriptor> boundaries_vertices;
halfedge_descriptor h = halfedge(src, trg, pm).first;
if(sign(sp_trg)!=EQUAL)
{
//split the first edge
auto pts = make_sorted_pair(get(vpm, src), get(vpm, trg));
typename GT::Point_3 ip = intersection_point(plane, pts.first, pts.second);
//visitor.before_edge_split(h, pm);
h = CGAL::Euler::split_edge(h, pm);
put(vpm, target(h, pm), ip);
}
vertex_descriptor v_start = target(h, pm);
halfedge_descriptor h_start=h;
do{
halfedge_descriptor h_previous = h;
CGAL_assertion(oriented_side(plane, get(vpm, source(h,pm)))==ON_POSITIVE_SIDE);
CGAL_assertion(oriented_side(plane, get(vpm, target(h,pm)))==ON_ORIENTED_BOUNDARY);
h = next(h, pm);
Oriented_side side_trg = oriented_side(plane, get(vpm, target(h,pm)));
if(side_trg != ON_NEGATIVE_SIDE){
// The face does not cross the plane
while(side_trg == ON_ORIENTED_BOUNDARY){
// The edge is along the plane, add it to boundaries
boundaries.emplace_back(h);
boundaries_vertices.emplace_back(target(h, pm));
set_halfedge(target(h, pm), h, pm);
h = next(h, pm);
side_trg=oriented_side(plane, get(vpm, target(h,pm)));
}
// continue on next face
h = opposite(h, pm);
continue;
}
// Search a crossing edge
h = next(h, pm);
side_trg=oriented_side(plane, get(vpm, target(h,pm)));
while(side_trg == ON_NEGATIVE_SIDE){
h = next(h,pm);
side_trg=oriented_side(plane, get(vpm, target(h,pm)));
}
if(side_trg != ON_ORIENTED_BOUNDARY){
// Split the edge
auto pts = make_sorted_pair(get(vpm, source(h,pm)), get(vpm, target(h,pm)));
typename GT::Point_3 ip = intersection_point(plane, pts.first, pts.second);
// visitor.before_edge_split(h, pm);
h = CGAL::Euler::split_edge(h, pm);
put(vpm, target(h, pm), ip);
// visitor.new_vertex_added(vid, target(h,pm), pm);
// visitor.edge_split(h, pm);
// visitor.after_edge_split();
}
// Split the face
halfedge_descriptor sh = CGAL::Euler::split_face(h_previous, h, pm);
boundaries.emplace_back(sh);
boundaries_vertices.emplace_back(target(sh, pm));
set_halfedge(target(sh, pm), sh, pm);
CGAL_assertion(target(sh, pm) == target(h, pm));
h = opposite(next(sh,pm), pm);
} while(target(h, pm)!=v_start || (boundaries.empty() && h!=h_start));
CGAL_assertion(is_valid_polygon_mesh(pm));
CGAL_assertion(!boundaries.empty());
// ________________________ Remove the negative side _________________________
std::set<vertex_descriptor> vertices_to_remove;
std::set<edge_descriptor> edges_to_remove;
std::set<face_descriptor> faces_to_remove;
std::sort(boundaries_vertices.begin(), boundaries_vertices.end());
// Go through to find all elements to delete
face_descriptor start_face(face(halfedge(src, pm), pm));
std::stack<face_descriptor> stack;
stack.push(start_face);
faces_to_remove.emplace(start_face);
while (!stack.empty()) {
face_descriptor f = stack.top();
stack.pop();
// Walk adjacent faces via halfedges
halfedge_descriptor h_start = halfedge(f, pm);
halfedge_descriptor h = h_start;
do {
if((std::find(boundaries.begin(), boundaries.end(), h)==boundaries.end()) &&
(std::find(boundaries.begin(), boundaries.end(), opposite(h,pm))==boundaries.end())){ // TODO avoid this linear operation
edges_to_remove.emplace(edge(h, pm));
if(!std::binary_search(boundaries_vertices.begin(), boundaries_vertices.end(), target(h, pm)))
vertices_to_remove.emplace(target(h, pm));
CGAL_assertion(oriented_side(plane, get(vpm, target(h, pm)))!=ON_NEGATIVE_SIDE);
face_descriptor fn = face(opposite(h, pm), pm);
if (faces_to_remove.find(fn)==faces_to_remove.end()) {
faces_to_remove.emplace(fn);
stack.push(fn);
}
}
h = next(h, pm);
} while (h != h_start);
}
for (vertex_descriptor v : vertices_to_remove){
remove_vertex(v, pm);
CGAL_assertion(oriented_side(plane, get(vpm, v))==ON_POSITIVE_SIDE);
}
for (edge_descriptor e : edges_to_remove)
remove_edge(e, pm);
for (face_descriptor f : faces_to_remove)
remove_face(f, pm);
// Reorder halfedges of the hole
for(size_t i=1; i<boundaries.size(); ++i)
set_next(boundaries[i-1], boundaries[i], pm);
set_next(boundaries.back(), boundaries[0], pm);
// Fill the hole
face_descriptor f=pm.add_face();
for(auto h: boundaries){
set_face(h, f, pm);
}
set_halfedge(f, boundaries[0], pm);
// std::ofstream("clip.off") << pm;
CGAL_assertion(is_valid_polygon_mesh(pm));
}
} } // CGAL::Polygon_mesh_processing
#endif // CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H

View File

@ -0,0 +1,663 @@
// Copyright (c) 2016-2025 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Sebastien Loriot
#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_KERNEL_H
#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_KERNEL_H
#include <CGAL/license/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/clip.h>
#include <CGAL/Polygon_mesh_processing/clip_convex.h>
#include <CGAL/Convex_hull_3/dual/halfspace_intersection_3.h>
#include <CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Homogeneous.h>
#include <CGAL/Exact_integer.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Cartesian_converter.h>
namespace CGAL {
namespace Polygon_mesh_processing {
namespace experimental {
template <class Kernel>
struct Three_point_cut_plane_traits
{
using FT = typename Kernel::FT;
using Plane_3 = std::array<typename Kernel::Point_3, 3>;
using Point_3 = typename Kernel::Point_3;
using Vector_3 = typename Kernel::Vector_3;
struct Does_not_support_CDT2{};
struct Oriented_side_3
{
Oriented_side operator()(const Plane_3& plane, const Point_3& p) const
{
return orientation(plane[0], plane[1], plane[2], p);
}
};
struct Construct_plane_line_intersection_point_3
{
Point_3 operator()(const Plane_3& plane, const Point_3& p, const Point_3& q)
{
typename Kernel::Construct_plane_line_intersection_point_3 construction;
return construction(plane[0], plane[1], plane[2], p, q);
}
};
struct Construct_orthogonal_vector_3{
Vector_3 operator()(const Plane_3& plane)
{
return typename Kernel::Plane_3(plane[0], plane[1], plane[2]).orthogonal_vector();
}
};
struct Compute_squared_distance_3
{
using Compute_scalar_product_3 = typename Kernel::Compute_scalar_product_3;
FT operator()(const Plane_3& plane, const Point_3& p)
{
typename Kernel::Plane_3 pl(plane[0], plane[1], plane[2]);
return Compute_scalar_product_3()(Vector_3(ORIGIN, p), pl.orthogonal_vector())+pl.d();
}
};
Oriented_side_3 oriented_side_3_object() const
{
return Oriented_side_3();
}
Construct_plane_line_intersection_point_3 construct_plane_line_intersection_point_3_object() const
{
return Construct_plane_line_intersection_point_3();
}
Construct_orthogonal_vector_3 construct_orthogonal_vector_3_object() const
{
return Construct_orthogonal_vector_3();
}
Compute_squared_distance_3 compute_squared_distance_3_object() const { return Compute_squared_distance_3(); }
#ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D
// for does self-intersect
using Segment_3 = typename Kernel::Segment_3;
using Triangle_3 = typename Kernel::Triangle_3;
using Construct_segment_3 = typename Kernel::Construct_segment_3;
using Construct_triangle_3 =typename Kernel::Construct_triangle_3;
using Do_intersect_3 = typename Kernel::Do_intersect_3;
Construct_segment_3 construct_segment_3_object() const { return Construct_segment_3(); }
Construct_triangle_3 construct_triangle_3_object() const { return Construct_triangle_3(); }
Do_intersect_3 do_intersect_3_object() const { return Do_intersect_3(); }
#endif
};
template <class TriangleMesh,
class NamedParameters = parameters::Default_named_parameters>
TriangleMesh
kernel(const TriangleMesh& pm,
const NamedParameters& np = parameters::default_values())
{
// TODO: bench if with EPECK we can directly use Kernel::Plane_3
// TODO: TriangleMesh as output is not correct since it is actually a PolygonMesh
using parameters::choose_parameter;
using parameters::get_parameter;
using GT = typename GetGeomTraits<TriangleMesh, NamedParameters>::type;
auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pm));
// GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
using Point_3 = typename GT::Point_3;
using parameters::choose_parameter;
using parameters::get_parameter;
//TODO: what do we do with a mesh that is not closed?
//TODO: what do we do if the input is not a triangle mesh?
if (vertices(pm).size() - edges(pm).size() + faces(pm).size() != 2)
return TriangleMesh();
CGAL::Bbox_3 bb3 = bbox(pm, np);
TriangleMesh kernel;
CGAL::make_hexahedron(Point_3(bb3.xmax(),bb3.ymin(),bb3.zmin()), Point_3(bb3.xmax(),bb3.ymax(),bb3.zmin()), Point_3(bb3.xmin(),bb3.ymax(),bb3.zmin()), Point_3(bb3.xmin(),bb3.ymin(),bb3.zmin()),
Point_3(bb3.xmin(),bb3.ymin(),bb3.zmax()), Point_3(bb3.xmax(),bb3.ymin(),bb3.zmax()), Point_3(bb3.xmax(),bb3.ymax(),bb3.zmax()), Point_3(bb3.xmin(),bb3.ymax(),bb3.zmax()),
kernel);
Three_point_cut_plane_traits<GT> kgt;
#ifdef CGAL_USE_CONCAVE_FACE_FIRST
auto compute_dihedral_angle=[](const TriangleMesh& mesh, const halfedge_descriptor h) {
// Calcul de l'angle dihédral entre les deux faces à partir des quatre points
return to_double(approximate_dihedral_angle(mesh.point(mesh.source(h)), mesh.point(mesh.target(h)),
mesh.point(mesh.target(mesh.next(h))), mesh.point(mesh.target(mesh.next(mesh.opposite(h))))));
};
std::vector<std::pair<face_descriptor, double>> faces_with_angles;
for(face_descriptor f : faces(pm)){
double max_dihedral_angle = 0.0;
auto h = pm.halfedge(f);
auto hf_circ = pm.halfedges_around_face(h);
for (auto he : hf_circ) {
if (!pm.is_border(he)) {
max_dihedral_angle = (std::max)(max_dihedral_angle,compute_dihedral_angle(pm, he));
}
faces_with_angles.push_back({f, max_dihedral_angle});
}
}
std::sort(faces_with_angles.begin(), faces_with_angles.end(),
[](const std::pair<face_descriptor, double>& a, const std::pair<face_descriptor, double>& b) {
return a.second > b.second;
});
for(auto pair : faces_with_angles)
{
auto h = halfedge(pair.first, pm);
auto plane = make_array( get(vpm,source(h, pm)),
get(vpm,target(h, pm)),
get(vpm,target(next(h, pm), pm)) );
#else
for (auto f : faces(pm))
{
auto h = halfedge(f, pm);
auto plane = make_array( get(vpm,source(h, pm)),
get(vpm,target(h, pm)),
get(vpm,target(next(h, pm), pm)) );
#endif
#ifdef CGAL_USE_OPTI_WITH_BBOX
auto pred = kgt.oriented_side_3_object();
bb3 = bbox(kernel);
std::array<Point_3, 8> corners = CGAL::make_array(Point_3(bb3.xmax(),bb3.ymin(),bb3.zmin()),
Point_3(bb3.xmax(),bb3.ymax(),bb3.zmin()),
Point_3(bb3.xmin(),bb3.ymax(),bb3.zmin()),
Point_3(bb3.xmin(),bb3.ymin(),bb3.zmin()),
Point_3(bb3.xmin(),bb3.ymin(),bb3.zmax()),
Point_3(bb3.xmax(),bb3.ymin(),bb3.zmax()),
Point_3(bb3.xmax(),bb3.ymax(),bb3.zmax()),
Point_3(bb3.xmin(),bb3.ymax(),bb3.zmax()));
int i=0;
auto first_ori=pred(plane, corners[i]);
while(++i!=8 && first_ori==ON_ORIENTED_BOUNDARY)
first_ori=pred(plane, corners[i]);
if (i==8) continue;
bool all_the_same=true;
for (;i<8;++i)
{
auto other_ori=pred(plane, corners[i]);
if (other_ori!=ON_ORIENTED_BOUNDARY && other_ori!=first_ori)
{
all_the_same=false;
break;
}
}
if (all_the_same)
{
if (first_ori==ON_NEGATIVE_SIDE) continue;
else
{
return TriangleMesh();
}
}
#endif
clip_convex(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(kgt).do_not_triangulate_faces(true).used_for_kernel(true));
if (is_empty(kernel)) break;
}
return kernel;
};
template <class TriangleMesh,
class NamedParameters = parameters::Default_named_parameters>
TriangleMesh
kernel_using_chull(const TriangleMesh& pm,
const NamedParameters& np = parameters::default_values())
{
using parameters::choose_parameter;
using parameters::get_parameter;
using GT = typename GetGeomTraits<TriangleMesh, NamedParameters>::type;
auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pm));
// GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
using parameters::choose_parameter;
using parameters::get_parameter;
std::vector<typename GT::Plane_3> planes;
planes.reserve(faces(pm).size());
for (auto f : faces(pm))
{
auto h = halfedge(f, pm);
planes.emplace_back( get(vpm,source(h, pm)),
get(vpm,target(h, pm)),
get(vpm,target(next(h, pm), pm)) );
}
auto origin_opt = CGAL::halfspace_intersection_interior_point_3(planes.begin(), planes.end());
TriangleMesh kernel;
if (origin_opt.has_value())
CGAL::halfspace_intersection_3(planes.begin(),
planes.end(),
kernel,
origin_opt);
return kernel;
}
template <class TriangleMesh,
class NamedParameters = parameters::Default_named_parameters>
TriangleMesh
kernel_using_chull_and_constructions(const TriangleMesh& pm,
const NamedParameters& np = parameters::default_values())
{
using parameters::choose_parameter;
using parameters::get_parameter;
using GT = typename GetGeomTraits<TriangleMesh, NamedParameters>::type;
auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pm));
// GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
using parameters::choose_parameter;
using parameters::get_parameter;
std::vector<typename GT::Plane_3> planes;
planes.reserve(faces(pm).size());
for (auto f : faces(pm))
{
auto h = halfedge(f, pm);
planes.emplace_back( get(vpm,source(h, pm)),
get(vpm,target(h, pm)),
get(vpm,target(next(h, pm), pm)) );
}
auto origin_opt = CGAL::halfspace_intersection_interior_point_3(planes.begin(), planes.end());
TriangleMesh kernel;
if (origin_opt.has_value())
CGAL::halfspace_intersection_with_constructions_3(planes.begin(),
planes.end(),
kernel,
origin_opt);
return kernel;
}
//__________________________________________________________________________________________________________________________________
using int256 = boost::multiprecision::int256_t;
using int512 = boost::multiprecision::int512_t;
using intexact = Exact_integer;
template <class Kernel>
struct Plane_based_traits
{
using FT = typename Kernel::FT;
using RT = typename Kernel::RT;
using Geometric_point_3 = typename Kernel::Point_3;
using Vector_3 = typename Kernel::Vector_3;
using plane_descriptor = std::size_t;
using Plane_3 = std::pair<typename Kernel::Plane_3, plane_descriptor>;
private:
using plane_range_pointer = std::shared_ptr<std::vector<Plane_3>>;
plane_range_pointer m_planes;
public:
struct Point_3 : public Geometric_point_3{
std::array<plane_descriptor, 3> supports;
mutable std::set<plane_descriptor> other_coplanar;
using Base = Geometric_point_3;
Point_3(){}
Point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c, const std::vector<Plane_3> &planes):
Base(*CGAL::Intersections::internal::intersection_point(planes[a].first, planes[b].first, planes[c].first, Kernel())),
supports({a,b,c}){
}
Point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p): Base(p), supports({a,b,c}){}
};
struct Construct_point_3{
plane_range_pointer m_planes;
Construct_point_3(plane_range_pointer planes) : m_planes(planes){}
Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c){
const std::vector<Plane_3> &planes = *m_planes;
auto res = intersection(planes[a].first, planes[b].first, planes[c].first);
return Point_3(a, b, c, std::get<Geometric_point_3>(*res));
}
Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p){
return Point_3(a, b, c, p);
}
};
struct Does_not_support_CDT2{};
struct Oriented_side_3
{
Oriented_side operator()(const Plane_3& plane, const Point_3& p) const
{
if((p.supports[0]==plane.second) || (p.supports[1]==plane.second) || (p.supports[2]==plane.second))
return COPLANAR;
Oriented_side ori = plane.first.oriented_side(p);
if(ori==COPLANAR)
p.other_coplanar.emplace(plane.second);
return ori;
}
};
struct Construct_plane_line_intersection_point_3
{
plane_range_pointer m_planes;
Construct_plane_line_intersection_point_3(plane_range_pointer planes) : m_planes(planes){}
Point_3 operator()(const Plane_3& plane, const Point_3& p, const Point_3& q)
{
Construct_point_3 point_3(m_planes);
const std::vector<Plane_3> &planes=*m_planes;
auto get_common_supports=[&](const Point_3& p, const Point_3 &q){
plane_descriptor first=-1;
plane_descriptor second=-1;
for(short int i=0; i!=3; ++i)
for(short int j=0; j!=3; ++j)
if(p.supports[i]==q.supports[j]){
if(first==-1){
first=p.supports[i];
break;
} else {
second=p.supports[i];
return std::make_pair(first, second);
}
}
for(plane_descriptor pd: p.other_coplanar)
for(short int j=0; j!=3; ++j)
if(pd==q.supports[j]){
if(first==-1){
first=pd;
break;
} else if(planes[first].first!=planes[pd].first){
second=pd;
return std::make_pair(first, second);
}
}
for(plane_descriptor pd: q.other_coplanar)
for(short int j=0; j!=3; ++j)
if(pd==p.supports[j]){
if(first==-1){
first=pd;
break;
} else if(planes[first].first!=planes[pd].first){
second=pd;
return std::make_pair(first, second);
}
}
for(plane_descriptor pd: p.other_coplanar)
for(plane_descriptor qd: q.other_coplanar)
if(pd==qd){
if(first==-1){
first=pd;
break;
} else if(planes[first].first!=planes[pd].first){
second=pd;
return std::make_pair(first, second);
}
}
// The two points do not shair a common support
std::cout << p.supports[0] << " " << p.supports[1] << " " << p.supports[2] << std::endl;
std::cout << q.supports[0] << " " << q.supports[1] << " " << q.supports[2] << std::endl;
CGAL_assertion_code(std::cout << "The two points do not share two common supporting planes" << std::endl;)
CGAL_assertion(0);
return std::make_pair(first, second);
};
std::pair<plane_descriptor, plane_descriptor> line_supports=get_common_supports(p, q);
assert((planes[line_supports.first].first != planes[line_supports.second].first) &&
(planes[line_supports.first].first != planes[line_supports.second].first.opposite()));
Point_3 res=point_3(plane.second, line_supports.first, line_supports.second);
namespace mp = boost::multiprecision;
CGAL_assertion_code(int256 max2E195=mp::pow(int256(2),195);)
CGAL_assertion_code(int256 max2E169=mp::pow(int256(2),169);)
// CGAL_assertion((mp::abs(p.hx())<=max2E195) && (mp::abs(p.hy())<=max2E195) && (mp::abs(p.hz())<=max2E195) && (mp::abs(p.hw())<=max2E169));
return res;
}
};
struct Construct_orthogonal_vector_3
{
Vector_3 operator()(const Plane_3& plane, const Point_3& p) const
{
if((p.supports[0]==plane.second) || (p.supports[1]==plane.second) || (p.supports[2]==plane.second))
return COPLANAR;
Oriented_side ori = plane.first.oriented_side(p);
if(ori==COPLANAR)
p.other_coplanar.emplace(plane.second);
return ori;
}
};
struct Compute_squared_distance_3
{
using Compute_scalar_product_3 = typename Kernel::Compute_scalar_product_3;
FT operator()(const Plane_3& plane, const Point_3& p)
{
return Compute_scalar_product_3()(Vector_3(ORIGIN, p), plane.first.orthogonal_vector())+plane.first.d();
}
};
template<class Mesh>
Plane_based_traits(const Mesh &m){
namespace mp = boost::multiprecision;
using face_descriptor = typename boost::graph_traits<Mesh>::face_descriptor;
using vertex_descriptor = typename boost::graph_traits<Mesh>::vertex_descriptor;
auto to_int=[](const typename Mesh::Point &p){
return Geometric_point_3(int(p.x()),int(p.y()),int(p.z()));
};
auto to_int_plane=[&](face_descriptor f){
auto pmap=boost::make_function_property_map<vertex_descriptor>([&](vertex_descriptor v){
return to_int(m.point(v));
});
auto pl = compute_face_normal(f, m, parameters::vertex_point_map(pmap));
return typename Kernel::Vector_3(pl.hx(),pl.hy(),pl.hz());
};
m_planes = std::make_shared<std::vector<Plane_3>>();
m_planes->reserve(faces(m).size());
// The ptr need to be initialized to get the functor
Construct_plane_3 plane_3 = construct_plane_3_object();
for(face_descriptor f : faces(m))
plane_3(typename Kernel::Plane_3(to_int(m.point(m.target(m.halfedge(f)))),
to_int_plane(f)));
}
struct Construct_plane_3{
plane_range_pointer m_planes;
Construct_plane_3(plane_range_pointer planes) : m_planes(planes){}
Plane_3 operator()(const typename Kernel::Plane_3 &pl){
namespace mp = boost::multiprecision;
m_planes->emplace_back(pl, m_planes->size());
CGAL_assertion_code(int256 max2E55=mp::pow(int256(2),55);)
CGAL_assertion_code(int256 max2E82=mp::pow(int256(2),82);)
// CGAL_assertion((mp::abs(pl.a())<=max2E55) && (mp::abs(pl.b())<=max2E55) && (mp::abs(pl.c())<=max2E82) && (mp::abs(pl.d())<=max2E82));
return m_planes->back();
}
Plane_3 operator()(const RT &a, const RT &b, const RT &c, const RT &d){
return (*this)(typename Kernel::Plane_3(a,b,c,d));
}
};
Oriented_side_3 oriented_side_3_object() const { return Oriented_side_3(); }
Construct_plane_3 construct_plane_3_object() { return Construct_plane_3(m_planes); }
Construct_point_3 construct_point_3_object() { return Construct_point_3(m_planes); }
Construct_orthogonal_vector_3 construct_orthogonal_vector_3_object() const { return Construct_orthogonal_vector_3(); }
Compute_squared_distance_3 compute_squared_distance_3_object() const { return Compute_squared_distance_3(); }
Construct_plane_line_intersection_point_3 construct_plane_line_intersection_point_3_object() const {
return Construct_plane_line_intersection_point_3(m_planes);
}
const std::vector<Plane_3>& planes(){
return *m_planes;
}
#ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D
// for does self-intersect
using Segment_3 = typename Kernel::Segment_3;
using Triangle_3 = typename Kernel::Triangle_3;
using Construct_segment_3 = typename Kernel::Construct_segment_3;
using Construct_triangle_3 =typename Kernel::Construct_triangle_3;
using Do_intersect_3 = typename Kernel::Do_intersect_3;
Construct_segment_3 construct_segment_3_object() const { return Construct_segment_3(); }
Construct_triangle_3 construct_triangle_3_object() const { return Construct_triangle_3(); }
Do_intersect_3 do_intersect_3_object() const { return Do_intersect_3(); }
#endif
};
template <class TriangleMesh,
class NamedParameters = parameters::Default_named_parameters>
Surface_mesh<Plane_based_traits<Homogeneous<int256>>::Point_3>
trettner_kernel(const TriangleMesh& pm,
const NamedParameters& np = parameters::default_values())
{
using parameters::choose_parameter;
using parameters::get_parameter;
using GT = Plane_based_traits<Homogeneous<int256>>;
using Point_3 = typename GT::Point_3;
using Plane_3 = typename GT::Plane_3;
using Construct_plane_3 = GT::Construct_plane_3;
using Construct_point_3 = GT::Construct_point_3;
using InternMesh = Surface_mesh<Point_3>;
GT gt(pm);
const std::vector<Plane_3> &planes=gt.planes();
Construct_plane_3 plane_3 = gt.construct_plane_3_object();
Construct_point_3 point_3 = gt.construct_point_3_object();
if (vertices(pm).size() - edges(pm).size() + faces(pm).size() != 2)
return Surface_mesh<Point_3>();
CGAL::Bbox_3 bb3 = bbox(pm, np);
InternMesh kernel;
Plane_3 xl=plane_3(1,0,0,int(-bb3.xmin()));
Plane_3 yl=plane_3(0,1,0,int(-bb3.ymin()));
Plane_3 zl=plane_3(0,0,1,int(-bb3.zmin()));
Plane_3 xr=plane_3(1,0,0,int(-bb3.xmax()));
Plane_3 yr=plane_3(0,1,0,int(-bb3.ymax()));
Plane_3 zr=plane_3(0,0,1,int(-bb3.zmax()));
CGAL::make_hexahedron(point_3(xl.second, yl.second, zl.second),
point_3(xl.second, yl.second, zr.second),
point_3(xl.second, yr.second, zr.second),
point_3(xl.second, yr.second, zl.second),
point_3(xr.second, yr.second, zl.second),
point_3(xr.second, yl.second, zl.second),
point_3(xr.second, yl.second, zr.second),
point_3(xr.second, yr.second, zr.second),
kernel);
for (auto plane: planes)
{
clip_convex(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true));
if (is_empty(kernel)) break;
}
return kernel;
}
template <class TriangleMesh,
class NamedParameters = parameters::Default_named_parameters>
Surface_mesh<Plane_based_traits<Exact_predicates_exact_constructions_kernel>::Point_3>
trettner_epeck_kernel(const TriangleMesh& pm,
const NamedParameters& np = parameters::default_values())
{
using parameters::choose_parameter;
using parameters::get_parameter;
// using GT = Plane_based_traits<Homogeneous<int256>>;
using GT = Plane_based_traits<Exact_predicates_exact_constructions_kernel>;
using Point_3 = typename GT::Point_3;
using Plane_3 = typename GT::Plane_3;
using Construct_plane_3 = GT::Construct_plane_3;
using Construct_point_3 = GT::Construct_point_3;
using InternMesh = Surface_mesh<Point_3>;
GT gt(pm);
const std::vector<Plane_3> &planes=gt.planes();
Construct_plane_3 plane_3 = gt.construct_plane_3_object();
Construct_point_3 point_3 = gt.construct_point_3_object();
if (vertices(pm).size() - edges(pm).size() + faces(pm).size() != 2)
return Surface_mesh<Point_3>();
CGAL::Bbox_3 bb3 = bbox(pm, np);
InternMesh kernel;
Plane_3 xl=plane_3(1,0,0,int(-bb3.xmin()));
Plane_3 yl=plane_3(0,1,0,int(-bb3.ymin()));
Plane_3 zl=plane_3(0,0,1,int(-bb3.zmin()));
Plane_3 xr=plane_3(1,0,0,int(-bb3.xmax()));
Plane_3 yr=plane_3(0,1,0,int(-bb3.ymax()));
Plane_3 zr=plane_3(0,0,1,int(-bb3.zmax()));
CGAL::make_hexahedron(point_3(xl.second, yl.second, zl.second),
point_3(xl.second, yl.second, zr.second),
point_3(xl.second, yr.second, zr.second),
point_3(xl.second, yr.second, zl.second),
point_3(xr.second, yr.second, zl.second),
point_3(xr.second, yl.second, zl.second),
point_3(xr.second, yl.second, zr.second),
point_3(xr.second, yr.second, zr.second),
kernel);
for (auto plane: planes)
{
clip_convex(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true));
if (is_empty(kernel)) break;
}
return kernel;
}
} } } // end of CGAL::Polygon_mesh_processing::experimental
#endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_KERNEL_H

View File

@ -24,10 +24,18 @@
#ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#endif
#include <boost/mpl/has_xxx.hpp>
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal
{
BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(Has_member_Does_not_support_CDT2,
Does_not_support_CDT2,
false)
} // internal namespace
#ifndef DOXYGEN_RUNNING
template <class PolygonMesh>
struct Default_cut_visitor
@ -74,6 +82,8 @@ struct Orthogonal_cut_plane_traits
using Plane_3 = std::pair<int, FT>;
using Point_3 = typename Kernel::Point_3;
struct Does_not_support_CDT2{};
struct Oriented_side_3
{
Oriented_side operator()(const Plane_3& plane, const Point_3& p) const

View File

@ -74,6 +74,7 @@ create_single_source_cgal_program("test_corefinement_nm_bo.cpp")
create_single_source_cgal_program("test_corefinement_cavities.cpp")
create_single_source_cgal_program("issue_8730.cpp")
create_single_source_cgal_program("issue_7164.cpp")
create_single_source_cgal_program("test_mesh_kernel.cpp")
# create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp")
find_package(Eigen3 QUIET)

View File

@ -0,0 +1,287 @@
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/clip.h>
#include <CGAL/Polygon_mesh_processing/internal/kernel.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Real_timer.h>
using K=CGAL::Exact_predicates_inexact_constructions_kernel;
using EK=CGAL::Exact_predicates_exact_constructions_kernel;
using Mesh=CGAL::Surface_mesh<K::Point_3>;
using EMesh=CGAL::Surface_mesh<EK::Point_3>;
namespace PMP=CGAL::Polygon_mesh_processing;
template<class Mesh>
void round_mesh(Mesh &m){
auto exp = [](const double v)
{
int n;
frexp(v, &n);
return n;
};
auto pow_2 = [](const int k)
{
return (k>=0)?std::pow(2,k):1./std::pow(2,-k);
};
CGAL::Bbox_3 bb;
for(auto vh : m.vertices()){
bb += m.point(vh).bbox();
}
size_t grid_size=26;
std::array<double, 3> max_abs{(std::max)(-bb.xmin(), bb.xmax()),
(std::max)(-bb.ymin(), bb.ymax()),
(std::max)(-bb.zmin(), bb.zmax())};
// Compute scale so that the exponent of max absolute value are 52-1.
std::array<double, 3> scale{pow_2(grid_size - exp(max_abs[0]) - 1),
pow_2(grid_size - exp(max_abs[1]) - 1),
pow_2(grid_size - exp(max_abs[2]) - 1)};
auto round=[&](typename Mesh::Point &p){
return typename Mesh::Point(std::ceil((CGAL::to_double(p.x()) * scale[0]) - 0.5) / scale[0],
std::ceil((CGAL::to_double(p.y()) * scale[1]) - 0.5) / scale[1],
std::ceil((CGAL::to_double(p.z()) * scale[2]) - 0.5) / scale[2]);
};
for(auto vh : m.vertices()){
auto &p = m.point(vh);
p = round(p);
}
}
template<class Mesh>
void to_integer_mesh(Mesh &m){
auto exp = [](const double v)
{
int n;
frexp(v, &n);
return n;
};
auto pow_2 = [](const int k)
{
return (k>=0)?std::pow(2,k):1./std::pow(2,-k);
};
CGAL::Bbox_3 bb;
for(auto vh : m.vertices()){
bb += m.point(vh).bbox();
}
size_t grid_size=26;
std::array<double, 3> max_abs{(std::max)(-bb.xmin(), bb.xmax()),
(std::max)(-bb.ymin(), bb.ymax()),
(std::max)(-bb.zmin(), bb.zmax())};
// Compute scale so that the exponent of max absolute value are 52-1.
std::array<double, 3> scale{pow_2(grid_size - exp(max_abs[0]) - 1),
pow_2(grid_size - exp(max_abs[1]) - 1),
pow_2(grid_size - exp(max_abs[2]) - 1)};
// auto round=[&](Point_3 &p){
auto round=[&](typename Mesh::Point &p){
return typename Mesh::Point(std::ceil((CGAL::to_double(p.x()) * scale[0]) - 0.5),
std::ceil((CGAL::to_double(p.y()) * scale[1]) - 0.5),
std::ceil((CGAL::to_double(p.z()) * scale[2]) - 0.5));
};
for(auto vh : m.vertices()){
auto &p = m.point(vh);
p = round(p);
}
}
void test_traits()
{
Mesh m;
std::ifstream(CGAL::data_file_path("meshes/elephant.off")) >> m;
std::pair p(0,0.05);
using Traits = PMP::Orthogonal_cut_plane_traits<K>;
PMP::clip(m, p, CGAL::parameters::geom_traits(Traits()));
std::ofstream("clipped.off") << m;
}
void elementary_test_kernel()
{
using Point_3 = typename EK::Point_3;
CGAL::Real_timer timer;
EMesh m;
// make_hexahedron(Point_3(0, 1, 0), Point_3(0, 1, 1), Point_3(1, 0, 1), Point_3(1, 0, 0),
// Point_3(0.75, 0.25, 0.25), Point_3(0.25, 0.75, 0.25), Point_3(0.25, 0.75, 0.75), Point_3(0.75, 0.25, 0.75),
// m);
// make_tetrahedron(Point_3(1, 0, 0), Point_3(0, 1, 0), Point_3(0, 0, 1), Point_3(0.5, 0.5, 0.5),
// m);
// make_tetrahedron(Point_3(1, 0, 0), Point_3(1, 1, 0), Point_3(0, 1, 0), Point_3(0, 0, 0),
// m);
make_hexahedron(Point_3(1,0,0), Point_3(1,1,0), Point_3(0,1,0), Point_3(0,0,0),
Point_3(0,0,1), Point_3(1,0,1), Point_3(1,1,1), Point_3(0,1,1),
m);
timer.start();
EMesh kernel = PMP::experimental::kernel(m);
timer.stop();
std::ofstream("kernel.off") << kernel;
std::cout << "test_kernel done in " << timer.time() << "\n";
}
void test_kernel(std::string fname)
{
CGAL::Real_timer timer;
Mesh m;
if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m))
{
std::cerr << "ERROR: cannot read " << fname << "\n";
exit(1);
}
timer.start();
Mesh kernel = PMP::experimental::kernel(m);
timer.stop();
std::ofstream("kernel.off") << kernel;
std::cout << "test_kernel done in " << timer.time() << "\n";
}
void test_kernel_with_rounding(std::string fname)
{
CGAL::Real_timer timer;
Mesh m;
if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m))
{
std::cerr << "ERROR: cannot read " << fname << "\n";
exit(1);
}
round_mesh(m);
timer.start();
Mesh kernel = PMP::experimental::kernel(m);
timer.stop();
std::ofstream("kernel.off") << kernel;
std::cout << "test_kernel done in " << timer.time() << "\n";
}
void test_exact_kernel(std::string fname)
{
CGAL::Real_timer timer;
EMesh m;
if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m))
{
std::cerr << "ERROR: cannot read " << fname << "\n";
exit(1);
}
timer.start();
EMesh kernel = PMP::experimental::kernel(m);
timer.stop();
std::ofstream("ekernel.off") << kernel;
std::cout << "test_exact_kernel done in " << timer.time() << "\n";
}
void test_exact_kernel_with_rounding(std::string fname)
{
CGAL::Real_timer timer;
EMesh m;
if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m))
{
std::cerr << "ERROR: cannot read " << fname << "\n";
exit(1);
}
// round_mesh(m);
to_integer_mesh(m);
timer.start();
EMesh kernel = PMP::experimental::kernel(m);
timer.stop();
std::ofstream("ekernel.off") << kernel;
std::cout << "test_exact_kernel done in " << timer.time() << "\n";
}
void test_kernel_with_chull(std::string fname)
{
CGAL::Real_timer timer;
Mesh m;
if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m))
{
std::cerr << "ERROR: cannot read " << fname << "\n";
exit(1);
}
timer.start();
Mesh kernel = PMP::experimental::kernel_using_chull(m);
timer.stop();
std::ofstream("kernel_with_chull.off") << kernel;
std::cout << "test_kernel_with_chull done in " << timer.time() << "\n";
}
void test_kernel_with_chull_and_constructions(std::string fname)
{
CGAL::Real_timer timer;
Mesh m;
if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m))
{
std::cerr << "ERROR: cannot read " << fname << "\n";
exit(1);
}
timer.start();
Mesh kernel = PMP::experimental::kernel_using_chull_and_constructions(m);
timer.stop();
std::ofstream("kernel_with_chull_and_constructions.off") << kernel;
std::cout << "test_kernel_with_chull_and_constructions done in " << timer.time() << "\n";
}
void test_trettner_kernel(std::string fname)
{
CGAL::Real_timer timer;
Mesh m;
if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m))
{
std::cerr << "ERROR: cannot read " << fname << "\n";
exit(1);
}
to_integer_mesh(m);
timer.start();
auto kernel = PMP::experimental::trettner_kernel(m);
timer.stop();
std::ofstream("tkernel.off") << kernel;
std::cout << "test_trettner_kernel done in " << timer.time() << "\n";
}
void test_trettner_epeck_kernel(std::string fname)
{
CGAL::Real_timer timer;
Mesh m;
if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m))
{
std::cerr << "ERROR: cannot read " << fname << "\n";
exit(1);
}
to_integer_mesh(m);
timer.start();
auto kernel = PMP::experimental::trettner_epeck_kernel(m);
timer.stop();
std::ofstream("tkernel.off") << kernel;
std::cout << "test_trettner_kernel with epeck done in " << timer.time() << "\n";
}
int main(int argc, char** argv)
{
// test_traits();
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
// elementary_test_kernel();
// test_kernel(filename);
// test_kernel_with_rounding(filename);
// test_exact_kernel(filename);
// test_exact_kernel_with_rounding(filename);
// test_trettner_kernel(filename);
test_trettner_epeck_kernel(filename);
// test_kernel_with_chull(filename);
// test_kernel_with_chull_and_constructions(filename);
}

View File

@ -200,6 +200,7 @@ CGAL_add_named_parameter(force_filtering_t, force_filtering, force_filtering)
//internal
CGAL_add_named_parameter(weight_calculator_t, weight_calculator, weight_calculator)
CGAL_add_named_parameter(use_bool_op_to_clip_surface_t, use_bool_op_to_clip_surface, use_bool_op_to_clip_surface)
CGAL_add_named_parameter(used_for_kernel_t, used_for_kernel, used_for_kernel)
// List of named parameters used in the Point Set Processing package
CGAL_add_named_parameter(query_point_t, query_point_map, query_point_map)