Skip to content

Commit

Permalink
GRIDEDIT-1566 Refactored triangulation and sample interpolators
Browse files Browse the repository at this point in the history
  • Loading branch information
@BillSenior
BillSenior committed Jan 6, 2025
1 parent ac960cd commit 018ba63
Show file tree
Hide file tree
Showing 12 changed files with 253 additions and 372 deletions.
66 changes: 0 additions & 66 deletions libs/MeshKernel/include/MeshKernel/Definitions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,70 +113,4 @@ namespace meshkernel
/// @brief Get the string representation of the CurvilinearDirection enumeration values.
const std::string& CurvilinearDirectionToString(CurvilinearDirection direction);

/// @brief The concept specifies that the array type must have an access operator returning the array element type or can be converted to one
template <typename ArrayType, typename ResultType>
concept ArrayConstAccessConcept = requires(const ArrayType& array, const size_t i) {
{
array[i]
} -> std::convertible_to<ResultType>;
};

/// @brief The concept specifies that the array type must have an access operator returning a reference to the array element type
template <typename ArrayType, typename ResultType>
concept ArrayNonConstAccessConcept = requires(ArrayType& array, const size_t i) {
{
array[i]
} -> std::same_as<ResultType&>;
};

/// @brief A concept that specifies that the array must have a size function return the number of elements in the array
template <typename ArrayType>
concept ArraySizeConcept = requires(const ArrayType& array) {
{
array.size()
} -> std::same_as<size_t>;
};

/// @brief A concept that specifies that the array must have a begin and end function.
///
/// Would like to also specify the return type here, but span needs some c++23 functionality here.
/// Then change all iterator usage to cbegin and cend returning a const_iterator
/// std::same_as<typename ArrayType::const_iterator>
template <typename ArrayType>
concept ArrayConstIteratorsConcept = requires(const ArrayType& array) {
{
array.begin()
};
{
array.end()
};
};

/// @brief A concept that specifies that the array must have a begin and end function.
///
/// Would like to also specify the return type here, but span needs some c++23 functionality here.
/// Then change all iterator usage to cbegin and cend returning a const_iterator
/// std::same_as<typename ArrayType::const_iterator>
template <typename ArrayType>
concept ArrayNonConstIteratorsConcept = requires(ArrayType& array) {
{
array.begin()
} -> std::same_as<typename ArrayType::iterator>;
{
array.end()
} -> std::same_as<typename ArrayType::iterator>;
};

/// @brief A concept that specifies all the functionality required to be usable as a constant array of doubles.
template <typename ArrayType>
concept ValidConstDoubleArray = ArrayConstAccessConcept<ArrayType, double> &&
ArrayConstIteratorsConcept<ArrayType> &&
ArraySizeConcept<ArrayType>;

/// @brief A concept that specifies all the functionality required to be usable as a constant array of doubles.
template <typename ArrayType>
concept ValidNonConstDoubleArray = ArrayNonConstAccessConcept<ArrayType, double> &&
ArrayNonConstIteratorsConcept<ArrayType> &&
ArraySizeConcept<ArrayType>;

} // namespace meshkernel
69 changes: 4 additions & 65 deletions libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
#include <array>
#include <functional>
#include <memory>
#include <span>
#include <string>
#include <vector>

Expand All @@ -49,9 +50,6 @@ namespace meshkernel
class BoundedArray
{
public:
/// @brief Constructor
BoundedArray() : m_size(0) {}

/// @brief Number of elements in the array
UInt size() const
{
Expand All @@ -66,7 +64,6 @@ namespace meshkernel
throw ConstraintError("array already at size limit: {}", Dimension);
}

// m_indices.push_back(index);
m_indices[m_size] = index;
++m_size;
}
Expand Down Expand Up @@ -117,14 +114,12 @@ namespace meshkernel
{
public:
/// @brief Constructor with array of points
template <meshkernel::ValidConstPointArray PointVector>
MeshTriangulation(const PointVector& nodes,
MeshTriangulation(const std::span<const Point> nodes,
const Projection projection);

/// @brief Constructor with separate arrays of x- and y-coordinates
template <meshkernel::ValidConstDoubleArray VectorType>
MeshTriangulation(const VectorType& xNodes,
const VectorType& yNodes,
MeshTriangulation(const std::span<const double> xNodes,
const std::span<const double> yNodes,
const Projection projection);

/// @brief Get the projection type used in the triangulation
Expand Down Expand Up @@ -199,62 +194,6 @@ inline meshkernel::Projection meshkernel::MeshTriangulation::GetProjection() con
return m_projection;
}

template <meshkernel::ValidConstPointArray PointVector>
meshkernel::MeshTriangulation::MeshTriangulation(const PointVector& nodes,
const Projection projection)
: m_nodes(nodes.begin(), nodes.end()),
m_projection(projection)
{
if (nodes.size() < constants::geometric::numNodesInTriangle)
{
throw ConstraintError("Not enough nodes for a single triangle: {}", nodes.size());
}

std::vector<double> xNodes(nodes.size());
std::vector<double> yNodes(nodes.size());

std::transform(nodes.begin(), nodes.end(), xNodes.begin(),
[](const Point& p)
{ return p.x; });
std::transform(nodes.begin(), nodes.end(), yNodes.begin(),
[](const Point& p)
{ return p.y; });

std::span<double> xNodesSpan(xNodes.data(), xNodes.size());
std::span<double> yNodesSpan(yNodes.data(), yNodes.size());

Compute(xNodesSpan, yNodesSpan);
}

template <meshkernel::ValidConstDoubleArray VectorType>
meshkernel::MeshTriangulation::MeshTriangulation(const VectorType& xNodes,
const VectorType& yNodes,
const Projection projection)
: m_nodes(xNodes.size()),
m_projection(projection)
{
if (xNodes.size() < constants::geometric::numNodesInTriangle)
{
throw ConstraintError("Not enough nodes for a single triangle: {}", xNodes.size());
}

if (xNodes.size() != yNodes.size())
{
throw ConstraintError("Size mismatch between x- and y-node values: {} /= {}",
xNodes.size(), yNodes.size());
}

for (size_t i = 0; i < xNodes.size(); ++i)
{
m_nodes[i] = Point{xNodes[i], yNodes[i]};
}

std::span<const double> xNodesSpan(xNodes.data(), xNodes.size());
std::span<const double> yNodesSpan(yNodes.data(), yNodes.size());

Compute(xNodesSpan, yNodesSpan);
}

inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfNodes() const
{
return static_cast<UInt>(m_nodes.size());
Expand Down
123 changes: 2 additions & 121 deletions libs/MeshKernel/include/MeshKernel/Operations.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -332,9 +332,8 @@ namespace meshkernel
/// @param[in] triangleNodes A series of node forming an open triangle
/// @param[in] projection The coordinate system projection.
/// @returns If point is inside the designated triangle
template <class PointVector>
[[nodiscard]] bool IsPointInTriangle(const Point& point,
const PointVector& triangleNodes,
const std::span<const Point> triangleNodes,
const Projection& projection);

/// @brief Computes three base components
Expand Down Expand Up @@ -595,8 +594,7 @@ namespace meshkernel
/// @param[in] points The point series.
/// @param[in] projection The projection to use.
/// @return The average coordinate.
template <ValidConstPointArray PointVector>
[[nodiscard]] Point ComputeAverageCoordinate(const PointVector& points, const Projection& projection);
[[nodiscard]] Point ComputeAverageCoordinate(const std::span<const Point> points, const Projection& projection);

/// @brief Cartesian projection of a point on a segment defined by other two points
/// @param firstNode The first node of the segment
Expand Down Expand Up @@ -675,120 +673,3 @@ inline meshkernel::Cartesian3DPoint meshkernel::operator*(const Cartesian3DPoint
{
return value * p;
}

namespace meshkernel
{
template <ValidConstPointArray PointVector>
[[nodiscard]] Point ComputeAverageCoordinate(const PointVector& points, const Projection& projection)
{
size_t validCount = std::ranges::count_if(points, [](const Point& p)
{ return p.IsValid(); });

if (projection == Projection::sphericalAccurate)
{

UInt firstValidPoint = 0;

if (validCount != points.size())
{
auto iterator = std::ranges::find_if(points, [](const Point& p)
{ return p.IsValid(); });
firstValidPoint = static_cast<UInt>(iterator - points.begin());
}

Cartesian3DPoint averagePoint3D{0.0, 0.0, 0.0};
for (const auto& point : points)
{
if (!point.IsValid())
{
continue;
}

const Cartesian3DPoint point3D{SphericalToCartesian3D(point)};
averagePoint3D.x += point3D.x;
averagePoint3D.y += point3D.y;
averagePoint3D.z += point3D.z;
}
averagePoint3D.x = averagePoint3D.x / static_cast<double>(validCount);
averagePoint3D.y = averagePoint3D.y / static_cast<double>(validCount);
averagePoint3D.z = averagePoint3D.z / static_cast<double>(validCount);

return Cartesian3DToSpherical(averagePoint3D, points[firstValidPoint].x);
}

auto result = std::accumulate(points.begin(), points.end(), Point{0.0, 0.0}, [](const Point& sum, const Point& current)
{ return current.IsValid() ? sum + current : sum; });
result.x = result.x / static_cast<double>(validCount);
result.y = result.y / static_cast<double>(validCount);
return result;
}
} // namespace meshkernel

template <class PointVector>
[[nodiscard]] bool meshkernel::IsPointInTriangle(const Point& point,
const PointVector& triangleNodes,
const Projection& projection)
{
if (triangleNodes.empty())
{
return true;
}

bool isInTriangle = false;

if (projection == Projection::cartesian || projection == Projection::spherical)
{
int windingNumber = 0;

for (UInt n = 0; n < constants::geometric::numNodesInTriangle; n++)
{
UInt endIndex = n == 2 ? 0 : n + 1;

const auto crossProductValue = crossProduct(triangleNodes[n], triangleNodes[endIndex], triangleNodes[n], point, Projection::cartesian);

if (IsEqual(crossProductValue, 0.0))
{
// point on the line
return true;
}

if (triangleNodes[n].y <= point.y) // an upward crossing
{
if (triangleNodes[endIndex].y > point.y && crossProductValue > 0.0)
{
++windingNumber; // have a valid up intersect
}
}
else
{
if (triangleNodes[endIndex].y <= point.y && crossProductValue < 0.0) // a downward crossing
{
--windingNumber; // have a valid down intersect
}
}
}

isInTriangle = windingNumber == 0 ? false : true;
}

if (projection == Projection::sphericalAccurate)
{
std::vector<Point> triangleCopy;
triangleCopy.reserve(constants::geometric::numNodesInTriangle + 1);
Point centre(0.0, 0.0);

for (UInt i = 0; i < constants::geometric::numNodesInTriangle; ++i)
{
centre += triangleNodes[i];
triangleCopy.emplace_back(triangleNodes[i]);
}

triangleCopy.emplace_back(triangleNodes[0]);

centre *= 1.0 / 3.0;

isInTriangle = IsPointInPolygonNodes(point, triangleCopy, projection, centre);
}

return isInTriangle;
}
6 changes: 0 additions & 6 deletions libs/MeshKernel/include/MeshKernel/Point.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,12 +265,6 @@ namespace meshkernel
/// @brief Get the delta-x and -y in Cartesian coordinate system
Vector GetDeltaCartesian(const Point& p1, const Point& p2);

/// @brief A concept that specifies all the functionality required to be usable as an array of points.
template <typename ArrayType>
concept ValidConstPointArray = ArrayConstAccessConcept<ArrayType, Point> &&
ArrayConstIteratorsConcept<ArrayType> &&
ArraySizeConcept<ArrayType>;

} // namespace meshkernel

inline void meshkernel::Point::SetInvalid()
Expand Down
Loading

0 comments on commit 018ba63

Please sign in to comment.