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
HtmPixelization.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 "lsst/sphgeom/curve.h"
30 
31 #include "PixelFinder.h"
32 
33 
34 namespace lsst {
35 namespace sphgeom {
36 
37 namespace {
38 
39 // `rootVertex` returns the i-th (0-3) root vertex of HTM root triangle r (0-8).
40 UnitVector3d const & rootVertex(int r, int i) {
41  static UnitVector3d const VERTICES[8][3] = {
50  };
51  return VERTICES[r][i];
52 }
53 
54 // `HtmPixelFinder` locates trixels that intersect a region.
55 template <typename RegionType, bool InteriorOnly>
56 class HtmPixelFinder: public detail::PixelFinder<
57  HtmPixelFinder<RegionType, InteriorOnly>, RegionType, InteriorOnly, 3>
58 {
59  using Base = detail::PixelFinder<
60  HtmPixelFinder<RegionType, InteriorOnly>, RegionType, InteriorOnly, 3>;
61  using Base::visit;
62 
63 public:
64  HtmPixelFinder(RangeSet & ranges,
65  RegionType const & region,
66  int level,
67  size_t maxRanges):
68  Base(ranges, region, level, maxRanges)
69  {}
70 
71  void operator()() {
72  UnitVector3d trixel[3];
73  // Loop over HTM root triangles.
74  for (uint64_t r = 0; r < 8; ++r) {
75  for (int v = 0; v < 3; ++v) {
76  trixel[v] = rootVertex(r, v);
77  }
78  visit(trixel, r + 8, 0);
79  }
80  }
81 
82  void subdivide(UnitVector3d const * trixel, uint64_t index, int level) {
83  UnitVector3d mid[3] = {
84  UnitVector3d(trixel[1] + trixel[2]),
85  UnitVector3d(trixel[2] + trixel[0]),
86  UnitVector3d(trixel[0] + trixel[1])
87  };
88  UnitVector3d child[3] = {trixel[0], mid[2], mid[1]};
89  index *= 4;
90  ++level;
91  visit(child, index, level);
92  child[0] = trixel[1];
93  child[1] = mid[0];
94  child[2] = mid[2];
95  ++index;
96  visit(child, index, level);
97  child[0] = trixel[2];
98  child[1] = mid[1];
99  child[2] = mid[0];
100  ++index;
101  visit(child, index, level);
102  ++index;
103  visit(mid, index, level);
104  }
105 };
106 
107 } // unnamed namespace
108 
109 
110 int HtmPixelization::level(uint64_t i) {
111  // An HTM index consists of 4 bits encoding the root triangle
112  // number (8 - 15), followed by 2l bits, where each of the l bit pairs
113  // encodes a child triangle number (0-3), and l is the desired level.
114  int j = log2(i);
115  // The level l is derivable from the index j of the MSB of i.
116  // For i to be valid, j must be an odd integer > 1.
117  if ((j & 1) == 0 || (j == 1)) {
118  return -1;
119  }
120  return (j - 3) >> 1;
121 }
122 
124  int l = level(i);
125  if (l < 0 || l > MAX_LEVEL) {
126  throw std::invalid_argument("Invalid HTM index");
127  }
128  l *= 2;
129  uint64_t r = (i >> l) & 7;
130  UnitVector3d v0 = rootVertex(r, 0);
131  UnitVector3d v1 = rootVertex(r, 1);
132  UnitVector3d v2 = rootVertex(r, 2);
133  for (l -= 2; l >= 0; l -= 2) {
134  int child = (i >> l) & 3;
135  UnitVector3d m12 = UnitVector3d(v1 + v2);
136  UnitVector3d m20 = UnitVector3d(v2 + v0);
137  UnitVector3d m01 = UnitVector3d(v0 + v1);
138  switch (child) {
139  case 0: v1 = m01; v2 = m20; break;
140  case 1: v0 = v1; v1 = m12; v2 = m01; break;
141  case 2: v0 = v2; v1 = m20; v2 = m12; break;
142  case 3: v0 = m12; v1 = m20; v2 = m01; break;
143  }
144  }
145  return ConvexPolygon(v0, v1, v2);
146 }
147 
149  char s[MAX_LEVEL + 2];
150  int l = level(i);
151  if (l < 0 || l > MAX_LEVEL) {
152  throw std::invalid_argument("Invalid HTM index");
153  }
154  // Print in base-4, from least to most significant digit.
155  char * p = s + (sizeof(s) - 1);
156  for (; l >= 0; --l, --p, i >>= 2) {
157  *p = '0' + (i & 3);
158  }
159  // The remaining bit corresponds to the hemisphere.
160  *p = (i & 1) == 0 ? 'S' : 'N';
161  return std::string(p, sizeof(s) - static_cast<size_t>(p - s));
162 }
163 
164 HtmPixelization::HtmPixelization(int level) : _level(level) {
165  if (level < 0 || level > MAX_LEVEL) {
166  throw std::invalid_argument("Invalid HTM subdivision level");
167  }
168 }
169 
170 uint64_t HtmPixelization::index(UnitVector3d const & v) const {
171  // Find the root triangle containing v.
172  uint64_t r;
173  if (v.z() < 0.0) {
174  // v is in the southern hemisphere (root triangle 0, 1, 2, or 3).
175  if (v.y() > 0.0) {
176  r = (v.x() > 0.0) ? 0 : 1;
177  } else if (v.y() == 0.0) {
178  r = (v.x() >= 0.0) ? 0 : 2;
179  } else {
180  r = (v.x() < 0.0) ? 2 : 3;
181  }
182  } else {
183  // v is in the northern hemisphere (root triangle 4, 5, 6, or 7).
184  if (v.y() > 0.0) {
185  r = (v.x() > 0.0) ? 7 : 6;
186  } else if (v.y() == 0.0) {
187  r = (v.x() >= 0.0) ? 7 : 5;
188  } else {
189  r = (v.x() < 0.0) ? 5 : 4;
190  }
191  }
192  UnitVector3d v0 = rootVertex(r, 0);
193  UnitVector3d v1 = rootVertex(r, 1);
194  UnitVector3d v2 = rootVertex(r, 2);
195  uint64_t i = r + 8;
196  for (int l = 0; l < _level; ++l) {
197  UnitVector3d m01 = UnitVector3d(v0 + v1);
198  UnitVector3d m20 = UnitVector3d(v2 + v0);
199  i <<= 2;
200  if (orientation(v, m01, m20) >= 0) {
201  v1 = m01; v2 = m20;
202  continue;
203  }
204  UnitVector3d m12 = UnitVector3d(v1 + v2);
205  if (orientation(v, m12, m01) >= 0) {
206  v0 = v1; v1 = m12; v2 = m01;
207  i += 1;
208  } else if (orientation(v, m20, m12) >= 0) {
209  v0 = v2; v1 = m20; v2 = m12;
210  i += 2;
211  } else {
212  v0 = m12; v1 = m20; v2 = m01;
213  i += 3;
214  }
215  }
216  return i;
217 }
218 
219 RangeSet HtmPixelization::_envelope(Region const & r, size_t maxRanges) const {
220  return detail::findPixels<HtmPixelFinder, false>(r, maxRanges, _level);
221 }
222 
223 RangeSet HtmPixelization::_interior(Region const & r, size_t maxRanges) const {
224  return detail::findPixels<HtmPixelFinder, true>(r, maxRanges, _level);
225 }
226 
227 }} // namespace lsst::sphgeom
static UnitVector3d X()
Definition: UnitVector3d.h:93
int orientation(UnitVector3d const &a, UnitVector3d const &b, UnitVector3d const &c)
orientation computes and returns the orientations of 3 unit vectors a, b and c.
Definition: orientation.cc:135
HtmPixelization(int level)
This constructor creates an HTM pixelization of the sphere with the given subdivision level...
static constexpr int MAX_LEVEL
MAX_LEVEL is the maximum supported HTM subdivision level.
static int level(uint64_t i)
level returns the subdivision level of the given HTM index.
static std::string asString(uint64_t i)
asString converts the given HTM index to a human readable string.
STL class.
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.
static UnitVector3d Y()
Definition: UnitVector3d.h:97
This file contains functions for space-filling curves.
ConvexPolygon is a closed convex polygon on the unit sphere.
Definition: ConvexPolygon.h:57
void visit(UnitVector3d const *pixel, uint64_t index, int level)
Definition: PixelFinder.h:81
static UnitVector3d Z()
Definition: UnitVector3d.h:101
A RangeSet is a set of unsigned 64 bit integers.
Definition: RangeSet.h:99
uint64_t index(UnitVector3d const &) const override
index computes the index of the pixel for v.
uint8_t log2(uint64_t x)
Definition: curve.h:98
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
This file declares functions for orienting points on the sphere.
static ConvexPolygon triangle(uint64_t i)
triangle returns the triangle corresponding to the given HTM index.
This file declares a Pixelization subclass for the HTM indexing scheme.