1 // Copyright (c) 2015,2016 GeometryFactory
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/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h $
7 // $Id: polylines_to_protect.h b69f643 2021-06-04T15:58:26+02:00 Jane Tournois
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Laurent Rineau
12
13 #ifndef CGAL_MESH_3_POLYLINES_TO_PROTECT_H
14 #define CGAL_MESH_3_POLYLINES_TO_PROTECT_H
15
16 #include <CGAL/license/Mesh_3.h>
17
18
19 #include <vector>
20 #include <map>
21 #include <utility> // std::swap
22 #include <CGAL/tuple.h>
23 #include <CGAL/Image_3.h>
24 #include <CGAL/boost/graph/split_graph_into_polylines.h>
25 #include <CGAL/internal/Mesh_3/Graph_manipulations.h>
26 #include <boost/graph/adjacency_list.hpp>
27 #include <CGAL/Labeled_mesh_domain_3.h> // for CGAL::Null_subdomain_index
28 #include <boost/utility.hpp> // for boost::prior
29 #include <boost/optional.hpp>
30
31 #include <CGAL/Search_traits_3.h>
32 #include <CGAL/Orthogonal_incremental_neighbor_search.h>
33
34 namespace CGAL {
35 namespace Mesh_3 {
36 namespace internal {
37
38 template <typename K, typename NT>
39 struct Returns_midpoint {
40 typedef typename K::Point_3 Point_3;
41
operatorReturns_midpoint42 Point_3 operator()(const Point_3& a,
43 const Point_3& b,
44 const NT,
45 const NT) const
46 {
47 typename K::Construct_midpoint_3 midpt
48 = K().construct_midpoint_3_object();
49 return a < b ? midpt(a, b) : midpt(b, a);
50 }
operatorReturns_midpoint51 double operator()(const NT,
52 const NT) const
53 {
54 return 0.5;
55 }
56 };
57
58 template <typename K, typename NT>
59 struct Linear_interpolator {
60 typedef typename K::Point_3 Point_3;
61
interpolateLinear_interpolator62 Point_3 interpolate(const Point_3& pa,
63 const Point_3& pb,
64 const NT a,
65 const NT b) const
66 {
67 return pa + (a / (a - b)) * ( pb - pa);
68 }
69
70
operatorLinear_interpolator71 Point_3 operator()(const Point_3& pa,
72 const Point_3& pb,
73 const NT a,
74 const NT b) const
75 {
76 return
77 (a < b) ?
78 interpolate(pa, pb, a, b) :
79 interpolate(pb, pa, b, a);
80 }
operatorLinear_interpolator81 double operator()(const NT a,
82 const NT b) const
83 {
84 return a / ( a - b );
85 }
86 };
87
88
89 class Isoline_equation {
90 double a;
91 double b;
92 double c;
93 double d;
94
95 public:
Isoline_equation(double a,double b,double c,double d)96 Isoline_equation(double a, double b, double c, double d)
97 : a(a)
98 , b(b)
99 , c(c)
100 , d(d)
101 {}
102
y(double x)103 double y(double x) const {
104 return ( x*(a-b)-a ) /
105 ( c-a+x*(a-b-c+d) );
106 // If (a+d) == (b+c), then the curve is a straight line.
107 }
x(double y)108 double x(double y) const {
109 return ( y*(a-c)-a ) /
110 ( b-a+y*(a-b-c+d) );
111 // If (a+d) == (b+c), then the curve is a straight line.
112 }
113 }; // end of class Isoline_equation
114
115 template <typename G_manip, typename vertex_descriptor, typename Geom_traits>
116 class Insert_curve_in_graph {
117 typedef typename Geom_traits::Point_3 Point;
118 typedef typename Geom_traits::Vector_3 Vector;
119
120 G_manip& g_manip;
121 const vertex_descriptor null_vertex;
122 const Geom_traits& gt;
123 const typename Geom_traits::Construct_translated_point_3 translate;
124 const typename Geom_traits::Construct_scaled_vector_3 scale;
125 const int prec;
126 const double max_squared_distance;
127
128 public:
129 Insert_curve_in_graph(G_manip& g_manip,
130 const Geom_traits& gt,
131 const int prec = 10,
132 const double max_squared_distance = 0)
g_manip(g_manip)133 : g_manip(g_manip)
134 , null_vertex(g_manip.null_vertex())
135 , gt(gt)
136 , translate(gt.construct_translated_point_3_object())
137 , scale(gt.construct_scaled_vector_3_object())
138 , prec(prec)
139 , max_squared_distance(max_squared_distance)
140 {}
141
142 struct Coords {
CoordsCoords143 Coords(double x, double y) : x(x), y(y) {}
144 double x;
145 double y;
146 };
147
insert_curve(Isoline_equation equation,vertex_descriptor start_v,vertex_descriptor end_v,Coords start,Coords end,Point p00,Vector vx,Vector vy)148 void insert_curve(Isoline_equation equation,
149 vertex_descriptor start_v,
150 vertex_descriptor end_v,
151 Coords start,
152 Coords end,
153 Point p00,
154 Vector vx,
155 Vector vy)
156 {
157 if(CGAL::abs(start.x - end.x) >= CGAL::abs(start.y - end.y)) {
158 insert_curve(equation, &Isoline_equation::y,
159 start_v, end_v,
160 start.x, end.x, p00, vx, vy);
161 } else {
162 insert_curve(equation, &Isoline_equation::x,
163 start_v, end_v,
164 start.y, end.y, p00, vy, vx);
165 }
166 }
167 private:
insert_curve(Isoline_equation equation,double (Isoline_equation::* f)(double)const,vertex_descriptor start_v,vertex_descriptor end_v,double start,double end,Point p00,Vector vx,Vector vy)168 void insert_curve(Isoline_equation equation,
169 double (Isoline_equation::* f)(double) const,
170 vertex_descriptor start_v,
171 vertex_descriptor end_v,
172 double start,
173 double end,
174 Point p00,
175 Vector vx,
176 Vector vy)
177 {
178 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
179 std::cerr << "New curve:\n"
180 << " base("
181 << p00 << " , "
182 << p00+vx << " , "
183 << p00+vy << ")\n"
184 << " vectors: "
185 << "( " << vx << " ) "
186 << " ( " << vy << " )\n";
187 #endif // CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
188 const double step = (end - start)/prec;
189 const double stop = end-step/2;
190 const bool step_is_positive = (step > 0);
191 vertex_descriptor old = start_v;
192 vertex_descriptor v_int = null_vertex;
193 for(double x = start + step;
194 (step_is_positive ? x < stop : x > stop);
195 x += step)
196 {
197 const double y = (equation.*f)(x);
198 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
199 std::cerr << " (" << x << ", " << y << ") -> ";
200 #endif // CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
201 const Point inter_p =
202 translate(translate(p00,
203 scale(vx, x)),
204 scale(vy, y));
205 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
206 std::cerr << "( " << inter_p << ")\n";
207 #endif // CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
208 v_int = g_manip.get_vertex(inter_p, false);
209 g_manip.try_add_edge(old, v_int);
210
211 CGAL_assertion_msg(max_squared_distance == 0 ||
212 CGAL::squared_distance(g_manip.g[old].point,
213 g_manip.g[v_int].point) <
214 max_squared_distance,
215 ([this, old, v_int] {
216 std::stringstream s;
217 s << "Problem at segment ("
218 << this->g_manip.g[old].point << ", "
219 << this->g_manip.g[v_int].point << ")";
220 return s.str();
221 }().c_str()));
222 old = v_int;
223 }
224 if(null_vertex != v_int) {
225 // v_int can be null if the curve is degenerated into one point.
226 g_manip.try_add_edge(v_int, end_v);
227 CGAL_assertion_msg(max_squared_distance == 0 ||
228 CGAL::squared_distance(g_manip.g[end_v].point,
229 g_manip.g[v_int].point) <
230 max_squared_distance,
231 ([this, end_v, v_int] {
232 std::stringstream s;
233 s << "Problem at segment ("
234 << this->g_manip.g[end_v].point << ", "
235 << this->g_manip.g[v_int].point << ")";
236 return s.str();
237 }().c_str()));
238 }
239 }
240 };
241
242 template <typename Pixel,
243 typename Point,
244 typename Domain_type,
245 typename Image_word_type>
246 struct Enriched_pixel {
247 Pixel pixel;
248 Point point;
249 Domain_type domain;
250 Image_word_type word;
251 bool on_edge_of_the_cube;
252 bool on_corner_of_the_cube;
253 }; // end struct template Enriched_pixel<Pix,P,D,C>
254
255 } // end namespace internal
256
257 template<typename P, typename G>
258 struct Polyline_visitor
259 {
260 std::vector<std::vector<P> >& polylines;
261 G& graph;
Polyline_visitorPolyline_visitor262 Polyline_visitor(typename std::vector<std::vector<P> >& lines, G& p_graph)
263 : polylines(lines), graph(p_graph)
264 {}
265
~Polyline_visitorPolyline_visitor266 ~Polyline_visitor()
267 {//DEBUG
268 #if CGAL_MESH_3_PROTECTION_DEBUG & 2
269 std::ofstream og("polylines_graph.polylines.txt");
270 og.precision(17);
271 for(const std::vector<P>& poly : polylines)
272 {
273 og << poly.size() << " ";
274 for(const P& p : poly)
275 og << p << " ";
276 og << std::endl;
277 }
278 #endif // CGAL_MESH_3_PROTECTION_DEBUG & 2
279 }
280
start_new_polylinePolyline_visitor281 void start_new_polyline()
282 {
283 std::vector<P> V;
284 polylines.push_back(V);
285 }
286
add_nodePolyline_visitor287 void add_node(typename boost::graph_traits<G>::vertex_descriptor vd)
288 {
289 std::vector<P>& polyline = polylines.back();
290 polyline.push_back(graph[vd].point);
291 }
292
end_polylinePolyline_visitor293 void end_polyline()
294 {
295 // ignore degenerated polylines
296 if(polylines.back().size() < 2)
297 polylines.resize(polylines.size() - 1);
298 }
299 };
300
301 template <typename Kernel>
302 struct Angle_tester
303 {
304 template <typename vertex_descriptor, typename Graph>
operatorAngle_tester305 bool operator()(vertex_descriptor& v, const Graph& g) const
306 {
307 typedef typename boost::graph_traits<Graph>::out_edge_iterator out_edge_iterator;
308 if (g[v].force_terminal || out_degree(v, g) != 2)
309 return true;
310 else
311 {
312 out_edge_iterator out_edge_it, out_edges_end;
313 boost::tie(out_edge_it, out_edges_end) = out_edges(v, g);
314
315 vertex_descriptor v1 = target(*out_edge_it++, g);
316 vertex_descriptor v2 = target(*out_edge_it++, g);
317 CGAL_assertion(out_edge_it == out_edges_end);
318
319 const typename Kernel::Point_3& p = g[v].point;
320 const typename Kernel::Point_3& p1 = g[v1].point;
321 const typename Kernel::Point_3& p2 = g[v2].point;
322
323 if(CGAL::angle(p1, p, p2) == CGAL::ACUTE) {
324 // const typename Kernel::Vector_3 e1 = p1 - p;
325 // const typename Kernel::Vector_3 e2 = p2 - p;
326 // std::cerr << "At point " << p << ": the angle is "
327 // << ( std::acos(e1 * e2
328 // / CGAL::sqrt(e1*e1)
329 // / CGAL::sqrt(e2*e2))
330 // * 180 / CGAL_PI ) << std::endl;
331 return true;
332 }
333 }
334 return false;
335 }
336 };
337
338 }//namespace Mesh_3
339
340 // this function is overloaded for when `PolylineInputIterator` is `int`.
341 template <typename K,
342 typename Graph,
343 typename PolylineInputIterator>
snap_graph_vertices(Graph & graph,const double vx,const double vy,const double vz,PolylineInputIterator existing_polylines_begin,PolylineInputIterator existing_polylines_end,K)344 void snap_graph_vertices(Graph& graph,
345 const double vx, const double vy, const double vz,
346 PolylineInputIterator existing_polylines_begin,
347 PolylineInputIterator existing_polylines_end,
348 K)
349 {
350 const double dist_bound = (std::min)(vx,
351 (std::min)(vy, vz)) / 256;
352 const double sq_dist_bound = dist_bound * dist_bound;
353
354 typedef CGAL::Search_traits_3<K> Tree_traits;
355 typedef CGAL::Orthogonal_incremental_neighbor_search<Tree_traits> NN_search;
356 typedef typename NN_search::Tree Tree;
357
358 Tree tree;
359 // insert the extremities of the polylines in the kd-tree
360 for(PolylineInputIterator poly_it = existing_polylines_begin;
361 poly_it != existing_polylines_end; ++poly_it)
362 {
363 if(poly_it->begin() != poly_it->end()) {
364 tree.insert(*poly_it->begin());
365 if(boost::next(poly_it->begin()) != poly_it->end()) {
366 tree.insert(*boost::prior(poly_it->end()));
367 }
368 }
369 }
370 if(tree.size() == 0) return;
371
372 for(typename boost::graph_traits<Graph>::vertex_descriptor v :
373 make_range(vertices(graph)))
374 {
375 const typename K::Point_3 p = graph[v].point;
376 NN_search nn(tree, p);
377 CGAL_assertion(nn.begin() != nn.end());
378 if(squared_distance(nn.begin()->first, p) < sq_dist_bound) {
379 graph[v].point = nn.begin()->first;
380 }
381 }
382 }
383
384 template <typename K,
385 typename Graph>
snap_graph_vertices(Graph &,double,double,double,int,int,K)386 void snap_graph_vertices(Graph&, double, double, double, int, int, K)
387 {}
388
389 template <typename Graph>
390 struct Less_for_Graph_vertex_descriptors
391 {
392 const Graph& graph;
Less_for_Graph_vertex_descriptorsLess_for_Graph_vertex_descriptors393 Less_for_Graph_vertex_descriptors(const Graph& graph) : graph(graph) {}
394
395 template <typename vertex_descriptor>
operatorLess_for_Graph_vertex_descriptors396 bool operator()(vertex_descriptor v1, vertex_descriptor v2) const {
397 return graph[v1].point < graph[v2].point;
398 }
399 }; // end of Less_for_Graph_vertex_descriptors<Graph>
400
401 namespace internal {
402 namespace polylines_to_protect_namespace {
403 template <typename Point>
404 struct Vertex_info {
405 Point point;
406 bool force_terminal;
Vertex_infoVertex_info407 Vertex_info() : force_terminal(false) {}
Vertex_infoVertex_info408 Vertex_info(const Point& p) : point(p), force_terminal(false) {}
409 };
410 } // end namespace polylines_to_protect_namespace
411 } // end namespace internal
412
413 template <typename P,
414 typename Image_word_type,
415 typename Null_subdomain_index,
416 typename DomainFunctor,
417 typename InterpolationFunctor,
418 typename PolylineInputIterator>
419 void
420 polylines_to_protect
421 (const CGAL::Image_3& cgal_image,
422 const double vx, const double vy, const double vz,
423 std::vector<std::vector<P> >& polylines,
424 Image_word_type*,
425 Null_subdomain_index null,
426 DomainFunctor domain_fct,
427 InterpolationFunctor interpolate,
428 PolylineInputIterator existing_polylines_begin,
429 PolylineInputIterator existing_polylines_end,
430 boost::optional<Image_word_type> scalar_interpolation_value = boost::none,
431 int prec = 10)
432 {
433 typedef typename DomainFunctor::result_type Domain_type;
434 typedef typename Kernel_traits<P>::Kernel K;
435 typedef P Point_3;
436
437 using CGAL::internal::polylines_to_protect_namespace::Vertex_info;
438 typedef boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS,
439 Vertex_info<Point_3> > Graph;
440
441 typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
442 // typedef typename boost::graph_traits<Graph>::edge_iterator edge_iterator;
443
444 const int xdim = static_cast<int>(cgal_image.xdim());
445 const int ydim = static_cast<int>(cgal_image.ydim());
446 const int zdim = static_cast<int>(cgal_image.zdim());
447 const int image_dims[3] = { xdim, ydim, zdim };
448
449 Graph graph;
450
451 using namespace CGAL::Mesh_3::internal;
452
453 typedef Graph_manipulations<Graph,
454 Point_3,
455 Image_word_type,
456 InterpolationFunctor> G_manip;
457
458 G_manip g_manip(graph, interpolate);
459
460 const float& tx = cgal_image.image()->tx;
461 const float& ty = cgal_image.image()->ty;
462 const float& tz = cgal_image.image()->tz;
463 double max_squared_distance =
464 (std::max)( (std::max)(vx, vy), vz );
465 max_squared_distance *= max_squared_distance;
466 max_squared_distance *= 2;
467
468 typedef Insert_curve_in_graph<G_manip, vertex_descriptor, K> Insert_c_in_g;
469 Insert_c_in_g insert_curve_in_graph(g_manip, K(), prec,
470 max_squared_distance);
471 typedef typename Insert_c_in_g::Coords Coords;
472
473 std::size_t
474 case4 = 0, // 4 colors
475 case211 = 0, case121 = 0, // 3 colors
476 case31 = 0, case22 = 0, // 2 colors
477 case1 = 0; // 1 color
478
479 typename K::Construct_midpoint_3 midpoint =
480 K().construct_midpoint_3_object();
481 typename K::Construct_vector_3 vector =
482 K().construct_vector_3_object();
483 typename K::Construct_translated_point_3 translate =
484 K().construct_translated_point_3_object();
485
486 for(int axis = 0; axis < 3; ++axis)
487 {
488 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
489 std::cerr << "axis = " << axis << "\n";
490 #endif // CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
491 for(int i = 0; i < xdim; i+= (axis == 0 ? (std::max)(1, xdim-1) : 1 ) )
492 for(int j = 0; j < ydim; j+= (axis == 1 ? (std::max)(1, ydim-1) : 1 ) )
493 for(int k = 0; k < zdim; k+= (axis == 2 ? (std::max)(1, zdim-1) : 1 ) )
494 {
495
496 using std::array;
497 using std::tuple;
498 using std::get;
499
500 typedef array<int, 3> Pixel;
501
502 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
503 std::cerr << "Pixel(" << i << ", " << j << ", " << k << ")\n";
504 #endif // CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
505 Pixel pix00 = {{i , j , k }},
506 pix10 = pix00, pix01 = pix00, pix11 = pix00;
507
508 const int axis_xx = (axis + 1) % 3;
509 const int axis_yy = (axis + 2) % 3;
510
511 ++pix10[axis_xx];
512 ++pix11[axis_xx]; ++pix11[axis_yy];
513 ++pix01[axis_yy];
514 if(pix11[0] >= xdim || pix11[1] >= ydim || pix11[2] >= zdim) {
515 // we have gone too far
516 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
517 std::cerr << " ... continue\n";
518 #endif // CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
519 continue;
520 }
521
522 typedef Enriched_pixel<Pixel,
523 Point_3,
524 Domain_type,
525 Image_word_type> Enriched_pixel;
526
527 array<array<Enriched_pixel, 2>, 2> square =
528 {{ {{ { pix00, Point_3(), Domain_type(), 0, false, false },
529 { pix01, Point_3(), Domain_type(), 0, false, false } }},
530 {{ { pix10, Point_3(), Domain_type(), 0, false, false },
531 { pix11, Point_3(), Domain_type(), 0, false, false } }} }};
532
533 std::map<Domain_type, int> pixel_values_set;
534 for(int ii = 0; ii < 2; ++ii) {
535 for(int jj = 0; jj < 2; ++jj) {
536 const Pixel& pixel = square[ii][jj].pixel;
537 double x = pixel[0] * vx + tx;
538 double y = pixel[1] * vy + ty;
539 double z = pixel[2] * vz + tz;
540 short sum_faces = ((0 == pixel[0] || (xdim - 1) == pixel[0]) ? 1 : 0)
541 + ((0 == pixel[1] || (ydim - 1) == pixel[1]) ? 1 : 0)
542 + ((0 == pixel[2] || (zdim - 1) == pixel[2]) ? 1 : 0);
543 square[ii][jj].on_edge_of_the_cube = (sum_faces > 1);
544 square[ii][jj].on_corner_of_the_cube = (sum_faces > 2);
545
546 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
547 if(square[ii][jj].on_edge_of_the_cube) {
548 std::cerr << " Pixel(" << pixel[0] << ", " << pixel[1] << ", "
549 << pixel[2] << ") is on edge\n";
550 }
551 if (square[ii][jj].on_corner_of_the_cube) {
552 std::cerr << " Pixel(" << pixel[0] << ", " << pixel[1] << ", "
553 << pixel[2] << ") is on corner\n";
554 }
555 #endif // CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
556
557 square[ii][jj].point = Point_3(x, y, z);
558 square[ii][jj].word =
559 CGAL::IMAGEIO::static_evaluate<Image_word_type>
560 (cgal_image.image(),
561 pixel[0],
562 pixel[1],
563 pixel[2]);
564 square[ii][jj].domain = domain_fct(square[ii][jj].word);
565 if(scalar_interpolation_value != boost::none) {
566 square[ii][jj].word =
567 Image_word_type(square[ii][jj].word -
568 (*scalar_interpolation_value));
569 }
570 ++pixel_values_set[square[ii][jj].domain];
571 }
572 }
573
574 const Point_3& p00 = square[0][0].point;
575 const Point_3& p10 = square[1][0].point;
576 const Point_3& p01 = square[0][1].point;
577 const Point_3& p11 = square[1][1].point;
578
579 const Image_word_type& v00 = square[0][0].word;
580 const Image_word_type& v10 = square[1][0].word;
581 const Image_word_type& v01 = square[0][1].word;
582 const Image_word_type& v11 = square[1][1].word;
583
584 {
585 bool out00 = null(square[0][0].domain);
586 bool out10 = null(square[1][0].domain);
587 bool out01 = null(square[0][1].domain);
588 bool out11 = null(square[1][1].domain);
589
590 bool is_corner00 = square[0][0].on_corner_of_the_cube;
591 bool is_corner10 = square[1][0].on_corner_of_the_cube;
592 bool is_corner01 = square[0][1].on_corner_of_the_cube;
593 bool is_corner11 = square[1][1].on_corner_of_the_cube;
594
595 //
596 // Protect the edges of the cube
597 //
598 if(pix00[axis_xx] == 0 &&
599 ! ( out00 && out01 ) )
600 {
601 g_manip.try_add_edge(g_manip.get_vertex(p00, is_corner00),
602 g_manip.get_vertex(p01, is_corner01));
603 }
604 if(pix11[axis_xx] == image_dims[axis_xx]-1 &&
605 ! ( out10 && out11 ) )
606 {
607 g_manip.try_add_edge(g_manip.get_vertex(p10, is_corner10),
608 g_manip.get_vertex(p11, is_corner11));
609 }
610 if(pix00[axis_yy] == 0 &&
611 ! ( out00 && out10 ) )
612 {
613 g_manip.try_add_edge(g_manip.get_vertex(p00, is_corner00),
614 g_manip.get_vertex(p10, is_corner10));
615 }
616 if(pix11[axis_yy] == image_dims[axis_yy]-1 &&
617 ! ( out01 && out11 ) )
618 {
619 g_manip.try_add_edge(g_manip.get_vertex(p01, is_corner01),
620 g_manip.get_vertex(p11, is_corner11));
621 }
622 } // end of scope for outIJ and on_edgeIJ
623
624 //
625 // Protect lines inside the square
626 //
627 switch(pixel_values_set.size()) {
628 case 4: {
629 CGAL_assertion(square[0][0].domain != square[0][1].domain);
630 CGAL_assertion(square[0][0].domain != square[1][0].domain);
631 CGAL_assertion(square[0][0].domain != square[1][1].domain);
632 CGAL_assertion(square[1][0].domain != square[1][1].domain);
633 CGAL_assertion(square[0][1].domain != square[1][1].domain);
634 CGAL_assertion(square[0][1].domain != square[1][0].domain);
635
636 case_4:
637 // case 4 or case 2-2
638 //
639 // ---------------
640 // | | |
641 // | | |
642 // |______|______|
643 // | | |
644 // | | |
645 // | | |
646 // ---------------
647 //
648 ++case4;
649 vertex_descriptor left = g_manip.split(square[0][0],
650 square[0][1], null);
651 vertex_descriptor right = g_manip.split(square[1][0],
652 square[1][1], null);
653 vertex_descriptor top = g_manip.split(square[0][1],
654 square[1][1], null);
655 vertex_descriptor bottom = g_manip.split(square[0][0],
656 square[1][0], null);
657
658 vertex_descriptor vmid = g_manip.split(graph[left].point,
659 graph[right].point,
660 v00, v10,
661 false, false,
662 false, false);
663 g_manip.try_add_edge(left , vmid);
664 g_manip.try_add_edge(right , vmid);
665 g_manip.try_add_edge(top , vmid);
666 g_manip.try_add_edge(bottom , vmid);
667 }
668 break;
669 case 3: {
670 if(square[0][0].domain == square[1][1].domain) {
671 // Diagonal, but the wrong one.
672 // Vertical swap
673 std::swap(square[0][1], square[0][0]);
674 std::swap(square[1][1], square[1][0]);
675 }
676 if(square[0][1].domain == square[1][0].domain) {
677 // diagonal case 1-2-1
678 //
679 // A-------------C
680 // | | |
681 // | \ |
682 // |___ \__|
683 // | \ |
684 // | \ |
685 // | | |
686 // B-------------A
687 //
688 // two curves: -1 ------------- 0 -1 ------------- 1
689 // | | | | |
690 // | | | \ |
691 // |___ | | \__|
692 // | \ | | |
693 // | \ | | |
694 // | | | | |
695 // 1 ------------- -1 0 ------------- -1
696 //
697 CGAL_assertion(square[0][1].domain == square[1][0].domain);
698 CGAL_assertion(square[1][1].domain != square[0][0].domain);
699 CGAL_assertion(square[0][1].domain != square[0][0].domain);
700 CGAL_assertion(square[0][1].domain != square[1][1].domain);
701 case_1_2_1:
702 ++case121;
703 vertex_descriptor left = g_manip.split(square[0][0],
704 square[0][1], null);
705 vertex_descriptor right = g_manip.split(square[1][0],
706 square[1][1], null);
707 vertex_descriptor top = g_manip.split(square[0][1],
708 square[1][1], null);
709 vertex_descriptor bottom = g_manip.split(square[0][0],
710 square[1][0], null);
711
712 Isoline_equation equation =
713 (scalar_interpolation_value == boost::none) ?
714 Isoline_equation(1, -1, -1, 0) :
715 Isoline_equation(v00, v10, v01, v11);
716 insert_curve_in_graph.insert_curve(equation,
717 left, bottom,
718 Coords(0., interpolate(v00, v01) ),
719 Coords(interpolate(v00, v10), 0. ),
720 p00,
721 p10 - p00,
722 p01 - p00);
723 if(scalar_interpolation_value == boost::none) {
724 equation = Isoline_equation(0, -1, -1, 1);
725 }
726 insert_curve_in_graph.insert_curve(equation,
727 top, right,
728 Coords(interpolate(v01, v11), 1. ),
729 Coords(1., interpolate(v10, v11) ),
730 p00,
731 p10 - p00,
732 p01 - p00);
733 } else {
734 // case 2-1-1
735 if(square[0][0].domain == square[1][0].domain) {
736 // Diagonal swap
737 std::swap(square[0][1], square[1][0]);
738 } else
739 if(square[0][1].domain == square[1][1].domain) {
740 // The other diagonal swap
741 std::swap(square[0][0], square[1][1]);
742 } else
743 if(square[1][0].domain == square[1][1].domain) {
744 // Vertical swap
745 std::swap(square[0][0], square[1][0]);
746 std::swap(square[0][1], square[1][1]);
747 }
748 CGAL_assertion(square[0][0].domain == square[0][1].domain);
749 CGAL_assertion(square[0][0].domain != square[1][0].domain);
750 CGAL_assertion(square[0][0].domain != square[1][1].domain);
751 CGAL_assertion(square[1][0].domain != square[1][1].domain);
752 ++case211;
753 // Note: this case cannot occur with non-segmented scalar
754 // images, because it needs three domains.
755 //
756 // A-------------B
757 // | | |
758 // | \ |
759 // | \__|
760 // | / |
761 // | / |
762 // | | |
763 // A-------------C
764 //
765 // Two curves portions:
766 //
767 // 1 ------------- 0 1 ------------- -1
768 // | /| | | |
769 // | / | | \ |
770 // | / | | \ |
771 // | / | | \ |
772 // | / | | \ |
773 // | | | | \|
774 // 1 ------------- -1 1 ------------- 0
775 //
776 // ...and a segment.
777 Point_3 midleft = midpoint(p00, p01);
778 Point_3 midright = midpoint(p10, p11);
779 Point_3 inter = translate(midleft
780 , (2./3) * vector(midleft, midright));
781 vertex_descriptor v_inter = g_manip.get_vertex(inter, false);
782 vertex_descriptor right = g_manip.split(square[1][0],
783 square[1][1], null);
784 vertex_descriptor top = g_manip.split(square[0][1],
785 square[1][1], null);
786 vertex_descriptor bottom = g_manip.split(square[0][0],
787 square[1][0], null);
788
789 insert_curve_in_graph.insert_curve(Isoline_equation(1, -1, 1, 0),
790 bottom, v_inter,
791 Coords( 0.5 , 0. ),
792 Coords( 2./3., 0.5 ),
793 p00,
794 p10 - p00,
795 p01 - p00);
796 insert_curve_in_graph.insert_curve(Isoline_equation(1, 0, 1, -1),
797 top, v_inter,
798 Coords( 0.5 , 1. ),
799 Coords( 2./3., 0.5 ),
800 p00,
801 p10 - p00,
802 p01 - p00);
803 g_manip.try_add_edge(right, v_inter);
804 } // end case 2-1-1
805 } // end `case 3:`
806 break;
807 case 2: {
808 if(pixel_values_set.begin()->second ==
809 pixel_values_set.rbegin()->second)
810 {
811 // Case of two colors with two pixels each.
812
813 if(square[0][0].domain==square[1][0].domain) {
814 // case 2-2, diagonal swap
815 std::swap(square[0][1], square[1][0]);
816 CGAL_assertion(square[0][0].domain==square[0][1].domain);
817 }
818 if(square[1][0].domain==square[1][1].domain) {
819 // case 2-2, vertical swap
820 std::swap(square[0][1], square[1][1]);
821 std::swap(square[0][0], square[1][0]);
822 CGAL_assertion(square[0][0].domain==square[0][1].domain);
823 }
824 if(square[0][1].domain==square[1][1].domain) {
825 // case 2-2, diagonal swap
826 std::swap(square[0][0], square[1][1]);
827 CGAL_assertion(square[0][0].domain==square[0][1].domain);
828 }
829
830 if(square[0][0].domain==square[0][1].domain) {
831 // vertical case 2-2
832 ++case22;
833 // case 2-2
834 //
835 // A-------------B
836 // | | |
837 // | | |
838 // | | |
839 // | | |
840 // | | |
841 // | | |
842 // A-------------B
843 //
844 CGAL_assertion(square[1][0].domain==square[1][1].domain);
845 CGAL_assertion(square[1][0].domain!=square[0][1].domain);
846 vertex_descriptor top = g_manip.split(square[0][1],
847 square[1][1], null);
848 vertex_descriptor bottom = g_manip.split(square[0][0],
849 square[1][0], null);
850 if(scalar_interpolation_value == boost::none) {
851 g_manip.try_add_edge(top, bottom);
852 } else {
853 insert_curve_in_graph.insert_curve(Isoline_equation(v00, v10,
854 v01, v11),
855 top, bottom,
856 Coords(interpolate(v01,v11), 1.),
857 Coords(interpolate(v00,v10), 0.),
858 p00,
859 p10 - p00,
860 p01 - p00);
861 }
862 } else {
863 // Else diagonal case case 2-2
864 // Same as the case with 4 colors
865 CGAL_assertion(square[0][0].domain==square[1][1].domain);
866 CGAL_assertion(square[1][0].domain==square[0][1].domain);
867 CGAL_assertion(square[0][0].domain!=square[0][1].domain);
868
869 if(scalar_interpolation_value != boost::none) {
870 // Compute the squared distance between the two branches of
871 // the hyperbola.
872 const double discrimant = double(v00) * v11 - double(v01) * v10;
873 const double squared_distance =
874 8 * CGAL::abs(discrimant) / CGAL::square(double(v00) - v10 - v01 + v11);
875 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
876 std::cerr << "squared_distance: " << squared_distance << "\n";
877 #endif // CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
878 if(CGAL::square(prec) * squared_distance > 1)
879 {
880 // In that case, the case 1-2-1 will be applied
881 if(discrimant > 0.) {
882 // Vertical swap
883 std::swap(square[0][1], square[0][0]);
884 std::swap(square[1][1], square[1][0]);
885 }
886 goto case_1_2_1;
887 }
888 } // scalar interpolation
889
890 goto case_4;
891 }
892 }
893 else {
894 // case of two colors with one pixel green and three red
895 Domain_type value_alone;
896 if(pixel_values_set.begin()->second == 1) {
897 value_alone = pixel_values_set.begin()->first;
898 } else {
899 CGAL_assertion(pixel_values_set.begin()->second == 3);
900 CGAL_assertion(pixel_values_set.rbegin()->second ==1);
901 value_alone = pixel_values_set.rbegin()->first;
902 }
903 if(square[0][1].domain == value_alone) {
904 // central symmetry
905 std::swap(square[0][1], square[1][0]);
906 std::swap(square[0][0], square[1][1]);
907 CGAL_assertion(square[1][0].domain == value_alone);
908 }
909 if(square[1][1].domain == value_alone) {
910 // vertical swap
911 std::swap(square[0][0], square[0][1]);
912 std::swap(square[1][0], square[1][1]);
913 CGAL_assertion(square[1][0].domain == value_alone);
914 }
915 if(square[0][0].domain == value_alone) {
916 // horizontal swap
917 std::swap(square[0][1], square[1][1]);
918 std::swap(square[0][0], square[1][0]);
919 CGAL_assertion(square[1][0].domain == value_alone);
920 }
921 ++case31;
922 //
923 // A-------------A 1 ------------- 1
924 // | | | |
925 // | | | |
926 // | ____| | ____|
927 // | / | | / |
928 // | / | | / |
929 // | | | | | |
930 // A-------------B 1 ------------- -1
931 //
932 CGAL_assertion(square[1][0].domain == value_alone);
933 CGAL_assertion(square[1][0].domain != square[0][0].domain);
934 CGAL_assertion(square[1][1].domain == square[0][0].domain);
935 CGAL_assertion(square[0][1].domain == square[0][0].domain);
936 vertex_descriptor bottom = g_manip.split(square[0][0],
937 square[1][0], null);
938 vertex_descriptor right = g_manip.split(square[1][0],
939 square[1][1], null);
940
941 Isoline_equation equation =
942 (scalar_interpolation_value == boost::none) ?
943 Isoline_equation(1, -1, 1, 1) :
944 Isoline_equation(v00, v10, v01, v11);
945
946 insert_curve_in_graph.insert_curve(equation,
947 bottom, right,
948 Coords(interpolate(v00, v10), 0. ),
949 Coords(1., interpolate(v10, v11) ),
950 p00,
951 p10 - p00,
952 p01 - p00);
953 }
954 }
955 break;
956 default: // case 1
957 ++case1;
958 // nothing to do
959 break;
960 }
961 }
962 }
963 // std::cerr << "case 4: " << case4 << std::endl;
964 // std::cerr << "case 2-1-1: " << case211 << std::endl;
965 // std::cerr << "case 1-2-1: " << case121 << std::endl;
966 // std::cerr << "case 3-1: " << case31 << std::endl;
967 // std::cerr << "case 2-2: " << case22 << std::endl;
968 // std::cerr << "case 1: " << case1 << std::endl;
969
970 const std::ptrdiff_t nb_facets =
971 case4 + case211 + case121 + case31 + case22 + case1;
972 const std::ptrdiff_t expected_nb_facets =
973 ((xdim != 1 && ydim != 1 && zdim != 1) ? 2 : 1)
974 *
975 ((xdim-1)*(ydim-1) + (ydim-1)*(zdim-1) + (xdim-1)*(zdim-1));
976
977 // std::cerr << "nb of facets: " << nb_facets << std::endl
978 // << " expected nb of facets: " << expected_nb_facets << std::endl;
979
980 CGAL_assertion(nb_facets == expected_nb_facets);
981 CGAL_USE(nb_facets); CGAL_USE(expected_nb_facets);
982
983 snap_graph_vertices(graph,
984 vx, vy, vz,
985 existing_polylines_begin, existing_polylines_end,
986 K());
987
988 Mesh_3::Polyline_visitor<Point_3, Graph> visitor(polylines, graph);
989 Less_for_Graph_vertex_descriptors<Graph> less(graph);
990 const Graph& const_graph = graph;
991 split_graph_into_polylines(const_graph, visitor,
992 Mesh_3::Angle_tester<K>(), less);
993 }
994
995 template <typename P,
996 typename PolylineInputIterator>
997 void
polylines_to_protect(std::vector<std::vector<P>> & polylines,PolylineInputIterator existing_polylines_begin,PolylineInputIterator existing_polylines_end)998 polylines_to_protect(std::vector<std::vector<P> >& polylines,
999 PolylineInputIterator existing_polylines_begin,
1000 PolylineInputIterator existing_polylines_end)
1001 {
1002 typedef P Point_3;
1003 typedef typename Kernel_traits<P>::Kernel K;
1004 using CGAL::internal::polylines_to_protect_namespace::Vertex_info;
1005 typedef boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS,
1006 Vertex_info<Point_3> > Graph;
1007 typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
1008 typedef typename std::iterator_traits<PolylineInputIterator>::value_type Polyline;
1009
1010 Graph graph;
1011 typedef Mesh_3::internal::Returns_midpoint<K, int> Midpoint_fct;
1012 Mesh_3::internal::Graph_manipulations<Graph,
1013 Point_3,
1014 int,
1015 Midpoint_fct> g_manip(graph);
1016
1017 for (PolylineInputIterator poly_it = existing_polylines_begin;
1018 poly_it != existing_polylines_end; ++poly_it)
1019 {
1020 Polyline polyline = *poly_it;
1021 if (polyline.size() < 2)
1022 continue;
1023
1024 typename Polyline::iterator pit = polyline.begin();
1025 while (boost::next(pit) != polyline.end())
1026 {
1027 vertex_descriptor v = g_manip.get_vertex(*pit, false);
1028 vertex_descriptor w = g_manip.get_vertex(*boost::next(pit), false);
1029 g_manip.try_add_edge(v, w);
1030 ++pit;
1031 }
1032 }
1033
1034 Mesh_3::Polyline_visitor<Point_3, Graph> visitor(polylines, graph);
1035 Less_for_Graph_vertex_descriptors<Graph> less(graph);
1036 const Graph& const_graph = graph;
1037 typedef typename Kernel_traits<P>::Kernel K;
1038 split_graph_into_polylines(const_graph, visitor,
1039 Mesh_3::Angle_tester<K>(), less);
1040 }
1041
1042 template <typename P, typename Image_word_type, typename Null_subdomain_index>
1043 void
polylines_to_protect(const CGAL::Image_3 & cgal_image,std::vector<std::vector<P>> & polylines,Image_word_type * word_type,Null_subdomain_index null)1044 polylines_to_protect(const CGAL::Image_3& cgal_image,
1045 std::vector<std::vector<P> >& polylines,
1046 Image_word_type* word_type,
1047 Null_subdomain_index null)
1048 {
1049 polylines_to_protect<P>
1050 (cgal_image,
1051 polylines,
1052 word_type,
1053 null,
1054 0,
1055 0);
1056 }
1057
1058 template <typename P, typename Image_word_type>
1059 void
polylines_to_protect(const CGAL::Image_3 & cgal_image,std::vector<std::vector<P>> & polylines)1060 polylines_to_protect(const CGAL::Image_3& cgal_image,
1061 std::vector<std::vector<P> >& polylines)
1062 {
1063 polylines_to_protect<P, Image_word_type>(cgal_image, polylines, 0, 0);
1064 }
1065
1066 template <typename P,
1067 typename Image_word_type,
1068 typename Null_subdomain_index,
1069 typename PolylineInputIterator>
1070 void
polylines_to_protect(const CGAL::Image_3 & cgal_image,std::vector<std::vector<P>> & polylines,Image_word_type * word_type,Null_subdomain_index null,PolylineInputIterator existing_polylines_begin,PolylineInputIterator existing_polylines_end)1071 polylines_to_protect(const CGAL::Image_3& cgal_image,
1072 std::vector<std::vector<P> >& polylines,
1073 Image_word_type* word_type,
1074 Null_subdomain_index null,
1075 PolylineInputIterator existing_polylines_begin,
1076 PolylineInputIterator existing_polylines_end)
1077 {
1078 typedef typename Kernel_traits<P>::Kernel K;
1079 polylines_to_protect<P>
1080 (cgal_image,
1081 cgal_image.vx(), cgal_image.vy(),cgal_image.vz(),
1082 polylines,
1083 word_type,
1084 null,
1085 CGAL::Identity<Image_word_type>(),
1086 Mesh_3::internal::Returns_midpoint<K, Image_word_type>(),
1087 existing_polylines_begin,
1088 existing_polylines_end);
1089 }
1090
1091 template <typename P,
1092 typename Image_word_type,
1093 typename PolylineInputIterator>
1094 void
polylines_to_protect(const CGAL::Image_3 & cgal_image,std::vector<std::vector<P>> & polylines,PolylineInputIterator existing_polylines_begin,PolylineInputIterator existing_polylines_end)1095 polylines_to_protect(const CGAL::Image_3& cgal_image,
1096 std::vector<std::vector<P> >& polylines,
1097 PolylineInputIterator existing_polylines_begin,
1098 PolylineInputIterator existing_polylines_end)
1099 {
1100 polylines_to_protect<P>
1101 (cgal_image,
1102 polylines,
1103 (Image_word_type*)0,
1104 CGAL::Null_subdomain_index(),
1105 existing_polylines_begin,
1106 existing_polylines_end);
1107 }
1108
1109 } // namespace CGAL
1110
1111 #endif // CGAL_MESH_3_POLYLINES_TO_PROTECT_H
1112