LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
PixelFinder.h
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2016 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 
23 #ifndef LSST_SPHGEOM_PIXELFINDER_H_
24 #define LSST_SPHGEOM_PIXELFINDER_H_
25 
28 
29 #include "lsst/sphgeom/RangeSet.h"
30 
31 #include "ConvexPolygonImpl.h"
32 
33 
34 namespace lsst {
35 namespace sphgeom {
36 namespace detail {
37 
38 // `PixelFinder` is a CRTP base class that locates pixels intersecting a
39 // region. It assumes a hierarchical pixelization, and that pixels are
40 // convex spherical polygons with a fixed number of vertices.
41 //
42 // The algorithm used is top-down tree traversal, implemented via recursion
43 // for simplicity. Subclasses must provide a method named `subdivide` with
44 // the following signature:
45 //
46 // void subdivide(UnitVector3d const * pixel,
47 // uint64_t index,
48 // int level);
49 //
50 // that subdivides a pixel into its children and then invokes visit() on
51 // each child. Children should be visited in ascending index order to keep
52 // RangeSet inserts efficient. The subclass is also responsible for
53 // implementing a top-level method that invokes visit() on each root pixel,
54 // or on some set of candidate pixels.
55 //
56 // The `RegionType` parameter avoids the need for virtual function calls to
57 // determine the spatial relationship between pixels and the input region. The
58 // boolean template parameter `InteriorOnly` is a flag that indicates whether
59 // to locate all pixels that intersect the input region, or only those that
60 // are entirely inside it. Finally, the `NumVertices` template parameter is
61 // the number of vertices in the polygonal representation of a pixel.
62 template <
63  typename Derived,
64  typename RegionType,
65  bool InteriorOnly,
66  size_t NumVertices
67 >
68 class PixelFinder {
69 public:
71  RegionType const & region,
72  int level,
73  size_t maxRanges):
74  _ranges{&ranges},
75  _region{&region},
76  _level{level},
77  _desiredLevel{level},
78  _maxRanges{maxRanges == 0 ? maxRanges - 1 : maxRanges}
79  {}
80 
81  void visit(UnitVector3d const * pixel,
82  uint64_t index,
83  int level)
84  {
85  if (level > _level) {
86  // Nothing to do - the subdivision level has been reduced
87  // or a pixel that completely contains the search region
88  // has been found.
89  return;
90  }
91  // Determine the relationship between the pixel and the search region.
92  Relationship r = detail::relate(pixel, pixel + NumVertices, *_region);
93  if ((r & DISJOINT) != 0) {
94  // The pixel is disjoint from the search region.
95  return;
96  }
97  if ((r & WITHIN) != 0) {
98  // The tree traversal has reached a pixel that is entirely within
99  // the search region.
100  _insert(index, level);
101  return;
102  } else if (level == _level) {
103  // The tree traversal has reached a leaf.
104  if (!InteriorOnly) {
105  _insert(index, level);
106  }
107  return;
108  }
109  static_cast<Derived *>(this)->subdivide(pixel, index, level);
110  }
111 
112 private:
113  RangeSet * _ranges;
114  RegionType const * _region;
115  int _level;
116  int const _desiredLevel;
117  size_t const _maxRanges;
118 
119  void _insert(uint64_t index, int level) {
120  int shift = 2 * (_desiredLevel - level);
121  _ranges->insert(index << shift, (index + 1) << shift);
122  while (_ranges->size() > _maxRanges) {
123  // Reduce the subdivision level.
124  --_level;
125  shift += 2;
126  // When looking for intersecting pixels, ranges are simplified
127  // by expanding them outwards, causing nearly adjacent small ranges
128  // to merge.
129  //
130  // When looking for interior pixels, ranges are simplified by
131  // shrinking them inwards, causing small ranges to disappear.
132  if (InteriorOnly) {
133  _ranges->complement();
134  }
135  _ranges->simplify(shift);
136  if (InteriorOnly) {
137  _ranges->complement();
138  }
139  }
140  }
141 };
142 
143 
144 // `findPixels` implements pixel-finding for an arbitrary Region, given a
145 // PixelFinder subclass for a specific pixelization.
146 template <
147  template <typename, bool> class Finder,
148  bool InteriorOnly
149 >
150 RangeSet findPixels(Region const & r, size_t maxRanges, int level) {
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 }
172 
173 }}} // namespace lsst::sphgeom::detail
174 
175 #endif // LSST_SPHGEOM_PIXELFINDER_H_
This file contains the meat of the ConvexPolygon implementation.
table::PointKey< int > pixel
This file provides a type for representing integer sets.
table::Key< int > b
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition: Box.h:54
Circle is a circular region on the unit sphere that contains its boundary.
Definition: Circle.h:46
ConvexPolygon is a closed convex polygon on the unit sphere.
Definition: ConvexPolygon.h:57
Ellipse is an elliptical region on the sphere.
Definition: Ellipse.h:169
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
Definition: Ellipse.cc:241
A RangeSet is a set of unsigned 64 bit integers.
Definition: RangeSet.h:99
void insert(uint64_t u)
Definition: RangeSet.h:285
RangeSet & complement()
complement replaces this set S with U ∖ S, where U is the universe of range sets, [0,...
Definition: RangeSet.h:328
size_t size() const
size returns the number of ranges in this set.
Definition: RangeSet.h:539
RangeSet & simplify(uint32_t n)
simplify simplifies this range set by "coarsening" its ranges.
Definition: RangeSet.cc:308
Region is a minimal interface for 2-dimensional regions on the unit sphere.
Definition: Region.h:79
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
void visit(UnitVector3d const *pixel, uint64_t index, int level)
Definition: PixelFinder.h:81
PixelFinder(RangeSet &ranges, RegionType const &region, int level, size_t maxRanges)
Definition: PixelFinder.h:70
RangeSet findPixels(Region const &r, size_t maxRanges, int level)
Definition: PixelFinder.h:150
Relationship relate(VertexIterator const begin, VertexIterator const end, Box const &b)
A base class for image defects.