LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Circle.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2014-2015 AURA/LSST.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <https://www.lsstcorp.org/LegalNotices/>.
21  */
22 
25 
26 #include "lsst/sphgeom/Circle.h"
27 
28 #include <ostream>
29 #include <stdexcept>
30 
31 #include "lsst/sphgeom/Box.h"
32 #include "lsst/sphgeom/Box3d.h"
34 #include "lsst/sphgeom/Ellipse.h"
35 #include "lsst/sphgeom/codec.h"
36 
37 
38 namespace lsst {
39 namespace sphgeom {
40 
42  if (a.asRadians() < 0.0) {
43  return -1.0;
44  }
45  if (a.asRadians() >= PI) {
46  return 4.0;
47  }
48  double s = sin(0.5 * a);
49  return 4.0 * s * s;
50 }
51 
52 Angle Circle::openingAngleFor(double squaredChordLength) {
53  // Note: the maximum error in the opening angle (and circle bounding box
54  // width) computations is ~ 2 * MAX_ASIN_ERROR.
55  if (squaredChordLength < 0.0) {
56  return Angle(-1.0);
57  }
58  if (squaredChordLength >= 4.0) {
59  return Angle(PI);
60  }
61  return Angle(2.0 * std::asin(0.5 * std::sqrt(squaredChordLength)));
62 }
63 
64 bool Circle::contains(Circle const & x) const {
65  if (isFull() || x.isEmpty()) {
66  return true;
67  }
68  if (isEmpty() || x.isFull()) {
69  return false;
70  }
71  NormalizedAngle cc(_center, x._center);
72  return _openingAngle >
73  cc + x._openingAngle + 4.0 * Angle(MAX_ASIN_ERROR);
74 }
75 
76 bool Circle::isDisjointFrom(Circle const & x) const {
77  if (isEmpty() || x.isEmpty()) {
78  return true;
79  }
80  if (isFull() || x.isFull()) {
81  return false;
82  }
83  NormalizedAngle cc(_center, x._center);
84  return cc > _openingAngle + x._openingAngle +
85  4.0 * Angle(MAX_ASIN_ERROR);
86 }
87 
89  *this = contains(x) ? Circle(x) : empty();
90  return *this;
91 }
92 
94  if (isEmpty() || x.isFull()) {
95  return *this;
96  }
97  if (isFull() || x.isEmpty()) {
98  *this = x;
99  return *this;
100  }
101  Angle a = _openingAngle;
102  Angle b = x._openingAngle;
103  NormalizedAngle cc(_center, x._center);
104  if (cc > a + b + 4.0 * Angle(MAX_ASIN_ERROR)) {
105  // This circle is disjoint from x.
106  *this = empty();
107  return *this;
108  }
109  // The circles (nearly) intersect, or one contains the other.
110  // For now, take the easy route and just use the smaller of
111  // the two circles as a bound on their intersection.
112  //
113  // TODO(smm): Compute the minimal bounding circle.
114  if (b < a) {
115  *this = x;
116  }
117  return *this;
118 }
119 
121  // For any circle c and unit vector x, c.expandTo(x).contains(x)
122  // should return true.
123  if (isEmpty()) {
124  *this = Circle(x);
125  } else if (!contains(x)) {
126  // Compute the normal vector for the plane defined by _center and x.
128  // The minimal bounding circle (MBC) includes unit vectors on the plane
129  // with normal n that span from _center.rotatedAround(n, -_openingAngle)
130  // to x. The MBC center is the midpoint of this interval.
131  NormalizedAngle cx(_center, x);
132  Angle o = 0.5 * (cx + _openingAngle);
133  Angle r = 0.5 * (cx - _openingAngle);
134  // Rotate _center by angle r around n to obtain the MBC center. This is
135  // done using Rodriques' formula, simplified by taking advantage of the
136  // orthogonality of _center and n.
137  _center = UnitVector3d(_center * cos(r) + n.cross(_center) * sin(r));
138  _squaredChordLength = squaredChordLengthFor(o + Angle(MAX_ASIN_ERROR));
139  _openingAngle = o + Angle(MAX_ASIN_ERROR);
140  }
141  return *this;
142 }
143 
145  if (isEmpty() || x.isFull()) {
146  *this = x;
147  return *this;
148  }
149  if (x.isEmpty() || isFull()) {
150  return *this;
151  }
152  NormalizedAngle cc(_center, x._center);
153  if (cc + x._openingAngle + 4.0 * Angle(MAX_ASIN_ERROR) <= _openingAngle) {
154  // This circle contains x.
155  return *this;
156  }
157  if (cc + _openingAngle + 4.0 * Angle(MAX_ASIN_ERROR) <= x._openingAngle) {
158  // x contains this circle.
159  *this = x;
160  return *this;
161  }
162  // The circles intersect or are disjoint.
163  Angle o = 0.5 * (cc + _openingAngle + x._openingAngle);
164  if (o + 2.0 * Angle(MAX_ASIN_ERROR) >= Angle(PI)) {
165  *this = full();
166  return *this;
167  }
168  // Compute the normal vector for the plane defined by the circle centers.
169  UnitVector3d n = UnitVector3d::orthogonalTo(_center, x._center);
170  // The minimal bounding circle (MBC) includes unit vectors on the plane
171  // with normal n that span from _center.rotatedAround(n, -_openingAngle)
172  // to x._center.rotatedAround(n, x._openingAngle). The MBC center is the
173  // midpoint of this interval.
174  Angle r = o - _openingAngle;
175  // Rotate _center by angle r around n to obtain the MBC center. This is
176  // done using Rodriques' formula, simplified by taking advantage of the
177  // orthogonality of _center and n.
178  _center = UnitVector3d(_center * cos(r) + n.cross(_center) * sin(r));
179  _squaredChordLength = squaredChordLengthFor(o + Angle(MAX_ASIN_ERROR));
180  _openingAngle = o + Angle(MAX_ASIN_ERROR);
181  return *this;
182 }
183 
185  if (!isEmpty() && !isFull() &&
186  (r.asRadians() > 0.0 || r.asRadians() < 0.0)) {
187  Angle o = _openingAngle + r;
188  _squaredChordLength = squaredChordLengthFor(o);
189  _openingAngle = o;
190  }
191  return *this;
192 }
193 
195  if (isEmpty()) {
196  // The complement of an empty circle is a full circle.
197  _squaredChordLength = 4.0;
198  _openingAngle = Angle(PI);
199  } else if (isFull()) {
200  // The complement of a full circle is an empty circle.
201  _squaredChordLength = -1.0;
202  _openingAngle = Angle(-1.0);
203  } else {
204  _center = -_center;
205  _squaredChordLength = 4.0 - _squaredChordLength;
206  _openingAngle = Angle(PI) - _openingAngle;
207  }
208  return *this;
209 }
210 
212  LonLat c(_center);
213  Angle h = _openingAngle + 2.0 * Angle(MAX_ASIN_ERROR);
216  return Box(c, w, h);
217 }
218 
220  static double const MAX_BOUNDARY_ERROR = 6.2e-16; // > 5.5ε, where ε = 2^-53
221  if (isEmpty()) {
222  return Box3d();
223  }
224  if (isFull()) {
225  return Box3d::aroundUnitSphere();
226  }
227  // Given circle center c and standard basis vector eᵢ, to check whether
228  // ±eᵢ is inside the circle we need to check that (c ∓ eᵢ)·(c ∓ eᵢ) ≤ s.
229  // Since c·c = 1, eᵢ·eᵢ = 1 (c and eᵢ are unit vectors) this is the
230  // same as checking that 2 ∓ 2c·eᵢ ≤ s, or 2 ∓ 2cᵢ ≤ s, where cᵢ is
231  // the i-th component of c.
232  //
233  // Besides any standard basis vectors inside the circle, the bounding box
234  // must also include the circle boundary. To find the extent of this
235  // boundary along a particular axis, note that we can write the i-th
236  // component of the circle center vector as the sine of a latitude angle
237  // (treating the i-th standard basis vector as "north"). So given a circle
238  // opening angle θ, the desired extent is
239  //
240  // [min(sin(asin(cᵢ) ± θ)), max(sin(asin(cᵢ) ± θ))]
241  //
242  // which can be simplified using the usual trigonometric identities to
243  // arrive at the code below.
244  Interval1d e[3];
245  double s = sin(_openingAngle);
246  double c = cos(_openingAngle);
247  for (int i = 0; i < 3; ++i) {
248  double ci = _center(i);
249  double di = std::sqrt(1.0 - ci * ci);
250  double bmin = 1.0, bmax = -1.0;
251  if (2.0 - 2.0 * ci <= _squaredChordLength) {
252  bmax = 1.0;
253  }
254  if (2.0 + 2.0 * ci <= _squaredChordLength) {
255  bmin = -1.0;
256  }
257  double b0 = ci * c + di * s;
258  bmax = std::max(bmax, b0 + MAX_BOUNDARY_ERROR);
259  bmin = std::min(bmin, b0 - MAX_BOUNDARY_ERROR);
260  double b1 = ci * c - di * s;
261  bmax = std::max(bmax, b1 + MAX_BOUNDARY_ERROR);
262  bmin = std::min(bmin, b1 - MAX_BOUNDARY_ERROR);
263  e[i] = Interval1d(std::max(-1.0, bmin), std::min(1.0, bmax));
264  }
265  return Box3d(e[0], e[1], e[2]);
266 }
267 
269  if (contains(v)) {
270  return CONTAINS;
271  } else if (isEmpty()) {
272  return DISJOINT | WITHIN;
273  }
274  return DISJOINT;
275 }
276 
278  // Box-Circle relations are implemented by Box.
279  return invert(b.relate(*this));
280 }
281 
283  if (isEmpty()) {
284  if (c.isEmpty()) {
285  return CONTAINS | DISJOINT | WITHIN;
286  }
287  return DISJOINT | WITHIN;
288  } else if (c.isEmpty()) {
289  return CONTAINS | DISJOINT;
290  }
291  // Neither circle is empty.
292  if (isFull()) {
293  if (c.isFull()) {
294  return CONTAINS | WITHIN;
295  }
296  return CONTAINS;
297  } else if (c.isFull()) {
298  return WITHIN;
299  }
300  // Neither circle is full.
301  NormalizedAngle cc(_center, c._center);
302  if (cc > _openingAngle + c._openingAngle + 4.0 * Angle(MAX_ASIN_ERROR)) {
303  return DISJOINT;
304  }
305  if (cc + c._openingAngle + 4.0 * Angle(MAX_ASIN_ERROR) <= _openingAngle) {
306  return CONTAINS;
307  } else if (cc + _openingAngle + 4.0 * Angle(MAX_ASIN_ERROR) <=
308  c._openingAngle) {
309  return WITHIN;
310  }
311  return INTERSECTS;
312 }
313 
315  // ConvexPolygon-Circle relations are implemented by ConvexPolygon.
316  return invert(p.relate(*this));
317 }
318 
320  // Ellipse-Circle relations are implemented by Ellipse.
321  return invert(e.relate(*this));
322 }
323 
325  std::vector<uint8_t> buffer;
326  uint8_t tc = TYPE_CODE;
327  buffer.reserve(ENCODED_SIZE);
328  buffer.push_back(tc);
329  encodeDouble(_center.x(), buffer);
330  encodeDouble(_center.y(), buffer);
331  encodeDouble(_center.z(), buffer);
332  encodeDouble(_squaredChordLength, buffer);
333  encodeDouble(_openingAngle.asRadians(), buffer);
334  return buffer;
335 }
336 
337 std::unique_ptr<Circle> Circle::decode(uint8_t const * buffer, size_t n) {
338  if (buffer == nullptr || n != ENCODED_SIZE || *buffer != TYPE_CODE) {
339  throw std::runtime_error("Byte-string is not an encoded Circle");
340  }
341  std::unique_ptr<Circle> circle(new Circle);
342  ++buffer;
343  double x = decodeDouble(buffer); buffer += 8;
344  double y = decodeDouble(buffer); buffer += 8;
345  double z = decodeDouble(buffer); buffer += 8;
346  double squaredChordLength = decodeDouble(buffer); buffer += 8;
347  double openingAngle = decodeDouble(buffer); buffer += 8;
348  circle->_center = UnitVector3d::fromNormalized(x, y, z);
349  circle->_squaredChordLength = squaredChordLength;
350  circle->_openingAngle = Angle(openingAngle);
351  return circle;
352 }
353 
355  char tail[32];
356  std::snprintf(tail, sizeof(tail), ", %.17g]}", c.getSquaredChordLength());
357  return os << "{\"Circle\": [" << c.getCenter() << tail;
358 }
359 
360 }} // namespace lsst::sphgeom
This file declares a class for representing axis-aligned bounding boxes in ℝ³.
double x
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
std::ostream * os
Definition: Schema.cc:557
int y
Definition: SpanSet.cc:48
table::Key< int > b
table::Key< int > a
T asin(T... args)
Angle represents an angle in radians.
Definition: Angle.h:43
double asRadians() const
asRadians returns the value of this angle in units of radians.
Definition: Angle.h:85
Box3d represents a box in ℝ³.
Definition: Box3d.h:42
static Box3d aroundUnitSphere()
aroundUnitSphere returns a minimal Box3d containing the unit sphere.
Definition: Box3d.h:56
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition: Box.h:54
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:43
Circle is a circular region on the unit sphere that contains its boundary.
Definition: Circle.h:46
static Angle openingAngleFor(double squaredChordLength)
openingAngleFor computes and returns the angular separation between points in S² that are separated b...
Definition: Circle.cc:52
bool isEmpty() const
Definition: Circle.h:109
Box getBoundingBox() const override
getBoundingBox returns a bounding-box for this region.
Definition: Circle.cc:211
static Circle empty()
Definition: Circle.h:50
Circle & clipTo(UnitVector3d const &x)
Definition: Circle.cc:88
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:184
static Circle full()
Definition: Circle.h:52
static std::unique_ptr< Circle > decode(std::vector< uint8_t > const &s)
Definition: Circle.h:251
bool isFull() const
Definition: Circle.h:114
std::vector< uint8_t > encode() const override
encode serializes this region into an opaque byte string.
Definition: Circle.cc:324
bool isDisjointFrom(UnitVector3d const &x) const
Definition: Circle.h:137
Relationship relate(UnitVector3d const &v) const
Definition: Circle.cc:268
double getSquaredChordLength() const
getSquaredChordLength returns the squared length of chords between the circle center and points on th...
Definition: Circle.h:123
UnitVector3d const & getCenter() const
getCenter returns the center of this circle as a unit vector.
Definition: Circle.h:118
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.
Definition: Circle.cc:219
virtual bool contains(UnitVector3d const &) const=0
contains tests whether the given unit vector is inside this region.
Circle & complement()
complement sets this circle to the closure of its complement.
Definition: Circle.cc:194
Circle & expandTo(UnitVector3d const &x)
Definition: Circle.cc:120
Circle()
This constructor creates an empty circle.
Definition: Circle.h:65
static double squaredChordLengthFor(Angle openingAngle)
squaredChordLengthFor computes and returns the squared chord length between points in S² that are sep...
Definition: Circle.cc:41
static constexpr uint8_t TYPE_CODE
Definition: Circle.h:48
ConvexPolygon is a closed convex polygon on the unit sphere.
Definition: ConvexPolygon.h:57
Relationship relate(Region const &r) const override
Ellipse is an elliptical region on the sphere.
Definition: Ellipse.h:170
Relationship relate(Region const &r) const override
Definition: Ellipse.h:289
Interval1d represents closed intervals of ℝ.
Definition: Interval1d.h:40
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition: LonLat.h:48
Angle getLat() const
Definition: LonLat.h:90
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.
Definition: UnitVector3d.h:55
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
Definition: UnitVector3d.h:82
static UnitVector3d orthogonalTo(Vector3d const &v)
orthogonalTo returns an arbitrary unit vector that is orthogonal to v.
Definition: UnitVector3d.cc:34
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this unit vector and v.
Definition: UnitVector3d.h:155
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)
lsst::geom::Angle Angle
Definition: misc.h:33
constexpr double MAX_ASIN_ERROR
Definition: constants.h:45
void encodeDouble(double item, std::vector< uint8_t > &buffer)
encode appends an IEEE double in little-endian byte order to the end of buffer.
Definition: codec.h:46
std::ostream & operator<<(std::ostream &, Angle const &)
Definition: Angle.cc:34
double sin(Angle const &a)
Definition: Angle.h:102
double decodeDouble(uint8_t const *buffer)
decode extracts an IEEE double from the 8 byte little-endian byte sequence in buffer.
Definition: codec.h:66
double cos(Angle const &a)
Definition: Angle.h:103
constexpr double PI
Definition: constants.h:36
Relationship invert(Relationship r)
Given the relationship between two sets A and B (i.e.
Definition: Relationship.h:55
A base class for image defects.
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