LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
Classes | Functions
lsst::sphgeom::detail Namespace Reference

Classes

class  PixelFinder
 

Functions

template<typename VertexIterator >
UnitVector3d centroid (VertexIterator const begin, VertexIterator const end)
 
template<typename VertexIterator >
Circle boundingCircle (VertexIterator const begin, VertexIterator const end)
 
template<typename VertexIterator >
Box boundingBox (VertexIterator const begin, VertexIterator const end)
 
template<typename VertexIterator >
Box3d boundingBox3d (VertexIterator const begin, VertexIterator const end)
 
template<typename VertexIterator >
bool contains (VertexIterator const begin, VertexIterator const end, UnitVector3d const &v)
 
template<typename VertexIterator >
Relationship relate (VertexIterator const begin, VertexIterator const end, Box const &b)
 
template<typename VertexIterator >
Relationship relate (VertexIterator const begin, VertexIterator const end, Circle const &c)
 
template<typename VertexIterator1 , typename VertexIterator2 >
Relationship relate (VertexIterator1 const begin1, VertexIterator1 const end1, VertexIterator2 const begin2, VertexIterator2 const end2)
 
template<typename VertexIterator >
Relationship relate (VertexIterator const begin, VertexIterator const end, ConvexPolygon const &p)
 
template<typename VertexIterator >
Relationship relate (VertexIterator const begin, VertexIterator const end, Ellipse const &e)
 
template<template< typename, bool > class Finder, bool InteriorOnly>
RangeSet findPixels (Region const &r, size_t maxRanges, int level)
 

Function Documentation

◆ boundingBox()

template<typename VertexIterator >
Box lsst::sphgeom::detail::boundingBox ( VertexIterator const  begin,
VertexIterator const  end 
)

Definition at line 85 of file ConvexPolygonImpl.h.

85  {
86  Angle const eps(5.0e-10); // ~ 0.1 milli-arcseconds
87  Box bbox;
88  VertexIterator i = std::prev(end);
89  VertexIterator j = begin;
90  bool haveCW = false;
91  bool haveCCW = false;
92  // Compute the bounding box for each vertex. When converting a Vector3d
93  // to a LonLat, the relative error on the longitude is about 4*2^-53,
94  // and the relative error on the latitude is about twice that (assuming
95  // std::atan2 and std::sqrt accurate to within 1 ulp). We convert each
96  // vertex to a conservative bounding box for its spherical coordinates,
97  // and compute a bounding box for the union of all these boxes.
98  //
99  // Furthermore, the latitude range of an edge can be greater than the
100  // latitude range of its endpoints - this occurs when the minimum or
101  // maximum latitude point on the great circle defined by the edge vertices
102  // lies in the edge interior.
103  for (; j != end; i = j, ++j) {
104  LonLat p(*j);
105  bbox.expandTo(Box(p, eps, eps));
106  if (!haveCW || !haveCCW) {
107  int o = orientationZ(*i, *j);
108  haveCCW = haveCCW || (o > 0);
109  haveCW = haveCW || (o < 0);
110  }
111  // Compute the plane normal for edge i, j.
112  Vector3d n = (*i).robustCross(*j);
113  // Compute a vector v with positive z component that lies on both the
114  // edge plane and on the plane defined by the z axis and the edge plane
115  // normal. This is the direction of maximum latitude for the great
116  // circle containing the edge, and -v is the direction of minimum
117  // latitude.
118  //
119  // TODO(smm): Do a proper error analysis.
120  Vector3d v(-n.x() * n.z(),
121  -n.y() * n.z(),
122  n.x() * n.x() + n.y() * n.y());
123  if (v != Vector3d()) {
124  // The plane defined by the z axis and n has normal
125  // (-n.y(), n.x(), 0.0). Compute the dot product of this plane
126  // normal with vertices i and j.
127  double zni = i->y() * n.x() - i->x() * n.y();
128  double znj = j->y() * n.x() - j->x() * n.y();
129  // Check if v or -v is in the edge interior.
130  if (zni > 0.0 && znj < 0.0) {
131  bbox = Box(bbox.getLon(), bbox.getLat().expandedTo(
132  LonLat::latitudeOf(v) + eps));
133  } else if (zni < 0.0 && znj > 0.0) {
134  bbox = Box(bbox.getLon(), bbox.getLat().expandedTo(
135  LonLat::latitudeOf(-v) - eps));
136  }
137  }
138  }
139  // If this polygon contains a pole, its bounding box must contain all
140  // longitudes.
141  if (!haveCW) {
142  Box northPole(Box::allLongitudes(), AngleInterval(Angle(0.5 * PI)));
143  bbox.expandTo(northPole);
144  } else if (!haveCCW) {
145  Box southPole(Box::allLongitudes(), AngleInterval(Angle(-0.5 * PI)));
146  bbox.expandTo(southPole);
147  }
148  return bbox;
149 }
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
int end
T begin(T... args)
lsst::geom::Angle Angle
Definition: misc.h:33
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
int orientationZ(UnitVector3d const &b, UnitVector3d const &c)
orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).
Definition: orientation.cc:237
T prev(T... args)

◆ boundingBox3d()

template<typename VertexIterator >
Box3d lsst::sphgeom::detail::boundingBox3d ( VertexIterator const  begin,
VertexIterator const  end 
)

Definition at line 152 of file ConvexPolygonImpl.h.

152  {
153  static double const maxError = 1.0e-14;
154  // Compute the extrema of all vertex coordinates.
155  VertexIterator j = begin;
156  double emin[3] = { j->x(), j->y(), j->z() };
157  double emax[3] = { j->x(), j->y(), j->z() };
158  for (++j; j != end; ++j) {
159  for (int i = 0; i < 3; ++i) {
160  double v = j->operator()(i);
161  emin[i] = std::min(emin[i], v);
162  emax[i] = std::max(emax[i], v);
163  }
164  }
165  // Compute the extrema of all edges.
166  //
167  // It can be shown that the great circle with unit normal vector
168  // n = (n₀, n₁, n₂) has extrema in x at:
169  //
170  // (∓√(1 - n₀²), ±n₁n₀/√(1 - n₀²), ±n₂n₀/√(1 - n₀²))
171  //
172  // in y at:
173  //
174  // (±n₀n₁/√(1 - n₁²), ∓√(1 - n₁²), ±n₂n₁/√(1 - n₁²))
175  //
176  // and in z at
177  //
178  // (±n₀n₂/√(1 - n₂²), ±n₁n₂/√(1 - n₂²), ∓√(1 - n₂²))
179  //
180  // Compute these vectors for each edge, determine whether they lie in
181  // the edge, and update the extrema if so. Rounding errors in these
182  // computations are compensated for by expanding the bounding box
183  // prior to returning it.
184  j = std::prev(end);
185  VertexIterator k = begin;
186  for (; k != end; j = k, ++k) {
187  UnitVector3d n(j->robustCross(*k));
188  for (int i = 0; i < 3; ++i) {
189  double ni = n(i);
190  double d = std::fabs(1.0 - ni * ni);
191  if (d > 0.0) {
192  Vector3d e(i == 0 ? -d : n.x() * ni,
193  i == 1 ? -d : n.y() * ni,
194  i == 2 ? -d : n.z() * ni);
195  // If e or -e lies in the lune defined by the half great
196  // circle passing through n and a and the half great circle
197  // passing through n and b, the edge contains an extremum.
198  Vector3d v = e.cross(n);
199  double vdj = v.dot(*j);
200  double vdk = v.dot(*k);
201  if (vdj >= 0.0 && vdk <= 0.0) {
202  emin[i] = std::min(emin[i], -std::sqrt(d));
203  }
204  if (vdj <= 0.0 && vdk >= 0.0) {
205  emax[i] = std::max(emax[i], std::sqrt(d));
206  }
207  }
208  }
209  }
210  // Check whether the standard basis vectors and their antipodes
211  // are inside this polygon.
212  bool a[3] = { true, true, true };
213  bool b[3] = { true, true, true };
214  j = std::prev(end);
215  k = begin;
216  for (; k != end; j = k, ++k) {
217  // Test the standard basis vectors against the plane defined by
218  // vertices (j, k). Note that orientation(-x, *j, *k) =
219  // -orientation(x, *j, *k).
220  int ox = orientationX(*j, *k);
221  a[0] = a[0] && (ox <= 0);
222  b[0] = b[0] && (ox >= 0);
223  int oy = orientationY(*j, *k);
224  a[1] = a[1] && (oy <= 0);
225  b[1] = b[1] && (oy >= 0);
226  int oz = orientationZ(*j, *k);
227  a[2] = a[2] && (oz <= 0);
228  b[2] = b[2] && (oz >= 0);
229  }
230  // At this point, b[i] is true iff the standard basis vector eᵢ
231  // is inside all the half spaces defined by the polygon edges.
232  // Similarly, a[i] is true iff -eᵢ is inside the same half spaces.
233  for (int i = 0; i < 3; ++i) {
234  emin[i] = a[i] ? -1.0 : std::max(-1.0, emin[i] - maxError);
235  emax[i] = b[i] ? 1.0 : std::min(1.0, emax[i] + maxError);
236  }
237  return Box3d(Interval1d(emin[0], emax[0]),
238  Interval1d(emin[1], emax[1]),
239  Interval1d(emin[2], emax[2]));
240 }
table::Key< int > b
table::Key< int > a
T fabs(T... args)
T max(T... args)
T min(T... args)
int orientationX(UnitVector3d const &b, UnitVector3d const &c)
orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).
Definition: orientation.cc:227
int orientationY(UnitVector3d const &b, UnitVector3d const &c)
orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).
Definition: orientation.cc:232
T sqrt(T... args)

◆ boundingCircle()

template<typename VertexIterator >
Circle lsst::sphgeom::detail::boundingCircle ( VertexIterator const  begin,
VertexIterator const  end 
)

Definition at line 69 of file ConvexPolygonImpl.h.

69  {
70  UnitVector3d c = centroid(begin, end);
71  // Compute the maximum squared chord length between the centroid and
72  // all vertices.
73  VertexIterator i = begin;
74  double cl2 = 0.0;
75  for (; i != end; ++i) {
76  cl2 = std::max(cl2, (*i - c).getSquaredNorm());
77  }
78  // Add double the maximum squared-chord-length error, so that the
79  // bounding circle we return also reliably CONTAINS this polygon.
80  return Circle(c, cl2 + 2.0 * MAX_SQUARED_CHORD_LENGTH_ERROR);
81 }
UnitVector3d centroid(VertexIterator const begin, VertexIterator const end)
constexpr double MAX_SQUARED_CHORD_LENGTH_ERROR
Definition: constants.h:50

◆ centroid()

template<typename VertexIterator >
UnitVector3d lsst::sphgeom::detail::centroid ( VertexIterator const  begin,
VertexIterator const  end 
)

Definition at line 48 of file ConvexPolygonImpl.h.

48  {
49  // The center of mass is obtained via trivial generalization of
50  // the formula for spherical triangles from:
51  //
52  // The centroid and inertia tensor for a spherical triangle
53  // John E. Brock
54  // 1974, Naval Postgraduate School, Monterey Calif.
55  Vector3d cm;
56  VertexIterator i = std::prev(end);
57  VertexIterator j = begin;
58  for (; j != end; i = j, ++j) {
59  Vector3d v = (*i).robustCross(*j);
60  double s = 0.5 * v.normalize();
61  double c = (*i).dot(*j);
62  double a = (s == 0.0 && c == 0.0) ? 0.0 : std::atan2(s, c);
63  cm += v * a;
64  }
65  return UnitVector3d(cm);
66 }
T atan2(T... args)

◆ contains()

template<typename VertexIterator >
bool lsst::sphgeom::detail::contains ( VertexIterator const  begin,
VertexIterator const  end,
UnitVector3d const &  v 
)

Definition at line 243 of file ConvexPolygonImpl.h.

246 {
247  VertexIterator i = std::prev(end);
248  VertexIterator j = begin;
249  for (; j != end; i = j, ++j) {
250  if (orientation(v, *i, *j) < 0) {
251  return false;
252  }
253  }
254  return true;
255 }
int orientation(UnitVector3d const &a, UnitVector3d const &b, UnitVector3d const &c)
orientation computes and returns the orientations of 3 unit vectors a, b and c.
Definition: orientation.cc:135

◆ findPixels()

template<template< typename, bool > class Finder, bool InteriorOnly>
RangeSet lsst::sphgeom::detail::findPixels ( Region const &  r,
size_t  maxRanges,
int  level 
)

Definition at line 150 of file PixelFinder.h.

150  {
151  RangeSet s;
152  Circle const * c = nullptr;
153  Ellipse const * e = nullptr;
154  Box const * b = nullptr;
155  if ((c = dynamic_cast<Circle const *>(&r))) {
156  Finder<Circle, InteriorOnly> find(s, *c, level, maxRanges);
157  find();
158  } else if ((e = dynamic_cast<Ellipse const *>(&r))) {
159  Finder<Circle, InteriorOnly> find(
160  s, e->getBoundingCircle(), level, maxRanges);
161  find();
162  } else if ((b = dynamic_cast<Box const *>(&r))) {
163  Finder<Box, InteriorOnly> find(s, *b, level, maxRanges);
164  find();
165  } else {
166  Finder<ConvexPolygon, InteriorOnly> find(
167  s, dynamic_cast<ConvexPolygon const &>(r), level, maxRanges);
168  find();
169  }
170  return s;
171 }
T find(T... args)

◆ relate() [1/5]

template<typename VertexIterator >
Relationship lsst::sphgeom::detail::relate ( VertexIterator const  begin,
VertexIterator const  end,
Box const &  b 
)

Definition at line 258 of file ConvexPolygonImpl.h.

261 {
262  // TODO(smm): be more accurate when computing box relations.
263  return boundingBox(begin, end).relate(b) & (DISJOINT | WITHIN);
264 }
Relationship relate(LonLat const &p) const
Definition: Box.h:297
Box boundingBox(VertexIterator const begin, VertexIterator const end)

◆ relate() [2/5]

template<typename VertexIterator >
Relationship lsst::sphgeom::detail::relate ( VertexIterator const  begin,
VertexIterator const  end,
Circle const &  c 
)

Definition at line 267 of file ConvexPolygonImpl.h.

270 {
271  if (c.isEmpty()) {
272  return CONTAINS | DISJOINT;
273  }
274  if (c.isFull()) {
275  return WITHIN;
276  }
277  // Determine whether or not the circle and polygon boundaries intersect.
278  // If the polygon vertices are not all inside or all outside of c, then the
279  // boundaries cross.
280  bool inside = false;
281  for (VertexIterator v = begin; v != end; ++v) {
282  double d = (*v - c.getCenter()).getSquaredNorm();
283  if (std::fabs(d - c.getSquaredChordLength()) <
285  // A polygon vertex is close to the circle boundary.
286  return INTERSECTS;
287  }
288  bool b = d < c.getSquaredChordLength();
289  if (v == begin) {
290  inside = b;
291  } else if (inside != b) {
292  // There are box vertices both inside and outside of c.
293  return INTERSECTS;
294  }
295  }
296  if (inside) {
297  // All polygon vertices are inside c. Look for points in the polygon
298  // edge interiors that are outside c.
299  for (VertexIterator a = std::prev(end), b = begin; b != end; a = b, ++b) {
300  Vector3d n = a->robustCross(*b);
301  double d = getMaxSquaredChordLength(c.getCenter(), *a, *b, n);
302  if (d > c.getSquaredChordLength() -
304  return INTERSECTS;
305  }
306  }
307  // The polygon boundary is conclusively inside c. It may still be the
308  // case that the circle punches a hole in the polygon. We check that
309  // the polygon does not contain the complement of c by testing whether
310  // or not it contains the anti-center of c.
311  if (contains(begin, end, -c.getCenter())) {
312  return INTERSECTS;
313  }
314  return WITHIN;
315  }
316  // All polygon vertices are outside c. Look for points in the polygon edge
317  // interiors that are inside c.
318  for (VertexIterator a = std::prev(end), b = begin; b != end; a = b, ++b) {
319  Vector3d n = a->robustCross(*b);
320  double d = getMinSquaredChordLength(c.getCenter(), *a, *b, n);
321  if (d < c.getSquaredChordLength() + MAX_SQUARED_CHORD_LENGTH_ERROR) {
322  return INTERSECTS;
323  }
324  }
325  // The polygon boundary is conclusively outside of c. If the polygon
326  // contains the circle center, then the polygon contains c. Otherwise, the
327  // polygon and circle are disjoint.
328  if (contains(begin, end, c.getCenter())) {
329  return CONTAINS;
330  }
331  return DISJOINT;
332 }
bool contains(VertexIterator const begin, VertexIterator const end, UnitVector3d const &v)
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

◆ relate() [3/5]

template<typename VertexIterator >
Relationship lsst::sphgeom::detail::relate ( VertexIterator const  begin,
VertexIterator const  end,
ConvexPolygon const &  p 
)

Definition at line 394 of file ConvexPolygonImpl.h.

397 {
398  return relate(begin, end, p.getVertices().begin(), p.getVertices().end());
399 }
Relationship relate(VertexIterator const begin, VertexIterator const end, Ellipse const &e)

◆ relate() [4/5]

template<typename VertexIterator >
Relationship lsst::sphgeom::detail::relate ( VertexIterator const  begin,
VertexIterator const  end,
Ellipse const &  e 
)

Definition at line 402 of file ConvexPolygonImpl.h.

405 {
406  return relate(begin, end, e.getBoundingCircle()) & (CONTAINS | DISJOINT);
407 }

◆ relate() [5/5]

template<typename VertexIterator1 , typename VertexIterator2 >
Relationship lsst::sphgeom::detail::relate ( VertexIterator1 const  begin1,
VertexIterator1 const  end1,
VertexIterator2 const  begin2,
VertexIterator2 const  end2 
)

Definition at line 336 of file ConvexPolygonImpl.h.

340 {
341  // TODO(smm): Make this more performant. Instead of the current quadratic
342  // implementation, it should be possible to determine whether the boundaries
343  // intersect by adapting the following method to the sphere:
344  //
345  // A new linear algorithm for intersecting convex polygons
346  // Computer Graphics and Image Processing, Volume 19, Issue 1, May 1982, Page 92
347  // Joseph O'Rourke, Chi-Bin Chien, Thomas Olson, David Naddor
348  //
349  // http://www.sciencedirect.com/science/article/pii/0146664X82900235
350  bool all1 = true;
351  bool any1 = false;
352  bool all2 = true;
353  bool any2 = false;
354  for (VertexIterator1 i = begin1; i != end1; ++i) {
355  bool b = contains(begin2, end2, *i);
356  all1 = b && all1;
357  any1 = b || any1;
358  }
359  for (VertexIterator2 j = begin2; j != end2; ++j) {
360  bool b = contains(begin1, end1, *j);
361  all2 = b && all2;
362  any2 = b || any2;
363  }
364  if (all1 || all2) {
365  // All vertices of one or both polygons are inside the other
366  return (all1 ? WITHIN : INTERSECTS) | (all2 ? CONTAINS : INTERSECTS);
367  }
368  if (any1 || any2) {
369  // The polygons have at least one point in common.
370  return INTERSECTS;
371  }
372  // No vertex of either polygon is inside the other. Consider all
373  // possible edge pairs and look for a crossing.
374  for (VertexIterator1 a = std::prev(end1), b = begin1;
375  b != end1; a = b, ++b) {
376  for (VertexIterator2 c = std::prev(end2), d = begin2;
377  d != end2; c = d, ++d) {
378  int acd = orientation(*a, *c, *d);
379  int bdc = orientation(*b, *d, *c);
380  if (acd == bdc && acd != 0) {
381  int cba = orientation(*c, *b, *a);
382  int dab = orientation(*d, *a, *b);
383  if (cba == dab && cba == acd) {
384  // Found a non-degenerate edge crossing
385  return INTERSECTS;
386  }
387  }
388  }
389  }
390  return DISJOINT;
391 }