|
BWAPI
|
00001 // Copyright (c) 2003,2004,2005,2006 INRIA Sophia-Antipolis (France) and 00002 // Notre Dame University (U.S.A.). All rights reserved. 00003 // 00004 // This file is part of CGAL (www.cgal.org); you may redistribute it under 00005 // the terms of the Q Public License version 1.0. 00006 // See the file LICENSE.QPL distributed with CGAL. 00007 // 00008 // Licensees holding a valid commercial license may use this file in 00009 // accordance with the commercial license agreement provided with the software. 00010 // 00011 // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE 00012 // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. 00013 // 00014 // $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.5-branch/Segment_Delaunay_graph_2/include/CGAL/Segment_Delaunay_graph_2/Voronoi_vertex_sqrt_field_C2.h $ 00015 // $Id: Voronoi_vertex_sqrt_field_C2.h 45094 2008-08-22 16:10:06Z spion $ 00016 // 00017 // 00018 // Author(s) : Menelaos Karavelas <mkaravel@cse.nd.edu> 00019 00020 00021 00022 00023 #ifndef CGAL_SEGMENT_DELAUNAY_GRAPH_2_VORONOI_VERTEX_SQRT_FIELD_C2_H 00024 #define CGAL_SEGMENT_DELAUNAY_GRAPH_2_VORONOI_VERTEX_SQRT_FIELD_C2_H 00025 00026 00027 00028 #include <CGAL/Segment_Delaunay_graph_2/Basic_predicates_C2.h> 00029 #include <CGAL/Segment_Delaunay_graph_2/Are_same_points_C2.h> 00030 #include <CGAL/Segment_Delaunay_graph_2/Are_same_segments_C2.h> 00031 00032 00033 CGAL_BEGIN_NAMESPACE 00034 00035 CGAL_SEGMENT_DELAUNAY_GRAPH_2_BEGIN_NAMESPACE 00036 00037 00038 template<class K> 00039 class Voronoi_vertex_sqrt_field_C2 00040 : public Basic_predicates_C2<K> 00041 { 00042 public: 00043 typedef Basic_predicates_C2<K> Base; 00044 00045 typedef enum {PPP = 0, PPS, PSS, SSS} vertex_t; 00046 struct PPP_Type {}; 00047 struct PPS_Type {}; 00048 struct PSS_Type {}; 00049 struct SSS_Type {}; 00050 00051 typedef typename Base::Point_2 Point_2; 00052 typedef typename Base::Segment_2 Segment_2; 00053 typedef typename Base::Line_2 Line_2; 00054 typedef typename Base::Site_2 Site_2; 00055 typedef typename Base::FT FT; 00056 typedef typename Base::RT RT; 00057 00058 typedef typename Base::Homogeneous_point_2 Homogeneous_point_2; 00059 00060 typedef typename Base::Orientation Orientation; 00061 typedef typename Base::Comparison_result Comparison_result; 00062 typedef typename Base::Oriented_side Oriented_side; 00063 typedef typename Base::Sign Sign; 00064 00065 private: 00066 typedef Are_same_points_C2<K> Are_same_points_2; 00067 typedef Are_same_segments_C2<K> Are_same_segments_2; 00068 00069 typedef typename K::Intersections_tag ITag; 00070 00071 Are_same_points_2 same_points; 00072 Are_same_segments_2 same_segments; 00073 00074 00075 private: 00076 //-------------------------------------------------------------------------- 00077 00078 void 00079 compute_ppp(const Site_2& sp, const Site_2& sq, const Site_2& sr) 00080 { 00081 CGAL_precondition( sp.is_point() && sq.is_point() && 00082 sr.is_point() ); 00083 00084 Point_2 p = sp.point(), q = sq.point(), r = sr.point(); 00085 00086 v_type = PPP; 00087 00088 FT np = CGAL::square(p.x()) + CGAL::square(p.y()); 00089 FT nq = CGAL::square(q.x()) + CGAL::square(q.y()); 00090 FT nr = CGAL::square(r.x()) + CGAL::square(r.y()); 00091 00092 ux = np * (q.y() - r.y()) + nq * (r.y() - p.y()) + nr * (p.y() - q.y()); 00093 uy = -(np * (q.x() - r.x()) + nq * (r.x() - p.x()) + nr * (p.x() - q.x())); 00094 uz = FT(2) * ( (q.x() * r.y() - r.x() * q.y()) + 00095 (r.x() * p.y() - p.x() * r.y()) + 00096 (p.x() * q.y() - q.x() * p.y()) ); 00097 } 00098 00099 //-------------------------------------------------------------------------- 00100 00101 void 00102 compute_pss(const Site_2& p, const Site_2& q, const Site_2& r) 00103 { 00104 CGAL_precondition( p.is_point() && q.is_segment() && 00105 r.is_segment() ); 00106 00107 v_type = PSS; 00108 00109 bool pq = 00110 same_points(p, q.source_site()) || same_points(p, q.target_site()); 00111 bool pr = 00112 same_points(p, r.source_site()) || same_points(p, r.target_site()); 00113 00114 Point_2 pp = p.point(); 00115 00116 if ( pq && pr ) { 00117 ux = pp.x(); 00118 uy = pp.y(); 00119 uz = FT(1); 00120 return; 00121 } 00122 00123 00124 FT a1, b1, c1, a2, b2, c2; 00125 compute_supporting_line(q.supporting_site(), a1, b1, c1); 00126 compute_supporting_line(r.supporting_site(), a2, b2, c2); 00127 00128 FT c1_ = a1 * pp.x() + b1 * pp.y() + c1; 00129 FT c2_ = a2 * pp.x() + b2 * pp.y() + c2; 00130 00131 if ( pq ) { 00132 c1_ = FT(0); 00133 } 00134 00135 if ( pr ) { 00136 c2_ = FT(0); 00137 } 00138 00139 Sign sgn_c1_ = CGAL::sign(c1_); 00140 Sign sgn_c2_ = CGAL::sign(c2_); 00141 00142 if ( sgn_c1_ == NEGATIVE ) { 00143 a1 = -a1; b1 = -b1; c1_ = -c1_; 00144 } else if ( sgn_c1_ == ZERO ) { 00145 00146 CGAL_assertion( pq ); 00147 if ( same_points(p, q.target_site()) ) { 00148 a1 = -a1; b1 = -b1; c1_ = -c1_; 00149 } 00150 } 00151 00152 if ( sgn_c2_ == NEGATIVE ) { 00153 a2 = -a2; b2 = -b2; c2_ = -c2_; 00154 } else if ( sgn_c2_ == ZERO ) { 00155 00156 CGAL_assertion( pr ); 00157 if ( same_points(p, r.source_site()) ) { 00158 a2 = -a2; b2 = -b2; c2_ = -c2_; 00159 } 00160 } 00161 00162 if ( pq ) { 00163 FT J = a1 * c2_; 00164 FT I = b1 * c2_; 00165 00166 FT n1 = CGAL::square(a1) + CGAL::square(b1); 00167 FT n2 = CGAL::square(a2) + CGAL::square(b2); 00168 00169 FT D1D2 = n1 * n2; 00170 00171 uz = -a1 * a2 - b1 * b2 + CGAL::sqrt(D1D2); 00172 00173 ux = J + pp.x() * uz; 00174 uy = I + pp.y() * uz; 00175 00176 } else if ( pr ) { 00177 FT J = a2 * c1_; 00178 FT I = b2 * c1_; 00179 00180 FT n1 = CGAL::square(a1) + CGAL::square(b1); 00181 FT n2 = CGAL::square(a2) + CGAL::square(b2); 00182 00183 FT D1D2 = n1 * n2; 00184 00185 uz = -a1 * a2 - b1 * b2 + CGAL::sqrt(D1D2); 00186 00187 ux = J + pp.x() * uz; 00188 uy = I + pp.y() * uz; 00189 00190 } else { 00191 Line_2 lq(a1, b1, c1_); 00192 Line_2 lr(a2, b2, c2_); 00193 compute_pll(pp, lq, lr); 00194 } 00195 } 00196 00197 00198 void 00199 compute_pll(const Point_2& p, const Line_2& lq, const Line_2& lr) 00200 { 00201 FT a1 = lq.a(), b1 = lq.b(), c1_ = lq.c(); 00202 FT a2 = lr.a(), b2 = lr.b(), c2_ = lr.c(); 00203 00204 CGAL_precondition( c1_ >= FT(0) ); 00205 CGAL_precondition( c2_ >= FT(0) ); 00206 00207 FT n1 = CGAL::square(a1) + CGAL::square(b1); 00208 FT n2 = CGAL::square(a2) + CGAL::square(b2); 00209 00210 FT I = b1 * c2_ + b2 * c1_; 00211 FT J = a1 * c2_ + a2 * c1_; 00212 00213 FT c1c2 = FT(2) * c1_ * c2_; 00214 FT a1a2 = a1 * a2; 00215 FT b1b2 = b1 * b2; 00216 00217 FT D1D2 = n1 * n2; 00218 00219 // compute sigma 00220 FT sigma_expr = b1 * CGAL::sqrt(n2) - b2 * CGAL::sqrt(n1); 00221 Sign s_sigma = CGAL::sign(sigma_expr); 00222 00223 int i_sigma(s_sigma); 00224 FT sigma(i_sigma); 00225 00226 // compute rho 00227 FT rho_expr = a1 * CGAL::sqrt(n2) - a2 * CGAL::sqrt(n1); 00228 Sign s_rho = CGAL::sign(rho_expr); 00229 00230 int i_rho(s_rho); 00231 FT rho(i_rho); 00232 00233 FT sqrt_D1D2 = CGAL::sqrt(D1D2); 00234 00235 FT A = a1a2 - b1b2; 00236 FT u1 = c1c2 * (sqrt_D1D2 + A); 00237 FT u2 = c1c2 * (sqrt_D1D2 - A); 00238 00239 uz = -a1a2 - b1b2 + CGAL::sqrt(D1D2); 00240 00241 ux = J + p.x() * uz + sigma * CGAL::sqrt(u1); 00242 uy = I + p.y() * uz - rho * CGAL::sqrt(u2); 00243 } 00244 00245 00246 //-------------------------------------------------------------------------- 00247 00248 00249 void 00250 compute_pps(const Site_2& p, const Site_2& q, const Site_2& r) 00251 { 00252 CGAL_precondition( p.is_point() && q.is_point() && 00253 r.is_segment() ); 00254 v_type = PPS; 00255 00256 FT a, b, c; 00257 compute_supporting_line(r.supporting_site(), a, b, c); 00258 00259 Point_2 pp = p.point(), qq = q.point(); 00260 00261 FT c_ = a * pp.x() + b * pp.y() + c; 00262 FT cq_ = a * qq.x() + b * qq.y() + c; 00263 00264 if ( same_points(p, r.source_site()) || 00265 same_points(p, r.target_site()) ) { 00266 c_ = FT(0); 00267 } 00268 if ( same_points(q, r.source_site()) || 00269 same_points(q, r.target_site()) ) { 00270 cq_ = FT(0); 00271 } 00272 00273 Sign s = CGAL::sign(c_); 00274 00275 if ( s == NEGATIVE ) { 00276 a = -a; b = -b; c = -c; c_ = -c_; cq_ = -cq_; 00277 } else if ( s == ZERO ) { 00278 Sign s1 = CGAL::sign(cq_); 00279 00280 CGAL_assertion( s1 != ZERO ); 00281 if ( s1 == NEGATIVE ) { 00282 a = -a; b = -b; c = -c; c_ = -c_; cq_ = -cq_; 00283 } 00284 } 00285 00286 FT nl = CGAL::square(a) + CGAL::square(b); 00287 00288 FT x_ = qq.x() - pp.x(); 00289 FT y_ = qq.y() - pp.y(); 00290 FT n_ = CGAL::square(x_) + CGAL::square(y_); 00291 00292 if ( c_ == cq_ ) { 00293 FT e1 = CGAL::square(c_); 00294 FT J = nl * (a * n_ + FT(4) * c_ * x_) - FT(4) * a * e1; 00295 FT I = nl * (b * n_ + FT(4) * c_ * y_) - FT(4) * b * e1; 00296 FT X = FT(8) * nl * c_; 00297 00298 ux = J + pp.x() * X; 00299 uy = I + pp.y() * X; 00300 uz = X; 00301 return; 00302 } 00303 00304 FT e1 = a * x_ + b * y_; 00305 FT e2 = b * x_ - a * y_; 00306 FT e3 = n_ * e1; 00307 FT e4 = FT(2) * c_ * e2; 00308 00309 FT X = FT(2) * CGAL::square(e1); 00310 FT I = b * e3 + x_ * e4; 00311 FT J = a * e3 - y_ * e4; 00312 FT sqrt_S = CGAL::sqrt(n_ * nl * c_ * cq_); 00313 00314 ux = J + pp.x() * X - FT(2) * y_ * sqrt_S; 00315 uy = I + pp.y() * X + FT(2) * x_ * sqrt_S; 00316 uz = X; 00317 } 00318 00319 00320 //-------------------------------------------------------------------------- 00321 bool check_if_exact(const Site_2& , unsigned int , 00322 const Tag_false&) const 00323 { 00324 return true; 00325 } 00326 00327 bool check_if_exact(const Site_2& s, unsigned int i, 00328 const Tag_true&) const 00329 { 00330 return s.is_input(i); 00331 } 00332 00333 // determines of the segment s is on the positive halfspace as 00334 // defined by the supporting line of the segment supp; the line l 00335 // is supposed to be the supporting line of the segment supp and we 00336 // pass it so that we do not have to recompute it 00337 bool 00338 is_on_positive_halfspace(const Site_2& supp, 00339 const Site_2& s, const Line_2& l) const 00340 { 00341 CGAL_precondition( supp.is_segment() && s.is_segment() ); 00342 00343 if ( same_segments(supp.supporting_site(), 00344 s.supporting_site()) ) { 00345 return false; 00346 } 00347 00348 if ( same_points(supp.source_site(), s.source_site()) || 00349 same_points(supp.target_site(), s.source_site()) ) { 00350 return oriented_side_of_line(l, s.target()) == ON_POSITIVE_SIDE; 00351 } 00352 00353 if ( same_points(supp.source_site(), s.target_site()) || 00354 same_points(supp.target_site(), s.target_site()) ) { 00355 return oriented_side_of_line(l, s.source()) == ON_POSITIVE_SIDE; 00356 } 00357 00358 ITag itag; 00359 00360 if ( !check_if_exact(s, 0, itag) && 00361 same_segments(supp.supporting_site(), 00362 s.crossing_site(0)) ) { 00363 return oriented_side_of_line(l, s.target()) == ON_POSITIVE_SIDE; 00364 } 00365 00366 if ( !check_if_exact(s, 1, itag) && 00367 same_segments(supp.supporting_site(), 00368 s.crossing_site(1)) ) { 00369 return oriented_side_of_line(l, s.source()) == ON_POSITIVE_SIDE; 00370 } 00371 00372 return Base::is_on_positive_halfspace(l, s.segment()); 00373 } 00374 00375 //-------------------------------------------------------------------------- 00376 00377 void 00378 orient_lines(const Site_2& sp, const Site_2& sq, 00379 const Site_2& sr, FT a[], FT b[], FT c[]) const 00380 { 00381 CGAL_precondition( sp.is_segment() && sq.is_segment() && 00382 sr.is_segment() ); 00383 00384 Line_2 l[3]; 00385 l[0] = compute_supporting_line(sp.supporting_site()); 00386 l[1] = compute_supporting_line(sq.supporting_site()); 00387 l[2] = compute_supporting_line(sr.supporting_site()); 00388 00389 bool is_oriented[3] = {false, false, false}; 00390 00391 if ( is_on_positive_halfspace(sp, sq, l[0]) || 00392 is_on_positive_halfspace(sp, sr, l[0]) ) { 00393 is_oriented[0] = true; 00394 } else { 00395 00396 l[0] = opposite_line(l[0]); 00397 if ( is_on_positive_halfspace(sp, sq, l[0]) || 00398 is_on_positive_halfspace(sp, sr, l[0]) ) { 00399 is_oriented[0] = true; 00400 } else { 00401 l[0] = opposite_line(l[0]); 00402 } 00403 } 00404 00405 if ( is_on_positive_halfspace(sq, sp, l[1]) || 00406 is_on_positive_halfspace(sq, sr, l[1]) ) { 00407 is_oriented[1] = true; 00408 } else { 00409 l[1] = opposite_line(l[1]); 00410 if ( is_on_positive_halfspace(sq, sp, l[1]) || 00411 is_on_positive_halfspace(sq, sr, l[1]) ) { 00412 is_oriented[1] = true; 00413 } else { 00414 l[1] = opposite_line(l[1]); 00415 } 00416 } 00417 00418 if ( is_on_positive_halfspace(sr, sp, l[2]) || 00419 is_on_positive_halfspace(sr, sq, l[2]) ) { 00420 is_oriented[2] = true; 00421 } else { 00422 l[2] = opposite_line(l[2]); 00423 if ( is_on_positive_halfspace(sr, sp, l[2]) || 00424 is_on_positive_halfspace(sr, sq, l[2]) ) { 00425 is_oriented[2] = true; 00426 } else { 00427 l[2] = opposite_line(l[2]); 00428 } 00429 } 00430 00431 if ( is_oriented[0] && is_oriented[1] && is_oriented[2] ) { 00432 for (int i = 0; i < 3; i++) { 00433 a[i] = l[i].a(); 00434 b[i] = l[i].b(); 00435 c[i] = l[i].c(); 00436 } 00437 return; 00438 } 00439 00440 int i_no(-1); 00441 for (int i = 0; i < 3; i++) { 00442 if ( !is_oriented[i] ) { 00443 i_no = i; 00444 CGAL_assertion( is_oriented[(i+1)%3] && is_oriented[(i+2)%3] ); 00445 break; 00446 } 00447 } 00448 00449 CGAL_assertion( i_no != -1 ); 00450 00451 FT sqrt_d[3]; 00452 for (int i = 0; i < 3; i++) { 00453 FT d1 = CGAL::square(l[i].a()) + CGAL::square(l[i].b()); 00454 sqrt_d[i] = CGAL::sqrt(d1); 00455 } 00456 00457 FT z[3]; 00458 for (int i = 0; i < 3; i++) { 00459 z[i] = l[(i+1)%3].a() * l[(i+2)%3].b() 00460 - l[(i+2)%3].a() * l[(i+1)%3].b(); 00461 } 00462 00463 00464 FT vz = z[0] * sqrt_d[0] + z[1] * sqrt_d[1] + z[2] * sqrt_d[2]; 00465 00466 Sign s_minus_vz = CGAL::sign(vz); 00467 00468 CGAL_assertion( s_minus_vz != ZERO ); 00469 00470 if ( s_minus_vz == NEGATIVE ) { 00471 l[i_no] = opposite_line(l[i_no]); 00472 00473 for (int i = 0; i < 3; i++) { 00474 a[i] = l[i].a(); 00475 b[i] = l[i].b(); 00476 c[i] = l[i].c(); 00477 } 00478 return; 00479 } 00480 00481 // now we have to check if the other orientation of l[i_no] 00482 // corresponds to a CCW triangle as well. 00483 z[(i_no+1)%3] = -z[(i_no+1)%3]; 00484 z[(i_no+2)%3] = -z[(i_no+2)%3]; 00485 00486 vz = z[0] * sqrt_d[0] + z[1] * sqrt_d[1] + z[2] * sqrt_d[2]; 00487 00488 Sign s_minus_vz_2 = CGAL::sign(vz); 00489 00490 CGAL_assertion( s_minus_vz_2 != ZERO ); 00491 00492 if ( s_minus_vz_2 == NEGATIVE ) { 00493 // the other orientation does not correspond to a CCW triangle. 00494 for (int i = 0; i < 3; i++) { 00495 a[i] = l[i].a(); 00496 b[i] = l[i].b(); 00497 c[i] = l[i].c(); 00498 } 00499 return; 00500 } 00501 00502 // now compute the Voronoi vertex; 00503 for (int i = 0; i < 3; i++) { 00504 a[i] = l[i].a(); 00505 b[i] = l[i].b(); 00506 c[i] = l[i].c(); 00507 } 00508 00509 FT x[3], y[3], w[3]; 00510 for (int i = 0; i < 3; i++) { 00511 x[i] = c[(i+1)%3] * b[(i+2)%3] - c[(i+2)%3] * b[(i+1)%3]; 00512 y[i] = -(c[(i+1)%3] * a[(i+2)%3] - c[(i+2)%3] * a[(i+1)%3]); 00513 w[i] = -(a[(i+1)%3] * b[(i+2)%3] - a[(i+2)%3] * b[(i+1)%3]); 00514 } 00515 00516 FT vx = FT(0), vy = FT(0), vw = FT(0); 00517 00518 for (int i = 0; i < 3; i++) { 00519 vx += x[i] * sqrt_d[i]; 00520 vy += y[i] * sqrt_d[i]; 00521 vw += w[i] * sqrt_d[i]; 00522 } 00523 00524 Sign s_vw = CGAL::sign(vw); 00525 00526 FT dist = 00527 a[(i_no+1)%3] * vx + b[(i_no+1)%3] * vy + c[(i_no+1)%3] * vw; 00528 00529 00530 Sign sgn_dist = s_vw * CGAL::sign(dist); 00531 00532 CGAL_assertion( sgn_dist != ZERO ); 00533 00534 if ( sgn_dist == NEGATIVE ) { 00535 a[i_no] = -a[i_no]; 00536 b[i_no] = -b[i_no]; 00537 c[i_no] = -c[i_no]; 00538 } 00539 } 00540 00541 00542 00543 void 00544 compute_sss(const Site_2& p, const Site_2& q, const Site_2& r) 00545 { 00546 CGAL_precondition( p.is_segment() && q.is_segment() && 00547 r.is_segment() ); 00548 00549 v_type = SSS; 00550 00551 FT a[3], b[3], c[3]; 00552 FT cx[3], cy[3], cz[3], sqrt_D[3]; 00553 00554 orient_lines(p, q, r, a, b, c); 00555 00556 for (int i = 0; i < 3; i++) { 00557 cx[i] = c[(i+1)%3] * b[(i+2)%3] - c[(i+2)%3] * b[(i+1)%3]; 00558 cy[i] = -(c[(i+1)%3] * a[(i+2)%3] - c[(i+2)%3] * a[(i+1)%3]); 00559 cz[i] = -(a[(i+1)%3] * b[(i+2)%3] - a[(i+2)%3] * b[(i+1)%3]); 00560 00561 FT d = CGAL::square(a[i]) + CGAL::square(b[i]); 00562 sqrt_D[i] = CGAL::sqrt(d); 00563 } 00564 00565 ux = cx[0] * sqrt_D[0] + cx[1] * sqrt_D[1] + cx[2] * sqrt_D[2]; 00566 uy = cy[0] * sqrt_D[0] + cy[1] * sqrt_D[1] + cy[2] * sqrt_D[2]; 00567 uz = cz[0] * sqrt_D[0] + cz[1] * sqrt_D[1] + cz[2] * sqrt_D[2]; 00568 } 00569 00570 //-------------------------------------------------------------------------- 00571 00572 void 00573 compute_vertex(const Site_2& s1, const Site_2& s2, const Site_2& s3) 00574 { 00575 if ( s1.is_point() && s2.is_point() && s3.is_point() ) { 00576 compute_ppp(s1, s2, s3); 00577 00578 } else if ( s1.is_segment() && s2.is_point() && s3.is_point() ) { 00579 compute_vertex(s2, s3, s1); 00580 pps_idx = 1; 00581 00582 } else if ( s1.is_point() && s2.is_segment() && s3.is_point() ) { 00583 compute_vertex(s3, s1, s2); 00584 pps_idx = 2; 00585 00586 } else if ( s1.is_point() && s2.is_point() && s3.is_segment() ) { 00587 compute_pps(s1, s2, s3); 00588 pps_idx = 0; 00589 00590 } else if ( s1.is_point() && s2.is_segment() && s3.is_segment() ) { 00591 compute_pss(s1, s2, s3); 00592 } else if ( s1.is_segment() && s2.is_point() && s3.is_segment() ) { 00593 compute_vertex(s2, s3, s1); 00594 } else if ( s1.is_segment() && s2.is_segment() && s3.is_point() ) { 00595 compute_vertex(s3, s1, s2); 00596 } else { 00597 compute_sss(s1, s2, s3); 00598 } 00599 00600 } 00601 00602 //-------------------------------------------------------------------------- 00603 00604 bool is_endpoint_of(const Site_2& p, const Site_2& s) const 00605 { 00606 CGAL_precondition( p.is_point() && s.is_segment() ); 00607 return ( same_points(p, s.source_site()) || 00608 same_points(p, s.target_site()) ); 00609 } 00610 00611 00612 //-------------------------------------------------------------------------- 00613 //-------------------------------------------------------------------------- 00614 // the incircle test 00615 //-------------------------------------------------------------------------- 00616 //-------------------------------------------------------------------------- 00617 00618 //-------------------------------------------------------------------------- 00619 // the incircle test when the fourth site is a point 00620 //-------------------------------------------------------------------------- 00621 00622 //-------------------------------------------------------------------------- 00623 00624 Sign 00625 check_easy_degeneracies(const Site_2& t, PPS_Type, 00626 bool& use_result) const 00627 { 00628 CGAL_precondition( t.is_point() ); 00629 00630 use_result = false; 00631 if ( ( p_.is_point() && same_points(p_, t) ) || 00632 ( q_.is_point() && same_points(q_, t) ) || 00633 ( r_.is_point() && same_points(r_, t) ) ) { 00634 use_result = true; 00635 return ZERO; 00636 } 00637 00638 if ( ( p_.is_segment() && is_endpoint_of(t, p_) ) || 00639 ( q_.is_segment() && is_endpoint_of(t, q_) ) || 00640 ( r_.is_segment() && is_endpoint_of(t, r_) ) ) { 00641 use_result = true; 00642 return POSITIVE; 00643 } 00644 00645 return ZERO; 00646 } 00647 00648 Sign 00649 check_easy_degeneracies(const Site_2& t, PSS_Type, 00650 bool& use_result) const 00651 { 00652 CGAL_precondition( t.is_point() ); 00653 00654 return check_easy_degeneracies(t, PPS_Type(), use_result); 00655 } 00656 00657 Sign 00658 check_easy_degeneracies(const Site_2& t, SSS_Type, 00659 bool& use_result) const 00660 { 00661 CGAL_precondition( t.is_point() ); 00662 00663 use_result = false; 00664 00665 // ADD THE CASES WHERE t IS AN ENDPOINT OF ONE OF THE SEGMENTS 00666 return ZERO; 00667 } 00668 00669 //-------------------------------------------------------------------------- 00670 00671 Sign incircle_p(const Site_2& st, PPP_Type) const 00672 { 00673 CGAL_precondition( st.is_point() ); 00674 00675 Point_2 t = st.point(); 00676 00677 Oriented_side os = 00678 side_of_oriented_circle(p_.point(), q_.point(), r_.point(), t); 00679 if ( os == ON_POSITIVE_SIDE ) { return NEGATIVE; } 00680 if ( os == ON_NEGATIVE_SIDE ) { return POSITIVE; } 00681 return ZERO; 00682 } 00683 00684 //-------------------------------------------------------------------------- 00685 00686 template<class Type> 00687 inline 00688 Sign incircle_p(const Site_2& st, Type type) const 00689 { 00690 CGAL_precondition( st.is_point() ); 00691 00692 bool use_result(false); 00693 Sign s = check_easy_degeneracies(st, type, use_result); 00694 if ( use_result ) { return s; } 00695 00696 return incircle_p_no_easy(st, type); 00697 } 00698 00699 template<class Type> 00700 Sign incircle_p_no_easy(const Site_2& st, Type) const 00701 { 00702 CGAL_precondition( st.is_point() ); 00703 00704 FT r2 = squared_radius(); 00705 00706 Point_2 t = st.point(); 00707 00708 FT d2 = CGAL::square(x() - t.x()) + 00709 CGAL::square(y() - t.y()); 00710 00711 return CGAL::compare(d2, r2); 00712 } 00713 00714 //-------------------------------------------------------------------------- 00715 00716 Sign incircle_p(const Site_2& t) const 00717 { 00718 if ( is_degenerate_Voronoi_circle() ) { 00719 return POSITIVE; 00720 } 00721 00722 Sign s(ZERO); 00723 switch ( v_type ) { 00724 case PPP: 00725 s = incircle_p(t, PPP_Type()); 00726 break; 00727 case PPS: 00728 s = incircle_p(t, PPS_Type()); 00729 break; 00730 case PSS: 00731 s = incircle_p(t, PSS_Type()); 00732 break; 00733 case SSS: 00734 s = incircle_p(t, SSS_Type()); 00735 break; 00736 } 00737 00738 return s; 00739 } 00740 00741 Sign incircle_p_no_easy(const Site_2& t) const 00742 { 00743 Sign s(ZERO); 00744 switch ( v_type ) { 00745 case PPP: 00746 s = incircle_p(t, PPP_Type()); 00747 break; 00748 case PPS: 00749 s = incircle_p_no_easy(t, PPS_Type()); 00750 break; 00751 case PSS: 00752 s = incircle_p_no_easy(t, PSS_Type()); 00753 break; 00754 case SSS: 00755 s = incircle_p_no_easy(t, SSS_Type()); 00756 break; 00757 } 00758 00759 return s; 00760 } 00761 00762 //-------------------------------------------------------------------------- 00763 00764 //-------------------------------------------------------------------------- 00765 // the incircle test when the fourth site is a segment 00766 //-------------------------------------------------------------------------- 00767 00768 //-------------------------------------------------------------------------- 00769 00770 00771 Oriented_side 00772 oriented_side(const Line_2& l, const Point_2& p) const 00773 { 00774 Line_2 l1(l.b(), -l.a(), l.a() * y() - l.b() * x()); 00775 00776 return oriented_side_of_line(l1, p); 00777 } 00778 00779 00780 Sign incircle(const Line_2& l) const 00781 { 00782 FT r2 = squared_radius(); 00783 00784 FT n2 = CGAL::square(l.a()) + CGAL::square(l.b()); 00785 00786 FT d2 = CGAL::square(l.a() * x() + l.b() * y() + l.c()); 00787 00788 return CGAL::compare(d2, r2 * n2); 00789 } 00790 00791 00792 00793 //-------------------------------------------------------------------------- 00794 //-------------------------------------------------------------------------- 00795 00796 Sign incircle_s(const Site_2& t, int i) const 00797 { 00798 CGAL_precondition( t.is_segment() ); 00799 00800 if ( v_type == PPP || v_type == PPS ) { 00801 if ( p_.is_point() && q_.is_point() && 00802 is_endpoint_of(p_, t) && is_endpoint_of(q_, t) ) { 00803 return NEGATIVE; 00804 } 00805 00806 if ( p_.is_point() && r_.is_point() && 00807 is_endpoint_of(p_, t) && is_endpoint_of(r_, t) ){ 00808 return NEGATIVE; 00809 } 00810 00811 if ( q_.is_point() && r_.is_point() && 00812 is_endpoint_of(q_, t) && is_endpoint_of(r_, t) ){ 00813 return NEGATIVE; 00814 } 00815 } 00816 00817 if ( v_type == PSS ) { 00818 if ( p_.is_segment() && 00819 same_segments(p_.supporting_site(), 00820 t.supporting_site()) ) { 00821 return POSITIVE; 00822 } 00823 if ( q_.is_segment() && 00824 same_segments(q_.supporting_site(), 00825 t.supporting_site()) ) { 00826 return POSITIVE; 00827 } 00828 if ( r_.is_segment() && 00829 same_segments(r_.supporting_site(), 00830 t.supporting_site()) ) { 00831 return POSITIVE; 00832 } 00833 } 00834 00835 return incircle_s_no_easy(t, i); 00836 } 00837 00838 00839 Sign incircle_s_no_easy(const Site_2& t, int) const 00840 { 00841 Sign d1, d2; 00842 if ( ( p_.is_point() && same_points(p_, t.source_site()) ) || 00843 ( q_.is_point() && same_points(q_, t.source_site()) ) || 00844 ( r_.is_point() && same_points(r_, t.source_site()) ) ) { 00845 d1 = ZERO; 00846 } else { 00847 d1 = incircle_p(t.source_site()); 00848 } 00849 //if ( d1 == NEGATIVE ) { return NEGATIVE; } 00850 if ( certainly(d1 == NEGATIVE) ) { return NEGATIVE; } 00851 if (! is_certain(d1 == NEGATIVE) ) { return indeterminate<Sign>(); } 00852 00853 if ( ( p_.is_point() && same_points(p_, t.target_site()) ) || 00854 ( q_.is_point() && same_points(q_, t.target_site()) ) || 00855 ( r_.is_point() && same_points(r_, t.target_site()) ) ) { 00856 d2 = ZERO; 00857 } else { 00858 d2 = incircle_p(t.target_site()); 00859 } 00860 //if ( d2 == NEGATIVE ) { return NEGATIVE; } 00861 if (certainly( d2 == NEGATIVE ) ) { return NEGATIVE; } 00862 if (! is_certain( d2 == NEGATIVE ) ) { return indeterminate<Sign>(); } 00863 00864 Line_2 l = compute_supporting_line(t.supporting_site()); 00865 Sign sl = incircle(l); 00866 00867 //if ( sl == POSITIVE ) { return sl; } 00868 if (certainly( sl == POSITIVE )) { return sl; } 00869 if (! is_certain( sl == POSITIVE )) { return indeterminate<Sign>(); } 00870 00871 if ( sl == ZERO && (d1 == ZERO || d2 == ZERO) ) { return ZERO; } 00872 00873 Oriented_side os1 = oriented_side(l, t.source()); 00874 Oriented_side os2 = oriented_side(l, t.target()); 00875 00876 if ( sl == ZERO ) { 00877 if (os1 == ON_ORIENTED_BOUNDARY || os2 == ON_ORIENTED_BOUNDARY) { 00878 return ZERO; 00879 } 00880 return ( os1 == os2 ) ? POSITIVE : ZERO; 00881 } 00882 00883 return (os1 == os2) ? POSITIVE : NEGATIVE; 00884 } 00885 00886 //-------------------------------------------------------------------------- 00887 00888 Sign incircle_s(const Site_2& t) const 00889 { 00890 CGAL_precondition( t.is_segment() ); 00891 00892 if ( is_degenerate_Voronoi_circle() ) { 00893 // case 1: the new segment is not adjacent to the center of the 00894 // degenerate Voronoi circle 00895 if ( !same_points( p_ref(), t.source_site() ) && 00896 !same_points( p_ref(), t.target_site() ) ) { 00897 return POSITIVE; 00898 } 00899 00900 CGAL_assertion( v_type == PSS ); 00901 00902 if ( p_.is_segment() && 00903 same_segments(p_.supporting_site(), 00904 t.supporting_site()) ) { 00905 return ZERO; 00906 } 00907 00908 if ( q_.is_segment() && 00909 same_segments(q_.supporting_site(), 00910 t.supporting_site()) ) { 00911 return ZERO; 00912 } 00913 00914 if ( r_.is_segment() && 00915 same_segments(r_.supporting_site(), 00916 t.supporting_site()) ) { 00917 return ZERO; 00918 } 00919 00920 Site_2 pr; 00921 Site_2 sp, sq; 00922 if ( p_.is_point() ) { 00923 CGAL_assertion( q_.is_segment() && r_.is_segment() ); 00924 pr = p_; 00925 sp = q_; 00926 sq = r_; 00927 } else if ( q_.is_point() ) { 00928 CGAL_assertion( r_.is_segment() && p_.is_segment() ); 00929 pr = q_; 00930 sp = r_; 00931 sq = p_; 00932 } else { 00933 CGAL_assertion( p_.is_segment() && q_.is_segment() ); 00934 pr = r_; 00935 sp = p_; 00936 sq = q_; 00937 } 00938 00939 Point_2 pq = sq.source(), pp = sp.source(), pt = t.source(); 00940 00941 if ( same_points(sp.source_site(), pr) ) { pp = sp.target(); } 00942 if ( same_points(sq.source_site(), pr) ) { pq = sq.target(); } 00943 if ( same_points( t.source_site(), pr) ) { pt = t.target(); } 00944 00945 Point_2 pr_ = pr.point(); 00946 00947 if ( CGAL::orientation(pr_, pp, pt) == LEFT_TURN && 00948 CGAL::orientation(pr_, pq, pt) == RIGHT_TURN ) { 00949 return NEGATIVE; 00950 } 00951 return ZERO; 00952 } // if ( is_degenerate_Voronoi_circle() ) 00953 00954 Sign s = incircle_s(t, 0); 00955 00956 return s; 00957 } 00958 00959 inline 00960 Sign incircle_s_no_easy(const Site_2& t) const 00961 { 00962 return incircle_s_no_easy(t, 0); 00963 } 00964 00965 //-------------------------------------------------------------------------- 00966 // subpredicates for the incircle test 00967 //-------------------------------------------------------------------------- 00968 00969 00970 public: 00971 bool is_degenerate_Voronoi_circle() const 00972 { 00973 if ( v_type != PSS ) { return false; } 00974 00975 if ( p_.is_point() ) { 00976 return ( is_endpoint_of(p_, q_) && is_endpoint_of(p_, r_) ); 00977 } else if ( q_.is_point() ) { 00978 return ( is_endpoint_of(q_, p_) && is_endpoint_of(q_, r_) ); 00979 } else { 00980 CGAL_assertion( r_.is_point() ); 00981 return ( is_endpoint_of(r_, p_) && is_endpoint_of(r_, q_) ); 00982 } 00983 } 00984 00985 00986 //-------------------------------------------------------------------------- 00987 00988 private: 00989 00990 //-------------------------------------------------------------------------- 00991 // the reference point (valid if v_type != SSS) 00992 //-------------------------------------------------------------------------- 00993 00994 Site_2 p_ref() const 00995 { 00996 CGAL_precondition ( v_type != SSS ); 00997 00998 if ( v_type == PPS ) { 00999 if ( pps_idx == 0 ) { return p_; } 01000 if ( pps_idx == 1 ) { return q_; } 01001 return r_; 01002 } 01003 01004 if ( p_.is_point() ) { 01005 return p_; 01006 } else if ( q_.is_point() ) { 01007 return q_; 01008 } else { 01009 CGAL_assertion( r_.is_point() ); 01010 return r_; 01011 } 01012 } 01013 01014 01015 public: 01016 //-------------------------------------------------------------------------- 01017 //-------------------------------------------------------------------------- 01018 // access methods 01019 //-------------------------------------------------------------------------- 01020 //-------------------------------------------------------------------------- 01021 01022 01023 FT x() const { return hx() / hw(); } 01024 FT y() const { return hy() / hw(); } 01025 01026 FT hx() const { 01027 return ux; 01028 } 01029 01030 FT hy() const { 01031 return uy; 01032 } 01033 01034 FT hw() const { 01035 return uz; 01036 } 01037 01038 FT squared_radius() const { 01039 switch (v_type) { 01040 case PPP: case PPS: case PSS: 01041 { 01042 Point_2 pref = p_ref().point(); 01043 FT dx2 = CGAL::square(x() - pref.x()); 01044 FT dy2 = CGAL::square(y() - pref.y()); 01045 return dx2 + dy2; 01046 } 01047 case SSS: 01048 { 01049 Line_2 l = compute_supporting_line(p_.supporting_site()); 01050 Homogeneous_point_2 q = compute_projection(l, point()); 01051 01052 FT dx2 = CGAL::square(x() - q.x()); 01053 FT dy2 = CGAL::square(y() - q.y()); 01054 return dx2 + dy2; 01055 } 01056 default: 01057 return FT(0); 01058 } 01059 } 01060 01061 01062 Point_2 point() const { 01063 if ( is_degenerate_Voronoi_circle() ) { 01064 return degenerate_point(); 01065 } 01066 01067 return Point_2(x(), y()); 01068 } 01069 01070 01071 Point_2 degenerate_point() const 01072 { 01073 CGAL_precondition( is_degenerate_Voronoi_circle() ); 01074 return p_ref().point(); 01075 } 01076 01077 01078 typename K::Circle_2 circle() const 01079 { 01080 typedef typename K::Circle_2 K_Circle_2; 01081 return K_Circle_2(point(), squared_radius()); 01082 } 01083 01084 vertex_t type() const { return v_type; } 01085 01086 public: 01087 Voronoi_vertex_sqrt_field_C2(const Site_2& p, 01088 const Site_2& q, 01089 const Site_2& r) 01090 : p_(p), q_(q), r_(r) 01091 { 01092 compute_vertex(p, q, r); 01093 } 01094 01095 //-------------------------------------------------------------------------- 01096 01097 Sign incircle(const Site_2& t) const 01098 { 01099 if ( t.is_point() ) { 01100 return incircle_p(t); 01101 } 01102 return incircle_s(t); 01103 } 01104 01105 Sign incircle_no_easy(const Site_2& t) const 01106 { 01107 Sign s; 01108 01109 if ( t.is_point() ) { 01110 s = incircle_p_no_easy(t); 01111 } else { 01112 s = incircle_s_no_easy(t); 01113 } 01114 01115 return s; 01116 } 01117 01118 //-------------------------------------------------------------------------- 01119 01120 01121 Orientation orientation(const Line_2& l) const 01122 { 01123 return CGAL::sign(l.a() * x() + l.b() * y() + l.c()); 01124 } 01125 01126 Oriented_side oriented_side(const Line_2& l) const 01127 { 01128 return orientation(l); 01129 } 01130 01131 //-------------------------------------------------------------------------- 01132 01133 private: 01134 const Site_2& p_, q_, r_; 01135 01136 vertex_t v_type; 01137 01138 // index that indicates the refence point for the case PPS 01139 short pps_idx; 01140 01141 FT ux, uy, uz; 01142 }; 01143 01144 CGAL_SEGMENT_DELAUNAY_GRAPH_2_END_NAMESPACE 01145 01146 CGAL_END_NAMESPACE 01147 01148 #endif // CGAL_SEGMENT_DELAUNAY_GRAPH_2_VORONOI_VEFTEX_SQRT_FIELD_C2_H
1.7.6.1
