update the implementation of clip plugin to use new corefinement PMP code

This commit is contained in:
Sébastien Loriot 2017-01-10 15:07:14 +01:00
parent c9bbc4e8d1
commit 59001acd92
6 changed files with 466 additions and 594 deletions

View File

@ -1,437 +0,0 @@
// Copyright (c) 2015 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
//
// Author(s) : Sebastien Loriot
#ifndef CGAL_INTERNAL_POLYHEDRON_PLANE_CLIPPING_3_H
#define CGAL_INTERNAL_POLYHEDRON_PLANE_CLIPPING_3_H
#include <CGAL/corefinement_operations.h>
#include <CGAL/iterator.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/array.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_2_projection_traits_3.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/convex_hull_3.h>
namespace CGAL{
namespace corefinement{
namespace internal{
template <class HDS,class T>
class Builder_from_T_2 : public CGAL::Modifier_base<HDS> {
typedef std::map<typename T::Vertex_handle,unsigned> Vertex_map;
const T& t;
template <class Builder>
static unsigned get_vertex_index( Vertex_map& vertex_map,
typename T::Vertex_handle vh,
Builder& builder,
unsigned& vindex)
{
std::pair<typename Vertex_map::iterator,bool>
res=vertex_map.insert(std::make_pair(vh,vindex));
if (res.second){
builder.add_vertex(vh->point());
++vindex;
}
return res.first->second;
}
public:
Builder_from_T_2(const T& t_):t(t_)
{
CGAL_assertion(t.dimension()==2);
}
void operator()( HDS& hds) {
// Postcondition: `hds' is a valid polyhedral surface.
CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
Vertex_map vertex_map;
//start the surface
B.begin_surface( t.number_of_vertices(), t.number_of_faces());
unsigned vindex=0;
for (typename T::Finite_faces_iterator it=t.finite_faces_begin();
it!=t.finite_faces_end();++it)
{
unsigned i0=get_vertex_index(vertex_map,it->vertex(0),B,vindex);
unsigned i1=get_vertex_index(vertex_map,it->vertex(1),B,vindex);
unsigned i2=get_vertex_index(vertex_map,it->vertex(2),B,vindex);
B.begin_facet();
B.add_vertex_to_facet( i0 );
B.add_vertex_to_facet( i1 );
B.add_vertex_to_facet( i2 );
B.end_facet();
}
B.end_surface();
}
};
} // end of namespace internal
template <class Polyhedron, class Plane_3>
Polyhedron clip_to_bbox(const Bbox_3& bbox, const Plane_3& plane)
{
typedef typename Polyhedron::Traits::Kernel Kernel;
typedef typename Kernel::Point_3 Point_3;
typedef typename Kernel::Segment_3 Segment_3;
cpp11::array<typename Kernel::Point_3,8> corners= {{
Point_3(bbox.xmin(),bbox.ymin(),bbox.zmin()),
Point_3(bbox.xmin(),bbox.ymax(),bbox.zmin()),
Point_3(bbox.xmax(),bbox.ymax(),bbox.zmin()),
Point_3(bbox.xmax(),bbox.ymin(),bbox.zmin()),
Point_3(bbox.xmin(),bbox.ymin(),bbox.zmax()),
Point_3(bbox.xmin(),bbox.ymax(),bbox.zmax()),
Point_3(bbox.xmax(),bbox.ymax(),bbox.zmax()),
Point_3(bbox.xmax(),bbox.ymin(),bbox.zmax())
}};
cpp11::array<CGAL::Oriented_side,8> orientations = {{
plane.oriented_side(corners[0]),
plane.oriented_side(corners[1]),
plane.oriented_side(corners[2]),
plane.oriented_side(corners[3]),
plane.oriented_side(corners[4]),
plane.oriented_side(corners[5]),
plane.oriented_side(corners[6]),
plane.oriented_side(corners[7])
}};
std::vector<Point_3> intersection_points;
// first look for intersections at corners
for (int i=0; i<8; ++i)
if (orientations[i]==ON_ORIENTED_BOUNDARY)
intersection_points.push_back(corners[i]);
// second look for intersections on edges
cpp11::array<int,24> edge_indices = {{ // 2 *12 edges
0,1, 1,2, 2,3, 3,0, // bottom face edges
4,5, 5,6, 6,7, 7,4, // top face edges
0,4, 1,5, 2,6, 3,7
}};
for (int i=0; i<12; ++i)
{
int i1=edge_indices[2*i], i2=edge_indices[2*i+1];
if (orientations[i1]==ON_ORIENTED_BOUNDARY) continue;
if (orientations[i2]==ON_ORIENTED_BOUNDARY) continue;
if (orientations[i1]!=orientations[i2])
intersection_points.push_back(
boost::get<Point_3>(
*CGAL::intersection(plane, Segment_3(corners[i1], corners[i2]) )
)
);
}
Polyhedron P;
//if less that 3 points there will be nothing to clipping.
if (intersection_points.size()<3) return P;
//triangulate the set of intersection points (I know it's overkill)
typedef CGAL::Triangulation_2_projection_traits_3<Kernel> P_traits;
typedef CGAL::Delaunay_triangulation_2<P_traits> DT;
DT dt(P_traits(plane.orthogonal_vector()));
dt.insert(intersection_points.begin(),
intersection_points.end());
// tangency with the bbox -> no intersection
if (dt.dimension()!=2) return P;
//now create the polyhedron from the triangulation
internal::Builder_from_T_2< typename Polyhedron::HalfedgeDS,DT > builder(dt);
P.delegate(builder);
return P;
}
template <class Polyhedron, class Plane_3>
Polyhedron clip_bbox(const Bbox_3& bbox, const Plane_3& plane)
{
typedef typename Polyhedron::Traits::Kernel Kernel;
typedef typename Kernel::Point_3 Point_3;
typedef typename Kernel::Segment_3 Segment_3;
cpp11::array<typename Kernel::Point_3,8> corners= {{
Point_3(bbox.xmin(),bbox.ymin(),bbox.zmin()),
Point_3(bbox.xmin(),bbox.ymax(),bbox.zmin()),
Point_3(bbox.xmax(),bbox.ymax(),bbox.zmin()),
Point_3(bbox.xmax(),bbox.ymin(),bbox.zmin()),
Point_3(bbox.xmin(),bbox.ymin(),bbox.zmax()),
Point_3(bbox.xmin(),bbox.ymax(),bbox.zmax()),
Point_3(bbox.xmax(),bbox.ymax(),bbox.zmax()),
Point_3(bbox.xmax(),bbox.ymin(),bbox.zmax())
}};
cpp11::array<CGAL::Oriented_side,8> orientations = {{
plane.oriented_side(corners[0]),
plane.oriented_side(corners[1]),
plane.oriented_side(corners[2]),
plane.oriented_side(corners[3]),
plane.oriented_side(corners[4]),
plane.oriented_side(corners[5]),
plane.oriented_side(corners[6]),
plane.oriented_side(corners[7])
}};
std::vector<Point_3> points;
// first look for intersections at corners
for (int i=0; i<8; ++i)
if (orientations[i]==ON_ORIENTED_BOUNDARY)
points.push_back(corners[i]);
// second look for intersections on edges
cpp11::array<int,24> edge_indices = {{ // 2 *12 edges
0,1, 1,2, 2,3, 3,0, // bottom face edges
4,5, 5,6, 6,7, 7,4, // top face edges
0,4, 1,5, 2,6, 3,7
}};
for (int i=0; i<12; ++i)
{
int i1=edge_indices[2*i], i2=edge_indices[2*i+1];
if (orientations[i1]==ON_ORIENTED_BOUNDARY) continue;
if (orientations[i2]==ON_ORIENTED_BOUNDARY) continue;
if (orientations[i1]!=orientations[i2])
points.push_back(
boost::get<Point_3>(
*CGAL::intersection(plane, Segment_3(corners[i1], corners[i2]) )
)
);
}
Polyhedron P;
//if less that 3 points there will be nothing to clipping.
if (points.size()<3) return P;
for (int i=0; i<8; ++i)
if (orientations[i]==ON_NEGATIVE_SIDE)
points.push_back(corners[i]);
// take the convex hull of the points on the negative side+intersection points
// overkill...
CGAL::convex_hull_3(points.begin(), points.end(), P);
return P;
}
template <class Polyhedron>
Polyhedron* clip_polyhedron(Polyhedron& P, Polyhedron& clipping_polyhedron)
{
std::pair <Polyhedron*,int> result;
typedef CGAL::Polyhedron_corefinement<Polyhedron> Corefinement;
Corefinement coref;
CGAL::Emptyset_iterator emptyset_iterator;
coref(P, clipping_polyhedron, emptyset_iterator,
&result, Corefinement::Intersection_tag);
return result.first;
}
template <class Polyhedron, class Plane_3>
Polyhedron* clip_polyhedron(const Polyhedron& P, const Plane_3& p)
{
if(P.empty()) return new Polyhedron();
CGAL::Bbox_3 bbox( CGAL::bbox_3(P.points_begin(), P.points_end()) );
//extend the bbox a bit to avoid border cases
double xd=(bbox.xmax()-bbox.xmin())/100;
double yd=(bbox.ymax()-bbox.ymin())/100;
double zd=(bbox.zmax()-bbox.zmin())/100;
bbox=CGAL::Bbox_3(bbox.xmin()-xd, bbox.ymin()-yd, bbox.zmin()-zd,
bbox.xmax()+xd, bbox.ymax()+yd, bbox.zmax()+zd);
Polyhedron clipping_polyhedron=clip_bbox<Polyhedron>(bbox, p);
if (clipping_polyhedron.empty()) //no intersection, result is all or nothing
{
if (p.oriented_side(*P.points_begin())==ON_POSITIVE_SIDE)
return new Polyhedron();
else
return new Polyhedron(P);
}
Polyhedron copy(P);
return clip_polyhedron(copy, clipping_polyhedron);
}
template <class Polyhedron, class Plane_3>
std::pair<Polyhedron*,Polyhedron*> split_polyhedron(const Polyhedron& P, const Plane_3& p)
{
std::pair<Polyhedron*, Polyhedron*> res;
if(P.empty()) {res.first=res.second= new Polyhedron(); return res;}
CGAL::Bbox_3 bbox( CGAL::bbox_3(P.points_begin(), P.points_end()) );
//extend the bbox a bit to avoid border cases
double xd=(bbox.xmax()-bbox.xmin())/100;
double yd=(bbox.ymax()-bbox.ymin())/100;
double zd=(bbox.zmax()-bbox.zmin())/100;
bbox=CGAL::Bbox_3(bbox.xmin()-xd, bbox.ymin()-yd, bbox.zmin()-zd,
bbox.xmax()+xd, bbox.ymax()+yd, bbox.zmax()+zd);
//First Polyhedron
Polyhedron clipping_polyhedron=clip_bbox<Polyhedron>(bbox, p);
if (clipping_polyhedron.empty()) //no intersection, result is all or nothing
{
if (p.oriented_side(*P.points_begin())==ON_POSITIVE_SIDE)
{
res.first = new Polyhedron(); res.second = new Polyhedron(P);
}
else
{
res.first = new Polyhedron(P); res.second = new Polyhedron();
}
return res;
}
Polyhedron copy(P);
res.first = clip_polyhedron(copy, clipping_polyhedron);
//Second Polyhedron
Plane_3 p_op = p.opposite();
clipping_polyhedron=clip_bbox<Polyhedron>(bbox, p_op);
copy = P;
res.second = clip_polyhedron(copy, clipping_polyhedron);
return res;
}
namespace internal{
template<class Polyhedron>
struct Edge_is_marked4coref{
std::set<typename Polyhedron::Halfedge_handle>& marked_halfedges;
typedef bool value_type;
typedef value_type reference;
typedef std::pair<typename Polyhedron::Halfedge_handle,Polyhedron*> key_type;
typedef boost::read_write_property_map_tag category;
Edge_is_marked4coref(std::set<typename Polyhedron::Halfedge_handle>& mh)
: marked_halfedges(mh)
{}
friend reference get(Edge_is_marked4coref& map,const key_type& key) {
return map.marked_halfedges.count(key.first)!=0;
}
friend void put(Edge_is_marked4coref& map,key_type key,value_type v) {
if (v) map.marked_halfedges.insert(key.first);
else map.marked_halfedges.erase(key.first);
}
};
template<class Polyhedron>
struct Edge_is_marked{
const std::set<typename Polyhedron::Halfedge_handle>* marked_halfedges;
typedef bool value_type;
typedef value_type reference;
typedef typename boost::graph_traits<Polyhedron>::edge_descriptor key_type;
typedef boost::readable_property_map_tag category;
Edge_is_marked(){}
Edge_is_marked(const std::set<typename Polyhedron::Halfedge_handle>& mh)
: marked_halfedges(&mh)
{}
friend reference get(const Edge_is_marked& map,const key_type& key) {
return map.marked_halfedges->count(key.halfedge())!=0;
}
};
} //end of internal namespace
template <class Polyhedron, class Plane_3>
void inplace_clip_open_polyhedron(Polyhedron& P, const Plane_3& p)
{
CGAL::Bbox_3 bbox( CGAL::bbox_3(P.points_begin(), P.points_end()) );
//extend the bbox a bit to avoid border cases
double xd=(bbox.xmax()-bbox.xmin())/100;
double yd=(bbox.ymax()-bbox.ymin())/100;
double zd=(bbox.zmax()-bbox.zmin())/100;
bbox=CGAL::Bbox_3(bbox.xmin()-xd, bbox.ymin()-yd, bbox.zmin()-zd,
bbox.xmax()+xd, bbox.ymax()+yd, bbox.zmax()+zd);
Polyhedron clipping_polyhedron=clip_to_bbox<Polyhedron>(bbox, p);
if (clipping_polyhedron.empty()) //no intersection, result is all or nothing
{
if (p.oriented_side(*P.points_begin())==ON_POSITIVE_SIDE)
P.clear();
return;
}
// set for marking edges of P intersected by the clipping plane
std::set<typename Polyhedron::Halfedge_handle> marked_halfedges;
internal::Edge_is_marked4coref<Polyhedron> cr_edge_is_marked(marked_halfedges);
typedef typename Polyhedron::Traits::Kernel K;
typedef CGAL::Node_visitor_refine_polyhedra<
Polyhedron,CGAL::Default,K,internal::Edge_is_marked4coref<Polyhedron> >
Split_visitor;
CGAL::Default_output_builder output_builder;
Split_visitor visitor(output_builder,
Default_polyhedron_ppmap<Polyhedron>(),
cr_edge_is_marked);
CGAL::Intersection_of_Polyhedra_3<Polyhedron,K,Split_visitor>
polyline_intersections(visitor);
CGAL::Emptyset_iterator emptyset_iterator;
// corefinement P and clipping_polyhedron
polyline_intersections(P,clipping_polyhedron,emptyset_iterator);
// extract connected components bounded by marked edges
internal::Edge_is_marked<Polyhedron> edge_is_marked(marked_halfedges);
namespace PMP=Polygon_mesh_processing;
std::map<typename Polyhedron::Face_handle, std::size_t> face_ccs;
std::size_t nb_cc=PMP::connected_components(P,
boost::make_assoc_property_map(face_ccs),
PMP::parameters::edge_is_constrained_map(edge_is_marked)
.face_index_map(get(boost::face_external_index,P))
);
// remove cc on the positive side of the plane
std::vector<bool> cc_handled(nb_cc, false);
std::vector<std::size_t> ccs_to_remove;
BOOST_FOREACH(typename Polyhedron::Face_handle f, faces(P))
{
std::size_t cc_id=face_ccs[f];
if (cc_handled[cc_id]) continue;
//look for a vertex not on the intersection
typename Polyhedron::Halfedge_handle h=f->halfedge();
for(int i=0;i<3;++i){
bool no_marked_edge=true;
BOOST_FOREACH(typename Polyhedron::Halfedge_handle h, halfedges_around_target(h, P))
if ( marked_halfedges.count(h) )
no_marked_edge=false;
if (no_marked_edge){
if ( p.oriented_side(h->vertex()->point())==ON_POSITIVE_SIDE )
ccs_to_remove.push_back(cc_id);
cc_handled[cc_id]=true;
if (--nb_cc==0) break;
break;
}
h=h->next();
}
}
//now remove the faces on the positive side
PMP::remove_connected_components(P,
ccs_to_remove,
boost::make_assoc_property_map(face_ccs),
PMP::parameters::vertex_index_map(get(boost::vertex_external_index,P))
);
}
} } // CGAL::corefinement
#endif // CGAL_INTERNAL_POLYHEDRON_PLANE_CLIPPING_3_H

View File

@ -0,0 +1,367 @@
// Copyright (c) 2016 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
//
// Author(s) : Sebastien Loriot
#ifndef CGAL_POLYGON_MESH_PROCESSING_CLIP_H
#define CGAL_POLYGON_MESH_PROCESSING_CLIP_H
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/convex_hull_3.h>
namespace CGAL{
namespace Polygon_mesh_processing {
namespace internal
{
template <class TriangleMesh,
class Ecm,
class NamedParameters1,
class NamedParameters2>
bool
clip_open_impl( TriangleMesh& tm,
TriangleMesh& clipper,
Ecm ecm,
const NamedParameters1& np_tm,
const NamedParameters2& np_c)
{
// first corefine the meshes
corefine(tm, clipper, np_tm, np_c);
typedef typename GetVertexPointMap<TriangleMesh,
NamedParameters1>::type Vpm;
typedef typename GetFaceIndexMap<TriangleMesh,
NamedParameters1>::type Fid_map;
typedef typename GetVertexIndexMap<TriangleMesh,
NamedParameters1>::type Vid_map;
typedef boost::graph_traits<TriangleMesh> GT;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters2>::type GeomTraits;
typedef typename GT::halfedge_descriptor halfedge_descriptor;
typedef typename GT::face_descriptor face_descriptor;
Fid_map fid_map = boost::choose_param(get_param(np_tm, face_index),
get_property_map(face_index, tm));
Vid_map vid_map = boost::choose_param(get_param(np_tm, boost::vertex_index),
get_property_map(boost::vertex_index, tm));
Vpm vpm1 = boost::choose_param(get_param(np_tm, vertex_point),
get_property_map(vertex_point, tm));
Vpm vpm_c = boost::choose_param(get_param(np_c, vertex_point),
get_property_map(vertex_point, clipper));
// init indices if needed
Corefinement::init_face_indices(tm, fid_map);
Corefinement::init_vertex_indices(tm, vid_map);
// set the connected component id of each face
std::vector<std::size_t> face_cc(num_faces(tm), std::size_t(-1));
std::size_t nb_cc =
connected_components(tm,
Corefinement::bind_maps(fid_map, make_property_map(face_cc)),
parameters::face_index_map(fid_map).
edge_is_constrained_map(ecm));
boost::dynamic_bitset<> cc_not_handled(nb_cc);
cc_not_handled.set();
std::vector <std::size_t> ccs_to_remove;
/// \todo clipper has been modified, this is not robust if inexact constructions are used
CGAL::Side_of_triangle_mesh<TriangleMesh, GeomTraits, Vpm>
side_of(clipper, vpm_c);
BOOST_FOREACH(face_descriptor f, faces(tm))
{
std::size_t cc_id = face_cc[ get(fid_map, f) ];
if ( !cc_not_handled.test(cc_id) ) continue;
halfedge_descriptor h=halfedge(f, tm);
for(int i=0;i<3;++i)
{
bool no_marked_edge=true;
BOOST_FOREACH(halfedge_descriptor h2, halfedges_around_target(h, tm))
if ( get(ecm, edge(h2, tm)) ){
no_marked_edge=false;
break;
}
if (no_marked_edge){
if ( side_of( get(vpm1, target(h, tm) ) ) == ON_UNBOUNDED_SIDE )
ccs_to_remove.push_back(cc_id);
cc_not_handled.reset(cc_id);
break;
}
h=next(h, tm);
}
if (!cc_not_handled.any()) break;
}
if (cc_not_handled.any())
{
///\todo handle using barycenters? won't work for coplanar faces
}
//now remove the cc
remove_connected_components(tm,
ccs_to_remove,
Corefinement::bind_maps(fid_map, make_property_map(face_cc)),
np_tm);
return true;
}
/// \TODO move this to property_map.h?
template <class Set>
struct Constrained_edge_map
{
typedef boost::read_write_property_map_tag category;
typedef bool value_type;
typedef bool reference;
typedef typename Set::value_type key_type;
Constrained_edge_map()
: edge_set(NULL)
{}
Constrained_edge_map(Set& set)
: edge_set(&set)
{}
friend bool get(const Constrained_edge_map<Set>& map, key_type k)
{
return map.edge_set->count(k)!=0;
}
friend void put(Constrained_edge_map<Set>& map, key_type k, bool b)
{
if (b)
map.edge_set->insert(k);
else
map.edge_set->erase(k);
}
private:
Set* edge_set;
};
template <class TriangleMesh,
class NamedParameters1,
class NamedParameters2>
bool
clip_open_impl( TriangleMesh& tm,
TriangleMesh& clipper,
boost::param_not_found,
const NamedParameters1& np_tm,
const NamedParameters2& np_c)
{
typedef typename boost::graph_traits<TriangleMesh>
::edge_descriptor edge_descriptor;
boost::unordered_set<edge_descriptor> constrained_edges;
Constrained_edge_map<boost::unordered_set<edge_descriptor> >
cst_map(constrained_edges);
return clip_open_impl(tm, clipper,
cst_map,
np_tm.edge_is_constrained_map(cst_map),
np_c);
}
} // end of internal namespace
#ifndef DOXYGEN_RUNNING
///\todo clipper const!
/// requires face_index_map, vertex_index_map for np_tm
/// requires face_index_map for np_c
/// if edge_is_constrained_map is not provided in np_tm a default one is
/// provided using boost::unordered_set<edge_descriptor>
template <class TriangleMesh,
class NamedParameters1,
class NamedParameters2>
bool
clip( TriangleMesh& tm,
/*const*/ TriangleMesh& clipper,
bool close,
const NamedParameters1& np_tm,
const NamedParameters2& np_c)
{
if (close && is_closed(tm))
return corefine_and_compute_intersection(tm, clipper, tm, np_tm, np_c);
return internal::clip_open_impl(tm, clipper,
get_param(np_tm, edge_is_constrained), np_tm, np_c);
}
/// \todo document me
template <class Plane_3,
class TriangleMesh,
class NamedParameters>
Oriented_side
clip_to_bbox(const Plane_3& plane,
const Bbox_3& bbox,
TriangleMesh& tm_out,
const NamedParameters& np )
{
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Geom_traits;
typedef typename Geom_traits::Point_3 Point_3;
typedef typename Geom_traits::Segment_3 Segment_3;
typedef typename GetVertexPointMap<TriangleMesh,
NamedParameters>::type Vpm;
Vpm vpm_out = boost::choose_param(get_param(np, vertex_point),
get_property_map(vertex_point, tm_out));
cpp11::array<Point_3,8> corners= {{
Point_3(bbox.xmin(),bbox.ymin(),bbox.zmin()),
Point_3(bbox.xmin(),bbox.ymax(),bbox.zmin()),
Point_3(bbox.xmax(),bbox.ymax(),bbox.zmin()),
Point_3(bbox.xmax(),bbox.ymin(),bbox.zmin()),
Point_3(bbox.xmin(),bbox.ymin(),bbox.zmax()),
Point_3(bbox.xmin(),bbox.ymax(),bbox.zmax()),
Point_3(bbox.xmax(),bbox.ymax(),bbox.zmax()),
Point_3(bbox.xmax(),bbox.ymin(),bbox.zmax())
}};
cpp11::array<CGAL::Oriented_side,8> orientations = {{
plane.oriented_side(corners[0]),
plane.oriented_side(corners[1]),
plane.oriented_side(corners[2]),
plane.oriented_side(corners[3]),
plane.oriented_side(corners[4]),
plane.oriented_side(corners[5]),
plane.oriented_side(corners[6]),
plane.oriented_side(corners[7])
}};
std::vector<Point_3> points;
// look for intersections on edges
cpp11::array<int,24> edge_indices = {{ // 2 *12 edges
0,1, 1,2, 2,3, 3,0, // bottom face edges
4,5, 5,6, 6,7, 7,4, // top face edges
0,4, 1,5, 2,6, 3,7
}};
for (int i=0; i<12; ++i)
{
int i1=edge_indices[2*i], i2=edge_indices[2*i+1];
if (orientations[i1]==ON_ORIENTED_BOUNDARY) continue;
if (orientations[i2]==ON_ORIENTED_BOUNDARY) continue;
if (orientations[i1]!=orientations[i2])
points.push_back(
boost::get<Point_3>(
*intersection(plane, Segment_3(corners[i1], corners[i2]) )
)
);
}
Oriented_side last_os = ON_ORIENTED_BOUNDARY;
for (int i=0; i<8; ++i)
if (orientations[i]!=ON_ORIENTED_BOUNDARY)
{
if (last_os==ON_ORIENTED_BOUNDARY)
last_os=orientations[i];
else
{
if(last_os!=orientations[i])
{
last_os=ON_ORIENTED_BOUNDARY;
break;
}
}
}
// the intersection is the full bbox
if (last_os!=ON_ORIENTED_BOUNDARY)
return last_os;
//add points on negative side and on the plane
for (int i=0; i<8; ++i)
if (orientations[i]!=ON_POSITIVE_SIDE)
points.push_back(corners[i]);
// take the convex hull of the points on the negative side+intersection points
// overkill...
Polyhedron_3<Geom_traits> P;
CGAL::convex_hull_3(points.begin(), points.end(), P);
copy_face_graph(P, tm_out,
Emptyset_iterator(), Emptyset_iterator(), Emptyset_iterator(),
get(vertex_point, P), vpm_out);
return ON_ORIENTED_BOUNDARY;
}
// convenience overload
template <class TriangleMesh,
class NamedParameters1>
bool
clip( TriangleMesh& tm,
/*const*/ TriangleMesh& clipper,
bool close,
const NamedParameters1& np_tm)
{
return clip(tm, clipper, close, np_tm, parameters::all_default());
}
// convenience overload
template <class TriangleMesh>
bool
clip( TriangleMesh& tm,
/*const*/ TriangleMesh& clipper,
bool close)
{
return clip(tm, clipper, close, parameters::all_default());
}
// works only with the default point map, for more complex use cases, use
// clip_to_bbox first and the other overload of clip with two meshes
/// \todo document me
template <class TriangleMesh,
class Plane_3>
void clip( TriangleMesh& tm,
const Plane_3& plane,
bool close)
{
if( boost::begin(faces(tm))==boost::end(faces(tm)) ) return;
CGAL::Bbox_3 bbox = bbox_3(tm);
//extend the bbox a bit to avoid border cases
double xd=(bbox.xmax()-bbox.xmin())/100;
double yd=(bbox.ymax()-bbox.ymin())/100;
double zd=(bbox.zmax()-bbox.zmin())/100;
bbox=CGAL::Bbox_3(bbox.xmin()-xd, bbox.ymin()-yd, bbox.zmin()-zd,
bbox.xmax()+xd, bbox.ymax()+yd, bbox.zmax()+zd);
TriangleMesh clipper;
Oriented_side os = clip_to_bbox(plane, bbox, clipper, parameters::all_default());
switch(os)
{
case ON_NEGATIVE_SIDE:
return; // nothing to clip, the full mesh is on the negative side
case ON_POSITIVE_SIDE:
clear(tm); // clear the mesh that is fully on the positive side
return;
default:
clip(tm, clipper, close);
}
}
#endif // !DOXYGEN_RUNNING
} } //end of namespace CGAL::Polygon_mesh_processing
#endif // CGAL_POLYGON_MESH_PROCESSING_CLIP_H

View File

@ -81,7 +81,7 @@ if (EIGEN3_FOUND)
create_single_source_cgal_program("test_corefinement_bool_op.cpp") create_single_source_cgal_program("test_corefinement_bool_op.cpp")
create_single_source_cgal_program("test_corefine.cpp") create_single_source_cgal_program("test_corefine.cpp")
create_single_source_cgal_program("test_does_bound_a_volume.cpp") create_single_source_cgal_program("test_does_bound_a_volume.cpp")
create_single_source_cgal_program("polyhedron_clipping.cpp") create_single_source_cgal_program("test_pmp_clip.cpp")
if(NOT (${EIGEN3_VERSION} VERSION_LESS 3.2.0)) if(NOT (${EIGEN3_VERSION} VERSION_LESS 3.2.0))
create_single_source_cgal_program("fairing_test.cpp") create_single_source_cgal_program("fairing_test.cpp")

View File

@ -1,67 +0,0 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/internal/Polyhedron_plane_clipping_3.h>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
int main(int argc, char** argv)
{
const char* fname = argc>1?argv[1]:"data/elephant.off";
std::ifstream input(fname);
Polyhedron P;
input >> P;
std::cout << "Polyhedron with " << P.size_of_vertices() << " vertices\n";
// test clipping by a plane
{
Polyhedron* P_clipped=
CGAL::corefinement::clip_polyhedron(P, Kernel::Plane_3(1,1,0,0));
std::ofstream output("clipped1.off");
output << *P_clipped;
delete P_clipped;
}
// test clipping with the opposite plane
{
Polyhedron* P_clipped=
CGAL::corefinement::clip_polyhedron(P, Kernel::Plane_3(-1,-1,0,0));
std::ofstream output("clipped2.off");
output << *P_clipped;
delete P_clipped;
}
// test clipping with a plane with no intersection having the elephant
// on its positive side, result should be empty
{
Polyhedron* P_clipped=
CGAL::corefinement::clip_polyhedron(P, Kernel::Plane_3(1,0,0,1));
if (!P_clipped->empty()){
std::cerr << "Error: Polyhedron should be empty!\n";
return 1;
}
delete P_clipped;
}
// test clipping with a plane with no intersection having the elephant
// on its negative side, result should be empty
{
Polyhedron* P_clipped=
CGAL::corefinement::clip_polyhedron(P, Kernel::Plane_3(-1,0,0,-1));
if (P_clipped->size_of_vertices()!=P.size_of_vertices()){
std::cerr << "Error: Polyhedron should be full!\n";
return 1;
}
delete P_clipped;
}
// clip with a plane but do not close the polyhedron
{
CGAL::corefinement::inplace_clip_open_polyhedron(P, Kernel::Plane_3(1,1,0,0));
std::ofstream output("open_clipped.off");
output << P;
}
return 0;
}

View File

@ -0,0 +1,65 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/internal/clip.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <iostream>
#include <fstream>
namespace PMP = CGAL::Polygon_mesh_processing;
namespace params = PMP::parameters;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Triangle_mesh;
typedef CGAL::Polyhedron_3<K> Polyhedron;
typedef Triangle_mesh::Property_map<Triangle_mesh::Edge_index,bool> Constrained_edge_map;
int main()
{
{
// test open clipping with Surface_mesh
Triangle_mesh tm1, tm2;
std::ifstream input("data-coref/elephant.off");
input >> tm1;
input.close();
input.open("data-coref/sphere.off");
input >> tm2;
input.close();
Constrained_edge_map ecm1 =
tm1.add_property_map<Triangle_mesh::Edge_index,bool>("e:cst", false).first;
PMP::clip(tm1, tm2, false, params::edge_is_constrained_map(ecm1));
std::ofstream output("clipped_opened.off");
output << tm1;
// test open clipping with Polyhedron
Polyhedron P, Q;
input.open("data-coref/elephant.off");
input >> P;
input.close();
input.open("data-coref/sphere.off");
input >> Q;
PMP::clip(P, Q, false,
params::face_index_map(get(CGAL::face_external_index, P)).
vertex_index_map(get(CGAL::vertex_external_index, P)),
params::face_index_map(get(CGAL::face_external_index, Q)));
assert(P.size_of_vertices() == tm1.number_of_vertices());
}
{
Triangle_mesh tm1, tm2;
std::ifstream input("data-coref/elephant.off");
input >> tm1;
input.close();
input.open("data-coref/sphere.off");
input >> tm2;
PMP::clip(tm1, tm2, true);
std::ofstream output("clipped_closed.off");
output << tm1;
}
return 0;
}

View File

@ -7,7 +7,8 @@
#include "Scene_plane_item.h" #include "Scene_plane_item.h"
#include <CGAL/Three/Viewer_interface.h> #include <CGAL/Three/Viewer_interface.h>
#include <CGAL/Three/Polyhedron_demo_plugin_interface.h> #include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
#include <CGAL/internal/Polyhedron_plane_clipping_3.h> #include <CGAL/Polygon_mesh_processing/internal/clip.h>
#include "ui_Clip_polyhedron_plugin.h" #include "ui_Clip_polyhedron_plugin.h"
#include "Viewer.h" #include "Viewer.h"
@ -187,100 +188,43 @@ public Q_SLOTS:
//apply the clipping function //apply the clipping function
Q_FOREACH(Scene_polyhedron_item* poly, polyhedra) Q_FOREACH(Scene_polyhedron_item* poly, polyhedra)
{ {
if (ui_widget.clip_radioButton->isChecked())
if(ui_widget.close_checkBox->isChecked() && poly->polyhedron()->is_closed())
{ {
std::pair<Polyhedron*, Polyhedron*>polyhedron; CGAL::Polygon_mesh_processing::clip(*(poly->polyhedron()),
if(ui_widget.clip_radioButton->isChecked()) plane->plane(),
{ ui_widget.close_checkBox->isChecked());
polyhedron.first = CGAL::corefinement::clip_polyhedron(*(poly->polyhedron()),plane->plane()); poly->invalidateOpenGLBuffers();
polyhedron.second = NULL; viewer->updateGL();
}
else
{
polyhedron = CGAL::corefinement::split_polyhedron(*(poly->polyhedron()),plane->plane());
}
if(polyhedron.first != NULL)
{
Scene_polyhedron_item* new_item = new Scene_polyhedron_item(polyhedron.first);
if(polyhedron.second != NULL)
new_item->setName(QString("%1 %2").arg(poly->name()).arg("1"));
else
new_item->setName(poly->name());
new_item->setColor(poly->color());
new_item->setRenderingMode(poly->renderingMode());
new_item->setVisible(poly->visible());
new_item->invalidateOpenGLBuffers();
new_item->setProperty("source filename", poly->property("source filename"));
scene->replaceItem(scene->item_id(poly),new_item);
new_item->invalidateOpenGLBuffers();
viewer->updateGL();
if(ui_widget.clip_radioButton->isChecked())
messages->information(QString("%1 clipped").arg(new_item->name()));
}
else
{
messages->information(QString("Could not clip %1 : returned polyhedron is null.").arg(poly->name()));
delete polyhedron.first;
}
if(polyhedron.second!= NULL)
{
Scene_polyhedron_item* new_item = new Scene_polyhedron_item(polyhedron.second);
new_item->setName(QString("%1 %2").arg(poly->name()).arg("2"));
new_item->setColor(poly->color());
new_item->setRenderingMode(poly->renderingMode());
new_item->setVisible(poly->visible());
new_item->invalidateOpenGLBuffers();
new_item->setProperty("source filename", poly->property("source filename"));
scene->addItem(new_item);
new_item->invalidateOpenGLBuffers();
viewer->updateGL();
if(!ui_widget.clip_radioButton->isChecked())
messages->information(QString("%1 splitted").arg(poly->name()));
}
if(polyhedron.first != NULL)
delete poly;
} }
else else
{ {
Scene_polyhedron_item* poly1 = new Scene_polyhedron_item(*(poly->polyhedron())); //part on the negative side
poly1->setProperty("source filename", poly->property("source filename")); Polyhedron* neg_side = new Polyhedron(*poly->polyhedron());
Scene_polyhedron_item* poly2 = NULL;
if(!ui_widget.clip_radioButton->isChecked())
{
poly2 = new Scene_polyhedron_item(*(poly->polyhedron()));
poly2->setProperty("source filename", poly->property("source filename"));
}
CGAL::Polygon_mesh_processing::clip(*neg_side,
plane->plane(),
ui_widget.close_checkBox->isChecked());
Scene_polyhedron_item* new_item = new Scene_polyhedron_item(neg_side);
new_item->setName(QString("%1 on %2").arg(poly->name()).arg("negative side"));
new_item->setColor(poly->color());
new_item->setRenderingMode(poly->renderingMode());
new_item->setVisible(poly->visible());
scene->addItem(new_item);
new_item->invalidateOpenGLBuffers();
// part on the positive side
Polyhedron* pos_side = new Polyhedron(*poly->polyhedron());
CGAL::Polygon_mesh_processing::clip(*pos_side,
plane->plane().opposite(),
ui_widget.close_checkBox->isChecked());
new_item = new Scene_polyhedron_item(pos_side);
new_item->setName(QString("%1 on %2").arg(poly->name()).arg("positive side"));
new_item->setColor(poly->color());
new_item->setRenderingMode(poly->renderingMode());
new_item->setVisible(poly->visible());
scene->addItem(new_item);
new_item->invalidateOpenGLBuffers();
CGAL::corefinement::inplace_clip_open_polyhedron(*(poly1->polyhedron()),plane->plane()); viewer->updateGL();
if(poly2 != NULL)
poly1->setName(QString("%1 %2").arg(poly->name()).arg("1"));
else
poly1->setName(poly->name());
poly1->setColor(poly->color());
poly1->setRenderingMode(poly->renderingMode());
poly1->setVisible(poly->visible());
scene->replaceItem(scene->item_id(poly),poly1);
poly1->invalidateOpenGLBuffers();
viewer->updateGL();
if(ui_widget.clip_radioButton->isChecked())
messages->information(QString("%1 clipped").arg(poly->name()));
if(poly2 != NULL)
{
CGAL::corefinement::inplace_clip_open_polyhedron(*(poly2->polyhedron()),plane->plane().opposite());
poly2->setName(QString("%1 %2").arg(poly->name()).arg("2"));
poly2->setColor(poly->color());
poly2->setRenderingMode(poly->renderingMode());
poly2->setVisible(poly->visible());
scene->addItem(poly2);
poly2->invalidateOpenGLBuffers();
viewer->updateGL();
messages->information(QString("%1 splitted").arg(poly->name()));
}
delete poly;
} }
} }
} }