|
BWAPI
|
00001 // Copyright (c) 2005 Tel-Aviv University (Israel). 00002 // 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/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_intersections_2.h $ 00015 // $Id: Conic_intersections_2.h 33081 2006-08-07 07:46:37Z wein $ 00016 // 00017 // 00018 // Author(s) : Ron Wein <wein@post.tau.ac.il> 00019 00020 #ifndef CGAL_CONIC_INTERSECTIONS_2_H 00021 #define CGAL_CONIC_INTERSECTIONS_2_H 00022 00027 CGAL_BEGIN_NAMESPACE 00028 00040 template <class Nt_traits> 00041 int _compute_resultant_roots (Nt_traits& nt_traits, 00042 const typename Nt_traits::Integer& r1, 00043 const typename Nt_traits::Integer& s1, 00044 const typename Nt_traits::Integer& t1, 00045 const typename Nt_traits::Integer& u1, 00046 const typename Nt_traits::Integer& v1, 00047 const typename Nt_traits::Integer& w1, 00048 const int& deg1, 00049 const typename Nt_traits::Integer& r2, 00050 const typename Nt_traits::Integer& s2, 00051 const typename Nt_traits::Integer& t2, 00052 const typename Nt_traits::Integer& u2, 00053 const typename Nt_traits::Integer& v2, 00054 const typename Nt_traits::Integer& w2, 00055 const int& deg2, 00056 typename Nt_traits::Algebraic *xs) 00057 { 00058 if (deg1 == 2 && deg2 == 1) 00059 { 00060 // If necessary, swap roles between the two curves, so that the first 00061 // curve always has the minimal degree. 00062 return (_compute_resultant_roots (nt_traits, 00063 r2, s2, t2, u2, v2, w2, 00064 deg2, 00065 r1, s1, t1, u1, v1, w1, 00066 deg1, 00067 xs)); 00068 } 00069 00070 // Act according to the degree of the first conic curve. 00071 const typename Nt_traits::Integer _two = 2; 00072 typename Nt_traits::Integer c[5]; 00073 unsigned int degree = 4; 00074 typename Nt_traits::Algebraic *xs_end; 00075 00076 if (deg1 == 1) 00077 { 00078 // The first curve has no quadratic coefficients, and represents a line. 00079 if (CGAL::sign (v1) == ZERO) 00080 { 00081 // The first line is u1*x + w1 = 0, therefore: 00082 xs[0] = nt_traits.convert(-w1) / nt_traits.convert(u1); 00083 return (1); 00084 } 00085 00086 // We can write the first curve as: y = -(u1*x + w1) / v1. 00087 if (deg2 == 1) 00088 { 00089 // The second curve is also a line. We therefore get the linear 00090 // equation c[1]*x + c[0] = 0: 00091 c[1] = v1*u2 - u1*v2; 00092 c[0] = v1*w2 - w1*v2; 00093 00094 if (CGAL::sign (c[1]) == ZERO) 00095 // The two lines are parallel: 00096 return (0); 00097 00098 xs[0] = nt_traits.convert(-c[0]) / nt_traits.convert(c[1]); 00099 return (1); 00100 } 00101 00102 // We substitute this expression into the equation of the second 00103 // conic, and get the quadratic equation c[2]*x^2 + c[1]*x + c[0] = 0: 00104 c[2] = u1*u1*s2 - u1*v1*t2 + v1*v1*r2; 00105 c[1] = _two*u1*w1*s2 - u1*v1*v2 - v1*w1*t2 + v1*v1*u2; 00106 c[0] = w1*w1*s2 - v1*w1*v2 + v1*v1*w2; 00107 00108 xs_end = nt_traits.solve_quadratic_equation (c[2], c[1], c[0], 00109 xs); 00110 return (xs_end - xs); 00111 } 00112 00113 // At this stage, both curves have degree 2. We obtain a qaurtic polynomial 00114 // whose roots are the x-coordinates of the intersection points. 00115 if (CGAL::sign (s1) == ZERO && CGAL::sign (s2) == ZERO) 00116 { 00117 // If both s1 and s2 are zero, we can write the two curves as: 00118 // C1: (t1*x + v1)*y + (r1*x^2 + u1*x + w1) = 0 00119 // C2: (t2*x + v2)*y + (r2*x^2 + u2*x + w2) = 0 00120 // By writing the resultant of these two polynomials we get a cubic 00121 // polynomial, whose coefficients are given by: 00122 c[3] = r2*t1 - r1*t2; 00123 c[2] = t1*u2 - t2*u1 + r2*v1 - r1*v2; 00124 c[1] = t1*w2 - t2*w1 + u2*v1 - u1*v2; 00125 c[0] = v1*w2 - v2*w1; 00126 00127 degree = 3; 00128 } 00129 else 00130 { 00131 // We can write the two curves as: 00132 // C1: (s1)*y^2 + (t1*x + v1)*y + (r1*x^2 + u1*x + w1) = 0 00133 // C2: (s2)*y^2 + (t2*x + v2)*y + (r2*x^2 + u2*x + w2) = 0 00134 // By writing the resultant of these two polynomials we get a quartic 00135 // polynomial, whose coefficients are given by: 00136 c[4] = -_two*s1*s2*r1*r2 + s1*t2*t2*r1 - s1*t2*t1*r2 + 00137 s1*s1*r2*r2 - s2*t1*r1*t2 + s2*t1*t1*r2 + s2*s2*r1*r1; 00138 00139 c[3] = -t2*r1*v1*s2 - u2*t1*t2*s1 - v2*r1*t1*s2 - 00140 r2*t1*v2*s1 - _two*s1*s2*r1*u2 - t2*u1*t1*s2 + u2*t1*t1*s2 - 00141 r2*v1*t2*s1 + u1*t2*t2*s1 + _two*v2*r1*t2*s1 + _two*u2*r2*s1*s1 + 00142 _two*r2*v1*t1*s2 + _two*u1*r1*s2*s2 - _two*s1*s2*u1*r2; 00143 00144 c[2] = -r2*v1*v2*s1 + u2*u2*s1*s1 + _two*w2*r2*s1*s1 + 00145 _two*u2*v1*t1*s2 - u2*v1*t2*s1 + w2*t1*t1*s2 - _two*s1*s2*u1*u2 - 00146 w2*t1*t2*s1 + v2*v2*r1*s1 + u1*u1*s2*s2 - v2*r1*v1*s2 + 00147 _two*w1*r1*s2*s2 - u2*t1*v2*s1 - t2*u1*v1*s2 - _two*s1*s2*r1*w2 - 00148 _two*s1*s2*w1*r2 + r2*v1*v1*s2 + w1*t2*t2*s1 - v2*u1*t1*s2 - 00149 t2*w1*t1*s2 + _two*v2*u1*t2*s1; 00150 00151 c[1] = _two*w2*u2*s1*s1 + _two*w2*v1*t1*s2 - w2*v1*t2*s1 + 00152 _two*v2*w1*t2*s1 + _two*w1*u1*s2*s2 - v2*u1*v1*s2 - _two*s1*s2*u1*w2 - 00153 v2*w1*t1*s2 + u2*v1*v1*s2 - t2*w1*v1*s2 - w2*t1*v2*s1 + 00154 v2*v2*u1*s1 - u2*v1*v2*s1 - _two*s1*s2*w1*u2; 00155 00156 c[0] = s2*v1*v1*w2 - s1*v2*v1*w2 - s2*v1*w1*v2 + s2*s2*w1*w1 - 00157 _two*s1*s2*w1*w2 + s1*w1*v2*v2 + s1*s1*w2*w2; 00158 00159 degree = 4; 00160 } 00161 00162 // Compute the roots of the resultant polynomial. 00163 typename Nt_traits::Polynomial poly = 00164 nt_traits.construct_polynomial (c, degree); 00165 00166 xs_end = nt_traits.compute_polynomial_roots (poly, 00167 xs); 00168 return (xs_end - xs); 00169 } 00170 00181 template <class Nt_traits> 00182 int _compute_resultant_roots (Nt_traits& nt_traits, 00183 const typename Nt_traits::Algebraic& r, 00184 const typename Nt_traits::Algebraic& s, 00185 const typename Nt_traits::Algebraic& t, 00186 const typename Nt_traits::Algebraic& u, 00187 const typename Nt_traits::Algebraic& v, 00188 const typename Nt_traits::Algebraic& w, 00189 const int& deg1, 00190 const typename Nt_traits::Algebraic& A, 00191 const typename Nt_traits::Algebraic& B, 00192 const typename Nt_traits::Algebraic& C, 00193 typename Nt_traits::Algebraic *xs) 00194 { 00195 if (deg1 == 1) 00196 { 00197 // We should actually compute the intersection of two line: 00198 // (u*x + v*y + w = 0) and (A*x + B*y + C = 0): 00199 const typename Nt_traits::Algebraic denom = A*v - B*u; 00200 00201 if (CGAL::sign (denom) == CGAL::ZERO) 00202 // The two lines are parallel and do not intersect. 00203 return (0); 00204 00205 xs[0] = (B*w - C*v) / denom; 00206 return (1); 00207 } 00208 00209 if (CGAL::sign (B) == CGAL::ZERO) 00210 { 00211 // The first line is A*x + C = 0, therefore: 00212 xs[0] = -C / A; 00213 return (1); 00214 } 00215 00216 // We can write the first curve as: y = -(A*x + C) / B. 00217 // We substitute this expression into the equation of the conic, and get 00218 // the quadratic equation c[2]*x^2 + c[1]*x + c[0] = 0: 00219 const typename Nt_traits::Algebraic _two = 2; 00220 typename Nt_traits::Algebraic c[3]; 00221 typename Nt_traits::Algebraic *xs_end; 00222 00223 c[2] = A*A*s - A*B*t + B*B*r; 00224 c[1] = _two*A*C*s - A*B*v - B*C*t + B*B*u; 00225 c[0] = C*C*s - B*C*v + B*B*w; 00226 00227 xs_end = nt_traits.solve_quadratic_equation (c[2], c[1], c[0], 00228 xs); 00229 return (xs_end - xs); 00230 } 00231 00232 CGAL_END_NAMESPACE 00233 00234 #endif
1.7.6.1
