LSST Applications  21.0.0+3c14b91618,21.0.0+9f51b1e3f7,21.0.0-1-ga51b5d4+6691386486,21.0.0-10-g2408eff+1a328412bf,21.0.0-10-g560fb7b+52fd22d7b4,21.0.0-10-g8d1d15d+2f9043cae0,21.0.0-10-gcf60f90+d15de71c48,21.0.0-11-g25eff31+d43066e4ef,21.0.0-15-g490e301a+a676f0d5cf,21.0.0-2-g103fe59+4758c8ef83,21.0.0-2-g1367e85+a9f57e981a,21.0.0-2-g45278ab+9f51b1e3f7,21.0.0-2-g5242d73+a9f57e981a,21.0.0-2-g7f82c8f+1bcc828e4f,21.0.0-2-g8f08a60+e6fd6d9ff9,21.0.0-2-ga326454+1bcc828e4f,21.0.0-2-gde069b7+66c51b65da,21.0.0-2-gecfae73+251b9830c3,21.0.0-2-gfc62afb+a9f57e981a,21.0.0-20-g09baf175d+e1e7d1c708,21.0.0-3-g357aad2+854c3902c3,21.0.0-3-g4be5c26+a9f57e981a,21.0.0-3-g65f322c+feaa1990e9,21.0.0-3-g7d9da8d+3c14b91618,21.0.0-3-ge02ed75+bda8df9b93,21.0.0-4-g591bb35+bda8df9b93,21.0.0-4-g65b4814+52fd22d7b4,21.0.0-4-g88306b8+656365ce3f,21.0.0-4-gccdca77+86bf7a300d,21.0.0-4-ge8a399c+b99e86088e,21.0.0-5-gd00fb1e+6a0dc09319,21.0.0-51-gd3b42663+bf9f0d1c8b,21.0.0-6-g2d4f3f3+9f51b1e3f7,21.0.0-7-g04766d7+ee2ae02087,21.0.0-7-g98eecf7+3609eddee2,21.0.0-7-gde99fe0+bda8df9b93,21.0.0-8-gd70ce6f+79e45e76e4,master-gac4afde19b+bda8df9b93,w.2021.10
LSST Data Management Base Package
Box.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/Box.h"
27 
28 #include <cmath>
29 #include <ostream>
30 #include <stdexcept>
31 
32 #include "lsst/sphgeom/Box3d.h"
33 #include "lsst/sphgeom/Circle.h"
35 #include "lsst/sphgeom/Ellipse.h"
36 #include "lsst/sphgeom/codec.h"
37 #include "lsst/sphgeom/utils.h"
38 
39 
40 namespace lsst {
41 namespace sphgeom {
42 
44  if (r <= Angle(0.0)) {
45  return NormalizedAngle(0.0);
46  }
47  // If a circle centered at the given latitude contains a pole, then
48  // its bounding box contains all possible longitudes.
49  if (abs(lat) + r >= Angle(0.5 * PI)) {
50  return NormalizedAngle(PI);
51  }
52  // Now, consider the circle with opening angle r > 0 centered at (0,δ)
53  // with r < π/2 and |δ| ≠ π/2. The circle center vector in ℝ³ is
54  // c = (cos δ, 0, sin δ). Its bounding box spans longitudes [-α,α], where
55  // α is the desired half-width. The plane corresponding to longitude α has
56  // normal vector (-sin α, cos α, 0) and is tangent to the circle at point
57  // p. The great circle segment between the center of the circle and the
58  // plane normal passes through p and has arc length π/2 + r, so that
59  //
60  // (cos δ, 0, sin δ) · (-sin α, cos α, 0) = cos (π/2 + r)
61  //
62  // Solving for α gives
63  //
64  // α = arcsin (sin r / cos δ)
65  //
66  // In the actual computation, there is an absolute value and an explicit
67  // arcsin domain check to cope with rounding errors. An alternate way to
68  // compute this is:
69  //
70  // α = arctan (sin r / √(cos(δ - r) cos(δ + r)))
71  double s = std::fabs(sin(r) / cos(lat));
72  if (s >= 1.0) {
73  return NormalizedAngle(0.5 * PI);
74  }
75  return NormalizedAngle(std::asin(s));
76 }
77 
79  // The basic idea is to compute the union of the bounding boxes for all
80  // circles of opening angle r with centers inside this box.
81  //
82  // The bounding box for a circle of opening angle r with center latitude
83  // |δ| ≤ π/2 - r has height 2r.
84  //
85  // Given fixed r, the width of the bounding box for the circle centered at
86  // latitude δ grows monotonically with |δ| - for justification, see the
87  // derivation in halfWidthForCircle(). The maximum width is therefore
88  // attained when the circle is centered at one of the latitude angle
89  // boundaries of this box. If max(|δ|) ≥ π/2 - r, it is 2π.
90  //
91  // Dilating the longitude interval of this box by the maximum width and
92  // the latitude interval by r gives the desired result.
93  if (isEmpty() || isFull() || r <= Angle(0.0)) {
94  return *this;
95  }
96  Angle maxAbsLatitude = std::max(abs(_lat.getA()), abs(_lat.getB()));
97  NormalizedAngle w = halfWidthForCircle(r, maxAbsLatitude);
98  return dilateBy(w, r);
99 }
100 
102  if (isEmpty() || isFull()) {
103  return *this;
104  }
105  _lon.dilateBy(w);
106  if (!h.isNan()) {
107  Angle a = (_lat.getA() > Angle(-0.5 * PI)) ? _lat.getA() - h : _lat.getA();
108  Angle b = (_lat.getB() < Angle(0.5 * PI)) ? _lat.getB() + h : _lat.getB();
109  _lat = AngleInterval(a, b);
110  }
111  _enforceInvariants();
112  return *this;
113 }
114 
115 double Box::getArea() const {
116  if (isEmpty()) {
117  return 0.0;
118  }
119  // Given a, b ∈ [-π/2, π/2] and a std::sin implementation that is not
120  // correctly rounded, b > a does not imply that std::sin(b) > std::sin(a).
121  // To avoid potentially returning a negative area, defensively take an
122  // absolute value.
123  double dz = sin(_lat.getB()) - sin(_lat.getA());
124  return std::fabs(_lon.getSize().asRadians() * dz);
125 }
126 
128  if (isEmpty()) {
129  return Box3d();
130  }
131  if (isFull()) {
132  return Box3d::aroundUnitSphere();
133  }
134  double slata = sin(_lat.getA()), clata = cos(_lat.getA());
135  double slatb = sin(_lat.getB()), clatb = cos(_lat.getB());
136  double slona = sin(_lon.getA()), clona = cos(_lon.getA());
137  double slonb = sin(_lon.getB()), clonb = cos(_lon.getB());
138  // Compute the minimum/maximum x/y values of the box vertices.
139  double xmin = std::min(std::min(clona * clata, clonb * clata),
140  std::min(clona * clatb, clonb * clatb)) - 2.5 * EPSILON;
141  double xmax = std::max(std::max(clona * clata, clonb * clata),
142  std::max(clona * clatb, clonb * clatb)) + 2.5 * EPSILON;
143  double ymin = std::min(std::min(slona * clata, slonb * clata),
144  std::min(slona * clatb, slonb * clatb)) - 2.5 * EPSILON;
145  double ymax = std::max(std::max(slona * clata, slonb * clata),
146  std::max(slona * clatb, slonb * clatb)) + 2.5 * EPSILON;
147  // Compute the maximum latitude cosine of points in the box.
148  double mlc;
149  if (_lat.contains(Angle(0.0))) {
150  mlc = 1.0;
151  // The box intersects the equator - the x or y extrema of the box may be
152  // at the intersection of the box edge meridians with the equator.
153  xmin = std::min(xmin, std::min(clona, clonb) - EPSILON);
154  xmax = std::max(xmax, std::max(clona, clonb) + EPSILON);
155  ymin = std::min(ymin, std::min(slona, slonb) - EPSILON);
156  ymax = std::max(ymax, std::max(slona, slonb) + EPSILON);
157  } else {
158  // Note that clata and clatb are positive.
159  mlc = std::max(clata, clatb) + EPSILON;
160  }
161  // Check for extrema on the box edges parallel to the equator.
162  if (_lon.contains(NormalizedAngle(0.0))) {
163  xmax = std::max(xmax, mlc);
164  }
165  if (_lon.contains(NormalizedAngle(0.5 * PI))) {
166  ymax = std::max(ymax, mlc);
167  }
168  if (_lon.contains(NormalizedAngle(PI))) {
169  xmin = std::min(xmin, -mlc);
170  }
171  if (_lon.contains(NormalizedAngle(1.5 * PI))) {
172  ymin = std::min(ymin, -mlc);
173  }
174  // Clamp x/y extrema to [-1, 1]
175  xmin = std::max(-1.0, xmin);
176  xmax = std::min(1.0, xmax);
177  ymin = std::max(-1.0, ymin);
178  ymax = std::min(1.0, ymax);
179  // Compute z extrema.
180  double zmin = std::max(-1.0, slata - EPSILON);
181  double zmax = std::min(1.0, slatb + EPSILON);
182  return Box3d(Interval1d(xmin, xmax),
183  Interval1d(ymin, ymax),
184  Interval1d(zmin, zmax));
185 }
186 
188  if (isEmpty()) {
189  return Circle::empty();
190  }
191  if (isFull()) {
192  return Circle::full();
193  }
195  // The minimal bounding circle center p lies on the meridian bisecting
196  // this box. Let δ₁ and δ₂ be the minimum and maximum box latitudes.
197  if (w.asRadians() <= PI) {
198  UnitVector3d p;
199  UnitVector3d boxVerts[4] = {
200  UnitVector3d(_lon.getA(), _lat.getA()),
201  UnitVector3d(_lon.getA(), _lat.getB()),
202  UnitVector3d(_lon.getB(), _lat.getA()),
203  UnitVector3d(_lon.getB(), _lat.getB())
204  };
205  // We take advantage of rotational symmetry to fix the bisecting
206  // meridian at a longitude of zero. The box vertices then have
207  // coordinates (±w/2, δ₁), (±w/2, δ₂), and p = (0, ϕ). Converting
208  // to Cartesian coordinates gives p = (cos ϕ, 0, sin ϕ), and box
209  // vertices at (cos w/2 cos δ₁, ±sin w/2 cos δ₁, sin δ₁) and
210  // (cos w/2 cos δ₂, ±sin w/2 cos δ₂, sin δ₂).
211  //
212  // The point p₁ on the meridian that has minimum angular separation
213  // to the vertices with latitude δ₁ lies on the plane they define.
214  // The sum of the two vertex vectors is on that plane and on the plane
215  // containing the meridian. Normalizing to obtain p₁, we have
216  //
217  // (cos ϕ₁, 0, sin ϕ₁) =
218  // λ ((cos w/2 cos δ₁, sin w/2 cos δ₁, sin δ₁) +
219  // (cos w/2 cos δ₁, -sin w/2 cos δ₁, sin δ₁))
220  //
221  // for some scaling factor λ. Simplifying, we get:
222  //
223  // cos ϕ₁ = λ cos w/2 cos δ₁
224  // sin ϕ₁ = λ sin δ₁
225  //
226  // so that
227  //
228  // tan ϕ₁ = sec w/2 tan δ₁
229  //
230  // Similarly, the point p₂ on the meridian that has minimum angular
231  // separation to the vertices with latitude δ₂ satisfies:
232  //
233  // tan ϕ₂ = sec w/2 tan δ₂
234  //
235  // where ϕ₁ ≤ ϕ₂ (since δ₁ ≤ δ₂). Finally, consider the point p₃
236  // separated from each box vertex by the same angle. The dot
237  // products of p₃ with the box vertices are all identical, so
238  //
239  // cos ϕ₃ cos w/2 cos δ₁ + sin ϕ₃ sin δ₁ =
240  // cos ϕ₃ cos w/2 cos δ₂ + sin ϕ₃ sin δ₂
241  //
242  // Rearranging gives:
243  //
244  // tan ϕ₃ = - cos w/2 (cos δ₁ - cos δ₂)/(sin δ₁ - sin δ₂)
245  //
246  // which can be simplified further using a tangent half-angle identity,
247  // yielding:
248  //
249  // tan ϕ₃ = cos w/2 tan (δ₁ + δ₂)/2
250  //
251  // Consider now the function f₁(ϕ) that gives the angular separation
252  // between p with latitude ϕ and the vertices at latitude δ₁. It has
253  // a line of symmetry at ϕ = ϕ₁, and increases monotonically with
254  // |ϕ - ϕ₁|. Similarly, f₂(ϕ) has a minimum at ϕ₂ and increases
255  // monotonically with |ϕ - ϕ₂|. The two functions cross at ϕ₃. The
256  // opening angle of the bounding circle centered at latitude ϕ is
257  // given by g = max(f₁, f₂), which we seek to minimize.
258  //
259  // If ϕ₁ ≤ ϕ₃ ≤ ϕ₂, then g is minimized at ϕ = ϕ₃. Otherwise, it
260  // is minimized at either ϕ₁ or ϕ₂.
261  double phi1, phi2, phi3;
262  double c = cos(0.5 * w);
263  if (c == 0.0) {
264  // This code should never execute. If it does, the implementation
265  // of std::cos is broken.
266  phi1 = ::copysign(0.5 * PI, _lat.getA().asRadians());
267  phi2 = ::copysign(0.5 * PI, _lat.getB().asRadians());
268  phi3 = 0.0;
269  } else {
270  phi1 = std::atan(tan(_lat.getA()) / c);
271  phi2 = std::atan(tan(_lat.getB()) / c);
272  phi3 = std::atan(c * tan(_lat.getCenter()));
273  }
274  if (phi1 <= phi3 && phi3 <= phi2) {
275  p = UnitVector3d(_lon.getCenter(), Angle(phi3));
276  } else {
277  UnitVector3d p1 = UnitVector3d(_lon.getCenter(), Angle(phi1));
278  UnitVector3d p2 = UnitVector3d(_lon.getCenter(), Angle(phi2));
279  if (p1.dot(boxVerts[0]) > p2.dot(boxVerts[1])) {
280  p = p2;
281  } else {
282  p = p1;
283  }
284  }
285  // Compute the maximum squared chord length between p and the box
286  // vertices, so that each one is guaranteed to lie in the bounding
287  // circle, regardless of numerical error in the above.
288  double cl2 = (p - boxVerts[0]).getSquaredNorm();
289  for (int i = 1; i < 4; ++i) {
290  cl2 = std::max(cl2, (p - boxVerts[i]).getSquaredNorm());
291  }
292  // Add double the maximum squared-chord-length error, so that the
293  // bounding circle we return also reliably CONTAINS this box.
294  return Circle(p, cl2 + 2.0 * MAX_SQUARED_CHORD_LENGTH_ERROR);
295  }
296  // The box spans more than π radians in longitude. First, pick the smaller
297  // of the bounding circles centered at the north and south pole.
298  Angle r;
299  UnitVector3d v;
300  if (abs(_lat.getA()) <= abs(_lat.getB())) {
301  v = UnitVector3d::Z();
302  r = Angle(0.5 * PI) - _lat.getA();
303  } else {
304  v = -UnitVector3d::Z();
305  r = _lat.getB() + Angle(0.5 * PI);
306  }
307  // If the box does not span all longitude angles, we also consider the
308  // equatorial bounding circle with center longitude equal to the longitude
309  // of the box center. The smaller of the polar and equatorial bounding
310  // circles is returned.
311  if (!_lon.isFull() && 0.5 * w < r) {
312  r = 0.5 * w;
313  v = UnitVector3d(_lon.getCenter(), Angle(0.0));
314  }
315  return Circle(v, r + 4.0 * Angle(MAX_ASIN_ERROR));
316 }
317 
318 Relationship Box::relate(Circle const & c) const {
319  if (isEmpty()) {
320  if (c.isEmpty()) {
321  return CONTAINS | DISJOINT | WITHIN;
322  }
323  return DISJOINT | WITHIN;
324  } else if (c.isEmpty()) {
325  return CONTAINS | DISJOINT;
326  }
327  if (isFull()) {
328  if (c.isFull()) {
329  return CONTAINS | WITHIN;
330  }
331  return CONTAINS;
332  } else if (c.isFull()) {
333  return WITHIN;
334  }
335  // Neither region is empty or full. We now determine whether or not the
336  // circle and box boundaries intersect.
337  //
338  // If the box vertices are not all inside or all outside of c, then the
339  // boundaries cross.
340  LonLat vertLonLat[4] = {
341  LonLat(_lon.getA(), _lat.getA()),
342  LonLat(_lon.getA(), _lat.getB()),
343  LonLat(_lon.getB(), _lat.getA()),
344  LonLat(_lon.getB(), _lat.getB())
345  };
346  UnitVector3d verts[4];
347  bool inside = false;
348  for (int i = 0; i < 4; ++i) {
349  verts[i] = UnitVector3d(vertLonLat[i]);
350  double d = (verts[i] - c.getCenter()).getSquaredNorm();
351  if (std::fabs(d - c.getSquaredChordLength()) <
353  // A box vertex is close to the circle boundary.
354  return INTERSECTS;
355  }
356  bool b = d < c.getSquaredChordLength();
357  if (i == 0) {
358  inside = b;
359  } else if (inside != b) {
360  // There are box vertices both inside and outside of c.
361  return INTERSECTS;
362  }
363  }
364  UnitVector3d norms[2] = {
367  };
368  if (inside) {
369  // All box vertices are inside c. Look for points in the box edge
370  // interiors that are outside c.
371  for (int i = 0; i < 2; ++i) {
372  double d = getMaxSquaredChordLength(
373  c.getCenter(), verts[2 * i + 1], verts[2 * i], norms[i]);
374  if (d > c.getSquaredChordLength() -
376  return INTERSECTS;
377  }
378  }
379  LonLat cc(-c.getCenter());
380  if (_lon.contains(cc.getLon())) {
381  // The points furthest from the center of c on the small circles
382  // defined by the box edges with constant latitude are in the box
383  // edge interiors. Find the largest squared chord length to either.
384  Angle a = std::min(getMinAngleToCircle(cc.getLat(), _lat.getA()),
385  getMinAngleToCircle(cc.getLat(), _lat.getB()));
386  double d = Circle::squaredChordLengthFor(Angle(PI) - a);
387  if (d > c.getSquaredChordLength() -
389  return INTERSECTS;
390  }
391  }
392  // The box boundary is completely inside c. However, the box is not
393  // necessarily within c: consider a circle with opening angle equal to
394  // π - ε. If a box contains the complement of such a circle, then
395  // intersecting it with that circle will punch a hole in the box. In
396  // this case each region contains the boundary of the other, but
397  // neither region contains the other.
398  //
399  // To handle this case, check that the box does not contain the
400  // complement of c - since the boundaries do not intersect, this is the
401  // case iff the box contains the center of the complement of c.
402  if (contains(cc)) {
403  return INTERSECTS;
404  }
405  return WITHIN;
406  }
407  // All box vertices are outside c. Look for points in the box edge
408  // interiors that are inside c.
409  for (int i = 0; i < 2; ++i) {
410  double d = getMinSquaredChordLength(
411  c.getCenter(), verts[2 * i + 1], verts[2 * i], norms[i]);
413  return INTERSECTS;
414  }
415  }
416  LonLat cc(c.getCenter());
417  if (_lon.contains(cc.getLon())) {
418  // The points closest to the center of c on the small circles
419  // defined by the box edges with constant latitude are in the box
420  // edge interiors. Find the smallest squared chord length to either.
421  Angle a = std::min(getMinAngleToCircle(cc.getLat(), _lat.getA()),
422  getMinAngleToCircle(cc.getLat(), _lat.getB()));
423  double d = Circle::squaredChordLengthFor(a);
425  return INTERSECTS;
426  }
427  }
428  // The box boundary is completely outside of c. If the box contains the
429  // circle center, then the box contains c. Otherwise, the box and circle
430  // are disjoint.
431  if (contains(cc)) {
432  return CONTAINS;
433  }
434  return DISJOINT;
435 }
436 
438  // ConvexPolygon-Box relations are implemented by ConvexPolygon.
439  return invert(p.relate(*this));
440 }
441 
442 Relationship Box::relate(Ellipse const & e) const {
443  // Ellipse-Box relations are implemented by Ellipse.
444  return invert(e.relate(*this));
445 }
446 
448  std::vector<uint8_t> buffer;
449  uint8_t tc = TYPE_CODE;
450  buffer.reserve(ENCODED_SIZE);
451  buffer.push_back(tc);
452  encodeDouble(_lon.getA().asRadians(), buffer);
453  encodeDouble(_lon.getB().asRadians(), buffer);
454  encodeDouble(_lat.getA().asRadians(), buffer);
455  encodeDouble(_lat.getB().asRadians(), buffer);
456  return buffer;
457 }
458 
459 std::unique_ptr<Box> Box::decode(uint8_t const * buffer, size_t n) {
460  if (buffer == nullptr || n != ENCODED_SIZE || *buffer != TYPE_CODE) {
461  throw std::runtime_error("Byte-string is not an encoded Box");
462  }
463  std::unique_ptr<Box> box(new Box);
464  ++buffer;
465  double a = decodeDouble(buffer); buffer += 8;
466  double b = decodeDouble(buffer); buffer += 8;
468  a = decodeDouble(buffer); buffer += 8;
469  b = decodeDouble(buffer); buffer += 8;
470  box->_lat = AngleInterval::fromRadians(a, b);
471  box->_enforceInvariants();
472  return box;
473 }
474 
476  return os << "{\"Box\": [" << b.getLon() << ", " << b.getLat() << "]}";
477 }
478 
479 }} // 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...
std::ostream * os
Definition: Schema.cc:746
int xmax
Definition: SpanSet.cc:49
int xmin
Definition: SpanSet.cc:49
T asin(T... args)
T atan(T... args)
Angle represents an angle in radians.
Definition: Angle.h:43
bool isNan() const
isNan returns true if the angle value is NaN.
Definition: Angle.h:91
double asRadians() const
asRadians returns the value of this angle in units of radians.
Definition: Angle.h:85
AngleInterval represents closed intervals of arbitrary angles.
Definition: AngleInterval.h:39
static AngleInterval fromRadians(double x, double y)
Definition: AngleInterval.h:48
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
std::vector< uint8_t > encode() const override
encode serializes this region into an opaque byte string.
Definition: Box.cc:447
NormalizedAngle getWidth() const
getWidth returns the width in longitude angle of this box.
Definition: Box.h:165
double getArea() const
getArea returns the area of this box in steradians.
Definition: Box.cc:115
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.
Definition: Box.cc:127
Relationship relate(LonLat const &p) const
Definition: Box.h:297
static std::unique_ptr< Box > decode(std::vector< uint8_t > const &s)
Definition: Box.h:338
bool isFull() const
isFull returns true if this box contains all points on the unit sphere.
Definition: Box.h:155
static constexpr uint8_t TYPE_CODE
Definition: Box.h:56
Box & dilateBy(Angle r)
dilateBy minimally expands this Box to include all points within angular separation r of its boundary...
Definition: Box.cc:78
bool contains(LonLat const &x) const
Definition: Box.h:174
bool isEmpty() const
isEmpty returns true if this box does not contain any points.
Definition: Box.h:151
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
Definition: Box.cc:187
Circle is a circular region on the unit sphere that contains its boundary.
Definition: Circle.h:46
bool isEmpty() const
Definition: Circle.h:109
static Circle empty()
Definition: Circle.h:50
static Circle full()
Definition: Circle.h:52
bool isFull() const
Definition: Circle.h:114
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
static double squaredChordLengthFor(Angle openingAngle)
squaredChordLengthFor computes and returns the squared chord length between points in S² that are sep...
Definition: Circle.cc:41
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:169
Relationship relate(Region const &r) const override
Definition: Ellipse.h:286
Interval1d represents closed intervals of ℝ.
Definition: Interval1d.h:39
Scalar getCenter() const
getCenter returns the center of this interval.
Definition: Interval.h:89
Scalar getA() const
getA returns the lower endpoint of this interval.
Definition: Interval.h:76
bool contains(Scalar x) const
Definition: Interval.h:98
Scalar getB() const
getB returns the upper endpoint of this interval.
Definition: Interval.h:80
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition: LonLat.h:48
Angle getLat() const
Definition: LonLat.h:90
NormalizedAngle getLon() const
Definition: LonLat.h:88
NormalizedAngle is an angle that lies in the range [0, 2π), with one exception - a NormalizedAngle ca...
double asRadians() const
asRadians returns the value of this angle in units of radians.
bool isFull() const
isFull returns true if this interval contains all normalized angles.
NormalizedAngle getA() const
getA returns the first endpoint of this interval.
NormalizedAngle getCenter() const
getCenter returns the center of this interval.
NormalizedAngle getSize() const
getSize returns the size (length, width) of this interval.
NormalizedAngle getB() const
getB returns the second endpoint of this interval.
NormalizedAngleInterval & dilateBy(Angle x)
static NormalizedAngleInterval fromRadians(double a, double b)
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
static UnitVector3d Z()
Definition: UnitVector3d.h:101
double dot(Vector3d const &v) const
dot returns the inner product of this unit vector and v.
Definition: UnitVector3d.h:152
static UnitVector3d orthogonalTo(Vector3d const &v)
orthogonalTo returns an arbitrary unit vector that is orthogonal to v.
Definition: UnitVector3d.cc:34
This file contains simple helper functions for encoding and decoding primitive types to/from byte str...
T fabs(T... args)
T max(T... args)
T min(T... args)
lsst::geom::Angle Angle
Definition: misc.h:33
Angle getMinAngleToCircle(Angle x, Angle c)
getMinAngleToCircle returns the minimum angular separation between a point at latitude x and the poin...
Definition: utils.h:62
constexpr double MAX_ASIN_ERROR
Definition: constants.h:45
Angle abs(Angle const &a)
Definition: Angle.h:106
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:38
std::ostream & operator<<(std::ostream &, Angle const &)
Definition: Angle.cc:34
double getMaxSquaredChordLength(Vector3d const &v, Vector3d const &a, Vector3d const &b, Vector3d const &n)
Let p be the unit vector furthest from v that lies on the plane with normal n in the direction of the...
Definition: utils.cc:58
double getMinSquaredChordLength(Vector3d const &v, Vector3d const &a, Vector3d const &b, Vector3d const &n)
Let p be the unit vector closest to v that lies on the plane with normal n in the direction of the cr...
Definition: utils.cc:36
double sin(Angle const &a)
Definition: Angle.h:102
constexpr double MAX_SQUARED_CHORD_LENGTH_ERROR
Definition: constants.h:50
double tan(Angle const &a)
Definition: Angle.h:104
constexpr double EPSILON
Definition: constants.h:54
double decodeDouble(uint8_t const *buffer)
decode extracts an IEEE double from the 8 byte little-endian byte sequence in buffer.
Definition: codec.h:59
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.
This file declares miscellaneous utility functions.
double w
Definition: CoaddPsf.cc:69
table::Key< int > b
table::Key< int > a