155 lines
5.2 KiB
C++
155 lines
5.2 KiB
C++
#include "core/algo/GeoVolumeBuilder.hpp"
|
||
|
||
#include <gtest/gtest.h>
|
||
|
||
#include <cmath>
|
||
#include <cstdint>
|
||
#include <filesystem>
|
||
#include <fstream>
|
||
#include <iomanip>
|
||
#include <string>
|
||
#include <vector>
|
||
|
||
#include "data/store/ChunkedVolumeStore.hpp"
|
||
|
||
namespace fs = std::filesystem;
|
||
using namespace geopro::core;
|
||
|
||
// ---- principalAxisAngle 纯函数 ----
|
||
TEST(GeoVolumeBuilder, PcaFindsRoadDirection) {
|
||
// 点沿 +Y（北向）排布 → 主轴角约 ±90°(±pi/2)。
|
||
std::vector<double> xs, ys;
|
||
for (int i = 0; i < 20; ++i) {
|
||
xs.push_back(0.01 * i); // 微小横向噪声
|
||
ys.push_back(static_cast<double>(i));
|
||
}
|
||
const double ang = principalAxisAngle(xs, ys);
|
||
EXPECT_NEAR(std::abs(ang), 3.14159265358979323846 / 2.0, 0.05);
|
||
}
|
||
|
||
TEST(GeoVolumeBuilder, PcaAlongEastIsZero) {
|
||
std::vector<double> xs, ys;
|
||
for (int i = 0; i < 20; ++i) {
|
||
xs.push_back(static_cast<double>(i));
|
||
ys.push_back(0.01 * i);
|
||
}
|
||
EXPECT_NEAR(principalAxisAngle(xs, ys), 0.0, 0.05);
|
||
}
|
||
|
||
TEST(GeoVolumeBuilder, PcaDegenerateReturnsZero) {
|
||
EXPECT_DOUBLE_EQ(principalAxisAngle({}, {}), 0.0);
|
||
EXPECT_DOUBLE_EQ(principalAxisAngle({1.0}, {1.0}), 0.0);
|
||
}
|
||
|
||
// ---- 合成小线建体 ----
|
||
namespace {
|
||
|
||
// 写一通道 .iprb + .iprh：值恒为 fixedVal（便于校验重叠均值）。
|
||
void writeChannel(const fs::path& iprb, int samples, int traces,
|
||
std::int16_t fixedVal) {
|
||
fs::path iprh = fs::path(iprb).replace_extension(".iprh");
|
||
std::ofstream h(iprh);
|
||
h << "SAMPLES: " << samples << "\n";
|
||
h << "LAST TRACE: " << (traces - 1) << "\n";
|
||
h << "CHANNELS: 2\n";
|
||
h << "TIMEWINDOW: 100.0\n";
|
||
h << "SOIL VELOCITY: 100.0\n"; // m/µs → 1e8 m/s
|
||
h << "DISTANCE INTERVAL: 0.05\n";
|
||
h.close();
|
||
|
||
std::ofstream b(iprb, std::ios::binary);
|
||
for (int t = 0; t < traces; ++t)
|
||
for (int s = 0; s < samples; ++s)
|
||
b.write(reinterpret_cast<const char*>(&fixedVal), sizeof(fixedVal));
|
||
}
|
||
|
||
// 写一条南北直线的 .gps（lat 从 lat0 递增到 lat1，lon 固定）。
|
||
void writeGps(const fs::path& path, double lat0, double lat1, double lon,
|
||
int pts) {
|
||
std::ofstream f(path);
|
||
for (int i = 0; i < pts; ++i) {
|
||
const double frac = pts > 1 ? static_cast<double>(i) / (pts - 1) : 0.0;
|
||
const double lat = lat0 + (lat1 - lat0) * frac;
|
||
f << "2023-06-03\t00:00:00:000\t" << std::fixed
|
||
<< std::setprecision(10) << lat << "\tN\t" << lon << "\tE\t9.0\tM\t4\n";
|
||
}
|
||
}
|
||
|
||
// 写两通道 .ord（横偏 -0.5 / +0.5，末列=1 有效）。
|
||
void writeOrd(const fs::path& path) {
|
||
std::ofstream f(path);
|
||
f << "0 -0.500000 -1.5 1\n";
|
||
f << "1 0.500000 -1.5 1\n";
|
||
}
|
||
|
||
GeoLineInput makeLine(const fs::path& dir, const std::string& tag, double lat0,
|
||
double lat1, double lon, int traces, std::int16_t val) {
|
||
const int samples = 8;
|
||
writeChannel(dir / (tag + "_A01.iprb"), samples, traces, val);
|
||
writeChannel(dir / (tag + "_A02.iprb"), samples, traces, val);
|
||
writeOrd(dir / (tag + ".ord"));
|
||
writeGps(dir / (tag + ".gps"), lat0, lat1, lon, traces);
|
||
GeoLineInput in;
|
||
in.iprb = {(dir / (tag + "_A01.iprb")).string(),
|
||
(dir / (tag + "_A02.iprb")).string()};
|
||
in.ord = (dir / (tag + ".ord")).string();
|
||
in.gps = (dir / (tag + ".gps")).string();
|
||
return in;
|
||
}
|
||
|
||
} // namespace
|
||
|
||
TEST(GeoVolumeBuilder, BuildsSyntheticLinesOverlapAveraged) {
|
||
const fs::path tmp = fs::temp_directory_path() / "geopro_geovol_test";
|
||
std::error_code ec;
|
||
fs::remove_all(tmp, ec);
|
||
fs::create_directories(tmp);
|
||
|
||
// 两条同位置南北线（lat 30.200→30.201，~111m），值不同 → 重叠 cell 取均值。
|
||
const double lat0 = 30.200, lat1 = 30.201, lon = 120.244;
|
||
std::vector<GeoLineInput> lines = {
|
||
makeLine(tmp, "synA_001", lat0, lat1, lon, /*traces=*/40, /*val=*/100),
|
||
makeLine(tmp, "synB_002", lat0, lat1, lon, /*traces=*/40, /*val=*/300),
|
||
};
|
||
|
||
const std::string store = (tmp / "store").string();
|
||
GeoGridSpec spec{/*cellXY=*/0.5, /*cellZ=*/0.1};
|
||
GeoBuildResult r = buildGeoVolume(lines, spec, store, /*pyramidLevels=*/1);
|
||
|
||
// 维度合理：约 111m 长 → nx 在数百量级；横路窄(~1m 阵列) → ny 小；nz>1。
|
||
EXPECT_GT(r.nx, 50);
|
||
EXPECT_GE(r.ny, 1);
|
||
EXPECT_GT(r.nz, 1);
|
||
EXPECT_GT(r.filled, 0);
|
||
EXPECT_LE(r.filled, r.total);
|
||
|
||
// store 可读，维度一致。
|
||
geopro::data::ChunkedVolumeStore s(store);
|
||
EXPECT_EQ(s.meta().nx, r.nx);
|
||
EXPECT_EQ(s.meta().ny, r.ny);
|
||
EXPECT_EQ(s.meta().nz, r.nz);
|
||
EXPECT_EQ(s.levels(), 2); // level0 + 1
|
||
|
||
// 重叠均值：两线值 100/300，命中同 cell → 均值 200。扫所有块找非 blank 体素，
|
||
// 其物理值应接近 200（量化误差内）。
|
||
const auto& m = s.meta();
|
||
bool foundNonBlank = false;
|
||
double sampleVal = 0.0;
|
||
for (int bz = 0; bz < s.bricksZ() && !foundNonBlank; ++bz)
|
||
for (int by = 0; by < s.bricksY() && !foundNonBlank; ++by)
|
||
for (int bx = 0; bx < s.bricksX() && !foundNonBlank; ++bx) {
|
||
auto vox = s.readBrick(bx, by, bz);
|
||
for (std::int16_t q : vox) {
|
||
if (q != geopro::core::ScalarVolumeI16::kBlank) {
|
||
sampleVal = m.quant.toPhys(q);
|
||
foundNonBlank = true;
|
||
break;
|
||
}
|
||
}
|
||
}
|
||
ASSERT_TRUE(foundNonBlank);
|
||
EXPECT_NEAR(sampleVal, 200.0, 2.0);
|
||
|
||
fs::remove_all(tmp, ec);
|
||
}
|