1 // Copyright (c) 2005 Tel-Aviv University (Israel). 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/Envelope_3/include/CGAL/Env_sphere_traits_3.h $ 7 // $Id: Env_sphere_traits_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 // Author(s) : Michal Meyerovitch <gorgymic@post.tau.ac.il> 11 // Baruch Zukerman <baruchzu@post.tau.ac.il> 12 // Ron Wein <wein@post.tau.ac.il> 13 // Efi Fogel <fogel@post.tau.ac.il> 14 15 #ifndef CGAL_ENV_SPHERE_TRAITS_3_H 16 #define CGAL_ENV_SPHERE_TRAITS_3_H 17 18 #include <CGAL/license/Envelope_3.h> 19 20 21 #include <CGAL/Object.h> 22 #include <CGAL/enum.h> 23 #include <CGAL/Bbox_3.h> 24 #include <CGAL/Sphere_3.h> 25 #include <CGAL/functions_on_signs.h> 26 #include <CGAL/Envelope_3/Envelope_base.h> 27 #include <CGAL/int.h> 28 29 namespace CGAL { 30 31 template <class ConicTraits_2> 32 class Env_sphere_traits_3 : public ConicTraits_2 33 { 34 public: 35 typedef ConicTraits_2 Traits_2; 36 typedef Env_sphere_traits_3<Traits_2> Self; 37 38 typedef typename Traits_2::Point_2 Point_2; 39 typedef typename Traits_2::Curve_2 Curve_2; 40 typedef typename Traits_2::X_monotone_curve_2 X_monotone_curve_2; 41 typedef typename Traits_2::Multiplicity Multiplicity; 42 43 typedef typename Traits_2::Rat_kernel Rat_kernel; 44 typedef typename Traits_2::Alg_kernel Alg_kernel; 45 typedef typename Traits_2::Nt_traits Nt_traits; 46 47 typedef typename Rat_kernel::FT Rational; 48 typedef typename Rat_kernel::Point_2 Rat_point_2; 49 typedef typename Rat_kernel::Segment_2 Rat_segment_2; 50 typedef typename Rat_kernel::Line_2 Rat_line_2; 51 typedef typename Rat_kernel::Circle_2 Rat_circle_2; 52 typedef typename Rat_kernel::Point_3 Rat_point_3; 53 54 typedef typename Alg_kernel::FT Algebraic; 55 typedef typename Alg_kernel::Point_2 Alg_point_2; 56 typedef typename Alg_kernel::Circle_2 Alg_circle_2; 57 58 typedef typename Rat_kernel::Sphere_3 Surface_3; 59 60 // here we refer to the lower part of the sphere only 61 typedef Surface_3 Xy_monotone_surface_3; 62 protected: 63 typedef std::pair<X_monotone_curve_2, 64 Multiplicity> Intersection_curve; 65 66 public: 67 class Make_xy_monotone_3 { 68 protected: 69 const Self & parent; 70 71 public: Make_xy_monotone_3(const Self * p)72 Make_xy_monotone_3(const Self * p) : parent(*p) 73 {} 74 75 // create xy-monotone surfaces from a general surface 76 // return a past-the-end iterator 77 template <class OutputIterator> operator()78 OutputIterator operator()(const Surface_3& s, 79 bool is_lower, 80 OutputIterator o) const 81 { 82 // our half sphere is of same type as our full sphere since we always 83 // need only the lower/upper part of each sphere 84 parent.m_is_lower = is_lower; 85 *o++ = s; 86 return o; 87 } 88 }; 89 90 /*! Get a Make_xy_monotone_3 functor object. */ 91 Make_xy_monotone_3 make_xy_monotone_3_object()92 make_xy_monotone_3_object() const 93 { 94 return Make_xy_monotone_3(this); 95 } 96 97 class Construct_projected_boundary_2 98 { 99 protected: 100 const Self & parent; 101 102 public: Construct_projected_boundary_2(const Self * p)103 Construct_projected_boundary_2(const Self * p) : parent(*p) 104 {} 105 106 // insert into the OutputIterator all the (2d) curves of the boundary of 107 // the vertical projection of the surface on the xy-plane 108 // the OutputIterator value type is X_monotone_curve_2 109 template <class OutputIterator> 110 OutputIterator operator()111 operator()(const Xy_monotone_surface_3& s, OutputIterator o) const 112 { 113 // the projected boundary in a circle, with a projected center, 114 // and same radius 115 Rat_point_2 proj_center = parent.project(s.center()); 116 Rat_circle_2 circ(proj_center, s.squared_radius()); 117 Curve_2 curve(circ); 118 Object objs[2]; 119 CGAL_assertion_code(Object *p = ) 120 parent.make_x_monotone_2_object()(curve, objs); 121 CGAL_assertion(p == objs + 2); 122 123 X_monotone_curve_2 cv1, cv2; 124 125 CGAL_assertion(assign(cv1, objs[0])); 126 CGAL_assertion(assign(cv2, objs[1])); 127 128 assign(cv1, objs[0]); 129 assign(cv2, objs[1]); 130 131 if(cv1.is_lower()) 132 { 133 CGAL_assertion(cv2.is_upper()); 134 *o++ = make_object(std::make_pair(cv1, ON_POSITIVE_SIDE)); 135 *o++ = make_object(std::make_pair(cv2, ON_NEGATIVE_SIDE)); 136 } 137 else 138 { 139 CGAL_assertion(cv2.is_lower()); 140 *o++ = make_object(std::make_pair(cv1, ON_NEGATIVE_SIDE)); 141 *o++ = make_object(std::make_pair(cv2, ON_POSITIVE_SIDE)); 142 } 143 144 return o; 145 } 146 }; 147 148 /*! Get a Construct_projected_boundary_2 functor object. */ 149 Construct_projected_boundary_2 construct_projected_boundary_2_object()150 construct_projected_boundary_2_object() const 151 { 152 return Construct_projected_boundary_2(this); 153 } 154 155 class Construct_projected_intersections_2 156 { 157 protected: 158 const Self & parent; 159 160 public: Construct_projected_intersections_2(const Self * p)161 Construct_projected_intersections_2(const Self * p) : parent(*p) 162 {} 163 164 // insert into OutputIterator all the (2d) projections on the xy plane of 165 // the intersection objects between the 2 surfaces 166 // the data type of OutputIterator is Object 167 template <class OutputIterator> 168 OutputIterator operator()169 operator()(const Xy_monotone_surface_3& s1, 170 const Xy_monotone_surface_3& s2, 171 OutputIterator o) const 172 { 173 Rat_point_3 p1 = s1.center(); 174 Rat_point_3 p2 = s2.center(); 175 const Rational a1 = p1.x(), b1 = p1.y(), c1 = p1.z(), 176 a2 = p2.x(), b2 = p2.y(), c2 = p2.z(); 177 const Rational sqr_r1 = s1.squared_radius(), 178 sqr_r2 = s2.squared_radius(); 179 180 // // the spheres intersect iff d(p1, p2) <= (r1+r2) 181 // // squaring this twice, we get the condition 182 // // sqr_d^2 + (1-2*sqr_d)(sqr_r1 + sqr_r2) - 2*sqr_r1*sqr_r2 <= 0 183 // // with only rational numbers involved. 184 // // todo: check if it helps 185 // Rat_kernel ratk; 186 // Rational sqr_d = ratk.compute_squared_distance_3_object()(p1, p2); 187 // Sign do_inter = CGAL_NTS sign(sqr_d*sqr_d + (1-2*sqr_d)*(sqr_r1+sqr_r2)-2*sqr_r1*sqr_r2); 188 // if (do_inter == POSITIVE) 189 // return o; 190 191 Nt_traits nt_traits; 192 193 // we check if the centers of the 2 spheres have same z coordinate - 194 // in this case the potential projected intersection is a segment 195 // (or point) 196 if (CGAL_NTS compare(c1, c2) == EQUAL) 197 { 198 if (CGAL_NTS compare(b1, b2) == EQUAL) 199 { 200 if (CGAL_NTS compare(a1, a2) == EQUAL) 201 // the same center, we have no intersection 202 // (we don't return overlappings as intersections) 203 { 204 return o; 205 } 206 207 // here we have c1 == c2, b1 == b2 208 209 // the intersection lies on the plane 210 // m 211 // (1) x = -------- 212 // 2(a1-a2) 213 // where m = (a1^2 - a2^2 + sqr_r2 - sqr_r1) 214 215 // which is orthogonal to the xy-plane 216 // we look at the intersection of this plane with the plane z = c1 217 // to get the projected segment of the spheres intersection 218 219 // we get the quadratic equation: 220 // A*y^2 + B*y + C = 0, where: 221 // A = 4(a1-a2)^2 222 // B = -8b1(a1-a2)^2 223 // C = m^2 - 4ma1(a1-a2) + 4(a1^2 + b1^2 - sqr_r1)(a1-a2)^2 224 // (we multiplied the equation by 4(a1-a2)^2 to get integer 225 // coefficients) 226 227 Rational a_diff = a1 - a2; 228 Rational sqr_a_diff = a_diff * a_diff; 229 Rational sqr_a1 = a1*a1, sqr_a2 = a2*a2; 230 Rational sqr_b1 = b1*b1; 231 Rational m = sqr_a1 - sqr_a2 + sqr_r2 - sqr_r1; 232 233 Rational A = 4*sqr_a_diff; 234 Rational B = -8*b1*sqr_a_diff; 235 Rational C = 4*sqr_a_diff*(sqr_a1+sqr_b1-sqr_r1) + m*m - 4*m*a1*a_diff; 236 237 238 Algebraic ys[2]; 239 Algebraic *ys_end; 240 std::ptrdiff_t n_ys; 241 242 ys_end = nt_traits.solve_quadratic_equation(A, B, C, ys); 243 n_ys = ys_end - ys; 244 245 if (n_ys == 0) 246 { 247 return o; // no intersection 248 } 249 250 // the x coordinate of the solution points 251 Algebraic xs = m / (2*a_diff); 252 253 if (n_ys == 1) 254 { 255 // intersection is a point 256 Point_2 inter_point(xs , ys[0]); 257 *o++ = make_object(inter_point); 258 return o; 259 } 260 261 CGAL_assertion(n_ys == 2); 262 // intersection is a segment, with non-rational endpoints 263 // so we construct a COLLINEAR conic (with equation as in (1)) 264 // with 2 endpoints 265 266 Alg_point_2 end1(xs, ys[0]); 267 Alg_point_2 end2(xs, ys[1]); 268 269 // equation (1) is: 270 // 2(a1-a2)x - m = 0 271 272 Curve_2 res(0,0,0, 2*a_diff, 0, -m, COLLINEAR, end1, end2); 273 274 parent.add_curve_to_output(res, o); 275 //*o++ = make_object(Intersection_curve(res, TRANSVERSAL)); 276 } 277 else 278 { 279 // here we have c1 == c2, b1 != b2. 280 // the intersection lies on the plane 281 // -2(a1-a2)x + m 282 // (1) y = ---------------- 283 // 2(b1-b2) 284 // where m = (a1^2 + b1^2 - a2^2 - b2^2 + sqr_r2 - sqr_r1) 285 286 // which is orthogonal to the xy-plane 287 // we look at the intersection of this plane with the plane z = c1 288 // to get the projected segment of the spheres intersection 289 290 // we get the quadratic equation: 291 // A*x^2 + B*x + C = 0 292 // where 293 // (a1-a2)^2 m(a1-a2) 2b1(a1-a2) 294 // A = 1 + --------- B = -2a1 - ---------- + ---------- 295 // (b1-b2)^2 (b1-b2)^2 (b1-b2) 296 // and m^2 b1*m 297 // C = a1^2 + b1^2 - sqr_r1 + ---------- - ------- 298 // 4(b1-b2)^2 (b1-b2) 299 300 // since we can solve only equations with integer coefficients we 301 // multiply everything by 4(b1 - b2)^2, and get: 302 // D*x^2 + E*x + F = 0 where 303 // D = 4(b1-b2)^2 + 4(a1-a2)^2 304 // E = -8a1(b1-b2)^2 - 4m(a1-a2) + 8b1(a1-a2)(b1-b2) 305 // F = 4(a1^2 + b1^2 - sqr_r1)(b1-b2)^2 + m^2 - 4mb1(b1-b2) 306 // TODO: in the new version coefficients can be rationals 307 Rational a_diff = a1 - a2; 308 Rational b_diff = b1 - b2; 309 Rational sqr_a_diff = a_diff * a_diff; 310 Rational sqr_b_diff = b_diff * b_diff; 311 Rational sqr_a1 = a1*a1, sqr_a2 = a2*a2; 312 Rational sqr_b1 = b1*b1, sqr_b2 = b2*b2; 313 Rational m = sqr_a1 + sqr_b1 - sqr_a2 - sqr_b2 + sqr_r2 - sqr_r1; 314 315 Rational D = 4*sqr_a_diff + 4*sqr_b_diff; 316 Rational E = -8*a1*sqr_b_diff - 4*m*a_diff + 8*b1*a_diff*b_diff; 317 Rational F = 4*sqr_b_diff*(sqr_a1+sqr_b1-sqr_r1) + m*m - 4*m*b1*b_diff; 318 319 Algebraic xs[2]; 320 Algebraic *xs_end; 321 std::ptrdiff_t n_xs; 322 323 xs_end = nt_traits.solve_quadratic_equation(D, E, F, xs); 324 n_xs = xs_end - xs; 325 326 if (n_xs == 0) 327 { 328 return o; // no intersection 329 } 330 if (n_xs == 1) 331 { 332 // intersection is a point 333 Point_2 inter_point(xs[0], (-2*a_diff*xs[0] + m)/(2*b_diff) ); 334 *o++ = make_object(inter_point); 335 return o; 336 } 337 338 CGAL_assertion(n_xs == 2); 339 // intersection is a segment, with non-rational endpoints 340 // so we construct a COLLINEAR conic (with equation as in (1)) 341 // with 2 endpoints 342 Algebraic ys[2]; 343 ys[0] = (-2*a_diff*xs[0] + m)/(2*b_diff); 344 ys[1] = (-2*a_diff*xs[1] + m)/(2*b_diff); 345 346 Alg_point_2 end1(xs[0], ys[0]); 347 Alg_point_2 end2(xs[1], ys[1]); 348 349 // equation (1) is: 350 // 2(a1-a2)x + 2(b1-b2)y - m = 0 351 Curve_2 res(0,0,0, 2*a_diff, 2*b_diff, -m, COLLINEAR, end1, end2); 352 parent.add_curve_to_output(res, o); 353 //*o++ = make_object(Intersection_curve(res, TRANSVERSAL)); 354 } 355 356 } 357 // now the potential intersection is (a part of) a circle, 358 // and the projection is (a part of) an ellipse 359 else 360 { 361 // here we have c1 != c2. 362 // the intersection lies on the plane: 363 // -2(a1-a2)x -2(b1-b2)y + m 364 // (*) z = -------------------------- 365 // 2(c1-c2) 366 // where m = a1^2 + b1^2 + c1^2 - a2^2 - b2^2 - c2^2 + sqr_r2 - sqr_r1 367 // 368 // (**) since we deal with only half sphere we are interested 369 // in the part below min(c1, c2) in the case of lower envelope 370 // and in the part above max(c1, c2) in the case of upper envelope 371 // 372 // substituting z in the sphere equation we get the ellipse equation: 373 // r*x^2 + s*y^2 + t*x*y + u*x + v*y + w = 0 where: 374 // 375 // (a1-a2)^2 (b1-b2)^2 2(a1-a2)(b1-b2) 376 // r = 1 + --------- s = 1 + --------- t = --------------- 377 // (c1-c2)^2 (c1-c2)^2 (c1-c2)^2 378 // 379 // 2c1(a1-a2) m(a1-a2) 380 // u = -2a1 + ---------- - ---------- 381 // (c1-c2) (c1-c2)^2 382 // 383 // 2c1(b1-b2) m(b1-b2) 384 // v = -2b1 + ---------- - ---------- // here we have c1 != c2 385 // (c1-c2) (c1-c2)^2 386 // 387 // m*c1 m^2 388 // w = a1^2 + b1^2 + c1^2 - sqr_r1 - ------- + ---------- 389 // (c1-c2) 4(c1-c2)^2 390 391 // since we can solve only equations with integer coefficients we 392 // multiply everything by 4(c1-c2)^2, and get: 393 // R*x^2 + S*y^2 + T*x*y + U*x + V*y + W = 0 where: 394 // R = 4(c1-c2)^2 + 4(a1-a2)^2 395 // S = 4(c1-c2)^2 + 4(b1-b2)^2 396 // T = 8(a1-a2)(b1-b2) 397 // U = -8a1(c1-c2)^2 + 8c1(a1-a2)(c1-c2) - 4m(a1-a2) 398 // V = -8b1(c1-c2)^2 + 8c1(b1-b2)(c1-c2) - 4m(b1-b2) 399 // W = 4(a1^2 + b1^2 + c1^2 - sqr_r1)(c1-c2)^2 - 4mc1(c1-c2) + m^2 400 // TODO: in the new version coefficients can be rationals 401 Rational a_diff = a1 - a2; 402 Rational b_diff = b1 - b2; 403 Rational c_diff = c1 - c2; 404 405 Rational sqr_a_diff = a_diff * a_diff; 406 Rational sqr_b_diff = b_diff * b_diff; 407 Rational sqr_c_diff = c_diff * c_diff; 408 409 Rational sqr_a1 = a1*a1, sqr_a2 = a2*a2; 410 Rational sqr_b1 = b1*b1, sqr_b2 = b2*b2; 411 Rational sqr_c1 = c1*c1, sqr_c2 = c2*c2; 412 413 Rational m = sqr_a1 + sqr_b1 + sqr_c1 - 414 sqr_a2 - sqr_b2 - sqr_c2 + sqr_r2 - sqr_r1; 415 416 Rational R = 4*sqr_c_diff + 4*sqr_a_diff; 417 Rational S = 4*sqr_c_diff + 4*sqr_b_diff; 418 Rational T = 8*a_diff*b_diff; 419 Rational U = -8*a1*sqr_c_diff + 8*c1*c_diff*a_diff - 4*m*a_diff; 420 Rational V = -8*b1*sqr_c_diff + 8*c1*c_diff*b_diff - 4*m*b_diff; 421 Rational W = 4*sqr_c_diff*(sqr_a1+sqr_b1+sqr_c1-sqr_r1) - 422 4*m*c1*c_diff + m*m; 423 424 // if the full spheres do not intersect, the equation we get has no 425 // real solution, so we should check it: 426 bool ellipse_is_point = false; 427 if (!parent.is_valid_conic_equation(R, S, T, U, V, W, 428 ellipse_is_point)) 429 { 430 return o; 431 } 432 433 // we need only a part of the ellipse (as stated in (**)) so we 434 // construct the cutting line, which is: 435 // equation (*) <= min(c1,c2) -- for lower envelope 436 // equation (*) >= max(c1,c2) -- for upper envelope 437 Rational z_plane; 438 if (parent.m_is_lower) 439 z_plane = ((c1 < c2) ? c1 : c2); 440 else 441 z_plane = ((c1 > c2) ? c1 : c2); 442 443 444 // we get (for lower envelope) 445 // -2(a1-a2)x -2(b1-b2)y + m 446 // (*) z = -------------------------- <= z_plane 447 // 2(c1-c2) 448 // and since we need integer coefficients, and also need to be in the 449 // positive side of the line (i.e., our halfplane equation should be of 450 // type ax+by+c >= 0), we have: 451 // sign_c_diff* [2(a1-a2)x + 2(b1-b2)y - m + 2(c1-c2)*z_plane] >= 0 452 453 // for upper envelope, we should multiply the line equation by -1 454 int envelope_coef = 1; 455 if (!parent.m_is_lower) 456 envelope_coef = -1; 457 458 Sign sign_c_diff = CGAL_NTS sign(c_diff); 459 Rational la = envelope_coef*2*a_diff*sign_c_diff; 460 Rational lb = envelope_coef*2*b_diff*sign_c_diff; 461 Rational lc = envelope_coef*sign_c_diff*(2*c_diff*z_plane - m); 462 463 if (ellipse_is_point) 464 { 465 // as specified in the is_valid_conic_equation method, the 466 // intersection point is: 467 // 468 Rational px = S*(4*U - T*V)/(T*T - 4*S*R); 469 px = px / 2; 470 Rational py = -(T*px + V)/(2*S); 471 // should check if the point is in the non-negative side of the 472 // line 473 if (CGAL_NTS sign(la*px + lb*py +lc) != NEGATIVE) 474 { 475 *o++ = make_object(Point_2(px, py)); 476 } 477 return o; 478 } 479 480 // EBEB 2012/06/29: Added because of 481 // no matching function for call to 'compare(CGAL::Env_sphere_traits_3<Conic_traits_2>::Rational&, int) 482 Rational zero(0); 483 484 // if (a1==a2) and (b1==b2) (*) is a plane parallel to the xy-plane 485 // and either all ellipse (which should be a circle) is the 486 // intersection - in which case lc >= 0 487 // or there is no intersection at all between the 2 half spheres - 488 // in which case lc < 0 489 if (CGAL_NTS compare(a_diff, zero) == EQUAL && 490 CGAL_NTS compare(b_diff, zero) == EQUAL) 491 { 492 Sign sign_lc = CGAL_NTS sign(lc); 493 if (sign_lc != NEGATIVE) 494 { 495 Curve_2 res(R, S, T, U, V, W); 496 parent.add_curve_to_output(res, o); 497 //*o++ = make_object(Intersection_curve(res, TRANSVERSAL)); 498 } 499 return o; 500 } 501 502 // find the intersection of the line 503 // la * x + lb * y + lc = 0 504 // with the conic 505 // R*x^2 + S*y^2 + T*xy + U*x + V*y + W = 0 506 Alg_point_2 source, target, pmid; 507 std::ptrdiff_t n_inter_points; 508 if (CGAL_NTS compare(lb, zero) != EQUAL) 509 { 510 // Find the x-coordinates of the intersection points of the conic 511 // curve and the line y = -(la*x + lc) / lb: 512 // we get a quadratic equation Ax^2 + Bx + C = 0 513 // where A = lb*lb*R + la*(la*S - lb*T) 514 // B = 2*la*lc*S - lb*(lc*T + la*V - lb*U) 515 // C = S*lc*lc + lb*(lb*W - lc*V) 516 Rational A = lb*lb*R + la*(la*S - lb*T), 517 B = 2*la*lc*S - lb*(lc*T + la*V - lb*U), 518 C = S*lc*lc + lb*(lb*W - lc*V); 519 520 Algebraic inter_xs[2]; 521 Algebraic *inter_xs_end; 522 523 inter_xs_end = nt_traits.solve_quadratic_equation(A, B, C, inter_xs); 524 n_inter_points = inter_xs_end - inter_xs; 525 526 if (n_inter_points > 0) 527 source = Alg_point_2(inter_xs[0], 528 -(la*inter_xs[0] + lc) / Algebraic(lb)); 529 if (n_inter_points == 2) 530 { 531 target = Alg_point_2(inter_xs[1], 532 -(la*inter_xs[1] + lc) / Algebraic(lb)); 533 534 // Get the conic points whose x-coordinate are in the middle of the 535 // two endpoints. 536 // since inter_xs[0] and inter_xs[1] are the roots of a quadratic 537 // equation Ax^2 + Bx + C = 0, their sum is -B/A which is rational 538 Algebraic x_mid = Algebraic(Rational(-B/(2*A))); 539 //Algebraic x_mid = (inter_xs[0] + inter_xs[1]) / 2; 540 Alg_point_2 x_mid_point(x_mid, 0); 541 542 CGAL_precondition_code(int x_mid_n_y_points;); 543 Alg_point_2 x_mid_y_points[2]; 544 545 Curve_2 inter_cv(R, S, T, U, V, W); 546 547 CGAL_precondition_code(x_mid_n_y_points = ) 548 inter_cv.points_at_x(x_mid_point, x_mid_y_points); 549 550 CGAL_precondition(x_mid_n_y_points > 0); 551 552 Algebraic y1 = x_mid_y_points[0].y(), y2 = x_mid_y_points[1].y(); 553 if (CGAL_NTS compare(Algebraic(la)*x_mid + Algebraic(lb)*y1 + Algebraic(lc), 554 Algebraic(0) 555 ) == LARGER) 556 { 557 pmid = Alg_point_2(x_mid, y1); 558 } 559 else 560 { 561 CGAL_assertion(CGAL_NTS compare(Algebraic(la)*x_mid + Algebraic(lb)*y2 + Algebraic(lc), 562 Algebraic(0)) == LARGER); 563 pmid = Alg_point_2(x_mid, y2); 564 } 565 } 566 } 567 else 568 { // lb == 0 569 CGAL_assertion(CGAL_NTS compare(la, zero) != EQUAL); 570 571 // Find the intersection of the vertical line x = -lc / la: 572 Rational inter_x = -lc/la; 573 // we should solve the quadratic equation A*y^2 + B*y + C = 0 574 // where A = S 575 // B = T*inter_x + V 576 // C = R*inter_x^2 + U*inter_x + W 577 Rational A = S, 578 B = T*inter_x + V, 579 C = R*inter_x*inter_x + U*inter_x + W; 580 581 Algebraic inter_points[2]; 582 Algebraic *inter_points_end; 583 584 inter_points_end = 585 nt_traits.solve_quadratic_equation(A, B, C, inter_points); 586 n_inter_points = inter_points_end - inter_points; 587 588 if (n_inter_points > 0) 589 source = Alg_point_2(Algebraic(inter_x), inter_points[0]); 590 if (n_inter_points == 2) 591 { 592 target = Alg_point_2(Algebraic(inter_x), inter_points[1]); 593 594 // Get the conic points whose y-coordinate are in the middle of the 595 // two endpoints. 596 // since inter_points[0] & inter_points[1] are roots of quadratic 597 // equation, their sum is -B/A, and mid_y is -B/2A 598 Algebraic y_mid = Algebraic(Rational(-B/(2*A))); 599 600 Alg_point_2 y_mid_point(0, y_mid); 601 Alg_point_2 y_mid_x_points[2]; 602 Curve_2 inter_cv(R, S, T, U, V, W); 603 604 CGAL_precondition_code(int y_mid_n_x_points =) 605 inter_cv.points_at_y(y_mid_point, y_mid_x_points); 606 607 CGAL_precondition(y_mid_n_x_points > 0); 608 609 Algebraic x1 = y_mid_x_points[0].x(), x2 = y_mid_x_points[1].x(); 610 if (CGAL_NTS compare( 611 Algebraic(la)*x1 + Algebraic(lb)*y_mid + Algebraic(lc), 612 Algebraic(0)) == LARGER) 613 { 614 pmid = Alg_point_2(x1, y_mid); 615 } 616 else 617 { 618 CGAL_assertion(CGAL_NTS compare ( 619 Algebraic(la)*x2 + Algebraic(lb)*y_mid + Algebraic(lc), 620 Algebraic(0)) == LARGER); 621 pmid = Alg_point_2(x2, Algebraic(y_mid)); 622 } 623 } 624 } 625 626 if (n_inter_points < 2) 627 { 628 // we should check whether the ellipse is in the positive side of the 629 // line - in which case we return the full ellipse 630 // or not - in which case there is no intersection if 631 // n_inter_points = 0, and a point intersection (equal to source) 632 // if n_inter_points = 1 633 634 // for this, we find a point inside the ellipse and substitute 635 // its coordinates in the line equation 636 637 Curve_2 inter_cv(R, S, T, U, V, W); 638 Alg_point_2 vtan_ps[2]; 639 640 CGAL_assertion_code(int n_vtan_ps =) 641 inter_cv.vertical_tangency_points(vtan_ps); 642 643 644 CGAL_assertion(n_vtan_ps == 2); 645 646 Algebraic lval = Algebraic(la)*vtan_ps[0].x() + 647 Algebraic(lb)*vtan_ps[0].y() + Algebraic(lc); 648 Sign lval_sign = CGAL_NTS sign(lval); 649 if (lval_sign == POSITIVE) 650 { 651 // the full ellipse is in the positive side 652 parent.add_curve_to_output(inter_cv, o); 653 //*o++ = make_object(Intersection_curve(inter_cv, TRANSVERSAL)); 654 return o; 655 } 656 else if (lval_sign == NEGATIVE) 657 { 658 // the full ellipse is in the negative side, except maybe the point 659 // source in the case n_inter_points = 1 (which lies on the line) 660 if (n_inter_points == 1) 661 *o++ = make_object(Point_2(source)); 662 return o; 663 } 664 665 CGAL_assertion(lval_sign == ZERO); 666 // in this case lval_sign lies on the line, so it must be that 667 // n_inter_points == 1 and source = vtan_ps[0] 668 CGAL_assertion(n_inter_points == 1 && source == vtan_ps[0]); 669 // so we try the other vertical tangency point 670 lval = Algebraic(la)*vtan_ps[1].x() + 671 Algebraic(lb)*vtan_ps[1].y() + Algebraic(lc); 672 lval_sign = CGAL_NTS sign(lval); 673 CGAL_assertion(lval_sign != ZERO); 674 675 if (lval_sign == POSITIVE) 676 parent.add_curve_to_output(inter_cv, o); 677 //*o++ = make_object(Intersection_curve(inter_cv, TRANSVERSAL)); 678 else 679 *o++ = make_object(Point_2(source)); 680 681 return o; 682 } 683 684 CGAL_assertion(n_inter_points == 2); 685 686 // find the correct orientation of the conic between the 2 endpoints 687 // it should lie on the positive side of line 688 // If the mid-point forms a left-turn with the source and the target 689 // points, the orientation is positive (going counterclockwise). 690 // Otherwise, it is negative (going clockwise). 691 Alg_kernel k; 692 typename Alg_kernel::Orientation_2 orient_f = k.orientation_2_object(); 693 Orientation orient; 694 if (orient_f(source, pmid, target) == LEFT_TURN) 695 orient = CGAL::COUNTERCLOCKWISE; 696 else 697 orient = CGAL::CLOCKWISE; 698 699 Curve_2 res(R, S, T, U, V, W, orient, source, target); 700 CGAL_assertion(res.is_valid()); 701 parent.add_curve_to_output(res, o); 702 //*o++ = make_object(Intersection_curve(res, TRANSVERSAL)); 703 } 704 705 return o; 706 } 707 }; 708 709 /*! Get a Construct_projected_intersections_2 functor object. */ 710 Construct_projected_intersections_2 construct_projected_intersections_2_object()711 construct_projected_intersections_2_object() const 712 { 713 return Construct_projected_intersections_2(this); 714 } 715 716 class Compare_z_at_xy_3 717 { 718 protected: 719 const Self & parent; 720 721 public: Compare_z_at_xy_3(const Self * p)722 Compare_z_at_xy_3(const Self * p) : parent(*p) 723 {} 724 725 // check which of the surfaces is closer to the envelope at the xy 726 // coordinates of point (i.e. lower if computing the lower envelope, or 727 // upper if computing the upper envelope) 728 // precondition: the surfaces are defined in point operator()729 Comparison_result operator()(const Point_2& p, 730 const Xy_monotone_surface_3& s1, 731 const Xy_monotone_surface_3& s2) const 732 { 733 Comparison_result c2 = compare_in_point_second_method(p, s1, s2); 734 735 CGAL_expensive_assertion_code( 736 Comparison_result c1 = compare_in_point_first_method(p, s1, s2); 737 ); 738 CGAL_expensive_assertion(c1 == c2); 739 740 return c2; 741 } 742 743 // check which of the surfaces is closer to the envelope at the xy 744 // coordinates of cv (i.e. lower if computing the lower envelope, or upper 745 // if computing the upper envelope) 746 // precondition: the surfaces are defined in all points of cv, and the 747 // answer is the same for each of these points operator()748 Comparison_result operator()(const X_monotone_curve_2& cv, 749 const Xy_monotone_surface_3& s1, 750 const Xy_monotone_surface_3& s2) const 751 { 752 // we compute a middle point on cv and use the previous function 753 Point_2 mid = parent.construct_middle_point(cv); 754 Comparison_result res = 755 parent.compare_z_at_xy_3_object()(mid, s1, s2); 756 return res; 757 } 758 759 protected: 760 // first method of compare in point, calculates the z value of both 761 // surfaces, and compares them 762 Comparison_result compare_in_point_first_method(const Point_2 & p,const Xy_monotone_surface_3 & s1,const Xy_monotone_surface_3 & s2)763 compare_in_point_first_method(const Point_2& p, 764 const Xy_monotone_surface_3& s1, 765 const Xy_monotone_surface_3& s2) const 766 { 767 // find the z coordinates of surface 1 over p 768 Algebraic z1 = parent.compute_envelope_z_in_point(p, s1); 769 // find the z coordinates of surface 2 over p 770 Algebraic z2 = parent.compute_envelope_z_in_point(p, s2); 771 772 Sign res = CGAL_NTS sign(z1 - z2); 773 if (parent.m_is_lower) 774 return res; 775 else 776 return -res; 777 } 778 779 // second method of compare in point 780 // p = (x1, y1) 781 // s1: (x-a1)^2 + (y-b1)^2 + (z-c1)^2 = r1^2 782 // s2: (x-a2)^2 + (y-b2)^2 + (z-c2)^2 = r2^2 783 // (both lower parts or upper parts) 784 // then in point p we get: 785 // s1(p): (z-c1)^2 = r1^2 - (x1-a1)^2 - (y1-b1)^2 = A1 786 // s2(p): (z-c2)^2 = r2^2 - (x1-a2)^2 - (y1-b2)^2 = A2 787 // so we get z - ci = +- sqrt(Ai) where -sqrt(Ai) is for the lower part 788 // and +sqrt(Ai) is for the upper part 789 // we now need to compute the sign of: 790 // c1 - sqrt(A1) - (c2 - sqrt(A2)) - for lower envelope 791 // c1 + sqrt(A1) - (c2 + sqrt(A2)) - for upper envelope 792 Comparison_result compare_in_point_second_method(const Point_2 & p,const Xy_monotone_surface_3 & s1,const Xy_monotone_surface_3 & s2)793 compare_in_point_second_method(const Point_2& p, 794 const Xy_monotone_surface_3& s1, 795 const Xy_monotone_surface_3& s2) const 796 { 797 Rat_point_3 p1 = s1.center(); 798 Rat_point_3 p2 = s2.center(); 799 const Rational a1 = p1.x(), b1 = p1.y(), c1 = p1.z(), 800 a2 = p2.x(), b2 = p2.y(), c2 = p2.z(); 801 const Rational sqr_r1 = s1.squared_radius(), 802 sqr_r2 = s2.squared_radius(); 803 const Algebraic x1 = p.x(), y1 = p.y(); 804 805 Rational c_diff = c1 - c2; 806 Algebraic x_diff1 = x1 - a1, y_diff1 = y1 - b1; 807 Algebraic x_diff2 = x1 - a2, y_diff2 = y1 - b2; 808 809 Algebraic A1 = sqr_r1 - x_diff1*x_diff1 - y_diff1*y_diff1; 810 Algebraic A2 = sqr_r2 - x_diff2*x_diff2 - y_diff2*y_diff2; 811 812 if (CGAL_NTS sign(A1) == NEGATIVE) 813 std::cout << "A1 = " << A1 << std::endl; 814 if (CGAL_NTS sign(A2) == NEGATIVE) 815 std::cout << "A2 = " << A2 << std::endl; 816 817 Sign res; 818 // sign_a_plus_b_x_sqrt_e_plus_c_x_sqrt_f is a CGAL method which 819 // computes the sign of quantity: a + b * sqrt(e) + c * sqrt(f) 820 821 res = CGAL::sign_a_plus_b_x_sqrt_e_plus_c_x_sqrt_f(Algebraic(c_diff), 822 Algebraic(-1), 823 Algebraic(1), 824 A1, 825 A2); 826 return res; 827 } 828 }; 829 830 /*! Get a Compare_z_at_xy_3 functor object. */ 831 Compare_z_at_xy_3 compare_z_at_xy_3_object()832 compare_z_at_xy_3_object() const 833 { 834 return Compare_z_at_xy_3(this); 835 } 836 837 class Compare_z_at_xy_above_3 838 { 839 protected: 840 const Self & parent; 841 842 public: Compare_z_at_xy_above_3(const Self * p)843 Compare_z_at_xy_above_3(const Self * p) : parent(*p) 844 {} 845 846 // check which of the surfaces is closer to the envelope on the points above 847 // the curve cv (i.e. lower if computing the lower envelope, or upper if 848 // computing the upper envelope) 849 // precondition: the surfaces are defined above cv 850 // the choise between s1 and s2 for the envelope is the same 851 // for every point in the infinitesimal region above cv 852 // the surfaces are EQUAL over the curve cv 853 Comparison_result operator()854 operator()(const X_monotone_curve_2& cv, 855 const Xy_monotone_surface_3& s1, 856 const Xy_monotone_surface_3& s2) const 857 { 858 Comparison_result res = parent.compare_on_side(cv, s1, s2, false); 859 return res; 860 } 861 }; 862 863 /*! Get a Compare_z_at_xy_above_3 functor object. */ 864 Compare_z_at_xy_above_3 compare_z_at_xy_above_3_object()865 compare_z_at_xy_above_3_object() const 866 { 867 return Compare_z_at_xy_above_3(this); 868 } 869 870 class Compare_z_at_xy_below_3 871 { 872 protected: 873 const Self & parent; 874 875 public: Compare_z_at_xy_below_3(const Self * p)876 Compare_z_at_xy_below_3(const Self * p) : parent(*p) 877 {} 878 879 Comparison_result operator()880 operator()(const X_monotone_curve_2& cv, 881 const Xy_monotone_surface_3& s1, 882 const Xy_monotone_surface_3& s2) const 883 { 884 Comparison_result res = parent.compare_on_side(cv, s1, s2, true); 885 return res; 886 } 887 }; 888 889 /*! Get a Compare_z_at_xy_below_3 functor object. */ 890 Compare_z_at_xy_below_3 compare_z_at_xy_below_3_object()891 compare_z_at_xy_below_3_object() const 892 { 893 return Compare_z_at_xy_below_3(this); 894 } 895 896 /***************************************************************************/ 897 898 // public method needed for testing 899 900 // checks if point is in the xy-range of surf 901 class Is_defined_over 902 { 903 protected: 904 const Self & parent; 905 906 public: Is_defined_over(const Self * p)907 Is_defined_over(const Self * p) : parent(*p) 908 {} 909 910 // checks if point is in the xy-range of surf operator()911 bool operator()(const Point_2& p, const Xy_monotone_surface_3& s) const 912 { 913 // project the surface on the plane 914 Rat_point_2 proj_center = parent.project(s.center()); 915 Rat_circle_2 boundary(proj_center, s.squared_radius()); 916 917 Nt_traits nt_traits; 918 Alg_kernel k; 919 Alg_point_2 aproj_center(proj_center.x(), proj_center.y()); 920 Alg_circle_2 aboundary(aproj_center, nt_traits.convert(s.squared_radius())); 921 922 // check if the projected point is inside the projected boundary 923 return (!k.has_on_unbounded_side_2_object()(aboundary, p)); 924 } 925 }; 926 927 /*! Get a Is_defined_over functor object. */ is_defined_over_object()928 Is_defined_over is_defined_over_object() const 929 { 930 return Is_defined_over(this); 931 } 932 933 /***************************************************************************/ 934 935 // helper methods 936 937 // compare the surfaces over the side (as specified in the compare_on_right 938 // parameter) of the curve, assuming they are defined there compare_on_side(const X_monotone_curve_2 & cv,const Xy_monotone_surface_3 & s1,const Xy_monotone_surface_3 & s2,bool compare_on_right)939 Comparison_result compare_on_side(const X_monotone_curve_2& cv, 940 const Xy_monotone_surface_3& s1, 941 const Xy_monotone_surface_3& s2, 942 bool compare_on_right) const 943 { 944 // cv(x,y) : r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 945 // let p be the leftmost endpoint of cv, p=(x0, y0) 946 // the tangence of cv at p is a line. on the infinitesimal region 947 // near p the relation between the surfaces to the right of cv is 948 // the same as the relation between the tangences of the surfaces 949 // in p to the right of this line (unless the tangence of the surface 950 // is vertical) 951 952 // we take a point in the internal of the curve, not an endpoint 953 // since we assume that such point represents better what is going 954 // on all internal curve points 955 Point_2 cv_point = construct_middle_point(cv); 956 Algebraic x0 = cv_point.x(), y0 = cv_point.y(); 957 958 // d(cv)/dx : 2r*x + 2s*y*dy/dx + t*y + t*x*dy/dx +u + v*dy/dx = 0 959 // in point p=(x0,y0) we get 960 // dy m 961 // -- = y' = - where m = -2rx0 -ty0 - u 962 // dx n n = 2sy0 + tx0 + v 963 // if n != 0 (if n = 0 we have a vertical line). 964 // 965 // So the tangence equation (in p) is: 966 // n = 0: x = x0 967 // n != 0: y - y0 = y'(x-x0) ==> -y'x + y + (y'x0 - y0) = 0 968 // and in general we have: 969 // -m*x + n*y + (m*x0 -n*y0) = 0 (with integer coordinates) 970 const Rational r = cv.r(), s = cv.s(), t = cv.t(), 971 u = cv.u(), v = cv.v(), w = cv.w(); 972 Algebraic m = -1 * (2*r*x0 + t*y0 + u); 973 Algebraic n = 2*s*y0 + t*x0 + v; 974 // line coefficients: A3, B3, C3 975 Algebraic A3 = -1*m, B3 = n, C3 = m*x0 - n*y0; 976 977 // the tangences of the spheres (in point (x0,y0,z0)): 978 Algebraic z0 = compute_envelope_z_in_point(cv_point, s1); 979 980 // we assume the surfaces are equal over cv: 981 CGAL_expensive_precondition_code( 982 Algebraic z0_2 = compute_envelope_z_in_point(cv_point, s2); 983 ) 984 // this test can be very time consuming ... 985 CGAL_expensive_precondition(CGAL_NTS compare(z0, z0_2) == EQUAL); 986 987 // the sphere i : fi(x,y,z) = (x-ai)^2 + (y-bi)^2 + (z-ci)^2 - ri^2 = 0 988 // dfi / dx = 2(x-ai) + 2(z-ci)*dz/dx = 0 989 // dfi / dy = 2(y-bi) + 2(z-ci)*dz/dy = 0 990 // if z = ci the tangent plane is vertical - if only one of the tangent 991 // planes is vertical, them its sphere wins (i.e. is on envelope). 992 // we assume not both are tangent, since this means that they are the 993 // same sphere 994 // if z != ci the tangent plane is: 995 996 // z-z0 = dz/dx (x-x0) + dz/dy (y-y0) 997 // ==> 998 // (x0-ai)(x-x0) + (y0-bi)(y-y0) + (z0-ci)(z-z0) = 0 999 // Ai*x + Bi*y + Ci*z + Di = 0 1000 // where Ai = (x0-ai) 1001 // Bi = (y0-bi) 1002 // Ci = (z0-ci) 1003 // Di = -(x0-ai)x0 - (y0-bi)y0 - (z0-ci)z0 1004 // 1005 // and we solve the problem as for triangles 1006 Rat_point_3 p1 = s1.center(); 1007 Rat_point_3 p2 = s2.center(); 1008 const Rational a1 = p1.x(), b1 = p1.y(), c1 = p1.z(), 1009 a2 = p2.x(), b2 = p2.y(), c2 = p2.z(); 1010 Algebraic A1 = x0 - a1, B1 = y0 - b1, C1 = z0 - c1; 1011 Algebraic A2 = x0 - a2, B2 = y0 - b2, C2 = z0 - c2; 1012 if (C1 != 0 && C2 != 0) 1013 { 1014 Sign sign1 = CGAL_NTS sign((A2*A3+B2*B3)/C2-(A1*A3+B1*B3)/C1); 1015 // to make sure the direction is correct, we take a second point on the 1016 // line: for vertical line we take (x0, y0+1) 1017 // otherwise we take (x0+1, y0+ m/n) 1018 Algebraic x1, y1; 1019 if (n == 0) 1020 { 1021 x1 = x0; 1022 y1 = y0+1; 1023 } 1024 else 1025 { 1026 x1 = x0+1; 1027 y1 = y0 + (m/n); 1028 } 1029 Sign sign2 = CGAL_NTS sign(-B3*x1+A3*y1-(-B3*x0+A3*y0)); 1030 1031 // the answer negates according to the side of the line we ask of 1032 Sign sign3 = (compare_on_right ? (CGAL_NTS sign(1)) : 1033 (CGAL_NTS sign(-1))); 1034 1035 return sign1 * sign2 * sign3; 1036 1037 } 1038 else if (C1 != 0 && C2 == 0) 1039 { 1040 // sphere 2 is on the envelope (both lower & upper) 1041 return LARGER; 1042 } 1043 else if (C1 == 0 && C2 != 0) 1044 1045 { 1046 // sphere 1 is on the envelope (both lower & upper) 1047 return SMALLER; 1048 } 1049 else 1050 CGAL_error(); 1051 1052 return EQUAL; 1053 } 1054 project(const Rat_point_3 & p)1055 Rat_point_2 project(const Rat_point_3& p) const 1056 { 1057 return Rat_point_2(p.x(), p.y()); 1058 } 1059 1060 // compute the z coordinate of the surface s in point p on the envelope 1061 // (i.e. take lower point if lower envelope, upper otherwise) 1062 // precondition: s is defined at p compute_envelope_z_in_point(const Point_2 & p,const Xy_monotone_surface_3 & s)1063 Algebraic compute_envelope_z_in_point(const Point_2& p, 1064 const Xy_monotone_surface_3& s) const 1065 { 1066 Algebraic res; 1067 1068 // the point coordinates 1069 const Algebraic x1 = p.x(), y1 = p.y(); 1070 1071 // the surface equations 1072 Rat_point_3 center = s.center(); 1073 const Rational a = center.x(), b = center.y(), c = center.z(); 1074 const Rational sqr_r = s.squared_radius(); 1075 1076 // we substitute x1 and y1 in the equation of s 1077 // (x-a)^2 + (y-b)^2 + (z-c)^2 = r^2 1078 // and get a quadratic equation of z: 1079 // z^2 - 2cz + [(x1-a)^2 + (y1-b)^2 + c^2 - r^2] = 0 1080 Algebraic x_diff = x1 - a, y_diff = y1 - b; 1081 // the coefficients are: 1082 Algebraic A = 1, 1083 B = -2*c, 1084 C = x_diff*x_diff + y_diff*y_diff + c*c - sqr_r; 1085 1086 Algebraic zs[2]; 1087 Algebraic *zs_end; 1088 std::ptrdiff_t n_zs; 1089 1090 Nt_traits nt_traits; 1091 zs_end = nt_traits.solve_quadratic_equation(A, B, C, zs); 1092 n_zs = zs_end - zs; 1093 1094 CGAL_precondition(n_zs > 0); 1095 1096 if (n_zs == 1) 1097 // only one point is defined at p, this is the result 1098 return zs[0]; 1099 CGAL_assertion(n_zs == 2); 1100 1101 Comparison_result comp = CGAL_NTS compare(zs[0], zs[1]); 1102 if (m_is_lower) 1103 res = ((comp == SMALLER) ? zs[0] : zs[1]); 1104 else 1105 res = ((comp == LARGER) ? zs[0] : zs[1]); 1106 1107 return res; 1108 } 1109 1110 // construct the point in the middle of cv construct_middle_point(const X_monotone_curve_2 & cv)1111 Point_2 construct_middle_point(const X_monotone_curve_2& cv) const 1112 { 1113 // get the x-value of the middle point 1114 Alg_kernel k; 1115 Alg_point_2 mid_x = k.construct_midpoint_2_object()(cv.source(), 1116 cv.target()); 1117 1118 // TODO_NEW_DESIGN - this is not implemented in X_monotone_curve_2, but maybe we want it there? 1119 // if (cv.is_segment()) 1120 // return mid_x; 1121 if (cv.is_vertical()) 1122 return Point_2(mid_x); 1123 1124 return Point_2(cv.point_at_x(mid_x)); 1125 } 1126 1127 1128 // for the test construct_middle_point(const Point_2 & p1,const Point_2 & p2)1129 Point_2 construct_middle_point(const Point_2& p1, const Point_2& p2) const 1130 { 1131 Alg_kernel k; 1132 return Point_2(k.construct_midpoint_2_object()(p1, p2)); 1133 } 1134 1135 // check if the equation 1136 // r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 1137 // has real solutions 1138 // is_point is set to true if the equation represents just one point 1139 template <class NT> is_valid_conic_equation(const NT & r,const NT & s,const NT & t,const NT & u,const NT & v,const NT & w,bool & is_point)1140 bool is_valid_conic_equation(const NT& r, const NT& s, const NT& t, 1141 const NT& u, const NT& v, const NT& w, 1142 bool &is_point) const 1143 { 1144 // initialize is_point to false, and will change it when we detect 1145 // that the equation represents a point 1146 is_point = false; 1147 // (*) r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 1148 // we fix x, and get a 1-variable quadratic equation: 1149 // (**) s*y^2 + (tx + v)*y + (rx^2 + ux + w) = 0 1150 // (*) has real solution (x,y) iff there exists x such that (**) has a 1151 // solution y, i.e. discriminant(**) >= 0 1152 // discriminant(**) = f(x) = (tx + v)^2 -4s(rx^2 + ux + w) 1153 // = (t^2 - 4sr)*x^2 + (stv - 4su)*x + (v^2 - 4sw) 1154 // = A*x^2 + B*x + C >= 0 1155 // where A = t^2 - 4sr 1156 // B = stv - 4su 1157 // C = v^2 - 4sw 1158 // so we should check if there exists x such that f(x) >= 0 1159 // if A > 0 we have a smiling parabula, and a positive answer 1160 // (the conic equation in this case represents hyperbola or 2 1161 // intersecting lines) 1162 Sign sign_A = CGAL_NTS sign(t*t - 4*s*r); 1163 if (sign_A == POSITIVE) 1164 return true; 1165 // if A < 0 we have a sad parabula, so we should check if it crosses the 1166 // x-axis, i.e. if the equation f(x) = 0 has a real solution x. 1167 // this means that discriminant(f(x)) >= 0 1168 // discriminant(f(x)) = B^2 - 4AC 1169 // = (2tv-4su)^2 - 4(t^2-4sr)(v^2-4sw) 1170 // = s(-tvu + su^2 + wt^2 + rv^2 - 4srw) 1171 if (sign_A == NEGATIVE) 1172 { 1173 // (in this case the conic equation represents ellipse, circle, point 1174 // or no curve) 1175 Sign sign_s = CGAL_NTS sign(s); 1176 Sign sign_eq = CGAL_NTS sign(-t*v*u + s*u*u + w*t*t + r*v*v - 4*s*r*w); 1177 // if sign_eq = 0 then discriminant(f(x))=0, and so we have only one x 1178 // solution for f(x), say x0. since we get f(x0)=0 and f(x)<0 forall 1179 // x!=x0, we have only one solution for (**). So the equation represents 1180 // a point with coordinates x0=-B/2A, y0=-(tx0 + v)/2s 1181 if (sign_eq == ZERO) 1182 is_point = true; 1183 1184 Sign sign_disc = CGAL_NTS sign(int(sign_s * sign_eq)); 1185 return (sign_disc != NEGATIVE); 1186 } 1187 // if A = 0 we get (***) f(x) = (stv - 4su)*x + (v^2 - 4sw) = B*x + C 1188 // if B != 0 we get a line equation, which always has x 1189 // such that f(x) >= 0 1190 // if B = 0 then f(x) = v^2 - 4sw = C and should check its sign 1191 CGAL_assertion(sign_A == ZERO); 1192 // (in this case the conic equation represents parabola, 2 parallel lines, 1193 // 1 line or no curve) 1194 Sign sign_B = CGAL_NTS sign(s*(t*v - 4*u)); 1195 if (sign_B != ZERO) 1196 return true; 1197 1198 Sign sign_C = CGAL_NTS sign(v*v - 4*s*w); 1199 1200 return (sign_C != NEGATIVE); 1201 } 1202 1203 // for the test: vertical_ray_shoot_2(const Point_2 & pt,const X_monotone_curve_2 & cv)1204 Point_2 vertical_ray_shoot_2(const Point_2& pt, 1205 const X_monotone_curve_2& cv) const 1206 { 1207 if (cv.is_vertical()) 1208 { 1209 Alg_kernel k; 1210 if (!k.less_y_2_object()(cv.left(), pt)) 1211 return cv.left(); 1212 else 1213 { 1214 CGAL_assertion(k.less_y_2_object()(cv.right(), pt)); 1215 return cv.right(); 1216 } 1217 } 1218 else 1219 return cv.point_at_x(pt); 1220 } 1221 1222 template <class OutputIterator> add_curve_to_output(const Curve_2 & c,OutputIterator oi)1223 OutputIterator add_curve_to_output(const Curve_2& c, OutputIterator oi) const 1224 { 1225 Object objs[2]; 1226 Object* p_obj = this->make_x_monotone_2_object()(c, objs); 1227 for(Object* o = objs; o != p_obj; ++o) 1228 { 1229 X_monotone_curve_2 cv; 1230 if(assign(cv, *o)) 1231 { 1232 *oi++ = make_object(Intersection_curve(cv, 1)); 1233 } 1234 else 1235 { 1236 Point_2 pt; 1237 CGAL_assertion(assign(pt, *o)); 1238 assign(pt, *o); 1239 *oi++ = make_object(pt); 1240 } 1241 } 1242 return oi; 1243 } 1244 1245 /*! Default constructor. */ Env_sphere_traits_3()1246 Env_sphere_traits_3() : m_is_lower(true) 1247 {} 1248 1249 protected: 1250 mutable bool m_is_lower; 1251 }; 1252 1253 /*! 1254 * Compare two spheres: first compare their center points in an 1255 * xyz-lexicographic order, then by their radii. 1256 */ 1257 template <class Kernel> 1258 bool operator< (const CGAL::Sphere_3<Kernel> & a, 1259 const CGAL::Sphere_3<Kernel> & b) 1260 { 1261 Kernel k; 1262 Comparison_result res = k.compare_xyz_3_object()(a.center(), b.center()); 1263 1264 if (res == EQUAL) 1265 { 1266 res = CGAL::compare (a.squared_radius(), b.squared_radius()); 1267 } 1268 1269 return (res == SMALLER); 1270 } 1271 1272 /*! 1273 * Compare two spheres for equality. 1274 */ 1275 template <class Kernel> 1276 bool operator== (const typename Kernel::Sphere_3& a, 1277 const typename Kernel::Sphere_3& b) 1278 { 1279 Kernel k; 1280 1281 if (! k.equal_3_object() (a.center(), b.center())) 1282 return (false); 1283 1284 return (CGAL::compare (a.squared_radius(), b.squared_radius()) == EQUAL); 1285 } 1286 1287 } //namespace CGAL 1288 1289 #endif // ENVELOPE_SPHERES_TRAITS_3_H 1290