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
PixelFinder.h
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
30#ifndef LSST_SPHGEOM_PIXELFINDER_H_
31#define LSST_SPHGEOM_PIXELFINDER_H_
32
35
37
38#include "ConvexPolygonImpl.h"
39
40
41namespace lsst {
42namespace sphgeom {
43namespace detail {
44
45// `PixelFinder` is a CRTP base class that locates pixels intersecting a
46// region. It assumes a hierarchical pixelization, and that pixels are
47// convex spherical polygons with a fixed number of vertices.
48//
49// The algorithm used is top-down tree traversal, implemented via recursion
50// for simplicity. Subclasses must provide a method named `subdivide` with
51// the following signature:
52//
53// void subdivide(UnitVector3d const * pixel,
54// uint64_t index,
55// int level);
56//
57// that subdivides a pixel into its children and then invokes visit() on
58// each child. Children should be visited in ascending index order to keep
59// RangeSet inserts efficient. The subclass is also responsible for
60// implementing a top-level method that invokes visit() on each root pixel,
61// or on some set of candidate pixels.
62//
63// The `RegionType` parameter avoids the need for virtual function calls to
64// determine the spatial relationship between pixels and the input region. The
65// boolean template parameter `InteriorOnly` is a flag that indicates whether
66// to locate all pixels that intersect the input region, or only those that
67// are entirely inside it. Finally, the `NumVertices` template parameter is
68// the number of vertices in the polygonal representation of a pixel.
69template <
70 typename Derived,
71 typename RegionType,
72 bool InteriorOnly,
73 size_t NumVertices
74>
76public:
78 RegionType const & region,
79 int level,
80 size_t maxRanges):
81 _ranges{&ranges},
82 _region{&region},
83 _level{level},
84 _desiredLevel{level},
85 _maxRanges{maxRanges == 0 ? maxRanges - 1 : maxRanges}
86 {}
87
88 void visit(UnitVector3d const * pixel,
89 uint64_t index,
90 int level)
91 {
92 if (level > _level) {
93 // Nothing to do - the subdivision level has been reduced
94 // or a pixel that completely contains the search region
95 // has been found.
96 return;
97 }
98 // Determine the relationship between the pixel and the search region.
99 Relationship r = detail::relate(pixel, pixel + NumVertices, *_region);
100 if ((r & DISJOINT) != 0) {
101 // The pixel is disjoint from the search region.
102 return;
103 }
104 if ((r & WITHIN) != 0) {
105 // The tree traversal has reached a pixel that is entirely within
106 // the search region.
107 _insert(index, level);
108 return;
109 } else if (level == _level) {
110 // The tree traversal has reached a leaf.
111 if (!InteriorOnly) {
112 _insert(index, level);
113 }
114 return;
115 }
116 static_cast<Derived *>(this)->subdivide(pixel, index, level);
117 }
118
119private:
120 RangeSet * _ranges;
121 RegionType const * _region;
122 int _level;
123 int const _desiredLevel;
124 size_t const _maxRanges;
125
126 void _insert(uint64_t index, int level) {
127 int shift = 2 * (_desiredLevel - level);
128 _ranges->insert(index << shift, (index + 1) << shift);
129 while (_ranges->size() > _maxRanges) {
130 // Reduce the subdivision level.
131 --_level;
132 shift += 2;
133 // When looking for intersecting pixels, ranges are simplified
134 // by expanding them outwards, causing nearly adjacent small ranges
135 // to merge.
136 //
137 // When looking for interior pixels, ranges are simplified by
138 // shrinking them inwards, causing small ranges to disappear.
139 if (InteriorOnly) {
140 _ranges->complement();
141 }
142 _ranges->simplify(shift);
143 if (InteriorOnly) {
144 _ranges->complement();
145 }
146 }
147 }
148};
149
150
151// `findPixels` implements pixel-finding for an arbitrary Region, given a
152// PixelFinder subclass for a specific pixelization.
153template <
154 template <typename, bool> class Finder,
155 bool InteriorOnly
156>
157RangeSet findPixels(Region const & r, size_t maxRanges, int level) {
158 RangeSet s;
159 Circle const * c = nullptr;
160 Ellipse const * e = nullptr;
161 Box const * b = nullptr;
162 if ((c = dynamic_cast<Circle const *>(&r))) {
163 Finder<Circle, InteriorOnly> find(s, *c, level, maxRanges);
164 find();
165 } else if ((e = dynamic_cast<Ellipse const *>(&r))) {
166 Finder<Circle, InteriorOnly> find(
167 s, e->getBoundingCircle(), level, maxRanges);
168 find();
169 } else if ((b = dynamic_cast<Box const *>(&r))) {
170 Finder<Box, InteriorOnly> find(s, *b, level, maxRanges);
171 find();
172 } else {
173 Finder<ConvexPolygon, InteriorOnly> find(
174 s, dynamic_cast<ConvexPolygon const &>(r), level, maxRanges);
175 find();
176 }
177 return s;
178}
179
180}}} // namespace lsst::sphgeom::detail
181
182#endif // LSST_SPHGEOM_PIXELFINDER_H_
This file contains the meat of the ConvexPolygon implementation.
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:61
Circle is a circular region on the unit sphere that contains its boundary.
Definition Circle.h:53
ConvexPolygon is a closed convex polygon on the unit sphere.
Ellipse is an elliptical region on the sphere.
Definition Ellipse.h:177
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
Definition Ellipse.cc:248
A RangeSet is a set of unsigned 64 bit integers.
Definition RangeSet.h:106
void insert(uint64_t u)
Definition RangeSet.h:292
RangeSet & complement()
complement replaces this set S with U ∖ S, where U is the universe of range sets, [0,...
Definition RangeSet.h:335
size_t size() const
size returns the number of ranges in this set.
Definition RangeSet.h:546
RangeSet & simplify(uint32_t n)
simplify simplifies this range set by "coarsening" its ranges.
Definition RangeSet.cc:315
Region is a minimal interface for 2-dimensional regions on the unit sphere.
Definition Region.h:86
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
void visit(UnitVector3d const *pixel, uint64_t index, int level)
Definition PixelFinder.h:88
PixelFinder(RangeSet &ranges, RegionType const &region, int level, size_t maxRanges)
Definition PixelFinder.h:77
RangeSet findPixels(Region const &r, size_t maxRanges, int level)
Relationship relate(VertexIterator const begin, VertexIterator const end, Box const &b)