Skip to content

Commit

Permalink
Add mkernel_mesh2d_refine_ridges_based_on_gridded_samples API function (
Browse files Browse the repository at this point in the history 
#287 | GRIDEDIT-873)
  • Loading branch information
@lucacarniato
lucacarniato authored Jan 25, 2024
1 parent 89f67c2 commit ce074dc
Show file tree
Hide file tree
Showing 11 changed files with 386 additions and 135 deletions.
8 changes: 2 additions & 6 deletions libs/MeshKernel/include/MeshKernel/BilinearInterpolationOnGriddedSamples.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,11 +38,7 @@ namespace meshkernel

/// @brief Defines the iterpolatable data types.
template <typename T>
concept InterpolatableType =
std::same_as<T, short> ||
std::same_as<T, int> ||
std::same_as<T, float> ||
std::same_as<T, double>;
concept InterpolatableType = std::same_as<T, short> || std::same_as<T, float>;

/// @brief A class for performing bilinear interpolation on gridded samples
template <InterpolatableType T>
Expand Down Expand Up @@ -261,7 +257,7 @@ namespace meshkernel
template <InterpolatableType T>
double BilinearInterpolationOnGriddedSamples<T>::GetGriddedValue(UInt columnIndex, UInt rowIndex) const
{
const auto index = rowIndex * m_numXCoord + columnIndex;
const auto index = (m_numYCoord - rowIndex - 1) * m_numXCoord + columnIndex;
return static_cast<double>(m_values[index]);
}

Expand Down
2 changes: 1 addition & 1 deletion libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,7 @@ namespace meshkernel
LandWater = 3
};

public:
/// @brief Enumerator describing the different refinement types
enum class RefinementType
{
Expand All @@ -89,7 +90,6 @@ namespace meshkernel
RidgeDetection = 3
};

public:
/// @brief The constructor for refining based on samples
/// @param[in] mesh The mesh to be refined
/// @param[in] interpolant The averaging interpolation to use
Expand Down
145 changes: 27 additions & 118 deletions libs/MeshKernel/tests/src/MeshRefinementTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,17 +3,18 @@

#include <gtest/gtest.h>

#include <MeshKernel/Mesh2D.hpp>
#include <MeshKernel/MeshRefinement.hpp>
#include <MeshKernel/Parameters.hpp>
#include <MeshKernel/Polygons.hpp>
#include <TestUtils/Definitions.hpp>
#include <TestUtils/MakeMeshes.hpp>
#include <TestUtils/SampleFileReader.hpp>
#include "MeshKernel/Mesh2D.hpp"
#include "MeshKernel/MeshRefinement.hpp"
#include "MeshKernel/Parameters.hpp"
#include "MeshKernel/Polygons.hpp"
#include "TestUtils/Definitions.hpp"
#include "TestUtils/MakeMeshes.hpp"
#include "TestUtils/SampleFileReader.hpp"
#include "TestUtils/SampleGenerator.hpp"

using namespace meshkernel;

TEST(MeshRefinement, FourByFourWithFourSamples)
TEST(MeshRefinement, MeshRefinementRefinementLevels_OnFourByFourWithFourSamples_ShouldRefinemesh)
{
auto mesh = MakeRectangularMeshForTesting(5, 5, 10.0, Projection::cartesian);

Expand Down Expand Up @@ -211,7 +212,7 @@ TEST(MeshRefinement, RefinementOnAFourByFourMeshWithSamplesShouldRefine)
ASSERT_EQ(10, mesh->m_edges[20].second);
}

TEST(MeshRefinement, SmallTriangualMeshTwoSamples)
TEST(MeshRefinement, MeshRefinementRefinementLevels_SmallTriangualMeshTwoSamples_ShouldRefinemesh)
{
// Prepare
auto mesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/SmallTriangularGrid_net.nc");
Expand Down Expand Up @@ -478,7 +479,7 @@ TEST(MeshRefinement, WindowOfRefinementFile)
ASSERT_EQ(326, mesh->m_edges[915].second);
}

TEST(MeshRefinement, WindowOfRefinementFileBasedOnLevels)
TEST(MeshRefinement, MeshRefinementRefinementLevels_OnWindowOfRefinementFile_ShouldRefinemesh)
{
// Prepare
auto mesh = MakeRectangularMeshForTesting(4, 4, 40.0, Projection::cartesian, {197253.0, 442281.0});
Expand Down Expand Up @@ -700,7 +701,7 @@ TEST(MeshRefinement, FourByFourWithFourSamplesSpherical)
ASSERT_EQ(6, mesh->m_edges[47].second);
}

TEST(MeshRefinement, Refine_SphericalMesh_ShouldRefine)
TEST(MeshRefinement, RefinementFileBasedOnLevels_OnSphericalMesh_ShouldRefine)
{
// Prepare
auto mesh = MakeRectangularMeshForTesting(6, 6, 0.0033, Projection::spherical, {41.1, 41.1});
Expand Down Expand Up @@ -819,9 +820,9 @@ TEST(MeshRefinement, BilinearInterpolationWithGriddedSamplesOnLandShouldNotRefin
// Setup
auto mesh = MakeRectangularMeshForTesting(2, 2, 10.0, Projection::cartesian);

std::vector values{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
std::vector<float> values{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
Point origin{-5.0, -5.0};
auto interpolator = std::make_unique<BilinearInterpolationOnGriddedSamples<double>>(*mesh, 2, 2, origin, 10.0, values);
auto interpolator = std::make_unique<BilinearInterpolationOnGriddedSamples<float>>(*mesh, 2, 2, origin, 10.0, values);

MeshRefinementParameters meshRefinementParameters;
meshRefinementParameters.max_num_refinement_iterations = 1;
Expand All @@ -846,9 +847,9 @@ TEST(MeshRefinement, BilinearInterpolationWithGriddedSamplesOnLandAndSeaShouldRe
// Setup
auto mesh = MakeRectangularMeshForTesting(2, 2, 10.0, Projection::cartesian);

std::vector values{-1.0, -2.0, 3.0, -4.0, -5.0, 6.0, 7.0, 8.0, 9.0};
std::vector<float> values{-1.0, -2.0, 3.0, -4.0, -5.0, 6.0, 7.0, 8.0, 9.0};
Point origin{-5.0, -5.0};
auto interpolator = std::make_unique<BilinearInterpolationOnGriddedSamples<double>>(*mesh, 3, 3, origin, 10.0, values);
auto interpolator = std::make_unique<BilinearInterpolationOnGriddedSamples<float>>(*mesh, 3, 3, origin, 10.0, values);

MeshRefinementParameters meshRefinementParameters;
meshRefinementParameters.max_num_refinement_iterations = 1;
Expand All @@ -873,9 +874,9 @@ TEST(MeshRefinement, BilinearInterpolationWithAllGriddedSamplesOnSeaShouldRefine
// Setup
auto mesh = MakeRectangularMeshForTesting(2, 2, 10.0, Projection::cartesian);

std::vector values{-1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0};
std::vector<float> values{-1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0};
Point origin{-5.0, -5.0};
auto interpolator = std::make_unique<BilinearInterpolationOnGriddedSamples<double>>(*mesh, 2, 2, origin, 10.0, values);
auto interpolator = std::make_unique<BilinearInterpolationOnGriddedSamples<float>>(*mesh, 2, 2, origin, 10.0, values);

MeshRefinementParameters meshRefinementParameters;
meshRefinementParameters.max_num_refinement_iterations = 1;
Expand All @@ -897,107 +898,16 @@ TEST(MeshRefinement, BilinearInterpolationWithAllGriddedSamplesOnSeaShouldRefine
ASSERT_EQ(4, mesh->GetNumEdges());
}

enum class RidgeRefinementTestCase
{
GaussianBump = 1,
GaussianWave = 2,
RidgeXDirection = 3,
ArctanFunction = 4,
};

std::vector<Sample> generateSampleData(RidgeRefinementTestCase testcase,
UInt nx = 10,
UInt ny = 10,
double deltaX = 10.0,
double deltaY = 10.0)
{
UInt start = 0;
UInt size = (nx - start) * (ny - start);
std::vector<Sample> sampleData(size);

std::vector sampleDataMatrix(ny, std::vector<double>(nx));

const double centreX = static_cast<double>((nx - 1) / 2) * deltaX;
const double centreY = static_cast<double>((ny - 1) / 2) * deltaY;

const double scale = ny / 4.0 * deltaY;

const double r = nx / 5 * deltaX;
const double maxx = (nx - 1) * deltaX;

std::function<double(double, double)> generateSample;
switch (testcase)
{
case RidgeRefinementTestCase::GaussianBump:
generateSample = [&](double x, double y)
{
const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY);
return 100.0 * std::exp(-0.025 * centre);
};
break;
case RidgeRefinementTestCase::GaussianWave:
generateSample = [&](double x, double y)
{
const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY);
const double factor = std::max(1e-6, std::exp(-0.00025 * centre));
return 100.0 * factor;
};
break;
case RidgeRefinementTestCase::RidgeXDirection:
generateSample = [&](double x, double y)
{
const double sinx = std::sin(x / maxx * M_PI * 4.0);
const double xxx = scale * sinx + centreY;
return 10 * (std::atan(20.0 * (xxx - y)) + M_PI / 2.0);
};
break;
case RidgeRefinementTestCase::ArctanFunction:
generateSample = [&](double x, double y)
{
const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY);
return 10 * (std::atan(20.0 * (r * r - centre)) + M_PI / 2.0);
};
break;
default:
throw std::invalid_argument("invalid ridge refinement test case");
}

for (int i = ny - 1; i >= 0; --i)
{

for (UInt j = start; j < nx; ++j)
{
const double y = deltaY * i;
const double x = deltaX * j;
sampleDataMatrix[ny - 1 - i][j] = generateSample(x, y);
}
}

UInt count = 0;
for (UInt j = start; j < nx; ++j)
{
for (int i = ny - 1; i >= 0; --i)
{
const double y = deltaY * i;
const double x = deltaX * j;
sampleData[count] = {x, y, sampleDataMatrix[ny - 1 - i][j]};
count++;
}
}

return sampleData;
}

class RidgeRefinementTestCases : public testing::TestWithParam<std::tuple<RidgeRefinementTestCase, UInt, UInt>>
class RidgeRefinementTestCases : public testing::TestWithParam<std::tuple<FunctionTestCase, UInt, UInt>>
{
public:
[[nodiscard]] static std::vector<std::tuple<RidgeRefinementTestCase, UInt, UInt>> GetData()
[[nodiscard]] static std::vector<std::tuple<FunctionTestCase, UInt, UInt>> GetData()
{
return std::vector{
std::make_tuple<RidgeRefinementTestCase, UInt, UInt>(RidgeRefinementTestCase::GaussianBump, 1165, 2344),
std::make_tuple<RidgeRefinementTestCase, UInt, UInt>(RidgeRefinementTestCase::GaussianWave, 5297, 10784),
std::make_tuple<RidgeRefinementTestCase, UInt, UInt>(RidgeRefinementTestCase::RidgeXDirection, 2618, 5694),
std::make_tuple<RidgeRefinementTestCase, UInt, UInt>(RidgeRefinementTestCase::ArctanFunction, 2309, 5028)};
std::make_tuple<FunctionTestCase, UInt, UInt>(FunctionTestCase::GaussianBump, 1165, 2344),
std::make_tuple<FunctionTestCase, UInt, UInt>(FunctionTestCase::GaussianWave, 5297, 10784),
std::make_tuple<FunctionTestCase, UInt, UInt>(FunctionTestCase::RidgeXDirection, 2618, 5694),
std::make_tuple<FunctionTestCase, UInt, UInt>(FunctionTestCase::ArctanFunction, 2309, 5028)};
}
};

Expand All @@ -1015,10 +925,9 @@ TEST_P(RidgeRefinementTestCases, expectedResults)
double dimX = (nx - 1) * deltaX;
double dimY = (ny - 1) * deltaY;

std::shared_ptr<Mesh2D> mesh = MakeRectangularMeshForTesting(nx, ny, dimX, dimY, Projection::cartesian);
auto mesh = MakeRectangularMeshForTesting(nx, ny, dimX, dimY, Projection::cartesian);

UInt superSample = 2;

UInt sampleNx = (nx - 1) * superSample + 1;
UInt sampleNy = (ny - 1) * superSample + 1;

Expand All @@ -1027,10 +936,10 @@ TEST_P(RidgeRefinementTestCases, expectedResults)

const auto sampleData = generateSampleData(testCase, sampleNx, sampleNy, sampleDeltaX, sampleDeltaY);

auto samples = SamplesHessianCalculator::ComputeSamplesHessian(sampleData, mesh->m_projection, 0, sampleNx, sampleNy);
auto samplesHessian = SamplesHessianCalculator::ComputeSamplesHessian(sampleData, mesh->m_projection, 0, sampleNx, sampleNy);

auto interpolator = std::make_unique<AveragingInterpolation>(*mesh,
samples,
samplesHessian,
AveragingInterpolation::Method::Max,
Location::Faces,
1.0,
Expand Down
19 changes: 19 additions & 0 deletions libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1231,6 +1231,25 @@ namespace meshkernelapi
int minimumNumSamples,
const meshkernel::MeshRefinementParameters& meshRefinementParameters);

/// @brief Refines a mesh2d based on samples with ridge refinement. This method automatically detects the ridges in a sample set.
///
/// The number of successive splits is indicated on the sample value.
/// For example a value of 0 means no split and hence no refinement, a value of 1 a single split (a quadrilateral face generates 4 faces),
/// a value of 2 two splits (a quadrilateral face generates 16 faces).
/// @param[in] meshKernelId The id of the mesh state
/// @param[in] griddedSamples The gridded samples
/// @param[in] relativeSearchRadius The relative search radius relative to the face size, used for some interpolation algorithms
/// @param[in] minimumNumSamples The minimum number of samples used for some averaging algorithms
/// @param[in] numberOfSmoothingIterations The number of smoothing iterations to apply to the input sample set
/// @param[in] meshRefinementParameters The mesh refinement parameters
/// @returns Error code
MKERNEL_API int mkernel_mesh2d_refine_ridges_based_on_gridded_samples(int meshKernelId,
const GriddedSamples& griddedSamples,
double relativeSearchRadius,
int minimumNumSamples,
int numberOfSmoothingIterations,
const meshkernel::MeshRefinementParameters& meshRefinementParameters);

/// @brief Remove any disconnected regions from a mesh2d.
///
/// The assumption is that the main region of interest has the largest number of elements.
Expand Down
69 changes: 68 additions & 1 deletion libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include <MeshKernel/Splines.hpp>

#include <MeshKernelApi/GeometryList.hpp>
#include <MeshKernelApi/GriddedSamples.hpp>
#include <MeshKernelApi/Mesh1D.hpp>
#include <MeshKernelApi/Mesh2D.hpp>

Expand Down Expand Up @@ -217,6 +218,72 @@ namespace meshkernelapi
}
}

/// @brief Computes the samples represented in gridded data in a vector of samples
/// @param[in] griddedSamples The gridded data to convert
/// @returns The converted vector of samples
template <meshkernel::InterpolatableType T>
static std::vector<meshkernel::Sample> ComputeGriddedDataSamples(const GriddedSamples& griddedSamples)
{
std::vector<meshkernel::Sample> result;
meshkernel::Point origin{griddedSamples.x_origin, griddedSamples.y_origin};
const auto numSamples = static_cast<size_t>(griddedSamples.num_x * griddedSamples.num_y);
result.resize(numSamples);
const T* valuePtr = static_cast<T*>(griddedSamples.values);
if (griddedSamples.x_coordinates == nullptr || griddedSamples.y_coordinates == nullptr)
{
meshkernel::UInt index = 0;

for (int j = 0; j < griddedSamples.num_x; ++j)
{
for (int i = griddedSamples.num_y - 1; i >= 0; --i)
{
const auto griddedIndex = griddedSamples.num_x * i + j;
result[index].x = origin.x + j * griddedSamples.cell_size;
result[index].y = origin.y + i * griddedSamples.cell_size;
result[index].value = static_cast<double>(valuePtr[griddedIndex]);
index++;
}
}
return result;
}

meshkernel::UInt index = 0;
for (int j = 0; j < griddedSamples.num_x; ++j)
{
for (int i = griddedSamples.num_y - 1; i >= 0; --i)
{
const auto griddedIndex = griddedSamples.num_x * i + j;
result[index].x = origin.x + griddedSamples.x_coordinates[griddedIndex];
result[index].y = origin.y + griddedSamples.y_coordinates[griddedIndex];
result[index].value = static_cast<double>(valuePtr[griddedIndex]);
index++;
}
}
return result;
}

/// @brief Converts the samples represented in gridded data in a vector of samples
/// @param[in] griddedSamples The gridded data to convert
/// @returns The converted vector of samples
static std::vector<meshkernel::Sample> ConvertGriddedData(const GriddedSamples& griddedSamples)
{
std::vector<meshkernel::Sample> result;
if (griddedSamples.num_x <= 0 || griddedSamples.num_y <= 0)
{
return result;
}

if (griddedSamples.value_type == static_cast<int>(meshkernel::InterpolationValues::shortType))
{
return ComputeGriddedDataSamples<short>(griddedSamples);
}
if (griddedSamples.value_type == static_cast<int>(meshkernel::InterpolationValues::floatType))
{
return ComputeGriddedDataSamples<float>(griddedSamples);
}
throw meshkernel::MeshKernelError("The value type for the gridded data samples is invalid.");
}

/// @brief Sets splines from a geometry list
/// @param[in] geometryListIn The input geometry list
/// @param[out] spline The spline which will be set
Expand Down Expand Up @@ -405,7 +472,7 @@ namespace meshkernelapi
const GriddedSamples& griddedSamples)
{
meshkernel::Point origin{griddedSamples.x_origin, griddedSamples.y_origin};
if (griddedSamples.x_coordinates == nullptr && griddedSamples.y_coordinates == nullptr)
if (griddedSamples.x_coordinates == nullptr || griddedSamples.y_coordinates == nullptr)
{
return std::make_unique<meshkernel::BilinearInterpolationOnGriddedSamples<T>>(mesh2d,
griddedSamples.num_x,
Expand Down
Loading

0 comments on commit ce074dc

Please sign in to comment.