1 // Copyright (c) 2000
2 // Utrecht University (The Netherlands),
3 // ETH Zurich (Switzerland),
4 // INRIA Sophia-Antipolis (France),
5 // Max-Planck-Institute Saarbruecken (Germany),
6 // and Tel-Aviv University (Israel). All rights reserved.
7 //
8 // This file is part of CGAL (www.cgal.org)
9 //
10 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Cartesian_kernel/include/CGAL/Cartesian/Circle_3.h $
11 // $Id: Circle_3.h fe5f655 2021-03-29T15:23:24+02:00 Maxime Gimeno
12 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
13 //
14 // Author(s) : Monique Teillaud, Pedro Machado, Sebastien Loriot
15
16 #ifndef CGAL_CARTESIAN_CIRCLEC3_H
17 #define CGAL_CARTESIAN_CIRCLEC3_H
18
19 #include <CGAL/Interval_nt.h>
20 #include "boost/tuple/tuple.hpp"
21
22 namespace CGAL {
23
24 template <class R_ >
25 class CircleC3 {
26 typedef typename R_::Sphere_3 Sphere_3;
27 typedef typename R_::Plane_3 Plane_3;
28 typedef typename R_::Point_3 Point_3;
29 typedef typename R_::Vector_3 Vector_3;
30 typedef typename R_::Direction_3 Direction_3;
31 typedef typename R_::FT FT;
32
33 //using a boost::tuple because std::pair and tuple cannot work with incomplete types.
34 typedef boost::tuple<Sphere_3, Plane_3> Rep;
35
36
37 typedef typename R_::template Handle<Rep>::type Base;
38 Base base;
39
40 public:
41 typedef R_ R;
42
CircleC3()43 CircleC3() {}
44
CircleC3(const Point_3 & center,const FT & squared_r,const Direction_3 & d)45 CircleC3(const Point_3& center, const FT& squared_r, const Direction_3& d)
46 {
47 CGAL_kernel_assertion(squared_r >= FT(0));
48 // non-degenerated Direction
49 CGAL_kernel_assertion((d.dx() != FT(0)) || (d.dy() != FT(0)) || (d.dz() != FT(0)));
50 base = Rep(Sphere_3(center,squared_r),
51 plane_from_point_direction(center, d));
52 }
53
CircleC3(const Point_3 & center,const FT & squared_r,const Vector_3 & normal)54 CircleC3(const Point_3& center, const FT& squared_r, const Vector_3& normal)
55 {
56 CGAL_kernel_assertion(squared_r >= FT(0));
57 // non-degenerated Vector
58 CGAL_kernel_assertion((normal.x() != FT(0)) ||
59 (normal.y() != FT(0)) ||
60 (normal.z() != FT(0)));
61 base = Rep(Sphere_3(center,squared_r),
62 Plane_3(center, normal.direction()));
63 }
64
CircleC3(const Point_3 & center,const FT & squared_r,const Plane_3 & p)65 CircleC3(const Point_3& center, const FT& squared_r, const Plane_3& p)
66 {
67 // the plane contains the center and it is not degenerate
68 CGAL_kernel_assertion(!R().is_degenerate_3_object()(p));
69 CGAL_kernel_assertion((p.a() * center.x() +
70 p.b() * center.y() +
71 p.c() * center.z() +
72 p.d()) == 0);
73 CGAL_kernel_assertion(squared_r >= FT(0));
74 base = Rep(Sphere_3(center,squared_r), p);
75 }
76
CircleC3(const Sphere_3 & s1,const Sphere_3 & s2)77 CircleC3(const Sphere_3 &s1, const Sphere_3 &s2) {
78 Object obj = R().intersect_3_object()(s1, s2);
79 // s1,s2 must intersect
80 CGAL_kernel_precondition(!(obj.is_empty()));
81 const typename R::Circle_3* circle_ptr=object_cast<typename R::Circle_3>(&obj);
82 if(circle_ptr!=nullptr)
83 base = Rep(circle_ptr->diametral_sphere(), circle_ptr->supporting_plane());
84 else {
85 const typename R::Point_3* point=object_cast<typename R::Point_3>(&obj);
86 CGAL_kernel_precondition(point!=nullptr);
87 CircleC3 circle = CircleC3(*point, FT(0), Vector_3(FT(1),FT(0),FT(0)));
88 base = Rep(circle.diametral_sphere(), circle.supporting_plane());
89 }
90 }
91
CircleC3(const Plane_3 & p,const Sphere_3 & s,int)92 CircleC3(const Plane_3 &p, const Sphere_3 &s, int) : base(s, p) {}
93
CircleC3(const Plane_3 & p,const Sphere_3 & s)94 CircleC3(const Plane_3 &p, const Sphere_3 &s) {
95 Object obj = R().intersect_3_object()(p, s);
96 // s1,s2 must intersect
97 CGAL_kernel_precondition(!(obj.is_empty()));
98 const typename R::Circle_3* circle_ptr=object_cast<typename R::Circle_3>(&obj);
99 if(circle_ptr!=nullptr)
100 base = Rep(circle_ptr->diametral_sphere(), circle_ptr->supporting_plane());
101 else {
102 const typename R::Point_3* point=object_cast<typename R::Point_3>(&obj);
103 CGAL_kernel_precondition(point!=nullptr);
104 CircleC3 circle = CircleC3(*point, FT(0), Vector_3(FT(1),FT(0),FT(0)));
105 base = Rep(circle.diametral_sphere(), circle.supporting_plane());
106 }
107 }
108
CircleC3(const Point_3 & p,const Point_3 & q,const Point_3 & r)109 CircleC3(const Point_3 &p, const Point_3 &q, const Point_3 &r) {
110 // p, q, r are not collinear
111 CGAL_kernel_precondition(!R().collinear_3_object()(p, q, r));
112 Plane_3 p1 = R().construct_plane_3_object()(p, q, r);
113 Plane_3 p2 = R().construct_bisector_3_object()(p, q);
114 Plane_3 p3 = R().construct_bisector_3_object()(p, r);
115 Object obj = R().intersect_3_object()(p1, p2, p3);
116 // must be a point, otherwise they are collinear
117 const Point_3& center=*object_cast<Point_3>(&obj);
118 FT sqr = R().compute_squared_distance_3_object()(center, r);
119 Sphere_3 s = R().construct_sphere_3_object()(center, sqr);
120 base = Rep(s, p1);
121 }
122
supporting_plane()123 const Plane_3& supporting_plane() const
124 {
125 return boost::get<1>(get_pointee_or_identity(base));
126 }
127
supporting_sphere()128 const Sphere_3& supporting_sphere() const
129 {
130 return diametral_sphere();
131 }
132
center()133 Point_3 center() const
134 {
135 return diametral_sphere().center();
136 }
137
squared_radius()138 FT squared_radius() const
139 {
140 return diametral_sphere().squared_radius();
141 }
142
diametral_sphere()143 const Sphere_3& diametral_sphere() const
144 {
145 return boost::get<0>(get_pointee_or_identity(base));
146 }
147
approximate_area()148 double approximate_area() const
149 {
150 return CGAL_PI * to_double(squared_radius());
151 }
152
approximate_squared_length()153 double approximate_squared_length() const
154 {
155 return CGAL_PI * CGAL_PI * 4.0 * to_double(squared_radius());
156 }
157
area_divided_by_pi()158 FT area_divided_by_pi() const
159 {
160 return squared_radius();
161 }
162
squared_length_divided_by_pi_square()163 FT squared_length_divided_by_pi_square() const
164 {
165 return 4 * squared_radius();
166 }
167
168 // this bbox function
169 // can be optimize by doing different cases
170 // for each variable = 0 (cases with is_zero)
bbox()171 CGAL::Bbox_3 bbox() const
172 {
173 typedef CGAL::Interval_nt<false> Interval;
174 CGAL::Interval_nt<false>::Protector ip;
175 const Sphere_3 &s = diametral_sphere();
176 const FT &sq_r = s.squared_radius();
177 const Point_3 &p = s.center();
178 if(sq_r == FT(0)) return p.bbox();
179 const Plane_3 &plane = supporting_plane();
180 const Interval a = CGAL::to_interval(plane.a());
181 const Interval b = CGAL::to_interval(plane.b());
182 const Interval c = CGAL::to_interval(plane.c());
183 const Interval x = CGAL::to_interval(p.x());
184 const Interval y = CGAL::to_interval(p.y());
185 const Interval z = CGAL::to_interval(p.z());
186 const Interval r2 = CGAL::to_interval(sq_r);
187 const Interval r = CGAL::sqrt(r2); // maybe we can work with r2
188 // in order to save this operation
189 // but if the coefficients are to high
190 // the multiplication would lead to inf
191 // results
192 const Interval a2 = CGAL::square(a);
193 const Interval b2 = CGAL::square(b);
194 const Interval c2 = CGAL::square(c);
195 const Interval sqr_sum = a2 + b2 + c2;
196 const Interval mx = r * CGAL::sqrt((sqr_sum - a2)/sqr_sum);
197 const Interval my = r * CGAL::sqrt((sqr_sum - b2)/sqr_sum);
198 const Interval mz = r * CGAL::sqrt((sqr_sum - c2)/sqr_sum);
199 return CGAL::Bbox_3((x-mx).inf(),(y-my).inf(),(z-mz).inf(),
200 (x+mx).sup(),(y+my).sup(),(z+mz).sup());
201 }
202
203 bool operator==(const CircleC3 &) const;
204 bool operator!=(const CircleC3 &) const;
205
206 bool has_on(const Point_3 &p) const;
207 bool has_on_bounded_side(const Point_3 &p) const;
208 bool has_on_unbounded_side(const Point_3 &p) const;
209 Bounded_side bounded_side(const Point_3 &p) const;
210
is_degenerate()211 bool is_degenerate() const
212 {
213 return diametral_sphere().is_degenerate();
214 }
215
216 };
217
218 template < class R >
219 inline
220 bool
221 CircleC3<R>::
has_on(const typename CircleC3<R>::Point_3 & p)222 has_on(const typename CircleC3<R>::Point_3 &p) const
223 {
224 return R().has_on_3_object()(diametral_sphere(),p) &&
225 R().has_on_3_object()(supporting_plane(),p);
226 }
227
228 template < class R >
229 inline
230 bool
231 CircleC3<R>::
has_on_bounded_side(const typename CircleC3<R>::Point_3 & p)232 has_on_bounded_side(const typename CircleC3<R>::Point_3 &p) const
233 {
234 CGAL_kernel_precondition(R().has_on_3_object()(supporting_plane(), p));
235 return squared_distance(center(),p) < squared_radius();
236 }
237
238 template < class R >
239 inline
240 bool
241 CircleC3<R>::
has_on_unbounded_side(const typename CircleC3<R>::Point_3 & p)242 has_on_unbounded_side(const typename CircleC3<R>::Point_3 &p) const
243 {
244 CGAL_kernel_precondition(R().has_on_3_object()(supporting_plane(), p));
245 return squared_distance(center(),p) > squared_radius();
246 }
247
248 template < class R >
249 CGAL_KERNEL_INLINE
250 Bounded_side
251 CircleC3<R>::
bounded_side(const typename CircleC3<R>::Point_3 & p)252 bounded_side(const typename CircleC3<R>::Point_3 &p) const
253 {
254 CGAL_kernel_precondition(is_degenerate() || R().has_on_3_object()(supporting_plane(), p));
255 return diametral_sphere().bounded_side(p);
256 }
257
258 template < class R >
259 CGAL_KERNEL_INLINE
260 bool
261 CircleC3<R>::operator==(const CircleC3<R> &t) const
262 {
263 if (CGAL::identical(base, t.base))
264 return true;
265 if(!(center() == t.center() &&
266 squared_radius() == t.squared_radius())) return false;
267
268 const typename R::Plane_3 p1 = supporting_plane();
269 const typename R::Plane_3 p2 = t.supporting_plane();
270
271 if(is_zero(p1.a())) {
272 if(!is_zero(p2.a())) return false;
273 if(is_zero(p1.b())) {
274 if(!is_zero(p2.b())) return false;
275 return p1.c() * p2.d() == p1.d() * p2.c();
276 }
277 return (p2.c() * p1.b() == p1.c() * p2.b()) &&
278 (p2.d() * p1.b() == p1.d() * p2.b());
279 }
280 return (p2.b() * p1.a() == p1.b() * p2.a()) &&
281 (p2.c() * p1.a() == p1.c() * p2.a()) &&
282 (p2.d() * p1.a() == p1.d() * p2.a());
283 }
284
285 template < class R >
286 CGAL_KERNEL_INLINE
287 bool
288 CircleC3<R>::operator!=(const CircleC3<R> &t) const
289 {
290 return !(*this == t);
291 }
292
293 } //namespace CGAL
294
295 #endif // CGAL_CARTESIAN_CIRCLEC3_H
296