Skip to content

Commit 2d893c4

Browse files
sgiraudotsgiraudot
committed
Partial handling of coplanar cases
1 parent a20fff9 commit 2d893c4

File tree

1 file changed

+144
-64
lines changed

1 file changed

+144
-64
lines changed

‎Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polyline_snapping.h‎

Lines changed: 144 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -65,12 +65,17 @@ struct Snapping_result
6565
template <typename Exact_kernel>
6666
Snapping_result<Exact_kernel>
6767
exact_snapping (const typename Exact_kernel::Segment_3& s0,
68-
const typename Exact_kernel::Segment_3& s1)
68+
const typename Exact_kernel::Segment_3& s1,
69+
const typename Exact_kernel::FT& tolerance)
6970
{
7071
using Point_3 = typename Exact_kernel::Point_3;
72+
using Point_2 = typename Exact_kernel::Point_2;
7173
using Segment_3 = typename Exact_kernel::Segment_3;
74+
using Segment_2 = typename Exact_kernel::Segment_2;
7275
using Line_3 = typename Exact_kernel::Line_3;
76+
using Line_2 = typename Exact_kernel::Line_2;
7377
using Vector_3 = typename Exact_kernel::Vector_3;
78+
using Vector_2 = typename Exact_kernel::Vector_2;
7479
using Plane_3 = typename Exact_kernel::Plane_3;
7580
using FT = typename Exact_kernel::FT;
7681
using Result = Snapping_result<Exact_kernel>;
@@ -85,34 +90,146 @@ exact_snapping (const typename Exact_kernel::Segment_3& s0,
8590
// Collinear segments
8691
if (normal == CGAL::NULL_VECTOR)
8792
{
88-
Point_3 proj0 = v1.supporting_line().projection(s0.source());
89-
Point_3 med = midpoint(s0.source(), proj0);
90-
Line_3 support (med, v0);
91-
Vector ref (med, support.projection(s0.target()));
93+
Line_3 l0 = s0.supporting_line();
94+
Line_3 l1 = s1.supporting_line();
95+
Vector_3 v0 = s0.to_vector();
96+
v0 = v0 / approximate_sqrt(v0 * v0);
97+
Vector_3 v1 = s1.to_vector();
98+
v1 = v1 / approximate_sqrt(v1 * v1);
9299

93-
std::vector<std::pair<FT, Point_3>> points_along_line;
100+
std::vector<std::pair<FT, FT>> points_along_lines;
94101

95102
auto create_projection = [&](const Point_3& p)
96103
{
97-
Point_3 proj = support.projection(p);
98-
Vector vec (med, proj);
99-
FT position = vec * ref;
100-
points_along_line.emplace_back(position, proj);
104+
Point_3 p0 = l0.projection(p);
105+
Vector_3 vec0 (s0.source(), p0);
106+
FT pos0 = vec0 * v0;
107+
108+
Point_3 p1 = l1.projection(p);
109+
Vector_3 vec1 (s1.source(), p1);
110+
FT pos1 = vec1 * v1;
111+
112+
points_along_lines.emplace_back(pos0, pos1);
101113
};
102114
create_projection (s0.source());
103115
create_projection (s0.target());
104116
create_projection (s1.source());
105117
create_projection (s1.target());
106118

107-
std::sort (points_along_line.begin(), points_along_line.end(),
119+
std::sort (points_along_lines.begin(), points_along_lines.end(),
108120
[](const auto& a, const auto& b) -> bool
109121
{
110122
return a.first < b.first;
111123
});
112124

125+
Segment_3 seg0 (s0.source() + points_along_lines[1].first * v0,
126+
s0.source() + points_along_lines[2].first * v0);
127+
Segment_3 seg1 (s1.source() + points_along_lines[1].second * v1,
128+
s1.source() + points_along_lines[2].second * v1);
129+
130+
131+
return Result(seg0, seg1);
132+
}
133+
134+
// Coplanar segments
135+
if (coplanar (s0.source(), s0.target(), s1.source(), s1.target()))
136+
{
137+
Plane_3 plane
138+
(s0.source(), s0.target(),
139+
collinear(s0.source(), s0.target(), s1.target()) ? s1.source() : s1.target());
140+
Vector_3 base1 = plane.base1();
141+
base1 = base1 / approximate_sqrt(base1 * base1);
142+
Vector_3 base2 = plane.base2();
143+
base2 = base2 / approximate_sqrt(base2 * base2);
144+
145+
auto to_2d = [&](const Point_3& p) -> Point_2
146+
{
147+
Vector_3 v (s0.source(), p);
148+
return Point_2 (v * base1, v * base2);
149+
};
150+
auto to_3d = [&](const Point_2& p) -> Point_3
151+
{
152+
return s0.source() + p.x() * base1 + p.y() * base2;
153+
};
154+
155+
auto coplanar_snapping
156+
= [&](const Segment_3& seg, const Segment_3& ref) -> std::pair<Point_3, typename Result::Type>
157+
{
158+
Segment_2 seg2 (to_2d (seg.source()), to_2d (seg.target()));
159+
Segment_2 ref2 (to_2d (ref.source()), to_2d (ref.target()));
160+
161+
auto seg2_i_ref2 = intersection(seg2, ref2);
162+
if (seg2_i_ref2)
163+
{
164+
const Point_2* point = boost::get<Point_2>(&*seg2_i_ref2);
165+
CGAL_assertion (point != nullptr);
166+
return std::make_pair(to_3d(*point), Result::POINT);
167+
}
168+
169+
Vector_2 vec2 = ref2.to_vector();
170+
Vector_2 ortho = vec2.perpendicular(CLOCKWISE);
171+
ortho = ortho / approximate_sqrt(ortho * ortho);
172+
173+
Segment_2 sc (ref2.source() + ortho * tolerance,
174+
ref2.target() + ortho * tolerance);
175+
Segment_2 scc (ref2.source() - ortho * tolerance,
176+
ref2.target() - ortho * tolerance);
177+
178+
Orientation source_sc = orientation (sc.source(), sc.target(), seg2.source());
179+
Orientation source_scc = orientation (scc.source(), scc.target(), seg2.source());
180+
Orientation target_sc = orientation (sc.source(), sc.target(), seg2.target());
181+
Orientation target_scc = orientation (scc.source(), scc.target(), seg2.target());
182+
183+
bool source_in = (source_sc != source_scc);
184+
bool target_in = (target_sc != target_scc);
185+
186+
// Whole segment is inside tolerance
187+
if (source_in && target_in)
188+
return std::make_pair(Point_3(), Result::SEGMENT);
189+
190+
191+
if (!source_in && !target_in)
192+
{
193+
Line_2 lc = sc.supporting_line();
194+
auto inter = intersection (seg2, lc);
195+
// Segment intersect both borders
196+
if (inter)
197+
return std::make_pair(Point_3(), Result::SEGMENT);
198+
199+
// Segment does not intersect at all
200+
return std::make_pair(Point_3(), Result::EMPTY);
201+
}
202+
203+
// One vertex only is in tolerance
204+
if (source_in)
205+
return std::make_pair(seg.source(), Result::POINT);
206+
// else
207+
return std::make_pair(seg.target(), Result::POINT);
208+
};
209+
210+
auto crop0 = coplanar_snapping (s0, s1);
211+
auto crop1 = coplanar_snapping (s1, s0);
212+
213+
if (crop0.second == Result::POINT && crop1.second == Result::POINT)
214+
return Result(crop0.first, crop1.first);
215+
216+
if (crop0.second == Result::POINT && crop1.second == Result::SEGMENT)
217+
return Result (crop0.first, s1.supporting_line().projection(crop0.first));
218+
219+
if (crop0.second == Result::SEGMENT && crop1.second == Result::POINT)
220+
return Result (crop1.first, s0.supporting_line().projection(crop1.first));
221+
222+
if (crop0.second == Result::SEGMENT && crop1.second == Result::SEGMENT)
223+
{
224+
// Note: there might be very specific cases where these segments should be cropped
225+
return Result (s0, s1);
226+
}
227+
113228
return Result();
114229
}
115230

231+
// General case (segments are not coplanar/collinear)
232+
116233
Plane_3 plane0 (s0.source(), normal);
117234
Plane_3 plane1 (s1.source(), normal);
118235

@@ -140,7 +257,6 @@ exact_snapping (const typename Exact_kernel::Segment_3& s0,
140257
return Result (*p0, *p1);
141258
}
142259

143-
// else coplanar segments
144260
return Result();
145261
}
146262

@@ -154,11 +270,13 @@ class Exact_snapping<Kernel,false>
154270
{
155271
//typedefs
156272
using Exact_kernel = Exact_predicates_exact_constructions_kernel;
273+
using Exact_FT = typename Exact_kernel::FT;
157274
using Exact_segment_3 = typename Exact_kernel::Segment_3;
158275
using Exact_point_3 = typename Exact_kernel::Point_3;
159276
using Exact_result = Snapping_result<Exact_kernel>;
160277

161278
using Inexact_kernel = Kernel;
279+
using Inexact_FT = typename Inexact_kernel::FT;
162280
using Inexact_segment_3 = typename Inexact_kernel::Segment_3;
163281
using Inexact_point_3 = typename Inexact_kernel::Point_3;
164282
using Inexact_result = Snapping_result<Inexact_kernel>;
@@ -171,9 +289,10 @@ class Exact_snapping<Kernel,false>
171289

172290
public:
173291

174-
Inexact_result intersection (const Inexact_segment_3& s0, const Inexact_segment_3& s1) const
292+
Inexact_result intersection (const Inexact_segment_3& s0, const Inexact_segment_3& s1,
293+
const Inexact_FT& tolerance) const
175294
{
176-
Exact_result result = exact_snapping<Exact_kernel> (to_exact(s0), to_exact(s1));
295+
Exact_result result = exact_snapping<Exact_kernel> (to_exact(s0), to_exact(s1), to_exact(tolerance));
177296
if (result.type == Exact_result::EMPTY)
178297
return Inexact_result();
179298
if (result.type == Exact_result::POINT)
@@ -189,15 +308,16 @@ class Exact_snapping<Kernel,false>
189308
template <class Kernel>
190309
class Exact_snapping<Kernel,true>
191310
{
311+
using FT = typename Kernel::FT;
192312
using Segment_3 = typename Kernel::Segment_3;
193313
using Point_3 = typename Kernel::Point_3;
194314
using Result = Snapping_result<Kernel>;
195315

196316
public:
197317

198-
Result intersection (const Segment_3& s0, const Segment_3& s1) const
318+
Result intersection (const Segment_3& s0, const Segment_3& s1, const FT& tolerance) const
199319
{
200-
return exact_snapping<Kernel>(s0,s1);
320+
return exact_snapping<Kernel>(s0,s1,tolerance);
201321
}
202322
};
203323

@@ -293,7 +413,7 @@ void polyline_snapping (const PolylineRange& polylines,
293413

294414
Point_3 new_point;
295415
bool okay;
296-
std::tie (new_point, okay) = exact_snapping.intersection (s1, s2);
416+
std::tie (new_point, okay) = exact_snapping.intersection (s1, s2, tolerance);
297417

298418
if (!okay)
299419
return;
@@ -413,6 +533,8 @@ void polyline_snapping (Graph& graph, PointMap point_map, double tolerance)
413533
using Vertex_position = std::pair<FT, vertex_descriptor>;
414534
std::map<edge_descriptor, std::vector<Vertex_position>> new_vertices;
415535

536+
std::map<Point_3, vertex_descriptor> map_p2v;
537+
416538
// Compute intersections
417539
box_self_intersection_d
418540
(boxes.begin(), boxes.end(),
@@ -437,7 +559,7 @@ void polyline_snapping (Graph& graph, PointMap point_map, double tolerance)
437559
return;
438560

439561
// Compute exact snapping point
440-
Result snapping = exact_snapping.intersection (sa, sb);
562+
Result snapping = exact_snapping.intersection (sa, sb, tolerance);
441563
if (snapping.type == Result::EMPTY)
442564
return;
443565

@@ -461,6 +583,10 @@ void polyline_snapping (Graph& graph, PointMap point_map, double tolerance)
461583
ma.first->second.emplace_back (ba, vd);
462584
mb.first->second.emplace_back (bb, vd);
463585
}
586+
else // snapping.type == Result::SEGMENT)
587+
{
588+
std::cerr << "Warning: coplanar segment snapping not handled" << std::endl;
589+
}
464590
},
465591
cutoff);
466592

@@ -485,52 +611,6 @@ void polyline_snapping (Graph& graph, PointMap point_map, double tolerance)
485611
add_edge (vertices[i].second, vertices[i+1].second, graph);
486612
}
487613

488-
#if 0
489-
// Detect edges smaller than threshold
490-
using vedge = std::pair<vertex_descriptor, vertex_descriptor>;
491-
std::vector<vedge> vedges;
492-
for (edge_descriptor ed : make_range(edges(graph).first, edges(graph).second))
493-
{
494-
vertex_descriptor vs = source (ed, graph);
495-
vertex_descriptor vt = target (ed, graph);
496-
vedges.emplace_back (vs, vt);
497-
}
498-
499-
// Map deleted vertex -> valid vertex
500-
std::map<vertex_descriptor, vertex_descriptor> map_v2v;
501-
// Map vertex -> barycenter weight
502-
std::map<vertex_descriptor, std::size_t> map_v2b;
503-
for (vertex_descriptor vd : make_range(vertices(graph).first, vertices(graph).second))
504-
{
505-
map_v2v.insert (std::make_pair(vd, vd));
506-
map_v2b.insert (std::make_pair(vd, 1));
507-
}
508-
509-
// Collapse edges
510-
for (vedge ve : vedges)
511-
{
512-
vertex_descriptor vs = map_v2v[ve.first];
513-
vertex_descriptor vt = map_v2v[ve.second];
514-
515-
if (vs == vt)
516-
continue;
517-
518-
const Point_3& ps = get (point_map, vs);
519-
const Point_3& pt = get (point_map, vt);
520-
521-
if (squared_distance(ps, ps) > tolerance * tolerance)
522-
continue;
523-
524-
std::size_t& bs = map_v2b[vs];
525-
std::size_t& bt = map_v2b[vt];
526-
Point_3 bary = barycenter (ps, bs, pt, bt);
527-
bs += bt;
528-
529-
put (point_map, vs, bary);
530-
531-
map_v2v[vt] = vs;
532-
}
533-
#endif
534614

535615
}
536616

0 commit comments

Comments
 (0)