LSST Applications g042eb84c57+730a74494b,g04e9c324dd+8c5ae1fdc5,g134cb467dc+1f1e3e7524,g199a45376c+0ba108daf9,g1fd858c14a+fa7d31856b,g210f2d0738+f66ac109ec,g262e1987ae+83a3acc0e5,g29ae962dfc+d856a2cb1f,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+a1e0c9f713,g47891489e3+0d594cb711,g4d44eb3520+c57ec8f3ed,g4d7b6aa1c5+f66ac109ec,g53246c7159+8c5ae1fdc5,g56a1a4eaf3+fd7ad03fde,g64539dfbff+f66ac109ec,g67b6fd64d1+0d594cb711,g67fd3c3899+f66ac109ec,g6985122a63+0d594cb711,g74acd417e5+3098891321,g786e29fd12+668abc6043,g81db2e9a8d+98e2ab9f28,g87389fa792+8856018cbb,g89139ef638+0d594cb711,g8d7436a09f+80fda9ce03,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+a8a29bda4b,g99822b682c+e3ec3c61f9,g9d5c6a246b+0d5dac0c3d,ga41d0fce20+9243b26dd2,gbf99507273+8c5ae1fdc5,gd7ef33dd92+0d594cb711,gdab6d2f7ff+3098891321,ge410e46f29+0d594cb711,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.38
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 (std::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, std::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
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 std::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
173 throw std::invalid_argument("Invalid HTM subdivision level");
174 }
175}
176
178 // Find the root triangle containing v.
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 std::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, universe());
228}
229
230RangeSet HtmPixelization::_interior(Region const & r, size_t maxRanges) const {
231 return detail::findPixels<HtmPixelFinder, true>(r, maxRanges, _level, universe());
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
static ConvexPolygon triangle(std::uint64_t i)
triangle returns the triangle corresponding to the given HTM index.
RangeSet universe() const override
universe returns the set of all pixel indexes for this pixelization.
static std::string asString(std::uint64_t i)
asString converts the given HTM index to a human readable string.
static int level(std::uint64_t i)
level returns the subdivision level of the given HTM index.
std::uint64_t index(UnitVector3d const &) const override
index computes the index of the pixel for v.
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:89
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
static UnitVector3d Z()
static UnitVector3d Y()
static UnitVector3d X()
This file contains functions for space-filling curves.
RangeSet findPixels(Region const &r, size_t maxRanges, int level, RangeSet const &universe)
std::uint8_t log2(std::uint64_t x)
Definition curve.h:105
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.
This file declares functions for orienting points on the sphere.