diff --git a/nanovdb/nanovdb/NanoVDB.h b/nanovdb/nanovdb/NanoVDB.h index 78955cd251..b1a661491d 100644 --- a/nanovdb/nanovdb/NanoVDB.h +++ b/nanovdb/nanovdb/NanoVDB.h @@ -1406,7 +1406,7 @@ struct Map double mTaperD; // 8B, placeholder for taper value /// @brief Default constructor for the identity map - __hostdev__ Map() + __hostdev__ constexpr Map() : mMatF{ 1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f} , mInvMatF{1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f} , mVecF{0.0f, 0.0f, 0.0f} @@ -1417,7 +1417,7 @@ struct Map , mTaperD{1.0} { } - __hostdev__ Map(double s, const Vec3d& t = Vec3d(0.0, 0.0, 0.0)) + __hostdev__ constexpr Map(double s, const Vec3d& t = Vec3d(0.0, 0.0, 0.0)) : mMatF{float(s), 0.0f, 0.0f, 0.0f, float(s), 0.0f, 0.0f, 0.0f, float(s)} , mInvMatF{1.0f / float(s), 0.0f, 0.0f, 0.0f, 1.0f / float(s), 0.0f, 0.0f, 0.0f, 1.0f / float(s)} , mVecF{float(t[0]), float(t[1]), float(t[2])} @@ -1449,7 +1449,7 @@ struct Map /// @param ijk 3D vector to be mapped - typically floating point index coordinates /// @return Forward mapping for affine transformation, i.e. (mat x ijk) + translation template - __hostdev__ Vec3T applyMap(const Vec3T& ijk) const { return math::matMult(mMatD, mVecD, ijk); } + __hostdev__ constexpr Vec3T applyMap(const Vec3T& ijk) const { return math::matMult(mMatD, mVecD, ijk); } /// @brief Apply the forward affine transformation to a vector using 32bit floating point arithmetics. /// @note Typically this operation is used for the scale, rotation and translation of index -> world mapping @@ -1457,7 +1457,7 @@ struct Map /// @param ijk 3D vector to be mapped - typically floating point index coordinates /// @return Forward mapping for affine transformation, i.e. (mat x ijk) + translation template - __hostdev__ Vec3T applyMapF(const Vec3T& ijk) const { return math::matMult(mMatF, mVecF, ijk); } + __hostdev__ constexpr Vec3T applyMapF(const Vec3T& ijk) const { return math::matMult(mMatF, mVecF, ijk); } /// @brief Apply the linear forward 3x3 transformation to an input 3d vector using 64bit floating point arithmetics, /// e.g. scale and rotation WITHOUT translation. @@ -1466,7 +1466,7 @@ struct Map /// @param ijk 3D vector to be mapped - typically floating point index coordinates /// @return linear forward 3x3 mapping of the input vector template - __hostdev__ Vec3T applyJacobian(const Vec3T& ijk) const { return math::matMult(mMatD, ijk); } + __hostdev__ constexpr Vec3T applyJacobian(const Vec3T& ijk) const { return math::matMult(mMatD, ijk); } /// @brief Apply the linear forward 3x3 transformation to an input 3d vector using 32bit floating point arithmetics, /// e.g. scale and rotation WITHOUT translation. @@ -1475,7 +1475,7 @@ struct Map /// @param ijk 3D vector to be mapped - typically floating point index coordinates /// @return linear forward 3x3 mapping of the input vector template - __hostdev__ Vec3T applyJacobianF(const Vec3T& ijk) const { return math::matMult(mMatF, ijk); } + __hostdev__ constexpr Vec3T applyJacobianF(const Vec3T& ijk) const { return math::matMult(mMatF, ijk); } /// @brief Apply the inverse affine mapping to a vector using 64bit floating point arithmetics. /// @note Typically this operation is used for the world -> index mapping @@ -1483,7 +1483,7 @@ struct Map /// @param xyz 3D vector to be mapped - typically floating point world coordinates /// @return Inverse affine mapping of the input @c xyz i.e. (xyz - translation) x mat^-1 template - __hostdev__ Vec3T applyInverseMap(const Vec3T& xyz) const + __hostdev__ constexpr Vec3T applyInverseMap(const Vec3T& xyz) const { return math::matMult(mInvMatD, Vec3T(xyz[0] - mVecD[0], xyz[1] - mVecD[1], xyz[2] - mVecD[2])); } @@ -1494,7 +1494,7 @@ struct Map /// @param xyz 3D vector to be mapped - typically floating point world coordinates /// @return Inverse affine mapping of the input @c xyz i.e. (xyz - translation) x mat^-1 template - __hostdev__ Vec3T applyInverseMapF(const Vec3T& xyz) const + __hostdev__ constexpr Vec3T applyInverseMapF(const Vec3T& xyz) const { return math::matMult(mInvMatF, Vec3T(xyz[0] - mVecF[0], xyz[1] - mVecF[1], xyz[2] - mVecF[2])); } @@ -1506,7 +1506,7 @@ struct Map /// @param ijk 3D vector to be mapped - typically floating point index coordinates /// @return linear inverse 3x3 mapping of the input vector i.e. xyz x mat^-1 template - __hostdev__ Vec3T applyInverseJacobian(const Vec3T& xyz) const { return math::matMult(mInvMatD, xyz); } + __hostdev__ constexpr Vec3T applyInverseJacobian(const Vec3T& xyz) const { return math::matMult(mInvMatD, xyz); } /// @brief Apply the linear inverse 3x3 transformation to an input 3d vector using 32bit floating point arithmetics, /// e.g. inverse scale and inverse rotation WITHOUT translation. @@ -1515,7 +1515,7 @@ struct Map /// @param ijk 3D vector to be mapped - typically floating point index coordinates /// @return linear inverse 3x3 mapping of the input vector i.e. xyz x mat^-1 template - __hostdev__ Vec3T applyInverseJacobianF(const Vec3T& xyz) const { return math::matMult(mInvMatF, xyz); } + __hostdev__ constexpr Vec3T applyInverseJacobianF(const Vec3T& xyz) const { return math::matMult(mInvMatF, xyz); } /// @brief Apply the transposed inverse 3x3 transformation to an input 3d vector using 64bit floating point arithmetics, /// e.g. inverse scale and inverse rotation WITHOUT translation. @@ -1524,12 +1524,12 @@ struct Map /// @param ijk 3D vector to be mapped - typically floating point index coordinates /// @return linear inverse 3x3 mapping of the input vector i.e. xyz x mat^-1 template - __hostdev__ Vec3T applyIJT(const Vec3T& xyz) const { return math::matMultT(mInvMatD, xyz); } + __hostdev__ constexpr Vec3T applyIJT(const Vec3T& xyz) const { return math::matMultT(mInvMatD, xyz); } template - __hostdev__ Vec3T applyIJTF(const Vec3T& xyz) const { return math::matMultT(mInvMatF, xyz); } + __hostdev__ constexpr Vec3T applyIJTF(const Vec3T& xyz) const { return math::matMultT(mInvMatF, xyz); } /// @brief Return a voxels size in each coordinate direction, measured at the origin - __hostdev__ Vec3d getVoxelSize() const { return this->applyMap(Vec3d(1)) - this->applyMap(Vec3d(0)); } + __hostdev__ constexpr Vec3d getVoxelSize() const { return this->applyMap(Vec3d(1)) - this->applyMap(Vec3d(0)); } }; // Map template diff --git a/nanovdb/nanovdb/math/Math.h b/nanovdb/nanovdb/math/Math.h index 19b1beee81..368ffbee43 100644 --- a/nanovdb/nanovdb/math/Math.h +++ b/nanovdb/nanovdb/math/Math.h @@ -17,6 +17,8 @@ #include // for __hostdev__ and lots of other utility functions +#include // for std::is_floating_point (matches long double in addition to float / double) + #if defined(__CUDA_ARCH__) #include // for ::cuda::std::numeric_limits #endif @@ -30,255 +32,316 @@ namespace math {// ============================================================= //@{ /// @brief Pi constant taken from Boost to match old behaviour template -inline __hostdev__ constexpr T pi() +inline __hostdev__ constexpr T pi() noexcept { return 3.141592653589793238462643383279502884e+00; } template<> -inline __hostdev__ constexpr float pi() +inline __hostdev__ constexpr float pi() noexcept { return 3.141592653589793238462643383279502884e+00F; } template<> -inline __hostdev__ constexpr double pi() +inline __hostdev__ constexpr double pi() noexcept { return 3.141592653589793238462643383279502884e+00; } template<> -inline __hostdev__ constexpr long double pi() +inline __hostdev__ constexpr long double pi() noexcept { return 3.141592653589793238462643383279502884e+00L; } //@} //@{ -/// Tolerance for floating-point comparison +/// @brief Per-type tolerance used by approximate floating-point comparisons. +/// @details Returned by @c Tolerance::value() and consumed by @c isApproxZero. +/// Specialized for @c float, @c double, and @c long double. Instantiating on +/// any other @c T is an "incomplete type" compile error by design — callers +/// must opt in by adding a specialization. template struct Tolerance; template<> struct Tolerance { - __hostdev__ static float value() { return 1e-8f; } + __hostdev__ [[nodiscard]] static constexpr float value() noexcept { return 1e-8f; } }; template<> struct Tolerance { - __hostdev__ static double value() { return 1e-15; } + __hostdev__ [[nodiscard]] static constexpr double value() noexcept { return 1e-15; } +}; +template<> +struct Tolerance +{ + __hostdev__ [[nodiscard]] static constexpr long double value() noexcept { return 1e-18L; } }; //@} //@{ -/// Delta for small floating-point offsets +/// @brief Per-type small floating-point offset, useful for nudging boundaries +/// off-edge in geometric tests. +/// @details Returned by @c Delta::value(). Specialized for @c float, +/// @c double, and @c long double; instantiating on any other @c T is a +/// compile error by design. template struct Delta; template<> struct Delta { - __hostdev__ static float value() { return 1e-5f; } + __hostdev__ [[nodiscard]] static constexpr float value() noexcept { return 1e-5f; } }; template<> struct Delta { - __hostdev__ static double value() { return 1e-9; } + __hostdev__ [[nodiscard]] static constexpr double value() noexcept { return 1e-9; } +}; +template<> +struct Delta +{ + __hostdev__ [[nodiscard]] static constexpr long double value() noexcept { return 1e-12L; } }; //@} //@{ -/// Maximum floating-point values +/// @brief Per-type largest representable value (mirrors @c std::numeric_limits::max). +/// @details Returned by @c Maximum::value(). The host fallback (`#else` +/// branch) forwards to @c std::numeric_limits so any @c T with a +/// specialization there works; the CUDA-device branch uses +/// @c ::cuda::std::numeric_limits; the HIP branch hard-codes a handful of +/// common types because @c std::numeric_limits is unavailable on the device +/// in some HIP toolchains. template struct Maximum; #if defined(__CUDA_ARCH__) template struct Maximum { - __hostdev__ static T value() { return ::cuda::std::numeric_limits::max(); } + __hostdev__ [[nodiscard]] static constexpr T value() noexcept { return ::cuda::std::numeric_limits::max(); } }; #elif defined(__HIP__) template<> struct Maximum { - __hostdev__ static int value() { return 2147483647; } + __hostdev__ [[nodiscard]] static constexpr int value() noexcept { return 2147483647; } }; template<> struct Maximum { - __hostdev__ static uint32_t value() { return 4294967295u; } + __hostdev__ [[nodiscard]] static constexpr uint32_t value() noexcept { return 4294967295u; } }; template<> struct Maximum { - __hostdev__ static float value() { return 1e+38f; } + __hostdev__ [[nodiscard]] static constexpr float value() noexcept { return 1e+38f; } }; template<> struct Maximum { - __hostdev__ static double value() { return 1e+308; } + __hostdev__ [[nodiscard]] static constexpr double value() noexcept { return 1e+308; } }; #else template struct Maximum { - static T value() { return std::numeric_limits::max(); } + [[nodiscard]] static constexpr T value() noexcept { return std::numeric_limits::max(); } }; #endif //@} +/// @brief Return @c true if @c |x| is below @c Tolerance::value(). template -__hostdev__ inline bool isApproxZero(const Type& x) +__hostdev__ [[nodiscard]] inline constexpr bool isApproxZero(const Type& x) noexcept { return !(x > Tolerance::value()) && !(x < -Tolerance::value()); } +/// @brief Same-type minimum (the primary template). template -__hostdev__ inline Type Min(Type a, Type b) +__hostdev__ [[nodiscard]] inline constexpr Type Min(Type a, Type b) noexcept { return (a < b) ? a : b; } -__hostdev__ inline int32_t Min(int32_t a, int32_t b) +// Mixed integer-type tiebreakers: the templated `Min` requires both +// args to deduce to the same `Type`, so calls like `Min(uint32_t, size_t)` +// would otherwise fall through to the float/double overloads with equal +// implicit-conversion cost — ambiguous. These integer overloads win. +__hostdev__ [[nodiscard]] inline constexpr int32_t Min(int32_t a, int32_t b) noexcept { - return int32_t(fminf(float(a), float(b))); + return a < b ? a : b; } -__hostdev__ inline uint32_t Min(uint32_t a, uint32_t b) +__hostdev__ [[nodiscard]] inline constexpr uint32_t Min(uint32_t a, uint32_t b) noexcept { - return uint32_t(fminf(float(a), float(b))); + return a < b ? a : b; } -__hostdev__ inline float Min(float a, float b) +__hostdev__ [[nodiscard]] inline float Min(float a, float b) noexcept { return fminf(a, b); } -__hostdev__ inline double Min(double a, double b) +__hostdev__ [[nodiscard]] inline double Min(double a, double b) noexcept { return fmin(a, b); } +/// @brief Same-type maximum (the primary template). template -__hostdev__ inline Type Max(Type a, Type b) +__hostdev__ [[nodiscard]] inline constexpr Type Max(Type a, Type b) noexcept { return (a > b) ? a : b; } - -__hostdev__ inline int32_t Max(int32_t a, int32_t b) +// See Min above: these integer overloads disambiguate mixed-int calls. +__hostdev__ [[nodiscard]] inline constexpr int32_t Max(int32_t a, int32_t b) noexcept { - return int32_t(fmaxf(float(a), float(b))); + return a > b ? a : b; } -__hostdev__ inline uint32_t Max(uint32_t a, uint32_t b) +__hostdev__ [[nodiscard]] inline constexpr uint32_t Max(uint32_t a, uint32_t b) noexcept { - return uint32_t(fmaxf(float(a), float(b))); + return a > b ? a : b; } -__hostdev__ inline float Max(float a, float b) +__hostdev__ [[nodiscard]] inline float Max(float a, float b) noexcept { return fmaxf(a, b); } -__hostdev__ inline double Max(double a, double b) +__hostdev__ [[nodiscard]] inline double Max(double a, double b) noexcept { return fmax(a, b); } -__hostdev__ inline float Clamp(float x, float a, float b) +//@{ +/// @brief Clamp @c x to the closed interval [@c a, @c b]. +__hostdev__ [[nodiscard]] inline float Clamp(float x, float a, float b) noexcept { return Max(Min(x, b), a); } -__hostdev__ inline double Clamp(double x, double a, double b) +__hostdev__ [[nodiscard]] inline double Clamp(double x, double a, double b) noexcept { return Max(Min(x, b), a); } +//@} -__hostdev__ inline float Fract(float x) +//@{ +/// @brief Fractional part of @c x, i.e. @c x @c - @c floor(x). Always non-negative. +__hostdev__ [[nodiscard]] inline float Fract(float x) noexcept { return x - floorf(x); } -__hostdev__ inline double Fract(double x) +__hostdev__ [[nodiscard]] inline double Fract(double x) noexcept { return x - floor(x); } +//@} -__hostdev__ inline int32_t Floor(float x) +//@{ +/// @brief Floor of @c x as a 32-bit signed integer (truncation toward -inf). +__hostdev__ [[nodiscard]] inline int32_t Floor(float x) noexcept { return int32_t(floorf(x)); } -__hostdev__ inline int32_t Floor(double x) +__hostdev__ [[nodiscard]] inline int32_t Floor(double x) noexcept { return int32_t(floor(x)); } +//@} -__hostdev__ inline int32_t Ceil(float x) +//@{ +/// @brief Ceil of @c x as a 32-bit signed integer (rounding toward +inf). +__hostdev__ [[nodiscard]] inline int32_t Ceil(float x) noexcept { return int32_t(ceilf(x)); } -__hostdev__ inline int32_t Ceil(double x) +__hostdev__ [[nodiscard]] inline int32_t Ceil(double x) noexcept { return int32_t(ceil(x)); } +//@} +/// @brief Return @c x * @c x. template -__hostdev__ inline T Pow2(T x) +__hostdev__ [[nodiscard]] inline constexpr T Pow2(T x) noexcept { return x * x; } +/// @brief Return @c x * @c x * @c x. template -__hostdev__ inline T Pow3(T x) +__hostdev__ [[nodiscard]] inline constexpr T Pow3(T x) noexcept { return x * x * x; } +/// @brief Return @c x raised to the fourth power. template -__hostdev__ inline T Pow4(T x) +__hostdev__ [[nodiscard]] inline constexpr T Pow4(T x) noexcept { return Pow2(x * x); } +/// @brief Absolute value (generic primary template, branchless for arithmetic types). +/// @note The @c float / @c double / @c int specializations below are not template -__hostdev__ inline T Abs(T x) +__hostdev__ [[nodiscard]] inline constexpr T Abs(T x) noexcept { return x < 0 ? -x : x; } template<> -__hostdev__ inline float Abs(float x) +__hostdev__ [[nodiscard]] inline float Abs(float x) noexcept { return fabsf(x); } template<> -__hostdev__ inline double Abs(double x) +__hostdev__ [[nodiscard]] inline double Abs(double x) noexcept { return fabs(x); } template<> -__hostdev__ inline int Abs(int x) +__hostdev__ [[nodiscard]] inline int Abs(int x) noexcept { return abs(x); } +/// @brief Round each component of @c xyz to its closest integer coordinate. +/// @details Forward declaration of the primary template — there is no +/// definition here. Callers resolve to one of the @c float / @c double +/// overloads below, both of which use the @c floor(x+0.5) rule +/// (round-half-toward-+inf), so a @c float and @c double input with the +/// same value yield the same integer coordinate. template class Vec3T> -__hostdev__ inline CoordT Round(const Vec3T& xyz); +__hostdev__ [[nodiscard]] inline CoordT Round(const Vec3T& xyz) noexcept; +/// @brief Round each component to its closest integer (round-half-toward-+inf) +/// using @c floor(x+0.5). Same rule is applied to both single and double +/// precision so float and double inputs yield the same integer coords. template class Vec3T> -__hostdev__ inline CoordT Round(const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline CoordT Round(const Vec3T& xyz) noexcept { - return CoordT(int32_t(rintf(xyz[0])), int32_t(rintf(xyz[1])), int32_t(rintf(xyz[2]))); - //return CoordT(int32_t(roundf(xyz[0])), int32_t(roundf(xyz[1])), int32_t(roundf(xyz[2])) ); - //return CoordT(int32_t(floorf(xyz[0] + 0.5f)), int32_t(floorf(xyz[1] + 0.5f)), int32_t(floorf(xyz[2] + 0.5f))); + return CoordT(int32_t(floorf(xyz[0] + 0.5f)), + int32_t(floorf(xyz[1] + 0.5f)), + int32_t(floorf(xyz[2] + 0.5f))); } +/// @brief Double-precision variant of @c Round — see the @c float overload above. template class Vec3T> -__hostdev__ inline CoordT Round(const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline CoordT Round(const Vec3T& xyz) noexcept { - return CoordT(int32_t(floor(xyz[0] + 0.5)), int32_t(floor(xyz[1] + 0.5)), int32_t(floor(xyz[2] + 0.5))); + return CoordT(int32_t(floor(xyz[0] + 0.5)), + int32_t(floor(xyz[1] + 0.5)), + int32_t(floor(xyz[2] + 0.5))); } +/// @brief Round each component of @c xyz down (toward -inf) into a @c CoordT. template class Vec3T> -__hostdev__ inline CoordT RoundDown(const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline CoordT RoundDown(const Vec3T& xyz) noexcept { return CoordT(Floor(xyz[0]), Floor(xyz[1]), Floor(xyz[2])); } //@{ /// Return the square root of a floating-point value. -__hostdev__ inline float Sqrt(float x) +__hostdev__ [[nodiscard]] inline float Sqrt(float x) noexcept { return sqrtf(x); } -__hostdev__ inline double Sqrt(double x) +__hostdev__ [[nodiscard]] inline double Sqrt(double x) noexcept { return sqrt(x); } @@ -286,13 +349,15 @@ __hostdev__ inline double Sqrt(double x) /// Return the sign of the given value as an integer (either -1, 0 or 1). template -__hostdev__ inline T Sign(const T& x) +__hostdev__ [[nodiscard]] inline constexpr T Sign(const T& x) noexcept { return ((T(0) < x) ? T(1) : T(0)) - ((x < T(0)) ? T(1) : T(0)); } +/// @brief Return the component index (0, 1, or 2) of the smallest element of @c v. +/// @note Ties resolve to the lowest index (i.e. x beats y, y beats z). template -__hostdev__ inline int MinIndex(const Vec3T& v) +__hostdev__ [[nodiscard]] inline constexpr int MinIndex(const Vec3T& v) noexcept { #if 0 static const int hashTable[8] = {2, 1, 9, 1, 2, 9, 0, 0}; //9 are dummy values @@ -308,8 +373,10 @@ __hostdev__ inline int MinIndex(const Vec3T& v) #endif } +/// @brief Return the component index (0, 1, or 2) of the largest element of @c v. +/// @note Ties resolve to the lowest index (i.e. x beats y, y beats z). template -__hostdev__ inline int MaxIndex(const Vec3T& v) +__hostdev__ [[nodiscard]] inline constexpr int MaxIndex(const Vec3T& v) noexcept { #if 0 static const int hashTable[8] = {2, 1, 9, 1, 2, 9, 0, 0}; //9 are dummy values @@ -329,7 +396,7 @@ __hostdev__ inline int MaxIndex(const Vec3T& v) /// /// @details both wordSize and byteSize are in byte units template -__hostdev__ inline uint64_t AlignUp(uint64_t byteCount) +__hostdev__ [[nodiscard]] inline constexpr uint64_t AlignUp(uint64_t byteCount) noexcept { const uint64_t r = byteCount % wordSize; return r ? byteCount - r + wordSize : byteCount; @@ -350,53 +417,60 @@ class Coord using IndexType = uint32_t; /// @brief Initialize all coordinates to zero. - __hostdev__ Coord() + __hostdev__ constexpr Coord() noexcept : mVec{0, 0, 0} { } /// @brief Initializes all coordinates to the given signed integer. - __hostdev__ explicit Coord(ValueType n) + __hostdev__ explicit constexpr Coord(ValueType n) noexcept : mVec{n, n, n} { } /// @brief Initializes coordinate to the given signed integers. - __hostdev__ Coord(ValueType i, ValueType j, ValueType k) + __hostdev__ constexpr Coord(ValueType i, ValueType j, ValueType k) noexcept : mVec{i, j, k} { } - __hostdev__ Coord(ValueType* ptr) + /// @brief Read three signed integers from @a ptr (no bounds check). + __hostdev__ constexpr Coord(ValueType* ptr) noexcept : mVec{ptr[0], ptr[1], ptr[2]} { } - __hostdev__ int32_t x() const { return mVec[0]; } - __hostdev__ int32_t y() const { return mVec[1]; } - __hostdev__ int32_t z() const { return mVec[2]; } + //@{ + /// @brief Named component accessors (x = mVec[0], y = mVec[1], z = mVec[2]). + __hostdev__ constexpr int32_t x() const noexcept { return mVec[0]; } + __hostdev__ constexpr int32_t y() const noexcept { return mVec[1]; } + __hostdev__ constexpr int32_t z() const noexcept { return mVec[2]; } - __hostdev__ int32_t& x() { return mVec[0]; } - __hostdev__ int32_t& y() { return mVec[1]; } - __hostdev__ int32_t& z() { return mVec[2]; } + __hostdev__ constexpr int32_t& x() noexcept { return mVec[0]; } + __hostdev__ constexpr int32_t& y() noexcept { return mVec[1]; } + __hostdev__ constexpr int32_t& z() noexcept { return mVec[2]; } + //@} - __hostdev__ static Coord max() { return Coord(int32_t((1u << 31) - 1)); } + /// @brief Largest representable @c Coord (all components = INT32_MAX). + __hostdev__ [[nodiscard]] static constexpr Coord max() noexcept { return Coord(int32_t((1u << 31) - 1)); } - __hostdev__ static Coord min() { return Coord(-int32_t((1u << 31) - 1) - 1); } + /// @brief Smallest representable @c Coord (all components = INT32_MIN). + __hostdev__ [[nodiscard]] static constexpr Coord min() noexcept { return Coord(-int32_t((1u << 31) - 1) - 1); } - __hostdev__ static size_t memUsage() { return sizeof(Coord); } + /// @brief Byte size of a @c Coord (always 12 on a 32-bit @c int32_t platform). + __hostdev__ [[nodiscard]] static constexpr size_t memUsage() noexcept { return sizeof(Coord); } /// @brief Return a const reference to the given Coord component. /// @warning The argument is assumed to be 0, 1, or 2. - __hostdev__ const ValueType& operator[](IndexType i) const { return mVec[i]; } + __hostdev__ constexpr const ValueType& operator[](IndexType i) const noexcept { NANOVDB_ASSERT(i < 3); return mVec[i]; } /// @brief Return a non-const reference to the given Coord component. /// @warning The argument is assumed to be 0, 1, or 2. - __hostdev__ ValueType& operator[](IndexType i) { return mVec[i]; } + __hostdev__ constexpr ValueType& operator[](IndexType i) noexcept { NANOVDB_ASSERT(i < 3); return mVec[i]; } /// @brief Assignment operator that works with openvdb::Coord template - __hostdev__ Coord& operator=(const CoordT& other) + __hostdev__ constexpr Coord& operator=(const CoordT& other) noexcept { static_assert(sizeof(Coord) == sizeof(CoordT), "Mis-matched sizeof"); mVec[0] = other[0]; @@ -406,16 +480,16 @@ class Coord } /// @brief Return a new instance with coordinates masked by the given unsigned integer. - __hostdev__ Coord operator&(IndexType n) const { return Coord(mVec[0] & n, mVec[1] & n, mVec[2] & n); } + __hostdev__ [[nodiscard]] constexpr Coord operator&(IndexType n) const noexcept { return Coord(mVec[0] & n, mVec[1] & n, mVec[2] & n); } - // @brief Return a new instance with coordinates left-shifted by the given unsigned integer. - __hostdev__ Coord operator<<(IndexType n) const { return Coord(mVec[0] << n, mVec[1] << n, mVec[2] << n); } + /// @brief Return a new instance with coordinates left-shifted by the given unsigned integer. + __hostdev__ [[nodiscard]] constexpr Coord operator<<(IndexType n) const noexcept { return Coord(mVec[0] << n, mVec[1] << n, mVec[2] << n); } - // @brief Return a new instance with coordinates right-shifted by the given unsigned integer. - __hostdev__ Coord operator>>(IndexType n) const { return Coord(mVec[0] >> n, mVec[1] >> n, mVec[2] >> n); } + /// @brief Return a new instance with coordinates right-shifted by the given unsigned integer. + __hostdev__ [[nodiscard]] constexpr Coord operator>>(IndexType n) const noexcept { return Coord(mVec[0] >> n, mVec[1] >> n, mVec[2] >> n); } /// @brief Return true if this Coord is lexicographically less than the given Coord. - __hostdev__ bool operator<(const Coord& rhs) const + __hostdev__ [[nodiscard]] constexpr bool operator<(const Coord& rhs) const noexcept { return mVec[0] < rhs[0] ? true : mVec[0] > rhs[0] ? false @@ -425,7 +499,7 @@ class Coord } /// @brief Return true if this Coord is lexicographically less or equal to the given Coord. - __hostdev__ bool operator<=(const Coord& rhs) const + __hostdev__ [[nodiscard]] constexpr bool operator<=(const Coord& rhs) const noexcept { return mVec[0] < rhs[0] ? true : mVec[0] > rhs[0] ? false @@ -434,8 +508,8 @@ class Coord : mVec[2] <=rhs[2] ? true : false; } - // @brief Return true if this Coord is lexicographically greater than the given Coord. - __hostdev__ bool operator>(const Coord& rhs) const + /// @brief Return true if this Coord is lexicographically greater than the given Coord. + __hostdev__ [[nodiscard]] constexpr bool operator>(const Coord& rhs) const noexcept { return mVec[0] > rhs[0] ? true : mVec[0] < rhs[0] ? false @@ -444,8 +518,8 @@ class Coord : mVec[2] > rhs[2] ? true : false; } - // @brief Return true if this Coord is lexicographically greater or equal to the given Coord. - __hostdev__ bool operator>=(const Coord& rhs) const + /// @brief Return true if this Coord is lexicographically greater or equal to the given Coord. + __hostdev__ [[nodiscard]] constexpr bool operator>=(const Coord& rhs) const noexcept { return mVec[0] > rhs[0] ? true : mVec[0] < rhs[0] ? false @@ -454,48 +528,58 @@ class Coord : mVec[2] >=rhs[2] ? true : false; } - // @brief Return true if the Coord components are identical. - __hostdev__ bool operator==(const Coord& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2]; } - __hostdev__ bool operator!=(const Coord& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2]; } - __hostdev__ Coord& operator&=(int n) + /// @brief Return true iff every component matches @a rhs. + __hostdev__ [[nodiscard]] constexpr bool operator==(const Coord& rhs) const noexcept { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2]; } + /// @brief Return true iff any component differs from @a rhs. + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Coord& rhs) const noexcept { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2]; } + /// @brief In-place component-wise bitwise AND with the given mask. + __hostdev__ constexpr Coord& operator&=(int n) noexcept { mVec[0] &= n; mVec[1] &= n; mVec[2] &= n; return *this; } - __hostdev__ Coord& operator<<=(uint32_t n) + /// @brief In-place left-shift of every component by @a n bits. + __hostdev__ constexpr Coord& operator<<=(uint32_t n) noexcept { mVec[0] <<= n; mVec[1] <<= n; mVec[2] <<= n; return *this; } - __hostdev__ Coord& operator>>=(uint32_t n) + /// @brief In-place right-shift of every component by @a n bits. + __hostdev__ constexpr Coord& operator>>=(uint32_t n) noexcept { mVec[0] >>= n; mVec[1] >>= n; mVec[2] >>= n; return *this; } - __hostdev__ Coord& operator+=(int n) + /// @brief Add a scalar to every component in place. + __hostdev__ constexpr Coord& operator+=(int n) noexcept { mVec[0] += n; mVec[1] += n; mVec[2] += n; return *this; } - __hostdev__ Coord operator+(const Coord& rhs) const { return Coord(mVec[0] + rhs[0], mVec[1] + rhs[1], mVec[2] + rhs[2]); } - __hostdev__ Coord operator-(const Coord& rhs) const { return Coord(mVec[0] - rhs[0], mVec[1] - rhs[1], mVec[2] - rhs[2]); } - __hostdev__ Coord operator-() const { return Coord(-mVec[0], -mVec[1], -mVec[2]); } - __hostdev__ Coord& operator+=(const Coord& rhs) + /// @brief Component-wise sum of two coordinates. + __hostdev__ [[nodiscard]] constexpr Coord operator+(const Coord& rhs) const noexcept { return Coord(mVec[0] + rhs[0], mVec[1] + rhs[1], mVec[2] + rhs[2]); } + /// @brief Component-wise difference of two coordinates. + __hostdev__ [[nodiscard]] constexpr Coord operator-(const Coord& rhs) const noexcept { return Coord(mVec[0] - rhs[0], mVec[1] - rhs[1], mVec[2] - rhs[2]); } + /// @brief Component-wise negation. + __hostdev__ [[nodiscard]] constexpr Coord operator-() const noexcept { return Coord(-mVec[0], -mVec[1], -mVec[2]); } + /// @brief In-place component-wise addition with another coordinate. + __hostdev__ constexpr Coord& operator+=(const Coord& rhs) noexcept { mVec[0] += rhs[0]; mVec[1] += rhs[1]; mVec[2] += rhs[2]; return *this; } - __hostdev__ Coord& operator-=(const Coord& rhs) + /// @brief In-place component-wise subtraction with another coordinate. + __hostdev__ constexpr Coord& operator-=(const Coord& rhs) noexcept { mVec[0] -= rhs[0]; mVec[1] -= rhs[1]; @@ -504,7 +588,7 @@ class Coord } /// @brief Perform a component-wise minimum with the other Coord. - __hostdev__ Coord& minComponent(const Coord& other) + __hostdev__ constexpr Coord& minComponent(const Coord& other) noexcept { if (other[0] < mVec[0]) mVec[0] = other[0]; @@ -516,7 +600,7 @@ class Coord } /// @brief Perform a component-wise maximum with the other Coord. - __hostdev__ Coord& maxComponent(const Coord& other) + __hostdev__ constexpr Coord& maxComponent(const Coord& other) noexcept { if (other[0] > mVec[0]) mVec[0] = other[0]; @@ -527,14 +611,16 @@ class Coord return *this; } #if defined(__CUDACC__) // the following functions only run on the GPU! - __device__ inline Coord& minComponentAtomic(const Coord& other) + /// @brief Device-only @c atomicMin component-wise. + __device__ inline Coord& minComponentAtomic(const Coord& other) noexcept { atomicMin(&mVec[0], other[0]); atomicMin(&mVec[1], other[1]); atomicMin(&mVec[2], other[2]); return *this; } - __device__ inline Coord& maxComponentAtomic(const Coord& other) + /// @brief Device-only @c atomicMax component-wise. + __device__ inline Coord& maxComponentAtomic(const Coord& other) noexcept { atomicMax(&mVec[0], other[0]); atomicMax(&mVec[1], other[1]); @@ -543,16 +629,18 @@ class Coord } #endif - __hostdev__ Coord offsetBy(ValueType dx, ValueType dy, ValueType dz) const + /// @brief Return a new @c Coord offset component-by-component by @a (dx, dy, dz). + __hostdev__ [[nodiscard]] constexpr Coord offsetBy(ValueType dx, ValueType dy, ValueType dz) const noexcept { return Coord(mVec[0] + dx, mVec[1] + dy, mVec[2] + dz); } - __hostdev__ Coord offsetBy(ValueType n) const { return this->offsetBy(n, n, n); } + /// @brief Return a new @c Coord offset by the same scalar @a n in every component. + __hostdev__ [[nodiscard]] constexpr Coord offsetBy(ValueType n) const noexcept { return this->offsetBy(n, n, n); } /// Return true if any of the components of @a a are smaller than the /// corresponding components of @a b. - __hostdev__ static inline bool lessThan(const Coord& a, const Coord& b) + __hostdev__ [[nodiscard]] static inline constexpr bool lessThan(const Coord& a, const Coord& b) noexcept { return (a[0] < b[0] || a[1] < b[1] || a[2] < b[2]); } @@ -560,7 +648,7 @@ class Coord /// @brief Return the largest integer coordinates that are not greater /// than @a xyz (node centered conversion). template - __hostdev__ static Coord Floor(const Vec3T& xyz) { return Coord(math::Floor(xyz[0]), math::Floor(xyz[1]), math::Floor(xyz[2])); } + __hostdev__ [[nodiscard]] static Coord Floor(const Vec3T& xyz) noexcept { return Coord(math::Floor(xyz[0]), math::Floor(xyz[1]), math::Floor(xyz[2])); } /// @brief Return a hash key derived from the existing coordinates. /// @details The hash function is originally taken from the SIGGRAPH paper: @@ -568,22 +656,22 @@ class Coord /// and the prime numbers are modified based on the ACM Transactions on Graphics paper: /// "Real-time 3D reconstruction at scale using voxel hashing" (the second number had a typo!) template - __hostdev__ uint32_t hash() const { return ((1 << Log2N) - 1) & (mVec[0] * 73856093 ^ mVec[1] * 19349669 ^ mVec[2] * 83492791); } + __hostdev__ [[nodiscard]] constexpr uint32_t hash() const noexcept { return ((1 << Log2N) - 1) & (mVec[0] * 73856093 ^ mVec[1] * 19349669 ^ mVec[2] * 83492791); } /// @brief Return the octant of this Coord - //__hostdev__ size_t octant() const { return (uint32_t(mVec[0])>>31) | ((uint32_t(mVec[1])>>31)<<1) | ((uint32_t(mVec[2])>>31)<<2); } - __hostdev__ uint8_t octant() const { return (uint8_t(bool(mVec[0] & (1u << 31)))) | + __hostdev__ [[nodiscard]] constexpr uint8_t octant() const noexcept { return (uint8_t(bool(mVec[0] & (1u << 31)))) | (uint8_t(bool(mVec[1] & (1u << 31))) << 1) | (uint8_t(bool(mVec[2] & (1u << 31))) << 2); } /// @brief Return a single precision floating-point vector of this coordinate - __hostdev__ inline Vec3 asVec3s() const; + __hostdev__ [[nodiscard]] inline constexpr Vec3 asVec3s() const noexcept; /// @brief Return a double precision floating-point vector of this coordinate - __hostdev__ inline Vec3 asVec3d() const; + __hostdev__ [[nodiscard]] inline constexpr Vec3 asVec3d() const noexcept; - // returns a copy of itself, so it mimics the behaviour of Vec3::round() - __hostdev__ inline Coord round() const { return *this; } + /// @brief Identity (Coord is already integer); provided so generic code + /// can call @c .round() uniformly on @c Coord and @c Vec3. + __hostdev__ [[nodiscard]] inline constexpr Coord round() const noexcept { return *this; } }; // Coord class @@ -603,51 +691,58 @@ class Coord2 using IndexType = uint32_t; /// @brief Initialize all coordinates to zero. - __hostdev__ Coord2() + __hostdev__ constexpr Coord2() noexcept : mVec{0, 0} { } /// @brief Initializes all coordinates to the given signed integer. - __hostdev__ explicit Coord2(ValueType n) + __hostdev__ explicit constexpr Coord2(ValueType n) noexcept : mVec{n, n} { } /// @brief Initializes coordinate to the given signed integers. - __hostdev__ Coord2(ValueType i, ValueType j) + __hostdev__ constexpr Coord2(ValueType i, ValueType j) noexcept : mVec{i, j} { } - __hostdev__ Coord2(ValueType* ptr) + /// @brief Read two signed integers from @a ptr (no bounds check). + __hostdev__ constexpr Coord2(ValueType* ptr) noexcept : mVec{ptr[0], ptr[1]} { } - __hostdev__ int32_t x() const { return mVec[0]; } - __hostdev__ int32_t y() const { return mVec[1]; } + //@{ + /// @brief Named component accessors (x = mVec[0], y = mVec[1]). + __hostdev__ constexpr int32_t x() const noexcept { return mVec[0]; } + __hostdev__ constexpr int32_t y() const noexcept { return mVec[1]; } - __hostdev__ int32_t& x() { return mVec[0]; } - __hostdev__ int32_t& y() { return mVec[1]; } + __hostdev__ constexpr int32_t& x() noexcept { return mVec[0]; } + __hostdev__ constexpr int32_t& y() noexcept { return mVec[1]; } + //@} - __hostdev__ static Coord2 max() { return Coord2(int32_t((1u << 31) - 1)); } + /// @brief Largest representable @c Coord2 (all components = INT32_MAX). + __hostdev__ [[nodiscard]] static constexpr Coord2 max() noexcept { return Coord2(int32_t((1u << 31) - 1)); } - __hostdev__ static Coord2 min() { return Coord2(-int32_t((1u << 31) - 1) - 1); } + /// @brief Smallest representable @c Coord2 (all components = INT32_MIN). + __hostdev__ [[nodiscard]] static constexpr Coord2 min() noexcept { return Coord2(-int32_t((1u << 31) - 1) - 1); } - __hostdev__ static size_t memUsage() { return sizeof(Coord2); } + /// @brief Byte size of a @c Coord2 (always 8 on a 32-bit @c int32_t platform). + __hostdev__ [[nodiscard]] static constexpr size_t memUsage() noexcept { return sizeof(Coord2); } /// @brief Return a const reference to the given Coord component. /// @warning The argument is assumed to be 0 or 1, - __hostdev__ const ValueType& operator[](IndexType i) const { return mVec[i]; } + __hostdev__ constexpr const ValueType& operator[](IndexType i) const noexcept { NANOVDB_ASSERT(i < 2); return mVec[i]; } /// @brief Return a non-const reference to the given Coord component. /// @warning The argument is assumed to be 0 or 1. - __hostdev__ ValueType& operator[](IndexType i) { return mVec[i]; } + __hostdev__ constexpr ValueType& operator[](IndexType i) noexcept { NANOVDB_ASSERT(i < 2); return mVec[i]; } /// @brief Assignment operator that works with openvdb::Coord template - __hostdev__ Coord2& operator=(const CoordT& other) + __hostdev__ constexpr Coord2& operator=(const CoordT& other) noexcept { static_assert(sizeof(Coord2) == sizeof(CoordT), "Mis-matched sizeof"); mVec[0] = other[0]; @@ -656,16 +751,16 @@ class Coord2 } /// @brief Return a new instance with coordinates masked by the given unsigned integer. - __hostdev__ Coord2 operator&(IndexType n) const { return Coord2(mVec[0] & n, mVec[1] & n); } + __hostdev__ [[nodiscard]] constexpr Coord2 operator&(IndexType n) const noexcept { return Coord2(mVec[0] & n, mVec[1] & n); } - // @brief Return a new instance with coordinates left-shifted by the given unsigned integer. - __hostdev__ Coord2 operator<<(IndexType n) const { return Coord2(mVec[0] << n, mVec[1] << n); } + /// @brief Return a new instance with coordinates left-shifted by the given unsigned integer. + __hostdev__ [[nodiscard]] constexpr Coord2 operator<<(IndexType n) const noexcept { return Coord2(mVec[0] << n, mVec[1] << n); } - // @brief Return a new instance with coordinates right-shifted by the given unsigned integer. - __hostdev__ Coord2 operator>>(IndexType n) const { return Coord2(mVec[0] >> n, mVec[1] >> n); } + /// @brief Return a new instance with coordinates right-shifted by the given unsigned integer. + __hostdev__ [[nodiscard]] constexpr Coord2 operator>>(IndexType n) const noexcept { return Coord2(mVec[0] >> n, mVec[1] >> n); } /// @brief Return true if this Coord is lexicographically less than the given Coord. - __hostdev__ bool operator<(const Coord2& rhs) const + __hostdev__ [[nodiscard]] constexpr bool operator<(const Coord2& rhs) const noexcept { return mVec[0] < rhs[0] ? true : mVec[0] > rhs[0] ? false @@ -673,66 +768,76 @@ class Coord2 } /// @brief Return true if this Coord is lexicographically less or equal to the given Coord. - __hostdev__ bool operator<=(const Coord2& rhs) const + __hostdev__ [[nodiscard]] constexpr bool operator<=(const Coord2& rhs) const noexcept { return mVec[0] < rhs[0] ? true : mVec[0] > rhs[0] ? false : mVec[1] <= rhs[1] ? true : false; } - // @brief Return true if this Coord is lexicographically greater than the given Coord. - __hostdev__ bool operator>(const Coord2& rhs) const + /// @brief Return true if this Coord is lexicographically greater than the given Coord. + __hostdev__ [[nodiscard]] constexpr bool operator>(const Coord2& rhs) const noexcept { return mVec[0] > rhs[0] ? true : mVec[0] < rhs[0] ? false : mVec[1] > rhs[1] ? true : false; } - // @brief Return true if this Coord is lexicographically greater or equal to the given Coord. - __hostdev__ bool operator>=(const Coord2& rhs) const + /// @brief Return true if this Coord is lexicographically greater or equal to the given Coord. + __hostdev__ [[nodiscard]] constexpr bool operator>=(const Coord2& rhs) const noexcept { return mVec[0] > rhs[0] ? true : mVec[0] < rhs[0] ? false : mVec[1] >= rhs[1] ? true : false; } - // @brief Return true if the Coord components are identical. - __hostdev__ bool operator==(const Coord2& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1]; } - __hostdev__ bool operator!=(const Coord2& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1]; } - __hostdev__ Coord2& operator&=(int n) + /// @brief Return true iff every component matches @a rhs. + __hostdev__ [[nodiscard]] constexpr bool operator==(const Coord2& rhs) const noexcept { return mVec[0] == rhs[0] && mVec[1] == rhs[1]; } + /// @brief Return true iff any component differs from @a rhs. + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Coord2& rhs) const noexcept { return mVec[0] != rhs[0] || mVec[1] != rhs[1]; } + /// @brief In-place component-wise bitwise AND with the given mask. + __hostdev__ constexpr Coord2& operator&=(int n) noexcept { mVec[0] &= n; mVec[1] &= n; return *this; } - __hostdev__ Coord2& operator<<=(uint32_t n) + /// @brief In-place left-shift of every component by @a n bits. + __hostdev__ constexpr Coord2& operator<<=(uint32_t n) noexcept { mVec[0] <<= n; mVec[1] <<= n; return *this; } - __hostdev__ Coord2& operator>>=(uint32_t n) + /// @brief In-place right-shift of every component by @a n bits. + __hostdev__ constexpr Coord2& operator>>=(uint32_t n) noexcept { mVec[0] >>= n; mVec[1] >>= n; return *this; } - __hostdev__ Coord2& operator+=(int n) + /// @brief Add a scalar to every component in place. + __hostdev__ constexpr Coord2& operator+=(int n) noexcept { mVec[0] += n; mVec[1] += n; return *this; } - __hostdev__ Coord2 operator+(const Coord2& rhs) const { return Coord2(mVec[0] + rhs[0], mVec[1] + rhs[1]); } - __hostdev__ Coord2 operator-(const Coord2& rhs) const { return Coord2(mVec[0] - rhs[0], mVec[1] - rhs[1]); } - __hostdev__ Coord2 operator-() const { return Coord2(-mVec[0], -mVec[1]); } - __hostdev__ Coord2& operator+=(const Coord2& rhs) + /// @brief Component-wise sum of two coordinates. + __hostdev__ [[nodiscard]] constexpr Coord2 operator+(const Coord2& rhs) const noexcept { return Coord2(mVec[0] + rhs[0], mVec[1] + rhs[1]); } + /// @brief Component-wise difference of two coordinates. + __hostdev__ [[nodiscard]] constexpr Coord2 operator-(const Coord2& rhs) const noexcept { return Coord2(mVec[0] - rhs[0], mVec[1] - rhs[1]); } + /// @brief Component-wise negation. + __hostdev__ [[nodiscard]] constexpr Coord2 operator-() const noexcept { return Coord2(-mVec[0], -mVec[1]); } + /// @brief In-place component-wise addition with another coordinate. + __hostdev__ constexpr Coord2& operator+=(const Coord2& rhs) noexcept { mVec[0] += rhs[0]; mVec[1] += rhs[1]; return *this; } - __hostdev__ Coord2& operator-=(const Coord2& rhs) + /// @brief In-place component-wise subtraction with another coordinate. + __hostdev__ constexpr Coord2& operator-=(const Coord2& rhs) noexcept { mVec[0] -= rhs[0]; mVec[1] -= rhs[1]; @@ -740,7 +845,7 @@ class Coord2 } /// @brief Perform a component-wise minimum with the other Coord. - __hostdev__ Coord2& minComponent(const Coord2& other) + __hostdev__ constexpr Coord2& minComponent(const Coord2& other) noexcept { if (other[0] < mVec[0]) mVec[0] = other[0]; @@ -750,7 +855,7 @@ class Coord2 } /// @brief Perform a component-wise maximum with the other Coord. - __hostdev__ Coord2& maxComponent(const Coord2& other) + __hostdev__ constexpr Coord2& maxComponent(const Coord2& other) noexcept { if (other[0] > mVec[0]) mVec[0] = other[0]; @@ -759,13 +864,15 @@ class Coord2 return *this; } #if defined(__CUDACC__) // the following functions only run on the GPU! - __device__ inline Coord2& minComponentAtomic(const Coord2& other) + /// @brief Device-only @c atomicMin component-wise. + __device__ inline Coord2& minComponentAtomic(const Coord2& other) noexcept { atomicMin(&mVec[0], other[0]); atomicMin(&mVec[1], other[1]); return *this; } - __device__ inline Coord2& maxComponentAtomic(const Coord2& other) + /// @brief Device-only @c atomicMax component-wise. + __device__ inline Coord2& maxComponentAtomic(const Coord2& other) noexcept { atomicMax(&mVec[0], other[0]); atomicMax(&mVec[1], other[1]); @@ -773,16 +880,18 @@ class Coord2 } #endif - __hostdev__ Coord2 offsetBy(ValueType dx, ValueType dy) const + /// @brief Return a new @c Coord2 offset component-by-component by @a (dx, dy). + __hostdev__ [[nodiscard]] constexpr Coord2 offsetBy(ValueType dx, ValueType dy) const noexcept { return Coord2(mVec[0] + dx, mVec[1] + dy); } - __hostdev__ Coord2 offsetBy(ValueType n) const { return this->offsetBy(n, n); } + /// @brief Return a new @c Coord2 offset by the same scalar @a n in both components. + __hostdev__ [[nodiscard]] constexpr Coord2 offsetBy(ValueType n) const noexcept { return this->offsetBy(n, n); } /// Return true if any of the components of @a a are smaller than the /// corresponding components of @a b. - __hostdev__ static inline bool lessThan(const Coord2& a, const Coord2& b) + __hostdev__ [[nodiscard]] static inline constexpr bool lessThan(const Coord2& a, const Coord2& b) noexcept { return (a[0] < b[0] || a[1] < b[1]); } @@ -790,219 +899,602 @@ class Coord2 /// @brief Return the largest integer coordinates that are not greater /// than @a xyz (node centered conversion). template - __hostdev__ static Coord2 Floor(const Vec2T& xy) { return Coord2(math::Floor(xy[0]), math::Floor(xy[1])); } + __hostdev__ [[nodiscard]] static Coord2 Floor(const Vec2T& xy) noexcept { return Coord2(math::Floor(xy[0]), math::Floor(xy[1])); } /// @brief Return a single precision floating-point vector of this coordinate - __hostdev__ inline Vec2 asVec2s() const; + __hostdev__ [[nodiscard]] inline constexpr Vec2 asVec2s() const noexcept; /// @brief Return a double precision floating-point vector of this coordinate - __hostdev__ inline Vec2 asVec2d() const; + __hostdev__ [[nodiscard]] inline constexpr Vec2 asVec2d() const noexcept; - // returns a copy of itself, so it mimics the behaviour of Vec3::round() - __hostdev__ inline Coord2 round() const { return *this; } + /// @brief Identity (Coord2 is already integer); provided so generic code + /// can call @c .round() uniformly on @c Coord2 and @c Vec2. + __hostdev__ [[nodiscard]] inline constexpr Coord2 round() const noexcept { return *this; } }; // Coord2 class -// ----------------------------> Vec2 <-------------------------------------- +// ----------------------------> VecBase <----------------------------------- -/// @brief A simple vector class with three components, similar to openvdb::math::Vec3 -template -class Vec2 +/// @brief Base class for fixed-size vectors. Provides shared element-wise +/// arithmetic, scalar arithmetic, equality, length/dot reductions, +/// component min/max, and integer-coord rounding. Derived classes provide +/// constructors and any dimension-specific extras (e.g. cross/outer on Vec3). +template +class VecBase { - T mVec[2]; +protected: + T mVec[N]; public: - static const int SIZE = 2; - static const int size = 2; // in openvdb::math::Tuple using ValueType = T; - Vec2() = default; - __hostdev__ explicit Vec2(T x) - : mVec{x, x} - { - } - __hostdev__ Vec2(T x, T y) - : mVec{x, y} - { - } - template class Vec2T, class T2> - __hostdev__ Vec2(const Vec2T& v) - : mVec{T(v[0]), T(v[1])} + static constexpr int SIZE = N; + + /// @brief Default-construct (mVec is left uninitialized for trivially-default-constructible @c T). + VecBase() noexcept = default; + +protected: + /// @brief Variadic component-init ctor — used only by derived ctors + /// (`Vec2(T x, T y) : Base(x, y) {}`). The member-init list directly + /// initializes @c mVec, avoiding the default-construct-then-assign that + /// a body-assignment ctor would impose on a non-fundamental @c T. + template + __hostdev__ explicit constexpr VecBase(Args... args) noexcept : mVec{T(args)...} { - static_assert(Vec2T::size == size, "expected Vec2T::size==2!"); + static_assert(sizeof...(Args) == N, "VecBase: wrong number of constructor arguments"); } - template - __hostdev__ explicit Vec2(const Vec2& v) - : mVec{T(v[0]), T(v[1])} - { + +public: + /// @brief Indexed element access. Asserts 0 <= i < N in debug builds. + __hostdev__ constexpr const T& operator[](int i) const noexcept { NANOVDB_ASSERT(i >= 0 && i < N); return mVec[i]; } + __hostdev__ constexpr T& operator[](int i) noexcept { NANOVDB_ASSERT(i >= 0 && i < N); return mVec[i]; } + + __hostdev__ constexpr const T& x() const noexcept { return mVec[0]; } + __hostdev__ constexpr T& x() noexcept { return mVec[0]; } + __hostdev__ constexpr const T& y() const noexcept { static_assert(N >= 2, "VecBase::y() requires N >= 2"); return mVec[1]; } + __hostdev__ constexpr T& y() noexcept { static_assert(N >= 2, "VecBase::y() requires N >= 2"); return mVec[1]; } + __hostdev__ constexpr const T& z() const noexcept { static_assert(N >= 3, "VecBase::z() requires N >= 3"); return mVec[2]; } + __hostdev__ constexpr T& z() noexcept { static_assert(N >= 3, "VecBase::z() requires N >= 3"); return mVec[2]; } + __hostdev__ constexpr const T& w() const noexcept { static_assert(N >= 4, "VecBase::w() requires N >= 4"); return mVec[3]; } + __hostdev__ constexpr T& w() noexcept { static_assert(N >= 4, "VecBase::w() requires N >= 4"); return mVec[3]; } + + /// @brief raw pointer to the underlying N-element storage + __hostdev__ constexpr T* asPointer() noexcept { return mVec; } + __hostdev__ constexpr const T* asPointer() const noexcept { return mVec; } + + // ---- generic element-wise helpers (taking/returning Derived) ---- + + /// @brief Return @c *this + @a rhs as a @c Derived. + template + __hostdev__ [[nodiscard]] constexpr Derived plus(const Derived& rhs) const noexcept { + Derived out{}; + for (int i = 0; i < N; ++i) out[i] = mVec[i] + rhs[i]; + return out; + } + /// @brief Return @c *this - @a rhs as a @c Derived. + template + __hostdev__ [[nodiscard]] constexpr Derived minus(const Derived& rhs) const noexcept { + Derived out{}; + for (int i = 0; i < N; ++i) out[i] = mVec[i] - rhs[i]; + return out; + } + /// @brief Return component-wise @c *this * @a rhs as a @c Derived (Hadamard product). + template + __hostdev__ [[nodiscard]] constexpr Derived mul(const Derived& rhs) const noexcept { + Derived out{}; + for (int i = 0; i < N; ++i) out[i] = mVec[i] * rhs[i]; + return out; + } + /// @brief Return component-wise @c *this / @a rhs as a @c Derived. + template + __hostdev__ [[nodiscard]] constexpr Derived div(const Derived& rhs) const noexcept { + Derived out{}; + for (int i = 0; i < N; ++i) out[i] = mVec[i] / rhs[i]; + return out; + } + /// @brief Return @c -(*this) as a @c Derived. + template + __hostdev__ [[nodiscard]] constexpr Derived negate() const noexcept { + Derived out{}; + for (int i = 0; i < N; ++i) out[i] = -mVec[i]; + return out; + } + /// @brief Return @c s * (*this) as a @c Derived (scalar broadcast). + template + __hostdev__ [[nodiscard]] constexpr Derived scale(const T& s) const noexcept { + Derived out{}; + for (int i = 0; i < N; ++i) out[i] = mVec[i] * s; + return out; + } + /// @brief Return @c (*this) / @a s element-wise as a @c Derived. Uses per-element + /// division (correct for integer @c T, unlike multiplying by 1/s). + template + __hostdev__ [[nodiscard]] constexpr Derived divideBy(const T& s) const noexcept { + Derived out{}; + for (int i = 0; i < N; ++i) out[i] = mVec[i] / s; + return out; + } + /// @brief Return a unit-length copy (@c *this divided by @c length()) as a @c Derived. + /// Const, non-mutating counterpart of the derived @c normalize(). Not @c constexpr — + /// @c length() calls @c Sqrt. + template + __hostdev__ [[nodiscard]] Derived normalized() const noexcept { + return this->template divideBy(this->length()); + } + + // ---- in-place compound assignment helpers ---- + + /// @brief Component-wise add @a rhs into @c *this and return @c *this. + __hostdev__ constexpr VecBase& addAssign(const VecBase& rhs) noexcept { + for (int i = 0; i < N; ++i) mVec[i] += rhs.mVec[i]; + return *this; } - __hostdev__ explicit Vec2(const Coord2& ijk) - : mVec{T(ijk[0]), T(ijk[1])} - { + /// @brief Component-wise subtract @a rhs from @c *this and return @c *this. + __hostdev__ constexpr VecBase& subAssign(const VecBase& rhs) noexcept { + for (int i = 0; i < N; ++i) mVec[i] -= rhs.mVec[i]; + return *this; } - __hostdev__ bool operator==(const Vec2& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1]; } - __hostdev__ bool operator!=(const Vec2& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1]; } - template class Vec2T, class T2> - __hostdev__ Vec2& operator=(const Vec2T& rhs) - { - static_assert(Vec2T::size == size, "expected Vec2T::size==2!"); - mVec[0] = rhs[0]; - mVec[1] = rhs[1]; + /// @brief Component-wise multiply @c *this by @a rhs (Hadamard product) and return @c *this. + __hostdev__ constexpr VecBase& mulAssign(const VecBase& rhs) noexcept { + for (int i = 0; i < N; ++i) mVec[i] *= rhs.mVec[i]; return *this; } - __hostdev__ const T& operator[](int i) const { return mVec[i]; } - __hostdev__ T& operator[](int i) { return mVec[i]; } - template - __hostdev__ T dot(const Vec2T& v) const { return mVec[0] * v[0] + mVec[1] * v[1]; } - __hostdev__ T lengthSqr() const - { - return mVec[0] * mVec[0] + mVec[1] * mVec[1]; // 3 flops - } - __hostdev__ T length() const { return Sqrt(this->lengthSqr()); } - __hostdev__ Vec2 operator-() const { return Vec2(-mVec[0], -mVec[1]); } - __hostdev__ Vec2 operator*(const Vec2& v) const { return Vec2(mVec[0] * v[0], mVec[1] * v[1]); } - __hostdev__ Vec2 operator/(const Vec2& v) const { return Vec2(mVec[0] / v[0], mVec[1] / v[1]); } - __hostdev__ Vec2 operator+(const Vec2& v) const { return Vec2(mVec[0] + v[0], mVec[1] + v[1]); } - __hostdev__ Vec2 operator-(const Vec2& v) const { return Vec2(mVec[0] - v[0], mVec[1] - v[1]); } - __hostdev__ Vec2 operator+(const Coord& ijk) const { return Vec2(mVec[0] + ijk[0], mVec[1] + ijk[1]); } - __hostdev__ Vec2 operator-(const Coord& ijk) const { return Vec2(mVec[0] - ijk[0], mVec[1] - ijk[1]); } - __hostdev__ Vec2 operator*(const T& s) const { return Vec2(s * mVec[0], s * mVec[1]); } - __hostdev__ Vec2 operator/(const T& s) const { return (T(1) / s) * (*this); } - __hostdev__ Vec2& operator+=(const Vec2& v) - { - mVec[0] += v[0]; - mVec[1] += v[1]; + /// @brief Component-wise divide @c *this by @a rhs and return @c *this. + __hostdev__ constexpr VecBase& divAssign(const VecBase& rhs) noexcept { + for (int i = 0; i < N; ++i) mVec[i] /= rhs.mVec[i]; return *this; } - __hostdev__ Vec2& operator+=(const Coord& ijk) - { - mVec[0] += T(ijk[0]); - mVec[1] += T(ijk[1]); + /// @brief Multiply every component by scalar @a s in place; return @c *this. + __hostdev__ constexpr VecBase& scaleAssign(const T& s) noexcept { + for (int i = 0; i < N; ++i) mVec[i] *= s; return *this; } - __hostdev__ Vec2& operator-=(const Vec2& v) - { - mVec[0] -= v[0]; - mVec[1] -= v[1]; + /// @brief Divide every component by scalar @a s in place (per-element, integer-safe); + /// return @c *this. + __hostdev__ constexpr VecBase& divideAssignScalar(const T& s) noexcept { + for (int i = 0; i < N; ++i) mVec[i] /= s; return *this; } - __hostdev__ Vec2& operator-=(const Coord& ijk) - { - mVec[0] -= T(ijk[0]); - mVec[1] -= T(ijk[1]); + + /// @brief Return @c true iff every component compares equal to @a rhs. + __hostdev__ [[nodiscard]] constexpr bool equals(const VecBase& rhs) const noexcept { + for (int i = 0; i < N; ++i) if (mVec[i] != rhs.mVec[i]) return false; + return true; + } + + // ---- reductions ---- + + /// @brief dot product. @a V must have @c operator[] valid for 0..N-1. + template + __hostdev__ [[nodiscard]] constexpr T dot(const V& v) const noexcept { + T s = T(0); + for (int i = 0; i < N; ++i) s += mVec[i] * v[i]; + return s; + } + /// @brief Squared L2 length (sum of squared components). Constexpr, no sqrt. + __hostdev__ [[nodiscard]] constexpr T lengthSqr() const noexcept { + T s = T(0); + for (int i = 0; i < N; ++i) s += mVec[i] * mVec[i]; + return s; + } + /// @brief L2 length (Euclidean norm). Not @c constexpr — calls @c Sqrt. + __hostdev__ [[nodiscard]] T length() const noexcept { return Sqrt(this->lengthSqr()); } + + /// @brief return the smallest of the N components + __hostdev__ [[nodiscard]] constexpr T smallestComponent() const noexcept { + T m = mVec[0]; + for (int i = 1; i < N; ++i) if (mVec[i] < m) m = mVec[i]; + return m; + } + /// @brief return the largest of the N components + __hostdev__ [[nodiscard]] constexpr T largestComponent() const noexcept { + T m = mVec[0]; + for (int i = 1; i < N; ++i) if (mVec[i] > m) m = mVec[i]; + return m; + } + + // ---- component-wise (mutating) min/max ---- + + /// @brief Component-wise take the @c min of @c *this and @a other (in place); return @c *this. + template + __hostdev__ constexpr VecBase& mergeMin(const V& other) noexcept { + for (int i = 0; i < N; ++i) if (other[i] < mVec[i]) mVec[i] = other[i]; return *this; } - __hostdev__ Vec2& operator*=(const T& s) - { - mVec[0] *= s; - mVec[1] *= s; + /// @brief Component-wise take the @c max of @c *this and @a other (in place); return @c *this. + template + __hostdev__ constexpr VecBase& mergeMax(const V& other) noexcept { + for (int i = 0; i < N; ++i) if (other[i] > mVec[i]) mVec[i] = other[i]; return *this; } - __hostdev__ Vec2& operator/=(const T& s) { return (*this) *= T(1) / s; } - __hostdev__ Vec2& normalize() { return (*this) /= this->length(); } - /// @brief Perform a component-wise minimum with the other Coord. - __hostdev__ Vec2& minComponent(const Vec2& other) + + // ---- integer rounding (toward -inf / +inf / nearest) ---- + // + + /// @brief floor-rounded components into the @c Result type (whose @c [i] + /// must accept int32_t). For integer @c T the value is passed through. + /// @note Only the integer-@c T specialization is usable as a constant + /// expression; the floating-point branch calls @c math::Floor, which + /// isn't constexpr until C++23. + template + __hostdev__ [[nodiscard]] constexpr Result floorAs() const noexcept { + Result r{}; + if constexpr (std::is_floating_point::value) { + for (int i = 0; i < N; ++i) r[i] = math::Floor(mVec[i]); + } else { + for (int i = 0; i < N; ++i) r[i] = static_cast(mVec[i]); + } + return r; + } + /// @brief ceil-rounded components into the @c Result type (whose @c [i] + /// must accept int32_t). For integer @c T the value is passed through. + /// @note Only the integer-@c T specialization is usable as a constant + /// expression; the floating-point branch calls @c math::Ceil, which + /// isn't constexpr until C++23. + template + __hostdev__ [[nodiscard]] constexpr Result ceilAs() const noexcept { + Result r{}; + if constexpr (std::is_floating_point::value) { + for (int i = 0; i < N; ++i) r[i] = math::Ceil(mVec[i]); + } else { + for (int i = 0; i < N; ++i) r[i] = static_cast(mVec[i]); + } + return r; + } + /// @brief nearest-integer rounding using floor(x + 0.5) for floating @c T + /// (round-half-toward-positive-infinity). Unifies behaviour between + /// float, double, and long double; pass-through for integer @c T. + /// @note See @c floorAs — only the integer-@c T branch is constexpr-usable. + template + __hostdev__ [[nodiscard]] constexpr Result roundAs() const noexcept { + Result r{}; + if constexpr (std::is_floating_point::value) { + const T half = T(0.5); + for (int i = 0; i < N; ++i) r[i] = math::Floor(mVec[i] + half); + } else { + for (int i = 0; i < N; ++i) r[i] = static_cast(mVec[i]); + } + return r; + } +}; // VecBase + +// ----------------------------> Vec2 <-------------------------------------- + +/// @brief A simple vector class with two components, similar to openvdb::math::Vec2 +/// +/// Aligned to 2*alignof(T) so the whole class fits in one SIMD-friendly +/// chunk (e.g. 8 bytes for Vec2, 16 bytes for Vec2) +template +class alignas(alignof(T) * 2) Vec2 final : public VecBase +{ + using Base = VecBase; + +public: + using ValueType = T; + static constexpr int size = 2; // openvdb::math::Tuple-compat alias of SIZE + + /// @brief Default-construct (components are left uninitialized for fundamental @c T). + Vec2() noexcept = default; + /// @brief Broadcast: set both components to @a x. + __hostdev__ explicit constexpr Vec2(T x) noexcept : Base(x, x) {} + /// @brief Component-wise construction. + __hostdev__ constexpr Vec2(T x, T y) noexcept : Base(x, y) {} + + /// @brief Cross-template converting ctor (e.g. from @c openvdb::Vec2). Implicit + /// to preserve foreign-type interop; same-class ctor below is @c explicit. + template class Vec2T, class T2> + __hostdev__ explicit constexpr Vec2(const Vec2T& v) noexcept : Base(v[0], v[1]) { - if (other[0] < mVec[0]) - mVec[0] = other[0]; - if (other[1] < mVec[1]) - mVec[1] = other[1]; + static_assert(Vec2T::size == 2, "expected Vec2T::size==2!"); + } + /// @brief Explicit cross-precision conversion within nanovdb (e.g. @c Vec2d → @c Vec2f). + template + __hostdev__ explicit constexpr Vec2(const Vec2& v) noexcept : Base(v[0], v[1]) {} + /// @brief Construct from a 2D integer coordinate. + __hostdev__ explicit constexpr Vec2(const Coord2& ijk) noexcept : Base(ijk[0], ijk[1]) {} + + /// @brief Assign from any 2-component vector type (foreign or nanovdb). + template class Vec2T, class T2> + __hostdev__ constexpr Vec2& operator=(const Vec2T& rhs) noexcept { + static_assert(Vec2T::size == 2, "expected Vec2T::size==2!"); + this->mVec[0] = rhs[0]; this->mVec[1] = rhs[1]; return *this; } - /// @brief Perform a component-wise maximum with the other Coord. - __hostdev__ Vec2& maxComponent(const Vec2& other) - { - if (other[0] > mVec[0]) - mVec[0] = other[0]; - if (other[1] > mVec[1]) - mVec[1] = other[1]; + // ---- element-wise (Vec & Vec) ---- + /// @brief Component-wise negation. + __hostdev__ [[nodiscard]] constexpr Vec2 operator-() const noexcept { return Base::template negate(); } + /// @brief Component-wise sum. + __hostdev__ [[nodiscard]] constexpr Vec2 operator+(const Vec2& v) const noexcept { return Base::template plus(v); } + /// @brief Component-wise difference. + __hostdev__ [[nodiscard]] constexpr Vec2 operator-(const Vec2& v) const noexcept { return Base::template minus(v); } + /// @brief Component-wise (Hadamard) product. + __hostdev__ [[nodiscard]] constexpr Vec2 operator*(const Vec2& v) const noexcept { return Base::template mul(v); } + /// @brief Component-wise division. + __hostdev__ [[nodiscard]] constexpr Vec2 operator/(const Vec2& v) const noexcept { return Base::template div(v); } + /// @brief In-place component-wise addition. + __hostdev__ constexpr Vec2& operator+=(const Vec2& v) noexcept { Base::addAssign(v); return *this; } + /// @brief In-place component-wise subtraction. + __hostdev__ constexpr Vec2& operator-=(const Vec2& v) noexcept { Base::subAssign(v); return *this; } + + // ---- mixed Vec2 / Coord2 ---- + /// @brief Add an integer @c Coord2 to this vector (component-wise). + __hostdev__ [[nodiscard]] constexpr Vec2 operator+(const Coord2& ijk) const noexcept { return Vec2(this->mVec[0] + ijk[0], this->mVec[1] + ijk[1]); } + /// @brief Subtract an integer @c Coord2 from this vector (component-wise). + __hostdev__ [[nodiscard]] constexpr Vec2 operator-(const Coord2& ijk) const noexcept { return Vec2(this->mVec[0] - ijk[0], this->mVec[1] - ijk[1]); } + /// @brief In-place component-wise addition of an integer @c Coord2. + __hostdev__ constexpr Vec2& operator+=(const Coord2& ijk) noexcept { + this->mVec[0] += T(ijk[0]); this->mVec[1] += T(ijk[1]); + return *this; + } + /// @brief In-place component-wise subtraction of an integer @c Coord2. + __hostdev__ constexpr Vec2& operator-=(const Coord2& ijk) noexcept { + this->mVec[0] -= T(ijk[0]); this->mVec[1] -= T(ijk[1]); return *this; } + + // ---- scalar ---- + /// @brief Component-wise multiply by scalar @a s. + __hostdev__ [[nodiscard]] constexpr Vec2 operator*(const T& s) const noexcept { return Base::template scale(s); } + /// @brief Component-wise divide by scalar @a s (integer-safe). + __hostdev__ [[nodiscard]] constexpr Vec2 operator/(const T& s) const noexcept { return Base::template divideBy(s); } + /// @brief In-place component-wise multiply by scalar @a s. + __hostdev__ constexpr Vec2& operator*=(const T& s) noexcept { Base::scaleAssign(s); return *this; } + /// @brief In-place component-wise divide by scalar @a s. + __hostdev__ constexpr Vec2& operator/=(const T& s) noexcept { Base::divideAssignScalar(s); return *this; } + /// @brief Normalize in place (divide by @c length()). Not @c constexpr — calls @c std::sqrt. + __hostdev__ Vec2& normalize() noexcept { return (*this) /= this->length(); } + /// @brief Return a normalized (unit-length) copy; const counterpart of @c normalize(). + __hostdev__ [[nodiscard]] Vec2 normalized() const noexcept { + return Base::template normalized(); + } + + // ---- equality ---- + /// @brief Component-wise equality. + __hostdev__ [[nodiscard]] constexpr bool operator==(const Vec2& rhs) const noexcept { return Base::equals(rhs); } + /// @brief Component-wise inequality. + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Vec2& rhs) const noexcept { return !Base::equals(rhs); } + + // ---- component-wise min/max ---- + /// @brief Take the component-wise minimum of @c *this and @a other in place. + __hostdev__ constexpr Vec2& minComponent(const Vec2& other) noexcept { Base::mergeMin(other); return *this; } + /// @brief Take the component-wise maximum of @c *this and @a other in place. + __hostdev__ constexpr Vec2& maxComponent(const Vec2& other) noexcept { Base::mergeMax(other); return *this; } + /// @brief Return the smallest vector component - __hostdev__ ValueType min() const + __hostdev__ [[nodiscard]] constexpr ValueType min() const noexcept { return Base::smallestComponent(); } + /// @brief Return the largest vector component + __hostdev__ [[nodiscard]] constexpr ValueType max() const noexcept { return Base::largestComponent(); } + + /// @brief Round each component down (toward negative infinity) + /// @return integer Coord2 + /// @note Only constexpr for integer @c T (floorAs uses non-constexpr math::Floor for floating point). + __hostdev__ [[nodiscard]] constexpr Coord2 floor() const noexcept { return Base::template floorAs(); } + /// @brief Round each component up (toward positive infinity) + /// @return integer Coord2 + /// @note Only constexpr for integer @c T (ceilAs uses non-constexpr math::Ceil for floating point). + __hostdev__ [[nodiscard]] constexpr Coord2 ceil() const noexcept { return Base::template ceilAs(); } + /// @brief Round each component to its closest integer value + /// @return integer Coord2 + /// @note Only constexpr for integer @c T (roundAs uses non-constexpr math::Floor for floating point). + __hostdev__ [[nodiscard]] constexpr Coord2 round() const noexcept { return Base::template roundAs(); } + + // ---- scalar * Vec / scalar / Vec (hidden friends — found only via ADL on Vec2, + // never participate in unrelated namespace-scope overload sets) ---- + + /// @brief Scalar-on-the-left multiplication (hidden friend). + template + __hostdev__ [[nodiscard]] friend constexpr Vec2 operator*(T1 scalar, const Vec2& vec) noexcept { - return mVec[0] < mVec[1] ? mVec[0] : mVec[1]; + return Vec2(scalar * vec[0], scalar * vec[1]); } - /// @brief Return the largest vector component - __hostdev__ ValueType max() const - { - return mVec[0] > mVec[1] ? mVec[0] : mVec[1]; - } - /// @brief Round each component if this Vec up to its integer value - /// @return Return an integer Coord - __hostdev__ Coord2 floor() const { return Coord2(Floor(mVec[0]), Floor(mVec[1])); } - /// @brief Round each component if this Vec down to its integer value - /// @return Return an integer Coord - __hostdev__ Coord2 ceil() const { return Coord2(Ceil(mVec[0]), Ceil(mVec[1])); } - /// @brief Round each component if this Vec to its closest integer value - /// @return Return an integer Coord - __hostdev__ Coord2 round() const - { - if constexpr(util::is_same::value) { - return Coord2(Floor(mVec[0] + 0.5f), Floor(mVec[1] + 0.5f)); - } else if constexpr(util::is_same::value) { - return Coord2(mVec[0], mVec[1]); - } else { - return Coord2(Floor(mVec[0] + 0.5), Floor(mVec[1] + 0.5)); - } + /// @brief Scalar-on-the-left division (hidden friend). + template + __hostdev__ [[nodiscard]] friend constexpr Vec2 operator/(T1 scalar, const Vec2& vec) noexcept + { + return Vec2(scalar / vec[0], scalar / vec[1]); } - - /// @brief return a non-const raw constant pointer to array of three vector components - __hostdev__ T* asPointer() { return mVec; } - /// @brief return a const raw constant pointer to array of three vector components - __hostdev__ const T* asPointer() const { return mVec; } }; // Vec2 -template -__hostdev__ inline Vec2 operator*(T1 scalar, const Vec2& vec) -{ - return Vec2(scalar * vec[0], scalar * vec[1]); -} -template -__hostdev__ inline Vec2 operator/(T1 scalar, const Vec2& vec) -{ - return Vec2(scalar / vec[0], scalar / vec[1]); -} - /// @brief Return a single precision floating-point vector of this coordinate -__hostdev__ inline Vec2 Coord2::asVec2s() const +__hostdev__ [[nodiscard]] inline constexpr Vec2 Coord2::asVec2s() const noexcept { return Vec2(float(mVec[0]), float(mVec[1])); } /// @brief Return a double precision floating-point vector of this coordinate -__hostdev__ inline Vec2 Coord2::asVec2d() const +__hostdev__ [[nodiscard]] inline constexpr Vec2 Coord2::asVec2d() const noexcept { return Vec2(double(mVec[0]), double(mVec[1])); } -// Matrix base class +/// @brief Base class for fixed-size matrices. Provides shared element-wise +/// arithmetic, scalar arithmetic, equality, transpose, mat*mat, and mat*vec. +/// Storage is a flat row-major array @c mData[ROWS * COLS]. Derived classes +/// supply constructors plus any dimension-specific extras (e.g. @c inverse() +/// on @c Mat2). template class MatBase { protected: T mData[ROWS * COLS]; // 1D array storage - static constexpr int rows() { return ROWS; } - static constexpr int cols() { return COLS; } - static constexpr int size() { return ROWS * COLS; } - public: - MatBase() = default; + using ValueType = T; + + /// @brief Compile-time row count. + [[nodiscard]] static constexpr int rows() noexcept { return ROWS; } + /// @brief Compile-time column count. + [[nodiscard]] static constexpr int cols() noexcept { return COLS; } + /// @brief Compile-time element count, i.e. @c ROWS * @c COLS. + [[nodiscard]] static constexpr int size() noexcept { return ROWS * COLS; } + + /// @brief Default-construct (entries are left uninitialized for trivially-default-constructible @c T). + MatBase() noexcept = default; + /// @brief Read @c ROWS*COLS elements from @a array in row-major order, converting each + /// element from @c S to @c T via @c static_cast. template - __hostdev__ MatBase(S* array) { - for (int i = 0; i < ROWS * COLS; ++i) { + __hostdev__ constexpr MatBase(S* array) noexcept { + for (int i = 0; i < size(); ++i) { mData[i] = static_cast(array[i]); } } - // 2D array access - __hostdev__ T* operator[](int row) { +protected: + /// @brief Variadic component-init ctor — used only by derived ctors + /// (`Mat2(T a, T b, T c, T d) : Base(a, b, c, d) {}`). The member-init + /// list directly initializes @c mData, avoiding the default-construct- + /// then-assign that a body-assignment ctor would impose on a + /// non-fundamental @c T. + template + __hostdev__ explicit constexpr MatBase(Args... args) noexcept : mData{T(args)...} + { + static_assert(sizeof...(Args) == ROWS * COLS, "MatBase: wrong number of constructor arguments"); + } + +public: + + /// @brief 2D row access. Returns a pointer to the start of @a row; the caller + /// indexes into it with the column. Asserts @c 0 <= @a row < @c ROWS in debug + /// builds; column index is the caller's responsibility (use @c 0 <= col < @c COLS). + __hostdev__ constexpr T* operator[](int row) noexcept { + NANOVDB_ASSERT(row >= 0 && row < ROWS); return &mData[row * COLS]; } - __hostdev__ const T* operator[](int row) const { + /// @brief Const overload of @c operator[](int); see non-const variant. + __hostdev__ constexpr const T* operator[](int row) const noexcept { + NANOVDB_ASSERT(row >= 0 && row < ROWS); return &mData[row * COLS]; } + + /// @brief return a raw pointer to the underlying 1D storage (row-major) + __hostdev__ constexpr T* data() noexcept { return mData; } + /// @brief return a const raw pointer to the underlying 1D storage (row-major) + __hostdev__ constexpr const T* data() const noexcept { return mData; } + + // ---- generic element-wise helpers ---- + + /// @brief return @c *this + @a rhs as a @c Derived + template + __hostdev__ [[nodiscard]] constexpr Derived plus(const Derived& rhs) const noexcept { + Derived out{}; + for (int i = 0; i < size(); ++i) out.data()[i] = mData[i] + rhs.data()[i]; + return out; + } + + /// @brief return @c *this - @a rhs as a @c Derived + template + __hostdev__ [[nodiscard]] constexpr Derived minus(const Derived& rhs) const noexcept { + Derived out{}; + for (int i = 0; i < size(); ++i) out.data()[i] = mData[i] - rhs.data()[i]; + return out; + } + + /// @brief return -(*this) as a @c Derived + template + __hostdev__ [[nodiscard]] constexpr Derived negate() const noexcept { + Derived out{}; + for (int i = 0; i < size(); ++i) out.data()[i] = -mData[i]; + return out; + } + + /// @brief return @a s * (*this) as a @c Derived + template + __hostdev__ [[nodiscard]] constexpr Derived scale(const T& s) const noexcept { + Derived out{}; + for (int i = 0; i < size(); ++i) out.data()[i] = mData[i] * s; + return out; + } + + /// @brief Element-wise add @a rhs into @c *this and return @c *this. + __hostdev__ constexpr MatBase& addAssign(const MatBase& rhs) noexcept { + for (int i = 0; i < size(); ++i) mData[i] += rhs.mData[i]; + return *this; + } + /// @brief Element-wise subtract @a rhs from @c *this and return @c *this. + __hostdev__ constexpr MatBase& subAssign(const MatBase& rhs) noexcept { + for (int i = 0; i < size(); ++i) mData[i] -= rhs.mData[i]; + return *this; + } + /// @brief Multiply every element by scalar @a s in place; return @c *this. + __hostdev__ constexpr MatBase& scaleAssign(const T& s) noexcept { + for (int i = 0; i < size(); ++i) mData[i] *= s; + return *this; + } + + /// @brief return (*this) / @a s element-wise as a @c Derived. Uses + /// per-element division (correct for integer @c T, unlike multiplying by 1/s). + template + __hostdev__ [[nodiscard]] constexpr Derived divideBy(const T& s) const noexcept { + Derived out{}; + for (int i = 0; i < size(); ++i) out.data()[i] = mData[i] / s; + return out; + } + /// @brief Divide every element by scalar @a s in place; return @c *this. + __hostdev__ constexpr MatBase& divideAssignScalar(const T& s) noexcept { + for (int i = 0; i < size(); ++i) mData[i] /= s; + return *this; + } + + /// @brief Return @c true iff every element compares equal to @a rhs. + __hostdev__ [[nodiscard]] constexpr bool equals(const MatBase& rhs) const noexcept { + for (int i = 0; i < size(); ++i) if (mData[i] != rhs.mData[i]) return false; + return true; + } + + // ---- generic transpose ---- + + /// @brief return the transpose of @c *this as the @c Result type. + /// @tparam Result a matrix type whose dimensions are (COLS, ROWS) + template + __hostdev__ [[nodiscard]] constexpr Result transposeAs() const noexcept { + static_assert(Result::rows() == COLS && Result::cols() == ROWS, + "transposeAs: result dims must be (COLS, ROWS)"); + Result r{}; + for (int i = 0; i < ROWS; ++i) + for (int j = 0; j < COLS; ++j) + r[j][i] = mData[i * COLS + j]; + return r; + } + + // ---- generic matrix * matrix ---- + + /// @brief return (*this) * @a rhs as the @c Result type. + /// @tparam Result a matrix type whose dimensions are (ROWS, Rhs::cols()) + /// @tparam Rhs a matrix type whose row count equals @c COLS + template + __hostdev__ [[nodiscard]] constexpr Result multiply(const Rhs& rhs) const noexcept { + static_assert(COLS == Rhs::rows(), "multiply: lhs.cols must equal rhs.rows"); + static_assert(Result::rows() == ROWS && Result::cols() == Rhs::cols(), + "multiply: result dims mismatch"); + Result r{}; + for (int i = 0; i < ROWS; ++i) { + for (int j = 0; j < Rhs::cols(); ++j) { + T sum = T(0); + for (int k = 0; k < COLS; ++k) + sum += mData[i * COLS + k] * rhs[k][j]; + r[i][j] = sum; + } + } + return r; + } + + // ---- generic matrix * vector ---- + + /// @brief return (*this) * @a v as the @c VecResult type. + /// @tparam VecResult a vector type whose @c SIZE equals @c ROWS + /// @tparam VecRhs a vector type whose @c SIZE equals @c COLS + template + __hostdev__ [[nodiscard]] constexpr VecResult multiplyVec(const VecRhs& v) const noexcept { + static_assert(VecRhs::SIZE == COLS && VecResult::SIZE == ROWS, + "multiplyVec: dim mismatch"); + VecResult r{}; + for (int i = 0; i < ROWS; ++i) { + T sum = T(0); + for (int k = 0; k < COLS; ++k) sum += mData[i * COLS + k] * v[k]; + r[i] = sum; + } + return r; + } }; // Forward declarations @@ -1011,115 +1503,136 @@ template class Mat2x3; template class Mat3x2; template class Mat3; template class Mat4; +template class Vec4; +/// @brief 2x2 row-major matrix. +/// @details Aligned to 4*alignof(T) — @c Mat2 stores 4 elements (2x2), which is +/// already a power-of-2 multiple of @c alignof(T), so this is free in size +/// and gives SIMD-friendly placement (16 bytes for @c Mat2, 32 bytes +/// for @c Mat2). template -class Mat2 : public MatBase { +class alignas(alignof(T) * 4) Mat2 final : public MatBase { using Base = MatBase; public: - Mat2() = default; + /// @brief Default-construct (entries left uninitialized for fundamental @c T). + Mat2() noexcept = default; /// @brief Constructor given individual array elements, the ordering is in row major form: /** @verbatim a b c d @endverbatim */ - __hostdev__ Mat2(T a, T b, T c, T d) { - this->mData[0] = a; this->mData[1] = b; - this->mData[2] = c; this->mData[3] = d; - } + __hostdev__ constexpr Mat2(T a, T b, T c, T d) noexcept : Base(a, b, c, d) {} /// @brief Constructor given array of elements, the ordering is in row major form template - __hostdev__ Mat2(Source* array) : Base(array) {} - - __hostdev__ Mat2 operator-() const { return Mat2(-(*this)[0][0], -(*this)[0][1], -(*this)[1][0], -(*this)[1][1]); } - - /// @brief Multiply by 2x2 matrix @a m and return the resulting matrix. - __hostdev__ Mat2 operator*(const Mat2& m) const { - return Mat2( - (*this)[0][0] * m[0][0] + (*this)[0][1] * m[1][0], - (*this)[0][0] * m[0][1] + (*this)[0][1] * m[1][1], - (*this)[1][0] * m[0][0] + (*this)[1][1] * m[1][0], - (*this)[1][0] * m[0][1] + (*this)[1][1] * m[1][1] - ); - } - - /// @brief Add each element of the given matrix to the corresponding element of this matrix. - __hostdev__ Mat2& operator+=(const Mat2& m) { - (*this)[0][0] += m[0][0]; - (*this)[0][1] += m[0][1]; - (*this)[1][0] += m[1][0]; - (*this)[1][1] += m[1][1]; - return *this; - } + __hostdev__ constexpr Mat2(Source* array) noexcept : Base(array) {} + + // ---- element-wise ---- + __hostdev__ [[nodiscard]] constexpr Mat2 operator-() const noexcept { return this->template negate(); } + __hostdev__ [[nodiscard]] constexpr Mat2 operator+(const Mat2& m) const noexcept { return this->template plus(m); } + __hostdev__ [[nodiscard]] constexpr Mat2 operator-(const Mat2& m) const noexcept { return this->template minus(m); } + __hostdev__ constexpr Mat2& operator+=(const Mat2& m) noexcept { Base::addAssign(m); return *this; } + __hostdev__ constexpr Mat2& operator-=(const Mat2& m) noexcept { Base::subAssign(m); return *this; } + + // ---- matrix * matrix / matrix * vector ---- + __hostdev__ [[nodiscard]] constexpr Mat2 operator*(const Mat2& m) const noexcept { return this->template multiply(m); } + __hostdev__ [[nodiscard]] constexpr Vec2 operator*(const Vec2& v) const noexcept { return this->template multiplyVec, Vec2>(v); } + + // ---- scalar ---- + __hostdev__ [[nodiscard]] constexpr Mat2 operator*(const T& s) const noexcept { return this->template scale(s); } + __hostdev__ [[nodiscard]] constexpr Mat2 operator/(const T& s) const noexcept { return this->template divideBy(s); } + __hostdev__ constexpr Mat2& operator*=(const T& s) noexcept { Base::scaleAssign(s); return *this; } + __hostdev__ constexpr Mat2& operator/=(const T& s) noexcept { Base::divideAssignScalar(s); return *this; } + + // ---- equality ---- + __hostdev__ [[nodiscard]] constexpr bool operator==(const Mat2& m) const noexcept { return Base::equals(m); } + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Mat2& m) const noexcept { return !Base::equals(m); } /// @brief returns transpose of this - __hostdev__ Mat2 transpose() const { - return Mat2((*this)[0][0], (*this)[1][0], (*this)[0][1], (*this)[1][1]); - } + __hostdev__ [[nodiscard]] constexpr Mat2 transpose() const noexcept { return this->template transposeAs(); } /// @brief returns inverse of this - __hostdev__ Mat2 inverse() const { + __hostdev__ [[nodiscard]] constexpr Mat2 inverse() const noexcept { T det = (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0]; if (isApproxZero(det)) { - return Mat2(); + return Mat2(T(0), T(0), T(0), T(0)); } - T invDet = 1.f / det; - return Mat2((*this)[1][1] * invDet, -(*this)[0][1] * invDet, -(*this)[1][0] * invDet, (*this)[0][0] * invDet); + T invDet = T(1) / det; + return Mat2((*this)[1][1] * invDet, -(*this)[0][1] * invDet, + -(*this)[1][0] * invDet, (*this)[0][0] * invDet); + } + + /// @brief scalar * Mat (hidden friend — found only via ADL on Mat2) + __hostdev__ [[nodiscard]] friend constexpr Mat2 operator*(const T& s, const Mat2& m) noexcept { + return m.template scale(s); } }; +/// @brief 2x3 row-major matrix. +/// @details Intentionally NOT @c alignas-elevated: its byte size +/// (6*sizeof(T)) is not a power-of-2 multiple of @c alignof(T), so any +/// @c alignas(N > alignof(T)) would force tail padding and break +/// packed-array layout plus on-disk format compatibility. template -class Mat2x3 : public MatBase { +class Mat2x3 final : public MatBase { using Base = MatBase; public: - Mat2x3() = default; + /// @brief Default-construct (entries left uninitialized for fundamental @c T). + Mat2x3() noexcept = default; /// @brief Constructor given individual array elements, the ordering is in row major form: /** @verbatim a b c d e f @endverbatim */ - __hostdev__ Mat2x3(T a, T b, T c, T d, T e, T f) { - this->mData[0] = a; this->mData[1] = b; this->mData[2] = c; - this->mData[3] = d; this->mData[4] = e; this->mData[5] = f; - } + __hostdev__ constexpr Mat2x3(T a, T b, T c, T d, T e, T f) noexcept : Base(a, b, c, d, e, f) {} /// @brief Constructor given array of elements, the ordering is in row major form template - __hostdev__ Mat2x3(Source* array) : Base(array) {} + __hostdev__ constexpr Mat2x3(Source* array) noexcept : Base(array) {} - /// @brief Add two matrices and return the resulting matrix. - __hostdev__ Mat2x3 operator+(const Mat2x3& m) const { - return Mat2x3( - (*this)[0][0] + m[0][0], (*this)[0][1] + m[0][1], (*this)[0][2] + m[0][2], - (*this)[1][0] + m[1][0], (*this)[1][1] + m[1][1], (*this)[1][2] + m[1][2] - ); - } + // ---- element-wise ---- + __hostdev__ [[nodiscard]] constexpr Mat2x3 operator-() const noexcept { return this->template negate(); } + __hostdev__ [[nodiscard]] constexpr Mat2x3 operator+(const Mat2x3& m) const noexcept { return this->template plus(m); } + __hostdev__ [[nodiscard]] constexpr Mat2x3 operator-(const Mat2x3& m) const noexcept { return this->template minus(m); } + __hostdev__ constexpr Mat2x3& operator+=(const Mat2x3& m) noexcept { Base::addAssign(m); return *this; } + __hostdev__ constexpr Mat2x3& operator-=(const Mat2x3& m) noexcept { Base::subAssign(m); return *this; } - /// @brief Add a 2x3 matrix to this matrix. - __hostdev__ Mat2x3& operator+=(const Mat2x3& m) { - (*this)[0][0] += m[0][0]; (*this)[0][1] += m[0][1]; (*this)[0][2] += m[0][2]; - (*this)[1][0] += m[1][0]; (*this)[1][1] += m[1][1]; (*this)[1][2] += m[1][2]; - return *this; - } + // ---- matrix * vector ---- + __hostdev__ [[nodiscard]] constexpr Vec2 operator*(const Vec3& v) const noexcept { return this->template multiplyVec, Vec3>(v); } + + // ---- scalar ---- + __hostdev__ [[nodiscard]] constexpr Mat2x3 operator*(const T& s) const noexcept { return this->template scale(s); } + __hostdev__ [[nodiscard]] constexpr Mat2x3 operator/(const T& s) const noexcept { return this->template divideBy(s); } + __hostdev__ constexpr Mat2x3& operator*=(const T& s) noexcept { Base::scaleAssign(s); return *this; } + __hostdev__ constexpr Mat2x3& operator/=(const T& s) noexcept { Base::divideAssignScalar(s); return *this; } + + // ---- equality ---- + __hostdev__ [[nodiscard]] constexpr bool operator==(const Mat2x3& m) const noexcept { return Base::equals(m); } + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Mat2x3& m) const noexcept { return !Base::equals(m); } /// @brief returns transpose of this - __hostdev__ Mat3x2 transpose() const { - return Mat3x2( - (*this)[0][0], (*this)[1][0], // First row - (*this)[0][1], (*this)[1][1], // Second row - (*this)[0][2], (*this)[1][2] // Third row - ); + __hostdev__ [[nodiscard]] constexpr Mat3x2 transpose() const noexcept { return this->template transposeAs>(); } + + /// @brief scalar * Mat (hidden friend — found only via ADL on Mat2x3) + __hostdev__ [[nodiscard]] friend constexpr Mat2x3 operator*(const T& s, const Mat2x3& m) noexcept { + return m.template scale(s); } }; +/// @brief 3x2 row-major matrix. +/// @details Intentionally NOT @c alignas-elevated: its byte size +/// (6*sizeof(T)) is not a power-of-2 multiple of @c alignof(T), so any +/// @c alignas(N > alignof(T)) would force tail padding and break +/// packed-array layout plus on-disk format compatibility. template -class Mat3x2: public MatBase { +class Mat3x2 final : public MatBase { using Base = MatBase; public: - Mat3x2() = default; + /// @brief Default-construct (entries left uninitialized for fundamental @c T). + Mat3x2() noexcept = default; /// @brief Constructor given individual array elements, the ordering is in row major form: /** @verbatim @@ -1127,31 +1640,52 @@ class Mat3x2: public MatBase { c d e f @endverbatim */ - template - __hostdev__ Mat3x2(Source a, Source b, Source c, Source d, Source e, Source f) - { - this->mData[0] = a; this->mData[1] = b; - this->mData[2] = c; this->mData[3] = d; - this->mData[4] = e; this->mData[5] = f; - } + __hostdev__ constexpr Mat3x2(T a, T b, T c, T d, T e, T f) noexcept : Base(a, b, c, d, e, f) {} /// @brief Constructor given array of elements, the ordering is in row major form template - __hostdev__ Mat3x2(Source *a): Base(a) {} + __hostdev__ constexpr Mat3x2(Source *a) noexcept : Base(a) {} + + // ---- element-wise ---- + __hostdev__ [[nodiscard]] constexpr Mat3x2 operator-() const noexcept { return this->template negate(); } + __hostdev__ [[nodiscard]] constexpr Mat3x2 operator+(const Mat3x2& m) const noexcept { return this->template plus(m); } + __hostdev__ [[nodiscard]] constexpr Mat3x2 operator-(const Mat3x2& m) const noexcept { return this->template minus(m); } + __hostdev__ constexpr Mat3x2& operator+=(const Mat3x2& m) noexcept { Base::addAssign(m); return *this; } + __hostdev__ constexpr Mat3x2& operator-=(const Mat3x2& m) noexcept { Base::subAssign(m); return *this; } + + // ---- matrix * vector ---- + __hostdev__ [[nodiscard]] constexpr Vec3 operator*(const Vec2& v) const noexcept { return this->template multiplyVec, Vec2>(v); } + + // ---- scalar ---- + __hostdev__ [[nodiscard]] constexpr Mat3x2 operator*(const T& s) const noexcept { return this->template scale(s); } + __hostdev__ [[nodiscard]] constexpr Mat3x2 operator/(const T& s) const noexcept { return this->template divideBy(s); } + __hostdev__ constexpr Mat3x2& operator*=(const T& s) noexcept { Base::scaleAssign(s); return *this; } + __hostdev__ constexpr Mat3x2& operator/=(const T& s) noexcept { Base::divideAssignScalar(s); return *this; } + + // ---- equality ---- + __hostdev__ [[nodiscard]] constexpr bool operator==(const Mat3x2& m) const noexcept { return Base::equals(m); } + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Mat3x2& m) const noexcept { return !Base::equals(m); } /// @brief returns transpose of this - __hostdev__ Mat2x3 transpose() const { - return Mat2x3((*this)[0][0], (*this)[1][0], (*this)[2][0], // First row - (*this)[0][1], (*this)[1][1], (*this)[2][1]); // Second row - } + __hostdev__ [[nodiscard]] constexpr Mat2x3 transpose() const noexcept { return this->template transposeAs>(); } + /// @brief scalar * Mat (hidden friend — found only via ADL on Mat3x2) + __hostdev__ [[nodiscard]] friend constexpr Mat3x2 operator*(const T& s, const Mat3x2& m) noexcept { + return m.template scale(s); + } }; +/// @brief 3x3 row-major matrix. +/// @details Intentionally NOT @c alignas-elevated: its byte size +/// (9*sizeof(T)) is not a power-of-2 multiple of @c alignof(T), so any +/// @c alignas(N > alignof(T)) would force tail padding and break +/// packed-array layout plus on-disk format compatibility. template -class Mat3 : public MatBase { +class Mat3 final : public MatBase { using Base = MatBase; public: - Mat3() = default; + /// @brief Default-construct (entries left uninitialized for fundamental @c T). + Mat3() noexcept = default; /// @brief Constructor given individual array elements, the ordering is in row major form: /** @verbatim @@ -1159,75 +1693,59 @@ class Mat3 : public MatBase { d e f g h i @endverbatim */ - template - __hostdev__ Mat3(Source a, Source b, Source c, - Source d, Source e, Source f, - Source g, Source h, Source i) - { - this->mData[0] = a; this->mData[1] = b; this->mData[2] = c; - this->mData[3] = d; this->mData[4] = e; this->mData[5] = f; - this->mData[6] = g; this->mData[7] = h; this->mData[8] = i; - } + __hostdev__ constexpr Mat3(T a, T b, T c, + T d, T e, T f, + T g, T h, T i) noexcept : Base(a, b, c, d, e, f, g, h, i) {} /// @brief Constructor given array of elements, the ordering is in row major form template - __hostdev__ Mat3(Source *a): Base(a) {} - - - /// @brief Add two matrices and return the resulting matrix. - __hostdev__ Mat3 operator+(const Mat3& m) const { - return Mat3( - (*this)[0][0] + m[0][0], (*this)[0][1] + m[0][1], (*this)[0][2] + m[0][2], - (*this)[1][0] + m[1][0], (*this)[1][1] + m[1][1], (*this)[1][2] + m[1][2], - (*this)[2][0] + m[2][0], (*this)[2][1] + m[2][1], (*this)[2][2] + m[2][2] - ); - } - - /// @brief Multiply by @a v and return the resulting vector. - __hostdev__ Vec3 operator*(const Vec3& v) const { - return Vec3( - (*this)[0][0] * v[0] + (*this)[0][1] * v[1] + (*this)[0][2] * v[2], - (*this)[1][0] * v[0] + (*this)[1][1] * v[1] + (*this)[1][2] * v[2], - (*this)[2][0] * v[0] + (*this)[2][1] * v[1] + (*this)[2][2] * v[2] - ); - } - - /// @brief Multiply by 3x3 matrix @a m and return the resulting matrix. - __hostdev__ Mat3 operator*(const Mat3& m) const { - return Mat3( - (*this)[0][0] * m[0][0] + (*this)[0][1] * m[1][0] + (*this)[0][2] * m[2][0], - (*this)[0][0] * m[0][1] + (*this)[0][1] * m[1][1] + (*this)[0][2] * m[2][1], - (*this)[0][0] * m[0][2] + (*this)[0][1] * m[1][2] + (*this)[0][2] * m[2][2], - (*this)[1][0] * m[0][0] + (*this)[1][1] * m[1][0] + (*this)[1][2] * m[2][0], - (*this)[1][0] * m[0][1] + (*this)[1][1] * m[1][1] + (*this)[1][2] * m[2][1], - (*this)[1][0] * m[0][2] + (*this)[1][1] * m[1][2] + (*this)[1][2] * m[2][2], - (*this)[2][0] * m[0][0] + (*this)[2][1] * m[1][0] + (*this)[2][2] * m[2][0], - (*this)[2][0] * m[0][1] + (*this)[2][1] * m[1][1] + (*this)[2][2] * m[2][1], - (*this)[2][0] * m[0][2] + (*this)[2][1] * m[1][2] + (*this)[2][2] * m[2][2] - ); - } - - /// @brief Add each element of the given matrix to the corresponding element of this matrix. - __hostdev__ Mat3& operator+=(const Mat3& m) { - (*this)[0][0] += m[0][0]; (*this)[0][1] += m[0][1]; (*this)[0][2] += m[0][2]; - (*this)[1][0] += m[1][0]; (*this)[1][1] += m[1][1]; (*this)[1][2] += m[1][2]; - (*this)[2][0] += m[2][0]; (*this)[2][1] += m[2][1]; (*this)[2][2] += m[2][2]; - return *this; - } + __hostdev__ constexpr Mat3(Source *a) noexcept : Base(a) {} + + + // ---- element-wise ---- + __hostdev__ [[nodiscard]] constexpr Mat3 operator-() const noexcept { return this->template negate(); } + __hostdev__ [[nodiscard]] constexpr Mat3 operator+(const Mat3& m) const noexcept { return this->template plus(m); } + __hostdev__ [[nodiscard]] constexpr Mat3 operator-(const Mat3& m) const noexcept { return this->template minus(m); } + __hostdev__ constexpr Mat3& operator+=(const Mat3& m) noexcept { Base::addAssign(m); return *this; } + __hostdev__ constexpr Mat3& operator-=(const Mat3& m) noexcept { Base::subAssign(m); return *this; } + + // ---- matrix * matrix / matrix * vector ---- + __hostdev__ [[nodiscard]] constexpr Mat3 operator*(const Mat3& m) const noexcept { return this->template multiply(m); } + __hostdev__ [[nodiscard]] constexpr Vec3 operator*(const Vec3& v) const noexcept { return this->template multiplyVec, Vec3>(v); } + + // ---- scalar ---- + __hostdev__ [[nodiscard]] constexpr Mat3 operator*(const T& s) const noexcept { return this->template scale(s); } + __hostdev__ [[nodiscard]] constexpr Mat3 operator/(const T& s) const noexcept { return this->template divideBy(s); } + __hostdev__ constexpr Mat3& operator*=(const T& s) noexcept { Base::scaleAssign(s); return *this; } + __hostdev__ constexpr Mat3& operator/=(const T& s) noexcept { Base::divideAssignScalar(s); return *this; } + + // ---- equality ---- + __hostdev__ [[nodiscard]] constexpr bool operator==(const Mat3& m) const noexcept { return Base::equals(m); } + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Mat3& m) const noexcept { return !Base::equals(m); } /// @brief returns transpose of this - __hostdev__ Mat3 transpose() const { - return Mat3((*this)[0][0], (*this)[1][0], (*this)[2][0], - (*this)[0][1], (*this)[1][1], (*this)[2][1], - (*this)[0][2], (*this)[1][2], (*this)[2][2]); + __hostdev__ [[nodiscard]] constexpr Mat3 transpose() const noexcept { return this->template transposeAs(); } + + /// @brief scalar * Mat (hidden friend — found only via ADL on Mat3) + __hostdev__ [[nodiscard]] friend constexpr Mat3 operator*(const T& s, const Mat3& m) noexcept { + return m.template scale(s); } }; +/// @brief 4x4 row-major matrix. +/// @details Aligned to 16*alignof(T) — @c Mat4 stores 16 elements (4x4), which is +/// already a power-of-2 multiple of @c alignof(T), so the alignment is free +/// in size. Whole-matrix alignment (64 bytes for @c Mat4, 128 bytes +/// for @c Mat4) is heavy compared to row-alignment (4*alignof(T)) +/// but matches the per-class "align to full size" rule used for @c Vec2 / +/// @c Vec4 / @c Mat2 above and lets a @c Mat4 load with a single +/// AVX-512 instruction. template -class Mat4 : public MatBase { +class alignas(alignof(T) * 16) Mat4 final : public MatBase { using Base = MatBase; public: - Mat4() = default; + /// @brief Default-construct (entries left uninitialized for fundamental @c T). + Mat4() noexcept = default; /// @brief Constructor given individual array elements, the ordering is in row major form: /** @verbatim @@ -1236,318 +1754,243 @@ class Mat4 : public MatBase { i j k l m n o p @endverbatim */ - template - __hostdev__ Mat4(Source a, Source b, Source c, Source d, - Source e, Source f, Source g, Source h, - Source i, Source j, Source k, Source l, - Source m, Source n, Source o, Source p) - { - this->mData[0] = a; this->mData[1] = b; this->mData[2] = c; this->mData[3] = d; - this->mData[4] = e; this->mData[5] = f; this->mData[6] = g; this->mData[7] = h; - this->mData[8] = i; this->mData[9] = j; this->mData[10] = k; this->mData[11] = l; - this->mData[12] = m; this->mData[13] = n; this->mData[14] = o; this->mData[15] = p; - } + __hostdev__ constexpr Mat4(T a, T b, T c, T d, + T e, T f, T g, T h, + T i, T j, T k, T l, + T m, T n, T o, T p) noexcept + : Base(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p) {} /// @brief Constructor given array of elements, the ordering is in row major form template - __hostdev__ Mat4(Source *a): Base(a) {} + __hostdev__ constexpr Mat4(Source *a) noexcept : Base(a) {} + + // ---- element-wise ---- + __hostdev__ [[nodiscard]] constexpr Mat4 operator-() const noexcept { return this->template negate(); } + __hostdev__ [[nodiscard]] constexpr Mat4 operator+(const Mat4& m) const noexcept { return this->template plus(m); } + __hostdev__ [[nodiscard]] constexpr Mat4 operator-(const Mat4& m) const noexcept { return this->template minus(m); } + __hostdev__ constexpr Mat4& operator+=(const Mat4& m) noexcept { Base::addAssign(m); return *this; } + __hostdev__ constexpr Mat4& operator-=(const Mat4& m) noexcept { Base::subAssign(m); return *this; } + + // ---- matrix * matrix / matrix * vector ---- + __hostdev__ [[nodiscard]] constexpr Mat4 operator*(const Mat4& m) const noexcept { return this->template multiply(m); } + __hostdev__ [[nodiscard]] constexpr Vec4 operator*(const Vec4& v) const noexcept { return this->template multiplyVec, Vec4>(v); } + + // ---- scalar ---- + __hostdev__ [[nodiscard]] constexpr Mat4 operator*(const T& s) const noexcept { return this->template scale(s); } + __hostdev__ [[nodiscard]] constexpr Mat4 operator/(const T& s) const noexcept { return this->template divideBy(s); } + __hostdev__ constexpr Mat4& operator*=(const T& s) noexcept { Base::scaleAssign(s); return *this; } + __hostdev__ constexpr Mat4& operator/=(const T& s) noexcept { Base::divideAssignScalar(s); return *this; } + + // ---- equality ---- + __hostdev__ [[nodiscard]] constexpr bool operator==(const Mat4& m) const noexcept { return Base::equals(m); } + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Mat4& m) const noexcept { return !Base::equals(m); } /// @brief returns transpose of this - __hostdev__ Mat4 transpose() const { - return Mat4((*this)[0][0], (*this)[1][0], (*this)[2][0], (*this)[3][0], - (*this)[0][1], (*this)[1][1], (*this)[2][1], (*this)[3][1], - (*this)[0][2], (*this)[1][2], (*this)[2][2], (*this)[3][2], - (*this)[0][3], (*this)[1][3], (*this)[2][3], (*this)[3][3]); + __hostdev__ [[nodiscard]] constexpr Mat4 transpose() const noexcept { return this->template transposeAs(); } + + /// @brief scalar * Mat (hidden friend — found only via ADL on Mat4) + __hostdev__ [[nodiscard]] friend constexpr Mat4 operator*(const T& s, const Mat4& m) noexcept { + return m.template scale(s); } }; -/// @brief Multiply a scalar by a 2x2 matrix, result is a 2x2 matrix -template -__hostdev__ Mat2 operator*(const T& s, const Mat2& m) { - return Mat2(m[0][0] * s, m[0][1] * s, m[1][0] * s, m[1][1] * s); -} /// @brief Multiply a 2x3 matrix by a 3x2 matrix, result is a 2x2 matrix template -__hostdev__ Mat2 operator*(const Mat2x3& lhs, const Mat3x2& rhs) { - return Mat2( - // First row - lhs[0][0] * rhs[0][0] + lhs[0][1] * rhs[1][0] + lhs[0][2] * rhs[2][0], // [0][0] - lhs[0][0] * rhs[0][1] + lhs[0][1] * rhs[1][1] + lhs[0][2] * rhs[2][1], // [0][1] - - // Second row - lhs[1][0] * rhs[0][0] + lhs[1][1] * rhs[1][0] + lhs[1][2] * rhs[2][0], // [1][0] - lhs[1][0] * rhs[0][1] + lhs[1][1] * rhs[1][1] + lhs[1][2] * rhs[2][1] // [1][1] - ); -} -/// @brief Multiply a 3x3 matrix by a 2x3 matrix, result is a 2x3 matrix -template -__hostdev__ Mat2x3 operator*(const Mat3& lhs, const Mat2x3& rhs) { - return Mat2x3( - // First row - lhs[0][0] * rhs[0][0] + lhs[0][1] * rhs[1][0], - lhs[0][0] * rhs[0][1] + lhs[0][1] * rhs[1][1], - lhs[0][0] * rhs[0][2] + lhs[0][1] * rhs[1][2], - - // Second row - lhs[1][0] * rhs[0][0] + lhs[1][1] * rhs[1][0], - lhs[1][0] * rhs[0][1] + lhs[1][1] * rhs[1][1], - lhs[1][0] * rhs[0][2] + lhs[1][1] * rhs[1][2] - ); +__hostdev__ [[nodiscard]] constexpr Mat2 operator*(const Mat2x3& lhs, const Mat3x2& rhs) noexcept { + return lhs.template multiply, Mat3x2>(rhs); } /// @brief Multiply a 2x3 matrix by a 3x3 matrix, result is a 2x3 matrix template -__hostdev__ Mat2x3 operator*(const Mat2x3& lhs, const Mat3& rhs) { - return Mat2x3( - // First row (3 elements) - lhs[0][0] * rhs[0][0] + lhs[0][1] * rhs[1][0] + lhs[0][2] * rhs[2][0], - lhs[0][0] * rhs[0][1] + lhs[0][1] * rhs[1][1] + lhs[0][2] * rhs[2][1], - lhs[0][0] * rhs[0][2] + lhs[0][1] * rhs[1][2] + lhs[0][2] * rhs[2][2], - - // Second row (3 elements) - lhs[1][0] * rhs[0][0] + lhs[1][1] * rhs[1][0] + lhs[1][2] * rhs[2][0], - lhs[1][0] * rhs[0][1] + lhs[1][1] * rhs[1][1] + lhs[1][2] * rhs[2][1], - lhs[1][0] * rhs[0][2] + lhs[1][1] * rhs[1][2] + lhs[1][2] * rhs[2][2] - ); +__hostdev__ [[nodiscard]] constexpr Mat2x3 operator*(const Mat2x3& lhs, const Mat3& rhs) noexcept { + return lhs.template multiply, Mat3>(rhs); } /// @brief Multiply a 3x2 matrix by a 2x2 matrix, result is a 3x2 matrix template -__hostdev__ Mat3x2 operator*(const Mat3x2& lhs, const Mat2& rhs) { - return Mat3x2( - lhs[0][0] * rhs[0][0] + lhs[0][1] * rhs[1][0], - lhs[0][0] * rhs[0][1] + lhs[0][1] * rhs[1][1], - lhs[1][0] * rhs[0][0] + lhs[1][1] * rhs[1][0], - lhs[1][0] * rhs[0][1] + lhs[1][1] * rhs[1][1], - lhs[2][0] * rhs[0][0] + lhs[2][1] * rhs[1][0], - lhs[2][0] * rhs[0][1] + lhs[2][1] * rhs[1][1] - ); +__hostdev__ [[nodiscard]] constexpr Mat3x2 operator*(const Mat3x2& lhs, const Mat2& rhs) noexcept { + return lhs.template multiply, Mat2>(rhs); } /// @brief Multiply a 2x2 matrix by a 2x3 matrix, result is a 2x3 matrix template -__hostdev__ Mat2x3 operator*(const Mat2& lhs, const Mat2x3& rhs) { - return Mat2x3( - // First row (3 elements) - lhs[0][0] * rhs[0][0] + lhs[0][1] * rhs[1][0], - lhs[0][0] * rhs[0][1] + lhs[0][1] * rhs[1][1], - lhs[0][0] * rhs[0][2] + lhs[0][1] * rhs[1][2], - - // Second row (3 elements) - lhs[1][0] * rhs[0][0] + lhs[1][1] * rhs[1][0], - lhs[1][0] * rhs[0][1] + lhs[1][1] * rhs[1][1], - lhs[1][0] * rhs[0][2] + lhs[1][1] * rhs[1][2] - ); +__hostdev__ [[nodiscard]] constexpr Mat2x3 operator*(const Mat2& lhs, const Mat2x3& rhs) noexcept { + return lhs.template multiply, Mat2x3>(rhs); } /// @brief Multiply a 3x2 matrix by a 2x3 matrix, result is a 3x3 matrix template -__hostdev__ Mat3 operator*(const Mat3x2& lhs, const Mat2x3& rhs) { - return Mat3( - lhs[0][0] * rhs[0][0] + lhs[0][1] * rhs[1][0], - lhs[0][0] * rhs[0][1] + lhs[0][1] * rhs[1][1], - lhs[0][0] * rhs[0][2] + lhs[0][1] * rhs[1][2], - - lhs[1][0] * rhs[0][0] + lhs[1][1] * rhs[1][0], - lhs[1][0] * rhs[0][1] + lhs[1][1] * rhs[1][1], - lhs[1][0] * rhs[0][2] + lhs[1][1] * rhs[1][2], - - lhs[2][0] * rhs[0][0] + lhs[2][1] * rhs[1][0], - lhs[2][0] * rhs[0][1] + lhs[2][1] * rhs[1][1], - lhs[2][0] * rhs[0][2] + lhs[2][1] * rhs[1][2] - ); +__hostdev__ [[nodiscard]] constexpr Mat3 operator*(const Mat3x2& lhs, const Mat2x3& rhs) noexcept { + return lhs.template multiply, Mat2x3>(rhs); +} +/// @brief Multiply a 3x3 matrix by a 3x2 matrix, result is a 3x2 matrix +template +__hostdev__ [[nodiscard]] constexpr Mat3x2 operator*(const Mat3& lhs, const Mat3x2& rhs) noexcept { + return lhs.template multiply, Mat3x2>(rhs); } // ----------------------------> Vec3 <-------------------------------------- /// @brief A simple vector class with three components, similar to openvdb::math::Vec3 +/// +/// Vec3 is intentionally NOT alignas-elevated: its byte size +/// (3*sizeof(T)) is not a power-of-2 multiple of alignof(T), so any +/// alignas(N > alignof(T)) would force tail padding and break +/// packed-array layout plus on-disk format compatibility. template -class Vec3 +class Vec3 final : public VecBase { - T mVec[3]; + using Base = VecBase; public: - static const int SIZE = 3; - static const int size = 3; // in openvdb::math::Tuple using ValueType = T; - Vec3() = default; - __hostdev__ explicit Vec3(T x) - : mVec{x, x, x} - { - } - __hostdev__ Vec3(T x, T y, T z) - : mVec{x, y, z} - { - } + static constexpr int size = 3; // openvdb::math::Tuple-compat alias of SIZE + + /// @brief Default-construct (components are left uninitialized for fundamental @c T). + Vec3() noexcept = default; + /// @brief Broadcast: set all three components to @a x. + __hostdev__ explicit constexpr Vec3(T x) noexcept : Base(x, x, x) {} + /// @brief Component-wise construction. + __hostdev__ constexpr Vec3(T x, T y, T z) noexcept : Base(x, y, z) {} + + /// @brief Cross-template converting ctor (e.g. from @c openvdb::Vec3). Implicit + /// to preserve foreign-type interop; same-class ctor below is @c explicit. template class Vec3T, class T2> - __hostdev__ Vec3(const Vec3T& v) - : mVec{T(v[0]), T(v[1]), T(v[2])} + __hostdev__ explicit constexpr Vec3(const Vec3T& v) noexcept : Base(v[0], v[1], v[2]) { - static_assert(Vec3T::size == size, "expected Vec3T::size==3!"); + static_assert(Vec3T::size == 3, "expected Vec3T::size==3!"); } + /// @brief Explicit cross-precision conversion within nanovdb (e.g. @c Vec3d → @c Vec3f). template - __hostdev__ explicit Vec3(const Vec3& v) - : mVec{T(v[0]), T(v[1]), T(v[2])} - { - } - __hostdev__ explicit Vec3(const Coord& ijk) - : mVec{T(ijk[0]), T(ijk[1]), T(ijk[2])} - { - } - __hostdev__ bool operator==(const Vec3& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2]; } - __hostdev__ bool operator!=(const Vec3& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2]; } + __hostdev__ explicit constexpr Vec3(const Vec3& v) noexcept : Base(v[0], v[1], v[2]) {} + /// @brief Construct from a 3D integer coordinate. + __hostdev__ explicit constexpr Vec3(const Coord& ijk) noexcept : Base(ijk[0], ijk[1], ijk[2]) {} + + /// @brief Assign from any 3-component vector type (foreign or nanovdb). template class Vec3T, class T2> - __hostdev__ Vec3& operator=(const Vec3T& rhs) - { - static_assert(Vec3T::size == size, "expected Vec3T::size==3!"); - mVec[0] = rhs[0]; - mVec[1] = rhs[1]; - mVec[2] = rhs[2]; + __hostdev__ constexpr Vec3& operator=(const Vec3T& rhs) noexcept { + static_assert(Vec3T::size == 3, "expected Vec3T::size==3!"); + this->mVec[0] = rhs[0]; this->mVec[1] = rhs[1]; this->mVec[2] = rhs[2]; return *this; } - __hostdev__ const T& operator[](int i) const { return mVec[i]; } - __hostdev__ T& operator[](int i) { return mVec[i]; } - template - __hostdev__ T dot(const Vec3T& v) const { return mVec[0] * v[0] + mVec[1] * v[1] + mVec[2] * v[2]; } - template - __hostdev__ Vec3 cross(const Vec3T& v) const - { - return Vec3(mVec[1] * v[2] - mVec[2] * v[1], - mVec[2] * v[0] - mVec[0] * v[2], - mVec[0] * v[1] - mVec[1] * v[0]); - } - /// @brief Outer product of a 3x1 vector and a 1x3 vector, result is a 3x3 matrix - template - __hostdev__ Mat3 outer(const Vec3T& v) const - { - return Mat3(mVec[0] * v[0], mVec[0] * v[1], mVec[0] * v[2], - mVec[1] * v[0], mVec[1] * v[1], mVec[1] * v[2], - mVec[2] * v[0], mVec[2] * v[1], mVec[2] * v[2]); - } - __hostdev__ T lengthSqr() const - { - return mVec[0] * mVec[0] + mVec[1] * mVec[1] + mVec[2] * mVec[2]; // 5 flops - } - __hostdev__ T length() const { return Sqrt(this->lengthSqr()); } - __hostdev__ Vec3 operator-() const { return Vec3(-mVec[0], -mVec[1], -mVec[2]); } - __hostdev__ Vec3 operator*(const Vec3& v) const { return Vec3(mVec[0] * v[0], mVec[1] * v[1], mVec[2] * v[2]); } - __hostdev__ Vec3 operator/(const Vec3& v) const { return Vec3(mVec[0] / v[0], mVec[1] / v[1], mVec[2] / v[2]); } - __hostdev__ Vec3 operator+(const Vec3& v) const { return Vec3(mVec[0] + v[0], mVec[1] + v[1], mVec[2] + v[2]); } - __hostdev__ Vec3 operator-(const Vec3& v) const { return Vec3(mVec[0] - v[0], mVec[1] - v[1], mVec[2] - v[2]); } - __hostdev__ Vec3 operator+(const Coord& ijk) const { return Vec3(mVec[0] + ijk[0], mVec[1] + ijk[1], mVec[2] + ijk[2]); } - __hostdev__ Vec3 operator-(const Coord& ijk) const { return Vec3(mVec[0] - ijk[0], mVec[1] - ijk[1], mVec[2] - ijk[2]); } - __hostdev__ Vec3 operator*(const T& s) const { return Vec3(s * mVec[0], s * mVec[1], s * mVec[2]); } - __hostdev__ Vec3 operator/(const T& s) const { return (T(1) / s) * (*this); } - __hostdev__ Vec3& operator+=(const Vec3& v) - { - mVec[0] += v[0]; - mVec[1] += v[1]; - mVec[2] += v[2]; - return *this; - } - __hostdev__ Vec3& operator+=(const Coord& ijk) - { - mVec[0] += T(ijk[0]); - mVec[1] += T(ijk[1]); - mVec[2] += T(ijk[2]); - return *this; - } - __hostdev__ Vec3& operator-=(const Vec3& v) - { - mVec[0] -= v[0]; - mVec[1] -= v[1]; - mVec[2] -= v[2]; + + // ---- element-wise (Vec & Vec) ---- + /// @brief Component-wise negation. + __hostdev__ [[nodiscard]] constexpr Vec3 operator-() const noexcept { return Base::template negate(); } + /// @brief Component-wise sum. + __hostdev__ [[nodiscard]] constexpr Vec3 operator+(const Vec3& v) const noexcept { return Base::template plus(v); } + /// @brief Component-wise difference. + __hostdev__ [[nodiscard]] constexpr Vec3 operator-(const Vec3& v) const noexcept { return Base::template minus(v); } + /// @brief Component-wise (Hadamard) product. + __hostdev__ [[nodiscard]] constexpr Vec3 operator*(const Vec3& v) const noexcept { return Base::template mul(v); } + /// @brief Component-wise division. + __hostdev__ [[nodiscard]] constexpr Vec3 operator/(const Vec3& v) const noexcept { return Base::template div(v); } + /// @brief In-place component-wise addition. + __hostdev__ constexpr Vec3& operator+=(const Vec3& v) noexcept { Base::addAssign(v); return *this; } + /// @brief In-place component-wise subtraction. + __hostdev__ constexpr Vec3& operator-=(const Vec3& v) noexcept { Base::subAssign(v); return *this; } + + // ---- mixed Vec3 / Coord (3D) ---- + /// @brief Add an integer @c Coord to this vector (component-wise). + __hostdev__ [[nodiscard]] constexpr Vec3 operator+(const Coord& ijk) const noexcept { return Vec3(this->mVec[0] + ijk[0], this->mVec[1] + ijk[1], this->mVec[2] + ijk[2]); } + /// @brief Subtract an integer @c Coord from this vector (component-wise). + __hostdev__ [[nodiscard]] constexpr Vec3 operator-(const Coord& ijk) const noexcept { return Vec3(this->mVec[0] - ijk[0], this->mVec[1] - ijk[1], this->mVec[2] - ijk[2]); } + /// @brief In-place component-wise addition of an integer @c Coord. + __hostdev__ constexpr Vec3& operator+=(const Coord& ijk) noexcept { + this->mVec[0] += T(ijk[0]); this->mVec[1] += T(ijk[1]); this->mVec[2] += T(ijk[2]); return *this; } - __hostdev__ Vec3& operator-=(const Coord& ijk) - { - mVec[0] -= T(ijk[0]); - mVec[1] -= T(ijk[1]); - mVec[2] -= T(ijk[2]); + /// @brief In-place component-wise subtraction of an integer @c Coord. + __hostdev__ constexpr Vec3& operator-=(const Coord& ijk) noexcept { + this->mVec[0] -= T(ijk[0]); this->mVec[1] -= T(ijk[1]); this->mVec[2] -= T(ijk[2]); return *this; } - __hostdev__ Vec3& operator*=(const T& s) - { - mVec[0] *= s; - mVec[1] *= s; - mVec[2] *= s; - return *this; + + // ---- scalar ---- + /// @brief Component-wise multiply by scalar @a s. + __hostdev__ [[nodiscard]] constexpr Vec3 operator*(const T& s) const noexcept { return Base::template scale(s); } + /// @brief Component-wise divide by scalar @a s (integer-safe). + __hostdev__ [[nodiscard]] constexpr Vec3 operator/(const T& s) const noexcept { return Base::template divideBy(s); } + /// @brief In-place component-wise multiply by scalar @a s. + __hostdev__ constexpr Vec3& operator*=(const T& s) noexcept { Base::scaleAssign(s); return *this; } + /// @brief In-place component-wise divide by scalar @a s. + __hostdev__ constexpr Vec3& operator/=(const T& s) noexcept { Base::divideAssignScalar(s); return *this; } + /// @brief Normalize in place (divide by @c length()). Not @c constexpr — calls @c std::sqrt. + __hostdev__ Vec3& normalize() noexcept { return (*this) /= this->length(); } + /// @brief Return a normalized (unit-length) copy; const counterpart of @c normalize(). + __hostdev__ [[nodiscard]] Vec3 normalized() const noexcept { + return Base::template normalized(); + } + + // ---- equality ---- + /// @brief Component-wise equality. + __hostdev__ [[nodiscard]] constexpr bool operator==(const Vec3& rhs) const noexcept { return Base::equals(rhs); } + /// @brief Component-wise inequality. + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Vec3& rhs) const noexcept { return !Base::equals(rhs); } + + // ---- component-wise min/max ---- + /// @brief Take the component-wise minimum of @c *this and @a other in place. + __hostdev__ constexpr Vec3& minComponent(const Vec3& other) noexcept { Base::mergeMin(other); return *this; } + /// @brief Take the component-wise maximum of @c *this and @a other in place. + __hostdev__ constexpr Vec3& maxComponent(const Vec3& other) noexcept { Base::mergeMax(other); return *this; } + + /// @brief Return the smallest vector component + __hostdev__ [[nodiscard]] constexpr ValueType min() const noexcept { return Base::smallestComponent(); } + /// @brief Return the largest vector component + __hostdev__ [[nodiscard]] constexpr ValueType max() const noexcept { return Base::largestComponent(); } + + /// @brief Round each component down (toward negative infinity) + /// @return integer Coord + /// @note Only constexpr for integer @c T (floorAs uses non-constexpr math::Floor for floating point). + __hostdev__ [[nodiscard]] constexpr Coord floor() const noexcept { return Base::template floorAs(); } + /// @brief Round each component up (toward positive infinity) + /// @return integer Coord + /// @note Only constexpr for integer @c T (ceilAs uses non-constexpr math::Ceil for floating point). + __hostdev__ [[nodiscard]] constexpr Coord ceil() const noexcept { return Base::template ceilAs(); } + /// @brief Round each component to its closest integer value + /// @return integer Coord + /// @note Only constexpr for integer @c T (roundAs uses non-constexpr math::Floor for floating point). + __hostdev__ [[nodiscard]] constexpr Coord round() const noexcept { return Base::template roundAs(); } + + // ---- 3D-specific ---- + + /// @brief cross product with another 3-vector + template + __hostdev__ [[nodiscard]] constexpr Vec3 cross(const Vec3T& v) const noexcept { + return Vec3(this->mVec[1] * v[2] - this->mVec[2] * v[1], + this->mVec[2] * v[0] - this->mVec[0] * v[2], + this->mVec[0] * v[1] - this->mVec[1] * v[0]); } - __hostdev__ Vec3& operator/=(const T& s) { return (*this) *= T(1) / s; } - __hostdev__ Vec3& normalize() { return (*this) /= this->length(); } - /// @brief Perform a component-wise minimum with the other Coord. - __hostdev__ Vec3& minComponent(const Vec3& other) - { - if (other[0] < mVec[0]) - mVec[0] = other[0]; - if (other[1] < mVec[1]) - mVec[1] = other[1]; - if (other[2] < mVec[2]) - mVec[2] = other[2]; - return *this; + + /// @brief Outer product of a 3x1 vector and a 1x3 vector, result is a 3x3 matrix + template + __hostdev__ [[nodiscard]] constexpr Mat3 outer(const Vec3T& v) const noexcept { + return Mat3(this->mVec[0] * v[0], this->mVec[0] * v[1], this->mVec[0] * v[2], + this->mVec[1] * v[0], this->mVec[1] * v[1], this->mVec[1] * v[2], + this->mVec[2] * v[0], this->mVec[2] * v[1], this->mVec[2] * v[2]); } - /// @brief Perform a component-wise maximum with the other Coord. - __hostdev__ Vec3& maxComponent(const Vec3& other) + // ---- scalar * Vec / scalar / Vec (hidden friends — found only via ADL on Vec3, + // never participate in unrelated namespace-scope overload sets) ---- + + /// @brief Scalar-on-the-left multiplication (hidden friend). + template + __hostdev__ [[nodiscard]] friend constexpr Vec3 operator*(T1 scalar, const Vec3& vec) noexcept { - if (other[0] > mVec[0]) - mVec[0] = other[0]; - if (other[1] > mVec[1]) - mVec[1] = other[1]; - if (other[2] > mVec[2]) - mVec[2] = other[2]; - return *this; + return Vec3(scalar * vec[0], scalar * vec[1], scalar * vec[2]); } - /// @brief Return the smallest vector component - __hostdev__ ValueType min() const + /// @brief Scalar-on-the-left division (hidden friend). + template + __hostdev__ [[nodiscard]] friend constexpr Vec3 operator/(T1 scalar, const Vec3& vec) noexcept { - return mVec[0] < mVec[1] ? (mVec[0] < mVec[2] ? mVec[0] : mVec[2]) : (mVec[1] < mVec[2] ? mVec[1] : mVec[2]); - } - /// @brief Return the largest vector component - __hostdev__ ValueType max() const - { - return mVec[0] > mVec[1] ? (mVec[0] > mVec[2] ? mVec[0] : mVec[2]) : (mVec[1] > mVec[2] ? mVec[1] : mVec[2]); - } - /// @brief Round each component if this Vec up to its integer value - /// @return Return an integer Coord - __hostdev__ Coord floor() const { return Coord(Floor(mVec[0]), Floor(mVec[1]), Floor(mVec[2])); } - /// @brief Round each component if this Vec down to its integer value - /// @return Return an integer Coord - __hostdev__ Coord ceil() const { return Coord(Ceil(mVec[0]), Ceil(mVec[1]), Ceil(mVec[2])); } - /// @brief Round each component if this Vec to its closest integer value - /// @return Return an integer Coord - __hostdev__ Coord round() const - { - if constexpr(util::is_same::value) { - return Coord(Floor(mVec[0] + 0.5f), Floor(mVec[1] + 0.5f), Floor(mVec[2] + 0.5f)); - } else if constexpr(util::is_same::value) { - return Coord(mVec[0], mVec[1], mVec[2]); - } else { - return Coord(Floor(mVec[0] + 0.5), Floor(mVec[1] + 0.5), Floor(mVec[2] + 0.5)); - } + return Vec3(scalar / vec[0], scalar / vec[1], scalar / vec[2]); } - - /// @brief return a non-const raw constant pointer to array of three vector components - __hostdev__ T* asPointer() { return mVec; } - /// @brief return a const raw constant pointer to array of three vector components - __hostdev__ const T* asPointer() const { return mVec; } }; // Vec3 -template -__hostdev__ inline Vec3 operator*(T1 scalar, const Vec3& vec) -{ - return Vec3(scalar * vec[0], scalar * vec[1], scalar * vec[2]); -} -template -__hostdev__ inline Vec3 operator/(T1 scalar, const Vec3& vec) -{ - return Vec3(scalar / vec[0], scalar / vec[1], scalar / vec[2]); -} - /// @brief Return a single precision floating-point vector of this coordinate -__hostdev__ inline Vec3 Coord::asVec3s() const +__hostdev__ [[nodiscard]] inline constexpr Vec3 Coord::asVec3s() const noexcept { return Vec3(float(mVec[0]), float(mVec[1]), float(mVec[2])); } /// @brief Return a double precision floating-point vector of this coordinate -__hostdev__ inline Vec3 Coord::asVec3d() const +__hostdev__ [[nodiscard]] inline constexpr Vec3 Coord::asVec3d() const noexcept { return Vec3(double(mVec[0]), double(mVec[1]), double(mVec[2])); } @@ -1555,144 +1998,136 @@ __hostdev__ inline Vec3 Coord::asVec3d() const // ----------------------------> Vec4 <-------------------------------------- /// @brief A simple vector class with four components, similar to openvdb::math::Vec4 +/// +/// Aligned to 4*alignof(T) so the whole class fits in one SIMD register +/// (16 bytes for Vec4, 32 bytes for Vec4), without any +/// tail padding because the byte size is already a power-of-2 multiple +/// of alignof(T). template -class Vec4 +class alignas(alignof(T) * 4) Vec4 final : public VecBase { - T mVec[4]; + using Base = VecBase; public: - static const int SIZE = 4; - static const int size = 4; using ValueType = T; - Vec4() = default; - __hostdev__ explicit Vec4(T x) - : mVec{x, x, x, x} - { - } - __hostdev__ Vec4(T x, T y, T z, T w) - : mVec{x, y, z, w} - { - } + static constexpr int size = 4; // openvdb::math::Tuple-compat alias of SIZE + + /// @brief Default-construct (components are left uninitialized for fundamental @c T). + Vec4() noexcept = default; + /// @brief Broadcast: set all four components to @a x. + __hostdev__ explicit constexpr Vec4(T x) noexcept : Base(x, x, x, x) {} + /// @brief Component-wise construction. + __hostdev__ constexpr Vec4(T x, T y, T z, T w) noexcept : Base(x, y, z, w) {} + + /// @brief Explicit cross-precision conversion within nanovdb (e.g. @c Vec4d → @c Vec4f). template - __hostdev__ explicit Vec4(const Vec4& v) - : mVec{T(v[0]), T(v[1]), T(v[2]), T(v[3])} - { - } + __hostdev__ explicit constexpr Vec4(const Vec4& v) noexcept : Base(v[0], v[1], v[2], v[3]) {} + /// @brief Cross-template converting ctor (e.g. from @c openvdb::Vec4). Implicit + /// to preserve foreign-type interop. template class Vec4T, class T2> - __hostdev__ Vec4(const Vec4T& v) - : mVec{T(v[0]), T(v[1]), T(v[2]), T(v[3])} + __hostdev__ explicit constexpr Vec4(const Vec4T& v) noexcept : Base(v[0], v[1], v[2], v[3]) { - static_assert(Vec4T::size == size, "expected Vec4T::size==4!"); + static_assert(Vec4T::size == 4, "expected Vec4T::size==4!"); } - __hostdev__ bool operator==(const Vec4& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2] && mVec[3] == rhs[3]; } - __hostdev__ bool operator!=(const Vec4& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2] || mVec[3] != rhs[3]; } + /// @brief Assign from any 4-component vector type (foreign or nanovdb). template class Vec4T, class T2> - __hostdev__ Vec4& operator=(const Vec4T& rhs) - { - static_assert(Vec4T::size == size, "expected Vec4T::size==4!"); - mVec[0] = rhs[0]; - mVec[1] = rhs[1]; - mVec[2] = rhs[2]; - mVec[3] = rhs[3]; + __hostdev__ constexpr Vec4& operator=(const Vec4T& rhs) noexcept { + static_assert(Vec4T::size == 4, "expected Vec4T::size==4!"); + this->mVec[0] = rhs[0]; this->mVec[1] = rhs[1]; this->mVec[2] = rhs[2]; this->mVec[3] = rhs[3]; return *this; } - __hostdev__ const T& operator[](int i) const { return mVec[i]; } - __hostdev__ T& operator[](int i) { return mVec[i]; } - template - __hostdev__ T dot(const Vec4T& v) const { return mVec[0] * v[0] + mVec[1] * v[1] + mVec[2] * v[2] + mVec[3] * v[3]; } - __hostdev__ T lengthSqr() const - { - return mVec[0] * mVec[0] + mVec[1] * mVec[1] + mVec[2] * mVec[2] + mVec[3] * mVec[3]; // 7 flops - } - __hostdev__ T length() const { return Sqrt(this->lengthSqr()); } - __hostdev__ Vec4 operator-() const { return Vec4(-mVec[0], -mVec[1], -mVec[2], -mVec[3]); } - __hostdev__ Vec4 operator*(const Vec4& v) const { return Vec4(mVec[0] * v[0], mVec[1] * v[1], mVec[2] * v[2], mVec[3] * v[3]); } - __hostdev__ Vec4 operator/(const Vec4& v) const { return Vec4(mVec[0] / v[0], mVec[1] / v[1], mVec[2] / v[2], mVec[3] / v[3]); } - __hostdev__ Vec4 operator+(const Vec4& v) const { return Vec4(mVec[0] + v[0], mVec[1] + v[1], mVec[2] + v[2], mVec[3] + v[3]); } - __hostdev__ Vec4 operator-(const Vec4& v) const { return Vec4(mVec[0] - v[0], mVec[1] - v[1], mVec[2] - v[2], mVec[3] - v[3]); } - __hostdev__ Vec4 operator*(const T& s) const { return Vec4(s * mVec[0], s * mVec[1], s * mVec[2], s * mVec[3]); } - __hostdev__ Vec4 operator/(const T& s) const { return (T(1) / s) * (*this); } - __hostdev__ Vec4& operator+=(const Vec4& v) - { - mVec[0] += v[0]; - mVec[1] += v[1]; - mVec[2] += v[2]; - mVec[3] += v[3]; - return *this; - } - __hostdev__ Vec4& operator-=(const Vec4& v) - { - mVec[0] -= v[0]; - mVec[1] -= v[1]; - mVec[2] -= v[2]; - mVec[3] -= v[3]; - return *this; - } - __hostdev__ Vec4& operator*=(const T& s) + // ---- element-wise (Vec & Vec) ---- + /// @brief Component-wise negation. + __hostdev__ [[nodiscard]] constexpr Vec4 operator-() const noexcept { return Base::template negate(); } + /// @brief Component-wise sum. + __hostdev__ [[nodiscard]] constexpr Vec4 operator+(const Vec4& v) const noexcept { return Base::template plus(v); } + /// @brief Component-wise difference. + __hostdev__ [[nodiscard]] constexpr Vec4 operator-(const Vec4& v) const noexcept { return Base::template minus(v); } + /// @brief Component-wise (Hadamard) product. + __hostdev__ [[nodiscard]] constexpr Vec4 operator*(const Vec4& v) const noexcept { return Base::template mul(v); } + /// @brief Component-wise division. + __hostdev__ [[nodiscard]] constexpr Vec4 operator/(const Vec4& v) const noexcept { return Base::template div(v); } + /// @brief In-place component-wise addition. + __hostdev__ constexpr Vec4& operator+=(const Vec4& v) noexcept { Base::addAssign(v); return *this; } + /// @brief In-place component-wise subtraction. + __hostdev__ constexpr Vec4& operator-=(const Vec4& v) noexcept { Base::subAssign(v); return *this; } + + // ---- scalar ---- + /// @brief Component-wise multiply by scalar @a s. + __hostdev__ [[nodiscard]] constexpr Vec4 operator*(const T& s) const noexcept { return Base::template scale(s); } + /// @brief Component-wise divide by scalar @a s (integer-safe). + __hostdev__ [[nodiscard]] constexpr Vec4 operator/(const T& s) const noexcept { return Base::template divideBy(s); } + /// @brief In-place component-wise multiply by scalar @a s. + __hostdev__ constexpr Vec4& operator*=(const T& s) noexcept { Base::scaleAssign(s); return *this; } + /// @brief In-place component-wise divide by scalar @a s. + __hostdev__ constexpr Vec4& operator/=(const T& s) noexcept { Base::divideAssignScalar(s); return *this; } + /// @brief Normalize in place (divide by @c length()). Not @c constexpr — calls @c std::sqrt. + __hostdev__ Vec4& normalize() noexcept { return (*this) /= this->length(); } + /// @brief Return a normalized (unit-length) copy; const counterpart of @c normalize(). + __hostdev__ [[nodiscard]] Vec4 normalized() const noexcept { + return Base::template normalized(); + } + + // ---- equality ---- + /// @brief Component-wise equality. + __hostdev__ [[nodiscard]] constexpr bool operator==(const Vec4& rhs) const noexcept { return Base::equals(rhs); } + /// @brief Component-wise inequality. + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Vec4& rhs) const noexcept { return !Base::equals(rhs); } + + // ---- component-wise min/max ---- + /// @brief Take the component-wise minimum of @c *this and @a other in place. + __hostdev__ constexpr Vec4& minComponent(const Vec4& other) noexcept { Base::mergeMin(other); return *this; } + /// @brief Take the component-wise maximum of @c *this and @a other in place. + __hostdev__ constexpr Vec4& maxComponent(const Vec4& other) noexcept { Base::mergeMax(other); return *this; } + + /// @brief Return the smallest vector component + __hostdev__ [[nodiscard]] constexpr ValueType min() const noexcept { return Base::smallestComponent(); } + /// @brief Return the largest vector component + __hostdev__ [[nodiscard]] constexpr ValueType max() const noexcept { return Base::largestComponent(); } + + /// @brief Round each component down (toward negative infinity) + /// @return Vec4 (NanoVDB has no Coord4) + /// @note Only constexpr for integer @c T (floorAs uses non-constexpr math::Floor for floating point). + __hostdev__ [[nodiscard]] constexpr Vec4 floor() const noexcept { return Base::template floorAs>(); } + /// @brief Round each component up (toward positive infinity) + /// @return Vec4 + /// @note Only constexpr for integer @c T (ceilAs uses non-constexpr math::Ceil for floating point). + __hostdev__ [[nodiscard]] constexpr Vec4 ceil() const noexcept { return Base::template ceilAs>(); } + /// @brief Round each component to its closest integer value + /// @return Vec4 + /// @note Only constexpr for integer @c T (roundAs uses non-constexpr math::Floor for floating point). + __hostdev__ [[nodiscard]] constexpr Vec4 round() const noexcept { return Base::template roundAs>(); } + + // ---- scalar * Vec / scalar / Vec (hidden friends — found only via ADL on Vec4, + // never participate in unrelated namespace-scope overload sets) ---- + + /// @brief Scalar-on-the-left multiplication (hidden friend). + template + __hostdev__ [[nodiscard]] friend constexpr Vec4 operator*(T1 scalar, const Vec4& vec) noexcept { - mVec[0] *= s; - mVec[1] *= s; - mVec[2] *= s; - mVec[3] *= s; - return *this; + return Vec4(scalar * vec[0], scalar * vec[1], scalar * vec[2], scalar * vec[3]); } - __hostdev__ Vec4& operator/=(const T& s) { return (*this) *= T(1) / s; } - __hostdev__ Vec4& normalize() { return (*this) /= this->length(); } - /// @brief Perform a component-wise minimum with the other Coord. - __hostdev__ Vec4& minComponent(const Vec4& other) + /// @brief Scalar-on-the-left division (hidden friend). + template + __hostdev__ [[nodiscard]] friend constexpr Vec4 operator/(T1 scalar, const Vec4& vec) noexcept { - if (other[0] < mVec[0]) - mVec[0] = other[0]; - if (other[1] < mVec[1]) - mVec[1] = other[1]; - if (other[2] < mVec[2]) - mVec[2] = other[2]; - if (other[3] < mVec[3]) - mVec[3] = other[3]; - return *this; - } - - /// @brief Perform a component-wise maximum with the other Coord. - __hostdev__ Vec4& maxComponent(const Vec4& other) - { - if (other[0] > mVec[0]) - mVec[0] = other[0]; - if (other[1] > mVec[1]) - mVec[1] = other[1]; - if (other[2] > mVec[2]) - mVec[2] = other[2]; - if (other[3] > mVec[3]) - mVec[3] = other[3]; - return *this; + return Vec4(scalar / vec[0], scalar / vec[1], scalar / vec[2], scalar / vec[3]); } }; // Vec4 - -template -__hostdev__ inline Vec4 operator*(T1 scalar, const Vec4& vec) -{ - return Vec4(scalar * vec[0], scalar * vec[1], scalar * vec[2], scalar * vec[3]); -} -template -__hostdev__ inline Vec4 operator/(T1 scalar, const Vec4& vec) -{ - return Vec4(scalar / vec[0], scalar / vec[1], scalar / vec[2], scalar / vec[3]); -} -/// @brief Return the matrix vector product of a 4x4 matrix and a 4d vector -/// @param m 4x4 matrix -/// @param v 4d vector -/// @return result of matrix-vector multiplication, i.e. m x v -template -__hostdev__ inline Vec4 operator*(const Mat4& m, const Vec4& v) { - return Vec4( - m[0][0] * v[0] + m[0][1] * v[1] + m[0][2] * v[2] + m[0][3] * v[3], - m[1][0] * v[0] + m[1][1] * v[1] + m[1][2] * v[2] + m[1][3] * v[3], - m[2][0] * v[0] + m[2][1] * v[1] + m[2][2] * v[2] + m[2][3] * v[3], - m[3][0] * v[0] + m[3][1] * v[1] + m[3][2] * v[2] + m[3][3] * v[3] - ); -} - // ----------------------------> matMult <-------------------------------------- +// +// All six matMult / matMultT overloads were originally written with +// fma / fmaf for the single-rounding precision benefit. Those stdlib +// functions are not constexpr in C++17, which transitively blocked +// Map::applyMap and BBox::transform from being +// constexpr. Switching to plain `a * b + c` form gives back the +// constexpr-eligibility (worth ~1 ulp of rounding accuracy in the +// worst case, well below NanoVDB's geometric precision) and the +// device-side codegen is unchanged in practice — nvcc contracts +// `a * b + c` back into a hardware FMA by default (-fmad=true). +// +// (C++23's constexpr fma will eventually obviate this trade-off.) /// @brief Multiply a 3x3 matrix and a 3d vector using 32bit floating point arithmetics /// @note This corresponds to a linear mapping, e.g. scaling, rotation etc. @@ -1701,11 +2136,14 @@ __hostdev__ inline Vec4 operator*(const Mat4& m, const Vec4& v) { /// @param xyz input vector to be multiplied by the matrix /// @return result of matrix-vector multiplication, i.e. mat x xyz template -__hostdev__ inline Vec3T matMult(const float* mat, const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline constexpr Vec3T matMult(const float* mat, const Vec3T& xyz) noexcept { - return Vec3T(fmaf(static_cast(xyz[0]), mat[0], fmaf(static_cast(xyz[1]), mat[1], static_cast(xyz[2]) * mat[2])), - fmaf(static_cast(xyz[0]), mat[3], fmaf(static_cast(xyz[1]), mat[4], static_cast(xyz[2]) * mat[5])), - fmaf(static_cast(xyz[0]), mat[6], fmaf(static_cast(xyz[1]), mat[7], static_cast(xyz[2]) * mat[8]))); // 6 fmaf + 3 mult = 9 flops + const float x = static_cast(xyz[0]); + const float y = static_cast(xyz[1]); + const float z = static_cast(xyz[2]); + return Vec3T(x * mat[0] + y * mat[1] + z * mat[2], + x * mat[3] + y * mat[4] + z * mat[5], + x * mat[6] + y * mat[7] + z * mat[8]); } /// @brief Multiply a 3x3 matrix and a 3d vector using 64bit floating point arithmetics @@ -1715,11 +2153,14 @@ __hostdev__ inline Vec3T matMult(const float* mat, const Vec3T& xyz) /// @param xyz input vector to be multiplied by the matrix /// @return result of matrix-vector multiplication, i.e. mat x xyz template -__hostdev__ inline Vec3T matMult(const double* mat, const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline constexpr Vec3T matMult(const double* mat, const Vec3T& xyz) noexcept { - return Vec3T(fma(static_cast(xyz[0]), mat[0], fma(static_cast(xyz[1]), mat[1], static_cast(xyz[2]) * mat[2])), - fma(static_cast(xyz[0]), mat[3], fma(static_cast(xyz[1]), mat[4], static_cast(xyz[2]) * mat[5])), - fma(static_cast(xyz[0]), mat[6], fma(static_cast(xyz[1]), mat[7], static_cast(xyz[2]) * mat[8]))); // 6 fmaf + 3 mult = 9 flops + const double x = static_cast(xyz[0]); + const double y = static_cast(xyz[1]); + const double z = static_cast(xyz[2]); + return Vec3T(x * mat[0] + y * mat[1] + z * mat[2], + x * mat[3] + y * mat[4] + z * mat[5], + x * mat[6] + y * mat[7] + z * mat[8]); } /// @brief Multiply a 3x3 matrix to a 3d vector and add another 3d vector using 32bit floating point arithmetics @@ -1730,11 +2171,14 @@ __hostdev__ inline Vec3T matMult(const double* mat, const Vec3T& xyz) /// @param xyz input vector to be multiplied by the matrix and a translated by @c vec /// @return result of affine transformation, i.e. (mat x xyz) + vec template -__hostdev__ inline Vec3T matMult(const float* mat, const float* vec, const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline constexpr Vec3T matMult(const float* mat, const float* vec, const Vec3T& xyz) noexcept { - return Vec3T(fmaf(static_cast(xyz[0]), mat[0], fmaf(static_cast(xyz[1]), mat[1], fmaf(static_cast(xyz[2]), mat[2], vec[0]))), - fmaf(static_cast(xyz[0]), mat[3], fmaf(static_cast(xyz[1]), mat[4], fmaf(static_cast(xyz[2]), mat[5], vec[1]))), - fmaf(static_cast(xyz[0]), mat[6], fmaf(static_cast(xyz[1]), mat[7], fmaf(static_cast(xyz[2]), mat[8], vec[2])))); // 9 fmaf = 9 flops + const float x = static_cast(xyz[0]); + const float y = static_cast(xyz[1]); + const float z = static_cast(xyz[2]); + return Vec3T(x * mat[0] + y * mat[1] + z * mat[2] + vec[0], + x * mat[3] + y * mat[4] + z * mat[5] + vec[1], + x * mat[6] + y * mat[7] + z * mat[8] + vec[2]); } /// @brief Multiply a 3x3 matrix to a 3d vector and add another 3d vector using 64bit floating point arithmetics @@ -1745,11 +2189,14 @@ __hostdev__ inline Vec3T matMult(const float* mat, const float* vec, const Vec3T /// @param xyz input vector to be multiplied by the matrix and a translated by @c vec /// @return result of affine transformation, i.e. (mat x xyz) + vec template -__hostdev__ inline Vec3T matMult(const double* mat, const double* vec, const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline constexpr Vec3T matMult(const double* mat, const double* vec, const Vec3T& xyz) noexcept { - return Vec3T(fma(static_cast(xyz[0]), mat[0], fma(static_cast(xyz[1]), mat[1], fma(static_cast(xyz[2]), mat[2], vec[0]))), - fma(static_cast(xyz[0]), mat[3], fma(static_cast(xyz[1]), mat[4], fma(static_cast(xyz[2]), mat[5], vec[1]))), - fma(static_cast(xyz[0]), mat[6], fma(static_cast(xyz[1]), mat[7], fma(static_cast(xyz[2]), mat[8], vec[2])))); // 9 fma = 9 flops + const double x = static_cast(xyz[0]); + const double y = static_cast(xyz[1]); + const double z = static_cast(xyz[2]); + return Vec3T(x * mat[0] + y * mat[1] + z * mat[2] + vec[0], + x * mat[3] + y * mat[4] + z * mat[5] + vec[1], + x * mat[6] + y * mat[7] + z * mat[8] + vec[2]); } /// @brief Multiply the transposed of a 3x3 matrix and a 3d vector using 32bit floating point arithmetics @@ -1759,11 +2206,14 @@ __hostdev__ inline Vec3T matMult(const double* mat, const double* vec, const Vec /// @param xyz input vector to be multiplied by the transposed matrix /// @return result of matrix-vector multiplication, i.e. mat^T x xyz template -__hostdev__ inline Vec3T matMultT(const float* mat, const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline constexpr Vec3T matMultT(const float* mat, const Vec3T& xyz) noexcept { - return Vec3T(fmaf(static_cast(xyz[0]), mat[0], fmaf(static_cast(xyz[1]), mat[3], static_cast(xyz[2]) * mat[6])), - fmaf(static_cast(xyz[0]), mat[1], fmaf(static_cast(xyz[1]), mat[4], static_cast(xyz[2]) * mat[7])), - fmaf(static_cast(xyz[0]), mat[2], fmaf(static_cast(xyz[1]), mat[5], static_cast(xyz[2]) * mat[8]))); // 6 fmaf + 3 mult = 9 flops + const float x = static_cast(xyz[0]); + const float y = static_cast(xyz[1]); + const float z = static_cast(xyz[2]); + return Vec3T(x * mat[0] + y * mat[3] + z * mat[6], + x * mat[1] + y * mat[4] + z * mat[7], + x * mat[2] + y * mat[5] + z * mat[8]); } /// @brief Multiply the transposed of a 3x3 matrix and a 3d vector using 64bit floating point arithmetics @@ -1773,52 +2223,81 @@ __hostdev__ inline Vec3T matMultT(const float* mat, const Vec3T& xyz) /// @param xyz input vector to be multiplied by the transposed matrix /// @return result of matrix-vector multiplication, i.e. mat^T x xyz template -__hostdev__ inline Vec3T matMultT(const double* mat, const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline constexpr Vec3T matMultT(const double* mat, const Vec3T& xyz) noexcept { - return Vec3T(fma(static_cast(xyz[0]), mat[0], fma(static_cast(xyz[1]), mat[3], static_cast(xyz[2]) * mat[6])), - fma(static_cast(xyz[0]), mat[1], fma(static_cast(xyz[1]), mat[4], static_cast(xyz[2]) * mat[7])), - fma(static_cast(xyz[0]), mat[2], fma(static_cast(xyz[1]), mat[5], static_cast(xyz[2]) * mat[8]))); // 6 fmaf + 3 mult = 9 flops + const double x = static_cast(xyz[0]); + const double y = static_cast(xyz[1]); + const double z = static_cast(xyz[2]); + return Vec3T(x * mat[0] + y * mat[3] + z * mat[6], + x * mat[1] + y * mat[4] + z * mat[7], + x * mat[2] + y * mat[5] + z * mat[8]); } +/// @brief Multiply the transpose of a 3x3 matrix by @a xyz and add @a vec (32-bit floats). +/// @note Corresponds to an inverse affine transform: @c (mat^T x xyz) + @c vec. +/// @tparam Vec3T Template type of the input and output 3d vectors template -__hostdev__ inline Vec3T matMultT(const float* mat, const float* vec, const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline constexpr Vec3T matMultT(const float* mat, const float* vec, const Vec3T& xyz) noexcept { - return Vec3T(fmaf(static_cast(xyz[0]), mat[0], fmaf(static_cast(xyz[1]), mat[3], fmaf(static_cast(xyz[2]), mat[6], vec[0]))), - fmaf(static_cast(xyz[0]), mat[1], fmaf(static_cast(xyz[1]), mat[4], fmaf(static_cast(xyz[2]), mat[7], vec[1]))), - fmaf(static_cast(xyz[0]), mat[2], fmaf(static_cast(xyz[1]), mat[5], fmaf(static_cast(xyz[2]), mat[8], vec[2])))); // 9 fmaf = 9 flops + const float x = static_cast(xyz[0]); + const float y = static_cast(xyz[1]); + const float z = static_cast(xyz[2]); + return Vec3T(x * mat[0] + y * mat[3] + z * mat[6] + vec[0], + x * mat[1] + y * mat[4] + z * mat[7] + vec[1], + x * mat[2] + y * mat[5] + z * mat[8] + vec[2]); } +/// @brief Multiply the transpose of a 3x3 matrix by @a xyz and add @a vec (64-bit floats). +/// @note Corresponds to an inverse affine transform: @c (mat^T x xyz) + @c vec. +/// @tparam Vec3T Template type of the input and output 3d vectors template -__hostdev__ inline Vec3T matMultT(const double* mat, const double* vec, const Vec3T& xyz) +__hostdev__ [[nodiscard]] inline constexpr Vec3T matMultT(const double* mat, const double* vec, const Vec3T& xyz) noexcept { - return Vec3T(fma(static_cast(xyz[0]), mat[0], fma(static_cast(xyz[1]), mat[3], fma(static_cast(xyz[2]), mat[6], vec[0]))), - fma(static_cast(xyz[0]), mat[1], fma(static_cast(xyz[1]), mat[4], fma(static_cast(xyz[2]), mat[7], vec[1]))), - fma(static_cast(xyz[0]), mat[2], fma(static_cast(xyz[1]), mat[5], fma(static_cast(xyz[2]), mat[8], vec[2])))); // 9 fma = 9 flops + const double x = static_cast(xyz[0]); + const double y = static_cast(xyz[1]); + const double z = static_cast(xyz[2]); + return Vec3T(x * mat[0] + y * mat[3] + z * mat[6] + vec[0], + x * mat[1] + y * mat[4] + z * mat[7] + vec[1], + x * mat[2] + y * mat[5] + z * mat[8] + vec[2]); } // ----------------------------> BBox <------------------------------------- -// Base-class for static polymorphism (cannot be constructed directly) +/// @brief Common base for floating-point and integer bounding boxes; not +/// constructible directly (only as a base of @c BBox). +/// @details Stores the closed-segment endpoints @c mCoord[0] (min corner) and +/// @c mCoord[1] (max corner). Whether the box is "min-inclusive max-inclusive" +/// or "min-inclusive max-exclusive" depends on the derived @c BBox +/// specialization (integer vs floating point). template struct BaseBBox { Vec3T mCoord[2]; - __hostdev__ bool operator==(const BaseBBox& rhs) const { return mCoord[0] == rhs.mCoord[0] && mCoord[1] == rhs.mCoord[1]; }; - __hostdev__ bool operator!=(const BaseBBox& rhs) const { return mCoord[0] != rhs.mCoord[0] || mCoord[1] != rhs.mCoord[1]; }; - __hostdev__ const Vec3T& operator[](int i) const { return mCoord[i]; } - __hostdev__ Vec3T& operator[](int i) { return mCoord[i]; } - __hostdev__ Vec3T& min() { return mCoord[0]; } - __hostdev__ Vec3T& max() { return mCoord[1]; } - __hostdev__ const Vec3T& min() const { return mCoord[0]; } - __hostdev__ const Vec3T& max() const { return mCoord[1]; } - __hostdev__ BaseBBox& translate(const Vec3T& xyz) + /// @brief Equality on the two corner points. + __hostdev__ [[nodiscard]] constexpr bool operator==(const BaseBBox& rhs) const noexcept { return mCoord[0] == rhs.mCoord[0] && mCoord[1] == rhs.mCoord[1]; }; + /// @brief Inequality on the two corner points. + __hostdev__ [[nodiscard]] constexpr bool operator!=(const BaseBBox& rhs) const noexcept { return mCoord[0] != rhs.mCoord[0] || mCoord[1] != rhs.mCoord[1]; }; + /// @brief Indexed corner access: @a i = 0 returns min, @a i = 1 returns max. + __hostdev__ constexpr const Vec3T& operator[](int i) const noexcept { NANOVDB_ASSERT(i >= 0 && i < 2); return mCoord[i]; } + /// @brief Mutable variant of @c operator[]. + __hostdev__ constexpr Vec3T& operator[](int i) noexcept { NANOVDB_ASSERT(i >= 0 && i < 2); return mCoord[i]; } + /// @brief Mutable accessor for the min corner. + __hostdev__ constexpr Vec3T& min() noexcept { return mCoord[0]; } + /// @brief Mutable accessor for the max corner. + __hostdev__ constexpr Vec3T& max() noexcept { return mCoord[1]; } + /// @brief Const accessor for the min corner. + __hostdev__ constexpr const Vec3T& min() const noexcept { return mCoord[0]; } + /// @brief Const accessor for the max corner. + __hostdev__ constexpr const Vec3T& max() const noexcept { return mCoord[1]; } + /// @brief Translate (rigid-shift) both corners by @a xyz; return @c *this. + __hostdev__ constexpr BaseBBox& translate(const Vec3T& xyz) noexcept { mCoord[0] += xyz; mCoord[1] += xyz; return *this; } /// @brief Expand this bounding box to enclose point @c xyz. - __hostdev__ BaseBBox& expand(const Vec3T& xyz) + __hostdev__ constexpr BaseBBox& expand(const Vec3T& xyz) noexcept { mCoord[0].minComponent(xyz); mCoord[1].maxComponent(xyz); @@ -1826,7 +2305,7 @@ struct BaseBBox } /// @brief Expand this bounding box to enclose the given bounding box. - __hostdev__ BaseBBox& expand(const BaseBBox& bbox) + __hostdev__ constexpr BaseBBox& expand(const BaseBBox& bbox) noexcept { mCoord[0].minComponent(bbox[0]); mCoord[1].maxComponent(bbox[1]); @@ -1834,18 +2313,16 @@ struct BaseBBox } /// @brief Intersect this bounding box with the given bounding box. - __hostdev__ BaseBBox& intersect(const BaseBBox& bbox) + __hostdev__ constexpr BaseBBox& intersect(const BaseBBox& bbox) noexcept { mCoord[0].maxComponent(bbox[0]); mCoord[1].minComponent(bbox[1]); return *this; } -// __hostdev__ BaseBBox expandBy(typename Vec3T::ValueType padding) const -// { -// return BaseBBox(mCoord[0].offsetBy(-padding),mCoord[1].offsetBy(padding)); -// } - __hostdev__ bool isInside(const Vec3T& xyz) + /// @brief Return @c true iff @a xyz lies in the closed box (each component + /// satisfies @c min[i] <= @c xyz[i] <= @c max[i]). + __hostdev__ [[nodiscard]] constexpr bool isInside(const Vec3T& xyz) const noexcept { if (xyz[0] < mCoord[0][0] || xyz[1] < mCoord[0][1] || xyz[2] < mCoord[0][2]) return false; @@ -1855,8 +2332,11 @@ struct BaseBBox } protected: - __hostdev__ BaseBBox() {} - __hostdev__ BaseBBox(const Vec3T& min, const Vec3T& max) + /// @brief Default-construct (leaves corners uninitialized — derived classes + /// supply a meaningful default). + __hostdev__ constexpr BaseBBox() noexcept {} + /// @brief Construct from @a min and @a max corners directly. + __hostdev__ constexpr BaseBBox(const Vec3T& min, const Vec3T& max) noexcept : mCoord{min, max} { } @@ -1878,37 +2358,47 @@ struct BBox : public BaseBBox using BaseT = BaseBBox; using BaseT::mCoord; /// @brief Default construction sets BBox to an empty bbox - __hostdev__ BBox() + __hostdev__ constexpr BBox() noexcept : BaseT(Vec3T( Maximum::value()), Vec3T(-Maximum::value())) { } - __hostdev__ BBox(const Vec3T& min, const Vec3T& max) + /// @brief Construct from explicit @a min / @a max corners. + __hostdev__ constexpr BBox(const Vec3T& min, const Vec3T& max) noexcept : BaseT(min, max) { } - __hostdev__ BBox(const Coord& min, const Coord& max) + /// @brief Convert an integer @c Coord box into this floating-point box. + /// @note @c max is exclusive in the floating-point convention, so the + /// integer @c max corner is incremented by 1 before conversion. + __hostdev__ constexpr BBox(const Coord& min, const Coord& max) noexcept : BaseT(Vec3T(ValueType(min[0]), ValueType(min[1]), ValueType(min[2])), Vec3T(ValueType(max[0] + 1), ValueType(max[1] + 1), ValueType(max[2] + 1))) { } - __hostdev__ static BBox createCube(const Coord& min, typename Coord::ValueType dim) + /// @brief Return a cube-shaped @c BBox with min corner @a min and edge length @a dim. + __hostdev__ [[nodiscard]] static constexpr BBox createCube(const Coord& min, typename Coord::ValueType dim) noexcept { return BBox(min, min.offsetBy(dim)); } - __hostdev__ BBox(const BaseBBox& bbox) + /// @brief Construct from a @c BaseBBox; delegates to the @c (min, @c max) ctor. + __hostdev__ constexpr BBox(const BaseBBox& bbox) noexcept : BBox(bbox[0], bbox[1]) { } - __hostdev__ bool empty() const { return mCoord[0][0] >= mCoord[1][0] || + /// @brief Return @c true if this bounding box is empty (max <= min in any axis). + __hostdev__ [[nodiscard]] constexpr bool empty() const noexcept { return mCoord[0][0] >= mCoord[1][0] || mCoord[0][1] >= mCoord[1][1] || mCoord[0][2] >= mCoord[1][2]; } - __hostdev__ operator bool() const { return mCoord[0][0] < mCoord[1][0] && + /// @brief Convert to bool: @c true iff this bounding box has positive volume. + __hostdev__ [[nodiscard]] constexpr operator bool() const noexcept { return mCoord[0][0] < mCoord[1][0] && mCoord[0][1] < mCoord[1][1] && mCoord[0][2] < mCoord[1][2]; } - __hostdev__ Vec3T dim() const { return *this ? this->max() - this->min() : Vec3T(0); } - __hostdev__ bool isInside(const Vec3T& p) const + /// @brief Return @c max - @c min when non-empty, else the zero vector. + __hostdev__ [[nodiscard]] constexpr Vec3T dim() const noexcept { return *this ? this->max() - this->min() : Vec3T(0); } + /// @brief Return @c true iff @a p is strictly inside the box (open interval — min is exclusive, max is exclusive). + __hostdev__ [[nodiscard]] constexpr bool isInside(const Vec3T& p) const noexcept { return p[0] > mCoord[0][0] && p[1] > mCoord[0][1] && p[2] > mCoord[0][2] && p[0] < mCoord[1][0] && p[1] < mCoord[1][1] && p[2] < mCoord[1][2]; @@ -1934,17 +2424,20 @@ struct BBox : public BaseBBox CoordT mPos; public: - __hostdev__ Iterator(const BBox& b) + /// @brief Construct an iterator positioned at @c b.min(). + __hostdev__ constexpr Iterator(const BBox& b) noexcept : mBBox(b) , mPos(b.min()) { } - __hostdev__ Iterator(const BBox& b, const Coord& p) + /// @brief Construct an iterator positioned at the given @a p inside @a b. + __hostdev__ constexpr Iterator(const BBox& b, const Coord& p) noexcept : mBBox(b) , mPos(p) { } - __hostdev__ Iterator& operator++() + /// @brief Pre-increment: advance to the next coordinate in row-major z-fastest order. + __hostdev__ constexpr Iterator& operator++() noexcept { if (mPos[2] < mBBox[1][2]) { // this is the most common case ++mPos[2];// increment z @@ -1958,49 +2451,64 @@ struct BBox : public BaseBBox } return *this; } - __hostdev__ Iterator operator++(int) + /// @brief Post-increment: advance to the next coordinate; return a copy of the previous state. + __hostdev__ constexpr Iterator operator++(int) noexcept { auto tmp = *this; ++(*this); return tmp; } - __hostdev__ bool operator==(const Iterator& rhs) const + /// @brief Iterator equality (asserts both iterators belong to the same @c BBox). + __hostdev__ [[nodiscard]] constexpr bool operator==(const Iterator& rhs) const noexcept { NANOVDB_ASSERT(mBBox == rhs.mBBox); return mPos == rhs.mPos; } - __hostdev__ bool operator!=(const Iterator& rhs) const + /// @brief Iterator inequality (asserts both iterators belong to the same @c BBox). + __hostdev__ [[nodiscard]] constexpr bool operator!=(const Iterator& rhs) const noexcept { NANOVDB_ASSERT(mBBox == rhs.mBBox); return mPos != rhs.mPos; } - __hostdev__ bool operator<(const Iterator& rhs) const + /// @brief Lexicographic less-than over the underlying position. + __hostdev__ [[nodiscard]] constexpr bool operator<(const Iterator& rhs) const noexcept { NANOVDB_ASSERT(mBBox == rhs.mBBox); return mPos < rhs.mPos; } - __hostdev__ bool operator<=(const Iterator& rhs) const + /// @brief Lexicographic less-than-or-equal over the underlying position. + __hostdev__ [[nodiscard]] constexpr bool operator<=(const Iterator& rhs) const noexcept { NANOVDB_ASSERT(mBBox == rhs.mBBox); return mPos <= rhs.mPos; } /// @brief Return @c true if the iterator still points to a valid coordinate. - __hostdev__ operator bool() const { return mPos <= mBBox[1]; } - __hostdev__ const CoordT& operator*() const { return mPos; } + __hostdev__ [[nodiscard]] constexpr operator bool() const noexcept { return mPos <= mBBox[1]; } + /// @brief Dereference to the current coordinate. + __hostdev__ [[nodiscard]] constexpr const CoordT& operator*() const noexcept { return mPos; } }; // Iterator - __hostdev__ Iterator begin() const { return Iterator{*this}; } - __hostdev__ Iterator end() const { return Iterator{*this, CoordT(mCoord[1][0]+1, mCoord[0][1], mCoord[0][2])}; } - __hostdev__ BBox() + /// @brief Iterator positioned at the min corner. + __hostdev__ [[nodiscard]] constexpr Iterator begin() const noexcept { return Iterator{*this}; } + /// @brief One-past-the-end iterator (max corner shifted by one in @c x). + __hostdev__ [[nodiscard]] constexpr Iterator end() const noexcept { return Iterator{*this, CoordT(mCoord[1][0]+1, mCoord[0][1], mCoord[0][2])}; } + /// @brief Default construct an empty bbox (@c min = @c CoordT::max(), @c max = @c CoordT::min()). + __hostdev__ constexpr BBox() noexcept : BaseT(CoordT::max(), CoordT::min()) { } - __hostdev__ BBox(const CoordT& min, const CoordT& max) + /// @brief Construct from explicit @a min / @a max integer corners. + __hostdev__ constexpr BBox(const CoordT& min, const CoordT& max) noexcept : BaseT(min, max) { } + /// @brief Splitting constructor (used by parallel range-based iteration). + /// @details Halves the box along its longest axis: @c other keeps the upper + /// half, the newly constructed @c BBox is the lower half. The @c SplitT tag + /// type is ignored; it exists so this overload doesn't collide with the + /// copy constructor. template - __hostdev__ BBox(BBox& other, const SplitT&) + __hostdev__ constexpr BBox(BBox& other, const SplitT&) noexcept : BaseT(other.mCoord[0], other.mCoord[1]) { NANOVDB_ASSERT(this->is_divisible()); @@ -2009,65 +2517,72 @@ struct BBox : public BaseBBox other.mCoord[0][n] = mCoord[1][n] + 1; } - __hostdev__ static BBox createCube(const CoordT& min, typename CoordT::ValueType dim) + /// @brief Return a cube-shaped @c BBox with min corner @a min and edge length @a dim. + /// @note Subtracts 1 from @a dim because the integer convention is min/max inclusive. + __hostdev__ [[nodiscard]] static constexpr BBox createCube(const CoordT& min, typename CoordT::ValueType dim) noexcept { return BBox(min, min.offsetBy(dim - 1)); } - __hostdev__ static BBox createCube(typename CoordT::ValueType min, typename CoordT::ValueType max) + /// @brief Return a cube-shaped @c BBox spanning @a min..@a max in every axis. + __hostdev__ [[nodiscard]] static constexpr BBox createCube(typename CoordT::ValueType min, typename CoordT::ValueType max) noexcept { return BBox(CoordT(min), CoordT(max)); } - __hostdev__ bool is_divisible() const { return mCoord[0][0] < mCoord[1][0] && + /// @brief Return @c true iff the box has more than one cell along every axis. + __hostdev__ [[nodiscard]] constexpr bool is_divisible() const noexcept { return mCoord[0][0] < mCoord[1][0] && mCoord[0][1] < mCoord[1][1] && mCoord[0][2] < mCoord[1][2]; } /// @brief Return true if this bounding box is empty, e.g. uninitialized - __hostdev__ bool empty() const { return mCoord[0][0] > mCoord[1][0] || + __hostdev__ [[nodiscard]] constexpr bool empty() const noexcept { return mCoord[0][0] > mCoord[1][0] || mCoord[0][1] > mCoord[1][1] || mCoord[0][2] > mCoord[1][2]; } /// @brief Convert this BBox to boolean true if it is not empty - __hostdev__ operator bool() const { return mCoord[0][0] <= mCoord[1][0] && + __hostdev__ [[nodiscard]] constexpr operator bool() const noexcept { return mCoord[0][0] <= mCoord[1][0] && mCoord[0][1] <= mCoord[1][1] && mCoord[0][2] <= mCoord[1][2]; } - __hostdev__ CoordT dim() const { return *this ? this->max() - this->min() + Coord(1) : Coord(0); } - __hostdev__ uint64_t volume() const + /// @brief Return the per-axis extent (max - min + 1 when non-empty, else 0). + __hostdev__ [[nodiscard]] constexpr CoordT dim() const noexcept { return *this ? this->max() - this->min() + Coord(1) : Coord(0); } + /// @brief Return the number of integer cells in the box, i.e. @c dim().x * @c dim().y * @c dim().z. + __hostdev__ [[nodiscard]] constexpr uint64_t volume() const noexcept { auto d = this->dim(); return uint64_t(d[0]) * uint64_t(d[1]) * uint64_t(d[2]); } - __hostdev__ bool isInside(const CoordT& p) const { return !(CoordT::lessThan(p, this->min()) || CoordT::lessThan(this->max(), p)); } + /// @brief Return @c true iff the integer coordinate @a p lies inside the closed box. + __hostdev__ [[nodiscard]] constexpr bool isInside(const CoordT& p) const noexcept { return !(CoordT::lessThan(p, this->min()) || CoordT::lessThan(this->max(), p)); } /// @brief Return @c true if the given bounding box is inside this bounding box. - __hostdev__ bool isInside(const BBox& b) const + __hostdev__ [[nodiscard]] constexpr bool isInside(const BBox& b) const noexcept { return !(CoordT::lessThan(b.min(), this->min()) || CoordT::lessThan(this->max(), b.max())); } /// @brief Return @c true if the given bounding box overlaps with this bounding box. - __hostdev__ bool hasOverlap(const BBox& b) const + __hostdev__ [[nodiscard]] constexpr bool hasOverlap(const BBox& b) const noexcept { return !(CoordT::lessThan(this->max(), b.min()) || CoordT::lessThan(b.max(), this->min())); } /// @warning This converts a CoordBBox into a floating-point bounding box which implies that max += 1 ! template - __hostdev__ BBox> asReal() const + __hostdev__ [[nodiscard]] constexpr BBox> asReal() const noexcept { static_assert(util::is_floating_point::value, "CoordBBox::asReal: Expected a floating point coordinate"); return BBox>(Vec3(RealT(mCoord[0][0]), RealT(mCoord[0][1]), RealT(mCoord[0][2])), Vec3(RealT(mCoord[1][0] + 1), RealT(mCoord[1][1] + 1), RealT(mCoord[1][2] + 1))); } /// @brief Return a new instance that is expanded by the specified padding. - __hostdev__ BBox expandBy(typename CoordT::ValueType padding) const + __hostdev__ [[nodiscard]] constexpr BBox expandBy(typename CoordT::ValueType padding) const noexcept { return BBox(mCoord[0].offsetBy(-padding), mCoord[1].offsetBy(padding)); } - /// @brief @brief transform this coordinate bounding box by the specified map + /// @brief Transform this coordinate bounding box by the specified map. /// @param map mapping of index to world coordinates /// @return world bounding box template - __hostdev__ auto transform(const Map& map) const + __hostdev__ [[nodiscard]] constexpr auto transform(const Map& map) const noexcept { using Vec3T = Vec3; const Vec3T tmp = map.applyMap(Vec3T(mCoord[0][0], mCoord[0][1], mCoord[0][2])); @@ -2083,19 +2598,22 @@ struct BBox : public BaseBBox } #if defined(__CUDACC__) // the following functions only run on the GPU! - __device__ inline BBox& expandAtomic(const CoordT& ijk) + /// @brief Device-only: atomically expand the box to enclose the integer point @a ijk. + __device__ inline BBox& expandAtomic(const CoordT& ijk) noexcept { mCoord[0].minComponentAtomic(ijk); mCoord[1].maxComponentAtomic(ijk); return *this; } - __device__ inline BBox& expandAtomic(const BBox& bbox) + /// @brief Device-only: atomically expand the box to enclose another @a bbox. + __device__ inline BBox& expandAtomic(const BBox& bbox) noexcept { mCoord[0].minComponentAtomic(bbox[0]); mCoord[1].maxComponentAtomic(bbox[1]); return *this; } - __device__ inline BBox& intersectAtomic(const BBox& bbox) + /// @brief Device-only: atomically intersect this box with @a bbox. + __device__ inline BBox& intersectAtomic(const BBox& bbox) noexcept { mCoord[0].maxComponentAtomic(bbox[0]); mCoord[1].minComponentAtomic(bbox[1]); @@ -2120,21 +2638,21 @@ class Rgba8 using ValueType = uint8_t; /// @brief Default copy constructor - Rgba8(const Rgba8&) = default; + Rgba8(const Rgba8&) noexcept = default; /// @brief Default move constructor - Rgba8(Rgba8&&) = default; + Rgba8(Rgba8&&) noexcept = default; /// @brief Default move assignment operator /// @return non-const reference to this instance - Rgba8& operator=(Rgba8&&) = default; + Rgba8& operator=(Rgba8&&) noexcept = default; /// @brief Default copy assignment operator /// @return non-const reference to this instance - Rgba8& operator=(const Rgba8&) = default; + Rgba8& operator=(const Rgba8&) noexcept = default; /// @brief Default ctor initializes all channels to zero - __hostdev__ Rgba8() + __hostdev__ constexpr Rgba8() noexcept : mData{{0, 0, 0, 0}} { static_assert(sizeof(uint32_t) == sizeof(Rgba8), "Unexpected sizeof"); @@ -2142,21 +2660,21 @@ class Rgba8 /// @brief integer r,g,b,a ctor where alpha channel defaults to opaque /// @note all values should be in the range 0u to 255u - __hostdev__ Rgba8(uint8_t r, uint8_t g, uint8_t b, uint8_t a = 255u) + __hostdev__ constexpr Rgba8(uint8_t r, uint8_t g, uint8_t b, uint8_t a = 255u) noexcept : mData{{r, g, b, a}} { } - /// @brief @brief ctor where all channels are initialized to the same value + /// @brief Ctor where all channels are initialized to the same value. /// @note value should be in the range 0u to 255u - explicit __hostdev__ Rgba8(uint8_t v) + explicit __hostdev__ constexpr Rgba8(uint8_t v) noexcept : mData{{v, v, v, v}} { } /// @brief floating-point r,g,b,a ctor where alpha channel defaults to opaque /// @note all values should be in the range 0.0f to 1.0f - __hostdev__ Rgba8(float r, float g, float b, float a = 1.0f) + __hostdev__ constexpr Rgba8(float r, float g, float b, float a = 1.0f) noexcept : mData{{static_cast(0.5f + r * 255.0f), // round floats to nearest integers static_cast(0.5f + g * 255.0f), // double {{}} is needed due to union static_cast(0.5f + b * 255.0f), @@ -2164,47 +2682,61 @@ class Rgba8 { } - /// @brief Vec3f r,g,b ctor (alpha channel it set to 1) + /// @brief Construct from a @c Vec3 with rgb components (alpha set to opaque). /// @note all values should be in the range 0.0f to 1.0f - __hostdev__ Rgba8(const Vec3& rgb) + __hostdev__ constexpr Rgba8(const Vec3& rgb) noexcept : Rgba8(rgb[0], rgb[1], rgb[2]) { } - /// @brief Vec4f r,g,b,a ctor + /// @brief Construct from a @c Vec4 with rgba components. /// @note all values should be in the range 0.0f to 1.0f - __hostdev__ Rgba8(const Vec4& rgba) + __hostdev__ constexpr Rgba8(const Vec4& rgba) noexcept : Rgba8(rgba[0], rgba[1], rgba[2], rgba[3]) { } - __hostdev__ bool operator< (const Rgba8& rhs) const { return mData.packed < rhs.mData.packed; } - __hostdev__ bool operator==(const Rgba8& rhs) const { return mData.packed == rhs.mData.packed; } - __hostdev__ float lengthSqr() const + /// @brief Compare on the 32-bit packed integer representation (lexicographic by byte). + __hostdev__ [[nodiscard]] bool operator< (const Rgba8& rhs) const noexcept { return mData.packed < rhs.mData.packed; } + /// @brief Equality on the 32-bit packed integer representation. + __hostdev__ [[nodiscard]] bool operator==(const Rgba8& rhs) const noexcept { return mData.packed == rhs.mData.packed; } + /// @brief Return the squared L2 length of the rgb channels in the [0, 1] range + /// (alpha is ignored). @c 0.0000153787005f is @c 1/255^2. + __hostdev__ [[nodiscard]] constexpr float lengthSqr() const noexcept { return 0.0000153787005f * (float(mData.c[0]) * mData.c[0] + float(mData.c[1]) * mData.c[1] + float(mData.c[2]) * mData.c[2]); //1/255^2 } - __hostdev__ float length() const { return sqrtf(this->lengthSqr()); } - /// @brief return n'th color channel as a float in the range 0 to 1 - __hostdev__ float asFloat(int n) const { return 0.003921569f*float(mData.c[n]); }// divide by 255 - __hostdev__ const uint8_t& operator[](int n) const { return mData.c[n]; } - __hostdev__ uint8_t& operator[](int n) { return mData.c[n]; } - __hostdev__ const uint32_t& packed() const { return mData.packed; } - __hostdev__ uint32_t& packed() { return mData.packed; } - __hostdev__ const uint8_t& r() const { return mData.c[0]; } - __hostdev__ const uint8_t& g() const { return mData.c[1]; } - __hostdev__ const uint8_t& b() const { return mData.c[2]; } - __hostdev__ const uint8_t& a() const { return mData.c[3]; } - __hostdev__ uint8_t& r() { return mData.c[0]; } - __hostdev__ uint8_t& g() { return mData.c[1]; } - __hostdev__ uint8_t& b() { return mData.c[2]; } - __hostdev__ uint8_t& a() { return mData.c[3]; } - __hostdev__ operator Vec3() const { + /// @brief L2 length of the rgb channels (alpha ignored). Not @c constexpr — calls @c sqrtf. + __hostdev__ [[nodiscard]] float length() const noexcept { return sqrtf(this->lengthSqr()); } + /// @brief Return the @a n'th color channel as a float in the range 0 to 1. + __hostdev__ [[nodiscard]] constexpr float asFloat(int n) const noexcept { return 0.003921569f*float(mData.c[n]); }// divide by 255 + /// @brief Indexed channel access. Asserts @c 0 <= @a n < 4 in debug builds. + __hostdev__ constexpr const uint8_t& operator[](int n) const noexcept { NANOVDB_ASSERT(n >= 0 && n < 4); return mData.c[n]; } + /// @brief Mutable variant of @c operator[]. + __hostdev__ constexpr uint8_t& operator[](int n) noexcept { NANOVDB_ASSERT(n >= 0 && n < 4); return mData.c[n]; } + /// @brief Const access to the 32-bit packed integer representation. + __hostdev__ const uint32_t& packed() const noexcept { return mData.packed; } + /// @brief Mutable access to the 32-bit packed integer representation. + __hostdev__ uint32_t& packed() noexcept { return mData.packed; } + //@{ + /// @brief Named channel accessors (r, g, b, a). Both const and non-const variants. + __hostdev__ constexpr const uint8_t& r() const noexcept { return mData.c[0]; } + __hostdev__ constexpr const uint8_t& g() const noexcept { return mData.c[1]; } + __hostdev__ constexpr const uint8_t& b() const noexcept { return mData.c[2]; } + __hostdev__ constexpr const uint8_t& a() const noexcept { return mData.c[3]; } + __hostdev__ constexpr uint8_t& r() noexcept { return mData.c[0]; } + __hostdev__ constexpr uint8_t& g() noexcept { return mData.c[1]; } + __hostdev__ constexpr uint8_t& b() noexcept { return mData.c[2]; } + __hostdev__ constexpr uint8_t& a() noexcept { return mData.c[3]; } + //@} + /// @brief Implicit conversion to a @c Vec3 with channels in [0, 1] (alpha dropped). + __hostdev__ [[nodiscard]] constexpr operator Vec3() const noexcept { return Vec3(this->asFloat(0), this->asFloat(1), this->asFloat(2)); } - __hostdev__ operator Vec4() const { + /// @brief Implicit conversion to a @c Vec4 with rgba channels in [0, 1]. + __hostdev__ [[nodiscard]] constexpr operator Vec4() const noexcept { return Vec4(this->asFloat(0), this->asFloat(1), this->asFloat(2), this->asFloat(3)); } }; // Rgba8 diff --git a/nanovdb/nanovdb/python/cuda/PySampleFromVoxels.cu b/nanovdb/nanovdb/python/cuda/PySampleFromVoxels.cu index aaebf66772..372ed2c87f 100644 --- a/nanovdb/nanovdb/python/cuda/PySampleFromVoxels.cu +++ b/nanovdb/nanovdb/python/cuda/PySampleFromVoxels.cu @@ -40,7 +40,9 @@ __global__ void sampleFromVoxels(unsigned int numPoints, const BuildT* points, c math::SampleFromVoxels sampler(d_grid->tree()); values[i] = sampler(indexPos); - Vec3T inv2Dx = (BuildT).5 / d_grid->voxelSize(); + // voxelSize() returns Vec3d; explicit conversion to Vec3T preserves the existing + // "compute reciprocal in double, then narrow to BuildT" semantics. + Vec3T inv2Dx = Vec3T((BuildT).5 / d_grid->voxelSize()); Vec3T gradient = Vec3T(sampler(indexPos + Vec3T(1, 0, 0)) - sampler(indexPos - Vec3T(1, 0, 0)), sampler(indexPos + Vec3T(0, 1, 0)) - sampler(indexPos - Vec3T(0, 1, 0)), sampler(indexPos + Vec3T(0, 0, 1)) - sampler(indexPos - Vec3T(0, 0, 1))) * diff --git a/nanovdb/nanovdb/tools/cuda/PointsToGrid.cuh b/nanovdb/nanovdb/tools/cuda/PointsToGrid.cuh index 5b22b32054..eaebead831 100644 --- a/nanovdb/nanovdb/tools/cuda/PointsToGrid.cuh +++ b/nanovdb/nanovdb/tools/cuda/PointsToGrid.cuh @@ -196,7 +196,7 @@ struct pointer_traits { template __hostdev__ inline static void worldToVoxel(Vec3u8 &voxel, const Vec3T &world, const Map &indexToWorld) { - const Vec3d ijk = indexToWorld.applyInverseMap(world);// world -> index + const Vec3d ijk = Vec3d(indexToWorld.applyInverseMap(world));// world -> index static constexpr double encode = double((1<<8) - 1); voxel[0] = uint8_t( encode*(ijk[0] - math::Floor(ijk[0] + 0.5) + 0.5) ); voxel[1] = uint8_t( encode*(ijk[1] - math::Floor(ijk[1] + 0.5) + 0.5) ); @@ -211,7 +211,7 @@ __hostdev__ inline static void worldToVoxel(Vec3u8 &voxel, const Vec3T &world, c template __hostdev__ inline static void worldToVoxel(Vec3u16 &voxel, const Vec3T &world, const Map &indexToWorld) { - const Vec3d ijk = indexToWorld.applyInverseMap(world);// world -> index + const Vec3d ijk = Vec3d(indexToWorld.applyInverseMap(world));// world -> index static constexpr double encode = double((1<<16) - 1); voxel[0] = uint16_t( encode*(ijk[0] - math::Floor(ijk[0] + 0.5) + 0.5) ); voxel[1] = uint16_t( encode*(ijk[1] - math::Floor(ijk[1] + 0.5) + 0.5) ); @@ -226,7 +226,7 @@ __hostdev__ inline static void worldToVoxel(Vec3u16 &voxel, const Vec3T &world, template __hostdev__ inline static void worldToVoxel(Vec3f &voxel, const Vec3T &world, const Map &indexToWorld) { - const Vec3d ijk = indexToWorld.applyInverseMap(world);// world -> index + const Vec3d ijk = Vec3d(indexToWorld.applyInverseMap(world));// world -> index voxel[0] = float( ijk[0] - math::Floor(ijk[0] + 0.5) ); voxel[1] = float( ijk[1] - math::Floor(ijk[1] + 0.5) ); voxel[2] = float( ijk[2] - math::Floor(ijk[2] + 0.5) ); diff --git a/nanovdb/nanovdb/unittest/TestNanoVDB.cc b/nanovdb/nanovdb/unittest/TestNanoVDB.cc index f26f5b225f..64d36c396c 100644 --- a/nanovdb/nanovdb/unittest/TestNanoVDB.cc +++ b/nanovdb/nanovdb/unittest/TestNanoVDB.cc @@ -1351,6 +1351,705 @@ TEST_F(TestNanoVDB, Vec4) EXPECT_NE(nanovdb::Vec4f(1, 2, 3, 4), nanovdb::Vec4f(1, 2, 3, 5)); }// Vec4 +// ===========================> Math.h refactor tests <=========================== +// The tests below cover behaviour added/changed by the MatBase/VecBase refactors: +// - integer Min/Max precision (no fminf/fmaxf trap) +// - math::Round float/double parity at half-integers +// - Vec2 (no prior test) plus extended Vec3/Vec4 operator coverage +// - All five matrix classes' constructors/operators/transpose/inverse +// - Cross-shape matrix products (Mat2x3*Mat3x2, Mat3*Mat3x2, Mat4*Mat4, ...) +// - Public shape/value introspection on MatBase and VecBase + +TEST_F(TestNanoVDB, IntegerMinMax) +{ + // 2^24 + 1 cannot be represented exactly as float; the old fminf-based + // overload would round both inputs to the same float and return a value + // that is neither of the inputs. + const int32_t a32 = 16777217; // 2^24 + 1 + const int32_t b32 = 16777218; // 2^24 + 2 + EXPECT_EQ(a32, nanovdb::math::Min(a32, b32)); + EXPECT_EQ(b32, nanovdb::math::Max(a32, b32)); + + // INT32_MAX through float can round up to 2^31 and convert back to INT_MIN + // (UB on x86). The ternary path returns the correct value. + const int32_t hi = (std::numeric_limits::max)(); + EXPECT_EQ(hi, nanovdb::math::Max(hi, 0)); + EXPECT_EQ(0, nanovdb::math::Min(hi, 0)); + + // Same risk on the unsigned side: float(UINT_MAX) rounds to 2^32. + const uint32_t uhi = (std::numeric_limits::max)(); + EXPECT_EQ(uhi, nanovdb::math::Max(uhi, 0u)); + EXPECT_EQ(0u, nanovdb::math::Min(uhi, 0u)); + + // Sanity: float/double overloads still work. + EXPECT_EQ(1.0f, nanovdb::math::Min(1.0f, 2.0f)); + EXPECT_EQ(2.0, nanovdb::math::Max(1.0, 2.0)); +}// IntegerMinMax + +TEST_F(TestNanoVDB, Round) +{ + // math::Round(Vec3) used rintf (round-half-to-even) while the + // double overload used floor(x + 0.5). They disagreed at -1.5: + // rintf(-1.5) -> -2, floor(-1) -> -1. + // After the unification both use floor(x + 0.5). + nanovdb::Vec3f vf(-1.5f, -0.5f, 1.5f); + nanovdb::Vec3d vd(-1.5, -0.5, 1.5); + + const auto rf = nanovdb::math::Round(vf); + const auto rd = nanovdb::math::Round(vd); + + EXPECT_EQ(rf[0], rd[0]); + EXPECT_EQ(rf[1], rd[1]); + EXPECT_EQ(rf[2], rd[2]); + + // floor(x + 0.5): -1.5 -> -1, -0.5 -> 0, 1.5 -> 2 + EXPECT_EQ(-1, rf[0]); + EXPECT_EQ( 0, rf[1]); + EXPECT_EQ( 2, rf[2]); +}// Round + +TEST_F(TestNanoVDB, Vec2) +{ + using Vec2d = nanovdb::math::Vec2; + using Vec2i = nanovdb::math::Vec2; + using Coord2 = nanovdb::math::Coord2; + + static_assert(Vec2d::SIZE == 2, ""); + static_assert(Vec2d::size == 2, ""); // openvdb::math::Tuple-compat alias + bool tt = nanovdb::util::is_same::value; + EXPECT_TRUE(tt); + + Vec2d a(1.0, 2.0); + EXPECT_EQ(1.0, a[0]); + EXPECT_EQ(2.0, a[1]); + + // construction from single scalar + EXPECT_EQ(Vec2d(5.0, 5.0), Vec2d(5.0)); + + // copy from Coord2 + Coord2 c(3, 4); + Vec2d fromC(c); + EXPECT_EQ(3.0, fromC[0]); + EXPECT_EQ(4.0, fromC[1]); + + // equality + EXPECT_EQ(Vec2d(1, 2), Vec2d(1, 2)); + EXPECT_NE(Vec2d(1, 2), Vec2d(1, 3)); + + // arithmetic + Vec2d b(3.0, 4.0); + EXPECT_EQ(Vec2d(4, 6), a + b); + EXPECT_EQ(Vec2d(-2, -2), a - b); + EXPECT_EQ(Vec2d(3, 8), a * b); + EXPECT_EQ(Vec2d(2, 2), b / Vec2d(1.5, 2.0)); + EXPECT_EQ(Vec2d(-1, -2), -a); + EXPECT_EQ(Vec2d(2, 4), a * 2.0); + EXPECT_EQ(Vec2d(0.5, 1.0), a / 2.0); + + Vec2d acc = a; acc += b; EXPECT_EQ(Vec2d(4, 6), acc); + acc = a; acc -= b; EXPECT_EQ(Vec2d(-2, -2), acc); + acc = a; acc *= 3.0; EXPECT_EQ(Vec2d(3, 6), acc); + acc = Vec2d(8, 4); acc /= 2.0; EXPECT_EQ(Vec2d(4, 2), acc); + + // mixed Vec2 / Coord2 (these used to take a 3D Coord -- bug #2) + EXPECT_EQ(Vec2d(4, 6), a + c); + EXPECT_EQ(Vec2d(-2, -2), a - c); + acc = a; acc += c; EXPECT_EQ(Vec2d(4, 6), acc); + acc = a; acc -= c; EXPECT_EQ(Vec2d(-2, -2), acc); + + // scalar*Vec free function + EXPECT_EQ(Vec2d(2, 4), 2.0 * a); + + // dot/length + EXPECT_EQ(1.0 + 4.0, a.lengthSqr()); + EXPECT_EQ(std::sqrt(5.0), a.length()); + EXPECT_EQ(1.0*3.0 + 2.0*4.0, a.dot(b)); + + // normalized() — const, non-mutating, unit length + const Vec2d e(3, 4); + EXPECT_NEAR(1.0, e.normalized().length(), 1e-12); + EXPECT_EQ(Vec2d(3, 4), e);// source unchanged + + // min/max scalar reductions + EXPECT_EQ(1.0, a.min()); + EXPECT_EQ(2.0, a.max()); + + // minComponent / maxComponent (mutating) + Vec2d p(5.0, 1.0); + Vec2d q(2.0, 7.0); + Vec2d mn = p; mn.minComponent(q); EXPECT_EQ(Vec2d(2, 1), mn); + Vec2d mx = p; mx.maxComponent(q); EXPECT_EQ(Vec2d(5, 7), mx); + + // floor / ceil / round + Vec2d r(-0.5, 1.5); + EXPECT_EQ(Coord2(-1, 1), r.floor()); + EXPECT_EQ(Coord2( 0, 2), r.ceil()); + // floor(x + 0.5): -0.5 -> 0, 1.5 -> 2 + EXPECT_EQ(Coord2( 0, 2), r.round()); + + // asPointer + Vec2d u(7.0, 9.0); + EXPECT_EQ(7.0, u.asPointer()[0]); + EXPECT_EQ(9.0, u.asPointer()[1]); + + // integer-T division must be per-element (not 1/s, which would give 0) + Vec2i vi(6, 8); + EXPECT_EQ(Vec2i(3, 4), vi / 2); + Vec2i vid = vi; vid /= 2; + EXPECT_EQ(Vec2i(3, 4), vid); +}// Vec2 + +TEST_F(TestNanoVDB, Vec3Ops) +{ + using Vec3d = nanovdb::math::Vec3; + using Vec3i = nanovdb::math::Vec3; + + Vec3d a(1.0, 2.0, 3.0); + Vec3d b(4.0, 5.0, 6.0); + + EXPECT_EQ(Vec3d(5,7,9), a + b); + EXPECT_EQ(Vec3d(3,3,3), b - a); + EXPECT_EQ(Vec3d(4,10,18), a * b); + EXPECT_EQ(Vec3d(2,1,1), b / Vec3d(2.0, 5.0, 6.0)); + EXPECT_EQ(Vec3d(-1,-2,-3), -a); + EXPECT_EQ(Vec3d(2,4,6), a * 2.0); + EXPECT_EQ(Vec3d(0.5,1,1.5), a / 2.0); + EXPECT_EQ(Vec3d(2,4,6), 2.0 * a); + + Vec3d acc = a; acc += b; EXPECT_EQ(Vec3d(5,7,9), acc); + acc = b; acc -= a; EXPECT_EQ(Vec3d(3,3,3), acc); + acc = a; acc *= 2.0; EXPECT_EQ(Vec3d(2,4,6), acc); + acc = Vec3d(8,4,2); acc /= 2.0; EXPECT_EQ(Vec3d(4,2,1), acc); + + EXPECT_EQ(1.0*4 + 2*5 + 3*6, a.dot(b)); + + // cross and outer (3D-only) + Vec3d ex(1,0,0), ey(0,1,0); + EXPECT_EQ(Vec3d(0,0,1), ex.cross(ey)); + auto M = ex.outer(ey); + EXPECT_EQ(0.0, M[0][0]); + EXPECT_EQ(1.0, M[0][1]); + EXPECT_EQ(0.0, M[2][2]); + + // min/max scalar reductions + Vec3d p(3, 1, 4); + EXPECT_EQ(1.0, p.min()); + EXPECT_EQ(4.0, p.max()); + + // minComponent / maxComponent + Vec3d q(2, 7, 1); + Vec3d mn = p; mn.minComponent(q); EXPECT_EQ(Vec3d(2,1,1), mn); + Vec3d mx = p; mx.maxComponent(q); EXPECT_EQ(Vec3d(3,7,4), mx); + + // floor / ceil + Vec3d r(-1.4, 2.7, 3.0); + EXPECT_EQ(nanovdb::Coord(-2, 2, 3), r.floor()); + EXPECT_EQ(nanovdb::Coord(-1, 3, 3), r.ceil()); + // floor(x + 0.5) + Vec3d s(-1.5, -0.5, 1.5); + EXPECT_EQ(nanovdb::Coord(-1, 0, 2), s.round()); + + // mixed Vec3/Coord (3D Coord is correct here) + nanovdb::Coord ijk(1, 2, 3); + EXPECT_EQ(Vec3d(2, 4, 6), a + ijk); + EXPECT_EQ(Vec3d(0, 0, 0), a - ijk); + Vec3d acc2 = a; acc2 += ijk; EXPECT_EQ(Vec3d(2, 4, 6), acc2); + acc2 = a; acc2 -= ijk; EXPECT_EQ(Vec3d(0, 0, 0), acc2); + + // asPointer + const double* dp = a.asPointer(); + EXPECT_EQ(1.0, dp[0]); + EXPECT_EQ(3.0, dp[2]); + + // integer-T division per-element + Vec3i vi(6, 8, 10); + EXPECT_EQ(Vec3i(3, 4, 5), vi / 2); + + // normalize + Vec3d n(3, 0, 4); + n.normalize(); + EXPECT_NEAR(1.0, n.length(), 1e-12); + + // normalized() — const, non-mutating counterpart; callable on a const vector + const Vec3d c(3, 0, 4); + const Vec3d u = c.normalized(); + EXPECT_NEAR(1.0, u.length(), 1e-12); + EXPECT_NEAR(0.6, u[0], 1e-12); + EXPECT_NEAR(0.0, u[1], 1e-12); + EXPECT_NEAR(0.8, u[2], 1e-12); + EXPECT_EQ(Vec3d(3, 0, 4), c);// source unchanged +}// Vec3Ops + +TEST_F(TestNanoVDB, Vec4Ops) +{ + using Vec4d = nanovdb::math::Vec4; + using Vec4i = nanovdb::math::Vec4; + + Vec4d a(1, 2, 3, 4); + Vec4d b(4, 5, 6, 7); + + EXPECT_EQ(Vec4d(5,7,9,11), a + b); + EXPECT_EQ(Vec4d(3,3,3,3), b - a); + EXPECT_EQ(Vec4d(4,10,18,28), a * b); + EXPECT_EQ(Vec4d(-1,-2,-3,-4), -a); + EXPECT_EQ(Vec4d(2,4,6,8), a * 2.0); + EXPECT_EQ(Vec4d(0.5,1,1.5,2), a / 2.0); + EXPECT_EQ(Vec4d(2,4,6,8), 2.0 * a); + + Vec4d acc = a; acc += b; EXPECT_EQ(Vec4d(5,7,9,11), acc); + acc = b; acc -= a; EXPECT_EQ(Vec4d(3,3,3,3), acc); + acc = a; acc *= 2.0; EXPECT_EQ(Vec4d(2,4,6,8), acc); + acc = Vec4d(8,4,2,1); acc /= 2.0; EXPECT_EQ(Vec4d(4,2,1,0.5), acc); + + // ---- Vec4 API gaps that the refactor closed ---- + + // min / max scalar reductions (Vec4 did not have these before) + Vec4d p(3, 1, 4, 1); + EXPECT_EQ(1.0, p.min()); + EXPECT_EQ(4.0, p.max()); + + // asPointer (was missing on Vec4) + const double* dp = a.asPointer(); + EXPECT_EQ(1.0, dp[0]); + EXPECT_EQ(4.0, dp[3]); + + // floor / ceil / round (was missing on Vec4) -> Vec4 + Vec4d r(-1.4, 2.7, 3.0, 0.5); + EXPECT_EQ(nanovdb::math::Vec4(-2, 2, 3, 0), r.floor()); + EXPECT_EQ(nanovdb::math::Vec4(-1, 3, 3, 1), r.ceil()); + // floor(x + 0.5) + EXPECT_EQ(nanovdb::math::Vec4(-1, 3, 3, 1), r.round()); + + // integer-T per-element division + Vec4i vi(6, 8, 10, 12); + EXPECT_EQ(Vec4i(3, 4, 5, 6), vi / 2); + + // normalize + Vec4d n(1, 0, 0, 0); + n.normalize(); + EXPECT_NEAR(1.0, n.length(), 1e-12); + + // normalized() — const, non-mutating counterpart; callable on a const vector + const Vec4d c(0, 3, 0, 4); + const Vec4d u = c.normalized(); + EXPECT_NEAR(1.0, u.length(), 1e-12); + EXPECT_EQ(Vec4d(0, 3, 0, 4), c);// source unchanged +}// Vec4Ops + +TEST_F(TestNanoVDB, Mat2) +{ + using Mat2d = nanovdb::math::Mat2; + + EXPECT_EQ(2, Mat2d::rows()); + EXPECT_EQ(2, Mat2d::cols()); + EXPECT_EQ(4, Mat2d::size()); + + Mat2d a(1, 2, + 3, 4); + EXPECT_EQ(1.0, a[0][0]); EXPECT_EQ(2.0, a[0][1]); + EXPECT_EQ(3.0, a[1][0]); EXPECT_EQ(4.0, a[1][1]); + + // negation, transpose, equality + Mat2d neg = -a; + EXPECT_EQ(-1.0, neg[0][0]); + EXPECT_EQ(-4.0, neg[1][1]); + + Mat2d t = a.transpose(); + EXPECT_EQ(1.0, t[0][0]); EXPECT_EQ(3.0, t[0][1]); + EXPECT_EQ(2.0, t[1][0]); EXPECT_EQ(4.0, t[1][1]); + + EXPECT_TRUE (a == Mat2d(1, 2, 3, 4)); + EXPECT_FALSE(a == Mat2d(1, 2, 3, 5)); + EXPECT_TRUE (a != Mat2d(0, 0, 0, 0)); + + // arithmetic + Mat2d b(5, 6, 7, 8); + Mat2d sum = a + b; EXPECT_EQ(Mat2d(6, 8, 10, 12), sum); + Mat2d dif = b - a; EXPECT_EQ(Mat2d(4, 4, 4, 4), dif); + Mat2d ac = a; ac += b; EXPECT_EQ(Mat2d(6, 8, 10, 12), ac); + ac = b; ac -= a; EXPECT_EQ(Mat2d(4, 4, 4, 4), ac); + + // scalar + Mat2d twoA = a * 2.0; EXPECT_EQ(Mat2d(2, 4, 6, 8), twoA); + Mat2d twoA2 = 2.0 * a; EXPECT_EQ(Mat2d(2, 4, 6, 8), twoA2); + Mat2d halfA = a / 2.0; EXPECT_EQ(Mat2d(0.5, 1, 1.5, 2), halfA); + Mat2d ac2 = a; ac2 *= 0.5; EXPECT_EQ(Mat2d(0.5, 1, 1.5, 2), ac2); + Mat2d ac3 = a; ac3 /= 2.0; EXPECT_EQ(Mat2d(0.5, 1, 1.5, 2), ac3); + + // Mat2 * Mat2 + // [1 2][5 6] = [1*5+2*7 1*6+2*8] = [19 22] + // [3 4][7 8] [3*5+4*7 3*6+4*8] [43 50] + Mat2d prod = a * b; + EXPECT_EQ(Mat2d(19, 22, 43, 50), prod); + + // Mat2 * Vec2 (newly added member) + nanovdb::math::Vec2 v(1, 2); + nanovdb::math::Vec2 mv = a * v; + // [1 2][1] = [5] + // [3 4][2] [11] + EXPECT_EQ(nanovdb::math::Vec2(5, 11), mv); + + // inverse: non-singular roundtrip + singular returns explicit zero + Mat2d inv = a.inverse(); + // det(a) = 1*4 - 2*3 = -2; invDet = -0.5 + // [a^-1] = [-2 1; 1.5 -0.5] + EXPECT_EQ(Mat2d(-2.0, 1.0, 1.5, -0.5), inv); + + Mat2d singular(1, 1, 1, 1); + Mat2d zinv = singular.inverse(); + EXPECT_EQ(Mat2d(0, 0, 0, 0), zinv); +}// Mat2 + +TEST_F(TestNanoVDB, Mat3) +{ + using Mat3d = nanovdb::math::Mat3; + using Vec3d = nanovdb::math::Vec3; + + EXPECT_EQ(3, Mat3d::rows()); + EXPECT_EQ(3, Mat3d::cols()); + EXPECT_EQ(9, Mat3d::size()); + + // Mixed-literal construction now works (Source->T unification) + Mat3d m(1, 2.0, 3, + 4.0, 5, 6.0, + 7, 8.0, 9); + EXPECT_EQ(1.0, m[0][0]); + EXPECT_EQ(5.0, m[1][1]); + EXPECT_EQ(9.0, m[2][2]); + + // transpose + Mat3d mt = m.transpose(); + EXPECT_EQ(1.0, mt[0][0]); EXPECT_EQ(4.0, mt[0][1]); EXPECT_EQ(7.0, mt[0][2]); + EXPECT_EQ(2.0, mt[1][0]); EXPECT_EQ(5.0, mt[1][1]); EXPECT_EQ(8.0, mt[1][2]); + + // ==, !=, +/-, scalar, +=, -=, *=, /= + EXPECT_TRUE (m == m); + EXPECT_FALSE(m != m); + + Mat3d twoM = m * 2.0; + EXPECT_EQ(2.0, twoM[0][0]); + EXPECT_EQ(18.0, twoM[2][2]); + + Mat3d sum = m + m; + EXPECT_EQ(twoM, sum); + + Mat3d dif = sum - m; + EXPECT_EQ(m, dif); + + Mat3d acc = m; acc += m; EXPECT_EQ(twoM, acc); + acc = sum; acc -= m; EXPECT_EQ(m, acc); + acc = m; acc *= 2.0; EXPECT_EQ(twoM, acc); + acc = twoM; acc /= 2.0; EXPECT_EQ(m, acc); + + Mat3d neg = -m; + EXPECT_EQ(-1.0, neg[0][0]); + EXPECT_EQ(-9.0, neg[2][2]); + + // scalar*Mat free function + EXPECT_EQ(twoM, 2.0 * m); + + // Mat3 * Mat3 + Mat3d id(1, 0, 0, 0, 1, 0, 0, 0, 1); + EXPECT_EQ(m, m * id); + Mat3d sq = m * m; + // Row 0, col 0: 1*1 + 2*4 + 3*7 = 30 + EXPECT_EQ(30.0, sq[0][0]); + + // Mat3 * Vec3 + Vec3d v(1, 2, 3); + Vec3d r = m * v; + // Row 0: 1+4+9 = 14, Row 1: 4+10+18 = 32, Row 2: 7+16+27 = 50 + EXPECT_EQ(Vec3d(14, 32, 50), r); +}// Mat3 + +TEST_F(TestNanoVDB, Mat4) +{ + using Mat4d = nanovdb::math::Mat4; + using Vec4d = nanovdb::math::Vec4; + + EXPECT_EQ(4, Mat4d::rows()); + EXPECT_EQ(4, Mat4d::cols()); + EXPECT_EQ(16, Mat4d::size()); + + // Mixed-literal construction (Source->T unification) + Mat4d m( 1, 2.0, 3, 4.0, + 5, 6.0, 7, 8.0, + 9, 10.0, 11, 12.0, + 13, 14.0, 15, 16.0); + EXPECT_EQ(1.0, m[0][0]); + EXPECT_EQ(16.0, m[3][3]); + + // transpose + Mat4d mt = m.transpose(); + EXPECT_EQ(1.0, mt[0][0]); EXPECT_EQ(5.0, mt[0][1]); + EXPECT_EQ(4.0, mt[3][0]); EXPECT_EQ(16.0, mt[3][3]); + + // ==/!=, +/-, +=, -=, scalar, *=, /=, unary minus + EXPECT_TRUE(m == m); + EXPECT_TRUE(m != Mat4d(0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,1)); + Mat4d twoM = m * 2.0; + EXPECT_EQ(32.0, twoM[3][3]); + EXPECT_EQ(twoM, 2.0 * m); + Mat4d sum = m + m; EXPECT_EQ(twoM, sum); + Mat4d dif = sum - m; EXPECT_EQ(m, dif); + Mat4d acc = m; acc += m; EXPECT_EQ(twoM, acc); + acc = sum; acc -= m; EXPECT_EQ(m, acc); + acc = m; acc *= 2.0; EXPECT_EQ(twoM, acc); + acc = twoM; acc /= 2.0; EXPECT_EQ(m, acc); + + Mat4d neg = -m; + EXPECT_EQ(-1.0, neg[0][0]); + + // Mat4 * Mat4 (newly added) + Mat4d id(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1); + EXPECT_EQ(m, m * id); + Mat4d sq = m * m; + // Row 0, col 0: 1*1 + 2*5 + 3*9 + 4*13 = 1 + 10 + 27 + 52 = 90 + EXPECT_EQ(90.0, sq[0][0]); + + // Mat4 * Vec4 (now a member, replacing the old standalone) + Vec4d v(1, 2, 3, 4); + Vec4d r = m * v; + // Row 0: 1+4+9+16 = 30 + // Row 1: 5+12+21+32 = 70 + // Row 2: 9+20+33+48 = 110 + // Row 3: 13+28+45+64 = 150 + EXPECT_EQ(Vec4d(30, 70, 110, 150), r); +}// Mat4 + +TEST_F(TestNanoVDB, Mat2x3_Mat3x2) +{ + using Mat2x3d = nanovdb::math::Mat2x3; + using Mat3x2d = nanovdb::math::Mat3x2; + using Vec2d = nanovdb::math::Vec2; + using Vec3d = nanovdb::math::Vec3; + + EXPECT_EQ(2, Mat2x3d::rows()); EXPECT_EQ(3, Mat2x3d::cols()); EXPECT_EQ(6, Mat2x3d::size()); + EXPECT_EQ(3, Mat3x2d::rows()); EXPECT_EQ(2, Mat3x2d::cols()); EXPECT_EQ(6, Mat3x2d::size()); + + Mat2x3d a(1, 2, 3, + 4, 5, 6); + Mat3x2d b(1, 2, + 3, 4, + 5, 6); + + // transpose round-trip between the two shapes + Mat3x2d at = a.transpose(); + EXPECT_EQ(1.0, at[0][0]); EXPECT_EQ(4.0, at[0][1]); + EXPECT_EQ(2.0, at[1][0]); EXPECT_EQ(5.0, at[1][1]); + EXPECT_EQ(3.0, at[2][0]); EXPECT_EQ(6.0, at[2][1]); + EXPECT_TRUE(a == at.transpose()); + + Mat2x3d bt = b.transpose(); + EXPECT_EQ(1.0, bt[0][0]); EXPECT_EQ(3.0, bt[0][1]); EXPECT_EQ(5.0, bt[0][2]); + EXPECT_EQ(2.0, bt[1][0]); EXPECT_EQ(4.0, bt[1][1]); EXPECT_EQ(6.0, bt[1][2]); + + // element-wise (only Mat2x3 had + and += previously; Mat3x2 now has them too) + Mat2x3d a2 = a + a; EXPECT_EQ(Mat2x3d(2,4,6,8,10,12), a2); + Mat2x3d ac = a; ac += a; EXPECT_EQ(a2, ac); + Mat2x3d ad = a; ad -= a; EXPECT_EQ(Mat2x3d(0,0,0,0,0,0), ad); + + Mat3x2d b2 = b + b; EXPECT_EQ(Mat3x2d(2,4,6,8,10,12), b2); + Mat3x2d bc = b; bc += b; EXPECT_EQ(b2, bc); + Mat3x2d bd = b; bd -= b; EXPECT_EQ(Mat3x2d(0,0,0,0,0,0), bd); + + // scalar + EXPECT_EQ(Mat2x3d(2,4,6,8,10,12), a * 2.0); + EXPECT_EQ(Mat2x3d(2,4,6,8,10,12), 2.0 * a); + EXPECT_EQ(a, (a * 2.0) / 2.0); + + EXPECT_EQ(Mat3x2d(2,4,6,8,10,12), b * 2.0); + EXPECT_EQ(Mat3x2d(2,4,6,8,10,12), 2.0 * b); + + // unary minus + EXPECT_EQ(Mat2x3d(-1,-2,-3,-4,-5,-6), -a); + EXPECT_EQ(Mat3x2d(-1,-2,-3,-4,-5,-6), -b); + + // ==, != + EXPECT_TRUE(a == Mat2x3d(1,2,3,4,5,6)); + EXPECT_TRUE(a != Mat2x3d(1,2,3,4,5,7)); + + // Mat2x3 * Vec3 (newly added member) -> Vec2 + Vec3d v3(1, 2, 3); + Vec2d r1 = a * v3; + // Row 0: 1+4+9 = 14, Row 1: 4+10+18 = 32 + EXPECT_EQ(Vec2d(14, 32), r1); + + // Mat3x2 * Vec2 (newly added member) -> Vec3 + Vec2d v2(1, 2); + Vec3d r2 = b * v2; + // Row 0: 1+4 = 5, Row 1: 3+8 = 11, Row 2: 5+12 = 17 + EXPECT_EQ(Vec3d(5, 11, 17), r2); +}// Mat2x3_Mat3x2 + +TEST_F(TestNanoVDB, MatMul) +{ + using namespace nanovdb::math; + using Mat2d = Mat2; + using Mat3d = Mat3; + using Mat4d = Mat4; + using Mat2x3d = Mat2x3; + using Mat3x2d = Mat3x2; + + Mat2d m22(1, 2, 3, 4); + Mat3d m33(1, 2, 3, 4, 5, 6, 7, 8, 9); + Mat4d m44(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16); + Mat2x3d m23(1, 2, 3, + 4, 5, 6); + Mat3x2d m32(1, 2, + 3, 4, + 5, 6); + + // Square * Square (member form) + EXPECT_EQ(Mat2d(7, 10, 15, 22), m22 * m22); + // [1 2 3][1 2 3] [30 36 42] + // [4 5 6][4 5 6] = [66 81 96] + // [7 8 9][7 8 9] [102 126 150] + Mat3d m33sq = m33 * m33; + EXPECT_EQ(30.0, m33sq[0][0]); + EXPECT_EQ(150.0, m33sq[2][2]); + + Mat4d m44sq = m44 * m44; + EXPECT_EQ(90.0, m44sq[0][0]); // 1+10+27+52 + EXPECT_EQ(600.0, m44sq[3][3]); // 13*4 + 14*8 + 15*12 + 16*16 + + // Rectangular: Mat2x3 * Mat3x2 -> Mat2 + // [1 2 3][1 2] [22 28] + // [4 5 6][3 4] = [49 64] + // [5 6] + EXPECT_EQ(Mat2d(22, 28, 49, 64), m23 * m32); + + // Mat3x2 * Mat2x3 -> Mat3 + // Each entry: lhs.row[i] dotted with rhs.col[j] + // Row 0, col 0: 1*1 + 2*4 = 9 + Mat3d m32x23 = m32 * m23; + EXPECT_EQ(9.0, m32x23[0][0]); + EXPECT_EQ(51.0, m32x23[2][2]); // 5*3 + 6*6 = 51 + + // Mat2 * Mat2x3 -> Mat2x3 (identity-passes test) + Mat2d id2(1, 0, 0, 1); + EXPECT_EQ(m23, id2 * m23); + + // Mat3x2 * Mat2 -> Mat3x2 + EXPECT_EQ(m32, m32 * id2); + + // Mat2x3 * Mat3 -> Mat2x3 + Mat3d id3(1, 0, 0, 0, 1, 0, 0, 0, 1); + EXPECT_EQ(m23, m23 * id3); + + // Mat3 * Mat3x2 -> Mat3x2 (the consistency-pass addition) + Mat3x2d r = m33 * m32; + // Row 0: [1*1+2*3+3*5 1*2+2*4+3*6] = [22 28] + EXPECT_EQ(22.0, r[0][0]); EXPECT_EQ(28.0, r[0][1]); + EXPECT_EQ(49.0, r[1][0]); EXPECT_EQ(64.0, r[1][1]); + EXPECT_EQ(76.0, r[2][0]); EXPECT_EQ(100.0, r[2][1]); + + // identity-on-the-left passes through for the rectangular shapes too + EXPECT_EQ(m32, id3 * m32); + + // Mat3 * Vec3 (member) + Vec3 v3(1, 2, 3); + EXPECT_EQ(Vec3(14, 32, 50), m33 * v3); + + // Mat4 * Vec4 (member after the consistency pass) + Vec4 v4(1, 2, 3, 4); + EXPECT_EQ(Vec4(30, 70, 110, 150), m44 * v4); + + // Mat2 * Vec2 (member, new) + EXPECT_EQ(Vec2(5, 11), m22 * Vec2(1, 2)); + + // Mat2x3 * Vec3 (member, new) -> Vec2 + EXPECT_EQ(Vec2(14, 32), m23 * v3); + + // Mat3x2 * Vec2 (member, new) -> Vec3 + EXPECT_EQ(Vec3(5, 11, 17), m32 * Vec2(1, 2)); +}// MatMul + +TEST_F(TestNanoVDB, MatVecIntrospection) +{ + using namespace nanovdb::math; + + // MatBase / Mat* shape and type introspection (public after refactor) + static_assert(Mat2::rows() == 2, ""); + static_assert(Mat2::cols() == 2, ""); + static_assert(Mat2::size() == 4, ""); + static_assert(Mat3::rows() == 3, ""); + static_assert(Mat3::cols() == 3, ""); + static_assert(Mat3::size() == 9, ""); + static_assert(Mat4::rows() == 4, ""); + static_assert(Mat4::cols() == 4, ""); + static_assert(Mat4::size() == 16, ""); + static_assert(Mat2x3::rows() == 2, ""); + static_assert(Mat2x3::cols() == 3, ""); + static_assert(Mat3x2::rows() == 3, ""); + static_assert(Mat3x2::cols() == 2, ""); + + bool tm = nanovdb::util::is_same::ValueType>::value; + EXPECT_TRUE(tm); + tm = nanovdb::util::is_same::ValueType>::value; + EXPECT_TRUE(tm); + + // VecBase / Vec* shape and type introspection + static_assert(Vec2::SIZE == 2, ""); + static_assert(Vec3::SIZE == 3, ""); + static_assert(Vec4::SIZE == 4, ""); + // openvdb::math::Tuple-compat lowercase alias + static_assert(Vec3::size == 3, ""); + + bool tv = nanovdb::util::is_same::ValueType>::value; + EXPECT_TRUE(tv); + tv = nanovdb::util::is_same::ValueType>::value; + EXPECT_TRUE(tv); + + // data() on Mat lays out row-major + Mat3 m(1, 2, 3, 4, 5, 6, 7, 8, 9); + EXPECT_EQ(1.0, m.data()[0]); + EXPECT_EQ(5.0, m.data()[4]); + EXPECT_EQ(9.0, m.data()[8]); + + // asPointer() on Vec lays out left-to-right + Vec4 v(10, 20, 30, 40); + EXPECT_EQ(10.0, v.asPointer()[0]); + EXPECT_EQ(40.0, v.asPointer()[3]); + + // Named component accessors x()/y()/z()/w(). x() is available on every + // Vec*; y() requires N >= 2; z() requires N >= 3; w() requires N >= 4. + // Each is exercised here at compile time (constexpr) and at runtime + // (read + write through the non-const overload). + { + constexpr Vec2 v2(11, 22); + static_assert(v2.x() == 11, ""); + static_assert(v2.y() == 22, ""); + } + { + constexpr Vec3 v3(31, 32, 33); + static_assert(v3.x() == 31, ""); + static_assert(v3.y() == 32, ""); + static_assert(v3.z() == 33, ""); + } + { + constexpr Vec4 v4(41, 42, 43, 44); + static_assert(v4.x() == 41, ""); + static_assert(v4.y() == 42, ""); + static_assert(v4.z() == 43, ""); + static_assert(v4.w() == 44, ""); + } + { + // Non-const overloads return mutable references. + Vec4 vw(0.0, 0.0, 0.0, 0.0); + vw.x() = 1.0; vw.y() = 2.0; vw.z() = 3.0; vw.w() = 4.0; + EXPECT_EQ(1.0, vw[0]); + EXPECT_EQ(2.0, vw[1]); + EXPECT_EQ(3.0, vw[2]); + EXPECT_EQ(4.0, vw[3]); + } + // The static_asserts in the bodies of y()/z()/w() guarantee that calling + // e.g. Vec2::w() fails with a clear compile-time error + // ("VecBase::w() requires N >= 4"); we don't exercise the negative + // case here because there's no way to assert "must fail to compile" + // inside a gtest TU. +}// MatVecIntrospection + TEST_F(TestNanoVDB, Map) { EXPECT_EQ(264u, sizeof(nanovdb::Map)); diff --git a/nanovdb/nanovdb/unittest/TestNanoVDB.cu b/nanovdb/nanovdb/unittest/TestNanoVDB.cu index f9d7a7f664..5a30356b74 100644 --- a/nanovdb/nanovdb/unittest/TestNanoVDB.cu +++ b/nanovdb/nanovdb/unittest/TestNanoVDB.cu @@ -2075,7 +2075,7 @@ TEST(TestNanoVDBCUDA, Sphere_CudaPointsToGrid_Voxel32) EXPECT_LE(voxel[0], 0.5f); EXPECT_LE(voxel[1], 0.5f); EXPECT_LE(voxel[2], 0.5f); - test = (begin[i] - nanovdb::voxelToWorld(voxel, ijk, grid->map())).length() < 1e-9; + test = (begin[i] - nanovdb::voxelToWorld(voxel, ijk, grid->map())).length() < 1e-9; } EXPECT_TRUE(test); } @@ -2202,7 +2202,7 @@ TEST(TestNanoVDBCUDA, Sphere_CudaPointsToGrid_Voxel16) EXPECT_LE(count, maxPointsPerVoxel); bool test = false; for (uint64_t j=0; test == false && jmap())).length() < 1e-6; + test = (begin[i] - nanovdb::voxelToWorld(start[j], ijk, grid->map())).length() < 1e-6; } } }); @@ -2331,7 +2331,7 @@ TEST(TestNanoVDBCUDA, Sphere_CudaPointsToGrid_Voxel8) EXPECT_LE(count, maxPointsPerVoxel); bool test = false; for (uint64_t j=0; test == false && jmap())).length() < 1e-2; + test = (begin[i] - nanovdb::voxelToWorld(start[j], ijk, grid->map())).length() < 1e-2; } EXPECT_TRUE(test); } diff --git a/nanovdb/nanovdb/unittest/TestOpenVDB.cc b/nanovdb/nanovdb/unittest/TestOpenVDB.cc index fd47a7e017..838b9f10ef 100644 --- a/nanovdb/nanovdb/unittest/TestOpenVDB.cc +++ b/nanovdb/nanovdb/unittest/TestOpenVDB.cc @@ -2243,9 +2243,9 @@ TEST_F(TestOpenVDB, NanoToOpenVDB_ValueOnIndex_Vec3f_SideCar) EXPECT_EQ("vec3s", openGrid->valueType()); EXPECT_EQ(2u, openGrid->activeVoxelCount()); auto openAcc = openGrid->getAccessor(); - EXPECT_EQ(nanovdb::Vec3f(0.0f), openAcc.getValue(openvdb::Coord(0, 0, 0))); - EXPECT_EQ(nanovdb::Vec3f(1.0f), openAcc.getValue(openvdb::Coord(1, 2, 3))); - EXPECT_EQ(nanovdb::Vec3f(2.0f), openAcc.getValue(openvdb::Coord(2, -2, 9))); + EXPECT_EQ(nanovdb::Vec3f(0.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(0, 0, 0)))); + EXPECT_EQ(nanovdb::Vec3f(1.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(1, 2, 3)))); + EXPECT_EQ(nanovdb::Vec3f(2.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(2, -2, 9)))); const auto nanoBBox = idxGrid->indexBBox(); const auto openBBox = openGrid->evalActiveVoxelBoundingBox(); @@ -2263,9 +2263,9 @@ TEST_F(TestOpenVDB, NanoToOpenVDB_ValueOnIndex_Vec3f_SideCar) EXPECT_EQ("vec3s", openGrid->valueType()); EXPECT_EQ(2u, openGrid->activeVoxelCount()); auto openAcc = openGrid->getAccessor(); - EXPECT_EQ(nanovdb::Vec3f(0.0f), openAcc.getValue(openvdb::Coord(0, 0, 0))); - EXPECT_EQ(nanovdb::Vec3f(1.0f), openAcc.getValue(openvdb::Coord(1, 2, 3))); - EXPECT_EQ(nanovdb::Vec3f(2.0f), openAcc.getValue(openvdb::Coord(2, -2, 9))); + EXPECT_EQ(nanovdb::Vec3f(0.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(0, 0, 0)))); + EXPECT_EQ(nanovdb::Vec3f(1.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(1, 2, 3)))); + EXPECT_EQ(nanovdb::Vec3f(2.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(2, -2, 9)))); const auto nanoBBox = idxGrid->indexBBox(); const auto openBBox = openGrid->evalActiveVoxelBoundingBox(); @@ -2323,9 +2323,9 @@ TEST_F(TestOpenVDB, NanoToOpenVDB_ValueIndex_Vec3f_SideCar) EXPECT_EQ("vec3s", openGrid->valueType()); EXPECT_EQ(2u, openGrid->activeVoxelCount()); auto openAcc = openGrid->getAccessor(); - EXPECT_EQ(nanovdb::Vec3f(0.0f), openAcc.getValue(openvdb::Coord(0, 0, 0))); - EXPECT_EQ(nanovdb::Vec3f(1.0f), openAcc.getValue(openvdb::Coord(1, 2, 3))); - EXPECT_EQ(nanovdb::Vec3f(2.0f), openAcc.getValue(openvdb::Coord(2, -2, 9))); + EXPECT_EQ(nanovdb::Vec3f(0.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(0, 0, 0)))); + EXPECT_EQ(nanovdb::Vec3f(1.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(1, 2, 3)))); + EXPECT_EQ(nanovdb::Vec3f(2.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(2, -2, 9)))); const auto nanoBBox = idxGrid->indexBBox(); const auto openBBox = openGrid->evalActiveVoxelBoundingBox(); @@ -2343,9 +2343,9 @@ TEST_F(TestOpenVDB, NanoToOpenVDB_ValueIndex_Vec3f_SideCar) EXPECT_EQ("vec3s", openGrid->valueType()); EXPECT_EQ(2u, openGrid->activeVoxelCount()); auto openAcc = openGrid->getAccessor(); - EXPECT_EQ(nanovdb::Vec3f(0.0f), openAcc.getValue(openvdb::Coord(0, 0, 0))); - EXPECT_EQ(nanovdb::Vec3f(1.0f), openAcc.getValue(openvdb::Coord(1, 2, 3))); - EXPECT_EQ(nanovdb::Vec3f(2.0f), openAcc.getValue(openvdb::Coord(2, -2, 9))); + EXPECT_EQ(nanovdb::Vec3f(0.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(0, 0, 0)))); + EXPECT_EQ(nanovdb::Vec3f(1.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(1, 2, 3)))); + EXPECT_EQ(nanovdb::Vec3f(2.0f), nanovdb::Vec3f(openAcc.getValue(openvdb::Coord(2, -2, 9)))); const auto nanoBBox = idxGrid->indexBBox(); const auto openBBox = openGrid->evalActiveVoxelBoundingBox(); diff --git a/pendingchanges/nanovdbmath.txt b/pendingchanges/nanovdbmath.txt new file mode 100644 index 0000000000..2a5a9a4896 --- /dev/null +++ b/pendingchanges/nanovdbmath.txt @@ -0,0 +1,20 @@ +NanoVDB: + Improvements: + - Refactored the math::Mat and math::Vec classes onto shared base classes + (MatBase and the new VecBase), giving every Mat/Vec type a complete and + consistent operator surface, bounds-checked element access, and constexpr + evaluation throughout Math.h and the Map API in NanoVDB.h. The + cross-template Vec conversion constructors are now explicit, so downstream + code relying on implicit Vec type conversions may need a one-line cast. + + Bug Fixes: + - Fixed math::Min/Max for 32-bit integers routing through fminf/fmaxf, which + lost precision above 2^24 and could overflow near INT_MAX. + - Fixed integer division in the Vec and Mat scalar operator/, which previously + computed a reciprocal first, so e.g. Vec3i(6,8,10) / 2 returned (0,0,0). + - Fixed Mat2::inverse returning uninitialized values on singular input; it now + returns zero. + - Fixed Vec2 arithmetic against Coord silently dropping the z component; the + overloads now take Coord2. + - Fixed math::Round(Vec3) disagreeing with the double overload at + half-integers; both now use floor(x + 0.5).