LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
Circle.cc
Go to the documentation of this file.
1/*
2 * This file is part of sphgeom.
3 *
4 * Developed for the LSST Data Management System.
5 * This product includes software developed by the LSST Project
6 * (http://www.lsst.org).
7 * See the COPYRIGHT file at the top-level directory of this distribution
8 * for details of code ownership.
9 *
10 * This software is dual licensed under the GNU General Public License and also
11 * under a 3-clause BSD license. Recipients may choose which of these licenses
12 * to use; please see the files gpl-3.0.txt and/or bsd_license.txt,
13 * respectively. If you choose the GPL option then the following text applies
14 * (but note that there is still no warranty even if you opt for BSD instead):
15 *
16 * This program is free software: you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation, either version 3 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program. If not, see <http://www.gnu.org/licenses/>.
28 */
29
32
33#include "lsst/sphgeom/Circle.h"
34
35#include <ostream>
36#include <stdexcept>
37
38#include "lsst/sphgeom/Box.h"
39#include "lsst/sphgeom/Box3d.h"
42#include "lsst/sphgeom/codec.h"
43
44
45namespace lsst {
46namespace sphgeom {
47
49 if (a.asRadians() < 0.0) {
50 return -1.0;
51 }
52 if (a.asRadians() >= PI) {
53 return 4.0;
54 }
55 double s = sin(0.5 * a);
56 return 4.0 * s * s;
57}
58
59Angle Circle::openingAngleFor(double squaredChordLength) {
60 // Note: the maximum error in the opening angle (and circle bounding box
61 // width) computations is ~ 2 * MAX_ASIN_ERROR.
62 if (squaredChordLength < 0.0) {
63 return Angle(-1.0);
64 }
65 if (squaredChordLength >= 4.0) {
66 return Angle(PI);
67 }
68 return Angle(2.0 * std::asin(0.5 * std::sqrt(squaredChordLength)));
69}
70
71bool Circle::contains(Circle const & x) const {
72 if (isFull() || x.isEmpty()) {
73 return true;
74 }
75 if (isEmpty() || x.isFull()) {
76 return false;
77 }
78 if (*this == x) {
79 return true;
80 }
81 NormalizedAngle cc(_center, x._center);
82 return _openingAngle >
83 cc + x._openingAngle + 4.0 * Angle(MAX_ASIN_ERROR);
84}
85
86bool Circle::isDisjointFrom(Circle const & x) const {
87 if (isEmpty() || x.isEmpty()) {
88 return true;
89 }
90 if (isFull() || x.isFull()) {
91 return false;
92 }
93 NormalizedAngle cc(_center, x._center);
94 return cc > _openingAngle + x._openingAngle +
95 4.0 * Angle(MAX_ASIN_ERROR);
96}
97
99 *this = contains(x) ? Circle(x) : empty();
100 return *this;
101}
102
104 if (isEmpty() || x.isFull()) {
105 return *this;
106 }
107 if (isFull() || x.isEmpty()) {
108 *this = x;
109 return *this;
110 }
111 Angle a = _openingAngle;
112 Angle b = x._openingAngle;
113 NormalizedAngle cc(_center, x._center);
114 if (cc > a + b + 4.0 * Angle(MAX_ASIN_ERROR)) {
115 // This circle is disjoint from x.
116 *this = empty();
117 return *this;
118 }
119 // The circles (nearly) intersect, or one contains the other.
120 // For now, take the easy route and just use the smaller of
121 // the two circles as a bound on their intersection.
122 //
123 // TODO(smm): Compute the minimal bounding circle.
124 if (b < a) {
125 *this = x;
126 }
127 return *this;
128}
129
131 // For any circle c and unit vector x, c.expandTo(x).contains(x)
132 // should return true.
133 if (isEmpty()) {
134 *this = Circle(x);
135 } else if (!contains(x)) {
136 // Compute the normal vector for the plane defined by _center and x.
138 // The minimal bounding circle (MBC) includes unit vectors on the plane
139 // with normal n that span from _center.rotatedAround(n, -_openingAngle)
140 // to x. The MBC center is the midpoint of this interval.
141 NormalizedAngle cx(_center, x);
142 Angle o = 0.5 * (cx + _openingAngle);
143 Angle r = 0.5 * (cx - _openingAngle);
144 // Rotate _center by angle r around n to obtain the MBC center. This is
145 // done using Rodriques' formula, simplified by taking advantage of the
146 // orthogonality of _center and n.
147 _center = UnitVector3d(_center * cos(r) + n.cross(_center) * sin(r));
148 _squaredChordLength = squaredChordLengthFor(o + Angle(MAX_ASIN_ERROR));
149 _openingAngle = o + Angle(MAX_ASIN_ERROR);
150 }
151 return *this;
152}
153
155 if (isEmpty() || x.isFull()) {
156 *this = x;
157 return *this;
158 }
159 if (x.isEmpty() || isFull()) {
160 return *this;
161 }
162 NormalizedAngle cc(_center, x._center);
163 if (cc + x._openingAngle + 4.0 * Angle(MAX_ASIN_ERROR) <= _openingAngle) {
164 // This circle contains x.
165 return *this;
166 }
167 if (cc + _openingAngle + 4.0 * Angle(MAX_ASIN_ERROR) <= x._openingAngle) {
168 // x contains this circle.
169 *this = x;
170 return *this;
171 }
172 // The circles intersect or are disjoint.
173 Angle o = 0.5 * (cc + _openingAngle + x._openingAngle);
174 if (o + 2.0 * Angle(MAX_ASIN_ERROR) >= Angle(PI)) {
175 *this = full();
176 return *this;
177 }
178 // Compute the normal vector for the plane defined by the circle centers.
179 UnitVector3d n = UnitVector3d::orthogonalTo(_center, x._center);
180 // The minimal bounding circle (MBC) includes unit vectors on the plane
181 // with normal n that span from _center.rotatedAround(n, -_openingAngle)
182 // to x._center.rotatedAround(n, x._openingAngle). The MBC center is the
183 // midpoint of this interval.
184 Angle r = o - _openingAngle;
185 // Rotate _center by angle r around n to obtain the MBC center. This is
186 // done using Rodriques' formula, simplified by taking advantage of the
187 // orthogonality of _center and n.
188 _center = UnitVector3d(_center * cos(r) + n.cross(_center) * sin(r));
189 _squaredChordLength = squaredChordLengthFor(o + Angle(MAX_ASIN_ERROR));
190 _openingAngle = o + Angle(MAX_ASIN_ERROR);
191 return *this;
192}
193
195 if (!isEmpty() && !isFull() &&
196 (r.asRadians() > 0.0 || r.asRadians() < 0.0)) {
197 Angle o = _openingAngle + r;
198 _squaredChordLength = squaredChordLengthFor(o);
199 _openingAngle = o;
200 }
201 return *this;
202}
203
205 if (isEmpty()) {
206 // The complement of an empty circle is a full circle.
207 _squaredChordLength = 4.0;
208 _openingAngle = Angle(PI);
209 } else if (isFull()) {
210 // The complement of a full circle is an empty circle.
211 _squaredChordLength = -1.0;
212 _openingAngle = Angle(-1.0);
213 } else {
214 _center = -_center;
215 _squaredChordLength = 4.0 - _squaredChordLength;
216 _openingAngle = Angle(PI) - _openingAngle;
217 }
218 return *this;
219}
220
222 LonLat c(_center);
223 Angle h = _openingAngle + 2.0 * Angle(MAX_ASIN_ERROR);
226 return Box(c, w, h);
227}
228
230 static double const MAX_BOUNDARY_ERROR = 6.2e-16; // > 5.5ε, where ε = 2^-53
231 if (isEmpty()) {
232 return Box3d();
233 }
234 if (isFull()) {
236 }
237 // Given circle center c and standard basis vector eᵢ, to check whether
238 // ±eᵢ is inside the circle we need to check that (c ∓ eᵢ)·(c ∓ eᵢ) ≤ s.
239 // Since c·c = 1, eᵢ·eᵢ = 1 (c and eᵢ are unit vectors) this is the
240 // same as checking that 2 ∓ 2c·eᵢ ≤ s, or 2 ∓ 2cᵢ ≤ s, where cᵢ is
241 // the i-th component of c.
242 //
243 // Besides any standard basis vectors inside the circle, the bounding box
244 // must also include the circle boundary. To find the extent of this
245 // boundary along a particular axis, note that we can write the i-th
246 // component of the circle center vector as the sine of a latitude angle
247 // (treating the i-th standard basis vector as "north"). So given a circle
248 // opening angle θ, the desired extent is
249 //
250 // [min(sin(asin(cᵢ) ± θ)), max(sin(asin(cᵢ) ± θ))]
251 //
252 // which can be simplified using the usual trigonometric identities to
253 // arrive at the code below.
254 Interval1d e[3];
255 double s = sin(_openingAngle);
256 double c = cos(_openingAngle);
257 for (int i = 0; i < 3; ++i) {
258 double ci = _center(i);
259 double di = std::sqrt(1.0 - ci * ci);
260 double bmin = 1.0, bmax = -1.0;
261 if (2.0 - 2.0 * ci <= _squaredChordLength) {
262 bmax = 1.0;
263 }
264 if (2.0 + 2.0 * ci <= _squaredChordLength) {
265 bmin = -1.0;
266 }
267 double b0 = ci * c + di * s;
268 bmax = std::max(bmax, b0 + MAX_BOUNDARY_ERROR);
269 bmin = std::min(bmin, b0 - MAX_BOUNDARY_ERROR);
270 double b1 = ci * c - di * s;
271 bmax = std::max(bmax, b1 + MAX_BOUNDARY_ERROR);
272 bmin = std::min(bmin, b1 - MAX_BOUNDARY_ERROR);
273 e[i] = Interval1d(std::max(-1.0, bmin), std::min(1.0, bmax));
274 }
275 return Box3d(e[0], e[1], e[2]);
276}
277
279 if (contains(v)) {
280 return CONTAINS;
281 } else if (isEmpty()) {
282 return DISJOINT | WITHIN;
283 }
284 return DISJOINT;
285}
286
288 // Box-Circle relations are implemented by Box.
289 return invert(b.relate(*this));
290}
291
293 if (isEmpty()) {
294 if (c.isEmpty()) {
295 return CONTAINS | DISJOINT | WITHIN;
296 }
297 return DISJOINT | WITHIN;
298 } else if (c.isEmpty()) {
299 return CONTAINS | DISJOINT;
300 }
301 // Neither circle is empty.
302 if (isFull()) {
303 if (c.isFull()) {
304 return CONTAINS | WITHIN;
305 }
306 return CONTAINS;
307 } else if (c.isFull()) {
308 return WITHIN;
309 }
310 // Special case equality, which can be missed by logic below due to
311 // round-off error.
312 if (*this == c) {
313 return CONTAINS | WITHIN;
314 }
315 // Neither circle is full.
316 NormalizedAngle cc(_center, c._center);
317 if (cc > _openingAngle + c._openingAngle + 4.0 * Angle(MAX_ASIN_ERROR)) {
318 return DISJOINT;
319 }
320 if (cc + c._openingAngle + 4.0 * Angle(MAX_ASIN_ERROR) <= _openingAngle) {
321 return CONTAINS;
322 } else if (cc + _openingAngle + 4.0 * Angle(MAX_ASIN_ERROR) <=
323 c._openingAngle) {
324 return WITHIN;
325 }
326 return INTERSECTS;
327}
328
330 // ConvexPolygon-Circle relations are implemented by ConvexPolygon.
331 return invert(p.relate(*this));
332}
333
335 // Ellipse-Circle relations are implemented by Ellipse.
336 return invert(e.relate(*this));
337}
338
341 uint8_t tc = TYPE_CODE;
342 buffer.reserve(ENCODED_SIZE);
343 buffer.push_back(tc);
344 encodeDouble(_center.x(), buffer);
345 encodeDouble(_center.y(), buffer);
346 encodeDouble(_center.z(), buffer);
347 encodeDouble(_squaredChordLength, buffer);
348 encodeDouble(_openingAngle.asRadians(), buffer);
349 return buffer;
350}
351
352std::unique_ptr<Circle> Circle::decode(uint8_t const * buffer, size_t n) {
353 if (buffer == nullptr || n != ENCODED_SIZE || *buffer != TYPE_CODE) {
354 throw std::runtime_error("Byte-string is not an encoded Circle");
355 }
357 ++buffer;
358 double x = decodeDouble(buffer); buffer += 8;
359 double y = decodeDouble(buffer); buffer += 8;
360 double z = decodeDouble(buffer); buffer += 8;
361 double squaredChordLength = decodeDouble(buffer); buffer += 8;
362 double openingAngle = decodeDouble(buffer); buffer += 8;
363 circle->_center = UnitVector3d::fromNormalized(x, y, z);
364 circle->_squaredChordLength = squaredChordLength;
365 circle->_openingAngle = Angle(openingAngle);
366 return circle;
367}
368
370 char tail[32];
371 std::snprintf(tail, sizeof(tail), ", %.17g]}", c.getSquaredChordLength());
372 return os << "{\"Circle\": [" << c.getCenter() << tail;
373}
374
375}} // namespace lsst::sphgeom
This file declares a class for representing axis-aligned bounding boxes in ℝ³.
This file declares a class for representing circular regions on the unit sphere.
This file declares a class for representing convex polygons with great circle edges on the unit spher...
double z
Definition Match.cc:44
int y
Definition SpanSet.cc:48
table::Key< int > b
T asin(T... args)
Angle represents an angle in radians.
Definition Angle.h:50
double asRadians() const
asRadians returns the value of this angle in units of radians.
Definition Angle.h:92
Box3d represents a box in ℝ³.
Definition Box3d.h:49
static Box3d aroundUnitSphere()
aroundUnitSphere returns a minimal Box3d containing the unit sphere.
Definition Box3d.h:63
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition Box.h:61
static NormalizedAngle halfWidthForCircle(Angle r, Angle lat)
halfWidthForCircle computes the half-width of bounding boxes for circles with radius r and centers at...
Definition Box.cc:50
Circle is a circular region on the unit sphere that contains its boundary.
Definition Circle.h:53
bool contains(Circle const &x) const
contains returns true if the intersection of this circle and x is equal to x.
Definition Circle.cc:71
static Angle openingAngleFor(double squaredChordLength)
openingAngleFor computes and returns the angular separation between points in S² that are separated b...
Definition Circle.cc:59
bool isEmpty() const
Definition Circle.h:116
Box getBoundingBox() const override
getBoundingBox returns a bounding-box for this region.
Definition Circle.cc:221
static Circle empty()
Definition Circle.h:57
Circle & clipTo(UnitVector3d const &x)
Definition Circle.cc:98
Circle & dilateBy(Angle r)
If r is positive, dilateBy increases the opening angle of this circle to include all points within an...
Definition Circle.cc:194
static Circle full()
Definition Circle.h:59
static std::unique_ptr< Circle > decode(std::vector< uint8_t > const &s)
Definition Circle.h:258
bool isFull() const
Definition Circle.h:121
std::vector< uint8_t > encode() const override
encode serializes this region into an opaque byte string.
Definition Circle.cc:339
bool isDisjointFrom(UnitVector3d const &x) const
Definition Circle.h:144
Relationship relate(UnitVector3d const &v) const
Definition Circle.cc:278
double getSquaredChordLength() const
getSquaredChordLength returns the squared length of chords between the circle center and points on th...
Definition Circle.h:130
UnitVector3d const & getCenter() const
getCenter returns the center of this circle as a unit vector.
Definition Circle.h:125
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.
Definition Circle.cc:229
Circle & complement()
complement sets this circle to the closure of its complement.
Definition Circle.cc:204
Circle & expandTo(UnitVector3d const &x)
Definition Circle.cc:130
Circle()
This constructor creates an empty circle.
Definition Circle.h:72
static double squaredChordLengthFor(Angle openingAngle)
squaredChordLengthFor computes and returns the squared chord length between points in S² that are sep...
Definition Circle.cc:48
static constexpr uint8_t TYPE_CODE
Definition Circle.h:55
ConvexPolygon is a closed convex polygon on the unit sphere.
Relationship relate(Region const &r) const override
Ellipse is an elliptical region on the sphere.
Definition Ellipse.h:177
Relationship relate(Region const &r) const override
Definition Ellipse.h:296
Interval1d represents closed intervals of ℝ.
Definition Interval1d.h:47
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition LonLat.h:55
Angle getLat() const
Definition LonLat.h:97
NormalizedAngle is an angle that lies in the range [0, 2π), with one exception - a NormalizedAngle ca...
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
static UnitVector3d orthogonalTo(Vector3d const &v)
orthogonalTo returns an arbitrary unit vector that is orthogonal to v.
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this unit vector and v.
This file contains simple helper functions for encoding and decoding primitive types to/from byte str...
T snprintf(T... args)
T max(T... args)
T min(T... args)
constexpr double MAX_ASIN_ERROR
Definition constants.h:52
void encodeDouble(double item, std::vector< uint8_t > &buffer)
encodeDouble appends an IEEE double in little-endian byte order to the end of buffer.
Definition codec.h:56
std::ostream & operator<<(std::ostream &, Angle const &)
Definition Angle.cc:41
double sin(Angle const &a)
Definition Angle.h:109
double decodeDouble(uint8_t const *buffer)
decodeDouble extracts an IEEE double from the 8 byte little-endian byte sequence in buffer.
Definition codec.h:76
double cos(Angle const &a)
Definition Angle.h:110
constexpr double PI
Definition constants.h:43
Relationship invert(Relationship r)
Given the relationship between two sets A and B (i.e.
T push_back(T... args)
T reserve(T... args)
This file declares a class for representing longitude/latitude angle boxes on the unit sphere.
This file declares a class for representing elliptical regions on the unit sphere.
T sqrt(T... args)
double w
Definition CoaddPsf.cc:69