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