1 // Copyright (c) 2014 INRIA Sophia-Antipolis (France). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h $ 7 // $Id: halfspace_intersection_3.h 278a26d 2020-05-18T18:01:56+02:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Jocelyn Meyron 12 // 13 14 #ifndef CGAL_HALFSPACE_INTERSECTION_3_H 15 #define CGAL_HALFSPACE_INTERSECTION_3_H 16 17 #include <CGAL/license/Convex_hull_3.h> 18 19 #include <CGAL/disable_warnings.h> 20 21 #include <CGAL/HalfedgeDS_default.h> 22 #include <CGAL/Convex_hull_3/dual/Convex_hull_traits_dual_3.h> 23 #include <CGAL/Origin.h> 24 #include <CGAL/convex_hull_3.h> 25 #include <CGAL/intersections.h> 26 #include <CGAL/assertions.h> 27 #include <CGAL/boost/graph/Euler_operations.h> 28 // For interior_polyhedron_3 29 #include <CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h> 30 #include <CGAL/internal/Exact_type_selector.h> 31 32 #include <boost/unordered_map.hpp> 33 #include <boost/type_traits/is_floating_point.hpp> 34 #include <deque> 35 36 namespace CGAL 37 { 38 namespace Convex_hull_3 39 { 40 namespace internal 41 { 42 // Build the primal polyhedron associated to a dual polyhedron 43 // We also need the `origin` which represents a point inside the primal polyhedron 44 template <class Polyhedron_dual, class Polyhedron, class Point_3> 45 void build_primal_polyhedron(Polyhedron & primal,const Polyhedron_dual & _dual,const Point_3 & origin)46 build_primal_polyhedron(Polyhedron& primal, 47 const Polyhedron_dual & _dual, 48 const Point_3& origin) 49 { 50 typedef typename Kernel_traits<Point_3>::Kernel Kernel; 51 typedef typename Kernel::RT RT; 52 53 // Typedefs for dual 54 55 typedef typename boost::graph_traits<Polyhedron_dual>::face_descriptor 56 Facet_const_handle; 57 typedef typename boost::graph_traits<Polyhedron_dual>::vertex_descriptor 58 Vertex_const_descriptor; 59 typedef typename boost::graph_traits<Polyhedron_dual>::halfedge_descriptor 60 Halfedge_const_descriptor; 61 62 // typedef and type for primal 63 typedef typename boost::graph_traits<Polyhedron>::vertex_descriptor vertex_descriptor; 64 typename boost::property_map<Polyhedron, vertex_point_t>::type vpm = 65 get(CGAL::vertex_point, primal); 66 67 // Typedefs for intersection 68 typedef typename Kernel::Plane_3 Plane_3; 69 typedef typename Kernel::Line_3 Line_3; 70 typedef boost::optional< boost::variant< Point_3, 71 Line_3, 72 Plane_3 > > result_inter; 73 74 75 boost::unordered_map <Facet_const_handle, vertex_descriptor> primal_vertices; 76 size_t n = 0; 77 78 // First, computing the primal vertices 79 for(Facet_const_handle fd : faces(_dual)){ 80 Halfedge_const_descriptor h = fd->halfedge(); 81 // Build the dual plane corresponding to the current facet 82 Plane_3 p1 = h->vertex()->point(); 83 Plane_3 p2 = h->next()->vertex()->point(); 84 Plane_3 p3 = h->next()->next()->vertex()->point(); 85 86 RT dp1 = p1.d() + origin.x() * p1.a() 87 + origin.y() * p1.b() + origin.z() * p1.c(); 88 RT dp2 = p2.d() + origin.x() * p2.a() 89 + origin.y() * p2.b() + origin.z() * p2.c(); 90 RT dp3 = p3.d() + origin.x() * p3.a() 91 + origin.y() * p3.b() + origin.z() * p3.c(); 92 93 Plane_3 pp1(p1.a(), p1.b(), p1.c(), dp1); 94 Plane_3 pp2(p2.a(), p2.b(), p2.c(), dp2); 95 Plane_3 pp3(p3.a(), p3.b(), p3.c(), dp3); 96 97 // Compute the intersection 98 result_inter result = CGAL::intersection(pp1, pp2, pp3); 99 CGAL_assertion_msg(bool(result), 100 "halfspace_intersection_3: no intersection"); 101 CGAL_assertion_msg(boost::get<Point_3>(& *result) != nullptr, 102 "halfspace_intersection_3: intersection is not a point"); 103 104 const Point_3* pp = boost::get<Point_3>(& *result); 105 106 // Primal vertex associated to the current dual plane 107 Point_3 ppp(origin.x() + pp->x(), 108 origin.y() + pp->y(), 109 origin.z() + pp->z()); 110 111 vertex_descriptor vd = add_vertex(primal); 112 primal_vertices[fd] = vd; 113 put(vpm, vd, ppp); 114 ++n; 115 } 116 117 // Then, add facets to the primal polyhedron 118 // To do this, for each dual vertex, we circulate around this vertex 119 // and we add an edge between each facet we encounter 120 121 for(Vertex_const_descriptor vd : vertices( _dual)) { 122 std::deque<vertex_descriptor> vertices; 123 for(Halfedge_const_descriptor hd : halfedges_around_target(vd, _dual)){ 124 vertices.push_front(primal_vertices[face(hd, _dual)]); 125 } 126 Euler::add_face(vertices,primal); 127 } 128 } 129 130 // Test if a point is inside a convex polyhedron 131 template <class Polyhedron, class Point> point_inside_convex_polyhedron(const Polyhedron & P,Point const & p)132 bool point_inside_convex_polyhedron (const Polyhedron &P, 133 Point const& p) { 134 // Compute the equations of the facets of the polyhedron 135 typedef typename boost::graph_traits<Polyhedron>::face_descriptor face_descriptor; 136 typedef typename boost::graph_traits<Polyhedron>::halfedge_descriptor halfedge_descriptor; 137 138 typename boost::property_map<Polyhedron, vertex_point_t>::const_type vpmap = get(CGAL::vertex_point, P); 139 140 for(face_descriptor fd : faces(P)) 141 { 142 halfedge_descriptor h = halfedge(fd,P), done(h); 143 Point const& p1 = get(vpmap, target(h,P)); 144 h = next(h,P); 145 Point const& p2 = get(vpmap, target(h,P)); 146 h = next(h,P); 147 Point p3 = get(vpmap, target(h,P)); 148 149 while( h!=done && collinear(p1,p2,p3) ) 150 { 151 h = next(h,P); 152 p3 = get(vpmap, target(h,P)); 153 } 154 if (h==done) continue; //degenerate facet, skip it 155 156 if ( orientation (p1, p2, p3, p) != CGAL::NEGATIVE ) return false; 157 } 158 159 return true; 160 } 161 162 template <class Plane> collinear_plane(Plane u,Plane v)163 bool collinear_plane (Plane u, Plane v) { 164 typedef typename Kernel_traits<Plane>::Kernel Kernel; 165 typedef typename Kernel::Vector_3 Vector; 166 167 Vector uu = u.orthogonal_vector(); 168 Vector vv = v.orthogonal_vector(); 169 170 return CGAL::cross_product(uu, vv) == Vector(0, 0, 0); 171 } 172 173 template <class Plane> coplanar_plane(Plane u,Plane v,Plane w)174 bool coplanar_plane (Plane u, Plane v, Plane w) { 175 typedef typename Kernel_traits<Plane>::Kernel Kernel; 176 typedef typename Kernel::Vector_3 Vector; 177 178 Vector uu = u.orthogonal_vector(); 179 Vector vv = v.orthogonal_vector(); 180 Vector ww = w.orthogonal_vector(); 181 182 return CGAL::orientation(uu, vv, ww) == CGAL::COPLANAR; 183 } 184 185 // Checks if the dimension of intersection of the planes 186 // is a polyhedron (dimension == 3) 187 template <class InputPlaneIterator> is_intersection_dim_3(InputPlaneIterator begin,InputPlaneIterator end)188 bool is_intersection_dim_3 (InputPlaneIterator begin, 189 InputPlaneIterator end) { 190 typedef typename std::iterator_traits<InputPlaneIterator>::value_type Plane; 191 typedef typename std::vector<Plane>::iterator PlaneIterator; 192 193 std::vector<Plane> planes(begin, end); 194 // Remove same planes 195 std::size_t size = planes.size(); 196 197 // At least 4 points 198 if (size < 4) 199 return false; 200 201 // Look for two non-parallel planes 202 PlaneIterator plane1_it = planes.begin(); 203 PlaneIterator plane2_it = std::next(planes.begin()); 204 205 while (plane2_it != planes.end() && 206 collinear_plane(*plane1_it, *plane2_it)) { 207 ++plane2_it; 208 } 209 210 if (plane2_it == planes.end()) return false; 211 212 PlaneIterator plane3_it = std::next(plane2_it); 213 214 // Look for a triple of planes intersecting in a point 215 while (plane3_it != planes.end() && 216 coplanar_plane(*plane1_it, *plane2_it, *plane3_it)) { 217 ++plane3_it; 218 } 219 220 if (plane3_it == planes.end()) return false; 221 222 return true; 223 } 224 } // namespace internal 225 } // namespace Convex_hull_3 226 227 // Compute the intersection of halfspaces. 228 // If the user gives an origin then the function used it, otherwise, it is 229 // computed using a linear program. 230 template <class PlaneIterator, class Polyhedron> 231 void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, 232 Polyhedron &P, 233 boost::optional<typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3> origin = boost::none) { 234 // Checks whether the intersection is a polyhedron 235 CGAL_assertion_msg(Convex_hull_3::internal::is_intersection_dim_3(begin, end), "halfspace_intersection_3: intersection not a polyhedron"); 236 237 // Types 238 typedef typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel K; 239 typedef Convex_hull_3::Convex_hull_traits_dual_3<K> Hull_traits_dual_3; 240 typedef HalfedgeDS_default<Hull_traits_dual_3, HalfedgeDS_items_3 > Polyhedron_dual_3; 241 242 // if a point inside is not provided find one using linear programming 243 if (!origin) { 244 // find a point inside the intersection 245 origin = halfspace_intersection_interior_point_3(begin, end); 246 247 CGAL_assertion_msg(origin!=boost::none, "halfspace_intersection_3: problem when determing a point inside the intersection"); 248 if (origin==boost::none) 249 return; 250 } 251 252 // make sure the origin is on the negative side of all the planes 253 CGAL_assertion_code(for(PlaneIterator pit=begin;pit!=end;++pit)) 254 CGAL_assertion(pit->has_on_negative_side(*origin)); 255 256 // compute the intersection of the half-space using the dual formulation 257 Hull_traits_dual_3 dual_traits(*origin); 258 Polyhedron_dual_3 dual_convex_hull; 259 CGAL::convex_hull_3(begin, end, dual_convex_hull, dual_traits); 260 Convex_hull_3::internal::build_primal_polyhedron(P, dual_convex_hull, *origin); 261 262 // Posterior check if the origin is inside the computed polyhedron 263 // The check is done only if the number type is not float or double because in that 264 // case we know the construction of dual points is not exact 265 CGAL_assertion_msg( 266 boost::is_floating_point<typename K::FT>::value || 267 Convex_hull_3::internal::point_inside_convex_polyhedron(P, *origin), 268 "halfspace_intersection_3: origin not in the polyhedron" 269 ); 270 } 271 272 #ifndef CGAL_NO_DEPRECATED_CODE 273 // Compute the intersection of halfspaces (an interior point is given.) 274 template <class PlaneIterator, class Polyhedron> halfspace_intersection_3(PlaneIterator begin,PlaneIterator end,Polyhedron & P,typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3 const & origin)275 void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, 276 Polyhedron &P, 277 typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3 const& origin) { 278 halfspace_intersection_3(begin, end , P, boost::make_optional(origin)); 279 } 280 #endif 281 } // namespace CGAL 282 283 #include <CGAL/enable_warnings.h> 284 285 #endif // CGAL_HALFSPACE_INTERSECTION_3_H 286 287