1 // Copyright (c) 2005-2006 INRIA Sophia-Antipolis (France). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // Partially supported by the IST Programme of the EU as a Shared-cost 7 // RTD (FET Open) Project under Contract No IST-2000-26473 8 // (ECG - Effective Computational Geometry for Curves and Surfaces) 9 // and a STREP (FET Open) Project under Contract No IST-006413 10 // (ACS -- Algorithms for Complex Shapes) 11 // 12 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Algebraic_kernel_for_spheres/include/CGAL/Root_for_spheres_2_3.h $ 13 // $Id: Root_for_spheres_2_3.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 14 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 15 // 16 // Author(s) : Monique Teillaud <Monique.Teillaud@sophia.inria.fr> 17 // Sylvain Pion 18 // Pedro Machado 19 // Julien Hazebrouck 20 // Damien Leroy 21 22 #ifndef CGAL_ROOT_FOR_SPHERES_2_3_H 23 #define CGAL_ROOT_FOR_SPHERES_2_3_H 24 25 #include <CGAL/license/Circular_kernel_3.h> 26 27 28 #include <iostream> 29 #include <CGAL/Polynomials_1_3.h> 30 #include <CGAL/Polynomials_2_3.h> 31 #include <CGAL/Polynomials_for_line_3.h> 32 #include <CGAL/Bbox_3.h> 33 #include <CGAL/Root_of_traits.h> 34 35 namespace CGAL { 36 37 template < typename RT_ > 38 class Root_for_spheres_2_3 { 39 40 typedef RT_ RT; 41 typedef typename Root_of_traits< RT >::RootOf_2 Root_of_2; 42 typedef typename Root_of_traits< RT >::RootOf_1 FT; 43 typedef CGAL::Polynomial_1_3< FT > Polynomial_1_3; 44 typedef CGAL::Polynomial_for_spheres_2_3< FT > Polynomial_for_spheres_2_3; 45 typedef CGAL::Polynomials_for_line_3< FT > Polynomials_for_line_3; 46 47 private: 48 Root_of_2 x_; 49 Root_of_2 y_; 50 Root_of_2 z_; 51 52 public: Root_for_spheres_2_3()53 Root_for_spheres_2_3(){} 54 55 Root_for_spheres_2_3(const Root_of_2 & r1,const Root_of_2 & r2,const Root_of_2 & r3)56 Root_for_spheres_2_3(const Root_of_2& r1, 57 const Root_of_2& r2, 58 const Root_of_2& r3) 59 : x_(r1), y_(r2), z_(r3) 60 { 61 // This assertion sont work if Root_of_2 is 62 // Interval_nt (and dont have is_rational, gamma, etc..) 63 /*CGAL_assertion( 64 ((r1.is_rational() && r2.is_rational()) || 65 (r1.is_rational() && r3.is_rational()) || 66 (r2.is_rational() && r3.is_rational()) || 67 ((r1.is_rational()) && (r2.gamma() == r3.gamma())) || 68 ((r2.is_rational()) && (r1.gamma() == r3.gamma())) || 69 ((r3.is_rational()) && (r1.gamma() == r2.gamma())) || 70 ((r1.gamma() == r2.gamma()) && (r2.gamma() == r3.gamma()))) 71 );*/ 72 } 73 x()74 const Root_of_2& x() const 75 { return x_; } 76 y()77 const Root_of_2& y() const 78 { return y_; } 79 z()80 const Root_of_2& z() const 81 { return z_; } 82 83 // On fait l'evaluation de (x,y,z) pour le plan 84 // aX + bY + cZ + d, donne evaluate(const Polynomial_1_3 & p)85 const Root_of_2 evaluate(const Polynomial_1_3 &p) const { 86 return (p.a() * x()) + (p.b() * y()) + (p.c() * z()) + p.d(); 87 } 88 89 // On fait l'evaluation de (x,y,z) pour le plan 90 // (X-a)^2 + (Y-b)^2 + (Z-c)^2 - r_sq, donne evaluate(const Polynomial_for_spheres_2_3 & p)91 const Root_of_2 evaluate(const Polynomial_for_spheres_2_3 &p) const { 92 return square(x() - p.a()) + 93 square(y() - p.b()) + 94 square(z() - p.c()) - 95 p.r_sq(); 96 } 97 98 // On verifie si (x,y,z) fait partie la ligne donne is_on_line(const Polynomials_for_line_3 & p)99 bool is_on_line(const Polynomials_for_line_3 &p) const { 100 Root_of_2 t; 101 bool already = false; 102 if(!is_zero(p.a1())) { 103 t = (x() - p.b1())/p.a1(); 104 already = true; 105 } else if(p.b1() != x()) return false; 106 if(!is_zero(p.a2())) { 107 if(!already) { 108 t = (y() - p.b2())/p.a2(); 109 already = true; 110 } 111 else if((p.a2() * t + p.b2()) != y()) return false; 112 } else if(p.b2() != y()) return false; 113 if(!is_zero(p.a3())) { 114 if(!already) return true; 115 else if((p.a3() * t + p.b3()) != z()) return false; 116 } else if(p.b3() != z()) return false; 117 return true; 118 } 119 bbox()120 CGAL::Bbox_3 bbox() const 121 { 122 const Root_of_2 &ox = x(); 123 const Root_of_2 &oy = y(); 124 const Root_of_2 &oz = z(); 125 126 CGAL::Interval_nt<> 127 ix=to_interval(ox), 128 iy=to_interval(oy), 129 iz=to_interval(oz); 130 return CGAL::Bbox_3(ix.inf(),iy.inf(),iz.inf(), 131 ix.sup(),iy.sup(),iz.sup()); 132 /* 133 // Note: This is a more efficient version 134 // but it won't work (in the future) 135 // with some Lazy_Curved_kernel_3 136 // because is_rational(), gamma(), etc.. is not defined 137 // for Interval_nt<false> data type 138 const Root_of_2 &ox = x(); 139 const Root_of_2 &oy = y(); 140 const Root_of_2 &oz = z(); 141 142 const bool x_rat = ox.is_rational(); 143 const bool y_rat = oy.is_rational(); 144 const bool z_rat = oz.is_rational(); 145 146 if(((x_rat?1:0) + (y_rat?1:0) +(z_rat?1:0)) > 1) { 147 CGAL::Interval_nt<> 148 ix=to_interval(ox), 149 iy=to_interval(oy), 150 iz=to_interval(oz); 151 return CGAL::Bbox_3(ix.inf(),iy.inf(),iz.inf(), 152 ix.sup(),iy.sup(),iz.sup()); 153 } 154 155 if(z_rat) { 156 const CGAL::Interval_nt<true> alpha1 = to_interval(ox.alpha()); 157 const CGAL::Interval_nt<true> beta1 = to_interval(ox.beta()); 158 const CGAL::Interval_nt<true> alpha2 = to_interval(oy.alpha()); 159 const CGAL::Interval_nt<true> beta2 = to_interval(oy.beta()); 160 const CGAL::Interval_nt<true> g = to_interval(ox.gamma()); 161 const CGAL::Interval_nt<true> sqrtg = CGAL::sqrt(g); 162 const CGAL::Interval_nt<true> ix = alpha1 + beta1 * sqrtg; 163 const CGAL::Interval_nt<true> iy = alpha2 + beta2 * sqrtg; 164 const CGAL::Interval_nt<true> iz = to_interval(oz); 165 return CGAL::Bbox_3(ix.inf(),iy.inf(),iz.inf(), 166 ix.sup(),iy.sup(),iz.sup()); 167 } 168 169 if(y_rat) { 170 const CGAL::Interval_nt<true> alpha1 = to_interval(ox.alpha()); 171 const CGAL::Interval_nt<true> beta1 = to_interval(ox.beta()); 172 const CGAL::Interval_nt<true> alpha2 = to_interval(oz.alpha()); 173 const CGAL::Interval_nt<true> beta2 = to_interval(oz.beta()); 174 const CGAL::Interval_nt<true> g = to_interval(ox.gamma()); 175 const CGAL::Interval_nt<true> sqrtg = CGAL::sqrt(g); 176 const CGAL::Interval_nt<true> ix = alpha1 + beta1 * sqrtg; 177 const CGAL::Interval_nt<true> iz = alpha2 + beta2 * sqrtg; 178 const CGAL::Interval_nt<true> iy = to_interval(oy); 179 return CGAL::Bbox_3(ix.inf(),iy.inf(),iz.inf(), 180 ix.sup(),iy.sup(),iz.sup()); 181 } 182 183 if(x_rat) { 184 const CGAL::Interval_nt<true> alpha1 = to_interval(oy.alpha()); 185 const CGAL::Interval_nt<true> beta1 = to_interval(oy.beta()); 186 const CGAL::Interval_nt<true> alpha2 = to_interval(oz.alpha()); 187 const CGAL::Interval_nt<true> beta2 = to_interval(oz.beta()); 188 const CGAL::Interval_nt<true> g = to_interval(oy.gamma()); 189 const CGAL::Interval_nt<true> sqrtg = CGAL::sqrt(g); 190 const CGAL::Interval_nt<true> iy = alpha1 + beta1 * sqrtg; 191 const CGAL::Interval_nt<true> iz = alpha2 + beta2 * sqrtg; 192 const CGAL::Interval_nt<true> ix = to_interval(ox); 193 return CGAL::Bbox_3(ix.inf(),iy.inf(),iz.inf(), 194 ix.sup(),iy.sup(),iz.sup()); 195 } 196 197 const CGAL::Interval_nt<true> alpha1 = to_interval(ox.alpha()); 198 const CGAL::Interval_nt<true> beta1 = to_interval(ox.beta()); 199 const CGAL::Interval_nt<true> alpha2 = to_interval(oy.alpha()); 200 const CGAL::Interval_nt<true> beta2 = to_interval(oy.beta()); 201 const CGAL::Interval_nt<true> alpha3 = to_interval(oz.alpha()); 202 const CGAL::Interval_nt<true> beta3 = to_interval(oz.beta()); 203 const CGAL::Interval_nt<true> g = to_interval(ox.gamma()); 204 const CGAL::Interval_nt<true> sqrtg = CGAL::sqrt(g); 205 const CGAL::Interval_nt<true> ix = alpha1 + beta1 * sqrtg; 206 const CGAL::Interval_nt<true> iy = alpha2 + beta2 * sqrtg; 207 const CGAL::Interval_nt<true> iz = alpha3 + beta3 * sqrtg; 208 return CGAL::Bbox_3(ix.inf(),iy.inf(),iz.inf(), 209 ix.sup(),iy.sup(),iz.sup()); 210 */ 211 } 212 213 }; 214 215 template < typename RT > 216 bool 217 operator == ( const Root_for_spheres_2_3<RT>& r1, 218 const Root_for_spheres_2_3<RT>& r2 ) 219 { return (r1.x() == r2.x()) && (r1.y() == r2.y()) && (r1.z() == r2.z()); } 220 221 template < typename RT > 222 std::ostream & 223 operator<<(std::ostream & os, const Root_for_spheres_2_3<RT> &r) 224 { return os << r.x() << " " << r.y() << " " << r.z() << " "; } 225 226 template < typename RT > 227 std::istream & 228 operator>>(std::istream & is, Root_for_spheres_2_3<RT> &r) 229 { 230 typedef typename Root_of_traits< RT >::RootOf_2 Root_of_2; 231 Root_of_2 x,y,z; 232 233 is >> x >> y >> z; 234 if(is) 235 r = Root_for_spheres_2_3<RT>(x,y,z); 236 return is; 237 } 238 239 } //namespace CGAL 240 241 #endif // CGAL_ROOT_FOR_SPHERES_2_3_H 242