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
Chunker.cc
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
32
34
35namespace lsst {
36namespace sphgeom {
37
38namespace {
39
40int32_t computeNumSegments(AngleInterval const & latitudes, Angle width) {
41 // TODO(smm): Document this.
42 if (width.asRadians() > PI) {
43 return 1;
44 }
45 Angle maxAbsLat = std::max(abs(latitudes.getA()), abs(latitudes.getB()));
46 if (maxAbsLat.asRadians() > 0.5 * PI - 4.85e-6) {
47 return 1;
48 }
49 double cosWidth = cos(width);
50 double sinLat = sin(maxAbsLat);
51 double cosLat = cos(maxAbsLat);
52 double x = cosWidth - sinLat * sinLat;
53 double u = cosLat * cosLat;
54 double y = std::sqrt(std::fabs(u * u - x * x));
55 return static_cast<int32_t>(
56 std::floor(2.0 * PI / std::fabs(std::atan2(y, x))));
57}
58
59constexpr double BOX_EPSILON = 5.0e-12; // ~1 micro-arcsecond
60
61} // unnamed namespace
62
63
64Chunker::Chunker(int32_t numStripes,
65 int32_t numSubStripesPerStripe) :
66 _numStripes(numStripes),
67 _numSubStripesPerStripe(numSubStripesPerStripe),
68 _numSubStripes(numStripes * numSubStripesPerStripe),
69 _maxSubChunksPerSubStripeChunk(0),
70 _subStripeHeight(Angle(PI) / _numSubStripes)
71{
72 if (numStripes < 1 || numSubStripesPerStripe < 1) {
73 throw std::runtime_error("The number of stripes and sub-stripes "
74 "per stripe must be positive");
75 }
76 if (numStripes * numSubStripesPerStripe > 180*3600) {
77 throw std::runtime_error("Sub-stripes are too small");
78 }
79 Angle const stripeHeight = Angle(PI) / _numStripes;
80 _stripes.reserve(_numStripes);
81 _subStripes.reserve(_numSubStripes);
82 for (int32_t s = 0; s < _numStripes; ++s) {
83 // Compute stripe latitude bounds.
84 AngleInterval sLat(s * stripeHeight - Angle(0.5 * PI),
85 (s + 1) * stripeHeight - Angle(0.5 * PI));
86 Stripe stripe;
87 int32_t const nc = computeNumSegments(sLat, stripeHeight);
88 stripe.chunkWidth = Angle(2.0 * PI) / nc;
89 stripe.numChunksPerStripe = nc;
90 int32_t ss = s * _numSubStripesPerStripe;
91 int32_t const ssEnd = ss + _numSubStripesPerStripe;
92 for (; ss < ssEnd; ++ss) {
93 // Compute sub-stripe latitude bounds.
94 AngleInterval ssLat(ss * _subStripeHeight - Angle(0.5 * PI),
95 (ss + 1) * _subStripeHeight - Angle(0.5 * PI));
96 SubStripe subStripe;
97 int32_t const nsc = computeNumSegments(ssLat, _subStripeHeight) / nc;
98 stripe.numSubChunksPerChunk += nsc;
99 subStripe.numSubChunksPerChunk = nsc;
100 if (nsc > _maxSubChunksPerSubStripeChunk) {
101 _maxSubChunksPerSubStripeChunk = nsc;
102 }
103 subStripe.subChunkWidth = Angle(2.0 * PI) / (nsc * nc);
104 _subStripes.push_back(subStripe);
105 }
106 _stripes.push_back(stripe);
107 }
108}
109
111 std::vector<int32_t> chunkIds;
112 // Find the stripes that intersect the bounding box of r.
113 Box b = r.getBoundingBox().dilatedBy(Angle(BOX_EPSILON));
114 double ya = std::floor((b.getLat().getA() + Angle(0.5 * PI)) / _subStripeHeight);
115 double yb = std::floor((b.getLat().getB() + Angle(0.5 * PI)) / _subStripeHeight);
116 int32_t minSS = std::min(static_cast<int32_t>(ya), _numSubStripes - 1);
117 int32_t maxSS = std::min(static_cast<int32_t>(yb), _numSubStripes - 1);
118 int32_t minS = minSS / _numSubStripesPerStripe;
119 int32_t maxS = maxSS / _numSubStripesPerStripe;
120 for (int32_t s = minS; s <= maxS; ++s) {
121 // Find the chunks of s that intersect the bounding box of r.
122 Angle chunkWidth = _stripes[s].chunkWidth;
123 int32_t nc = _stripes[s].numChunksPerStripe;
124 double xa = std::floor(b.getLon().getA() / chunkWidth);
125 double xb = std::floor(b.getLon().getB() / chunkWidth);
126 int32_t ca = std::min(static_cast<int32_t>(xa), nc - 1);
127 int32_t cb = std::min(static_cast<int32_t>(xb), nc - 1);
128 if (ca == cb && b.getLon().wraps()) {
129 ca = 0;
130 cb = nc - 1;
131 }
132 // Examine each chunk overlapping the bounding box of r.
133 if (ca <= cb) {
134 for (int32_t c = ca; c <= cb; ++c) {
135 if ((r.relate(getChunkBoundingBox(s, c)) & DISJOINT) == 0) {
136 chunkIds.push_back(_getChunkId(s, c));
137 }
138 }
139 } else {
140 for (int32_t c = 0; c <= cb; ++c) {
141 if ((r.relate(getChunkBoundingBox(s, c)) & DISJOINT) == 0) {
142 chunkIds.push_back(_getChunkId(s, c));
143 }
144 }
145 for (int32_t c = ca; c < nc; ++c) {
146 if ((r.relate(getChunkBoundingBox(s, c)) & DISJOINT) == 0) {
147 chunkIds.push_back(_getChunkId(s, c));
148 }
149 }
150 }
151 }
152 return chunkIds;
153}
154
156 Region const & r) const
157{
159 // Find the stripes that intersect the bounding box of r.
160 Box b = r.getBoundingBox().dilatedBy(Angle(BOX_EPSILON));
161 double ya = std::floor((b.getLat().getA() + Angle(0.5 * PI)) / _subStripeHeight);
162 double yb = std::floor((b.getLat().getB() + Angle(0.5 * PI)) / _subStripeHeight);
163 int32_t minSS = std::min(static_cast<int32_t>(ya), _numSubStripes - 1);
164 int32_t maxSS = std::min(static_cast<int32_t>(yb), _numSubStripes - 1);
165 int32_t minS = minSS / _numSubStripesPerStripe;
166 int32_t maxS = maxSS / _numSubStripesPerStripe;
167 for (int32_t s = minS; s <= maxS; ++s) {
168 // Find the chunks of s that intersect the bounding box of r.
169 Angle chunkWidth = _stripes[s].chunkWidth;
170 int32_t nc = _stripes[s].numChunksPerStripe;
171 double xa = std::floor(b.getLon().getA() / chunkWidth);
172 double xb = std::floor(b.getLon().getB() / chunkWidth);
173 int32_t ca = std::min(static_cast<int32_t>(xa), nc - 1);
174 int32_t cb = std::min(static_cast<int32_t>(xb), nc - 1);
175 if (ca == cb && b.getLon().wraps()) {
176 ca = 0;
177 cb = nc - 1;
178 }
179 // Examine sub-chunks for each chunk overlapping the bounding box of r.
180 if (ca <= cb) {
181 for (int32_t c = ca; c <= cb; ++c) {
182 _getSubChunks(chunks, r, b.getLon(), s, c, minSS, maxSS);
183 }
184 } else {
185 for (int32_t c = 0; c <= cb; ++c) {
186 _getSubChunks(chunks, r, b.getLon(), s, c, minSS, maxSS);
187 }
188 for (int32_t c = ca; c < nc; ++c) {
189 _getSubChunks(chunks, r, b.getLon(), s, c, minSS, maxSS);
190 }
191 }
192 }
193 return chunks;
194}
195
196void Chunker::_getSubChunks(std::vector<SubChunks> & chunks,
197 Region const & r,
198 NormalizedAngleInterval const & lon,
199 int32_t stripe,
200 int32_t chunk,
201 int32_t minSS,
202 int32_t maxSS) const
203{
204 SubChunks subChunks;
205 subChunks.chunkId = _getChunkId(stripe, chunk);
206 if ((r.relate(getChunkBoundingBox(stripe, chunk)) & CONTAINS) != 0) {
207 // r contains the entire chunk, so there is no need to test sub-chunks
208 // for intersection with r.
209 subChunks.subChunkIds = getAllSubChunks(subChunks.chunkId);
210 } else {
211 // Find the sub-stripes to iterate over.
212 minSS = std::max(minSS, stripe * _numSubStripesPerStripe);
213 maxSS = std::min(maxSS, (stripe + 1) * _numSubStripesPerStripe - 1);
214 int32_t const nc = _stripes[stripe].numChunksPerStripe;
215 for (int32_t ss = minSS; ss <= maxSS; ++ss) {
216 // Find the sub-chunks of ss to iterate over.
217 Angle subChunkWidth = _subStripes[ss].subChunkWidth;
218 int32_t const nsc = _subStripes[ss].numSubChunksPerChunk;
219 double xa = std::floor(lon.getA() / subChunkWidth);
220 double xb = std::floor(lon.getB() / subChunkWidth);
221 int32_t sca = std::min(static_cast<int32_t>(xa), nc * nsc - 1);
222 int32_t scb = std::min(static_cast<int32_t>(xb), nc * nsc - 1);
223 if (sca == scb && lon.wraps()) {
224 sca = 0;
225 scb = nc * nsc - 1;
226 }
227 int32_t minSC = chunk * nsc;
228 int32_t maxSC = (chunk + 1) * nsc - 1;
229 // Test each sub-chunk against r, and record those that intersect.
230 if (sca <= scb) {
231 minSC = std::max(sca, minSC);
232 maxSC = std::min(scb, maxSC);
233 for (int32_t sc = minSC; sc <= maxSC; ++sc) {
234 if ((r.relate(getSubChunkBoundingBox(ss, sc)) & DISJOINT) == 0) {
235 subChunks.subChunkIds.push_back(
236 _getSubChunkId(stripe, ss, chunk, sc));
237 }
238 }
239 } else {
240 sca = std::max(sca, minSC);
241 scb = std::min(scb, maxSC);
242 for (int32_t sc = sca; sc <= maxSC; ++sc) {
243 if ((r.relate(getSubChunkBoundingBox(ss, sc)) & DISJOINT) == 0) {
244 subChunks.subChunkIds.push_back(
245 _getSubChunkId(stripe, ss, chunk, sc));
246 }
247 }
248 for (int32_t sc = minSC; sc <= scb; ++sc) {
249 if ((r.relate(getSubChunkBoundingBox(ss, sc)) & DISJOINT) == 0) {
250 subChunks.subChunkIds.push_back(
251 _getSubChunkId(stripe, ss, chunk, sc));
252 }
253 }
254 }
255 }
256 }
257 // If any sub-chunks of this chunk intersect r,
258 // append them to the result vector.
259 if (!subChunks.subChunkIds.empty()) {
260 chunks.push_back(SubChunks());
261 chunks.back().swap(subChunks);
262 }
263}
264
266 std::vector<int32_t> chunkIds;
267 for (int32_t s = 0; s < _numStripes; ++s) {
268 int32_t nc = _stripes[s].numChunksPerStripe;
269 for (int32_t c = 0; c < nc; ++c) {
270 chunkIds.push_back(_getChunkId(s, c));
271 }
272 }
273 return chunkIds;
274}
275
277 std::vector<int32_t> subChunkIds;
278 int32_t s = getStripe(chunkId);
279 subChunkIds.reserve(_stripes.at(s).numSubChunksPerChunk);
280 int32_t const ssBeg = s * _numSubStripesPerStripe;
281 int32_t const ssEnd = ssBeg + _numSubStripesPerStripe;
282 for (int32_t ss = ssBeg; ss < ssEnd; ++ss) {
283 int32_t const scEnd = _subStripes[ss].numSubChunksPerChunk;
284 int32_t const subChunkIdBase = _maxSubChunksPerSubStripeChunk * (ss - ssBeg);
285 for (int32_t sc = 0; sc < scEnd; ++sc) {
286 subChunkIds.push_back(subChunkIdBase + sc);
287 }
288 }
289 return subChunkIds;
290}
291
292bool Chunker::valid(int32_t chunkId) const {
293 int32_t const s = getStripe(chunkId);
294 return s >= 0 and s < _numStripes and
295 getChunk(chunkId, s) < _stripes.at(s).numChunksPerStripe;
296}
297
298Box Chunker::getChunkBoundingBox(int32_t stripe, int32_t chunk) const {
299 Angle chunkWidth = _stripes[stripe].chunkWidth;
300 NormalizedAngleInterval lon(chunkWidth * chunk,
301 chunkWidth * (chunk + 1));
302 int32_t ss = stripe * _numSubStripesPerStripe;
303 int32_t ssEnd = ss + _numSubStripesPerStripe;
304 AngleInterval lat(ss * _subStripeHeight - Angle(0.5 * PI),
305 ssEnd * _subStripeHeight - Angle(0.5 * PI));
306 return Box(lon, lat).dilatedBy(Angle(BOX_EPSILON));
307}
308
309Box Chunker::getSubChunkBoundingBox(int32_t subStripe, int32_t subChunk) const {
310 Angle subChunkWidth = _subStripes[subStripe].subChunkWidth;
311 NormalizedAngleInterval lon(subChunkWidth * subChunk,
312 subChunkWidth * (subChunk + 1));
313 AngleInterval lat(subStripe * _subStripeHeight - Angle(0.5 * PI),
314 (subStripe + 1) * _subStripeHeight - Angle(0.5 * PI));
315 return Box(lon, lat).dilatedBy(Angle(BOX_EPSILON));
316}
317
318}} // namespace lsst::sphgeom
This file declares a class for partitioning the sky into chunks and sub-chunks.
int y
Definition SpanSet.cc:48
table::Key< int > b
T at(T... args)
T atan2(T... args)
T back(T... args)
Angle represents an angle in radians.
Definition Angle.h:50
AngleInterval represents closed intervals of arbitrary angles.
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition Box.h:61
Box dilatedBy(Angle r) const
Definition Box.h:280
std::vector< int32_t > getAllChunks() const
getAllChunks returns the complete set of chunk IDs for the unit sphere.
Definition Chunker.cc:265
int32_t getStripe(int32_t chunkId) const
Return the stripe for the specified chunkId.
Definition Chunker.h:123
Box getSubChunkBoundingBox(int32_t subStripe, int32_t subChunk) const
Definition Chunker.cc:309
int32_t getChunk(int32_t chunkId, int32_t stripe) const
Return the chunk for the given chunkId and stripe.
Definition Chunker.h:128
Chunker(int32_t numStripes, int32_t numSubStripesPerStripe)
Definition Chunker.cc:64
std::vector< int32_t > getAllSubChunks(int32_t chunkId) const
getAllSubChunks returns the complete set of sub-chunk IDs for the given chunk.
Definition Chunker.cc:276
std::vector< SubChunks > getSubChunksIntersecting(Region const &r) const
getSubChunksIntersecting returns all the sub-chunks that potentially intersect the given region.
Definition Chunker.cc:155
Box getChunkBoundingBox(int32_t stripe, int32_t chunk) const
Definition Chunker.cc:298
bool valid(int32_t chunkId) const
Return 'true' if the specified chunk number is valid.
Definition Chunker.cc:292
std::vector< int32_t > getChunksIntersecting(Region const &r) const
getChunksIntersecting returns all the chunks that potentially intersect the given region.
Definition Chunker.cc:110
NormalizedAngleInterval represents closed intervals of normalized angles, i.e.
bool wraps() const
wraps returns true if the interval "wraps" around the 0/2π angle discontinuity, i....
NormalizedAngle getA() const
getA returns the first endpoint of this interval.
NormalizedAngle getB() const
getB returns the second endpoint of this interval.
Region is a minimal interface for 2-dimensional regions on the unit sphere.
Definition Region.h:86
virtual Relationship relate(Region const &) const =0
virtual Box getBoundingBox() const =0
getBoundingBox returns a bounding-box for this region.
T empty(T... args)
T fabs(T... args)
T floor(T... args)
T max(T... args)
T min(T... args)
Angle abs(Angle const &a)
Definition Angle.h:113
double sin(Angle const &a)
Definition Angle.h:109
double cos(Angle const &a)
Definition Angle.h:110
constexpr double PI
Definition constants.h:43
T push_back(T... args)
T reserve(T... args)
T sqrt(T... args)
SubChunks represents a set of sub-chunks of a particular chunk.
Definition Chunker.h:50
std::vector< int32_t > subChunkIds
Definition Chunker.h:52