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