33 int32_t computeNumSegments(AngleInterval
const & latitudes,
Angle width) {
35 if (width.asRadians() >
PI) {
39 if (maxAbsLat.asRadians() > 0.5 *
PI - 4.85e-6) {
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;
48 return static_cast<int32_t
>(
52 constexpr
double BOX_EPSILON = 5.0e-12;
58 int32_t numSubStripesPerStripe) :
59 _numStripes(numStripes),
60 _numSubStripesPerStripe(numSubStripesPerStripe),
61 _numSubStripes(numStripes * numSubStripesPerStripe),
62 _maxSubChunksPerSubStripeChunk(0),
63 _subStripeHeight(
Angle(
PI) / _numSubStripes)
65 if (numStripes < 1 || numSubStripesPerStripe < 1) {
67 "per stripe must be positive");
69 if (numStripes * numSubStripesPerStripe > 180*3600) {
74 _subStripes.
reserve(_numSubStripes);
75 for (int32_t s = 0; s < _numStripes; ++s) {
78 (s + 1) * stripeHeight -
Angle(0.5 *
PI));
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) {
88 (ss + 1) * _subStripeHeight -
Angle(0.5 *
PI));
90 int32_t
const nsc = computeNumSegments(ssLat, _subStripeHeight) / nc;
91 stripe.numSubChunksPerChunk += nsc;
92 subStripe.numSubChunksPerChunk = nsc;
93 if (nsc > _maxSubChunksPerSubStripeChunk) {
94 _maxSubChunksPerSubStripeChunk = nsc;
96 subStripe.subChunkWidth =
Angle(2.0 *
PI) / (nsc * nc);
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) {
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()) {
127 for (int32_t c = ca; c <= cb; ++c) {
133 for (int32_t c = 0; c <= cb; ++c) {
138 for (int32_t c = ca; c < nc; ++c) {
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) {
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()) {
174 for (int32_t c = ca; c <= cb; ++c) {
175 _getSubChunks(chunks, r,
b.getLon(), s, c, minSS, maxSS);
178 for (int32_t c = 0; c <= cb; ++c) {
179 _getSubChunks(chunks, r,
b.getLon(), s, c, minSS, maxSS);
181 for (int32_t c = ca; c < nc; ++c) {
182 _getSubChunks(chunks, r,
b.getLon(), s, c, minSS, maxSS);
198 subChunks.
chunkId = _getChunkId(stripe, chunk);
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) {
210 Angle subChunkWidth = _subStripes[ss].subChunkWidth;
211 int32_t
const nsc = _subStripes[ss].numSubChunksPerChunk;
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()) {
220 int32_t minSC = chunk * nsc;
221 int32_t maxSC = (chunk + 1) * nsc - 1;
226 for (int32_t sc = minSC; sc <= maxSC; ++sc) {
229 _getSubChunkId(stripe, ss, chunk, sc));
235 for (int32_t sc = sca; sc <= maxSC; ++sc) {
238 _getSubChunkId(stripe, ss, chunk, sc));
241 for (int32_t sc = minSC; sc <= scb; ++sc) {
244 _getSubChunkId(stripe, ss, chunk, sc));
254 chunks.
back().swap(subChunks);
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) {
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);
287 return s >= 0 and s < _numStripes and
288 getChunk(chunkId, s) < _stripes.
at(s).numChunksPerStripe;
292 Angle chunkWidth = _stripes[stripe].chunkWidth;
294 chunkWidth * (chunk + 1));
295 int32_t ss = stripe * _numSubStripesPerStripe;
296 int32_t ssEnd = ss + _numSubStripesPerStripe;
298 ssEnd * _subStripeHeight -
Angle(0.5 *
PI));
303 Angle subChunkWidth = _subStripes[subStripe].subChunkWidth;
305 subChunkWidth * (subChunk + 1));
307 (subStripe + 1) * _subStripeHeight -
Angle(0.5 *
PI));
This file declares a class for partitioning the sky into chunks and sub-chunks.
Angle represents an angle in radians.
AngleInterval represents closed intervals of arbitrary angles.
Box represents a rectangle in spherical coordinate space that contains its boundary.
Box dilatedBy(Angle r) const
std::vector< int32_t > getAllChunks() const
getAllChunks returns the complete set of chunk IDs for the unit sphere.
int32_t getStripe(int32_t chunkId) const
Return the stripe for the specified chunkId.
Box getSubChunkBoundingBox(int32_t subStripe, int32_t subChunk) const
int32_t getChunk(int32_t chunkId, int32_t stripe) const
Return the chunk for the given chunkId and stripe.
Chunker(int32_t numStripes, int32_t numSubStripesPerStripe)
std::vector< int32_t > getAllSubChunks(int32_t chunkId) const
getAllSubChunks returns the complete set of sub-chunk IDs for the given chunk.
std::vector< SubChunks > getSubChunksIntersecting(Region const &r) const
getSubChunksIntersecting returns all the sub-chunks that potentially intersect the given region.
Box getChunkBoundingBox(int32_t stripe, int32_t chunk) const
bool valid(int32_t chunkId) const
Return 'true' if the specified chunk number is valid.
std::vector< int32_t > getChunksIntersecting(Region const &r) const
getChunksIntersecting returns all the chunks that potentially intersect the given region.
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.
virtual Relationship relate(Region const &) const =0
virtual Box getBoundingBox() const =0
getBoundingBox returns a bounding-box for this region.
Angle abs(Angle const &a)
double sin(Angle const &a)
double cos(Angle const &a)
A base class for image defects.
SubChunks represents a set of sub-chunks of a particular chunk.
std::vector< int32_t > subChunkIds