LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
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
Angle abs(Angle const &a)
Definition: Angle.h:106
T empty(T... args)
Scalar getA() const
getA returns the lower endpoint of this interval.
Definition: Interval.h:76
NormalizedAngle getB() const
getB returns the second endpoint of this interval.
AngleInterval represents closed intervals of arbitrary angles.
Definition: AngleInterval.h:39
table::Key< int > b
bool valid(int32_t chunkId) const
Return &#39;true&#39; if the specified chunk number is valid.
Definition: Chunker.cc:285
int y
Definition: SpanSet.cc:49
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< int32_t > subChunkIds
Definition: Chunker.h:45
T floor(T... args)
double sin(Angle const &a)
Definition: Angle.h:102
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition: Box.h:54
virtual Box getBoundingBox() const =0
getBoundingBox returns a bounding-box for this region.
double cos(Angle const &a)
Definition: Angle.h:103
T atan2(T... args)
bool wraps() const
wraps returns true if the interval "wraps" around the 0/2π angle discontinuity, i.e.
T min(T... args)
T at(T... args)
T push_back(T... args)
Box dilatedBy(Angle r) const
Definition: Box.h:273
A base class for image defects.
Region is a minimal interface for 2-dimensional regions on the unit sphere.
Definition: Region.h:79
NormalizedAngleInterval const & getLon() const
getLon returns the longitude interval of this box.
Definition: Box.h:145
constexpr double PI
Definition: constants.h:36
virtual Relationship relate(Region const &) const =0
Chunker(int32_t numStripes, int32_t numSubStripesPerStripe)
Definition: Chunker.cc:57
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.
T fabs(T... args)
T max(T... args)
solver_t * s
std::vector< SubChunks > getSubChunksIntersecting(Region const &r) const
getSubChunksIntersecting returns all the sub-chunks that potentially intersect the given region...
Definition: Chunker.cc:148
double x
Scalar getB() const
getB returns the upper endpoint of this interval.
Definition: Interval.h:80
Angle represents an angle in radians.
Definition: Angle.h:43
AngleInterval const & getLat() const
getLat returns the latitude interval of this box.
Definition: Box.h:148
NormalizedAngle getA() const
getA returns the first endpoint of this interval.
T back(T... args)
lsst::geom::Angle Angle
Definition: misc.h:33
T sqrt(T... args)
std::vector< int32_t > getAllChunks() const
getAllChunks returns the complete set of chunk IDs for the unit sphere.
Definition: Chunker.cc:258
This file declares a class for partitioning the sky into chunks and sub-chunks.
T reserve(T... args)
SubChunks represents a set of sub-chunks of a particular chunk.
Definition: Chunker.h:43