diff --git a/.mailmap b/.mailmap index fdd425c3a1..5f5e8e3031 100644 --- a/.mailmap +++ b/.mailmap @@ -15,10 +15,12 @@ Benjamin Curtice Corbett Benjamin Corbett Ben Corbett <32752943+corbett5@users.noreply.github.com> Brian Manh Hien Han Brian Han Brian Manh Hien Han Brian Manh Hien Han +Brian T.N. Gunney Brian Gunney <45609916+gunney1@users.noreply.github.com> Chris White Christopher A. White Cyrus D. Harrison Cyrus Harrison Cyrus D. Harrison Cyrus Harrison Cyrus D. Harrison Cyrus +Daniel Taller Danny Taller <66029857+dtaller@users.noreply.github.com> Esteban Pauli Esteban Pauli <40901502+estebanpauli@users.noreply.github.com> Evan Taylor Desantola Evan Taylor DeSantola George Zagaris George Zagaris diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index b01e8460c2..d175c020f1 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -88,6 +88,10 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ - Adds an overload to quest's `SignedDistance` query to return the closest point on the surface to the query point and the surface normal at that point. Also exposes this functionality in quest's signed_distance C API. +- Adds utility function for linear interpolation (`lerp`) of two numbers +- Adds utility function to compute binomial coefficients +- Adds a `CurvedPolygon` class to primal representing a polygon with `BezierCurves` as edges +- Adds functions to compute the moments (area, centroid) of a `CurvedPolygon` ### Changed - Axom now requires C++14 and will default to that if not specified via `BLT_CXX_STD`. @@ -140,6 +144,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ and added a new derived class `sidre::IndexedCollection` - Spin: `BVH::findPoints/Rays/BoundingBoxes()` candidate search methods now accept an `axom::ArrayView` for the `offsets` and `counts` output arrays, and return `candidates` as an `axom::Array`. +- Renamed `primal::Polygon::centroid()` to `primal::Polygon::vertexMean()` because it was not actually computing the centroid. ### Fixed - Fixed a bug relating to swap and assignment operations for multidimensional `axom::Array`s diff --git a/src/axom/core/numerics/Matrix.hpp b/src/axom/core/numerics/Matrix.hpp index c2a071a7ad..901148c8d0 100644 --- a/src/axom/core/numerics/Matrix.hpp +++ b/src/axom/core/numerics/Matrix.hpp @@ -117,6 +117,11 @@ template class Matrix { public: + /*! + * \brief Default constructor + */ + Matrix() : m_rows(0), m_cols(0), m_data(nullptr), m_usingExternal(false) { } + /*! * \brief Constructor, creates a Matrix with the given rows and columns. * @@ -527,12 +532,6 @@ class Matrix /// @} private: - /*! - * \brief Default constructor. Does nothing. - * \note Made private to prevent host-code from calling this. - */ - Matrix() : m_rows(0), m_cols(0), m_data(nullptr) {}; - /// \name Private Helper Methods /// @{ diff --git a/src/axom/core/tests/utils_utilities.hpp b/src/axom/core/tests/utils_utilities.hpp index 7bf48e025e..287cabacf4 100644 --- a/src/axom/core/tests/utils_utilities.hpp +++ b/src/axom/core/tests/utils_utilities.hpp @@ -158,3 +158,139 @@ TEST(utils_utilities, floor_ceil) EXPECT_EQ(-5.0, axom::utilities::ceil(val)); } } + +//------------------------------------------------------------------------------ +TEST(utils_utilities, binomial_coefficient) +{ + std::cout << "Testing binomial coefficient function." << std::endl; + + // test n less than zero + { + const int n = -1; + const int exp = 0; + for(int k = -1; k < 10; ++k) + { + auto binom_k_n = axom::utilities::binomialCoefficient(n, k); + EXPECT_EQ(exp, binom_k_n); + } + } + + // test n := 0 + { + const int n = 0; + + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 0)); + + EXPECT_EQ(0, axom::utilities::binomialCoefficient(n, -1)); + EXPECT_EQ(0, axom::utilities::binomialCoefficient(n, 1)); + } + + // test n := 1 + { + const int n = 1; + + EXPECT_EQ(0, axom::utilities::binomialCoefficient(n, -1)); + + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 0)); + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 1)); + + EXPECT_EQ(0, axom::utilities::binomialCoefficient(n, 2)); + } + + // test n := 2 + { + const int n = 2; + + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 0)); + EXPECT_EQ(2, axom::utilities::binomialCoefficient(n, 1)); + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 2)); + } + + // test n := 3 + { + const int n = 3; + + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 0)); + EXPECT_EQ(3, axom::utilities::binomialCoefficient(n, 1)); + EXPECT_EQ(3, axom::utilities::binomialCoefficient(n, 2)); + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 3)); + } + + // test n := 4 + { + const int n = 4; + + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 0)); + EXPECT_EQ(4, axom::utilities::binomialCoefficient(n, 1)); + EXPECT_EQ(6, axom::utilities::binomialCoefficient(n, 2)); + EXPECT_EQ(4, axom::utilities::binomialCoefficient(n, 3)); + EXPECT_EQ(1, axom::utilities::binomialCoefficient(n, 4)); + } + + // test recurrence relation nCk = (n-1)C(k-1) + (n-1)C(k) + { + for(int n = 1; n < 10; ++n) + { + for(int k = 1; k <= n; ++k) + { + auto binom_n_k = axom::utilities::binomialCoefficient(n, k); + auto binom_n1_k1 = axom::utilities::binomialCoefficient(n - 1, k - 1); + auto binom_n1_k = axom::utilities::binomialCoefficient(n - 1, k); + + EXPECT_EQ(binom_n_k, binom_n1_k1 + binom_n1_k); + } + } + } +} + +TEST(utils_utilities, lerp) +{ + std::cout << "Testing lerp function." << std::endl; + + { + double A = 0.; + double B = 1.; + + EXPECT_DOUBLE_EQ(0., axom::utilities::lerp(A, B, 0.)); + EXPECT_DOUBLE_EQ(1., axom::utilities::lerp(A, B, 1.)); + EXPECT_DOUBLE_EQ(.5, axom::utilities::lerp(A, B, .5)); + + for(double t = -2.0; t < 2.0; t += 0.05) + { + EXPECT_DOUBLE_EQ(t, axom::utilities::lerp(A, B, t)); + } + } + + // Test endpoint interpolation + { + const double lower = -.23; // Arbitrary end points + const double upper = 5.73; + for(int i = 0; i < 100; ++i) + { + double A = axom::utilities::random_real(lower, upper); + double B = axom::utilities::random_real(lower, upper); + + EXPECT_DOUBLE_EQ(A, axom::utilities::lerp(A, B, 0.)); + EXPECT_DOUBLE_EQ(B, axom::utilities::lerp(B, A, 0.)); + EXPECT_DOUBLE_EQ(B, axom::utilities::lerp(A, B, 1.)); + EXPECT_DOUBLE_EQ(A, axom::utilities::lerp(B, A, 1.)); + } + } + + // Compute using different form + { + const double lower = -.23; // Arbitrary end points + const double upper = 5.73; + for(int i = 0; i < 100; ++i) + { + double A = axom::utilities::random_real(lower, upper); + double B = axom::utilities::random_real(lower, upper); + + // Test interpolation and also extrapolation beyond endpoints. + double t = axom::utilities::random_real(-1.5, 1.5); + + double exp = A + (B - A) * t; + EXPECT_NEAR(exp, axom::utilities::lerp(A, B, t), 1e-12); + } + } +} diff --git a/src/axom/core/utilities/Utilities.cpp b/src/axom/core/utilities/Utilities.cpp index d64b28363b..88234d4072 100644 --- a/src/axom/core/utilities/Utilities.cpp +++ b/src/axom/core/utilities/Utilities.cpp @@ -4,8 +4,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) /*! - * - * \file + * \file Utilities.cpp * * \brief Implementation file for utility functions. * @@ -39,5 +38,29 @@ void processAbort() #endif } +int binomialCoefficient(int n, int k) +{ + if(k > n || k < 0) // check if out-of-bounds + { + return 0; + } + if(k == n || k == 0) // early return + { + return 1; + } + if(k > n - k) // exploit symmetry to reduce work + { + k = n - k; + } + + int val = 1; + for(int i = 1; i <= k; ++i) + { + val *= (n - k + i); + val /= i; + } + return val; +} + } // end namespace utilities } // end namespace axom diff --git a/src/axom/core/utilities/Utilities.hpp b/src/axom/core/utilities/Utilities.hpp index 1c7b86bcbf..b74234b921 100644 --- a/src/axom/core/utilities/Utilities.hpp +++ b/src/axom/core/utilities/Utilities.hpp @@ -109,6 +109,15 @@ inline AXOM_HOST_DEVICE void swap(T& a, T& b) b = tmp; } +/*! + * \brief returns the linear interpolation of \a A and \a B at \a t. i.e. (1-t)A+tB + */ +template +inline AXOM_HOST_DEVICE T lerp(T A, T B, T t) +{ + return (1 - t) * A + t * B; +} + /*! * \brief Returns the base 2 logarithm of the input. * \param [in] val The input value @@ -163,6 +172,14 @@ inline AXOM_HOST_DEVICE T clampLower(T val, T lower) return val < lower ? lower : val; } +/*! + * \brief Computes the binomial coefficient `n choose k` + * + * \return \f$ {n\choose k} = n! / (k! * (n-k)!)\f$ + * when \f$ n \ge k \ge 0 \f$, 0 otherwise. + */ +int binomialCoefficient(int n, int k); + /*! * \brief Returns a random real number within the specified interval * diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index 83a7a175eb..920fd4a3b8 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -20,6 +20,7 @@ set( primal_headers ## geometry geometry/BezierCurve.hpp geometry/BoundingBox.hpp + geometry/CurvedPolygon.hpp geometry/OrientedBoundingBox.hpp geometry/OrientationResult.hpp geometry/NumericArray.hpp @@ -42,14 +43,16 @@ set( primal_headers operators/orientation.hpp operators/squared_distance.hpp operators/compute_bounding_box.hpp + operators/compute_moments.hpp operators/in_sphere.hpp operators/split.hpp operators/detail/clip_impl.hpp + operators/detail/compute_moments_impl.hpp operators/detail/intersect_bezier_impl.hpp - operators/detail/intersect_ray_impl.hpp operators/detail/intersect_bounding_box_impl.hpp operators/detail/intersect_impl.hpp + operators/detail/intersect_ray_impl.hpp ## utils utils/ZipIndexable.hpp diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index 0ce440dfd7..76d1a39781 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/src/axom/primal/geometry/BezierCurve.hpp @@ -12,7 +12,7 @@ #ifndef AXOM_PRIMAL_BEZIERCURVE_HPP_ #define AXOM_PRIMAL_BEZIERCURVE_HPP_ -#include "axom/core/utilities/Utilities.hpp" +#include "axom/core.hpp" #include "axom/slic.hpp" #include "axom/primal/geometry/NumericArray.hpp" @@ -50,7 +50,6 @@ std::ostream& operator<<(std::ostream& os, const BezierCurve& bCurve); * The curve is approximated by the control points, * parametrized from t=0 to t=1. */ - template class BezierCurve { @@ -63,6 +62,12 @@ class BezierCurve using BoundingBoxType = BoundingBox; using OrientedBoundingBoxType = OrientedBoundingBox; + AXOM_STATIC_ASSERT_MSG((NDIMS == 2) || (NDIMS == 3), + "A Bezier Curve object may be defined in 2-D or 3-D"); + AXOM_STATIC_ASSERT_MSG( + std::is_arithmetic::value, + "A Bezier Curve must be defined using an arithmetic type"); + public: /*! * \brief Constructor for a Bezier Curve that reserves space for @@ -73,12 +78,6 @@ class BezierCurve */ explicit BezierCurve(int ord = -1) { - AXOM_STATIC_ASSERT_MSG( - (NDIMS == 2) || (NDIMS == 3), - "A Bezier Curve object may be defined in 2-D or 3-D"); - AXOM_STATIC_ASSERT_MSG( - std::is_arithmetic::value, - "A Bezier Curve must be defined using an arithmetic type"); SLIC_ASSERT(ord >= -1); const int sz = utilities::max(-1, ord + 1); m_controlPoints.resize(sz); @@ -98,12 +97,6 @@ class BezierCurve */ BezierCurve(T* pts, int ord) { - AXOM_STATIC_ASSERT_MSG( - (NDIMS == 2) || (NDIMS == 3), - "A Bezier Curve object may be defined in 2-D or 3-D"); - AXOM_STATIC_ASSERT_MSG( - std::is_arithmetic::value, - "A Bezier Curve must be defined using an arithmetic type"); SLIC_ASSERT(pts != nullptr); SLIC_ASSERT(ord >= 0); @@ -128,15 +121,8 @@ class BezierCurve * \pre order is greater than or equal to zero * */ - BezierCurve(PointType* pts, int ord) { - AXOM_STATIC_ASSERT_MSG( - (NDIMS == 2) || (NDIMS == 3), - "A Bezier Curve object may be defined in 2-D or 3-D"); - AXOM_STATIC_ASSERT_MSG( - std::is_arithmetic::value, - "A Bezier Curve must be defined using an arithmetic type"); SLIC_ASSERT(pts != nullptr); SLIC_ASSERT(ord >= 0); @@ -149,13 +135,30 @@ class BezierCurve } } - /*! Sets the order of the Bezier Curve*/ + /*! + * \brief Constructor for a Bezier Curve from an vector of coordinates + * + * \param [in] pts a vector with ord+1 control points + * \param [in] ord The Curve's polynomial order + * \pre order is greater than or equal to zero + * + */ + BezierCurve(const std::vector& pts, int ord) + { + SLIC_ASSERT(ord >= 0); + + const int sz = utilities::max(0, ord + 1); + m_controlPoints.resize(sz); + m_controlPoints = pts; + } + + /// Sets the order of the Bezier Curve void setOrder(int ord) { m_controlPoints.resize(ord + 1); } - /*! Returns the order of the Bezier Curve*/ + /// Returns the order of the Bezier Curve int getOrder() const { return static_cast(m_controlPoints.size()) - 1; } - /*! Clears the list of control points*/ + /// Clears the list of control points void clear() { const int ord = getOrder(); @@ -165,13 +168,13 @@ class BezierCurve } } - /*! Retrieves the control point at index \a idx */ + /// Retrieves the control point at index \a idx PointType& operator[](int idx) { return m_controlPoints[idx]; } - /*! Retrieves the control point at index \a idx */ + /// Retrieves the control point at index \a idx const PointType& operator[](int idx) const { return m_controlPoints[idx]; } - /* Checks equality of two Bezier Curve */ + /// Checks equality of two Bezier Curve friend inline bool operator==(const BezierCurve& lhs, const BezierCurve& rhs) { @@ -184,17 +187,28 @@ class BezierCurve return !(lhs == rhs); } - /*! Returns a copy of the Bezier curve's control points */ + /// Returns a copy of the Bezier curve's control points CoordsVec getControlPoints() const { return m_controlPoints; } - /*! Returns an axis-aligned bounding box containing the Bezier curve */ + /// Reverses the order of the Bezier curve's control points + void reverseOrientation() + { + const int ord = getOrder(); + const int mid = (ord + 1) / 2; + for(int i = 0; i < mid; ++i) + { + axom::utilities::swap(m_controlPoints[i], m_controlPoints[ord - i]); + } + } + + /// Returns an axis-aligned bounding box containing the Bezier curve BoundingBoxType boundingBox() const { return BoundingBoxType(m_controlPoints.data(), static_cast(m_controlPoints.size())); } - /*! Returns an oriented bounding box containing the Bezier curve */ + /// Returns an oriented bounding box containing the Bezier curve OrientedBoundingBoxType orientedBoundingBox() const { return OrientedBoundingBoxType(m_controlPoints.data(), @@ -209,9 +223,10 @@ class BezierCurve * * \note We typically evaluate the curve at \a t between 0 and 1 */ - PointType evaluate(T t) const { + using axom::utilities::lerp; + PointType ptval; const int ord = getOrder(); @@ -230,7 +245,7 @@ class BezierCurve const int end = ord - p; for(int k = 0; k <= end; ++k) { - dCarray[k] = (1 - t) * dCarray[k] + t * dCarray[k + 1]; + dCarray[k] = lerp(dCarray[k], dCarray[k + 1], t); } } ptval[i] = dCarray[0]; @@ -240,11 +255,52 @@ class BezierCurve } /*! - * \brief Splits a Bezier curve into two Bezier curves at particular parameter - * value between 0 and 1 + * \brief Computes the tangent of a Bezier curve at a particular parameter value \a t + * + * \param [in] t parameter value at which to compute tangent + * \return p the tangent vector of the Bezier curve at t + * + * \note We typically find the tangent of the curve at \a t between 0 and 1 + */ + VectorType dt(T t) const + { + using axom::utilities::lerp; + VectorType val; + + const int ord = getOrder(); + std::vector dCarray(ord + 1); + + // Run de Casteljau algorithm on each dimension + for(int i = 0; i < NDIMS; ++i) + { + for(int p = 0; p <= ord; ++p) + { + dCarray[p] = m_controlPoints[p][i]; + } + + // stop one step early and take difference of last two values + for(int p = 1; p <= ord - 1; ++p) + { + const int end = ord - p; + for(int k = 0; k <= end; ++k) + { + dCarray[k] = lerp(dCarray[k], dCarray[k + 1], t); + } + } + val[i] = ord * (dCarray[1] - dCarray[0]); + } + + return val; + } + + /*! + * \brief Splits a Bezier curve into two Bezier curves at a given parameter value * * \param [in] t parameter value between 0 and 1 at which to evaluate - * \param [out] c1, c2 Bezier curves that split the original + * \param [out] c1 First output Bezier curve + * \param [out] c2 Second output Bezier curve + * + * \pre Parameter \a t must be between 0 and 1 */ void split(T t, BezierCurve& c1, BezierCurve& c2) const { @@ -265,12 +321,7 @@ class BezierCurve const int end = ord - p; for(int k = 0; k <= end; ++k) { - PointType& pt1 = c2[k]; - const PointType& pt2 = c2[k + 1]; - for(int i = 0; i < NDIMS; ++i) - { - pt1[i] = (1 - t) * pt1[i] + t * pt2[i]; - } + c2[k] = PointType::lerp(c2[k], c2[k + 1], t); } c1[p] = c2[0]; } @@ -331,7 +382,7 @@ class BezierCurve }; //------------------------------------------------------------------------------ -/// Free functions implementing BezierCurve's operators +/// Free functions related to BezierCurve //------------------------------------------------------------------------------ template std::ostream& operator<<(std::ostream& os, const BezierCurve& bCurve) diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp new file mode 100644 index 0000000000..5805572c40 --- /dev/null +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -0,0 +1,254 @@ +// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file CurvedPolygon.hpp + * + * \brief A polygon primitive whose edges are Bezier curves + */ + +#ifndef AXOM_PRIMAL_CURVEDPOLYGON_HPP_ +#define AXOM_PRIMAL_CURVEDPOLYGON_HPP_ + +#include "axom/slic.hpp" + +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/Vector.hpp" +#include "axom/primal/geometry/NumericArray.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/BoundingBox.hpp" + +#include +#include + +namespace axom +{ +namespace primal +{ +// Forward declare the templated classes and operator functions +template +class CurvedPolygon; + +/*! \brief Overloaded output operator for polygons */ +template +std::ostream& operator<<(std::ostream& os, const CurvedPolygon& poly); + +/*! + * \class CurvedPolygon + * + * \brief Represents a polygon with curved edges defined by BezierCurves + * \tparam T the coordinate type, e.g., double, float, etc. + * \tparam NDIMS the number of dimensions + * \note The component curves should be ordered in a counter clockwise + * orientation with respect to the polygon's normal vector + */ +template +class CurvedPolygon +{ +public: + using PointType = Point; + using VectorType = Vector; + using NumArrayType = NumericArray; + using BezierCurveType = BezierCurve; + using BoundingBoxType = typename BezierCurveType::BoundingBoxType; + +public: + /// Default constructor for an empty polygon + CurvedPolygon() = default; + + /*! + * \brief Constructor for an empty CurvedPolygon that reserves space for + * the given number of Edges + * + * \param [in] numExpectedEdges number of edges for which to reserve space + * \pre numExpectedEdges is at least 1 + * + */ + explicit CurvedPolygon(int nEdges) + { + SLIC_ASSERT(nEdges >= 1); + m_edges.reserve(nEdges); + m_edges.resize(nEdges); + } + + /// Constructor from an array of \a nEdges curves + CurvedPolygon(BezierCurveType* curves, int nEdges) + { + SLIC_ASSERT(curves != nullptr); + SLIC_ASSERT(nEdges >= 1); + + m_edges.reserve(nEdges); + + for(int e = 0; e < nEdges; ++e) + { + this->addEdge(curves[e]); + } + } + + /// Clears the list of edges + void clear() { m_edges.clear(); } + + /// Returns true if the polygon has no edges + bool empty() const { return m_edges.empty(); } + + /// \name Operations on edges + /// @{ + + /// Return the number of edges in the polygon + int numEdges() const { return m_edges.size(); } + + void setNumEdges(int ngon) + { + SLIC_ASSERT(ngon >= 0); + m_edges.resize(ngon); + } + + /// Appends a BezierCurve to the list of edges + void addEdge(const BezierCurveType& c1) { m_edges.push_back(c1); } + + /// Splits an edge "in place" + void splitEdge(int idx, T t) + { + SLIC_ASSERT(idx < static_cast(m_edges.size())); + + m_edges.insert(m_edges.begin() + idx + 1, 1, m_edges[idx]); + auto& csplit = m_edges[idx]; + csplit.split(t, m_edges[idx], m_edges[idx + 1]); + } + + std::vector> getEdges() const { return m_edges; } + + /// @} + + /*! Retrieves the Bezier Curve at index idx */ + BezierCurveType& operator[](int idx) { return m_edges[idx]; } + /*! Retrieves the vertex at index idx */ + const BezierCurveType& operator[](int idx) const { return m_edges[idx]; } + + /// Tests equality of two CurvedPolygons + friend inline bool operator==(const CurvedPolygon& lhs, + const CurvedPolygon& rhs) + { + return lhs.m_edges == rhs.m_edges; + } + + /// Tests inequality of two CurvedPolygons + friend inline bool operator!=(const CurvedPolygon& lhs, + const CurvedPolygon& rhs) + { + return !(lhs == rhs); + } + + /*! + * \brief Simple formatted print of a CurvedPolygon instance + * + * \param os The output stream to write to + * \return A reference to the modified ostream + */ + std::ostream& print(std::ostream& os) const + { + const int sz = numEdges(); + + os << "{" << sz << "-sided Bezier polygon:"; + for(int i = 0; i < sz - 1; ++i) + { + os << m_edges[i] << ","; + } + if(sz >= 2) + { + os << m_edges[sz - 1]; + } + os << "}"; + + return os; + } + + /*! + * \brief Check closedness of a CurvedPolygon + * + * A CurvedPolygon is closed when the endpoint of each edge coincides with startpoint of next edge + * \return \a true, if the polygon is closed, \a false otherwise + */ + bool isClosed(double tol = 1e-5) const + { + using axom::utilities::isNearlyEqual; + + const double sq_tol = tol * tol; + const int nEdges = numEdges(); + + // initial basic check: no edges, or one edge or linear or quadratic order cannot be closed + if(nEdges < 1 || (nEdges == 1 && m_edges[0].getOrder() <= 2)) + { + return false; + } + + // foreach edge: check last vertex of previous edge against first vertex of current edge + for(int cur = 0, prev = nEdges - 1; cur < nEdges; prev = cur++) + { + const auto ord = m_edges[prev].getOrder(); + const auto& lastPrev = m_edges[prev][ord]; + const auto& firstCur = m_edges[cur][0]; + if(!isNearlyEqual(squared_distance(lastPrev, firstCur), 0., sq_tol)) + { + return false; + } + } + + return true; + } + + /// \brief Reverses orientation of a CurvedPolygon + void reverseOrientation() + { + const int ngon = numEdges(); + const int mid = ngon >> 1; + const bool isOdd = (ngon & 1) != 0; + + // Swap matching left/right cases, unmatched center is dealt with below + for(int i = 0; i < mid; ++i) + { + const int left = i; + const int right = ngon - i - 1; + m_edges[left].reverseOrientation(); + m_edges[right].reverseOrientation(); + axom::utilities::swap(m_edges[left], m_edges[right]); + } + + // Handle unmatched center curve, if necessary + if(isOdd) + { + m_edges[mid].reverseOrientation(); + } + } + + /// Returns an axis-aligned bounding box containing the CurvedPolygon + BoundingBoxType boundingBox() const + { + BoundingBoxType bbox; + for(const auto& cp : m_edges) + { + bbox.addBox(cp.boundingBox()); + } + return bbox; + } + +private: + std::vector> m_edges; +}; + +//------------------------------------------------------------------------------ +/// Free functions implementing Polygon's operators +//------------------------------------------------------------------------------ +template +std::ostream& operator<<(std::ostream& os, const CurvedPolygon& poly) +{ + poly.print(os); + return os; +} + +} // namespace primal +} // namespace axom + +#endif // AXOM_PRIMAL_CURVEDPOLYGON_HPP_ diff --git a/src/axom/primal/geometry/Polygon.hpp b/src/axom/primal/geometry/Polygon.hpp index 0b47173729..b90c1822a1 100644 --- a/src/axom/primal/geometry/Polygon.hpp +++ b/src/axom/primal/geometry/Polygon.hpp @@ -12,12 +12,10 @@ #ifndef AXOM_PRIMAL_POLYGON_HPP_ #define AXOM_PRIMAL_POLYGON_HPP_ +#include "axom/core/Array.hpp" #include "axom/primal/geometry/Point.hpp" -#include "axom/primal/geometry/Vector.hpp" -#include "axom/primal/geometry/NumericArray.hpp" -#include -#include // for std::ostream +#include namespace axom { @@ -45,11 +43,6 @@ class Polygon { public: using PointType = Point; - using VectorType = Vector; - using NumArrayType = NumericArray; - -private: - using Coords = std::vector; public: /*! Default constructor for an empty polygon */ @@ -84,26 +77,25 @@ class Polygon const PointType& operator[](int idx) const { return m_vertices[idx]; } /*! - * \brief Computes the centroid as the average of the polygon's vertex - * positions - * - * \return The centroid of the polygon's vertices + * \brief Computes the average of the polygon's vertex positions * + * \return A point at the mean of the polygon's vertices * \pre polygon.isValid() is true */ - PointType centroid() const + PointType vertexMean() const { SLIC_ASSERT(isValid()); - NumArrayType sum; + PointType sum; - for(int i = 0; i < numVertices(); ++i) + const int sz = numVertices(); + for(int i = 0; i < sz; ++i) { - sum += m_vertices[i].array(); + sum.array() += m_vertices[i].array(); } - sum /= numVertices(); + sum.array() /= sz; - return PointType(sum); + return sum; } /*! @@ -139,7 +131,7 @@ class Polygon bool isValid() const { return m_vertices.size() >= 3; } private: - Coords m_vertices; + axom::Array m_vertices; }; //------------------------------------------------------------------------------ diff --git a/src/axom/primal/operators/compute_moments.hpp b/src/axom/primal/operators/compute_moments.hpp new file mode 100644 index 0000000000..88f258737e --- /dev/null +++ b/src/axom/primal/operators/compute_moments.hpp @@ -0,0 +1,146 @@ +// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#ifndef AXOM_PRIMAL_COMPUTE_MOMENTS_HPP_ +#define AXOM_PRIMAL_COMPUTE_MOMENTS_HPP_ + +/*! + * \file compute_moments.hpp + * + * \brief Consists of a set of methods to compute areas/volumes and centroids + * for Polygons and CurvedPolygons composed of BezierCurves + */ + +#include "axom/config.hpp" +#include "axom/core.hpp" +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/CurvedPolygon.hpp" +#include "axom/primal/operators/detail/compute_moments_impl.hpp" + +#include +#include + +namespace axom +{ +namespace primal +{ +/*! + * \brief Calculates the sector area of a planar Bezier Curve + * + * The sector area is the area between the curve and the origin. + * The equation and derivation is described in: + * Ueda, K. "Signed area of sectors between spline curves and the origin" + * IEEE International Conference on Information Visualization, 1999. + */ +template +T sector_area(const primal::BezierCurve& curve) +{ + // Weights for each polynomial order are precomputed and memoized + static detail::MemoizedSectorAreaWeights s_weights; + + T A = 0; + const int ord = curve.getOrder(); + const auto& weights = s_weights.getWeights(ord); + + for(int p = 0; p <= ord; ++p) + { + for(int q = 0; q <= ord; ++q) + { + A += weights(p, q) * curve[p][1] * curve[q][0]; + } + } + return A; +} + +/*! + * \brief Calculates the sector centroid of a planar Bezier Curve + * + * This is the centroid of the region between the curve and the origin. + * The equation and derivation are generalizations of: + * Ueda, K. "Signed area of sectors between spline curves and the origin" + * IEEE International Conference on Information Visualization, 1999. + */ +template +primal::Point sector_centroid(const primal::BezierCurve& curve) +{ + // Weights for each polynomial order's centroid are precomputed and memoized + static detail::MemoizedSectorCentroidWeights s_weights; + + T Mx = 0; + T My = 0; + const int ord = curve.getOrder(); + for(int r = 0; r <= ord; ++r) + { + const auto& weights_r = s_weights.getWeights(ord, r); + for(int p = 0; p <= ord; ++p) + { + for(int q = 0; q <= ord; ++q) + { + Mx += weights_r(p, q) * curve[p][1] * curve[q][0] * curve[r][0]; + My += weights_r(p, q) * curve[p][1] * curve[q][0] * curve[r][1]; + } + } + } + return primal::Point {Mx, My}; +} + +/// \brief Returns the area enclosed by the CurvedPolygon +template +T area(const primal::CurvedPolygon& poly, double tol = 1e-8) +{ + const int ngon = poly.numEdges(); + T A = 0.0; + if(!poly.isClosed(1e3 * tol)) + { + SLIC_DEBUG( + "Warning! The area is 0 because the curved polygon is not closed."); + return A; + } + else + { + for(int ed = 0; ed < ngon; ++ed) + { + A += primal::sector_area(poly[ed]); + } + return A; + } +} + +/// \brief Returns the centroid of the CurvedPolygon +template +primal::Point centroid(const primal::CurvedPolygon& poly, + double tol = 1e-8) +{ + using PointType = primal::Point; + + const int ngon = poly.numEdges(); + PointType M; + + if(!poly.isClosed(1e3 * tol)) + { + SLIC_DEBUG( + "Warning! The moments are 0 because the curved polygon is not closed."); + return M; + } + else + { + const T A = area(poly, tol); + if(A != 0.) + { + for(int ed = 0; ed < ngon; ++ed) + { + M.array() += primal::sector_centroid(poly[ed]).array(); + } + M.array() /= A; + } + return M; + } +} + +} // namespace primal +} // namespace axom + +#endif // AXOM_PRIMAL_COMPUTE_MOMENTS_HPP_ diff --git a/src/axom/primal/operators/detail/compute_moments_impl.hpp b/src/axom/primal/operators/detail/compute_moments_impl.hpp new file mode 100644 index 0000000000..b3df6060cb --- /dev/null +++ b/src/axom/primal/operators/detail/compute_moments_impl.hpp @@ -0,0 +1,204 @@ +// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#ifndef AXOM_PRIMAL_COMPUTE_MOMENTS_IMPL_HPP_ +#define AXOM_PRIMAL_COMPUTE_MOMENTS_IMPL_HPP_ + +/*! + * \file compute_moments_impl.hpp + * + * \brief Consists of implementation helpers for computing areas/volumes and centroids + * for Polygons and CurvedPolygons composed of BezierCurves + */ + +#include "axom/config.hpp" +#include "axom/core.hpp" +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/CurvedPolygon.hpp" + +#include +#include + +namespace axom +{ +namespace primal +{ +namespace detail +{ +/// Utility class that caches precomputed coefficient matrices for sector_area computation +template +class MemoizedSectorAreaWeights +{ +public: + using SectorWeights = numerics::Matrix; + + MemoizedSectorAreaWeights() = default; + + ~MemoizedSectorAreaWeights() + { + for(auto& p : m_sectorWeightsMap) + { + delete[] p.second->data(); // delete the matrix's data + delete p.second; // delete the matrix + p.second = nullptr; + } + m_sectorWeightsMap.clear(); + } + + /// Returns a memoized matrix of coeficients for sector area computation + const SectorWeights& getWeights(int order) const + { + // Compute and cache the weights if they are not already available + if(m_sectorWeightsMap.find(order) == m_sectorWeightsMap.end()) + { + SectorWeights* weights = generateBezierCurveSectorWeights(order); + m_sectorWeightsMap[order] = weights; + } + + return *(m_sectorWeightsMap[order]); + } + +private: + /*! + * \brief Computes the weights for BezierCurve's sector_area() function + * + * \param order The polynomial order of the curve + * \return An anti-symmetric matrix with (order+1)*{order+1) entries + * containing the integration weights for entry (i,j) + * + * The derivation is provided in: + * Ueda, K. "Signed area of sectors between spline curves and the origin" + * IEEE International Conference on Information Visualization, 1999. + */ + SectorWeights* generateBezierCurveSectorWeights(int ord) const + { + const bool memoryIsExternal = true; + const int SZ = ord + 1; + SectorWeights* weights = + new SectorWeights(SZ, SZ, new T[SZ * SZ], memoryIsExternal); + + T binom_2n_n = static_cast(utilities::binomialCoefficient(2 * ord, ord)); + for(int i = 0; i <= ord; ++i) + { + (*weights)(i, i) = 0.; // zero on the diagonal + for(int j = i + 1; j <= ord; ++j) + { + double val = 0.; + if(i != j) + { + T binom_ij_i = static_cast(utilities::binomialCoefficient(i + j, i)); + T binom_2nij_nj = static_cast( + utilities::binomialCoefficient(2 * ord - i - j, ord - j)); + + val = ((j - i) * ord) / binom_2n_n * + (binom_ij_i / static_cast(i + j)) * + (binom_2nij_nj / (2. * ord - j - i)); + } + (*weights)(i, j) = val; // antisymmetric + (*weights)(j, i) = -val; + } + } + return weights; + } + +private: + mutable std::map m_sectorWeightsMap; +}; + +/// Utility class that caches precomputed coefficient matrices for sector_centroid() computation +template +class MemoizedSectorCentroidWeights +{ +public: + using SectorWeights = numerics::Matrix; + + MemoizedSectorCentroidWeights() = default; + + ~MemoizedSectorCentroidWeights() + { + for(auto& p : m_sectorWeightsMap) // for each matrix of weights + { + delete[] p.second->data(); // delete the matrix's data + delete p.second; // delete the matrix + p.second = nullptr; + } + m_sectorWeightsMap.clear(); + } + + /// Returns a memoized matrix of sector moment coeficients for component \a dim or order \a order + const SectorWeights& getWeights(int order, int dim) const + { + // Compute and cache the weights if they are not already available + if(m_sectorWeightsMap.find(std::make_pair(order, dim)) == + m_sectorWeightsMap.end()) + { + auto vec = generateBezierCurveSectorCentroidWeights(order); + for(int d = 0; d <= order; ++d) + { + m_sectorWeightsMap[std::make_pair(order, d)] = vec[d]; + } + } + + return *(m_sectorWeightsMap[std::make_pair(order, dim)]); + } + + /*! + * \brief Computes the weights for BezierCurve's sectorCentroid() function + * + * \param order The polynomial order of the curve + * \return An anti-symmetric matrix with (order+1)*{order+1) entries + * containing the integration weights for entry (i,j) + * + * The derivation is provided in: + * Ueda, K. "Signed area of sectors between spline curves and the origin" + * IEEE International Conference on Information Visualization, 1999. + */ + std::vector generateBezierCurveSectorCentroidWeights(int ord) const + { + const bool memoryIsExternal = true; + const int SZ = ord + 1; + + std::vector weights; + weights.resize(SZ); + for(int k = 0; k <= ord; ++k) + { + SectorWeights* weights_k = + new SectorWeights(SZ, SZ, new T[SZ * SZ], memoryIsExternal); + for(int i = 0; i <= ord; ++i) + { + (*weights_k)(i, i) = 0.; // zero on the diagonal + for(int j = i + 1; j <= ord; ++j) + { + double val = 0.; + if(i != j) + { + T binom_n_i = static_cast(utilities::binomialCoefficient(ord, i)); + T binom_n_j = static_cast(utilities::binomialCoefficient(ord, j)); + T binom_n_k = static_cast(utilities::binomialCoefficient(ord, k)); + T binom_3n2_ijk1 = static_cast( + utilities::binomialCoefficient(3 * ord - 2, i + j + k - 1)); + + val = (1. * (j - i)) / (3. * (3 * ord - 1)) * + (1. * binom_n_i * binom_n_j * binom_n_k / (1. * binom_3n2_ijk1)); + } + (*weights_k)(i, j) = val; // antisymmetric + (*weights_k)(j, i) = -val; + } + } + weights[k] = weights_k; + } + return weights; + } + +private: + mutable std::map, SectorWeights*> m_sectorWeightsMap; +}; + +} // namespace detail +} // namespace primal +} // namespace axom + +#endif // AXOM_PRIMAL_COMPUTE_MOMENTS_IMPL_HPP_ \ No newline at end of file diff --git a/src/axom/primal/operators/detail/intersect_bezier_impl.hpp b/src/axom/primal/operators/detail/intersect_bezier_impl.hpp index a89bbb96f4..66925002a8 100644 --- a/src/axom/primal/operators/detail/intersect_bezier_impl.hpp +++ b/src/axom/primal/operators/detail/intersect_bezier_impl.hpp @@ -57,9 +57,9 @@ namespace detail * \return True if the two curves intersect, False otherwise * \sa intersect_bezier */ -template -bool intersect_bezier_curves(const BezierCurve &c1, - const BezierCurve &c2, +template +bool intersect_bezier_curves(const BezierCurve &c1, + const BezierCurve &c2, std::vector &sp, std::vector &tp, double sq_tol, @@ -97,19 +97,19 @@ bool intersect_bezier_curves(const BezierCurve &c1, * \note This function does not properly handle collinear lines */ -template -bool intersect_2d_linear(const Point &a, - const Point &b, - const Point &c, - const Point &d, +template +bool intersect_2d_linear(const Point &a, + const Point &b, + const Point &c, + const Point &d, T &s, T &t); //------------------------------ IMPLEMENTATIONS ------------------------------ -template -bool intersect_bezier_curves(const BezierCurve &c1, - const BezierCurve &c2, +template +bool intersect_bezier_curves(const BezierCurve &c1, + const BezierCurve &c2, std::vector &sp, std::vector &tp, double sq_tol, @@ -120,8 +120,7 @@ bool intersect_bezier_curves(const BezierCurve &c1, double t_offset, double t_scale) { - using BCurve = BezierCurve; - SLIC_ASSERT(NDIMS == 2); + using BCurve = BezierCurve; // Check bounding boxes to short-circuit the intersection if(!intersect(c1.boundingBox(), c2.boundingBox())) @@ -152,6 +151,7 @@ bool intersect_bezier_curves(const BezierCurve &c1, s_scale *= scaleFac; + // Note: we want to find all intersections, so don't short-circuit if(intersect_bezier_curves(c2, c3, tp, @@ -185,11 +185,11 @@ bool intersect_bezier_curves(const BezierCurve &c1, return foundIntersection; } -template -bool intersect_2d_linear(const Point &a, - const Point &b, - const Point &c, - const Point &d, +template +bool intersect_2d_linear(const Point &a, + const Point &b, + const Point &c, + const Point &d, T &s, T &t) { @@ -199,18 +199,16 @@ bool intersect_2d_linear(const Point &a, // Note: Uses exact floating point comparisons since the subdivision algorithm // provides both sides of the line segments for interior curve points. - AXOM_STATIC_ASSERT(NDIMS == 2); - // compute signed areas of endpoints of segment (c,d) w.r.t. segment (a,b) - auto area1 = twoDcross(a, b, c); - auto area2 = twoDcross(a, b, d); + const auto area1 = twoDcross(a, b, c); + const auto area2 = twoDcross(a, b, d); // early return if both have same orientation, or if d is collinear w/ (a,b) if(area2 == 0. || (area1 * area2) > 0.) return false; // compute signed areas of endpoints of segment (a,b) w.r.t. segment (c,d) - auto area3 = twoDcross(c, d, a); - auto area4 = area3 + area1 - area2; // equivalent to twoDcross(c,d,b) + const auto area3 = twoDcross(c, d, a); + const auto area4 = area3 + area1 - area2; // equivalent to twoDcross(c,d,b) // early return if both have same orientation, or if b is collinear w/ (c,d) if(area4 == 0. || (area3 * area4) > 0.) return false; diff --git a/src/axom/primal/operators/intersect.hpp b/src/axom/primal/operators/intersect.hpp index e6e7baac13..f20eca17d0 100644 --- a/src/axom/primal/operators/intersect.hpp +++ b/src/axom/primal/operators/intersect.hpp @@ -12,6 +12,7 @@ #ifndef AXOM_PRIMAL_INTERSECT_HPP_ #define AXOM_PRIMAL_INTERSECT_HPP_ +#include "axom/config.hpp" #include "axom/core/Macros.hpp" #include "axom/core/utilities/Utilities.hpp" @@ -505,9 +506,9 @@ bool intersect(const OrientedBoundingBox& b1, * contain their first endpoint, but not their last endpoint. Thus, the * curves do not intersect at \f$ s==1 \f$ or at \f$ t==1 \f$. */ -template -bool intersect(const BezierCurve& c1, - const BezierCurve& c2, +template +bool intersect(const BezierCurve& c1, + const BezierCurve& c2, std::vector& sp, std::vector& tp, double tol = 1E-8) diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index b8d8667634..6d246266bb 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -13,6 +13,8 @@ set( primal_tests primal_clip.cpp primal_closest_point.cpp primal_compute_bounding_box.cpp + primal_compute_moments.cpp + primal_curved_polygon.cpp primal_in_sphere.cpp primal_intersect.cpp primal_intersect_impl.cpp diff --git a/src/axom/primal/tests/primal_bezier_curve.cpp b/src/axom/primal/tests/primal_bezier_curve.cpp index 9f95645ec5..e355424043 100644 --- a/src/axom/primal/tests/primal_bezier_curve.cpp +++ b/src/axom/primal/tests/primal_bezier_curve.cpp @@ -3,12 +3,15 @@ // // SPDX-License-Identifier: (BSD-3-Clause) -/* /file primal_bezier_curve.cpp - * /brief This file tests primal's Bezier curve functionality +/*! + * \file primal_bezier_curve.cpp + * \brief This file tests primal's Bezier curve functionality */ #include "gtest/gtest.h" +#include "axom/slic.hpp" + #include "axom/primal/geometry/BezierCurve.hpp" #include "axom/primal/operators/squared_distance.hpp" @@ -56,8 +59,8 @@ TEST(primal_beziercurve, set_order) EXPECT_EQ(-1, bCurve.getOrder()); const int order = 1; - PointType controlPoints[] = {PointType::make_point(0.6, 1.2, 1.0), - PointType::make_point(0.0, 1.6, 1.8)}; + PointType controlPoints[2] = {PointType {0.6, 1.2, 1.0}, + PointType {0.0, 1.6, 1.8}}; bCurve.setOrder(order); EXPECT_EQ(order, bCurve.getOrder()); @@ -85,8 +88,8 @@ TEST(primal_beziercurve, point_array_constructor) using PointType = primal::Point; using BezierCurveType = primal::BezierCurve; - PointType controlPoints[2] = {PointType::make_point(0.6, 1.2, 1.0), - PointType::make_point(0.0, 1.6, 1.8)}; + PointType controlPoints[2] = {PointType {0.6, 1.2, 1.0}, + PointType {0.0, 1.6, 1.8}}; BezierCurveType bCurve(controlPoints, 1); @@ -140,14 +143,14 @@ TEST(primal_beziercurve, evaluate) using BezierCurveType = primal::BezierCurve; const int order = 3; - PointType data[order + 1] = {PointType::make_point(0.6, 1.2, 1.0), - PointType::make_point(1.3, 1.6, 1.8), - PointType::make_point(2.9, 2.4, 2.3), - PointType::make_point(3.2, 3.5, 3.0)}; + PointType data[order + 1] = {PointType {0.6, 1.2, 1.0}, + PointType {1.3, 1.6, 1.8}, + PointType {2.9, 2.4, 2.3}, + PointType {3.2, 3.5, 3.0}}; BezierCurveType b2Curve(data, order); - PointType midtval = PointType::make_point(2.05, 2.0875, 2.0375); + PointType midtval {2.05, 2.0875, 2.0375}; // Evaluate the curve at several parameter values // Curve should interpolate endpoints @@ -163,6 +166,43 @@ TEST(primal_beziercurve, evaluate) } } +//------------------------------------------------------------------------------ +TEST(primal_beziercurve_, tangent) +{ + SLIC_INFO("Testing Bezier tangent calculation"); + + const int DIM = 3; + using CoordType = double; + using PointType = primal::Point; + using VectorType = primal::Vector; + using BezierCurveType = primal::BezierCurve; + + const int order = 3; + PointType data[order + 1] = {PointType {0.6, 1.2, 1.0}, + PointType {1.3, 1.6, 1.8}, + PointType {2.9, 2.4, 2.3}, + PointType {3.2, 3.5, 3.0}}; + + BezierCurveType b2Curve(data, order); + + VectorType midtval = VectorType {3.15, 2.325, 1.875}; + VectorType starttval = VectorType {2.1, 1.2, 2.4}; + VectorType endtval = VectorType {.9, 3.3, 2.1}; + + // Evaluate the curve at several parameter values + // Curve should be tangent to control net at endpoints + VectorType eval0 = b2Curve.dt(0.0); + VectorType eval1 = b2Curve.dt(1.0); + VectorType evalMid = b2Curve.dt(0.5); + + for(int i = 0; i < DIM; ++i) + { + EXPECT_NEAR(starttval[i], eval0[i], 1e-15); + EXPECT_NEAR(endtval[i], eval1[i], 1e-15); + EXPECT_NEAR(midtval[i], evalMid[i], 1e-15); + } +} + //------------------------------------------------------------------------------ TEST(primal_beziercurve, split_cubic) { @@ -174,10 +214,10 @@ TEST(primal_beziercurve, split_cubic) using BezierCurveType = primal::BezierCurve; const int order = 3; - PointType data[order + 1] = {PointType::make_point(0.6, 1.2, 1.0), - PointType::make_point(1.3, 1.6, 1.8), - PointType::make_point(2.9, 2.4, 2.3), - PointType::make_point(3.2, 3.5, 3.0)}; + PointType data[order + 1] = {PointType {0.6, 1.2, 1.0}, + PointType {1.3, 1.6, 1.8}, + PointType {2.9, 2.4, 2.3}, + PointType {3.2, 3.5, 3.0}}; BezierCurveType b2Curve(data, order); BezierCurveType b3Curve(order); // Checks split with order constructor @@ -242,8 +282,7 @@ TEST(primal_beziercurve, split_linear) using BezierCurveType = primal::BezierCurve; const int order = 1; - PointType data[order + 1] = {PointType::make_point(-1, -5), - PointType::make_point(1, 5)}; + PointType data[order + 1] = {PointType {-1, -5}, PointType {1, 5}}; BezierCurveType b(data, order); { @@ -292,9 +331,9 @@ TEST(primal_beziercurve, split_quadratic) const int order = 2; // Control points for the three levels of the quadratic de Casteljau algorithm - PointType lev0[3] = {PointType::make_point(1.1, 1.1), - PointType::make_point(5.5, 5.5), - PointType::make_point(9.9, 2.2)}; + PointType lev0[3] = {PointType {1.1, 1.1}, + PointType {5.5, 5.5}, + PointType {9.9, 2.2}}; PointType lev1[2] = {PointType::lerp(lev0[0], lev0[1], t), PointType::lerp(lev0[1], lev0[2], t)}; @@ -318,9 +357,11 @@ TEST(primal_beziercurve, split_quadratic) BezierCurveType c1, c2; b.split(t, c1, c2); - SLIC_INFO("" - << "Original quadratic: " << b << "\nCurves after splitting at t = " - << t << "\n\t c1: " << c1 << "\n\t c2: " << c2); + SLIC_INFO("" // + << "Original quadratic: " << b // + << "\nCurves after splitting at t = " << t // + << "\n\t c1: " << c1 // + << "\n\t c2: " << c2); // Check values for(int p = 0; p <= order; ++p) @@ -360,8 +401,8 @@ TEST(primal_beziercurve, isLinear) auto curve = BezierCurveType(order); EXPECT_TRUE(curve.isLinear()); - curve[0] = PointType::make_point(1., 1.8); - curve[1] = PointType::make_point(-12., 3.5); + curve[0] = PointType {1., 1.8}; + curve[1] = PointType {-12., 3.5}; EXPECT_TRUE(curve.isLinear()); } @@ -372,14 +413,14 @@ TEST(primal_beziercurve, isLinear) EXPECT_TRUE(curve.isLinear()); // straight line - curve[0] = PointType::make_point(1, 1); - curve[1] = PointType::make_point(2, 2); - curve[2] = PointType::make_point(3, 3); + curve[0] = PointType {1, 1}; + curve[1] = PointType {2, 2}; + curve[2] = PointType {3, 3}; EXPECT_TRUE(curve.isLinear()); // move middle point and check linearity with different tolerances VectorType v(curve[2], curve[0]); - auto normal = VectorType::make_vector(-v[1], v[0]); + auto normal = VectorType {-v[1], v[0]}; curve[1].array() += 0.005 * normal.array(); SLIC_INFO("Updated curve: " << curve); @@ -402,6 +443,55 @@ TEST(primal_beziercurve, isLinear) } } +TEST(primal_beziercurve, reverseOrientation) +{ + SLIC_INFO("Testing reverseOrientation() on Bezier curves"); + + { + const int DIM = 2; + using CoordType = double; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + + // test different orders + for(int order = 0; order <= 10; ++order) + { + // control points for curve monotonically increase + axom::Array pts(order + 1); + for(int i = 0; i <= order; ++i) + { + pts[i] = PointType(i); + } + BezierCurveType curve(pts.data(), order); + + for(int i = 1; i <= order; ++i) + { + EXPECT_GT(curve[i][0], curve[i - 1][0]); + } + + // create a reversed curve and check that it monotonically decreases + BezierCurveType reversed = curve; + reversed.reverseOrientation(); + + for(int i = 1; i <= order; ++i) + { + EXPECT_LT(reversed[i][0], reversed[i - 1][0]); + } + + // Check that the control points are actually reversed + for(int i = 0; i <= order; ++i) + { + EXPECT_EQ(curve[i], reversed[order - i]); + } + + // check that reversing again reverts to the original + BezierCurveType reversedAgain = reversed; + reversedAgain.reverseOrientation(); + EXPECT_EQ(curve, reversedAgain); + } + } +} + //------------------------------------------------------------------------------ int main(int argc, char* argv[]) @@ -409,7 +499,8 @@ int main(int argc, char* argv[]) int result = 0; ::testing::InitGoogleTest(&argc, argv); - axom::slic::SimpleLogger logger; // create & initialize test logger, + + axom::slic::SimpleLogger logger; result = RUN_ALL_TESTS(); diff --git a/src/axom/primal/tests/primal_bezier_intersect.cpp b/src/axom/primal/tests/primal_bezier_intersect.cpp index 87cbea79dc..2988f01b7c 100644 --- a/src/axom/primal/tests/primal_bezier_intersect.cpp +++ b/src/axom/primal/tests/primal_bezier_intersect.cpp @@ -3,8 +3,9 @@ // // SPDX-License-Identifier: (BSD-3-Clause) -/* /file primal_bezier_intersect.cpp - * /brief This file tests the Bezier curve intersection routines +/*! + * \file primal_bezier_intersect.cpp + * \brief This file tests the Bezier curve intersection routines */ #include "gtest/gtest.h" @@ -31,15 +32,16 @@ namespace primal = axom::primal; * Param \a shouldPrintIntersections is used for debugging and for generating * the initial array of expected intersections. */ -template -void checkIntersections(const primal::BezierCurve& curve1, - const primal::BezierCurve& curve2, +template +void checkIntersections(const primal::BezierCurve& curve1, + const primal::BezierCurve& curve2, const std::vector& exp_s, const std::vector& exp_t, double eps, double test_eps, bool shouldPrintIntersections = false) { + constexpr int DIM = 2; using Array = std::vector; // Check validity of input data exp_s and exp_t. @@ -127,13 +129,11 @@ TEST(primal_bezier_inter, linear_bezier) { SCOPED_TRACE("linear bezier simple"); - PointType data1[order + 1] = {PointType::make_point(0.0, 0.0), - PointType::make_point(1.0, 1.0)}; + PointType data1[order + 1] = {PointType {0.0, 0.0}, PointType {1.0, 1.0}}; BezierCurveType curve1(data1, order); - PointType data2[order + 1] = {PointType::make_point(0.0, 1.0), - PointType::make_point(1.0, 0.0)}; + PointType data2[order + 1] = {PointType {0.0, 1.0}, PointType {1.0, 0.0}}; BezierCurveType curve2(data2, order); std::vector exp_intersections1 = {0.5}; @@ -152,13 +152,11 @@ TEST(primal_bezier_inter, linear_bezier) { SCOPED_TRACE("linear bezier endpoints"); - PointType data1[order + 1] = {PointType::make_point(1.0, 1.0), - PointType::make_point(0.0, 0.0)}; + PointType data1[order + 1] = {PointType {1.0, 1.0}, PointType {0.0, 0.0}}; BezierCurveType curve1(data1, order); - PointType data2[order + 1] = {PointType::make_point(1.0, 1.0), - PointType::make_point(2.0, 0.0)}; + PointType data2[order + 1] = {PointType {1.0, 1.0}, PointType {2.0, 0.0}}; BezierCurveType curve2(data2, order); std::vector exp_intersections1 = {0.0}; @@ -177,13 +175,11 @@ TEST(primal_bezier_inter, linear_bezier) { SCOPED_TRACE("linear bezier non-midpoint"); - PointType data1[order + 1] = {PointType::make_point(0.0, 0.0), - PointType::make_point(4.0, 2.0)}; + PointType data1[order + 1] = {PointType {0.0, 0.0}, PointType {4.0, 2.0}}; BezierCurveType curve1(data1, order); - PointType data2[order + 1] = {PointType::make_point(-2.0, 2.0), - PointType::make_point(2.0, 0.0)}; + PointType data2[order + 1] = {PointType {-2.0, 2.0}, PointType {2.0, 0.0}}; BezierCurveType curve2(data2, order); std::vector exp_intersections1 = {.25}; @@ -201,7 +197,7 @@ TEST(primal_bezier_inter, linear_bezier) TEST(primal_bezier_inter, linear_bezier_interp_params) { - static const int DIM = 2; + constexpr int DIM = 2; using CoordType = double; using PointType = primal::Point; @@ -225,12 +221,10 @@ TEST(primal_bezier_inter, linear_bezier_interp_params) sstr << "linear bezier perpendicular (s,t) = (" << s << "," << t << ")"; SCOPED_TRACE(sstr.str()); - PointType data1[order + 1] = {PointType::make_point(0.0, s), - PointType::make_point(1.0, s)}; + PointType data1[order + 1] = {PointType {0.0, s}, PointType {1.0, s}}; BezierCurveType curve1(data1, order); - PointType data2[order + 1] = {PointType::make_point(t, 0.0), - PointType::make_point(t, 1.0)}; + PointType data2[order + 1] = {PointType {t, 0.0}, PointType {t, 1.0}}; BezierCurveType curve2(data2, order); std::vector exp_intersections1 = {t}; @@ -269,18 +263,18 @@ TEST(primal_bezier_inter, no_intersections_bezier) const int order = 3; // cubic line - PointType data1[order + 1] = {PointType::make_point(0.0, 0.0), - PointType::make_point(1.0, 0.0), - PointType::make_point(2.0, 0.0), - PointType::make_point(3.0, 0.0)}; + PointType data1[order + 1] = {PointType {0.0, 0.0}, + PointType {1.0, 0.0}, + PointType {2.0, 0.0}, + PointType {3.0, 0.0}}; BezierCurveType curve1(data1, order); // Cubic curve - PointType data2[order + 1] = {PointType::make_point(0.0, 0.5), - PointType::make_point(1.0, 1.0), - PointType::make_point(2.0, 3.0), - PointType::make_point(3.0, 1.5)}; + PointType data2[order + 1] = {PointType {0.0, 0.5}, + PointType {1.0, 1.0}, + PointType {2.0, 3.0}, + PointType {3.0, 1.5}}; BezierCurveType curve2(data2, order); std::vector exp_intersections; @@ -316,10 +310,10 @@ TEST(primal_bezier_inter, cubic_quadratic_bezier) curve1[0] = data1; // Cubic curve - PointType data2[order2 + 1] = {PointType::make_point(0.0, 0.5), - PointType::make_point(1.0, -1.0), - PointType::make_point(2.0, 1.0), - PointType::make_point(3.0, -0.5)}; + PointType data2[order2 + 1] = {PointType {0.0, 0.5}, + PointType {1.0, -1.0}, + PointType {2.0, 1.0}, + PointType {3.0, -0.5}}; BezierCurveType curve2(data2, order2); // Note: same intersection params for curve and line @@ -337,7 +331,7 @@ TEST(primal_bezier_inter, cubic_quadratic_bezier) { curve1[i][0] = curve1[i][0] * (otherorder - 1) / (1.0 * otherorder); } - curve1[otherorder] = PointType::make_point(3.0, 0); + curve1[otherorder] = PointType {3.0, 0}; SLIC_INFO("Testing w/ order 3 and " << otherorder); std::stringstream sstr; @@ -366,18 +360,18 @@ TEST(primal_bezier_inter, cubic_bezier_varying_eps) const int order = 3; // cubic line - PointType data1[order + 1] = {PointType::make_point(0.0, 0.0), - PointType::make_point(1.0, 0.0), - PointType::make_point(2.0, 0.0), - PointType::make_point(3.0, 0.0)}; + PointType data1[order + 1] = {PointType {0.0, 0.0}, + PointType {1.0, 0.0}, + PointType {2.0, 0.0}, + PointType {3.0, 0.0}}; BezierCurveType curve1(data1, order); // Cubic curve - PointType data2[order + 1] = {PointType::make_point(0.0, 0.5), - PointType::make_point(1.0, -1.0), - PointType::make_point(2.0, 1.0), - PointType::make_point(3.0, -0.5)}; + PointType data2[order + 1] = {PointType {0.0, 0.5}, + PointType {1.0, -1.0}, + PointType {2.0, 1.0}, + PointType {3.0, -0.5}}; BezierCurveType curve2(data2, order); // Note: same intersection params for curve and line @@ -415,17 +409,17 @@ TEST(primal_bezier_inter, cubic_bezier_nine_intersections) // A configuration of two cubic bezier curves that intersect at nine points const int order = 3; - PointType data1[order + 1] = {PointType::make_point(100, 90), - PointType::make_point(125, 260), - PointType::make_point(125, 0), - PointType::make_point(140, 145)}; + PointType data1[order + 1] = {PointType {100, 90}, + PointType {125, 260}, + PointType {125, 0}, + PointType {140, 145}}; BezierCurveType curve1(data1, order); - PointType data2[order + 1] = {PointType::make_point(75, 110), - PointType::make_point(265, 120), - PointType::make_point(0, 130), - PointType::make_point(145, 135)}; + PointType data2[order + 1] = {PointType {75, 110}, + PointType {265, 120}, + PointType {0, 130}, + PointType {145, 135}}; BezierCurveType curve2(data2, order); const double eps = 1E-16; diff --git a/src/axom/primal/tests/primal_clip.cpp b/src/axom/primal/tests/primal_clip.cpp index 578f206cb9..7b87081dbd 100644 --- a/src/axom/primal/tests/primal_clip.cpp +++ b/src/axom/primal/tests/primal_clip.cpp @@ -43,16 +43,18 @@ TEST(primal_clip, simple_clip) bbox.addPoint(PointType::ones()); PointType points[] = { - PointType::make_point(2, 2, 2), - PointType::make_point(2, 2, 4), - PointType::make_point(2, 4, 2), - PointType::make_point(-100, -100, 0.5), - PointType::make_point(-100, 100, 0.5), - PointType::make_point(100, 0, 0.5), - PointType::make_point(0.25, 0.25, 0.5), - PointType::make_point(0.75, 0.25, 0.5), - PointType::make_point(0.66, 0.5, 0.5), - PointType::make_point(1.5, 0.5, 0.5), + PointType {2, 2, 2}, + PointType {2, 2, 4}, + PointType {2, 4, 2}, + + PointType {-100, -100, 0.5}, + PointType {-100, 100, 0.5}, + PointType {100, 0, 0.5}, + + PointType {0.25, 0.25, 0.5}, + PointType {0.75, 0.25, 0.5}, + PointType {0.66, 0.5, 0.5}, + PointType {1.5, 0.5, 0.5}, }; { @@ -78,7 +80,7 @@ TEST(primal_clip, simple_clip) PolygonType poly = axom::primal::clip(tri, bbox); EXPECT_EQ(4, poly.numVertices()); - EXPECT_EQ(PointType(.5), poly.centroid()); + EXPECT_EQ(PointType(.5), poly.vertexMean()); SLIC_INFO("Intersection of triangle " << tri << " and bounding box " << bbox << " is polygon" << poly); @@ -101,12 +103,12 @@ TEST(primal_clip, unit_simplex) double delta = 1e-5; // Test the "unit simplex", and a jittered version - PointType points[] = {PointType::make_point(1, 0, 0), - PointType::make_point(0, 1, 0), - PointType::make_point(0, 0, 1), - PointType::make_point(1 + delta, delta, delta), - PointType::make_point(delta, 1 + delta, delta), - PointType::make_point(delta, delta, 1 + delta)}; + PointType points[] = {PointType {1, 0, 0}, + PointType {0, 1, 0}, + PointType {0, 0, 1}, + PointType {1 + delta, delta, delta}, + PointType {delta, 1 + delta, delta}, + PointType {delta, delta, 1 + delta}}; BoundingBoxType bbox; bbox.addPoint(PointType::zero()); @@ -150,25 +152,25 @@ TEST(primal_clip, boundingBoxOptimization) PointType midpoint = PointType::zero(); PointType points[] = { - PointType::make_point(VAL1, VAL2, 0), - PointType::make_point(-VAL1, VAL2, 0), - PointType::make_point(VAL1, -VAL2, 0), - PointType::make_point(-VAL1, -VAL2, 0), - - PointType::make_point(VAL1, 0, VAL2), - PointType::make_point(-VAL1, 0, VAL2), - PointType::make_point(VAL1, 0, -VAL2), - PointType::make_point(-VAL1, 0, -VAL2), - - PointType::make_point(0, VAL2, VAL1), - PointType::make_point(0, VAL2, -VAL1), - PointType::make_point(0, -VAL2, VAL1), - PointType::make_point(0, -VAL2, -VAL1), - - PointType::make_point(0, VAL1, VAL2), - PointType::make_point(0, -VAL1, VAL2), - PointType::make_point(0, VAL1, -VAL2), - PointType::make_point(0, -VAL1, -VAL2), + PointType {VAL1, VAL2, 0}, + PointType {-VAL1, VAL2, 0}, + PointType {VAL1, -VAL2, 0}, + PointType {-VAL1, -VAL2, 0}, + + PointType {VAL1, 0, VAL2}, + PointType {-VAL1, 0, VAL2}, + PointType {VAL1, 0, -VAL2}, + PointType {-VAL1, 0, -VAL2}, + + PointType {0, VAL2, VAL1}, + PointType {0, VAL2, -VAL1}, + PointType {0, -VAL2, VAL1}, + PointType {0, -VAL2, -VAL1}, + + PointType {0, VAL1, VAL2}, + PointType {0, -VAL1, VAL2}, + PointType {0, VAL1, -VAL2}, + PointType {0, -VAL1, -VAL2}, }; for(int i = 0; i < 16; i += 2) @@ -187,20 +189,20 @@ TEST(primal_clip, experimentalData) const double EPS = 1e-8; // Triangle 248 from sphere mesh - TriangleType tri(PointType::make_point(0.405431, 3.91921, 3.07821), - PointType::make_point(1.06511, 3.96325, 2.85626), - PointType::make_point(0.656002, 4.32465, 2.42221)); + TriangleType tri(PointType {0.405431, 3.91921, 3.07821}, + PointType {1.06511, 3.96325, 2.85626}, + PointType {0.656002, 4.32465, 2.42221}); // Block index {grid pt: (19,29,24); level: 5} from InOutOctree - BoundingBoxType box12(PointType::make_point(0.937594, 4.06291, 2.50025), - PointType::make_point(1.25012, 4.37544, 2.81278)); + BoundingBoxType box12(PointType {0.937594, 4.06291, 2.50025}, + PointType {1.25012, 4.37544, 2.81278}); PolygonType poly = axom::primal::clip(tri, box12); EXPECT_EQ(3, poly.numVertices()); SLIC_INFO("Intersection of triangle " << tri << " \n\t and bounding box " << box12 << " \n\t is polygon" - << poly << " with centroid " << poly.centroid()); + << poly << " with centroid " << poly.vertexMean()); // Check that the polygon vertices are on the triangle for(int i = 0; i < poly.numVertices(); ++i) @@ -226,7 +228,7 @@ TEST(primal_clip, experimentalData) // Check that the polygon centroid is on the triangle { - PointType centroid = poly.centroid(); + PointType centroid = poly.vertexMean(); PointType bary = tri.physToBarycentric(centroid); PointType reconstructed = tri.baryToPhysical(bary); diff --git a/src/axom/primal/tests/primal_compute_moments.cpp b/src/axom/primal/tests/primal_compute_moments.cpp new file mode 100644 index 0000000000..0fb0d018f0 --- /dev/null +++ b/src/axom/primal/tests/primal_compute_moments.cpp @@ -0,0 +1,251 @@ +// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file primal_compute_moments.cpp + * \brief This file tests primal's functionality related to computing moments + */ + +#include "gtest/gtest.h" + +#include "axom/core.hpp" +#include "axom/slic.hpp" + +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/CurvedPolygon.hpp" +#include "axom/primal/operators/compute_moments.hpp" +#include "axom/primal/operators/detail/compute_moments_impl.hpp" + +namespace primal = axom::primal; + +const double EPS = 2e-15; + +//------------------------------------------------------------------------------ +TEST(primal_compute_moments, sector_area_cubic) +{ + const int DIM = 2; + using T = double; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + using axom::utilities::isNearlyEqual; + + { + SLIC_INFO("Testing Bezier sector area calculation for a cubic"); + const int order = 3; + PointType data[order + 1] = {PointType {0.6, 1.2}, + PointType {1.3, 1.6}, + PointType {2.9, 2.4}, + PointType {3.2, 3.5}}; + + BezierCurveType bCurve(data, order); + const T area = primal::sector_area(bCurve); + + EXPECT_NEAR(.1455, area, EPS); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_compute_moments, sector_moment_cubic) +{ + const int DIM = 2; + using CoordType = double; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + + { + SLIC_INFO("Testing Bezier sector moment calculation for a cubic"); + const int order = 3; + PointType data[order + 1] = {PointType {0.6, 1.2}, + PointType {1.3, 1.6}, + PointType {2.9, 2.4}, + PointType {3.2, 3.5}}; + + BezierCurveType bCurve(data, order); + PointType M = primal::sector_centroid(bCurve); + EXPECT_NEAR(-.429321428571429, M[0], EPS); + EXPECT_NEAR(-.354010714285715, M[1], EPS); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_compute_moments, sector_area_point) +{ + const int DIM = 2; + using CoordType = double; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + + { + SLIC_INFO("Testing Bezier sector area calculation for a point"); + const int order = 0; + PointType data[order + 1] = {PointType {0.6, 1.2}}; + + BezierCurveType bCurve(data, order); + EXPECT_DOUBLE_EQ(0., primal::sector_area(bCurve)); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_compute_moments, sector_moment_point) +{ + const int DIM = 2; + using CoordType = double; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + + { + SLIC_INFO("Testing Bezier sector moment calculation for a point"); + const int order = 0; + PointType data[order + 1] = {PointType {0.6, 1.2}}; + + BezierCurveType bCurve(data, order); + PointType M = primal::sector_centroid(bCurve); + EXPECT_DOUBLE_EQ(M[0], 0.0); + EXPECT_DOUBLE_EQ(M[1], 0.0); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_compute_moments, sector_weights) +{ + SLIC_INFO("Testing weights for BezierCurve::sectorArea()"); + + // NOTE: Expected weights are provided in the reference paper [Ueda99] + // See doxygen comment for primal::sector_area(BezierCurve) + + using CoordType = double; + primal::detail::MemoizedSectorAreaWeights memoizedSectorWeights; + + // order 1 + { + const int ord = 1; + auto weights = memoizedSectorWeights.getWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(2, 1); + axom::numerics::Matrix exp(ord + 1, ord + 1); + // clang-format off + exp(0,0) = 0; exp(0,1) = 1; + exp(1,0) = -1; exp(1,1) = 0; + // clang-format on + + for(int i = 0; i <= ord; ++i) + { + for(int j = 0; j <= ord; ++j) + { + EXPECT_DOUBLE_EQ(exp(i, j) * binomInv, weights(i, j)); + } + } + } + + // order 2 + { + const int ord = 2; + auto weights = memoizedSectorWeights.getWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(4, 2); + axom::numerics::Matrix exp(ord + 1, ord + 1); + // clang-format off + exp(0,0) = 0; exp(0,1) = 2; exp(0,2) = 1; + exp(1,0) = -2; exp(1,1) = 0; exp(1,2) = 2; + exp(2,0) = -1; exp(2,1) = -2; exp(2,2) = 0; + // clang-format on + + for(int i = 0; i <= ord; ++i) + { + for(int j = 0; j <= ord; ++j) + { + EXPECT_DOUBLE_EQ(exp(i, j) * binomInv, weights(i, j)); + } + } + } + + // order 3 + { + const int ord = 3; + auto weights = memoizedSectorWeights.getWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(6, 3); + axom::numerics::Matrix exp(ord + 1, ord + 1); + // clang-format off + exp(0,0) = 0; exp(0,1) = 6; exp(0,2) = 3; exp(0,3) = 1; + exp(1,0) = -6; exp(1,1) = 0; exp(1,2) = 3; exp(1,3) = 3; + exp(2,0) = -3; exp(2,1) = -3; exp(2,2) = 0; exp(2,3) = 6; + exp(3,0) = -1; exp(3,1) = -3; exp(3,2) = -6; exp(3,3) = 0; + // clang-format on + + for(int i = 0; i <= ord; ++i) + { + for(int j = 0; j <= ord; ++j) + { + EXPECT_DOUBLE_EQ(exp(i, j) * binomInv, weights(i, j)); + } + } + } + + // order 4 + { + const int ord = 4; + auto weights = memoizedSectorWeights.getWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(8, 4); + axom::numerics::Matrix exp(ord + 1, ord + 1); + // clang-format off + exp(0,0) = 0; exp(0,1) = 20; exp(0,2) = 10; exp(0,3) = 4; exp(0,4) = 1; + exp(1,0) =-20; exp(1,1) = 0; exp(1,2) = 8; exp(1,3) = 8; exp(1,4) = 4; + exp(2,0) =-10; exp(2,1) = -8; exp(2,2) = 0; exp(2,3) = 8; exp(2,4) = 10; + exp(3,0) = -4; exp(3,1) = -8; exp(3,2) = -8; exp(3,3) = 0; exp(3,4) = 20; + exp(4,0) = -1; exp(4,1) = -4; exp(4,2) =-10; exp(4,3) =-20; exp(4,4) = 0; + // clang-format on + + for(int i = 0; i <= ord; ++i) + { + for(int j = 0; j <= ord; ++j) + { + EXPECT_DOUBLE_EQ(exp(i, j) * binomInv, weights(i, j)); + } + } + } + + // order 5 + { + const int ord = 5; + auto weights = memoizedSectorWeights.getWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(10, 5); + axom::numerics::Matrix exp(ord + 1, ord + 1); + // clang-format off + exp(0,0) = 0; exp(0,1) = 70; exp(0,2) = 35; exp(0,3) = 15; exp(0,4) = 5; exp(0,5) = 1; + exp(1,0) =-70; exp(1,1) = 0; exp(1,2) = 25; exp(1,3) = 25; exp(1,4) = 15; exp(1,5) = 5; + exp(2,0) =-35; exp(2,1) =-25; exp(2,2) = 0; exp(2,3) = 20; exp(2,4) = 25; exp(2,5) = 15; + exp(3,0) =-15; exp(3,1) =-25; exp(3,2) =-20; exp(3,3) = 0; exp(3,4) = 25; exp(3,5) = 35; + exp(4,0) = -5; exp(4,1) =-15; exp(4,2) =-25; exp(4,3) =-25; exp(4,4) = 0; exp(4,5) = 70; + exp(5,0) = -1; exp(5,1) = -5; exp(5,2) =-15; exp(5,3) =-35; exp(5,4) =-70; exp(5,5) = 0; + // clang-format on + + for(int i = 0; i <= ord; ++i) + { + for(int j = 0; j <= ord; ++j) + { + EXPECT_DOUBLE_EQ(exp(i, j) * binomInv, weights(i, j)); + } + } + } +} + +//------------------------------------------------------------------------------ + +int main(int argc, char* argv[]) +{ + int result = 0; + + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::SimpleLogger logger; + + result = RUN_ALL_TESTS(); + + return result; +} diff --git a/src/axom/primal/tests/primal_curved_polygon.cpp b/src/axom/primal/tests/primal_curved_polygon.cpp new file mode 100644 index 0000000000..ef0ef95373 --- /dev/null +++ b/src/axom/primal/tests/primal_curved_polygon.cpp @@ -0,0 +1,526 @@ +// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file primal_curved_polygon.cpp + * \brief This file tests the CurvedPolygon class + */ + +#include "gtest/gtest.h" + +#include "axom/config.hpp" +#include "axom/slic.hpp" +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/Segment.hpp" +#include "axom/primal/geometry/CurvedPolygon.hpp" +#include "axom/primal/geometry/OrientationResult.hpp" +#include "axom/primal/operators/intersect.hpp" +#include "axom/primal/operators/compute_moments.hpp" +#include "axom/primal/operators/orientation.hpp" + +#include +#include + +namespace primal = axom::primal; + +/*! + * Helper function to compute the area and centroid of a curved polygon and to check that they match expectations, + * stored in \a expArea and \a expCentroid. Areas and Moments are computed within tolerance \a eps and checks use \a test_eps. + */ +template +void checkMoments(const primal::CurvedPolygon& bPolygon, + const CoordType expArea, + const primal::Point& expMoment, + double eps, + double test_eps) +{ + EXPECT_NEAR(expArea, primal::area(bPolygon, eps), test_eps); + + const auto centroid = primal::centroid(bPolygon, eps); + for(int i = 0; i < 2; ++i) + { + EXPECT_NEAR(expMoment[i], centroid[i], test_eps); + } +} + +/*! + * Helper function to create a CurvedPolygon from a list of control points and a list of orders of component curves. + * Control points should be given as a list of Points in order of orientation with no duplicates except that + * the first control point should also be the last control point (if the polygon is closed). + * Orders should be given as a list of ints in order of orientation, representing the orders of the component curves. + */ +template +primal::CurvedPolygon createPolygon( + const std::vector> ControlPoints, + const std::vector orders) +{ + using PointType = primal::Point; + using CurvedPolygonType = primal::CurvedPolygon; + using BezierCurveType = primal::BezierCurve; + + const int num_edges = orders.size(); + const int num_unique_control_points = ControlPoints.size(); + + //checks if the orders and control points given will yield a valid curved polygon + { + const int sum_of_orders = std::accumulate(orders.begin(), orders.end(), 0); + EXPECT_EQ(sum_of_orders + 1, num_unique_control_points); + } + + //Converts the control points to BezierCurves of specified orders and stores them in a CurvedPolygon object. + CurvedPolygonType bPolygon; + int iter = 0; + for(int j = 0; j < num_edges; ++j) + { + std::vector subCP; + subCP.assign(ControlPoints.begin() + iter, + ControlPoints.begin() + iter + orders[j] + 1); + BezierCurveType addCurve(subCP, orders[j]); + bPolygon.addEdge(addCurve); + iter += (orders[j]); + } + return bPolygon; +} + +//------------------------------------------------------------------------------ +TEST(primal_curvedpolygon, constructor) +{ + const int DIM = 3; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using BezierCurveType = primal::BezierCurve; + + { + SLIC_INFO("Testing default CurvedPolygon constructor "); + CurvedPolygonType bPolygon; + + int expNumEdges = 0; + EXPECT_EQ(expNumEdges, bPolygon.numEdges()); + EXPECT_EQ(expNumEdges, bPolygon.getEdges().size()); + EXPECT_EQ(std::vector(), bPolygon.getEdges()); + } + + { + SLIC_INFO("Testing CurvedPolygon numEdges constructor "); + + CurvedPolygonType bPolygon(1); + int expNumEdges = 1; + EXPECT_EQ(expNumEdges, bPolygon.numEdges()); + EXPECT_EQ(expNumEdges, static_cast(bPolygon.getEdges().size())); + } +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, add_edges) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + + SLIC_INFO("Test adding edges to empty CurvedPolygon"); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + PointType controlPoints[2] = {PointType {0.6, 1.2}, PointType {0.0, 1.6}}; + + BezierCurveType bCurve(controlPoints, 1); + + bPolygon.addEdge(bCurve); + bPolygon.addEdge(bCurve); + + EXPECT_EQ(2, bPolygon.numEdges()); + for(int p = 0; p < bPolygon.numEdges(); ++p) + { + BezierCurveType& bc = bPolygon[p]; + for(int sz = 0; sz <= bc.getOrder(); ++sz) + { + auto& pt = bc[sz]; + for(int i = 0; i < DIM; ++i) + { + EXPECT_DOUBLE_EQ(controlPoints[sz][i], pt[i]); + } + } + } +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, isClosed) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test checking if CurvedPolygon is closed."); + + { + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + EXPECT_FALSE(bPolygon.isClosed()); + } + + std::vector CP = {PointType {0.6, 1.2}, + PointType {0.3, 2.0}, + PointType {0.0, 1.6}, + PointType {0.6, 1.2}}; + std::vector orders = {1, 1, 1}; + + { + std::vector subCP = {PointType {0.6, 1.2}, PointType {0.3, 2.0}}; + std::vector suborders = {1}; + CurvedPolygonType subPolygon = createPolygon(subCP, suborders); + EXPECT_FALSE(subPolygon.isClosed()); + } + + { + CurvedPolygonType bPolygon = createPolygon(CP, orders); + EXPECT_EQ(3, bPolygon.numEdges()); + EXPECT_TRUE(bPolygon.isClosed()); + + bPolygon[2][1][0] -= 2e-15; + EXPECT_FALSE(bPolygon.isClosed(1e-15)); + } + + { + CurvedPolygonType bPolygon = createPolygon(CP, orders); + + bPolygon[1][0][0] = 5; + EXPECT_FALSE(bPolygon.isClosed(1e-15)); + } +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, isClosed_BiGon) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test checking if CurvedPolygon is closed for a Bi-Gon"); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + EXPECT_FALSE(bPolygon.isClosed()); + + // Bi-gon defined by a quadratic edge and a straight line + std::vector CP = {PointType {0.8, .25}, + PointType {2.0, .50}, + PointType {0.8, .75}, + PointType {0.8, .25}}; + std::vector orders = {2, 1}; + + CurvedPolygonType poly = createPolygon(CP, orders); + EXPECT_TRUE(poly.isClosed()); + + // modify a vertex of the quadratic and check again + CurvedPolygonType poly2 = poly; + poly2[0][2] = PointType {0.8, 1.0}; + EXPECT_FALSE(poly2.isClosed()); +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, split_edge) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + + SLIC_INFO("Test checking CurvedPolygon edge split."); + + std::vector CP = {PointType {0.6, 1.2}, + PointType {0.3, 2.0}, + PointType {0.0, 1.6}, + PointType {0.6, 1.2}}; + + std::vector orders32 = {1, 1, 1}; + CurvedPolygonType bPolygon32 = createPolygon(CP, orders32); + std::cout << "Got here!! " << std::endl; + std::vector subCP; + + subCP.assign(CP.begin(), CP.begin() + 2); + BezierCurveType bCurve(subCP, 1); + bPolygon32.splitEdge(0, .5); + + BezierCurveType bCurve2; + BezierCurveType bCurve3; + bCurve.split(.5, bCurve2, bCurve3); + + EXPECT_EQ(bPolygon32.numEdges(), 4); + for(int i = 0; i < bPolygon32[0].getOrder(); ++i) + { + for(int dimi = 0; dimi < DIM; ++dimi) + { + EXPECT_EQ(bPolygon32[0][i][dimi], bCurve2[i][dimi]); + EXPECT_EQ(bPolygon32[1][i][dimi], bCurve3[i][dimi]); + } + } +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, moments_triangle_degenerate) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + + SLIC_INFO("Testing area computation of degenerate triangles"); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + EXPECT_EQ(0.0, primal::area(bPolygon)); + PointType origin {0.0, 0.0}; + + PointType controlPoints[2] = {PointType {0.6, 1.2}, PointType {0.3, 2.0}}; + + PointType controlPoints2[2] = {PointType {0.3, 2.0}, PointType {0.0, 1.6}}; + PointType controlPoints3[2] = {PointType {0.0, 1.6}, PointType {0.6, 1.2}}; + + BezierCurveType bCurve(controlPoints, 1); + bPolygon.addEdge(bCurve); + checkMoments(bPolygon, 0.0, origin, 1e-14, 1e-15); + + BezierCurveType bCurve2(controlPoints2, 1); + bPolygon.addEdge(bCurve2); + checkMoments(bPolygon, 0.0, origin, 1e-14, 1e-15); + + BezierCurveType bCurve3(controlPoints3, 1); + bPolygon.addEdge(bCurve3); + + bPolygon[2][1][0] -= 1e-11; + checkMoments(bPolygon, 0.0, origin, 1e-14, 1e-15); +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, moments_triangle_linear) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test moment computation of a linear triangle"); + std::vector CP = {PointType {0.6, 1.2}, + PointType {0.3, 2.0}, + PointType {0.0, 1.6}, + PointType {0.6, 1.2}}; + + std::vector orders = {1, 1, 1}; + CurvedPolygonType bPolygon = createPolygon(CP, orders); + + CoordType trueA = -.18; + PointType trueC = PointType::make_point(0.3, 1.6); + + checkMoments(bPolygon, trueA, trueC, 1e-14, 1e-15); +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, moments_triangle_quadratic) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test moment computation of quadratic triangle"); + + std::vector CP = {PointType {0.6, 1.2}, + PointType {0.4, 1.3}, + PointType {0.3, 2.0}, + PointType {0.27, 1.5}, + PointType {0.0, 1.6}, + PointType {0.1, 1.5}, + PointType {0.6, 1.2}}; + + std::vector orders = {2, 2, 2}; + CurvedPolygonType bPolygon = createPolygon(CP, orders); + + CoordType trueA = -0.097333333333333; + PointType trueC {.294479452054794, 1.548219178082190}; + + checkMoments(bPolygon, trueA, trueC, 1e-15, 1e-14); +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, moments_triangle_mixed_order) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO( + "Test moment computation for curved triangle with mixed order edges"); + + std::vector CP = {PointType {0.6, 1.2}, + PointType {0.4, 1.3}, + PointType {0.3, 2.0}, + PointType {0.27, 1.5}, + PointType {0.0, 1.6}, + PointType {0.6, 1.2}}; + + std::vector orders = {2, 2, 1}; + CurvedPolygonType bPolygon = createPolygon(CP, orders); + + CoordType trueA = -.0906666666666666666666; + PointType trueC {.2970147058823527, 1.55764705882353}; + + checkMoments(bPolygon, trueA, trueC, 1e-14, 1e-15); +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, moments_quad_all_orders) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test moment computation for quads of different orders"); + std::vector CPorig = {PointType {0.0, 0.0}, + PointType {0.0, 1.0}, + PointType {1.0, 1.0}, + PointType {1.0, 0.0}, + PointType {0.0, 0.0}}; + + std::vector orders = {1, 1, 1, 1}; + CurvedPolygonType bPolygon = createPolygon(CPorig, orders); + + CoordType trueA = 1.0; + PointType trueC = PointType::make_point(0.5, 0.5); + + checkMoments(bPolygon, trueA, trueC, 1e-14, 1e-15); + for(int p = 2; p < 11; ++p) + { + std::vector CP = CPorig; + for(int side = 0; side < 4; ++side) + { + for(int i = 1; i < p; ++i) + { + const int offset = i + (side * p); + const double t = static_cast(i) / p; + switch(side) + { + case 0: + CP.insert(CP.begin() + offset, PointType {0., t}); + break; + case 1: + CP.insert(CP.begin() + offset, PointType {t, 1.}); + break; + case 2: + CP.insert(CP.begin() + offset, PointType {1., 1. - t}); + break; + case 3: + CP.insert(CP.begin() + offset, PointType {1. - t, 0.}); + break; + } + } + orders[side] += 1; + } + // for (int i=0; i; + using PointType = primal::Point; + using SegmentType = primal::Segment; + using BezierCurveType = primal::BezierCurve; + + // Test several n-gons discretizing the unit circle + const int MAX_SEG = 10; + const PointType origin; + for(int nseg = 3; nseg < MAX_SEG; ++nseg) + { + // Create an n-gon with line segments going CCW along the unit circle + CurvedPolygonType poly(nseg); + axom::Array pts(nseg + 1); + for(int i = 0; i < nseg; ++i) + { + const double theta = i / static_cast(nseg) * (2. * M_PI); + pts[i] = PointType {cos(theta), sin(theta)}; + } + pts[nseg] = pts[0]; + + for(int i = 0; i < nseg; ++i) + { + poly[i] = BezierCurveType(&pts[i], order); + } + EXPECT_TRUE(poly.isClosed()); + + // Perform some checks on the polygon + for(int i = 0; i < nseg; ++i) + { + // check that the end point of each segment is equal to the start of the next + auto& currentEnd = poly[i][order]; + auto& nextStart = poly[(i + 1) % nseg][0]; + EXPECT_EQ(currentEnd, nextStart); + + // check that the orientation of segment midpoints goes in the same direction + SegmentType seg(poly[i].evaluate(0.5), poly[(i + 1) % nseg].evaluate(0.5)); + EXPECT_EQ(primal::ON_NEGATIVE_SIDE, primal::orientation(origin, seg)); + } + + // Create a polygon with reversed orientation + CurvedPolygonType reversed = poly; + reversed.reverseOrientation(); + + EXPECT_EQ(poly.numEdges(), reversed.numEdges()); + + // Perform some checks on the reversed polygon + for(int i = 0; i < nseg; ++i) + { + // check that order of each segment stayed the same + EXPECT_EQ(poly[i].getOrder(), reversed[i].getOrder()); + + // check that the end point of each segment is equal to the start of the next + auto& currentEnd = reversed[i][order]; + auto& nextStart = reversed[(i + 1) % nseg][0]; + EXPECT_EQ(currentEnd, nextStart); + + // check that segment midpoints are oriented in same direction (opposite of origin) + SegmentType seg(reversed[i].evaluate(0.5), + reversed[(i + 1) % nseg].evaluate(0.5)); + EXPECT_EQ(primal::ON_POSITIVE_SIDE, primal::orientation(origin, seg)); + } + + // Check that reversing twice yields the original + CurvedPolygonType reversedAgain = reversed; + reversedAgain.reverseOrientation(); + EXPECT_EQ(poly, reversedAgain); + } +} + +//---------------------------------------------------------------------------------- +int main(int argc, char* argv[]) +{ + int result = 0; + + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::SimpleLogger logger; + + result = RUN_ALL_TESTS(); + + return result; +} diff --git a/src/axom/quest/InOutOctree.hpp b/src/axom/quest/InOutOctree.hpp index 0111786f8e..9bf9fd3fc8 100644 --- a/src/axom/quest/InOutOctree.hpp +++ b/src/axom/quest/InOutOctree.hpp @@ -1123,7 +1123,7 @@ typename std::enable_if::type InOutOctree::withinGrayBlock } } - triPt = poly.centroid(); + triPt = poly.vertexMean(); /// Use a ray from the query point to the triangle point to find an /// intersection. Note: We have to check all triangles to ensure that diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index 65cb8423e4..591a307d7c 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -8,6 +8,8 @@ set(quest_example_depends axom fmt cli11) +blt_list_append(TO quest_example_depends ELEMENTS mpi IF ENABLE_MPI) +blt_list_append(TO quest_example_depends ELEMENTS openmp IF ENABLE_OPENMP) blt_list_append(TO quest_example_depends ELEMENTS cuda IF ENABLE_CUDA) blt_list_append(TO quest_example_depends ELEMENTS blt::hip IF ENABLE_HIP) blt_list_append(TO quest_example_depends ELEMENTS RAJA IF RAJA_FOUND) diff --git a/src/axom/quest/readers/C2CReader.cpp b/src/axom/quest/readers/C2CReader.cpp index 3bfc04c2c0..91c98981df 100644 --- a/src/axom/quest/readers/C2CReader.cpp +++ b/src/axom/quest/readers/C2CReader.cpp @@ -17,9 +17,6 @@ namespace axom { namespace quest { -/// returns the linear interpolation of \a A and \a B at \a t. i.e. (1-t)A+tB -inline double lerp(double A, double B, double t) { return (1 - t) * A + t * B; } - /*! * \brief Helper class for interpolating points on a NURBS curve * @@ -282,6 +279,8 @@ void C2CReader::log() void C2CReader::getLinearMesh(mint::UnstructuredMesh* mesh, int segmentsPerKnotSpan) { + using axom::utilities::lerp; + // Sanity checks SLIC_ERROR_IF(mesh == nullptr, "supplied mesh is null!"); SLIC_ERROR_IF(mesh->getDimension() != 2, "C2C reader expects a 2D mesh!"); diff --git a/src/axom/quest/tests/CMakeLists.txt b/src/axom/quest/tests/CMakeLists.txt index 46e0bd292e..4666415b43 100644 --- a/src/axom/quest/tests/CMakeLists.txt +++ b/src/axom/quest/tests/CMakeLists.txt @@ -31,6 +31,7 @@ set(quest_tests_depends gtest ) +blt_list_append( TO quest_tests_depends ELEMENTS mpi IF ENABLE_MPI ) blt_list_append( TO quest_tests_depends ELEMENTS cuda IF ENABLE_CUDA ) blt_list_append( TO quest_tests_depends ELEMENTS blt::hip IF ENABLE_HIP )