1 // Copyright (c) 2009 GeometryFactory (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/Triangulation_2/include/CGAL/Triangulation_2/internal/Triangulation_2_projection_traits_base_3.h $ 7 // $Id: Triangulation_2_projection_traits_base_3.h 1a040c8 2021-02-19T07:59:59+01:00 Laurent Rineau 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Laurent Rineau 12 13 14 #ifndef CGAL_INTERNAL_TRIANGULATION_2_PROJECTION_TRAITS_BASE_3_H 15 #define CGAL_INTERNAL_TRIANGULATION_2_PROJECTION_TRAITS_BASE_3_H 16 17 #include <CGAL/license/Triangulation_2.h> 18 19 20 #include <CGAL/Profile_timer.h> 21 #include <CGAL/intersections.h> 22 #include <CGAL/predicates/sign_of_determinant.h> 23 24 namespace CGAL { 25 26 namespace TriangulationProjectionTraitsCartesianFunctors { 27 28 template <class Traits> 29 class Projected_orientation_with_normal_3 30 { 31 // private members 32 typename Traits::Vector_3 normal; 33 34 // private type aliases 35 typedef typename Traits::K K; 36 typedef typename Traits::Point_2 Point; 37 typedef typename Traits::Vector_3 Vector_3; 38 public: 39 typedef typename K::Orientation Orientation; 40 typedef Orientation result_type; 41 Projected_orientation_with_normal_3(const Vector_3 & normal_)42 Projected_orientation_with_normal_3(const Vector_3& normal_) 43 : normal(normal_) 44 { 45 CGAL_PROFILER("Construct Projected_orientation_with_normal_3.") 46 CGAL_TIME_PROFILER("Construct Projected_orientation_with_normal_3") 47 } 48 operator()49 Orientation operator()(const Point& p, 50 const Point& q, 51 const Point& r) const 52 { 53 CGAL_PROFILER("Projected_orientation_with_normal_3::operator()"); 54 CGAL_TIME_PROFILER("Projected_orientation_with_normal_3::operator()"); 55 return orientation(q-p, r-p, normal); 56 } 57 }; // end class Projected_orientation_with_normal_3<Traits> 58 59 template <class Traits> 60 class Projected_side_of_oriented_circle_with_normal_3 61 { 62 // private members 63 typename Traits::Vector_3 normal; 64 65 // private types aliases 66 typedef typename Traits::K K; 67 typedef typename Traits::Point_2 Point; 68 typedef typename Traits::Vector_3 Vector_3; 69 typedef typename Traits::FT FT; 70 71 typedef Projected_side_of_oriented_circle_with_normal_3<Traits> Self; 72 73 public: 74 typedef typename K::Oriented_side Oriented_side; 75 typedef Oriented_side result_type; 76 Projected_side_of_oriented_circle_with_normal_3(const Vector_3 & normal_)77 Projected_side_of_oriented_circle_with_normal_3(const Vector_3& normal_) 78 : normal(normal_) 79 { 80 CGAL_PROFILER("Construct Projected_side_of_oriented_circle_with_normal_3.") 81 CGAL_TIME_PROFILER("Construct Projected_side_of_oriented_circle_with_normal_3.") 82 } 83 operator()84 Oriented_side operator()(const Point& p, 85 const Point& q, 86 const Point& r, 87 const Point& t) const 88 { 89 CGAL_PROFILER("Projected_side_of_oriented_circle_with_normal_3::operator()") 90 CGAL_TIME_PROFILER("Projected_side_of_oriented_circle_with_normal_3::operator()") 91 const Vector_3& u = normal; 92 93 const Vector_3 tp = p - t; 94 const Vector_3 tq = q - t; 95 const Vector_3 tr = r - t; 96 97 const FT tp2 = tp * tp; 98 const FT tq2 = tq * tq; 99 const FT tr2 = tr * tr; 100 const FT u2 = u * u; 101 102 const FT k_p = tp * u; 103 const FT k_q = tq * u; 104 const FT k_r = tr * u; 105 106 return sign_of_determinant<FT>( 107 tp.x(), tp.y(), tp.z(), (tp2 + k_p) * u2 - k_p * k_p, 108 tr.x(), tr.y(), tr.z(), (tr2 + k_r) * u2 - k_r * k_r, 109 tq.x(), tq.y(), tq.z(), (tq2 + k_q) * u2 - k_q * k_q, 110 u.x(), u.y(), u.z(), u2 * u2); 111 // Note that q and r have been swapped in the determinant above, to 112 // inverse its sign. 113 } 114 }; // end class Projected_side_of_oriented_circle_with_normal_3 115 116 template <class Traits> 117 class Projected_squared_distance_with_normal_3 118 { 119 // private members 120 typename Traits::Vector_3 normal; 121 122 // private types aliases 123 typedef typename Traits::K K; 124 typedef typename Traits::Point_2 Point; 125 typedef typename Traits::Line_2 Line; 126 typedef typename Traits::Vector_3 Vector_3; 127 typedef typename Traits::FT FT; 128 129 public: Projected_squared_distance_with_normal_3(const Vector_3 & normal_)130 Projected_squared_distance_with_normal_3(const Vector_3& normal_) 131 : normal(normal_) 132 { 133 CGAL_PROFILER("Construct Projected_squared_distance_with_normal_3.") 134 CGAL_TIME_PROFILER("Construct Projected_squared_distance_with_normal_3") 135 } 136 operator()137 FT operator()(const Point& p, const Point& q) 138 { 139 return squared_distance(p, q); 140 } 141 operator()142 FT operator()(const Line& line, const Point& p) 143 { 144 CGAL_PROFILER("Projected_squared_distance_with_normal_3::operator()") 145 CGAL_TIME_PROFILER("Projected_squared_distance_with_normal_3::operator()") 146 const Vector_3& vl = line.to_vector(); 147 const Point& pl = line.point(); 148 const Vector_3 v = cross_product(normal, 149 vl); 150 if(v == NULL_VECTOR) { 151 // den == 0 if the line is vertical 152 // In that case, the distance is the distance to the line 153 const Vector_3 w = cross_product(pl - p, 154 vl); 155 return (w * w) / (vl * vl); 156 } 157 else { 158 const FT det = determinant(normal, 159 vl, 160 pl - p); 161 return (det * det) / ( v * v ); 162 } 163 } 164 }; // end class Projected_squared_distance_with_normal_3 165 166 template <class Traits> 167 class Projected_intersect_3 168 { 169 // private members 170 typename Traits::Vector_3 normal; 171 172 // private types aliases 173 typedef typename Traits::K K; 174 typedef typename Traits::Point_2 Point; 175 typedef typename Traits::Line_2 Line; 176 typedef typename Traits::Segment_2 Segment; 177 typedef typename K::Plane_3 Plane_3; 178 typedef typename Traits::Vector_3 Vector_3; 179 typedef typename Traits::FT FT; 180 public: Projected_intersect_3(const Vector_3 & normal_)181 Projected_intersect_3(const Vector_3& normal_) 182 : normal(normal_) 183 { 184 CGAL_PROFILER("Construct Projected_intersect_3") 185 CGAL_TIME_PROFILER("Construct Projected_intersect_3") 186 } 187 188 boost::optional<boost::variant<Point,Segment> > operator()189 operator()(const Segment& s1, const Segment& s2) 190 { 191 CGAL_PROFILER("Projected_intersect_3::operator()") 192 CGAL_TIME_PROFILER("Projected_intersect_3::operator()") 193 const Vector_3 u1 = cross_product(s1.to_vector(), normal); 194 if(u1 == NULL_VECTOR) 195 return K().intersect_3_object()(s1.supporting_line(), s2); 196 197 const Vector_3 u2 = cross_product(s2.to_vector(), normal); 198 if(u2 == NULL_VECTOR) 199 return K().intersect_3_object()(s1, s2.supporting_line()); 200 201 const Plane_3 plane_1(s1.source(), u1); 202 const Plane_3 plane_2(s2.source(), u2); 203 204 auto planes_intersection = intersection(plane_1, plane_2); 205 if(! planes_intersection) { 206 std::cerr << "planes_intersection is empty\n"; 207 return boost::none; 208 } 209 if(const Line* line = boost::get<Line>(&*planes_intersection)) 210 { 211 const Point& pi = line->point(0); 212 if(cross_product(normal, pi - s1.source()) 213 * cross_product(normal, pi - s1.target()) > FT(0) 214 || 215 cross_product(normal, pi - s2.source()) 216 * cross_product(normal, pi - s2.target()) > FT(0) ) 217 { 218 // the intersection of the lines is not inside the segments 219 std::cerr << "intersection not inside\n"; 220 return boost::none; 221 } 222 else 223 { 224 // Let the plane passing through s1.source() and with normal 225 // the cross product of s1.to_vector() and s2.to_vector(). That 226 // plane should intersect *l, now. 227 auto inter = intersection(*line, Plane_3(s1.source(), 228 cross_product(s1.to_vector(), 229 s2.to_vector()))); 230 if(! inter){ 231 return boost::none; 232 } 233 if(const Point* point = boost::get<Point>(&*inter)){ 234 typedef boost::variant<Point, Segment> variant_type; 235 return boost::make_optional(variant_type(*point)); 236 } 237 } 238 } 239 if(boost::get<Plane_3>(&*planes_intersection)) 240 { 241 std::cerr << "coplanar lines\n"; 242 CGAL_error(); 243 return boost::none; 244 } 245 return boost::none; 246 } 247 }; // end class Projected_intersect_3 248 249 250 template <class Traits> 251 class Less_along_axis 252 { 253 // private members 254 typedef typename Traits::Vector_3 Vector_3; 255 typedef typename Traits::Point_2 Point; 256 Vector_3 base; 257 public: Less_along_axis(const Vector_3 & base)258 Less_along_axis(const Vector_3& base) : base(base) 259 { 260 CGAL_PROFILER("Construct Less_along_axis") 261 CGAL_TIME_PROFILER("Construct Less_along_axis") 262 } 263 264 typedef bool result_type; 265 operator()266 bool operator() (const Point &p, const Point &q) const { 267 return base * (p - q) < 0; 268 } 269 }; // end class Less_along_axis 270 271 template <class Traits> 272 class Compare_along_axis 273 { 274 // private members 275 typedef typename Traits::Vector_3 Vector_3; 276 typedef typename Traits::Point_2 Point; 277 Vector_3 base; 278 public: Compare_along_axis(const Vector_3 & base)279 Compare_along_axis(const Vector_3& base) : base(base) 280 { 281 CGAL_PROFILER("Construct Compare_along_axis") 282 CGAL_TIME_PROFILER("Construct Compare_along_axis") 283 } 284 285 typedef Comparison_result result_type; 286 operator()287 Comparison_result operator() (const Point &p, const Point &q) const { 288 return compare(base * (p - q), 0); 289 } 290 }; // end class Compare_along_axis 291 292 template <class Traits> 293 class Less_xy_along_axis 294 { 295 // private members 296 typedef typename Traits::Vector_3 Vector_3; 297 typedef typename Traits::Point_2 Point; 298 Vector_3 base1, base2; 299 public: Less_xy_along_axis(const Vector_3 & base1,const Vector_3 & base2)300 Less_xy_along_axis(const Vector_3& base1, const Vector_3& base2) : base1(base1), base2(base2) 301 { 302 CGAL_PROFILER("Construct Less_xy_along_axis") 303 CGAL_TIME_PROFILER("Construct Less_xy_along_axis") 304 } 305 306 typedef bool result_type; 307 operator()308 bool operator() (const Point &p, const Point &q) const { 309 310 Compare_along_axis<Traits> cx(base1); 311 Comparison_result crx = cx(p, q); 312 if (crx == SMALLER) { return true; } 313 if (crx == LARGER) { return false; } 314 Less_along_axis<Traits> ly(base2); 315 return ly(p, q); 316 } 317 }; // end class Less_xy_along_axis 318 319 } // end namespace TriangulationProjectionTraitsCartesianFunctors 320 321 322 template < class Kernel > 323 class Triangulation_2_projection_traits_base_3 324 { 325 typedef Triangulation_2_projection_traits_base_3<Kernel> Self; 326 327 typename Kernel::Vector_3 n, b1, b2; 328 329 public: 330 typedef typename Kernel::Vector_3 Vector_3; 331 332 Triangulation_2_projection_traits_base_3(const Vector_3 & n_)333 explicit Triangulation_2_projection_traits_base_3(const Vector_3& n_) 334 : n(n_) 335 { 336 typedef typename Kernel::FT FT; 337 typedef typename Kernel::Vector_3 Vector_3; 338 339 const FT& nx = n.x(); 340 const FT& ny = n.y(); 341 const FT& nz = n.z(); 342 if(CGAL::abs(nz) >= CGAL::abs(ny)) { 343 b1 = Vector_3(nz, 0, -nx); 344 } 345 else { 346 b1 = Vector_3(ny, -nx, 0); 347 } 348 b2 = cross_product(n, b1); 349 } 350 normal()351 const Vector_3& normal() const 352 { 353 return n; 354 } 355 base1()356 const Vector_3& base1() const{ 357 return b1; 358 } 359 base2()360 const Vector_3& base2() const{ 361 return b2; 362 } 363 364 typedef Kernel K; 365 typedef typename K::FT FT; 366 typedef typename K::Point_3 Point_2; 367 typedef typename K::Segment_3 Segment_2; 368 typedef typename K::Vector_3 Vector_2; 369 typedef typename K::Triangle_3 Triangle_2; 370 typedef typename K::Line_3 Line_2; 371 372 typedef typename K::Angle_3 Angle_2; 373 374 typedef TriangulationProjectionTraitsCartesianFunctors:: 375 Compare_along_axis<Self> Compare_x_2; 376 typedef TriangulationProjectionTraitsCartesianFunctors:: 377 Compare_along_axis<Self> Compare_y_2; 378 379 typedef TriangulationProjectionTraitsCartesianFunctors:: 380 Less_along_axis<Self> Less_x_2; 381 typedef TriangulationProjectionTraitsCartesianFunctors:: 382 Less_along_axis<Self> Less_y_2; 383 typedef TriangulationProjectionTraitsCartesianFunctors:: 384 Less_xy_along_axis<Self> Less_xy_2; 385 386 typedef TriangulationProjectionTraitsCartesianFunctors:: 387 Projected_orientation_with_normal_3<Self> Orientation_2; 388 389 typedef TriangulationProjectionTraitsCartesianFunctors:: 390 Projected_side_of_oriented_circle_with_normal_3<Self> Side_of_oriented_circle_2; 391 392 typedef TriangulationProjectionTraitsCartesianFunctors:: 393 Projected_squared_distance_with_normal_3<Self> Compute_squared_distance_2; 394 395 typedef TriangulationProjectionTraitsCartesianFunctors:: 396 Projected_intersect_3<Self> Intersect_2; 397 398 typedef typename K::Construct_point_3 Construct_point_2; 399 typedef typename K::Construct_weighted_point_3 Construct_weighted_point_2; 400 typedef typename K::Construct_segment_3 Construct_segment_2; 401 typedef typename K::Construct_vector_3 Construct_vector_2; 402 typedef typename K::Construct_line_3 Construct_line_2; 403 typedef typename K::Construct_triangle_3 Construct_triangle_2; 404 405 typedef typename K::Construct_scaled_vector_3 Construct_scaled_vector_2; 406 typedef typename K::Construct_translated_point_3 Construct_translated_point_2; 407 typedef typename K::Construct_midpoint_3 Construct_midpoint_2; 408 typedef typename K::Construct_circumcenter_3 Construct_circumcenter_2; 409 410 typedef typename K::Compute_area_3 Compute_area_2; 411 typedef typename K::Construct_bbox_3 Construct_bbox_2; 412 413 Less_x_2 less_x_2_object()414 less_x_2_object() const 415 { 416 return Less_x_2(this->base1()); 417 } 418 419 Less_y_2 less_y_2_object()420 less_y_2_object() const 421 { 422 return Less_y_2(this->base2()); 423 } 424 425 Less_xy_2 less_xy_2_object()426 less_xy_2_object() const 427 { 428 return Less_xy_2(this->base1(), this->base2()); 429 } 430 431 Compare_x_2 compare_x_2_object()432 compare_x_2_object() const 433 { 434 return Compare_x_2(this->base1()); 435 } 436 437 Compare_y_2 compare_y_2_object()438 compare_y_2_object() const 439 { 440 return Compare_y_2(this->base2()); 441 } 442 443 Orientation_2 orientation_2_object()444 orientation_2_object() const 445 { 446 return Orientation_2(this->normal()); 447 } 448 449 Side_of_oriented_circle_2 side_of_oriented_circle_2_object()450 side_of_oriented_circle_2_object() const 451 { 452 return Side_of_oriented_circle_2(this->normal()); 453 } 454 455 Compute_squared_distance_2 compute_squared_distance_2_object()456 compute_squared_distance_2_object() const 457 { 458 return Compute_squared_distance_2(this->normal()); 459 } 460 461 Intersect_2 intersect_2_object()462 intersect_2_object () const 463 { 464 return Intersect_2(this->normal()); 465 } 466 angle_2_object()467 Angle_2 angle_2_object() const 468 {return Angle_2();} 469 construct_point_2_object()470 Construct_point_2 construct_point_2_object() const 471 {return Construct_point_2();} 472 construct_weighted_point_2_object()473 Construct_weighted_point_2 construct_weighted_point_2_object() const 474 {return Construct_weighted_point_2();} 475 construct_segment_2_object()476 Construct_segment_2 construct_segment_2_object() const 477 {return Construct_segment_2();} 478 construct_vector_2_object()479 Construct_vector_2 construct_vector_2_object() const 480 {return Construct_vector_2();} 481 construct_scaled_vector_2_object()482 Construct_scaled_vector_2 construct_scaled_vector_2_object() const 483 {return Construct_scaled_vector_2();} 484 construct_midpoint_2_object()485 Construct_midpoint_2 construct_midpoint_2_object() const 486 {return Construct_midpoint_2();} 487 construct_circumcenter_2_object()488 Construct_circumcenter_2 construct_circumcenter_2_object() const 489 {return Construct_circumcenter_2();} 490 construct_translated_point_2_object()491 Construct_translated_point_2 construct_translated_point_2_object() const 492 {return Construct_translated_point_2();} 493 construct_line_2_object()494 Construct_line_2 construct_line_2_object() const 495 {return Construct_line_2();} 496 construct_triangle_2_object()497 Construct_triangle_2 construct_triangle_2_object() const 498 {return Construct_triangle_2();} 499 compute_area_2_object()500 Compute_area_2 compute_area_2_object() const 501 {return Compute_area_2();} 502 503 construct_bbox_2_object()504 Construct_bbox_2 construct_bbox_2_object() const 505 {return Construct_bbox_2();} 506 507 508 // Special functor, not in the Kernel concept 509 class Projection_to_plan { 510 // Remeber: Point_2 is K::Point_3 511 const Point_2& plane_point; 512 const Vector_3& normal; 513 public: 514 // Return the projection of a point to a plane passing through 515 // the point 'plane_point' and with orthogonal vector normal(). Projection_to_plan(const Point_2 & plane_point_,const Self & self)516 Projection_to_plan(const Point_2& plane_point_, const Self& self) 517 : plane_point(plane_point_), 518 normal(self.normal()) 519 {} 520 operator()521 Point_2 operator()(const Point_2& point) const 522 { 523 return point + 524 ( ( (plane_point - point) * normal ) / (normal * normal) ) * normal; 525 } 526 }; // end Projection_to_plan 527 projection_to_plan_object(const Point_2 & plane_point)528 Projection_to_plan projection_to_plan_object(const Point_2& plane_point) const 529 { 530 return Projection_to_plan(plane_point, *this); 531 } 532 533 }; // end class Triangulation_2_projection_traits_base_3<Kernel> 534 535 } // end namespace CGAL 536 537 #endif // CGAL_INTERNAL_TRIANGULATION_2_PROJECTION_TRAITS_BASE_3_H 538