feat(core): IDW 插值器(IInterpolator->ScalarVolume, 含 maxDist 包络裁剪)

This commit is contained in:
gaozheng 2026-06-07 19:53:22 +08:00
parent 868c49ca2c
commit e5a48c5af7
6 changed files with 88 additions and 0 deletions

View File

@ -3,6 +3,7 @@ find_package(Eigen3 CONFIG REQUIRED)
add_library(geopro_core STATIC add_library(geopro_core STATIC
geo/LocalFrame.cpp geo/LocalFrame.cpp
model/ColorScale.cpp model/ColorScale.cpp
algo/IdwInterpolator.cpp
) )
target_include_directories(geopro_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) target_include_directories(geopro_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

View File

@ -0,0 +1,21 @@
#pragma once
#include <vector>
#include "model/Field.hpp"
namespace geopro::core {
struct PointSet { std::vector<double> 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

View File

@ -0,0 +1,31 @@
#include "algo/IdwInterpolator.hpp"
#include <cmath>
#include <limits>
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<double>::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<double>::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

View File

@ -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

View File

@ -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_local_frame.cpp)
target_sources(geopro_tests PRIVATE core/test_model.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_color_scale.cpp)
target_sources(geopro_tests PRIVATE core/test_idw.cpp)
target_link_libraries(geopro_tests PRIVATE geopro_core) target_link_libraries(geopro_tests PRIVATE geopro_core)
add_subdirectory(spike) # spike S3: banded contour add_subdirectory(spike) # spike S3: banded contour

26
tests/core/test_idw.cpp Normal file
View File

@ -0,0 +1,26 @@
#include <gtest/gtest.h>
#include <cmath>
#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)
}