diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 9a4a31c..a0d2a72 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -3,6 +3,7 @@ find_package(Eigen3 CONFIG REQUIRED) add_library(geopro_core STATIC geo/LocalFrame.cpp model/ColorScale.cpp + algo/IdwInterpolator.cpp ) target_include_directories(geopro_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/algo/IInterpolator.hpp b/src/core/algo/IInterpolator.hpp new file mode 100644 index 0000000..cea17c1 --- /dev/null +++ b/src/core/algo/IInterpolator.hpp @@ -0,0 +1,21 @@ +#pragma once +#include +#include "model/Field.hpp" +namespace geopro::core { + +struct PointSet { std::vector x, y, z, v; }; +struct GridSpec { + int nx, ny, nz; + double ox, oy, oz; // 原点 + double dx, dy, dz; // 步长 + double power; // IDW 幂 + double maxDist; // 超过则 blank(NaN),约束插值域(设计 §10) +}; + +class IInterpolator { +public: + virtual ~IInterpolator() = default; + virtual ScalarVolume interpolate(const PointSet& pts, const GridSpec& spec) const = 0; +}; + +} // namespace geopro::core diff --git a/src/core/algo/IdwInterpolator.cpp b/src/core/algo/IdwInterpolator.cpp new file mode 100644 index 0000000..c2e59a3 --- /dev/null +++ b/src/core/algo/IdwInterpolator.cpp @@ -0,0 +1,31 @@ +#include "algo/IdwInterpolator.hpp" +#include +#include +namespace geopro::core { + +ScalarVolume IdwInterpolator::interpolate(const PointSet& pts, const GridSpec& s) const { + ScalarVolume vol(s.nx, s.ny, s.nz); + const double nan = std::numeric_limits::quiet_NaN(); + const size_t n = pts.v.size(); + for (int k = 0; k < s.nz; ++k) + for (int j = 0; j < s.ny; ++j) + for (int i = 0; i < s.nx; ++i) { + const double gx = s.ox + i * s.dx, gy = s.oy + j * s.dy, gz = s.oz + k * s.dz; + double wsum = 0.0, vsum = 0.0, nearest = std::numeric_limits::max(); + bool hit = false; double hitVal = 0.0; + for (size_t p = 0; p < n; ++p) { + const double ddx = gx - pts.x[p], ddy = gy - pts.y[p], ddz = gz - pts.z[p]; + const double d = std::sqrt(ddx * ddx + ddy * ddy + ddz * ddz); + if (d < nearest) nearest = d; + if (d < 1e-12) { hit = true; hitVal = pts.v[p]; break; } + const double w = 1.0 / std::pow(d, s.power); + wsum += w; vsum += w * pts.v[p]; + } + if (hit) vol.at(i, j, k) = hitVal; + else if (nearest > s.maxDist || wsum == 0.0) vol.at(i, j, k) = nan; + else vol.at(i, j, k) = vsum / wsum; + } + return vol; +} + +} // namespace geopro::core diff --git a/src/core/algo/IdwInterpolator.hpp b/src/core/algo/IdwInterpolator.hpp new file mode 100644 index 0000000..fe16723 --- /dev/null +++ b/src/core/algo/IdwInterpolator.hpp @@ -0,0 +1,8 @@ +#pragma once +#include "algo/IInterpolator.hpp" +namespace geopro::core { +class IdwInterpolator : public IInterpolator { +public: + ScalarVolume interpolate(const PointSet& pts, const GridSpec& spec) const override; +}; +} // namespace geopro::core diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2533525..4d7024a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -12,6 +12,7 @@ gtest_discover_tests(geopro_tests) target_sources(geopro_tests PRIVATE core/test_local_frame.cpp) target_sources(geopro_tests PRIVATE core/test_model.cpp) target_sources(geopro_tests PRIVATE core/test_color_scale.cpp) +target_sources(geopro_tests PRIVATE core/test_idw.cpp) target_link_libraries(geopro_tests PRIVATE geopro_core) add_subdirectory(spike) # spike S3: banded contour 渲染验证 diff --git a/tests/core/test_idw.cpp b/tests/core/test_idw.cpp new file mode 100644 index 0000000..4799a30 --- /dev/null +++ b/tests/core/test_idw.cpp @@ -0,0 +1,26 @@ +#include +#include +#include "algo/IdwInterpolator.hpp" +using namespace geopro::core; + +TEST(Idw, ReproducesSampleAtNode) { + PointSet pts; + pts.x = {0.0, 10.0}; pts.y = {0.0, 0.0}; pts.z = {0.0, 0.0}; pts.v = {100.0, 200.0}; + GridSpec spec{ /*nx*/3, /*ny*/1, /*nz*/1, /*ox*/0, /*oy*/0, /*oz*/0, /*dx*/5, /*dy*/1, /*dz*/1, + /*power*/2.0, /*maxDist*/1e9 }; + IdwInterpolator idw; + ScalarVolume vol = idw.interpolate(pts, spec); + EXPECT_NEAR(vol.at(0, 0, 0), 100.0, 1e-6); // x=0 命中点1 + EXPECT_NEAR(vol.at(2, 0, 0), 200.0, 1e-6); // x=10 命中点2 + EXPECT_GT(vol.at(1, 0, 0), 100.0); + EXPECT_LT(vol.at(1, 0, 0), 200.0); +} + +TEST(Idw, BlanksOutsideMaxDist) { + PointSet pts; pts.x = {0.0}; pts.y = {0.0}; pts.z = {0.0}; pts.v = {50.0}; + GridSpec spec{2,1,1, 0,0,0, 100,1,1, 2.0, /*maxDist*/1.0}; + IdwInterpolator idw; + ScalarVolume vol = idw.interpolate(pts, spec); + EXPECT_NEAR(vol.at(0,0,0), 50.0, 1e-6); // 命中 + EXPECT_TRUE(std::isnan(vol.at(1,0,0))); // 超距 -> NaN(blank) +}