Skip to content

Commit 01d9632

Browse files
MelchiorSchuhMelchiorSchuh
authored
Merge pull request #1142 from Geode-solutions/fix/hilbert
fix(Geometry): add Hilbert sorting of sets of points
2 parents a1e8dd5 + e77754b commit 01d9632

File tree

2 files changed

+75
-24
lines changed

2 files changed

+75
-24
lines changed

‎include/geode/geometry/points_sort.hpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,4 +42,8 @@ namespace geode
4242
template < index_t dimension >
4343
[[nodiscard]] std::vector< index_t > morton_mapping(
4444
absl::Span< const Point< dimension > > points );
45+
46+
template < index_t dimension >
47+
[[nodiscard]] std::vector< index_t > hilbert_mapping(
48+
absl::Span< const Point< dimension > > points );
4549
} // namespace geode

‎src/geode/geometry/points_sort.cpp‎

Lines changed: 71 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,28 @@ namespace
5757
};
5858
ALIAS_1D_AND_2D_AND_3D( Morton_cmp );
5959

60+
template < geode::index_t dimension >
61+
class Hilbert_cmp
62+
{
63+
public:
64+
Hilbert_cmp( absl::Span< const geode::Point< dimension > > points,
65+
geode::local_index_t coord )
66+
: points_( points ), coord_( coord )
67+
{
68+
}
69+
70+
bool operator()( geode::index_t box1, geode::index_t box2 ) const
71+
{
72+
return points_[box1].value( coord_ )
73+
> points_[box2].value( coord_ );
74+
}
75+
76+
private:
77+
absl::Span< const geode::Point< dimension > > points_;
78+
geode::local_index_t coord_;
79+
};
80+
ALIAS_1D_AND_2D_AND_3D( Hilbert_cmp );
81+
6082
/**
6183
* \brief Splits a sequence into two ordered halves.
6284
* \details The algorithm shuffles the sequence and
@@ -88,7 +110,8 @@ namespace
88110
* In CGAL User and Reference Manual. CGAL Editorial Board,
89111
* 3.9 edition, 2011
90112
*/
91-
template < geode::local_index_t COORDX >
113+
template < geode::local_index_t COORDX,
114+
template < geode::index_t > typename Comparator >
92115
void morton_mapping( absl::Span< const geode::Point3D > points,
93116
const itr& begin,
94117
const itr& end )
@@ -100,9 +123,9 @@ namespace
100123
constexpr geode::local_index_t COORDY = COORDX == 2 ? 0 : COORDX + 1;
101124
constexpr geode::local_index_t COORDZ = COORDY == 2 ? 0 : COORDY + 1;
102125

103-
const Morton_cmp3D compX{ points, COORDX };
104-
const Morton_cmp3D compY{ points, COORDY };
105-
const Morton_cmp3D compZ{ points, COORDZ };
126+
const Comparator< 3 > compX{ points, COORDX };
127+
const Comparator< 3 > compY{ points, COORDY };
128+
const Comparator< 3 > compZ{ points, COORDZ };
106129

107130
const auto m0 = begin;
108131
const auto m8 = end;
@@ -113,17 +136,18 @@ namespace
113136
const auto m6 = split_container( m4, m8, compY );
114137
const auto m5 = split_container( m4, m6, compZ );
115138
const auto m7 = split_container( m6, m8, compZ );
116-
morton_mapping< COORDZ >( points, m0, m1 );
117-
morton_mapping< COORDY >( points, m1, m2 );
118-
morton_mapping< COORDY >( points, m2, m3 );
119-
morton_mapping< COORDX >( points, m3, m4 );
120-
morton_mapping< COORDX >( points, m4, m5 );
121-
morton_mapping< COORDY >( points, m5, m6 );
122-
morton_mapping< COORDY >( points, m6, m7 );
123-
morton_mapping< COORDZ >( points, m7, m8 );
139+
morton_mapping< COORDZ, Comparator >( points, m0, m1 );
140+
morton_mapping< COORDY, Comparator >( points, m1, m2 );
141+
morton_mapping< COORDY, Comparator >( points, m2, m3 );
142+
morton_mapping< COORDX, Comparator >( points, m3, m4 );
143+
morton_mapping< COORDX, Comparator >( points, m4, m5 );
144+
morton_mapping< COORDY, Comparator >( points, m5, m6 );
145+
morton_mapping< COORDY, Comparator >( points, m6, m7 );
146+
morton_mapping< COORDZ, Comparator >( points, m7, m8 );
124147
}
125148

126-
template < geode::local_index_t COORDX >
149+
template < geode::local_index_t COORDX,
150+
template < geode::index_t > typename Comparator >
127151
void morton_mapping( absl::Span< const geode::Point2D > points,
128152
const itr& begin,
129153
const itr& end )
@@ -134,21 +158,22 @@ namespace
134158
}
135159
constexpr geode::local_index_t COORDY = COORDX == 1 ? 0 : COORDX + 1;
136160

137-
const Morton_cmp2D compX{ points, COORDX };
138-
const Morton_cmp2D compY{ points, COORDY };
161+
const Comparator< 2 > compX{ points, COORDX };
162+
const Comparator< 2 > compY{ points, COORDY };
139163

140164
const auto m0 = begin;
141165
const auto m4 = end;
142166
const auto m2 = split_container( m0, m4, compX );
143167
const auto m1 = split_container( m0, m2, compY );
144168
const auto m3 = split_container( m2, m4, compY );
145-
morton_mapping< COORDY >( points, m0, m1 );
146-
morton_mapping< COORDX >( points, m1, m2 );
147-
morton_mapping< COORDX >( points, m2, m3 );
148-
morton_mapping< COORDY >( points, m3, m4 );
169+
morton_mapping< COORDY, Comparator >( points, m0, m1 );
170+
morton_mapping< COORDX, Comparator >( points, m1, m2 );
171+
morton_mapping< COORDX, Comparator >( points, m2, m3 );
172+
morton_mapping< COORDY, Comparator >( points, m3, m4 );
149173
}
150174

151-
template < geode::local_index_t COORDX >
175+
template < geode::local_index_t COORDX,
176+
template < geode::index_t > typename Comparator >
152177
void morton_mapping( absl::Span< const geode::Point1D > points,
153178
const itr& begin,
154179
const itr& end )
@@ -157,13 +182,13 @@ namespace
157182
{
158183
return;
159184
}
160-
const Morton_cmp1D compX{ points, COORDX };
185+
const Comparator< 1 > compX{ points, COORDX };
161186

162187
const auto m0 = begin;
163188
const auto m2 = end;
164189
const auto m1 = split_container( m0, m2, compX );
165-
morton_mapping< COORDX >( points, m0, m1 );
166-
morton_mapping< COORDX >( points, m1, m2 );
190+
morton_mapping< COORDX, Comparator >( points, m0, m1 );
191+
morton_mapping< COORDX, Comparator >( points, m1, m2 );
167192
}
168193

169194
/*
@@ -214,7 +239,22 @@ namespace geode
214239
[&mapping]( index_t i ) {
215240
mapping[i] = i;
216241
} );
217-
::morton_mapping< 0_uc >( points, mapping.begin(), mapping.end() );
242+
::morton_mapping< 0_uc, Morton_cmp >(
243+
points, mapping.begin(), mapping.end() );
244+
return mapping;
245+
}
246+
247+
template < index_t dimension >
248+
std::vector< index_t > hilbert_mapping(
249+
absl::Span< const Point< dimension > > points )
250+
{
251+
std::vector< index_t > mapping( points.size() );
252+
async::parallel_for( async::irange( size_t{ 0 }, mapping.size() ),
253+
[&mapping]( index_t i ) {
254+
mapping[i] = i;
255+
} );
256+
::morton_mapping< 0_uc, Hilbert_cmp >(
257+
points, mapping.begin(), mapping.end() );
218258
return mapping;
219259
}
220260

@@ -231,4 +271,11 @@ namespace geode
231271
absl::Span< const Point< 2 > > );
232272
template std::vector< index_t > opengeode_geometry_api morton_mapping(
233273
absl::Span< const Point< 3 > > );
274+
275+
template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
276+
absl::Span< const Point< 1 > > );
277+
template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
278+
absl::Span< const Point< 2 > > );
279+
template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
280+
absl::Span< const Point< 3 > > );
234281
} // namespace geode

0 commit comments

Comments
 (0)