Skip to content

fix(Geometry): add Hilbert sorting of sets of points #1142

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
MelchiorSchuh merged 1 commit into from
May 13, 2025
Merged

fix(Geometry): add Hilbert sorting of sets of points #1142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions include/geode/geometry/points_sort.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,4 +42,8 @@ namespace geode
template < index_t dimension >
[[nodiscard]] std::vector< index_t > morton_mapping(
absl::Span< const Point< dimension > > points );

template < index_t dimension >
[[nodiscard]] std::vector< index_t > hilbert_mapping(
absl::Span< const Point< dimension > > points );
} // namespace geode
95 changes: 71 additions & 24 deletions src/geode/geometry/points_sort.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,28 @@ namespace
};
ALIAS_1D_AND_2D_AND_3D( Morton_cmp );

template < geode::index_t dimension >
class Hilbert_cmp
{
public:
Hilbert_cmp( absl::Span< const geode::Point< dimension > > points,
geode::local_index_t coord )
: points_( points ), coord_( coord )
{
}

bool operator()( geode::index_t box1, geode::index_t box2 ) const
{
return points_[box1].value( coord_ )
> points_[box2].value( coord_ );
}

private:
absl::Span< const geode::Point< dimension > > points_;
geode::local_index_t coord_;
};
ALIAS_1D_AND_2D_AND_3D( Hilbert_cmp );

/**
* \brief Splits a sequence into two ordered halves.
* \details The algorithm shuffles the sequence and
Expand Down Expand Up @@ -88,7 +110,8 @@ namespace
* In CGAL User and Reference Manual. CGAL Editorial Board,
* 3.9 edition, 2011
*/
template < geode::local_index_t COORDX >
template < geode::local_index_t COORDX,
template < geode::index_t > typename Comparator >
void morton_mapping( absl::Span< const geode::Point3D > points,
const itr& begin,
const itr& end )
Expand All @@ -100,9 +123,9 @@ namespace
constexpr geode::local_index_t COORDY = COORDX == 2 ? 0 : COORDX + 1;
constexpr geode::local_index_t COORDZ = COORDY == 2 ? 0 : COORDY + 1;

const Morton_cmp3D compX{ points, COORDX };
const Morton_cmp3D compY{ points, COORDY };
const Morton_cmp3D compZ{ points, COORDZ };
const Comparator< 3 > compX{ points, COORDX };
const Comparator< 3 > compY{ points, COORDY };
const Comparator< 3 > compZ{ points, COORDZ };

const auto m0 = begin;
const auto m8 = end;
Expand All @@ -113,17 +136,18 @@ namespace
const auto m6 = split_container( m4, m8, compY );
const auto m5 = split_container( m4, m6, compZ );
const auto m7 = split_container( m6, m8, compZ );
morton_mapping< COORDZ >( points, m0, m1 );
morton_mapping< COORDY >( points, m1, m2 );
morton_mapping< COORDY >( points, m2, m3 );
morton_mapping< COORDX >( points, m3, m4 );
morton_mapping< COORDX >( points, m4, m5 );
morton_mapping< COORDY >( points, m5, m6 );
morton_mapping< COORDY >( points, m6, m7 );
morton_mapping< COORDZ >( points, m7, m8 );
morton_mapping< COORDZ, Comparator >( points, m0, m1 );
morton_mapping< COORDY, Comparator >( points, m1, m2 );
morton_mapping< COORDY, Comparator >( points, m2, m3 );
morton_mapping< COORDX, Comparator >( points, m3, m4 );
morton_mapping< COORDX, Comparator >( points, m4, m5 );
morton_mapping< COORDY, Comparator >( points, m5, m6 );
morton_mapping< COORDY, Comparator >( points, m6, m7 );
morton_mapping< COORDZ, Comparator >( points, m7, m8 );
}

template < geode::local_index_t COORDX >
template < geode::local_index_t COORDX,
template < geode::index_t > typename Comparator >
void morton_mapping( absl::Span< const geode::Point2D > points,
const itr& begin,
const itr& end )
Expand All @@ -134,21 +158,22 @@ namespace
}
constexpr geode::local_index_t COORDY = COORDX == 1 ? 0 : COORDX + 1;

const Morton_cmp2D compX{ points, COORDX };
const Morton_cmp2D compY{ points, COORDY };
const Comparator< 2 > compX{ points, COORDX };
const Comparator< 2 > compY{ points, COORDY };

const auto m0 = begin;
const auto m4 = end;
const auto m2 = split_container( m0, m4, compX );
const auto m1 = split_container( m0, m2, compY );
const auto m3 = split_container( m2, m4, compY );
morton_mapping< COORDY >( points, m0, m1 );
morton_mapping< COORDX >( points, m1, m2 );
morton_mapping< COORDX >( points, m2, m3 );
morton_mapping< COORDY >( points, m3, m4 );
morton_mapping< COORDY, Comparator >( points, m0, m1 );
morton_mapping< COORDX, Comparator >( points, m1, m2 );
morton_mapping< COORDX, Comparator >( points, m2, m3 );
morton_mapping< COORDY, Comparator >( points, m3, m4 );
}

template < geode::local_index_t COORDX >
template < geode::local_index_t COORDX,
template < geode::index_t > typename Comparator >
void morton_mapping( absl::Span< const geode::Point1D > points,
const itr& begin,
const itr& end )
Expand All @@ -157,13 +182,13 @@ namespace
{
return;
}
const Morton_cmp1D compX{ points, COORDX };
const Comparator< 1 > compX{ points, COORDX };

const auto m0 = begin;
const auto m2 = end;
const auto m1 = split_container( m0, m2, compX );
morton_mapping< COORDX >( points, m0, m1 );
morton_mapping< COORDX >( points, m1, m2 );
morton_mapping< COORDX, Comparator >( points, m0, m1 );
morton_mapping< COORDX, Comparator >( points, m1, m2 );
}

/*
Expand Down Expand Up @@ -214,7 +239,22 @@ namespace geode
[&mapping]( index_t i ) {
mapping[i] = i;
} );
::morton_mapping< 0_uc >( points, mapping.begin(), mapping.end() );
::morton_mapping< 0_uc, Morton_cmp >(
points, mapping.begin(), mapping.end() );
return mapping;
}

template < index_t dimension >
std::vector< index_t > hilbert_mapping(
absl::Span< const Point< dimension > > points )
{
std::vector< index_t > mapping( points.size() );
async::parallel_for( async::irange( size_t{ 0 }, mapping.size() ),
[&mapping]( index_t i ) {
mapping[i] = i;
} );
::morton_mapping< 0_uc, Hilbert_cmp >(
points, mapping.begin(), mapping.end() );
return mapping;
}

Expand All @@ -231,4 +271,11 @@ namespace geode
absl::Span< const Point< 2 > > );
template std::vector< index_t > opengeode_geometry_api morton_mapping(
absl::Span< const Point< 3 > > );

template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
absl::Span< const Point< 1 > > );
template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
absl::Span< const Point< 2 > > );
template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
absl::Span< const Point< 3 > > );
} // namespace geode