|
BWAPI
|
00001 // Copyright (c) 2000,2001 Utrecht University (The Netherlands), 00002 // ETH Zurich (Switzerland), Freie Universitaet Berlin (Germany), 00003 // INRIA Sophia-Antipolis (France), Martin-Luther-University Halle-Wittenberg 00004 // (Germany), Max-Planck-Institute Saarbruecken (Germany), RISC Linz (Austria), 00005 // and Tel-Aviv University (Israel). All rights reserved. 00006 // 00007 // This file is part of CGAL (www.cgal.org); you can redistribute it and/or 00008 // modify it under the terms of the GNU Lesser General Public License as 00009 // published by the Free Software Foundation; version 2.1 of the License. 00010 // See the file LICENSE.LGPL distributed with CGAL. 00011 // 00012 // Licensees holding a valid commercial license may use this file in 00013 // accordance with the commercial license agreement provided with the software. 00014 // 00015 // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE 00016 // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. 00017 // 00018 // $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.5-branch/Conic_2/include/CGAL/ConicHPA2.h $ 00019 // $Id: ConicHPA2.h 34970 2006-10-28 13:07:32Z hemmer $ 00020 // 00021 // 00022 // Author(s) : Bernd Gaertner, Sven Schoenherr <sven@inf.ethz.ch> 00023 00024 #ifndef CGAL_CONICHPA2_H 00025 #define CGAL_CONICHPA2_H 00026 00027 // includes 00028 #include <CGAL/Conic_misc.h> 00029 #include <CGAL/kernel_assertions.h> 00030 00031 CGAL_BEGIN_NAMESPACE 00032 00033 template < class PT, class DA> 00034 class ConicHPA2; 00035 00036 template < class PT, class DA> 00037 class _Min_ellipse_2_adapterH2__Ellipse; 00038 00039 00040 template < class _PT, class _DA> 00041 class ConicHPA2 00042 { 00043 public: 00044 // types 00045 typedef _PT PT; 00046 typedef _DA DA; 00047 typedef typename _DA::RT RT; 00048 00049 //private: 00050 //friend class Conic_2< CGAL::Homogeneous<RT> >; 00051 friend class _Min_ellipse_2_adapterH2__Ellipse<PT,DA>; 00052 00053 DA dao; 00054 RT _r, _s, _t, _u, _v, _w; 00055 Conic_type type; 00056 CGAL::Orientation o; 00057 bool empty, trivial, degenerate; 00058 00059 00060 void 00061 set_linear_combination (const RT& a1, const ConicHPA2<PT,DA>& c1, 00062 const RT& a2, const ConicHPA2<PT,DA>& c2) 00063 { 00064 _r = a1 * c1.r() + a2 * c2.r(); 00065 _s = a1 * c1.s() + a2 * c2.s(); 00066 _t = a1 * c1.t() + a2 * c2.t(); 00067 _u = a1 * c1.u() + a2 * c2.u(); 00068 _v = a1 * c1.v() + a2 * c2.v(); 00069 _w = a1 * c1.w() + a2 * c2.w(); 00070 } 00071 00072 static void set_two_linepairs (const PT& p1, 00073 const PT& p2, 00074 const PT& p3, 00075 const PT& p4, 00076 ConicHPA2<PT,DA>& pair1, 00077 ConicHPA2<PT,DA>& pair2) 00078 { 00079 RT x1, y1, h1, x2, y2, h2, x3, y3, h3, x4, y4, h4; 00080 const DA& da = pair1.da(); 00081 da.get (p1, x1, y1, h1); 00082 da.get (p2, x2, y2, h2); 00083 da.get (p3, x3, y3, h3); 00084 da.get (p4, x4, y4, h4); 00085 00086 CGAL::Orientation side1_24 = (CGAL::Orientation)(CGAL_NTS sign 00087 (-h2*x1*y4+h1*x2*y4 00088 +h2*x4*y1-h4*x2*y1 00089 +h4*x1*y2-h1*x4*y2)), 00090 side3_24 = (CGAL::Orientation)(CGAL_NTS sign 00091 (-h2*x3*y4+h3*x2*y4 00092 +h2*x4*y3-h4*x2*y3 00093 +h4*x3*y2-h3*x4*y2)); 00094 if (side1_24 != side3_24) { 00095 // (counter)clockwise order 00096 pair1.set_linepair (p1, p2, p3, p4); 00097 pair2.set_linepair (p2, p3, p4, p1); 00098 } else { 00099 CGAL::Orientation side1_32 = (CGAL::Orientation)(CGAL_NTS sign 00100 (-h3*x1*y2+h1*x3*y2 00101 +h3*x2*y1-h2*x3*y1 00102 +h2*x1*y3-h1*x2*y3)); 00103 if (side1_32 != side3_24) { 00104 // p1, p2 need to be swapped 00105 pair1.set_linepair (p2, p1, p3, p4); 00106 pair2.set_linepair (p1, p3, p4, p2); 00107 } else { 00108 // p2, p3 need to be swapped 00109 pair1.set_linepair (p1, p3, p2, p4); 00110 pair2.set_linepair (p3, p2, p4, p1); 00111 } 00112 } 00113 } 00114 00115 void set_ellipse (const ConicHPA2<PT,DA>& pair1, 00116 const ConicHPA2<PT,DA>& pair2) 00117 { 00118 RT b = RT(2) * (pair1.r() * pair2.s() + pair1.s() * pair2.r()) - 00119 pair1.t() * pair2.t(); 00120 set_linear_combination (pair2.det()-b, pair1, 00121 pair1.det()-b, pair2); 00122 } 00123 00124 void set (const ConicHPA2<PT,DA>& c1, 00125 const ConicHPA2<PT,DA>& c2, 00126 const PT& p) 00127 { 00128 set_linear_combination (c2.evaluate(p), c1, -c1.evaluate(p), c2); 00129 } 00130 00131 CGAL::Sign vol_derivative (RT dr, RT ds, RT dt, 00132 RT du, RT dv, RT dw) const 00133 { 00134 RT a1 = RT(4)*r()*ds+RT(4)*dr*s()-RT(2)*t()*dt, 00135 a0 = RT(4)*r()*s()-t()*t(), 00136 b1 = (RT(4)*r()*s()-t()*t())*dw+(RT(4)*r()*ds+RT(4)*dr*s()- 00137 RT(2)*t()*dt)*w()-u()*u()*ds - 00138 RT(2)*u()*du*s()-v()*v()*dr-RT(2)*v()*dv*r()+u()*v()*dt+ 00139 (u()*dv+du*v())*t(), 00140 b0 = (RT(4)*r()*s()-t()*t())*w() 00141 -u()*u()*s()-v()*v()*r()+u()*v()*t(), 00142 c0 = -RT(2)*a0*b1 + RT(3)*a1*b0; 00143 00144 return CGAL_NTS sign ((int)-CGAL_NTS sign (c0)*o); 00145 } 00146 00147 double vol_minimum (RT dr, RT ds, RT dt, RT du, RT dv, RT dw) const 00148 { 00149 RT a2 = RT(4)*dr*ds-dt*dt, 00150 a1 = RT(4)*r()*ds+RT(4)*dr*s()-RT(2)*t()*dt, 00151 a0 = RT(4)*r()*s()-t()*t(), 00152 b3 = (RT(4)*dr*ds-dt*dt)*dw-du*du*ds-dv*dv*dr+du*dv*dt, 00153 b2 = (RT(4)*r()*ds+RT(4)*dr*s()-RT(2)*t()*dt)*dw+ 00154 (RT(4)*dr*ds-dt*dt)*w()-RT(2)*u()*du*ds-du*du*s()- 00155 RT(2)*v()*dv*dr-dv*dv*r()+(u()*dv+du*v())*dt+du*dv*t(), 00156 b1 = (RT(4)*r()*s()-t()*t())*dw+(RT(4)*r()*ds+RT(4)*dr*s()- 00157 RT(2)*t()*dt)*w()-u()*u()*ds - 00158 RT(2)*u()*du*s()-v()*v()*dr-RT(2)*v()*dv*r()+u()*v()*dt+ 00159 (u()*dv+du*v())*t(), 00160 b0 = (RT(4)*r()*s()-t()*t())*w() 00161 -u()*u()*s()-v()*v()*r()+u()*v()*t(), 00162 c3 = -RT(3)*a1*b3 + RT(2)*a2*b2, 00163 c2 = -RT(6)*a0*b3 - a1*b2 + RT(4)*a2*b1, 00164 c1 = -RT(4)*a0*b2 + a1*b1 + RT(6)*a2*b0, 00165 c0 = -RT(2)*a0*b1 + RT(3)*a1*b0; 00166 00167 double roots[3]; 00168 int nr_roots = solve_cubic 00169 (CGAL::to_double(c3), CGAL::to_double(c2), 00170 CGAL::to_double(c1), CGAL::to_double(c0), 00171 roots[0], roots[1], roots[2]); 00172 CGAL_kernel_precondition (nr_roots > 0); // minimum exists 00173 return best_value (roots, nr_roots, 00174 CGAL::to_double(a2), CGAL::to_double(a1), 00175 CGAL::to_double(a0), CGAL::to_double(b3), 00176 CGAL::to_double(b2), CGAL::to_double(b1), 00177 CGAL::to_double(b0)); 00178 } 00179 00180 00181 00182 protected: 00183 RT det () const 00184 { 00185 return RT(4)*s()*r() - t()*t(); 00186 } 00187 00188 void analyse( ) 00189 { 00190 RT d = det(); 00191 type = (Conic_type)(CGAL_NTS sign(d)); 00192 switch (type) { 00193 case HYPERBOLA: 00194 { 00195 trivial = empty = false; 00196 RT z_prime = d*w() - u()*u()*s() - v()*v()*r() + u()*v()*t(); 00197 o = (CGAL::Orientation)(CGAL_NTS sign (z_prime)); 00198 degenerate = (o == CGAL::ZERO); 00199 00200 00201 } 00202 break; 00203 case PARABOLA: 00204 { 00205 if (!CGAL_NTS is_zero (r())) { 00206 trivial = false; 00207 degenerate = (t()*u() == RT(2)*r()*v()); 00208 if (degenerate) { 00209 CGAL::Sign discr = (CGAL::Sign) 00210 CGAL_NTS sign(u()*u()-RT(4)*r()*w()); 00211 switch (discr) { 00212 case CGAL::NEGATIVE: 00213 empty = true; 00214 o = (CGAL::Orientation)(CGAL_NTS sign (w())); 00215 break; 00216 case CGAL::ZERO: 00217 empty = false; 00218 o = (CGAL::Orientation)(CGAL_NTS sign (r())); 00219 break; 00220 case CGAL::POSITIVE: 00221 empty = false; 00222 o = CGAL::ZERO; 00223 break; 00224 } 00225 } else { 00226 empty = false; 00227 o = (CGAL::Orientation)(-CGAL_NTS sign (r())); 00228 } 00229 } else if (!CGAL_NTS is_zero (s())) { 00230 trivial = false; 00231 degenerate = (t()*v() == RT(2)*s()*u()); 00232 if (degenerate) { 00233 CGAL::Sign discr = (CGAL::Sign) 00234 CGAL_NTS sign(v()*v()-RT(4)*s()*w()); 00235 switch (discr) { 00236 case CGAL::NEGATIVE: 00237 empty = true; 00238 o = (CGAL::Orientation)(CGAL_NTS sign (w())); 00239 break; 00240 case CGAL::ZERO: 00241 empty = false; 00242 o = (CGAL::Orientation)(CGAL_NTS sign (s())); 00243 break; 00244 case CGAL::POSITIVE: 00245 empty = false; 00246 o = CGAL::ZERO; 00247 break; 00248 } 00249 } else { 00250 empty = false; 00251 o = (CGAL::Orientation)(-CGAL_NTS sign (s())); 00252 } 00253 } else { // r=0, s=0 00254 degenerate = true; 00255 bool uv_zero = CGAL_NTS is_zero (u()) 00256 && CGAL_NTS is_zero (v()); 00257 trivial = uv_zero && CGAL_NTS is_zero (w()); 00258 empty = uv_zero && !trivial; 00259 if (empty) 00260 o = (CGAL::Orientation)(CGAL_NTS sign (w())); 00261 else if (trivial) 00262 o = CGAL::POSITIVE; 00263 else 00264 o = CGAL::ZERO; 00265 } 00266 00267 00268 } 00269 break; 00270 case ELLIPSE: 00271 { 00272 trivial = false; 00273 RT z_prime = d*w() - u()*u()*s() - v()*v()*r() + u()*v()*t(); 00274 if (CGAL_NTS is_positive (r())) { 00275 empty = (CGAL::POSITIVE == CGAL_NTS sign (z_prime)); 00276 empty ? o = CGAL::POSITIVE : o = CGAL::NEGATIVE; 00277 } else { 00278 empty = (CGAL::NEGATIVE == CGAL_NTS sign (z_prime)); 00279 empty ? o = CGAL::NEGATIVE : o = CGAL::POSITIVE; 00280 } 00281 degenerate = empty || CGAL_NTS is_zero (z_prime); 00282 } 00283 break; 00284 } 00285 } 00286 00287 RT evaluate (const PT& p) const 00288 { 00289 RT x, y, h; 00290 dao.get (p, x, y, h); 00291 return r()*x*x + s()*y*y + t()*x*y + u()*x*h + v()*y*h + w()*h*h; 00292 } 00293 00294 00295 00296 public: 00297 ConicHPA2 ( const DA& da = DA()) : dao( da) { } 00298 00299 ConicHPA2 (RT r, RT s, RT t, RT u, RT v, RT w, const DA& da = DA()) 00300 : dao( da), _r(r), _s(s), _t(t), _u(u), _v(v), _w(w) 00301 { 00302 analyse(); 00303 } 00304 00305 const DA& da() const 00306 { 00307 return dao; 00308 } 00309 00310 RT r() const { return _r;} 00311 RT s() const { return _s;} 00312 RT t() const { return _t;} 00313 RT u() const { return _u;} 00314 RT v() const { return _v;} 00315 RT w() const { return _w;} 00316 00317 PT center () const 00318 { 00319 CGAL_kernel_precondition (type != PARABOLA); 00320 // PT p; 00321 // replaced previous line by following hack (no idea 00322 // why original version doesn't work) 00323 typename DA::Point_2 p; 00324 RT two = RT(2); 00325 dao.set( p, two*s()*u() - t()*v(), two*r()*v() - t()*u(), -det()); 00326 return p; 00327 } 00328 00329 Conic_type conic_type () const 00330 { 00331 return type; 00332 } 00333 00334 bool is_hyperbola () const 00335 { 00336 return (type == HYPERBOLA); 00337 } 00338 00339 bool is_parabola () const 00340 { 00341 return (type == PARABOLA); 00342 } 00343 00344 bool is_ellipse () const 00345 { 00346 return (type == ELLIPSE); 00347 } 00348 00349 bool is_circle () const 00350 { 00351 return (type == ELLIPSE && (r()==s()) && CGAL_NTS is_zero (t())); 00352 } 00353 00354 bool is_empty () const 00355 { 00356 return empty; 00357 } 00358 00359 bool is_trivial () const 00360 { 00361 return trivial; 00362 } 00363 00364 bool is_degenerate () const 00365 { 00366 return degenerate; 00367 } 00368 00369 CGAL::Orientation orientation () const 00370 { 00371 return o; 00372 } 00373 00374 CGAL::Oriented_side oriented_side (const PT& p) const 00375 { 00376 return (CGAL::Oriented_side)(CGAL_NTS sign (evaluate (p))); 00377 } 00378 00379 bool has_on_positive_side (const PT& p) const 00380 { 00381 return (CGAL_NTS is_positive (evaluate(p))); 00382 } 00383 00384 bool has_on_negative_side (const PT& p) const 00385 { 00386 return (CGAL_NTS is_negative (evaluate(p))); 00387 } 00388 00389 bool has_on_boundary (const PT& p) const 00390 { 00391 return (CGAL_NTS is_zero (evaluate(p))); 00392 } 00393 00394 bool has_on (const PT& p) const 00395 { 00396 return (CGAL_NTS is_zero (evaluate(p))); 00397 } 00398 00399 Convex_side convex_side (const PT& p) const 00400 { 00401 switch (o) { 00402 case CGAL::POSITIVE: 00403 return (Convex_side)( CGAL_NTS sign (evaluate (p))); 00404 case CGAL::NEGATIVE: 00405 return (Convex_side)(-CGAL_NTS sign (evaluate (p))); 00406 case CGAL::ZERO: 00407 return (Convex_side)( CGAL_NTS sign (CGAL_NTS abs (evaluate(p)))); 00408 } 00409 // keeps g++ happy 00410 return( Convex_side( 0)); 00411 } 00412 00413 bool has_on_convex_side (const PT& p) const 00414 { 00415 return (convex_side (p) == ON_CONVEX_SIDE); 00416 } 00417 00418 bool has_on_nonconvex_side (const PT& p) const 00419 { 00420 return (convex_side (p) == ON_NONCONVEX_SIDE); 00421 } 00422 00423 bool operator == ( const ConicHPA2<_PT,_DA>& c) const 00424 { 00425 // find coefficient != 0 00426 RT factor1(0); 00427 if ( ! CGAL_NTS is_zero( r())) factor1 = r(); else 00428 if ( ! CGAL_NTS is_zero( s())) factor1 = s(); else 00429 if ( ! CGAL_NTS is_zero( t())) factor1 = t(); else 00430 if ( ! CGAL_NTS is_zero( u())) factor1 = u(); else 00431 if ( ! CGAL_NTS is_zero( v())) factor1 = v(); else 00432 if ( ! CGAL_NTS is_zero( w())) factor1 = w(); else 00433 CGAL_kernel_assertion_msg( false, "all coefficients zero"); 00434 00435 // find coefficient != 0 00436 RT factor2(0); 00437 if ( ! CGAL_NTS is_zero( c.r())) factor2 = c.r(); else 00438 if ( ! CGAL_NTS is_zero( c.s())) factor2 = c.s(); else 00439 if ( ! CGAL_NTS is_zero( c.t())) factor2 = c.t(); else 00440 if ( ! CGAL_NTS is_zero( c.u())) factor2 = c.u(); else 00441 if ( ! CGAL_NTS is_zero( c.v())) factor2 = c.v(); else 00442 if ( ! CGAL_NTS is_zero( c.w())) factor2 = c.w(); else 00443 CGAL_kernel_assertion_msg( false, "all coefficients zero"); 00444 00445 return( ( r()*factor2 == c.r()*factor1) 00446 && ( s()*factor2 == c.s()*factor1) 00447 && ( t()*factor2 == c.t()*factor1) 00448 && ( u()*factor2 == c.u()*factor1) 00449 && ( v()*factor2 == c.v()*factor1) 00450 && ( w()*factor2 == c.w()*factor1)); 00451 } 00452 00453 void set (RT r_, RT s_, RT t_, RT u_, RT v_, RT w_) 00454 { 00455 _r = r_; _s = s_; _t = t_; _u = u_; _v = v_; _w = w_; 00456 analyse(); 00457 } 00458 00459 void set_opposite () 00460 { 00461 _r = -r(); _s = -s(); _t = -t(); _u = -u(); _v = -v(); _w = -w(); 00462 o = CGAL::opposite(orientation()); 00463 } 00464 00465 void set_circle (const PT& p1, const PT& p2, const PT& p3) 00466 { 00467 // the circle will have r = s = det*h1*h2*h3, t=0 00468 RT x1, y1, h1, x2, y2, h2, x3, y3, h3; 00469 dao.get (p1, x1, y1, h1); 00470 dao.get (p2, x2, y2, h2); 00471 dao.get (p3, x3, y3, h3); 00472 00473 // precondition: p1, p2, p3 not collinear 00474 RT det = -h1*x3*y2+h3*x1*y2+h1*x2*y3-h2*x1*y3+h2*x3*y1-h3*x2*y1; 00475 CGAL_kernel_precondition (!CGAL_NTS is_zero (det)); 00476 00477 // Cramer's rule 00478 RT sqr1 = (-x1*x1 - y1*y1)*h2*h3; 00479 RT sqr2 = (-x2*x2 - y2*y2)*h1*h3; 00480 RT sqr3 = (-x3*x3 - y3*y3)*h1*h2; 00481 00482 _u = -h1*sqr3*y2+h3*sqr1*y2+h1*sqr2*y3-h2*sqr1*y3+h2*sqr3*y1-h3*sqr2*y1; 00483 _v = -h1*x3*sqr2+h3*x1*sqr2+h1*x2*sqr3-h2*x1*sqr3+h2*x3*sqr1-h3*x2*sqr1; 00484 _w = -sqr1*x3*y2+sqr3*x1*y2+sqr1*x2*y3-sqr2*x1*y3+sqr2*x3*y1-sqr3*x2*y1; 00485 _r = det*h1*h2*h3; 00486 _s = _r; 00487 _t = RT(0); 00488 00489 analyse(); 00490 CGAL_kernel_postcondition(is_circle()); 00491 CGAL_kernel_postcondition(has_on_boundary(p1)); 00492 CGAL_kernel_postcondition(has_on_boundary(p2)); 00493 CGAL_kernel_postcondition(has_on_boundary(p3)); 00494 } 00495 00496 void set_linepair (const PT& p1, const PT& p2, const PT& p3, 00497 const PT& p4, const DA& da = DA()) 00498 { 00499 RT x1, y1, h1, x2, y2, h2, x3, y3, h3, x4, y4, h4; 00500 da.get (p1, x1, y1, h1); 00501 da.get (p2, x2, y2, h2); 00502 da.get (p3, x3, y3, h3); 00503 da.get (p4, x4, y4, h4); 00504 00505 // precondition: p1 != p2, p3 != p4 00506 CGAL_kernel_precondition 00507 ( ((x1*h2 != x2*h1) || (y1*h2 != y2*h1)) && 00508 ((x3*h4 != x4*h3) || (y3*h4 != y4*h3)) ); 00509 00510 RT h1x2_x1h2 = h1*x2-x1*h2; 00511 RT h3x4_x3h4 = h3*x4-x3*h4; 00512 RT y1h2_h1y2 = y1*h2-h1*y2; 00513 RT y3h4_h3y4 = y3*h4-h3*y4; 00514 RT x1y2_y1x2 = x1*y2-y1*x2; 00515 RT x3y4_y3x4 = x3*y4-y3*x4; 00516 00517 _r = y1h2_h1y2 * y3h4_h3y4; 00518 _s = h1x2_x1h2 * h3x4_x3h4; 00519 _t = h1x2_x1h2 * y3h4_h3y4 + y1h2_h1y2 * h3x4_x3h4; 00520 _u = x1y2_y1x2 * y3h4_h3y4 + y1h2_h1y2 * x3y4_y3x4; 00521 _v = x1y2_y1x2 * h3x4_x3h4 + h1x2_x1h2 * x3y4_y3x4; 00522 _w = x1y2_y1x2 * x3y4_y3x4; 00523 00524 analyse(); 00525 } 00526 00527 void set_ellipse (const PT& p1, const PT& p2, const PT& p3) 00528 { 00529 RT x1, y1, h1, x2, y2, h2, x3, y3, h3; 00530 dao.get (p1, x1, y1, h1); 00531 dao.get (p2, x2, y2, h2); 00532 dao.get (p3, x3, y3, h3); 00533 00534 // precondition: p1, p2, p3 not collinear 00535 RT det = -h1*x3*y2+h3*x1*y2+h1*x2*y3-h2*x1*y3+h2*x3*y1-h3*x2*y1; 00536 CGAL_kernel_precondition (!CGAL_NTS is_zero (det)); 00537 00538 RT x1x1 = x1*x1, y1y1 = y1*y1, 00539 x2x2 = x2*x2, y2y2 = y2*y2, 00540 x3x3 = x3*x3, y3y3 = y3*y3, // x_i^2, y_i^2 00541 x1h1 = x1*h1, y1h1 = y1*h1, 00542 x2h2 = x2*h2, y2h2 = y2*h2, 00543 x3h3 = x3*h3, y3h3 = y3*h3, // x_i h_i, y_i h_i 00544 h1h1 = h1*h1, 00545 h2h2 = h2*h2, 00546 h3h3 = h3*h3, // h_i^2 00547 two = RT(2); // 2 00548 00549 _r = y1y1*h2h2*h3h3 - y1h1*y2h2*h3h3 - y1h1*h2h2*y3h3 + 00550 h1h1*y2y2*h3h3 - h1h1*y2h2*y3h3 + h1h1*h2h2*y3y3; 00551 00552 _s = x1x1*h2h2*h3h3 - x1h1*x2h2*h3h3 - x1h1*h2h2*x3h3 + 00553 h1h1*x2x2*h3h3 - h1h1*x2h2*x3h3 + h1h1*h2h2*x3x3; 00554 00555 _t = -two*x1*y1*h2h2*h3h3 + x1h1*y2h2*h3h3 + x1h1*h2h2*y3h3 + 00556 y1h1*x2h2*h3h3 -two*h1h1*x2*y2*h3h3 + h1h1*x2h2*y3h3 + 00557 y1h1*h2h2*x3h3 + h1h1*y2h2*x3h3 -two*h1h1*h2h2*x3*y3; 00558 00559 _u = -(h1h1*y2y2*x3h3 - h1h1*x2*y2*y3h3 - h1h1*y2h2*x3*y3 + 00560 x1h1*h2h2*y3y3 + h1h1*x2h2*y3y3 +y1y1*x2h2*h3h3 + 00561 y1y1*h2h2*x3h3 - x1*y1*y2h2*h3h3 - y1h1*x2*y2*h3h3 - 00562 x1*y1*h2h2*y3h3 - y1h1*h2h2*x3*y3 + x1h1*y2y2*h3h3); 00563 00564 _v = -(h1h1*x2x2*y3h3 - h1h1*x2*y2*x3h3 - h1h1*x2h2*x3*y3 + 00565 y1h1*h2h2*x3x3 + h1h1*y2h2*x3x3 +x1x1*y2h2*h3h3 + 00566 x1x1*h2h2*y3h3 - x1*y1*x2h2*h3h3 - x1h1*x2*y2*h3h3 - 00567 x1*y1*h2h2*x3h3 - x1h1*h2h2*x3*y3 + y1h1*x2x2*h3h3); 00568 00569 _w = y1y1*x2h2*x3h3 - x1*y1*y2h2*x3h3 - y1h1*x2*y2*x3h3 + 00570 y1h1*y2h2*x3x3 - x1*y1*x2h2*y3h3 + y1h1*x2x2*y3h3 - 00571 y1h1*x2h2*x3*y3 + x1h1*y2y2*x3h3 + x1x1*y2h2*y3h3 - 00572 x1h1*x2*y2*y3h3 - x1h1*y2h2*x3*y3 + x1h1*x2h2*y3y3; 00573 00574 type = ELLIPSE; 00575 degenerate = trivial = empty = false; 00576 o = CGAL::NEGATIVE; 00577 if (CGAL_NTS is_positive (det)) set_opposite (); 00578 00579 } 00580 00581 void set_ellipse (const PT& p1, const PT& p2, 00582 const PT& p3, const PT& p4, 00583 CGAL::Orientation _o = POSITIVE) 00584 { 00585 ConicHPA2<PT,DA> pair1, pair2; 00586 set_two_linepairs (p1, p2, p3, p4, pair1, pair2); 00587 set_ellipse (pair1, pair2); 00588 analyse(); 00589 if (o != _o) set_opposite(); 00590 } 00591 00592 void set (const PT& p1, const PT& p2, const PT& p3, const PT& p4, 00593 const PT& p5, CGAL::Orientation _o = POSITIVE) 00594 { 00595 ConicHPA2<PT,DA> c1; c1.set_linepair (p1, p2, p3, p4); 00596 ConicHPA2<PT,DA> c2; c2.set_linepair (p1, p4, p2, p3); 00597 set_linear_combination (c2.evaluate (p5), c1, 00598 -c1.evaluate (p5), c2); 00599 analyse(); 00600 // precondition: all points distinct <=> conic nontrivial 00601 CGAL_kernel_precondition (!is_trivial()); 00602 if (o != _o) set_opposite(); 00603 } 00604 00605 00606 }; 00607 00608 #ifndef CGAL_NO_OSTREAM_INSERT_CONICHPA2 00609 template< class _PT, class _DA> 00610 std::ostream& operator << ( std::ostream& os, const ConicHPA2<_PT,_DA>& c) 00611 { 00612 return( os << c.r() << ' ' << c.s() << ' ' << c.t() << ' ' 00613 << c.u() << ' ' << c.v() << ' ' << c.w()); 00614 } 00615 00616 template< class _PT, class _DA> 00617 std::istream& operator >> ( std::istream& is, ConicHPA2<_PT,_DA>& c) 00618 { 00619 typedef ConicHPA2<_PT,_DA> Conic; 00620 typedef typename _DA::RT RT; 00621 00622 RT r, s, t, u, v, w; 00623 is >> r >> s >> t >> u >> v >> w; 00624 c.set( r, s, t, u, v, w); 00625 00626 return( is); 00627 } 00628 #endif // CGAL_NO_OSTREAM_INSERT_CONICHPA2 00629 00630 CGAL_END_NAMESPACE 00631 00632 #endif // CGAL_CONICHPA2_H 00633 00634 // ===== EOF ==================================================================
1.7.6.1
