diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 1bcba77..e41bbdd 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -8,6 +8,7 @@ add_library(geopro_core STATIC model/ColorScale.cpp algo/IdwInterpolator.cpp algo/VolumeBuilder.cpp + algo/GprVolumeBuilder.cpp ) target_include_directories(geopro_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/algo/GprVolumeBuilder.cpp b/src/core/algo/GprVolumeBuilder.cpp new file mode 100644 index 0000000..b06ceda --- /dev/null +++ b/src/core/algo/GprVolumeBuilder.cpp @@ -0,0 +1,112 @@ +#include "algo/GprVolumeBuilder.hpp" + +#include +#include + +namespace geopro::core { + +namespace { + +// 把世界坐标轴值映射到最近整数索引；落在 [0, n-1] 外返回 -1（标记越界）。 +int nearestIndex(double world, double origin, double step, int n) { + if (n <= 0 || step == 0.0) return -1; + long idx = std::lround((world - origin) / step); + if (idx < 0 || idx >= n) return -1; + return static_cast(idx); +} + +} // namespace + +BuiltI16 buildGprVolume(const GprSurvey& s, const GridSpec& spec) { + BuiltI16 out; + out.origin = {spec.ox, spec.oy, spec.oz}; + out.spacing = {spec.dx, spec.dy, spec.dz}; + + // 1. 物理值域（扫 values，跳过 NaN）。 + double vmin = std::numeric_limits::infinity(); + double vmax = -std::numeric_limits::infinity(); + for (double v : s.values) { + if (std::isnan(v)) continue; + if (v < vmin) vmin = v; + if (v > vmax) vmax = v; + } + if (!(vmin <= vmax)) { // 无有效值：退化为 [0,0]。 + vmin = 0.0; + vmax = 0.0; + } + out.vminPhys = vmin; + out.vmaxPhys = vmax; + + // 2. 量化：scale=(vmax-vmin)/64000 把整个物理值域铺满 int16 的 ~64000 个码位 + // （-32000..+32000），故 offset 取值域中点使量化对称——既用满 64000 码位， + // 又两端各留 ~700 余量不撞 int16 边界(±32767)与 kBlank(INT16_MIN)。 + // vmax==vmin 时 scale=1。 + out.quant.scale = (vmax > vmin) ? (vmax - vmin) / 64000.0 : 1.0; + out.quant.offset = 0.5 * (vmin + vmax); + + // 3. 分配体（构造即填 0），origin/spacing 已设。 + out.vol = ScalarVolumeI16(spec.nx, spec.ny, spec.nz); + + const int nChan = static_cast(s.channelY.size()); + + // 4. 逐网格点落值。 + for (int gk = 0; gk < spec.nz; ++gk) { + const double worldZ = spec.oz + gk * spec.dz; + const int sIdx = nearestIndex(worldZ, s.z0, s.dz, s.samples); + + for (int gj = 0; gj < spec.ny; ++gj) { + const double worldY = spec.oy + gj * spec.dy; + + for (int gi = 0; gi < spec.nx; ++gi) { + const double worldX = spec.ox + gi * spec.dx; + const int tIdx = nearestIndex(worldX, s.x0, s.dx, s.ntraces); + + // X / Z 越界 → blank。 + if (tIdx < 0 || sIdx < 0 || nChan == 0) { + out.vol.at(gi, gj, gk) = ScalarVolumeI16::kBlank; + continue; + } + + // Y → 跨通道 1D 线性插值（channelY 升序）。 + double phys = 0.0; + bool blank = false; + if (worldY <= s.channelY.front()) { + // 边界外不外推：在 maxDist 内 clamp 到首通道，否则 blank。 + if (s.channelY.front() - worldY > spec.maxDist) { + blank = true; + } else { + phys = s.at(0, tIdx, sIdx); + } + } else if (worldY >= s.channelY.back()) { + if (worldY - s.channelY.back() > spec.maxDist) { + blank = true; + } else { + phys = s.at(nChan - 1, tIdx, sIdx); + } + } else { + // 定位 worldY 落在哪两相邻通道间，线性插值。 + int lo = 0; + while (lo + 1 < nChan && s.channelY[lo + 1] < worldY) ++lo; + const int hi = lo + 1; + const double yLo = s.channelY[lo]; + const double yHi = s.channelY[hi]; + const double span = yHi - yLo; + const double w = (span > 0.0) ? (worldY - yLo) / span : 0.0; + const double vLo = s.at(lo, tIdx, sIdx); + const double vHi = s.at(hi, tIdx, sIdx); + phys = vLo + (vHi - vLo) * w; + } + + if (blank || std::isnan(phys)) { + out.vol.at(gi, gj, gk) = ScalarVolumeI16::kBlank; + } else { + out.vol.at(gi, gj, gk) = out.quant.toQ(phys); + } + } + } + } + + return out; +} + +} // namespace geopro::core diff --git a/src/core/algo/GprVolumeBuilder.hpp b/src/core/algo/GprVolumeBuilder.hpp new file mode 100644 index 0000000..d4e2623 --- /dev/null +++ b/src/core/algo/GprVolumeBuilder.hpp @@ -0,0 +1,27 @@ +#ifndef GEOPRO_CORE_ALGO_GPRVOLUMEBUILDER_HPP +#define GEOPRO_CORE_ALGO_GPRVOLUMEBUILDER_HPP + +#include + +#include "algo/IInterpolator.hpp" // GridSpec +#include "model/GprSurvey.hpp" +#include "model/ScalarVolumeI16.hpp" // ScalarVolumeI16, Quant + +namespace geopro::core { + +// 结构化建体产物：int16 量化体 + 量化映射 + 几何（origin/spacing）+ 物理值域。 +struct BuiltI16 { + ScalarVolumeI16 vol{0, 0, 0}; + Quant quant; + std::array origin{{0, 0, 0}}; + std::array spacing{{0, 0, 0}}; + double vminPhys = 0, vmaxPhys = 0; +}; + +// 结构化建体：X/Z 直接落格（取最近道/采样）+ 仅 Y 向跨通道 1D 线性插值。 +// 超 X/Z 范围或 Y 越界且超 maxDist 的网格点 → ScalarVolumeI16::kBlank（不外推）。 +BuiltI16 buildGprVolume(const GprSurvey& s, const GridSpec& spec); + +} // namespace geopro::core + +#endif // GEOPRO_CORE_ALGO_GPRVOLUMEBUILDER_HPP diff --git a/src/core/model/GprSurvey.hpp b/src/core/model/GprSurvey.hpp new file mode 100644 index 0000000..c54f1e1 --- /dev/null +++ b/src/core/model/GprSurvey.hpp @@ -0,0 +1,34 @@ +#ifndef GEOPRO_CORE_MODEL_GPRSURVEY_HPP +#define GEOPRO_CORE_MODEL_GPRSURVEY_HPP + +#include +#include + +namespace geopro::core { + +// 规则化建体输入（GPR 三维体）。雷达数据沿测线(X)、深度(Z)规则密采样， +// 仅横向/跨通道(Y)稀疏。结构化建体据此做 X/Z 落格 + 仅 Y 向 1D 线性插值。 +// +// 放 core/model（geopro::core）以保持 geopro_core 自洽：buildGprVolume 在 core， +// 若 GprSurvey 落 io/gpr 则 core 需反向 include io，故置于 core。 +// 真实数据 → GprSurvey 的装配在后续 POC 台（Task 9）完成。 +struct GprSurvey { + int ntraces = 0; // 沿测线道数（X） + int samples = 0; // 每道采样数（Z 深度） + double x0 = 0, dx = 1; // 沿测线轴：第 t 道 X = x0 + t*dx + double z0 = 0, dz = 1; // 深度轴：第 s 采样 Z = z0 + s*dz + + std::vector channelY; // 各通道横向位置（Y，稀疏，需按 Y 升序） + + // 值：channelY.size() × ntraces × samples，布局 [(c*ntraces + t)*samples + s]。 + std::vector values; + + // 64 位索引访问：values[(size_t(c)*ntraces + t)*samples + s]。 + double at(int c, int t, int s) const { + return values[(static_cast(c) * ntraces + t) * samples + s]; + } +}; + +} // namespace geopro::core + +#endif // GEOPRO_CORE_MODEL_GPRSURVEY_HPP diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b4b28de..2a02e8f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -33,6 +33,8 @@ target_sources(geopro_tests PRIVATE core/test_crs_transform.cpp) target_sources(geopro_tests PRIVATE core/test_model_data.cpp) target_sources(geopro_tests PRIVATE core/test_geo_frame.cpp) target_sources(geopro_tests PRIVATE core/test_scalar_volume_i16.cpp) +# GprVolumeBuilder：结构化建体（X/Z 落格 + 仅 Y 向 1D 线性插值 → int16 量化体）。 +target_sources(geopro_tests PRIVATE core/test_gpr_volume_builder.cpp) target_link_libraries(geopro_tests PRIVATE geopro_core) target_sources(geopro_tests PRIVATE data/test_parsers.cpp) diff --git a/tests/core/test_gpr_volume_builder.cpp b/tests/core/test_gpr_volume_builder.cpp new file mode 100644 index 0000000..3ecf116 --- /dev/null +++ b/tests/core/test_gpr_volume_builder.cpp @@ -0,0 +1,63 @@ +#include "core/algo/GprVolumeBuilder.hpp" +#include "core/algo/IInterpolator.hpp" // GridSpec +#include + +using namespace geopro::core; + +// GprSurvey 放 core/model（geopro::core 命名空间），保持 core 自洽，避免 core->io 反向依赖。 +TEST(GprVolumeBuilder, InterpolatesAcrossChannelsOnly) { + geopro::core::GprSurvey s; + s.ntraces = 1; + s.samples = 1; + s.x0 = 0; + s.dx = 1; + s.z0 = 0; + s.dz = 1; + s.channelY = {0.0, 1.0}; // 两通道：Y=0 值0、Y=1 值100 + s.values = {0.0, 100.0}; // [(c*1+0)*1+0] => c=0->0, c=1->100 + GridSpec spec{}; + spec.nx = 1; + spec.ny = 3; + spec.nz = 1; + spec.ox = 0; + spec.oy = 0; + spec.oz = 0; + spec.dx = 1; + spec.dy = 0.5; + spec.dz = 1; + spec.power = 2; + spec.maxDist = 9.9; + auto b = buildGprVolume(s, spec); + EXPECT_EQ(b.vol.nx(), 1); + EXPECT_EQ(b.vol.ny(), 3); + EXPECT_EQ(b.vol.nz(), 1); + EXPECT_NEAR(b.quant.toPhys(b.vol.at(0, 0, 0)), 0.0, 1.0); // Y=0 -> ch0 + EXPECT_NEAR(b.quant.toPhys(b.vol.at(0, 1, 0)), 50.0, 1.5); // Y=0.5 -> 中点≈50 + EXPECT_NEAR(b.quant.toPhys(b.vol.at(0, 2, 0)), 100.0, 1.0); // Y=1.0 -> ch1 +} + +TEST(GprVolumeBuilder, OutOfRangeBecomesBlank) { + geopro::core::GprSurvey s; + s.ntraces = 1; + s.samples = 1; + s.x0 = 0; + s.dx = 1; + s.z0 = 0; + s.dz = 1; + s.channelY = {0.0}; + s.values = {5.0}; + GridSpec spec{}; + spec.nx = 3; + spec.ny = 1; + spec.nz = 1; // X 超出仅 1 道的范围 + spec.ox = 0; + spec.oy = 0; + spec.oz = 0; + spec.dx = 10; + spec.dy = 1; + spec.dz = 1; + spec.power = 2; + spec.maxDist = 0.5; + auto b = buildGprVolume(s, spec); + EXPECT_EQ(b.vol.at(2, 0, 0), ScalarVolumeI16::kBlank); // 远端无覆盖->blank +}