Skip to content

Commit

Permalink
Update codebase to use new Vector
Browse files Browse the repository at this point in the history
  • Loading branch information
@aprokop
aprokop committed Apr 26, 2024
1 parent bcd7b52 commit 7770d93
Show file tree
Hide file tree
Showing 5 changed files with 61 additions and 115 deletions.
2 changes: 1 addition & 1 deletion examples/raytracing/example_raytracing.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -269,7 +269,7 @@ int main(int argc, char *argv[])
g, 2.f * Kokkos::numbers::pi_v<float>);
float theta =
acos(1 - 2 * Kokkos::rand<GeneratorType, float>::draw(g));
ArborX::Experimental::Vector direction{
typename ArborX::Experimental::Ray::Vector direction{
cos(upsilon) * sin(theta), sin(upsilon) * sin(theta), cos(theta)};

rays(j + i * rays_per_box) =
Expand Down
51 changes: 16 additions & 35 deletions src/geometry/ArborX_DetailsAlgorithms.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include <ArborX_Box.hpp>
#include <ArborX_DetailsKokkosExtArithmeticTraits.hpp>
#include <ArborX_DetailsKokkosExtMinMaxOperations.hpp> // min, max
#include <ArborX_DetailsVector.hpp>
#include <ArborX_GeometryTraits.hpp>

#include <Kokkos_Assert.hpp> // KOKKOS_ASSERT
Expand Down Expand Up @@ -252,24 +253,6 @@ struct distance<PointTag, TriangleTag, Point, Triangle>

static_assert(DIM == 2 || DIM == 3);

struct Vector : private Point
{
using Point::Point;
using Point::operator[];
KOKKOS_FUNCTION Vector(Point const &a, Point const &b)
{
for (int d = 0; d < DIM; ++d)
(*this)[d] = b[d] - a[d];
}
};

KOKKOS_FUNCTION static auto dot_product(Vector const &v, Vector const &w)
{
auto r = v[0] * w[0];
for (int d = 1; d < DIM; ++d)
r += v[d] * w[d];
return r;
}
KOKKOS_FUNCTION static auto combine(Point const &a, Point const &b,
Point const &c, Coordinate u,
Coordinate v, Coordinate w)
Expand All @@ -294,24 +277,24 @@ struct distance<PointTag, TriangleTag, Point, Triangle>
| |
*/

Vector ab(a, b);
Vector ac(a, c);
Vector ap(a, p);
auto const ab = b - a;
auto const ac = c - a;
auto const ap = p - a;

auto const d1 = dot_product(ab, ap);
auto const d2 = dot_product(ac, ap);
auto const d1 = ab.dot(ap);
auto const d2 = ac.dot(ap);
if (d1 <= 0 && d2 <= 0) // zone 1
return a;

Vector bp(b, p);
auto const d3 = dot_product(ab, bp);
auto const d4 = dot_product(ac, bp);
auto const bp = p - b;
auto const d3 = ab.dot(bp);
auto const d4 = ac.dot(bp);
if (d3 >= 0 && d4 <= d3) // zone 2
return b;

Vector cp(c, p);
auto const d5 = dot_product(ab, cp);
auto const d6 = dot_product(ac, cp);
auto const cp = p - c;
auto const d5 = ab.dot(cp);
auto const d6 = ac.dot(cp);
if (d6 >= 0 && d5 <= d6) // zone 3
return c;

Expand Down Expand Up @@ -610,18 +593,16 @@ struct intersects<BoxTag, TriangleTag, Box, Triangle>
}

using Point = decltype(a);
Point vector_ab{b[0] - a[0], b[1] - a[1], b[2] - a[2]};
Point vector_ac{c[0] - a[0], c[1] - a[1], c[2] - a[2]};
auto const vector_ab = b - a;
auto const vector_ac = c - a;
Point extents{(max_corner[0] - min_corner[0]) / 2,
(max_corner[1] - min_corner[1]) / 2,
(max_corner[2] - min_corner[2]) / 2};

// test normal of the triangle
// check if the projection of the triangle its normal lies in the interval
// defined by the projecting of the box onto the same vector
Point normal{{vector_ab[1] * vector_ac[2] - vector_ab[2] * vector_ac[1],
vector_ab[2] * vector_ac[0] - vector_ab[0] * vector_ac[2],
vector_ab[0] * vector_ac[1] - vector_ab[1] * vector_ac[0]}};
auto normal = vector_ab.cross(vector_ac);
auto radius = extents[0] * Kokkos::abs(normal[0]) +
extents[1] * Kokkos::abs(normal[1]) +
extents[2] * Kokkos::abs(normal[2]);
Expand All @@ -630,7 +611,7 @@ struct intersects<BoxTag, TriangleTag, Box, Triangle>
return false;

// Test crossproducts in a similar way as the triangle's normal above
Point vector_bc{c[0] - b[0], c[1] - b[1], c[2] - b[2]};
auto const vector_bc = c - b;

// e_x x vector_ab = (0, -vector_ab[2], vector_ab[1])
{
Expand Down
78 changes: 21 additions & 57 deletions src/geometry/ArborX_Ray.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include <ArborX_DetailsAlgorithms.hpp> // equal
#include <ArborX_DetailsKokkosExtArithmeticTraits.hpp>
#include <ArborX_DetailsKokkosExtSwap.hpp>
#include <ArborX_DetailsVector.hpp>
#include <ArborX_HyperTriangle.hpp>
#include <ArborX_Point.hpp>
#include <ArborX_Sphere.hpp>
Expand All @@ -27,54 +28,10 @@
namespace ArborX::Experimental
{

struct Vector : private Point
{
using Point::Point;
using Point::operator[];
friend KOKKOS_FUNCTION constexpr bool operator==(Vector const &v,
Vector const &w)
{
return v[0] == w[0] && v[1] == w[1] && v[2] == w[2];
}
};

template <typename Point1, typename Point2>
KOKKOS_INLINE_FUNCTION constexpr std::enable_if_t<
GeometryTraits::is_point<Point1>::value &&
GeometryTraits::is_point<Point2>::value &&
GeometryTraits::dimension_v<Point1> == 3 &&
GeometryTraits::dimension_v<Point2> == 3,
Vector>
makeVector(Point1 const &begin, Point2 const &end)
{
Vector v;
for (int d = 0; d < 3; ++d)
{
v[d] = end[d] - begin[d];
}
return v;
}

KOKKOS_INLINE_FUNCTION constexpr auto dotProduct(Vector const &v,
Vector const &w)
{
return v[0] * w[0] + v[1] * w[1] + v[2] * w[2];
}

KOKKOS_INLINE_FUNCTION constexpr Vector crossProduct(Vector const &v,
Vector const &w)
{
return {v[1] * w[2] - v[2] * w[1], v[2] * w[0] - v[0] * w[2],
v[0] * w[1] - v[1] * w[0]};
}

KOKKOS_INLINE_FUNCTION constexpr bool equals(Vector const &v, Vector const &w)
{
return v == w;
}

struct Ray
{
using Vector = ArborX::Details::Vector<3>;

Point _origin = {};
Vector _direction = {};

Expand Down Expand Up @@ -133,7 +90,7 @@ KOKKOS_INLINE_FUNCTION
constexpr bool equals(Ray const &l, Ray const &r)
{
using ArborX::Details::equals;
return equals(l.origin(), r.origin()) && equals(l.direction(), r.direction());
return equals(l.origin(), r.origin()) && l.direction() == r.direction();
}

KOKKOS_INLINE_FUNCTION
Expand Down Expand Up @@ -206,7 +163,7 @@ bool intersects(Ray const &ray, Box const &box)

// The function returns the index of the largest
// component of the direction vector.
KOKKOS_INLINE_FUNCTION int findLargestComp(Vector const &dir)
KOKKOS_INLINE_FUNCTION int findLargestComp(typename Ray::Vector const &dir)
{
int kz = 0;

Expand Down Expand Up @@ -303,6 +260,8 @@ bool intersection(Ray const &ray,
{
namespace KokkosExt = Details::KokkosExt;

using Vector = typename Ray::Vector;

auto dir = ray.direction();
// normalize the direction vector by its largest component.
auto kz = findLargestComp(dir);
Expand All @@ -319,15 +278,20 @@ bool intersection(Ray const &ray,
s[1] = dir[ky] * s[2];

// calculate vertices relative to ray origin
Vector const oA = makeVector(ray.origin(), triangle.a);
Vector const oB = makeVector(ray.origin(), triangle.b);
Vector const oC = makeVector(ray.origin(), triangle.c);
constexpr int dim = GeometryTraits::dimension_v<Point>;
using Float = GeometryTraits::coordinate_type_t<Point>;
auto const &o =
reinterpret_cast<ExperimentalHyperGeometry::Point<dim, Float> const &>(
ray.origin());
auto const oA = triangle.a - o;
auto const oB = triangle.b - o;
auto const oC = triangle.c - o;

// oA, oB, oB need to be normalized, otherwise they
// will scale with the problem size.
float const mag_oA = std::sqrt(dotProduct(oA, oA));
float const mag_oB = std::sqrt(dotProduct(oB, oB));
float const mag_oC = std::sqrt(dotProduct(oC, oC));
auto const mag_oA = oA.norm();
auto const mag_oB = oB.norm();
auto const mag_oC = oC.norm();

auto mag_bar = 3.0 / (mag_oA + mag_oB + mag_oC);

Expand Down Expand Up @@ -526,11 +490,11 @@ KOKKOS_INLINE_FUNCTION bool intersection(Ray const &ray, Sphere const &sphere,
auto const &r = sphere.radius();

// Vector oc = (origin_of_ray - center_of_sphere)
Vector const oc = makeVector(sphere.centroid(), ray.origin());
auto const oc = ray.origin() - sphere.centroid();

float const a2 = 1.f; // directions are normalized
float const a1 = 2.f * dotProduct(ray.direction(), oc);
float const a0 = dotProduct(oc, oc) - r * r;
float const a1 = 2.f * oc.dot(ray.direction());
float const a0 = oc.dot(oc) - r * r;

if (solveQuadratic(a2, a1, a0, tmin, tmax))
{
Expand Down
9 changes: 5 additions & 4 deletions test/tstQueryTreeRay.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@

#include "ArborXTest_StdVectorToKokkosView.hpp"
#include <ArborX.hpp>
#include <ArborX_DetailsVector.hpp>
#include <ArborX_Ray.hpp>

#include <boost/test/unit_test.hpp>
Expand Down Expand Up @@ -43,7 +44,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_ray_box_nearest, DeviceType,
ArborX::BVH<memory_space> const tree(exec_space, device_boxes);

ArborX::Experimental::Ray ray{ArborX::Point{0, 0, 0},
ArborX::Experimental::Vector{.15, .1, 0.}};
ArborX::Details::Vector{.15f, .1f, 0.f}};
Kokkos::View<ArborX::Experimental::Ray *, DeviceType> device_rays("rays", 1);
Kokkos::deep_copy(exec_space, device_rays, ray);

Expand Down Expand Up @@ -73,7 +74,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_ray_box_intersection, DeviceType,
ArborX::BVH<memory_space> const tree(exec_space, device_boxes);

ArborX::Experimental::Ray ray{ArborX::Point{0, 0, 0},
ArborX::Experimental::Vector{.1, .1, .1}};
ArborX::Details::Vector{.1f, .1f, .1f}};
Kokkos::View<ArborX::Experimental::Ray *, DeviceType> device_rays("rays", 1);
Kokkos::deep_copy(exec_space, device_rays, ray);

Expand Down Expand Up @@ -149,10 +150,10 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_ray_box_intersection_new, DeviceType,
2);
host_rays(0) = ArborX::Experimental::Ray{
ArborX::Point{0, 0, 0},
ArborX::Experimental::Vector{1.f / n, 1.f / n, 1.f / n}};
ArborX::Details::Vector{1.f / n, 1.f / n, 1.f / n}};
host_rays(1) = ArborX::Experimental::Ray{
ArborX::Point{n, n, n},
ArborX::Experimental::Vector{-1.f / n, -1.f / n, -1.f / n}};
ArborX::Details::Vector{-1.f / n, -1.f / n, -1.f / n}};
auto device_rays = Kokkos::create_mirror_view_and_copy(
typename DeviceType::memory_space{}, host_rays);
Kokkos::View<int *[2], DeviceType> device_ordered_intersections(
Expand Down
36 changes: 18 additions & 18 deletions test/tstRay.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -617,32 +617,32 @@ BOOST_AUTO_TEST_CASE(ray_triangle_intersection,

BOOST_AUTO_TEST_CASE(make_euclidean_vector)
{
using ArborX::Experimental::makeVector;
using ArborX::Experimental::Vector;
using Vector = ArborX::Details::Vector<3>;
using ArborX::ExperimentalHyperGeometry::Point;
static_assert(makeVector(Point{0, 0, 0}, Point{1, 2, 3}) == Vector{1, 2, 3});
static_assert(makeVector(Point{1, 2, 3}, Point{4, 5, 6}) == Vector{3, 3, 3});
static_assert(Point{1, 2, 3} - Point{0, 0, 0} == Vector{1, 2, 3});
static_assert(Point{4, 5, 6} - Point{1, 2, 3} == Vector{3, 3, 3});
}

BOOST_AUTO_TEST_CASE(dot_product)
{
using ArborX::Experimental::dotProduct;
static_assert(dotProduct({1, 0, 0}, {1, 0, 0}) == 1);
static_assert(dotProduct({1, 0, 0}, {0, 1, 0}) == 0);
static_assert(dotProduct({1, 0, 0}, {0, 0, 1}) == 0);
static_assert(dotProduct({1, 1, 1}, {1, 1, 1}) == 3);
using ArborX::Details::Vector;
static_assert(Vector{1, 0, 0}.dot(Vector{1, 0, 0}) == 1);
static_assert(Vector{1, 0, 0}.dot(Vector{0, 1, 0}) == 0);
static_assert(Vector{1, 0, 0}.dot(Vector{0, 0, 1}) == 0);
static_assert(Vector{1, 1, 1}.dot(Vector{1, 1, 1}) == 3);
}

BOOST_AUTO_TEST_CASE(cross_product)
{
using ArborX::Experimental::crossProduct;
using ArborX::Experimental::Vector;
static_assert(crossProduct({1, 0, 0}, {1, 0, 0}) == Vector{0, 0, 0});
static_assert(crossProduct({1, 0, 0}, {0, 1, 0}) == Vector{0, 0, 1});
static_assert(crossProduct({1, 0, 0}, {0, 0, 1}) == Vector{0, -1, 0});
static_assert(crossProduct({0, 1, 0}, {1, 0, 0}) == Vector{0, 0, -1});
static_assert(crossProduct({0, 1, 0}, {0, 1, 0}) == Vector{0, 0, 0});
static_assert(crossProduct({0, 1, 0}, {0, 0, 1}) == Vector{1, 0, 0});
static_assert(crossProduct({1, 1, 1}, {1, 1, 1}) == Vector{0, 0, 0});
using ArborX::Details::Vector;
// clang-format off
static_assert(Vector{1, 0, 0}.cross(Vector{1, 0, 0}) == Vector{0, 0, 0});
static_assert(Vector{1, 0, 0}.cross(Vector{0, 1, 0}) == Vector{0, 0, 1});
static_assert(Vector{1, 0, 0}.cross(Vector{0, 0, 1}) == Vector{0, -1, 0});
static_assert(Vector{0, 1, 0}.cross(Vector{1, 0, 0}) == Vector{0, 0, -1});
static_assert(Vector{0, 1, 0}.cross(Vector{0, 1, 0}) == Vector{0, 0, 0});
static_assert(Vector{0, 1, 0}.cross(Vector{0, 0, 1}) == Vector{1, 0, 0});
static_assert(Vector{1, 1, 1}.cross(Vector{1, 1, 1}) == Vector{0, 0, 0});
// clang-format on
}
#undef static_assert

0 comments on commit 7770d93

Please sign in to comment.