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
Q3cPixelization.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 3, 3, 3, 3, UNUSED, 0, UNUSED, UNUSED,
57 UNUSED, UNUSED, 5, UNUSED, UNUSED, UNUSED, UNUSED, UNUSED,
58 UNUSED, UNUSED, UNUSED, 2, UNUSED, 0, UNUSED, 2,
59 UNUSED, UNUSED, 5, 2, UNUSED, UNUSED, UNUSED, 2,
60 4, UNUSED, UNUSED, UNUSED, 4, 0, UNUSED, UNUSED,
61 4, UNUSED, 5, UNUSED, 4, UNUSED, UNUSED, UNUSED,
62 UNUSED, UNUSED, UNUSED, UNUSED, UNUSED, 0, UNUSED, UNUSED,
63 UNUSED, UNUSED, 5, UNUSED, 1, 1, 1, 1
64};
65
66uint8_t const FACE_COMP[6][4] = {
67 {1, 0, 2, UNUSED}, {1, 2, 0, UNUSED}, {0, 2, 1, UNUSED},
68 {1, 2, 0, UNUSED}, {0, 2, 1, UNUSED}, {1, 0, 2, UNUSED}
69};
70
71alignas(16) double 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// The amount by which pixel boundaries are dilated (in u and v) prior to
81// being mapped to unit vectors. This dilation ensures that the polygonal
82// representation of a pixel contains all unit vectors that map to that
83// pixel.
84//
85// The value is slightly more than the maximum absolute error incurred by
86// mapping face coordinates (u,v) to the unit sphere and back.
87constexpr double DILATION = 1.0e-15;
88
89// `wrapIndex` returns the Q3C index for grid coordinates (face, s, t) at
90// the given level. Both s and t may underflow or overflow by 1, i.e. wrap
91// to an adjacent face.
92uint64_t wrapIndex(int level,
93 int face,
94 uint32_t s,
95 uint32_t t)
96{
97 uint32_t const stMax = (static_cast<uint32_t>(1) << level) - 1;
98 // Wrap until no more underflow or overflow is detected.
99 while (true) {
100 if (s == static_cast<uint32_t>(-1)) {
101 switch (face) {
102 case 0: face = 4; s = stMax - t; t = stMax; break;
103 case 1: face = 4; s = stMax; break;
104 case 2: face = 1; s = stMax; break;
105 case 3: face = 2; s = stMax; break;
106 case 4: face = 3; s = stMax; break;
107 case 5: face = 4; s = t; t = 0; break;
108 default: break;
109 }
110 continue;
111 } else if (s > stMax) {
112 switch (face) {
113 case 0: face = 2; s = t; t = stMax; break;
114 case 1: face = 2; s = 0; break;
115 case 2: face = 3; s = 0; break;
116 case 3: face = 4; s = 0; break;
117 case 4: face = 1; s = 0; break;
118 case 5: face = 2; s = stMax - t; t = 0; break;
119 default: break;
120 }
121 continue;
122 } else if (t == static_cast<uint32_t>(-1)) {
123 switch (face) {
124 case 0: face = 1; t = stMax; break;
125 case 1: face = 5; t = stMax; break;
126 case 2: face = 5; t = stMax - s; s = stMax; break;
127 case 3: face = 5; t = 0; s = stMax - s; break;
128 case 4: face = 5; t = s; s = 0; break;
129 case 5: face = 3; t = 0; s = stMax - s; break;
130 default: break;
131 }
132 continue;
133 } else if (t > stMax) {
134 switch (face) {
135 case 0: face = 3; t = stMax; s = stMax - s; break;
136 case 1: face = 0; t = 0; break;
137 case 2: face = 0; t = s; s = stMax; break;
138 case 3: face = 0; t = stMax; s = stMax - s; break;
139 case 4: face = 0; t = stMax - s; s = 0; break;
140 case 5: face = 1; t = 0; break;
141 default: break;
142 }
143 continue;
144 }
145 break;
146 }
147 return (static_cast<uint64_t>(face) << (2 * level)) | mortonIndex(s, t);
148}
149
150int findNeighborhood(int level, uint64_t i, uint64_t * dst) {
151 uint64_t const mask = (static_cast<uint64_t>(1) << (2 * level)) - 1;
152 int const face = static_cast<int>(i >> (2 * level));
153 uint32_t s, t;
154 std::tie(s, t) = mortonIndexInverse(i & mask);
155 dst[0] = wrapIndex(level, face, s - 1, t - 1);
156 dst[1] = wrapIndex(level, face, s , t - 1);
157 dst[2] = wrapIndex(level, face, s + 1, t - 1);
158 dst[3] = wrapIndex(level, face, s - 1, t);
159 dst[4] = i;
160 dst[5] = wrapIndex(level, face, s + 1, t);
161 dst[6] = wrapIndex(level, face, s - 1, t + 1);
162 dst[7] = wrapIndex(level, face, s , t + 1);
163 dst[8] = wrapIndex(level, face, s + 1, t + 1);
164 std::sort(dst, dst + 9);
165 return static_cast<int>(std::unique(dst, dst + 9) - dst);
166}
167
168#if defined(NO_SIMD) || !defined(__x86_64__)
169 void makeQuad(uint64_t i, int level, UnitVector3d * verts) {
170 uint64_t const mask = (static_cast<uint64_t>(1) << (2 * level)) - 1;
171 int const face = static_cast<int>(i >> (2 * level));
172 double const faceScale = FACE_SCALE[level];
173 double u0, v0;
174 uint32_t s, t;
175 std::tie(s, t) = mortonIndexInverse(i & mask);
176 std::tie(u0, v0) = gridToFace(
177 level, static_cast<int32_t>(s), static_cast<int32_t>(t));
178 double u1 = (u0 + faceScale) + DILATION;
179 double v1 = (v0 + faceScale) + DILATION;
180 u0 -= DILATION;
181 v0 -= DILATION;
182 verts[0] = faceToSphere(face, u0, v0, FACE_COMP, FACE_CONST);
183 verts[1] = faceToSphere(face, u1, v0, FACE_COMP, FACE_CONST);
184 verts[2] = faceToSphere(face, u1, v1, FACE_COMP, FACE_CONST);
185 verts[3] = faceToSphere(face, u0, v1, FACE_COMP, FACE_CONST);
186 }
187#else
188 void makeQuad(uint64_t i, int level, UnitVector3d * verts) {
189 uint64_t const mask = (static_cast<uint64_t>(1) << (2 * level)) - 1;
190 int const face = static_cast<int>(i >> (2 * level));
191 __m128d faceScale = _mm_set1_pd(FACE_SCALE[level]);
192 __m128d dilation = _mm_set1_pd(DILATION);
193 __m128d u0v0 = gridToFace(level, mortonIndexInverseSimd(i & mask));
194 __m128d u1v1 = _mm_add_pd(u0v0, faceScale);
195 u0v0 = _mm_sub_pd(u0v0, dilation);
196 u1v1 = _mm_add_pd(u1v1, dilation);
197 verts[0] = faceToSphere(face, u0v0, FACE_COMP, FACE_CONST);
198 verts[1] = faceToSphere(face, _mm_shuffle_pd(u1v1, u0v0, 2),
199 FACE_COMP, FACE_CONST);
200 verts[2] = faceToSphere(face, u1v1, FACE_COMP, FACE_CONST);
201 verts[3] = faceToSphere(face, _mm_shuffle_pd(u0v0, u1v1, 2),
202 FACE_COMP, FACE_CONST);
203 }
204#endif
205
206
207// `Q3cPixelFinder` locates Q3C pixels that intersect a region.
208//
209// For now, we always begin with a loop over the root cube faces. For small
210// regions, this could be made significantly faster by computing the Q3C
211// index of the region bounding circle center, and looping over that pixel
212// and its neighbors.
213//
214// The subdivision level for the initial index computation would have to be
215// chosen such that the 8 or 9 pixel neighborhood of the center pixel is
216// guaranteed to contain the bounding circle. The minimum angular pixel width
217// could be precomputed per level. Alternatively, the minimum angle (and
218// chord length) between two points separated by at least one pixel can be
219// shown to be greater than √2/3 * 2^-L at level L. Given the bounding circle
220// radius R, the subdivision level L of the initial neighborhood is the binary
221// exponent of √2/3/R (and can be extracted by calling std::frexp).
222//
223// This is left as a future optimization.
224template <typename RegionType, bool InteriorOnly>
225class Q3cPixelFinder: public detail::PixelFinder<
226 Q3cPixelFinder<RegionType, InteriorOnly>, RegionType, InteriorOnly, 4>
227{
228private:
229 using Base = detail::PixelFinder<
230 Q3cPixelFinder<RegionType, InteriorOnly>, RegionType, InteriorOnly, 4>;
231 using Base::visit;
232
233public:
234 Q3cPixelFinder(RangeSet & ranges,
235 RegionType const & region,
236 int level,
237 size_t maxRanges):
238 Base(ranges, region, level, maxRanges)
239 {}
240
241 void operator()() {
242 UnitVector3d pixel[4];
243 // Loop over cube faces
244 for (uint64_t f = 0; f < 6; ++f) {
245 makeQuad(f, 0, pixel);
246 visit(pixel, f, 0);
247 }
248 }
249
250 void subdivide(UnitVector3d const *, uint64_t i, int level) {
251 UnitVector3d pixel[4];
252 ++level;
253 for (uint64_t c = i * 4; c != i * 4 + 4; ++c) {
254 makeQuad(c, level, pixel);
255 visit(pixel, c, level);
256 }
257 }
258};
259
260} // unnamed namespace
261
262
263Q3cPixelization::Q3cPixelization(int level) : _level{level} {
264 if (level < 0 || level > MAX_LEVEL) {
265 throw std::invalid_argument("Q3C subdivision level not in [0, 30]");
266 }
267}
268
270 if (i >= static_cast<uint64_t>(6) << (2 * _level)) {
271 throw std::invalid_argument("Invalid Q3C index");
272 }
273 UnitVector3d verts[4];
274 makeQuad(i, _level, verts);
275 return ConvexPolygon(verts[0], verts[1], verts[2], verts[3]);
276}
277
279 if (i >= static_cast<uint64_t>(6) << (2 * _level)) {
280 throw std::invalid_argument("Invalid Q3C index");
281 }
282 uint64_t indexes[9];
283 int n = findNeighborhood(_level, i, indexes);
284 return std::vector<uint64_t>(indexes, indexes + n);
285}
286
288 static char const FACE_NORM[6][2] = {
289 {'+', 'Z'}, {'+', 'X'}, {'+', 'Y'},
290 {'-', 'X'}, {'-', 'Y'}, {'-', 'Z'},
291 };
292 char s[MAX_LEVEL + 2];
293 if (i >= static_cast<uint64_t>(6) << (2 * _level)) {
294 throw std::invalid_argument("Invalid Q3C index");
295 }
296 // Print in base-4, from least to most significant digit.
297 char * p = s + (sizeof(s) - 1);
298 for (int l = _level; l > 0; --l, --p, i >>= 2) {
299 *p = '0' + (i & 3);
300 }
301 // The remaining bits correspond to the cube face.
302 --p;
303 p[0] = FACE_NORM[i][0];
304 p[1] = FACE_NORM[i][1];
305 return std::string(p, sizeof(s) - static_cast<size_t>(p - s));
306}
307
309 if (i >= static_cast<uint64_t>(6) << (2 * _level)) {
310 throw std::invalid_argument("Invalid Q3C index");
311 }
312 UnitVector3d verts[4];
313 makeQuad(i, _level, verts);
315 new ConvexPolygon(verts[0], verts[1], verts[2], verts[3]));
316}
317
318#if defined(NO_SIMD) || !defined(__x86_64__)
319 uint64_t Q3cPixelization::index(UnitVector3d const & p) const {
320 int face = faceNumber(p, FACE_NUM);
321 double w = std::fabs(p(FACE_COMP[face][2]));
322 double u = (p(FACE_COMP[face][0]) / w) * FACE_CONST[face][0];
323 double v = (p(FACE_COMP[face][1]) / w) * FACE_CONST[face][1];
324 std::tuple<int32_t, int32_t> g = faceToGrid(_level, u, v);
325 uint64_t z = mortonIndex(static_cast<uint32_t>(std::get<0>(g)),
326 static_cast<uint32_t>(std::get<1>(g)));
327 return (static_cast<uint64_t>(face) << (2 * _level)) | z;
328 }
329#else
330 uint64_t Q3cPixelization::index(UnitVector3d const & p) const {
331 int face = faceNumber(p, FACE_NUM);
332 __m128d ww = _mm_set1_pd(p(FACE_COMP[face][2]));
333 __m128d uv = _mm_set_pd(p(FACE_COMP[face][1]), p(FACE_COMP[face][0]));
334 uv = _mm_mul_pd(
335 _mm_div_pd(uv, _mm_andnot_pd(_mm_set_pd(-0.0, -0.0), ww)),
336 _mm_set_pd(FACE_CONST[face][1], FACE_CONST[face][0])
337 );
338 __m128i st = faceToGrid(_level, uv);
339 return (static_cast<uint64_t>(face) << (2 * _level)) | mortonIndex(st);
340 }
341#endif
342
343RangeSet Q3cPixelization::_envelope(Region const & r, size_t maxRanges) const {
344 return detail::findPixels<Q3cPixelFinder, false>(r, maxRanges, _level);
345}
346
347RangeSet Q3cPixelization::_interior(Region const & r, size_t maxRanges) const {
348 return detail::findPixels<Q3cPixelFinder, true>(r, maxRanges, _level);
349}
350
351}} // namespace lsst::sphgeom
This file declares a class for representing convex polygons with great circle edges on the unit spher...
afw::table::Key< afw::table::Array< MaskPixelT > > mask
double z
Definition Match.cc:44
This file provides a base class for pixel finders.
This file declares a Pixelization subclass for the Q3C indexing scheme.
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.
std::string toString(uint64_t i) const override
toString converts the given Q3C index to a human readable string.
RangeSet _interior(Region const &r, size_t maxRanges) const override
uint64_t index(UnitVector3d const &v) const override
index computes the index of the pixel for v.
std::vector< uint64_t > neighborhood(uint64_t i) const
neighborhood returns the indexes of all pixels that share a vertex with pixel i (including i itself).
Q3cPixelization(int level)
This constructor creates a Q3C pixelization of the sphere with the given subdivision level.
static constexpr int MAX_LEVEL
The maximum supported cube-face grid resolution is 2^30 by 2^30.
ConvexPolygon quad(uint64_t i) const
quad returns the quadrilateral corresponding to the Q3C pixel with index i.
std::unique_ptr< Region > pixel(uint64_t i) const override
pixel returns the spherical region corresponding to the pixel with index i.
RangeSet _envelope(Region const &r, size_t maxRanges) const override
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 > mortonIndexInverse(uint64_t z)
mortonIndexInverse separates the even and odd bits of z.
Definition curve.h:202
uint64_t mortonIndex(uint32_t x, uint32_t y)
mortonIndex interleaves the bits of x and y.
Definition curve.h:155
T sort(T... args)
double w
Definition CoaddPsf.cc:69
T tie(T... args)
T unique(T... args)