1 // Copyright (c) 2000,2001
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/Kernel_d/include/CGAL/Kernel_d/Sphere_d.h $
11 // $Id: Sphere_d.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
12 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
13 //
14 //
15 // Author(s) : Michael Seel
16
17 #ifndef CGAL_SPHERE_D_H
18 #define CGAL_SPHERE_D_H
19
20 #include <CGAL/basic.h>
21 #include <vector>
22 #include <CGAL/Dimension.h>
23 #include <CGAL/use.h>
24 #include <CGAL/Handle.h>
25 #include <CGAL/Kernel_d/Tuple_d.h>
26 #include <CGAL/Kernel_d/Aff_transformation_d.h>
27 #include <CGAL/Kernel_d/Vector_d.h>
28
29 namespace CGAL {
30
31 template <class R> class Sphere_d;
32 template <class R> bool equal_as_sets(const Sphere_d<R>&, const Sphere_d<R>&);
33
34 template <class R>
35 class Sphere_d_rep {
36
37 typedef typename R::Point_d Point_d;
38
39 friend class Sphere_d<R>;
40 friend bool equal_as_sets <>
41 (const Sphere_d<R>&, const Sphere_d<R>&);
42
43 std::vector< Point_d > P; // d+1 defining points, index range 0-d
44 Orientation orient; // orientation(P)
45 Point_d* cp; // pointer to center (lazy calculated)
46
47 public:
Sphere_d_rep()48 Sphere_d_rep() : cp(0) {}
Sphere_d_rep(int d)49 Sphere_d_rep(int d) : P(d), cp(0) {}
50
51 template <class ForwardIterator>
Sphere_d_rep(int d,ForwardIterator first,ForwardIterator last)52 Sphere_d_rep(int d, ForwardIterator first, ForwardIterator last) :
53 P(first,last), cp(0)
54 { TUPLE_DIM_CHECK(P.begin(),P.end(),Sphere_d);
55 CGAL_assertion(d+1==int(P.size()));
56 CGAL_USE(d);
57 typename R::Orientation_d orientation_;
58 orient = orientation_(P.begin(),P.end()); }
59
~Sphere_d_rep()60 ~Sphere_d_rep() { if (cp) delete cp; }
61
62 }; // Sphere_d_rep<R>
63
64 /*{\Manpage {Sphere_d}{R}{Simple Spheres}{S}}*/
65
66 template <class R_>
67 class Sphere_d : public Handle_for< Sphere_d_rep<R_> > {
68
69 /*{\Mdefinition
70 An instance $S$ of the data type |Sphere_d| is an oriented sphere in
71 some $d$-dimensional space. A sphere is defined by $d+1$ points with
72 rational coordinates (class |Point_d<R>|). We use $A$ to denote the
73 array of the defining points. A set $A$ of defining points is
74 \emph{legal} if either the points are affinely independent or if the
75 points are all equal. Only a legal set of points defines a sphere in
76 the geometric sense and hence many operations on spheres require the
77 set of defining points to be legal. The orientation of $S$ is equal
78 to the orientation of the defining points, i.e., |orientation(A)|. }*/
79
80 public:
81 typedef CGAL::Dynamic_dimension_tag Ambient_dimension;
82 typedef CGAL::Dynamic_dimension_tag Feature_dimension;
83
84 /*{\Mtypes 4}*/
85
86 typedef Sphere_d_rep<R_> Rep;
87 typedef Handle_for<Rep> Base;
88 typedef Sphere_d<R_> Self;
89 typedef typename R_::Point_d Point_d;
90
91 using Base::ptr;
92
Sphere_d(const Base & b)93 Sphere_d(const Base& b) : Base(b) {}
94
95 typedef R_ R;
96 /*{\Mtypemember the representation type.}*/
97
98 typedef typename R::RT RT;
99 /*{\Mtypemember the ring type.}*/
100
101 typedef typename R::FT FT;
102 /*{\Mtypemember the field type.}*/
103
104 typedef typename R::LA LA;
105 /*{\Mtypemember the linear algebra layer.}*/
106
107 typedef typename std::vector< Point_d >::const_iterator point_iterator;
108 /*{\Mtypemember a read-only iterator for points defining the sphere.}*/
109
110 /*{\Mcreation 4}*/
111
112 Sphere_d(int d = 0) : Base( Rep(d+1) )
113 /*{\Mcreate introduces a variable |\Mvar| of type |\Mname|. |\Mvar|
114 is initialized to the empty sphere centered at the origin of
115 $d$-dimensional space. }*/
116 {
117 Point_d p(d);
118 for (int i = 0; i <= d; i++) ptr()->P[i] = p;
119 ptr()->orient = ZERO;
120 ptr()->cp = new Point_d(p);
121 }
122
123 template <class ForwardIterator>
Sphere_d(int d,ForwardIterator first,ForwardIterator last)124 Sphere_d(int d, ForwardIterator first, ForwardIterator last) :
125 /*{\Mcreate introduces a variable |\Mvar| of type |\Mname|. |\Mvar| is
126 initialized to the sphere through the points in |A = set [first,last)|.
127 \precond $A$ consists of $d+1$ $d$-dimensional points.}*/
128 Base( Rep(d,first,last) ) {}
129
130
131 /*{\Moperations 4 3}*/
132
dimension()133 int dimension() const
134 /*{\Mop returns the dimension of |\Mvar|.}*/
135 { return static_cast<int>(ptr()->P.size()) - 1; }
136
point(int i)137 Point_d point(int i) const
138 /*{\Mop returns the $i$-th defining point. \precond $0 \le i \le |dim|$.}*/
139 { CGAL_assertion_msg((0<=i && i<=dimension()),
140 "Sphere_d::point(): index out of range.");
141 return ptr()->P[i];
142 }
143
points_begin()144 point_iterator points_begin() const { return ptr()->P.begin(); }
145 /*{\Mop returns an iterator pointing to the first defining point.}*/
146
points_end()147 point_iterator points_end() const { return ptr()->P.end(); }
148 /*{\Mop returns an iterator pointing beyond the last defining point.}*/
149
is_degenerate()150 bool is_degenerate() const { return (ptr()->orient == CGAL::ZERO); }
151 /*{\Mop returns true iff the defining points are not full dimenional.}*/
152
is_legal()153 bool is_legal() const
154 /*{\Mop returns true iff the set of defining points is legal.
155 A set of defining points is legal iff their orientation is
156 non-zero or if they are all equal.}*/
157 { if (ptr()->orient != ZERO ) return true;
158 const std::vector< Point_d >& A = ptr()->P;
159 Point_d po = A[0];
160 for (int i = 1; i < int(A.size()); ++i)
161 if (A[i] != po) return false;
162 return true;
163 }
164
center()165 Point_d center() const
166 /*{\Mop returns the center of |\Mvar|. \precond The orientation
167 of |\Mvar| is non-zero. }*/
168 {
169 if (ptr()->cp == 0) {
170 if (ptr()->orient == 0) {
171 const std::vector< Point_d >& A = ptr()->P;
172 Point_d po = A[0];
173 for (int i = 1; i < int(A.size()); ++i)
174 if (A[i] != po)
175 CGAL_error_msg("Sphere_d::center(): points are illegal.");
176 const_cast<Self&>(*this).ptr()->cp = new Point_d(A[0]);
177 return *(ptr()->cp);
178 }
179 typename R::Center_of_sphere_d center_of_sphere_;
180 const_cast<Self&>(*this).ptr()->cp =
181 new Point_d(center_of_sphere_(points_begin(),points_end()));
182 }
183 return *(ptr()->cp);
184 }
185
186
187
squared_radius()188 FT squared_radius() const
189 /*{\Mop returns the squared radius of the sphere.}*/
190 { if (is_degenerate()) return 0;
191 return (point(0)-center()).squared_length();
192 }
193
orientation()194 Orientation orientation() const { return ptr()->orient; }
195 /*{\Mop returns the orientation of |\Mvar|.}*/
196
oriented_side(const Point_d & p)197 Oriented_side oriented_side(const Point_d& p) const
198 /*{\Mop returns either the constant |ON_ORIENTED_BOUNDARY|,
199 |ON_POSITIVE_SIDE|, or |ON_NEGATIVE_SIDE|, iff p lies on the boundary,
200 properly on the positive side, or properly on the negative side of
201 circle, resp.}*/
202 { typename R::Side_of_oriented_sphere_d side;
203 return side(points_begin(),points_end(),p); }
204
bounded_side(const Point_d & p)205 Bounded_side bounded_side(const Point_d& p) const
206 /*{\Mop returns |ON_BOUNDED_SIDE|, |ON_BOUNDARY|, or
207 |ON_UNBOUNDED_SIDE| iff p lies properly inside, on
208 the boundary, or properly outside of circle, resp.}*/
209 { typename R::Side_of_bounded_sphere_d side;
210 return side(points_begin(),points_end(),p); }
211
has_on_positive_side(const Point_d & p)212 bool has_on_positive_side (const Point_d& p) const
213 /*{\Mop returns |\Mvar.oriented_side(p)==ON_POSITIVE_SIDE|.}*/
214 { return oriented_side(p) == ON_POSITIVE_SIDE; }
215
has_on_negative_side(const Point_d & p)216 bool has_on_negative_side (const Point_d& p) const
217 /*{\Mop returns |\Mvar.oriented_side(p)==ON_NEGATIVE_SIDE|.}*/
218 { return oriented_side(p) == ON_NEGATIVE_SIDE; }
219
has_on_boundary(const Point_d & p)220 bool has_on_boundary (const Point_d& p) const
221 /*{\Mop returns |\Mvar.oriented_side(p)==ON_ORIENTED_BOUNDARY|,
222 which is the same as |\Mvar.bounded_side(p)==ON_BOUNDARY|.}*/
223 { return oriented_side(p) == ON_ORIENTED_BOUNDARY; }
224
has_on_bounded_side(const Point_d & p)225 bool has_on_bounded_side (const Point_d& p) const
226 /*{\Mop returns |\Mvar.bounded_side(p)==ON_BOUNDED_SIDE|.}*/
227 { return (int(ptr()->orient) * int(oriented_side(p))) > 0 ; }
228
has_on_unbounded_side(const Point_d & p)229 bool has_on_unbounded_side (const Point_d& p) const
230 /*{\Mop returns |\Mvar.bounded_side(p)==ON_UNBOUNDED_SIDE|.}*/
231 { return (int(ptr()->orient) * int(oriented_side(p))) < 0; }
232
opposite()233 Sphere_d<R> opposite() const
234 /*{\Mop returns the sphere with the same center and squared
235 radius as |\Mvar| but with opposite orientation.}*/
236 { CGAL_assertion(dimension()>1);
237 if (is_degenerate()) return *this;
238 std::vector< Point_d > V(points_begin(),points_end());
239 std::swap(V[0],V[1]);
240 return Sphere_d<R>(dimension(),V.begin(),V.end());
241 }
242
transform(const Aff_transformation_d<R> & t)243 Sphere_d<R> transform(const Aff_transformation_d<R>& t) const
244 /*{\Mopl returns $t(s)$. }*/
245 { std::vector< Point_d > B(points_begin(),points_end());
246 typename std::vector< Point_d >::iterator it;
247 for (it = B.begin(); it != B.end(); ++it)
248 *it = it->transform(t);
249 return Sphere_d<R>(dimension(),B.begin(),B.end());
250 }
251
252 Sphere_d<R> operator+(const Vector_d<R>& v) const
253 /*{\Mbinop returns the sphere translated by |v|. }*/
254 { std::vector< Point_d > B(points_begin(),points_end());
255 typename std::vector< Point_d >::iterator it;
256 for (it = B.begin(); it != B.end(); ++it) it->operator+=(v);
257 return Sphere_d<R>(dimension(),B.begin(),B.end());
258 }
259
260 bool operator==(const Sphere_d<R>& D) const
261 { if (this->identical(D)) return true;
262 if (dimension() != D.dimension()) return false;
263 return (center()==D.center() &&
264 squared_radius() == D.squared_radius() &&
265 orientation() == D.orientation());
266 }
267
268 bool operator!=(const Sphere_d<R>& D) const
269 { return !operator==(D); }
270
271
272 }; // end of class Sphere_d
273
274 /*{\Mtext \headerline{Non-Member Functions} }*/
275 template <class R>
weak_equality(const Sphere_d<R> & s1,const Sphere_d<R> & s2)276 bool weak_equality(const Sphere_d<R>& s1, const Sphere_d<R>& s2)
277 /*{\Mfunc Test for equality as unoriented spheres.}*/
278 { if (s1.identical(s2)) return true;
279 if (s1.dimension() != s2.dimension()) return false;
280 return (s1.center()==s2.center() &&
281 s1.squared_radius() == s2.squared_radius());
282 }
283
284 /*{\Mimplementation Spheres are implemented by a vector of points as
285 an item type. All operations like creation, initialization, tests,
286 input and output of a sphere $s$ take time
287 $O(|s.dimension()|)$. |dimension()|, point access take constant time.
288 The space requirement for spheres is $O(|s.dimension()|)$
289 neglecting the storage room of the points.}*/
290
291 template <class R>
292 std::ostream& operator<<(std::ostream& os, const CGAL::Sphere_d<R>& s)
293 { typedef typename Sphere_d<R>::point_iterator iterator;
294 os << s.dimension()+1 << " ";
295 for (iterator it=s.points_begin(); it!=s.points_end(); ++it)
296 os << *it << " ";
297 return os;
298 }
299
300 template <class R> std::istream&
301 operator>>(std::istream& is, CGAL::Sphere_d<R>& s)
302 { int d; is >> d;
303 std::vector< Point_d<R> > V(d);
304 Point_d<R> p;
305 while ( d-- ) {
306 if (!(is >> p)) return is;
307 V[d] = p;
308 }
309 s = Sphere_d<R>(static_cast<int>(V.size())-1, V.begin(), V.end() );
310 return is;
311 }
312
313
314 /*
315 The center is only defined if the set of defining points are
316 legal. If the defining points are all equal the sphere is trivial. So
317 assume otherwise. Then the center $c$ is the unique point with equal
318 distance to all the defining points. A point $c$ has equal distance to
319 point $p_0$ and $p_i$ if it lies on the perpendicual bisector of $p_d$
320 and $p_i$, i.e., if it satiesfies the plane equation $(p_i - p_d)^T c
321 = (p_i - p_d) (p_i + p_d)/2$. Note that $p_i - p_d$ is the normal
322 vector of the bisector hyperplane and $(p_i + p_d)/2$ is the midpoint
323 of $p_d$ and $p_i$. The equation above translates into the equation \[
324 \sum_{0 \le j < d} 2*p_{dd}p_{id}(p_{ij}p_{dd} - p_{dj}p_{id})c_j/c_d
325 = \sum_{0 \le j < d} (p_{ij}p_{dd} - p_{dj}p_{id})(p_{ij}p_{dd} +
326 p_{dj}p_{id}) \] for the homogeneous coordinates of the points and the
327 center. We may tentatively assume that $c_d = 1$, solve the
328 corresponding linear system, and then define the center.
329 */
330
331 } //namespace CGAL
332
333 #endif // CGAL_SPHERE_D_H
334 //----------------------- end of file ----------------------------------
335