1 // Copyright (c) 2006-2007 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/Surface_mesher/include/CGAL/Surface_mesher/Sphere_oracle_3.h $ 7 // $Id: Sphere_oracle_3.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Laurent RINEAU 12 13 #ifndef CGAL_SURFACE_MESHER_SPHERE_ORACLE_3_H 14 #define CGAL_SURFACE_MESHER_SPHERE_ORACLE_3_H 15 16 #include <CGAL/license/Surface_mesher.h> 17 18 #include <CGAL/disable_warnings.h> 19 20 #include <CGAL/Surface_mesher/Null_oracle_visitor.h> 21 #include <CGAL/point_generators_3.h> 22 #include <CGAL/number_utils.h> 23 #include <CGAL/Origin.h> 24 25 #include <boost/tuple/tuple.hpp> 26 27 namespace CGAL { 28 29 namespace Surface_mesher { 30 31 template < 32 class GT, 33 class Point_creator = Creator_uniform_3<typename GT::FT, 34 typename GT::Point_3>, 35 class Visitor = Null_oracle_visitor 36 > 37 class Sphere_oracle_3 38 { 39 // private types 40 typedef Sphere_oracle_3<GT, Point_creator, Visitor> Self; 41 42 typedef typename GT::Point_3 Point; 43 typedef typename GT::FT FT; 44 typedef typename GT::Sphere_3 Sphere_3; 45 46 public: 47 48 // Public types 49 typedef GT Geom_traits; 50 typedef typename GT::Point_3 Point_3; 51 typedef typename GT::Segment_3 Segment_3; 52 typedef typename GT::Ray_3 Ray_3; 53 typedef typename GT::Line_3 Line_3; 54 55 typedef Sphere_3 Surface_3; 56 57 typedef Point Intersection_point; 58 private: 59 // Private members 60 Visitor visitor; // a visitor that can modify a point, before returning it. 61 62 public: 63 64 // Constructors 65 Sphere_oracle_3 (Visitor visitor_ = Visitor() ) : visitor(visitor_)66 visitor(visitor_) 67 { 68 #ifdef CGAL_SURFACE_MESHER_DEBUG_CONSTRUCTORS 69 # ifndef CGAL_SURFACE_MESHER_IMPLICIT_SURFACE_ORACLE_3_H 70 std::cerr << "CONS: Sphere_oracle_3\n"; 71 # endif 72 #endif 73 } 74 get_visitor()75 const Visitor& get_visitor() const 76 { 77 return visitor; 78 } 79 80 // Predicates and Constructions 81 is_in_volume(const Surface_3 & sphere,const Point & p)82 bool is_in_volume(const Surface_3& sphere, const Point& p) const 83 { 84 typename GT::Has_on_bounded_side_3 on_bounded_side_of_sphere = 85 GT().has_on_bounded_side_3_object(); 86 87 return on_bounded_side_of_sphere(sphere, p); 88 } 89 90 class Intersect_3 91 { 92 const Self& oracle; 93 94 boost::tuple<int, FT, FT> intersection_line_sphere_lambda(const Surface_3 & sphere,const Point & a,const Point & b)95 intersection_line_sphere_lambda(const Surface_3& sphere, 96 const Point& a, 97 const Point& b) const 98 { 99 /* 100 Let the vectorial line equation: 101 m = a + lambda * ( b - a ) 102 (a, b, and m are points, and lambda if a real.) 103 104 Let c be the center of the sphere, of radius r. 105 The intersection of the line and the sphere is given by: 106 (c-m)^2 = r^2 107 That is: 108 ((c-a)^2 - r^2) 109 - 2 lambda (c-a)*(b-a) 110 + lambda^2 (b-a)^2 == 0 111 112 (second degre equation) 113 114 deltaprime = delta/4 = ((c-a)(b-a))^2 - (b-a)^2 * ( (c-a)^2 -r^2 ) 115 116 if delta > 0, root_1 = ((c-a)(b-a) - \sqrt(delta/4)) / (b-a)^2 117 root_2 = ((c-a)(b-a) + \sqrt(delta/4)) / (b-a)^2 118 (root_1 < root_2) 119 */ 120 121 typedef typename GT::Vector_3 Vector_3; 122 123 typename GT::Construct_vector_3 vector = 124 GT().construct_vector_3_object(); 125 typename GT::Compute_scalar_product_3 scalar_product = 126 GT().compute_scalar_product_3_object(); 127 typename GT::Compute_squared_distance_3 squared_distance = 128 GT().compute_squared_distance_3_object(); 129 typename GT::Construct_center_3 center = 130 GT().construct_center_3_object(); 131 typename GT::Compute_squared_radius_3 squared_radius = 132 GT().compute_squared_radius_3_object(); 133 134 const Point c = center(sphere); 135 const Vector_3 ab = vector(a, b); 136 const Vector_3 ac = vector(a, c); 137 const FT ab_ac = scalar_product(ab, ac); 138 const FT ab2 = squared_distance(a, b); 139 const FT ac2 = squared_distance(a, c); 140 const FT r2 = squared_radius(sphere); 141 const FT deltaprime = ab_ac * ab_ac - ab2 * ( ac2 - r2 ); 142 143 switch( CGAL::sign(deltaprime) ) 144 { 145 case ZERO: 146 return boost::make_tuple(1, ab_ac / ab2, 0); 147 case POSITIVE: 148 { 149 const FT sqrt_deltaprime = CGAL::sqrt(deltaprime); 150 return boost::make_tuple(2, 151 (ab_ac - sqrt_deltaprime) / ab2, 152 (ab_ac + sqrt_deltaprime) / ab2); 153 } 154 case NEGATIVE: 155 break; 156 } 157 return boost::make_tuple(0, 0, 0); 158 } //end intersection_line_sphere_lambda 159 160 template <class Assert_on_lambda> private_intersection(const Surface_3 & sphere,const Point & a,const Point & b,Assert_on_lambda test)161 Object private_intersection(const Surface_3& sphere, 162 const Point& a, 163 const Point& b, 164 Assert_on_lambda test) const 165 { 166 typedef typename GT::Vector_3 Vector; 167 168 typename GT::Construct_vector_3 vector = 169 GT().construct_vector_3_object(); 170 typename GT::Construct_scaled_vector_3 scaled_vector = 171 GT().construct_scaled_vector_3_object(); 172 typename GT::Construct_translated_point_3 translated_point = 173 GT().construct_translated_point_3_object(); 174 175 int number_of_roots; 176 FT root_1, root_2; 177 boost::tie(number_of_roots, root_1, root_2) = 178 intersection_line_sphere_lambda(sphere, a, b); 179 180 const Vector ab = vector(a, b); 181 if(number_of_roots > 0 && test(root_1)) 182 { 183 Point p = translated_point(a, scaled_vector(ab, root_1)); 184 oracle.get_visitor().new_point(p); 185 return make_object(p); 186 } 187 else if (number_of_roots > 1 && test(root_2)) 188 { 189 Point p = translated_point(a, scaled_vector(ab, root_2)); 190 oracle.get_visitor().new_point(p); 191 return make_object(p); 192 } 193 // else 194 return Object(); 195 } // end private_intersection 196 197 struct Lambda_between_0_and_1 : public CGAL::cpp98::unary_function<FT, bool> 198 { operatorLambda_between_0_and_1199 bool operator()(const FT x) const 200 { 201 return FT(0) <= x && x <= FT(1); 202 } 203 }; 204 205 struct Lambda_positive : public CGAL::cpp98::unary_function<FT, bool> 206 { operatorLambda_positive207 bool operator()(const FT x) const 208 { 209 return FT(0) <= x; 210 } 211 }; 212 213 struct Always_true : public CGAL::cpp98::unary_function<FT, bool> 214 { operatorAlways_true215 bool operator()(const FT) const 216 { 217 return true; 218 } 219 }; 220 221 public: Intersect_3(const Self & oracle)222 Intersect_3(const Self& oracle) : oracle(oracle) 223 { 224 } 225 operator()226 Object operator()(const Surface_3& sphere, Segment_3 s) const 227 { 228 typename GT::Construct_point_on_3 point_on = 229 GT().construct_point_on_3_object(); 230 231 const Point& a = point_on(s, 0); 232 const Point& b = point_on(s, 1); 233 234 return private_intersection(sphere, a, b, Lambda_between_0_and_1()); 235 } // end operator()(Surface_3, Segment_3) 236 operator()237 Object operator()(const Surface_3& sphere, const Ray_3& r) const { 238 typename GT::Construct_point_on_3 point_on = 239 GT().construct_point_on_3_object(); 240 241 const Point& a = point_on(r, 0); 242 const Point& b = point_on(r, 1); 243 244 return private_intersection(sphere, a, b, Lambda_positive()); 245 } // end operator()(Surface_3, Ray_3) 246 operator()247 Object operator()(const Surface_3& sphere, const Line_3& l) const { 248 typename GT::Construct_point_on_3 point_on = 249 GT().construct_point_on_3_object(); 250 251 const Point& a = point_on(l, 0); 252 const Point& b = point_on(l, 1); 253 254 return private_intersection(sphere, a, b, Always_true()); 255 } // end operator()(Surface_3, Line_3) 256 257 /** Modifies s = [a, b] by clipping it to sphere. 258 Return false iff s is outside sphere. */ clip_segment(const Surface_3 & sphere,Point_3 & a,Point_3 & b)259 bool clip_segment(const Surface_3& sphere, 260 Point_3& a, 261 Point_3& b) const 262 { 263 typedef typename GT::Vector_3 Vector; 264 265 typename GT::Has_on_bounded_side_3 on_bounded_side_of_sphere = 266 GT().has_on_bounded_side_3_object(); 267 typename GT::Construct_vector_3 vector = 268 GT().construct_vector_3_object(); 269 typename GT::Construct_scaled_vector_3 scaled_vector = 270 GT().construct_scaled_vector_3_object(); 271 typename GT::Construct_translated_point_3 translated_point = 272 GT().construct_translated_point_3_object(); 273 274 const bool a_in_sphere = on_bounded_side_of_sphere(sphere, a); 275 const bool b_in_sphere = on_bounded_side_of_sphere(sphere, b); 276 277 if( a_in_sphere && b_in_sphere ) 278 return true; 279 280 int number_of_roots; 281 FT root_1, root_2; 282 283 boost::tie(number_of_roots, root_1, root_2) = 284 intersection_line_sphere_lambda(sphere, a, b); 285 286 #ifdef CGAL_SURFACE_MESHER_DEBUG_IMPLICIT_ORACLE 287 std::cerr << "Clip segment. Roots=(" 288 << root_1 << ", " << root_2 << ")\n"; 289 #endif 290 if( number_of_roots < 2 ) 291 return false; 292 293 if( root_1 > FT(1) ) // root_x \in ]1,\infinity[ 294 return false; // no intersection 295 296 if( root_1 >= FT(0) ) // root_1 \in [0,1[ 297 { // move point a 298 const Point original_a = a; 299 const Vector ab = vector(a, b); 300 a = translated_point(original_a, scaled_vector(ab, root_1)); 301 if( root_2 <= FT(1) ) /// move b iif root_2 <=1 302 { 303 b = translated_point(original_a, scaled_vector(ab, root_2)); 304 } 305 return true; 306 } 307 else // root_1 in ]-\infinity, 0[ 308 { // do not move point a 309 if( root_2 < FT(0) ) // root_x in ]-\infinity, 0[ 310 return false; // no intersection 311 else 312 { 313 const Vector ab = vector(a, b); 314 if( root_2 <= FT(1) ) 315 b = translated_point(a, scaled_vector(ab, root_2)); 316 return true; 317 } 318 } 319 } 320 321 /** The return value s is r clipped to sphere. 322 Return false iff r does not intersect sphere. */ clip_ray(const Surface_3 & sphere,const Ray_3 & r,Point_3 & a,Point_3 & b)323 bool clip_ray(const Surface_3& sphere, 324 const Ray_3& r, 325 Point_3& a, 326 Point_3& b) const 327 { 328 typedef typename GT::Vector_3 Vector; 329 330 typename GT::Construct_point_on_3 point_on = 331 GT().construct_point_on_3_object(); 332 typename GT::Construct_vector_3 vector = 333 GT().construct_vector_3_object(); 334 typename GT::Construct_scaled_vector_3 scaled_vector = 335 GT().construct_scaled_vector_3_object(); 336 typename GT::Construct_translated_point_3 translated_point = 337 GT().construct_translated_point_3_object(); 338 339 a = point_on(r, 0); 340 b = point_on(r, 1); 341 342 int number_of_roots; 343 FT root_1, root_2; 344 345 boost::tie(number_of_roots, root_1, root_2) = 346 intersection_line_sphere_lambda(sphere, a, b); 347 348 if( number_of_roots == 2 && root_2 > FT(0) ) 349 { 350 const Vector ab = vector(a, b); 351 b = translated_point(a, scaled_vector(ab, root_2)); 352 if(root_1 > FT(0)) 353 a = translated_point(a, scaled_vector(ab, root_1)); 354 // if root_1 <= 0, a is in the ball 355 return true; 356 } 357 // else r does not intersect the sphere 358 return false; 359 } // end clip_ray 360 361 /** The return value s=(ab) is l clipped to sphere. 362 Return false iff l does not intersect sphere. */ clip_line(const Surface_3 & sphere,const Line_3 & l,Point & a,Point & b)363 bool clip_line(const Surface_3& sphere, const Line_3& l, 364 Point& a, 365 Point& b) const 366 { 367 typedef typename GT::Vector_3 Vector; 368 369 typename GT::Construct_point_on_3 point_on = 370 GT().construct_point_on_3_object(); 371 typename GT::Construct_vector_3 vector = 372 GT().construct_vector_3_object(); 373 typename GT::Construct_scaled_vector_3 scaled_vector = 374 GT().construct_scaled_vector_3_object(); 375 typename GT::Construct_translated_point_3 translated_point = 376 GT().construct_translated_point_3_object(); 377 378 a = point_on(l, 0); 379 b = point_on(l, 1); 380 381 int number_of_roots; 382 FT root_1, root_2; 383 384 boost::tie(number_of_roots, root_1, root_2) = 385 intersection_line_sphere_lambda(sphere, a, b); 386 387 if( number_of_roots == 2 ) 388 { 389 const Point original_a = a; 390 const Vector ab = vector(a, b); 391 a = translated_point(original_a, scaled_vector(ab, root_1)); 392 b = translated_point(original_a, scaled_vector(ab, root_2)); 393 return true; 394 } 395 // else l does not intersect the sphere 396 return false; 397 } // end clip_line 398 399 }; // end nested class Intersect_3 400 401 class Construct_initial_points 402 { 403 const Self& oracle; 404 public: Construct_initial_points(const Self & oracle)405 Construct_initial_points(const Self& oracle) : oracle(oracle) 406 { 407 } 408 409 // Random points 410 template <typename OutputIteratorPoints> operator()411 OutputIteratorPoints operator() (const Surface_3& sphere, 412 OutputIteratorPoints out, 413 int n = 20) const // WARNING: why 20? 414 { 415 const Point center = 416 GT().construct_center_3_object()(sphere); 417 const FT squared_radius = 418 GT().compute_squared_radius_3_object()(sphere); 419 const double radius_in_double = 420 CGAL::sqrt(CGAL::to_double(squared_radius)); 421 422 typename CGAL::Random_points_on_sphere_3<Point, 423 Point_creator> random_point_on_sphere(radius_in_double); 424 typename GT::Construct_vector_3 vector_3 = 425 GT().construct_vector_3_object(); 426 typename GT::Construct_translated_point_3 translate = 427 GT().construct_translated_point_3_object(); 428 429 while (n-->0) 430 { 431 Point p = translate(*random_point_on_sphere++, 432 vector_3(CGAL::ORIGIN, center)); 433 oracle.get_visitor().new_point(p); 434 *out++ = p; 435 } 436 return out; 437 } 438 }; // end nested class Construct_initial_points 439 construct_initial_points_object()440 Construct_initial_points construct_initial_points_object() const 441 { 442 return Construct_initial_points(*this); 443 } 444 intersect_3_object()445 Intersect_3 intersect_3_object() const 446 { 447 return Intersect_3(*this); 448 } 449 }; // end Sphere_oracle_3 450 451 } // namespace Surface_mesher 452 453 } // namespace CGAL 454 455 #include <CGAL/enable_warnings.h> 456 457 #endif // CGAL_SURFACE_MESHER_SPHERE_ORACLE_3_H 458