1 /**
2  *   SFCGAL
3  *
4  *   Copyright (C) 2012-2013 Oslandia <infos@oslandia.com>
5  *   Copyright (C) 2012-2013 IGN (http://www.ign.fr)
6  *
7  *   This library is free software; you can redistribute it and/or
8  *   modify it under the terms of the GNU Library General Public
9  *   License as published by the Free Software Foundation; either
10  *   version 2 of the License, or (at your option) any later version.
11  *
12  *   This library is distributed in the hope that it will be useful,
13  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15  *   Library General Public License for more details.
16 
17  *   You should have received a copy of the GNU Library General Public
18  *   License along with this library; if not, see <http://www.gnu.org/licenses/>.
19  */
20 
21 
22 #include <CGAL/intersections.h>
23 
24 #include <CGAL/Polygon_mesh_processing/corefinement.h>
25 #include <CGAL/Polygon_mesh_processing/clip.h>
26 #include <CGAL/Polygon_mesh_processing/connected_components.h>
27 
28 #include <SFCGAL/detail/GeometrySet.h>
29 #include <SFCGAL/algorithm/intersects.h>
30 #include <SFCGAL/algorithm/intersection.h>
31 #include <SFCGAL/detail/triangulate/triangulateInGeometrySet.h>
32 
33 #include <CGAL/IO/Polyhedron_iostream.h>
34 
35 #include <CGAL/Side_of_triangle_mesh.h>
36 
37 using namespace SFCGAL::detail;
38 
39 namespace SFCGAL {
40 namespace algorithm {
41 
_intersection_solid_segment(const PrimitiveHandle<3> & pa,const PrimitiveHandle<3> & pb,GeometrySet<3> & output)42 void _intersection_solid_segment( const PrimitiveHandle<3>& pa, const PrimitiveHandle<3>& pb, GeometrySet<3>& output )
43 {
44     //typedef CGAL::Polyhedral_mesh_domain_3<MarkedPolyhedron, Kernel> Mesh_domain;
45 
46     const MarkedPolyhedron* ext_poly = pa.as<MarkedPolyhedron>();
47     BOOST_ASSERT( ext_poly->is_closed() );
48     const CGAL::Segment_3<Kernel>* segment = pb.as<CGAL::Segment_3<Kernel> >();
49 
50     MarkedPolyhedron* ext_poly_nc = const_cast<MarkedPolyhedron*>( ext_poly );
51     CGAL::Side_of_triangle_mesh<MarkedPolyhedron, Kernel> is_in_ext( *ext_poly_nc );
52 
53     GeometrySet<3> triangles;
54     GeometrySet<3> spoint( segment->source() );
55     GeometrySet<3> tpoint( segment->target() );
56     triangulate::triangulate( *ext_poly, triangles );
57 
58     bool source_inside = ( is_in_ext( segment->source() ) != CGAL::ON_UNBOUNDED_SIDE ) || intersects( triangles, spoint );
59     bool target_inside = ( is_in_ext( segment->target() ) != CGAL::ON_UNBOUNDED_SIDE ) || intersects( triangles, tpoint );
60 
61     if ( source_inside && target_inside ) {
62         // the entire segment intersects the volume, return the segment
63         output.addPrimitive( pb );
64     }
65     else {
66         GeometrySet<3> triangles, g;
67         triangulate::triangulate( *ext_poly, triangles );
68         g.addPrimitive( pb );
69         GeometrySet<3> inter;
70 
71         // call recursively on triangulated polyhedron
72         intersection( g, triangles, inter );
73 
74         if ( ! inter.points().empty() ) {
75             // the intersection is a point, build a segment from that point to the other end
76             if ( !source_inside && target_inside ) {
77                 CGAL::Segment_3<Kernel> interSeg( inter.points().begin()->primitive(),
78                                                   segment->target() );
79 
80                 if ( interSeg.source() == interSeg.target() ) {
81                     output.addPrimitive( segment->target() );
82                 }
83                 else {
84                     output.addPrimitive( interSeg );
85                 }
86             }
87             else if ( source_inside && !target_inside ) {
88                 CGAL::Segment_3<Kernel> interSeg( segment->source(),
89                                                   inter.points().begin()->primitive() );
90 
91                 if ( interSeg.source() == interSeg.target() ) {
92                     output.addPrimitive( segment->source() );
93                 }
94                 else {
95                     output.addPrimitive( interSeg );
96                 }
97             }
98             else { // !source_inside && !target_inside => intersection on a point
99                 output.addPrimitive( inter.points().begin()->primitive() );
100             }
101         }
102 
103         if ( ! inter.segments().empty() ) {
104             // the intersection is a segment
105             output.addPrimitive( inter.segments().begin()->primitive() );
106         }
107     }
108 }
109 
110 typedef std::vector<Kernel::Point_3> Polyline_3;
111 
112 struct Is_not_marked {
operator ()SFCGAL::algorithm::Is_not_marked113     bool operator()( MarkedPolyhedron::Halfedge_const_handle h ) const {
114         return !h->mark;
115     }
116 };
117 
118 
_intersection_solid_triangle(const MarkedPolyhedron & pa,const CGAL::Triangle_3<Kernel> & tri,GeometrySet<3> & output)119 void _intersection_solid_triangle( const MarkedPolyhedron& pa, const CGAL::Triangle_3<Kernel>& tri, GeometrySet<3>& output )
120 {
121     BOOST_ASSERT( pa.is_closed() );
122 
123     MarkedPolyhedron polyb;
124     polyb.make_triangle( tri.vertex( 0 ), tri.vertex( 1 ), tri.vertex( 2 ) );
125 
126     MarkedPolyhedron polya =  pa;
127     CGAL::Side_of_triangle_mesh<MarkedPolyhedron, Kernel> side_of_tm(polya);
128 
129     std::list<Polyline_3> polylines;
130     CGAL::Polygon_mesh_processing::surface_intersection(polya, polyb, std::back_inserter(polylines));
131     if ( polylines.size() == 0 ) {
132         // no surface intersection
133         // if one of the point of the triangle is inside the polyhedron,
134         // the triangle is inside
135         if ( side_of_tm( tri.vertex( 0 ) ) != CGAL::ON_UNBOUNDED_SIDE )
136             output.addPrimitive( tri );
137         return;
138     }
139 
140     CGAL::Polygon_mesh_processing::clip(polyb, polya, CGAL::parameters::use_compact_clipper(true));
141 
142     bool hasSurface = false;
143 
144     std::vector<MarkedPolyhedron> ccs;
145     CGAL::Polygon_mesh_processing::split_connected_components(polyb, ccs);
146 
147     for (const MarkedPolyhedron& mp : ccs)
148     {
149       // check if all vertices are on polya
150       bool all_on = true;
151       for (auto v : vertices(mp))
152         if ( side_of_tm(v->point()) != CGAL::ON_BOUNDARY )
153         {
154           all_on=false;
155           break;
156         }
157       if (all_on)
158         output.addPrimitive(mp);
159       else
160       {
161         hasSurface=true;
162         output.addPrimitive(mp, FLAG_IS_PLANAR );
163       }
164     }
165 
166     if ( hasSurface ) {
167         return;
168     }
169 
170     for ( std::list<Polyline_3>::const_iterator lit = polylines.begin(); lit != polylines.end(); ++lit ) {
171         if ( lit->size() == 1 ) {
172             // it's a point
173             output.addPrimitive( ( *lit )[0] );
174         }
175         else {
176             for ( size_t k = 1; k < lit->size(); ++k ) {
177                 CGAL::Segment_3<Kernel> seg( ( *lit )[k-1], ( *lit )[k] );
178                 output.addPrimitive( seg );
179             }
180         }
181     }
182 }
183 
_intersection_solid_solid(const MarkedPolyhedron & pa,const MarkedPolyhedron & pb,GeometrySet<3> & output)184 void _intersection_solid_solid( const MarkedPolyhedron& pa, const MarkedPolyhedron& pb, GeometrySet<3>& output )
185 {
186     // 1. find intersections on surfaces
187     // CGAL corefinement or polyhedra_intersection do not return polygon intersections between two solids
188     // they only return points, lines and volumes (but no surfaces ...)
189     {
190         // call CGAL::intersection on triangles
191         GeometrySet<3> gsa, gsb;
192         // convert polyhedra to geometry sets
193         // (no actual triangulation is done if the polyhedra are pure_triangle()
194         triangulate::triangulate( pa, gsa );
195         triangulate::triangulate( pb, gsb );
196         // "recurse" call on triangles
197         algorithm::intersection( gsa, gsb, output );
198     }
199 
200     // 2. find intersections in volumes
201     {
202         MarkedPolyhedron polya = pa, polyb = pb;
203         if ( CGAL::Polygon_mesh_processing::corefine_and_compute_intersection(polya, polyb, polya) )
204           if (std::next(vertices(polya).first)!=vertices(polya).second)
205             output.addPrimitive(polya);
206     }
207 }
208 
209 //
210 // must be called with pa's dimension larger than pb's
intersection(const PrimitiveHandle<3> & pa,const PrimitiveHandle<3> & pb,GeometrySet<3> & output,dim_t<3>)211 void intersection( const PrimitiveHandle<3>& pa, const PrimitiveHandle<3>& pb,
212                    GeometrySet<3>& output, dim_t<3> )
213 {
214     // everything vs a point
215     if ( pb.handle.which() == PrimitivePoint ) {
216         if ( algorithm::intersects( pa, pb ) ) {
217             output.addPrimitive( *boost::get<const TypeForDimension<3>::Point*>( pb.handle ) );
218         }
219     }
220     else if ( pa.handle.which() == PrimitiveSegment && pb.handle.which() == PrimitiveSegment ) {
221         const CGAL::Segment_3<Kernel>* seg1 = pa.as<CGAL::Segment_3<Kernel> >();
222         const CGAL::Segment_3<Kernel>* seg2 = pb.as<CGAL::Segment_3<Kernel> >();
223         CGAL::Object interObj = CGAL::intersection( *seg1, *seg2 );
224         output.addPrimitive( interObj );
225     }
226     else if ( pa.handle.which() == PrimitiveSurface ) {
227         const CGAL::Triangle_3<Kernel>* tri1 = pa.as<CGAL::Triangle_3<Kernel> >();
228 
229         if ( pb.handle.which() == PrimitiveSegment ) {
230             const CGAL::Segment_3<Kernel>* seg2 = pb.as<CGAL::Segment_3<Kernel> >();
231             CGAL::Object interObj = CGAL::intersection( *tri1, *seg2 );
232             output.addPrimitive( interObj );
233         }
234         else if ( pb.handle.which() == PrimitiveSurface ) {
235             const CGAL::Triangle_3<Kernel>* tri2 = pb.as<CGAL::Triangle_3<Kernel> >();
236             CGAL::Object interObj = CGAL::intersection( *tri1, *tri2 );
237             output.addPrimitive( interObj, /* pointsAsRing */ true );
238         }
239     }
240     else if ( pa.handle.which() == PrimitiveVolume ) {
241         if ( pb.handle.which() == PrimitiveSegment ) {
242             _intersection_solid_segment( pa, pb, output );
243         }
244         else if ( pb.handle.which() == PrimitiveSurface ) {
245             _intersection_solid_triangle( *pa.as<MarkedPolyhedron>(),
246                                           *pb.as<CGAL::Triangle_3<Kernel> >(),
247                                           output );
248         }
249         else if ( pb.handle.which() == PrimitiveVolume ) {
250             const MarkedPolyhedron& sa = *pa.as<MarkedPolyhedron>();
251             const MarkedPolyhedron& sb = *pb.as<MarkedPolyhedron>();
252             BOOST_ASSERT( sa.is_closed() );
253             BOOST_ASSERT( sb.is_closed() );
254             _intersection_solid_solid( sa, sb, output );
255         }
256     }
257 }
258 
259 } // algorithm
260 } // SFCGAL
261