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