AABB.hpp

Go to the documentation of this file.

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//
 
#pragma once
 
#include "shambase/SourceLocation.hpp"
#include "shambase/assert.hpp"
#include "shambackends/math.hpp"
#include "shambackends/vec.hpp"
#include <limits>
 
namespace shammath {
 
    template<class T>
    struct Ray {
        using T_prop = shambase::VectorProperties<T>;
        using Tscal  = typename T_prop::component_type;
 
        T origin;
        T direction;
        T inv_direction;
 
        inline Ray(T origin, T direction)
            : origin(origin), direction(direction), inv_direction(Tscal{1.} / direction) {
 
            Tscal f = sycl::length(direction);
            SHAM_ASSERT(f > 0);
 
            this->direction /= f;
            this->inv_direction *= f;
        }
    };
 
    template<class T>
    struct RingRay {
        using T_prop = shambase::VectorProperties<T>;
        using Tscal  = typename T_prop::component_type;
 
        T center;
        Tscal radius;
        T e_x;
        T e_y;
 
        T get_ez() { return sycl::cross(e_x, e_y); }
        inline RingRay(T center, Tscal radius, T e_x, T e_y)
            : center(center), radius(radius), e_x(e_x), e_y(e_y) {}
    };
 
    template<class T>
    struct AABB {
 
        using T_prop = shambase::VectorProperties<T>;   
        using Tscal  = typename T_prop::component_type; 
 
        T lower; 
        T upper; 
 
        inline AABB() = default;
 
        inline AABB(T lower, T upper)
            : lower(lower), upper(upper) {
 
              };
 
        inline AABB(std::tuple<T, T> range)
            : lower(std::get<0>(range)), upper(std::get<1>(range)) {}
 
        inline AABB(std::pair<T, T> range) : lower(std::get<0>(range)), upper(std::get<1>(range)) {}
 
        inline T delt() const { return upper - lower; }
 
        inline Tscal get_radius() const { return sycl::length(delt()) / 2; }
 
        inline Tscal get_volume() { return sham::product_accumulate(upper - lower); }
 
        inline T get_center() const noexcept { return (lower + upper) / 2; }
 
        inline T sum_bounds() const noexcept { return lower + upper; }
 
        inline AABB expand_all(Tscal value) { return AABB{lower - value, upper + value}; }
 
        template<class Tb>
        inline AABB<Tb> convert() {
            using Tb_prop = shambase::VectorProperties<Tb>;
            static_assert(
                Tb_prop::dim == T_prop::dim, "you cannot change the dimension in convert");
 
            return {
                lower.template convert<Tb_prop::component_type>(),
                upper.template convert<Tb_prop::component_type>()};
        }
 
        inline AABB get_intersect(AABB other) const noexcept {
            return {sham::max(lower, other.lower), sham::min(upper, other.upper)};
        }
 
        inline bool contains(AABB other) const noexcept {
            // return lower <= other.lower && upper >= other.upper;
            return sham::vec_compare_leq(lower, other.lower)
                   && sham::vec_compare_geq(upper, other.upper);
        }
 
        inline bool contains_asymmetric(T point) const noexcept {
            return sham::vec_compare_leq(lower, point) && sham::vec_compare_g(upper, point);
        }
 
        [[nodiscard]] inline bool is_not_empty() const noexcept {
            return sham::vec_compare_geq(upper, lower);
        }
 
        [[nodiscard]] inline bool is_volume_not_null() const noexcept {
            return sham::vec_compare_g(upper, lower);
        }
 
        [[nodiscard]] inline bool is_surface() const noexcept {
            return sham::component_have_only_one_zero(delt()) && (is_not_empty());
        }
 
        [[nodiscard]] inline bool is_surface_or_volume() const noexcept {
            return sham::component_have_at_most_one_zero(delt()) && (is_not_empty());
        }
 
        [[nodiscard]] inline T clamp_coord(T coord) const noexcept {
            return sycl::clamp(coord, lower, upper);
        }
 
        [[nodiscard]] inline bool intersect_ray(Ray<T> ray) const noexcept;
 
        [[nodiscard]] inline bool intersect_ring_ray_approx(RingRay<T> ring_ray) const noexcept;
 
        inline bool operator==(const AABB<T> &other) const noexcept {
            return sham::equals(lower, other.lower) && sham::equals(upper, other.upper);
        }
 
        inline bool operator!=(const AABB<T> &other) const noexcept { return !(*this == other); }
    };
 
    template<class T>
    [[nodiscard]] inline bool AABB<T>::intersect_ray(Ray<T> ray) const noexcept {
        Tscal tmin = -shambase::get_infty<Tscal>(), tmax = shambase::get_infty<Tscal>();
 
        Tscal tx1 = (lower.x() - ray.origin.x()) * ray.inv_direction.x();
        Tscal tx2 = (upper.x() - ray.origin.x()) * ray.inv_direction.x();
 
        tmin = sycl::max(tmin, sycl::min(tx1, tx2));
        tmax = sycl::min(tmax, sycl::max(tx1, tx2));
 
        Tscal ty1 = (lower.y() - ray.origin.y()) * ray.inv_direction.y();
        Tscal ty2 = (upper.y() - ray.origin.y()) * ray.inv_direction.y();
 
        tmin = sycl::max(tmin, sycl::min(ty1, ty2));
        tmax = sycl::min(tmax, sycl::max(ty1, ty2));
 
        Tscal tz1 = (lower.z() - ray.origin.z()) * ray.inv_direction.z();
        Tscal tz2 = (upper.z() - ray.origin.z()) * ray.inv_direction.z();
 
        tmin = sycl::max(tmin, sycl::min(tz1, tz2));
        tmax = sycl::min(tmax, sycl::max(tz1, tz2));
 
        return tmax >= tmin;
    }
 
    template<class T>
    [[nodiscard]] inline bool AABB<T>::intersect_ring_ray_approx(
        RingRay<T> ring_ray) const noexcept {
        T aabb_center     = get_center();
        Tscal aabb_radius = get_radius();
 
        T r_center = ring_ray.center - aabb_center;
 
        Tscal x_val = sycl::dot(r_center, ring_ray.e_x);
        Tscal y_val = sycl::dot(r_center, ring_ray.e_y);
        Tscal z_val = sycl::dot(r_center, ring_ray.get_ez());
 
        Tscal r_val     = sycl::sqrt(x_val * x_val + y_val * y_val);
        Tscal delta_r   = r_val - ring_ray.radius;
        Tscal rab2_ring = z_val * z_val + delta_r * delta_r;
 
        return rab2_ring <= aabb_radius * aabb_radius;
    }
 
} // namespace shammath