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
Mq3cPixelization.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
35#include <stdexcept>
36
38#include "lsst/sphgeom/curve.h"
40
41#include "PixelFinder.h"
42#include "Q3cPixelizationImpl.h"
43
44
45namespace lsst {
46namespace sphgeom {
47
48namespace {
49
50// See commentary in Q3cPixelizationImpl.h for an explanation of
51// these lookup tables.
52
53constexpr uint8_t UNUSED = 255;
54
55alignas(64) uint8_t const FACE_NUM[64] = {
56 4, 4, 4, 4, UNUSED, 3, UNUSED, UNUSED,
57 UNUSED, UNUSED, 0, UNUSED, UNUSED, UNUSED, UNUSED, UNUSED,
58 UNUSED, UNUSED, UNUSED, 2, UNUSED, 3, UNUSED, 2,
59 UNUSED, UNUSED, 0, 2, UNUSED, UNUSED, UNUSED, 2,
60 5, UNUSED, UNUSED, UNUSED, 5, 3, UNUSED, UNUSED,
61 5, UNUSED, 0, UNUSED, 5, UNUSED, UNUSED, UNUSED,
62 UNUSED, UNUSED, UNUSED, UNUSED, UNUSED, 3, UNUSED, UNUSED,
63 UNUSED, UNUSED, 0, UNUSED, 1, 1, 1, 1
64};
65
66uint8_t const FACE_COMP[6][4] = {
67 {0, 1, 2, UNUSED}, {1, 2, 0, UNUSED}, {2, 0, 1, UNUSED},
68 {0, 1, 2, UNUSED}, {1, 2, 0, UNUSED}, {2, 0, 1, UNUSED}
69};
70
71double const FACE_CONST[6][4] = {
72 { 1.0, 1.0, -1.0, 0.0},
73 { 1.0, 1.0, 1.0, 0.0},
74 { 1.0, -1.0, 1.0, 0.0},
75 {-1.0, -1.0, 1.0, 0.0},
76 {-1.0, -1.0, -1.0, 0.0},
77 {-1.0, 1.0, -1.0, 0.0}
78};
79
80// TODO: Fix and document this constant!
81constexpr double DILATION = 1.0e-15;
82
83// `wrapIndex` returns the modified-Q3C index for grid coordinates (face, s, t)
84// at the given level. Both s and t may underflow or overflow by 1, i.e. wrap
85// to an adjacent face.
86uint64_t wrapIndex(int level,
87 int face,
88 uint32_t s,
89 uint32_t t)
90{
91 uint32_t const stMax = (static_cast<uint32_t>(1) << level) - 1;
92 // Wrap until no more underflow or overflow is detected.
93 while (true) {
94 if (s == static_cast<uint32_t>(-1)) {
95 face = (face + 4) % 6;
96 s = stMax - t;
97 t = stMax;
98 continue;
99 } else if (s > stMax) {
100 face = (face + 1) % 6;
101 s = t;
102 t = 0;
103 continue;
104 } else if (t == static_cast<uint32_t>(-1)) {
105 face = (face + 5) % 6;
106 t = s;
107 s = stMax;
108 continue;
109 } else if (t > stMax) {
110 face = (face + 2) % 6;
111 t = stMax - s;
112 s = 0;
113 continue;
114 }
115 break;
116 }
117 return (static_cast<uint64_t>(face + 10) << (2 * level)) |
118 hilbertIndex(s, t, level);
119}
120
121int findNeighborhood(int level, uint64_t i, uint64_t * dst) {
122 int const face = static_cast<int>(i >> (2 * level)) - 10;
123 uint32_t s, t;
124 std::tie(s, t) = hilbertIndexInverse(i, level);
125 dst[0] = wrapIndex(level, face, s - 1, t - 1);
126 dst[1] = wrapIndex(level, face, s , t - 1);
127 dst[2] = wrapIndex(level, face, s + 1, t - 1);
128 dst[3] = wrapIndex(level, face, s - 1, t);
129 dst[4] = i;
130 dst[5] = wrapIndex(level, face, s + 1, t);
131 dst[6] = wrapIndex(level, face, s - 1, t + 1);
132 dst[7] = wrapIndex(level, face, s , t + 1);
133 dst[8] = wrapIndex(level, face, s + 1, t + 1);
134 std::sort(dst, dst + 9);
135 return static_cast<int>(std::unique(dst, dst + 9) - dst);
136}
137
138#if defined(NO_SIMD) || !defined(__x86_64__)
139 void makeQuad(uint64_t i, int level, UnitVector3d * verts) {
140 int const face = static_cast<int>(i >> (2 * level)) - 10;
141 double const faceScale = FACE_SCALE[level];
142 double u0, v0;
143 uint32_t s, t;
144 std::tie(s, t) = hilbertIndexInverse(i, level);
145 std::tie(u0, v0) = gridToFace(
146 level, static_cast<int32_t>(s), static_cast<int32_t>(t));
147 double u1 = (u0 + faceScale) + DILATION;
148 double v1 = (v0 + faceScale) + DILATION;
149 u0 -= DILATION;
150 v0 -= DILATION;
151 std::tie(u0, v0) = atanApproxInverse(u0, v0);
152 std::tie(u1, v1) = atanApproxInverse(u1, v1);
153 verts[0] = faceToSphere(face, u0, v0, FACE_COMP, FACE_CONST);
154 verts[1] = faceToSphere(face, u1, v0, FACE_COMP, FACE_CONST);
155 verts[2] = faceToSphere(face, u1, v1, FACE_COMP, FACE_CONST);
156 verts[3] = faceToSphere(face, u0, v1, FACE_COMP, FACE_CONST);
157 // Even face numbers have right-handed uv coordinate systems,
158 // odd face numbers have left-handed ones. This has to be taken
159 // into account when generating pixel vertices, since convex
160 // polygon vertices must be specified in counter-clockwise order.
161 if ((face & 1) == 0) {
162 std::swap(verts[1], verts[3]);
163 }
164 }
165#else
166 void makeQuad(uint64_t i, int level, UnitVector3d * verts) {
167 int const face = static_cast<int>(i >> (2 * level)) - 10;
168 __m128d faceScale = _mm_set1_pd(FACE_SCALE[level]);
169 __m128d dilation = _mm_set1_pd(DILATION);
170 __m128d u0v0 = gridToFace(level, hilbertIndexInverseSimd(i, level));
171 __m128d u1v1 = _mm_add_pd(u0v0, faceScale);
172 u0v0 = atanApproxInverse(_mm_sub_pd(u0v0, dilation));
173 u1v1 = atanApproxInverse(_mm_add_pd(u1v1, dilation));
174 verts[0] = faceToSphere(face, u0v0, FACE_COMP, FACE_CONST);
175 verts[1] = faceToSphere(face, _mm_shuffle_pd(u1v1, u0v0, 2),
176 FACE_COMP, FACE_CONST);
177 verts[2] = faceToSphere(face, u1v1, FACE_COMP, FACE_CONST);
178 verts[3] = faceToSphere(face, _mm_shuffle_pd(u0v0, u1v1, 2),
179 FACE_COMP, FACE_CONST);
180 if ((face & 1) == 0) {
181 std::swap(verts[1], verts[3]);
182 }
183 }
184#endif
185
186
187// `Mq3cPixelFinder` locates modified-Q3C pixels that intersect a region.
188//
189// For now, we always begin with a loop over the root cube faces. For small
190// regions, this could be made significantly faster by computing the modified
191// Q3C index of the region bounding circle center, and looping over that pixel
192// and its neighbors.
193//
194// The subdivision level for the initial index computation would have to be
195// chosen such that the 8 or 9 pixel neighborhood of the center pixel is
196// guaranteed to contain the bounding circle. The minimum angular pixel width
197// could be precomputed per level. Alternatively, there is some constant W
198// such that the minimum angle between two points separated by at least one
199// pixel is greater than W * 2^-L at level L. Given the bounding circle
200// radius R, the subdivision level L of the initial neighborhood is the binary
201// exponent of W/R (and can be extracted by calling std::frexp).
202//
203// Finding W and implementing the above is left as a future optimization.
204template <typename RegionType, bool InteriorOnly>
205class Mq3cPixelFinder: public detail::PixelFinder<
206 Mq3cPixelFinder<RegionType, InteriorOnly>, RegionType, InteriorOnly, 4>
207{
208private:
209 using Base = detail::PixelFinder<
210 Mq3cPixelFinder<RegionType, InteriorOnly>, RegionType, InteriorOnly, 4>;
211 using Base::visit;
212
213public:
214 Mq3cPixelFinder(RangeSet & ranges,
215 RegionType const & region,
216 int level,
217 size_t maxRanges):
218 Base(ranges, region, level, maxRanges)
219 {}
220
221 void operator()() {
222 UnitVector3d pixel[4];
223 // Loop over cube faces
224 for (uint64_t f = 10; f < 16; ++f) {
225 makeQuad(f, 0, pixel);
226 visit(pixel, f, 0);
227 }
228 }
229
230 void subdivide(UnitVector3d const *, uint64_t i, int level) {
231 UnitVector3d pixel[4];
232 ++level;
233 for (uint64_t c = i * 4; c != i * 4 + 4; ++c) {
234 makeQuad(c, level, pixel);
235 visit(pixel, c, level);
236 }
237 }
238};
239
240} // unnamed namespace
241
242
243int Mq3cPixelization::level(uint64_t i) {
244 // A modified Q3C index consists of 4 bits encoding the root cube face
245 // (10 - 15), followed by 2l bits, where each of the l bit pairs encodes
246 // a child index (0-3), and l is the desired level.
247 int j = log2(i);
248 // The level l is derivable from the index j of the MSB of i. For i to
249 // be valid, j must be an odd integer > 1 and the upper 4 bits of i must
250 // be at least 10.
251 if ((j & 1) == 0 || (j == 1) || ((i >> (j - 3)) < 10)) {
252 return -1;
253 }
254 return (j - 3) >> 1;
255}
256
258 int l = level(i);
259 if (l < 0 || l > MAX_LEVEL) {
260 throw std::invalid_argument("Invalid modified-Q3C index");
261 }
262 UnitVector3d verts[4];
263 makeQuad(i, l, verts);
264 return ConvexPolygon(verts[0], verts[1], verts[2], verts[3]);
265}
266
268 int l = level(i);
269 if (l < 0 || l > MAX_LEVEL) {
270 throw std::invalid_argument("Invalid modified-Q3C index");
271 }
272 uint64_t indexes[9];
273 int n = findNeighborhood(l, i, indexes);
274 return std::vector<uint64_t>(indexes, indexes + n);
275}
276
278 static char const FACE_NORM[6][2] = {
279 {'-', 'Z'}, {'+', 'X'}, {'+', 'Y'},
280 {'+', 'Z'}, {'-', 'X'}, {'-', 'Y'},
281 };
282 char s[MAX_LEVEL + 2];
283 int l = level(i);
284 if (l < 0 || l > MAX_LEVEL) {
285 throw std::invalid_argument("Invalid modified-Q3C index");
286 }
287 // Print in base-4, from least to most significant digit.
288 char * p = s + (sizeof(s) - 1);
289 for (; l > 0; --l, --p, i >>= 2) {
290 *p = '0' + (i & 3);
291 }
292 // The remaining bits correspond to the cube face.
293 --p;
294 p[0] = FACE_NORM[i - 10][0];
295 p[1] = FACE_NORM[i - 10][1];
296 return std::string(p, sizeof(s) - static_cast<size_t>(p - s));
297}
298
299Mq3cPixelization::Mq3cPixelization(int level) : _level{level} {
300 if (level < 0 || level > MAX_LEVEL) {
302 "Modified-Q3C subdivision level not in [0, 30]");
303 }
304}
305
307 uint64_t f = i >> (2 * _level);
308 if (f < 10 || f > 15) {
309 throw std::invalid_argument("Invalid modified-Q3C index");
310 }
311 UnitVector3d verts[4];
312 makeQuad(i, _level, verts);
314 new ConvexPolygon(verts[0], verts[1], verts[2], verts[3]));
315}
316
317#if defined(NO_SIMD) || !defined(__x86_64__)
318 uint64_t Mq3cPixelization::index(UnitVector3d const & p) const {
319 int face = faceNumber(p, FACE_NUM);
320 double w = std::fabs(p(FACE_COMP[face][2]));
321 double u = (p(FACE_COMP[face][0]) / w) * FACE_CONST[face][0];
322 double v = (p(FACE_COMP[face][1]) / w) * FACE_CONST[face][1];
323 std::tie(u, v) = atanApprox(u, v);
324 std::tuple<int32_t, int32_t> g = faceToGrid(_level, u, v);
325 uint64_t h = hilbertIndex(static_cast<uint32_t>(std::get<0>(g)),
326 static_cast<uint32_t>(std::get<1>(g)),
327 _level);
328 return (static_cast<uint64_t>(face + 10) << (2 * _level)) | h;
329 }
330#else
331 uint64_t Mq3cPixelization::index(UnitVector3d const & p) const {
332 int face = faceNumber(p, FACE_NUM);
333 __m128d ww = _mm_set1_pd(p(FACE_COMP[face][2]));
334 __m128d uv = _mm_set_pd(p(FACE_COMP[face][1]), p(FACE_COMP[face][0]));
335 uv = _mm_mul_pd(
336 _mm_div_pd(uv, _mm_andnot_pd(_mm_set_pd(-0.0, -0.0), ww)),
337 _mm_set_pd(FACE_CONST[face][1], FACE_CONST[face][0])
338 );
339 __m128i st = faceToGrid(_level, atanApprox(uv));
340 uint64_t h = hilbertIndex(st, _level);
341 return (static_cast<uint64_t>(face + 10) << (2 * _level)) | h;
342 }
343#endif
344
345RangeSet Mq3cPixelization::_envelope(Region const & r, size_t maxRanges) const {
346 return detail::findPixels<Mq3cPixelFinder, false>(r, maxRanges, _level);
347}
348
349RangeSet Mq3cPixelization::_interior(Region const & r, size_t maxRanges) const {
350 return detail::findPixels<Mq3cPixelFinder, true>(r, maxRanges, _level);
351}
352
353}} // namespace lsst::sphgeom
This file declares a class for representing convex polygons with great circle edges on the unit spher...
This file declares a Pixelization subclass for the modified Q3C indexing scheme.
This file provides a base class for pixel finders.
This file contains functions used by Q3C pixelization implementations.
This file declares a class for representing unit vectors in ℝ³.
ConvexPolygon is a closed convex polygon on the unit sphere.
static int level(uint64_t i)
level returns the subdivision level of the given modified Q3C index.
static constexpr int MAX_LEVEL
The maximum supported cube-face grid resolution is 2^30 by 2^30.
static ConvexPolygon quad(uint64_t i)
quad returns the quadrilateral corresponding to the modified Q3C pixel with index i.
RangeSet _envelope(Region const &r, size_t maxRanges) const override
RangeSet _interior(Region const &r, size_t maxRanges) const override
static std::string asString(uint64_t i)
toString converts the given modified-Q3C index to a human readable string.
uint64_t index(UnitVector3d const &v) const override
index computes the index of the pixel for v.
Mq3cPixelization(int level)
This constructor creates a modified Q3C pixelization of the sphere with the given subdivision level.
static std::vector< uint64_t > neighborhood(uint64_t i)
neighborhood returns the indexes of all pixels that share a vertex with pixel i (including i itself).
std::unique_ptr< Region > pixel(uint64_t i) const override
pixel returns the spherical region corresponding to the pixel with index i.
A RangeSet is a set of unsigned 64 bit integers.
Definition RangeSet.h:106
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
This file contains functions for space-filling curves.
T fabs(T... args)
std::tuple< uint32_t, uint32_t > hilbertIndexInverse(uint64_t h, int m)
hilbertIndexInverse returns the point (x, y) with Hilbert index h, where x and y are m bit integers.
Definition curve.h:368
uint8_t log2(uint64_t x)
Definition curve.h:105
uint64_t hilbertIndex(uint32_t x, uint32_t y, int m)
hilbertIndex returns the index of (x, y) in a 2-D Hilbert curve.
Definition curve.h:356
T sort(T... args)
double w
Definition CoaddPsf.cc:69
T swap(T... args)
T tie(T... args)
T unique(T... args)