Skip to content
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

Add ellipsoid geometry #1221

Merged
merged 7 commits into from
Mar 11, 2025
Merged

Add ellipsoid geometry #1221

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
0399d5f
Add EllipoidTag
aprokop Mar 8, 2025
d020222
Add Ellipsoid
aprokop Mar 8, 2025
288238e
Add centroid for Ellipsoid
aprokop Mar 8, 2025
a1c6690
Add intersects ellipsoid-point
aprokop Mar 8, 2025
a6109b2
Add intersects ellipsoid-segment
aprokop Mar 8, 2025
13de0b4
Add intersects ellipse-box
aprokop Mar 8, 2025
3667fb5
Try to work around MSVC
aprokop Mar 11, 2025
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
102 changes: 102 additions & 0 deletions src/geometry/ArborX_Ellipsoid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
/****************************************************************************
* Copyright (c) 2025, ArborX authors *
* All rights reserved. *
* *
* This file is part of the ArborX library. ArborX is *
* distributed under a BSD 3-clause license. For the licensing terms see *
* the LICENSE file in the top-level directory. *
* *
* SPDX-License-Identifier: BSD-3-Clause *
****************************************************************************/
#ifndef ARBORX_ELLIPSOID_HPP
#define ARBORX_ELLIPSOID_HPP

#include <ArborX_GeometryTraits.hpp>
#include <ArborX_Point.hpp>

#include <Kokkos_Assert.hpp>
#include <Kokkos_Macros.hpp>

#include <initializer_list>
#include <type_traits>

namespace ArborX::Experimental
{

template <int DIM, class Coordinate = float>
struct Ellipsoid
{
KOKKOS_DEFAULTED_FUNCTION
Ellipsoid() = default;

KOKKOS_FUNCTION
constexpr Ellipsoid(Point<DIM, Coordinate> const &centroid,
Coordinate const rmt[DIM][DIM])
: _centroid(centroid)
{
for (int i = 0; i < DIM; ++i)
for (int j = 0; j < DIM; ++j)
_rmt[i][j] = rmt[i][j];
}

KOKKOS_FUNCTION
constexpr Ellipsoid(
Point<DIM, Coordinate> const &centroid,
std::initializer_list<std::initializer_list<Coordinate>> const rmt)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we want this constructor?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For unit tests to allow {{a, b}, {b, c}} initialization.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could(should?) have used array of arrays

: _centroid(centroid)
{
KOKKOS_ASSERT(rmt.size() == DIM);
int i = 0;
for (auto const &row : rmt)
{
KOKKOS_ASSERT(row.size() == DIM);
int j = 0;
for (auto const &value : row)
_rmt[i][j++] = value;
++i;
}
}

KOKKOS_FUNCTION
constexpr auto &centroid() { return _centroid; }

KOKKOS_FUNCTION
constexpr auto const &centroid() const { return _centroid; }

KOKKOS_FUNCTION
constexpr auto const &rmt() const { return _rmt; }
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's RMT?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Riemann metric tensor

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment...


Point<DIM, Coordinate> _centroid = {};
Coordinate _rmt[DIM][DIM] = {};
};

template <typename T, std::size_t N>
#if KOKKOS_VERSION >= 40400
KOKKOS_DEDUCTION_GUIDE
#else
KOKKOS_FUNCTION
#endif
Ellipsoid(T const (&)[N], T const (&)[N][N]) -> Ellipsoid<N, T>;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we have this deduction guide and a constructor from an initializer lists instead of having a constructor from array of arrays?

Copy link
Contributor Author

@aprokop aprokop Mar 11, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was trying to make things work well enough. If you have a better approach that would simplify things, it would be nice. I want to be able to construct ellipsoid as

Ellipsoid ellipse{{3, 1}, {1, 3}}; // testing
float rmt[2][2];
// initialize rmt
Ellipsoid ellipse1{Point{0,0}, rmt};


} // namespace ArborX::Experimental

template <int DIM, class Coordinate>
struct ArborX::GeometryTraits::dimension<
ArborX::Experimental::Ellipsoid<DIM, Coordinate>>
{
static constexpr int value = DIM;
};
template <int DIM, class Coordinate>
struct ArborX::GeometryTraits::tag<
ArborX::Experimental::Ellipsoid<DIM, Coordinate>>
{
using type = EllipsoidTag;
};
template <int DIM, class Coordinate>
struct ArborX::GeometryTraits::coordinate_type<
ArborX::Experimental::Ellipsoid<DIM, Coordinate>>
{
using type = Coordinate;
};

#endif
3 changes: 2 additions & 1 deletion src/geometry/ArborX_GeometryTraits.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,14 +65,15 @@ DEFINE_GEOMETRY(kdop, KDOPTag);
DEFINE_GEOMETRY(tetrahedron, TetrahedronTag);
DEFINE_GEOMETRY(ray, RayTag);
DEFINE_GEOMETRY(segment, SegmentTag);
DEFINE_GEOMETRY(ellipsoid, EllipsoidTag);
#undef DEFINE_GEOMETRY

template <typename Geometry>
inline constexpr bool is_valid_geometry =
(is_point_v<Geometry> || is_box_v<Geometry> || is_sphere_v<Geometry> ||
is_kdop_v<Geometry> || is_triangle_v<Geometry> ||
is_tetrahedron_v<Geometry> || is_ray_v<Geometry> ||
is_segment_v<Geometry>);
is_segment_v<Geometry> || is_ellipsoid_v<Geometry>);

template <typename Geometry>
using DimensionNotSpecializedArchetypeAlias =
Expand Down
9 changes: 9 additions & 0 deletions src/geometry/algorithms/ArborX_Centroid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,15 @@ struct centroid<SegmentTag, Segment>
}
};

template <typename Ellipsoid>
struct centroid<EllipsoidTag, Ellipsoid>
{
KOKKOS_FUNCTION static auto apply(Ellipsoid const &ellipsoid)
{
return ellipsoid.centroid();
}
};

} // namespace Dispatch

} // namespace ArborX::Details
Expand Down
120 changes: 120 additions & 0 deletions src/geometry/algorithms/ArborX_Intersects.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,11 @@
#include "ArborX_Distance.hpp"
#include "ArborX_Expand.hpp"
#include <ArborX_GeometryTraits.hpp>
#include <ArborX_Segment.hpp>
#include <misc/ArborX_Vector.hpp>

#include <Kokkos_Array.hpp>
#include <Kokkos_Clamp.hpp>
#include <Kokkos_MathematicalFunctions.hpp>

namespace ArborX::Details
Expand Down Expand Up @@ -512,6 +515,123 @@ struct intersects<BoxTag, SegmentTag, Box, Segment>
}
};

namespace
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why anonymous space here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I try to use anonymous spaces when I know it's not needed anywhere outside of this file.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a header file though

{
// Computes x^t R y
template <int DIM, typename Coordinate>
KOKKOS_INLINE_FUNCTION auto rmt_multiply(Details::Vector<DIM, Coordinate> x,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any good reason not to spell out that it returns Coordinate instead of auto?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No good reason. I find either acceptable.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In that case it is slightly less readable

Coordinate const (&rmt)[DIM][DIM],
Details::Vector<DIM, Coordinate> y)
{
Coordinate r = 0;
for (int i = 0; i < DIM; ++i)
for (int j = 0; j < DIM; ++j)
r += x[i] * rmt[i][j] * y[j];
return r;
}
} // namespace

template <typename Ellipsoid, typename Point>
struct intersects<EllipsoidTag, PointTag, Ellipsoid, Point>
{
KOKKOS_FUNCTION static constexpr bool apply(Ellipsoid const &ellipsoid,
Point const &point)
{
auto d = point - ellipsoid.centroid();
return rmt_multiply(d, ellipsoid.rmt(), d) <= 1;
}
};

template <typename Point, typename Ellipsoid>
struct intersects<PointTag, EllipsoidTag, Point, Ellipsoid>
{
KOKKOS_FUNCTION static constexpr bool apply(Point const &point,
Ellipsoid const &ellipsoid)
{
return Details::intersects(ellipsoid, point);
}
};
Comment on lines +545 to +553
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't we generate intersects(geom2, geom1) from instersects(geom1, geom2)?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't. If you know of an easy way to do it, I'm all ears. Boost geometry does it, but I can't quite recall why we couldn't do it that way.


template <typename Ellipsoid, typename Segment>
struct intersects<EllipsoidTag, SegmentTag, Ellipsoid, Segment>
{
KOKKOS_FUNCTION static constexpr bool apply(Ellipsoid const &ellipsoid,
Segment const &segment)
{
// Preliminaries:
// - parametric segment formula: a + t(b-a)
// - ellipsoid formula: (x-c)^T R (x-c) <= 1
//
// Steps:
// - shift coordinates so that ellipsoid center is at origin
// new segment formula: a-c + t(b-a)
// new ellipsoid formula: x^T R x <= 1
// - substitute segment parametric equation into ellipsoid formula
// - find the value of t minimizing it
// t = -(R(b-a), a-c) / (R(b-a), b-a)
// - clamp t to [0, 1]
// - Plug the resulting point into ellipsoid equation
auto const &rmt = ellipsoid.rmt();

auto ab = segment.b - segment.a;
auto ca = segment.a - ellipsoid.centroid();

// At^2 + 2B^t + C
auto A = rmt_multiply(ab, rmt, ab);
auto B = rmt_multiply(ca, rmt, ab);
auto C = rmt_multiply(ca, rmt, ca);
auto t = -B / A;

using Float = coordinate_type_t<Segment>;
t = Kokkos::clamp(t, (Float)0, (Float)1);

return A * t * t + 2 * B * t + C <= 1;
}
};

template <typename Segment, typename Ellipsoid>
struct intersects<SegmentTag, EllipsoidTag, Segment, Ellipsoid>
{
KOKKOS_FUNCTION static constexpr bool apply(Segment const &segment,
Ellipsoid const &ellipsoid)
{
return Details::intersects(ellipsoid, segment);
}
};

template <typename Ellipsoid, typename Box>
struct intersects<EllipsoidTag, BoxTag, Ellipsoid, Box>
{
KOKKOS_FUNCTION static constexpr bool apply(Ellipsoid const &ellipsoid,
Box const &box)
{
static_assert(GeometryTraits::dimension_v<Box> == 2,
"Ellipsoid-box intersection is only implemented for 2D");

auto min_corner = box.minCorner();
auto max_corner = box.maxCorner();
using ::ArborX::Experimental::Segment;
// clang-format off
return
Details::intersects(ellipsoid.centroid(), box) ||
Details::intersects(ellipsoid, Segment{min_corner, {min_corner[0], max_corner[1]}}) ||
Details::intersects(ellipsoid, Segment{min_corner, {max_corner[0], min_corner[1]}}) ||
Details::intersects(ellipsoid, Segment{max_corner, {min_corner[0], max_corner[1]}}) ||
Details::intersects(ellipsoid, Segment{max_corner, {max_corner[0], min_corner[1]}});
// clang-format on
}
};

template <typename Box, typename Ellipsoid>
struct intersects<BoxTag, EllipsoidTag, Box, Ellipsoid>
{
KOKKOS_FUNCTION static constexpr bool apply(Box const &box,
Ellipsoid const &ellipsoid)
{
return Details::intersects(ellipsoid, box);
}
};

} // namespace Dispatch

} // namespace ArborX::Details
Expand Down
42 changes: 42 additions & 0 deletions test/tstDetailsAlgorithms.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
****************************************************************************/

#include <ArborX_Box.hpp>
#include <ArborX_Ellipsoid.hpp>
#include <ArborX_Point.hpp>
#include <ArborX_Segment.hpp>
#include <ArborX_Sphere.hpp>
Expand Down Expand Up @@ -341,6 +342,42 @@ BOOST_AUTO_TEST_CASE(intersects)
BOOST_TEST(!intersects(Segment2{{0.5, 1.6}, {2, 0}}, box2));
BOOST_TEST(intersects(box2, Segment2{{-1, 2}, {2, -1}}));
BOOST_TEST(!intersects(Segment2{{0.5, 1.6}, {2, 0}}, box2));

// ellipsoid [2x^2 - 3xy + 2y^2 <= 1]
using Ellipse = ArborX::Experimental::Ellipsoid<2>;
constexpr Ellipse ellipse{{1.f, 1.f}, {{2.f, -1.5f}, {-1.5f, 2.f}}};
BOOST_TEST(intersects(ellipse, Point2{1.f, 1.f}));
BOOST_TEST(intersects(ellipse, Point2{0.f, 0.f}));
BOOST_TEST(!intersects(ellipse, Point2{-0.01f, -0.01f}));
BOOST_TEST(intersects(ellipse, Point2{0.f, 0.5f}));
BOOST_TEST(!intersects(ellipse, Point2{0.f, 0.6f}));
BOOST_TEST(intersects(ellipse, Point2{2.f, 2.f}));
BOOST_TEST(!intersects(ellipse, Point2{0.5f, 1.5f}));
BOOST_TEST(intersects(ellipse, Point2{1.f, 0.3f}));
BOOST_TEST(!intersects(ellipse, Point2{1.f, 0.29f}));
BOOST_TEST(intersects(ellipse, Point2{1.f, 1.69f}));
BOOST_TEST(intersects(ellipse, Point2{1.f, 1.70f}));

BOOST_TEST(intersects(ellipse, Segment2{{-1, 1}, {1, -1}}));
BOOST_TEST(!intersects(ellipse, Segment2{{-1.1, 1}, {1, -1}}));
BOOST_TEST(intersects(ellipse, Segment2{{0, 0}, {0, 1}}));
BOOST_TEST(intersects(ellipse, Segment2{{0.5, 0.5}, {1.5, 1.5}}));
BOOST_TEST(intersects(ellipse, Segment2{{0.0, 1.9}, {3.0, 1.9}}));
BOOST_TEST(!intersects(ellipse, Segment2{{2.1, 0}, {2.1, 3}}));

using Box2 = ArborX::Box<2>;
BOOST_TEST(intersects(ellipse, Box2{{-10, -10}, {10, 10}}));
BOOST_TEST(intersects(ellipse, Box2{{0.5, 0.5}, {1.0, 1.0}}));
BOOST_TEST(intersects(ellipse, Box2{{-1, -1}, {0, 0}}));
BOOST_TEST(intersects(ellipse, Box2{{2, 2}, {3, 3}}));
BOOST_TEST(intersects(ellipse, Box2{{-1, -1}, {0, 2}}));
BOOST_TEST(intersects(ellipse, Box2{{-1, -1}, {2, 0}}));
BOOST_TEST(intersects(ellipse, Box2{{2, 1}, {3, 3}}));
BOOST_TEST(intersects(ellipse, Box2{{1, 2}, {3, 3}}));
BOOST_TEST(!intersects(ellipse, Box2{{1.5, 0}, {2, 0.5}}));
BOOST_TEST(!intersects(ellipse, Box2{{-1, -1}, {-0.1, -0.1}}));
BOOST_TEST(!intersects(ellipse, Box2{{0, 1.5}, {0.5, 2}}));
BOOST_TEST(!intersects(ellipse, Box2{{2.1, 2.1}, {3, 3}}));
}

BOOST_AUTO_TEST_CASE(equals)
Expand Down Expand Up @@ -456,6 +493,11 @@ BOOST_AUTO_TEST_CASE(centroid)
Segment segment{{-1.f, -1.f}, {3.f, 3.f}};
auto seg_center = returnCentroid(segment);
BOOST_TEST(equals(seg_center, ArborX::Point{1.f, 1.f}));

using ArborX::Experimental::Ellipsoid;
Ellipsoid ellipse{{1.f, 0.f}, {{2.f, 1.f}, {1.f, 2.f}}};
auto ell_center = returnCentroid(ellipse);
BOOST_TEST(equals(ell_center, ArborX::Point{1.f, 0.f}));
}

BOOST_AUTO_TEST_CASE(is_valid)
Expand Down