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