258 lines
9.3 KiB
C++
258 lines
9.3 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;
|
||
}
|
||
|
||
// 写一条 L 形弯轨迹的 .gps：先沿北(纬增,lon 固定)走半程，再沿东(经增,lat 固定)走半程。
|
||
void writeLGps(const fs::path& path, double lat0, double lon0, double dLatHalf,
|
||
double dLonHalf, int pts) {
|
||
std::ofstream f(path);
|
||
const int half = pts / 2;
|
||
for (int i = 0; i < pts; ++i) {
|
||
double lat, lon;
|
||
if (i <= half) {
|
||
const double frac = half > 0 ? static_cast<double>(i) / half : 0.0;
|
||
lat = lat0 + dLatHalf * frac; // 北段：纬增
|
||
lon = lon0;
|
||
} else {
|
||
const double frac =
|
||
(pts - 1 - half) > 0
|
||
? static_cast<double>(i - half) / (pts - 1 - half)
|
||
: 0.0;
|
||
lat = lat0 + dLatHalf; // 拐角后纬固定
|
||
lon = lon0 + dLonHalf * 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";
|
||
}
|
||
}
|
||
|
||
GeoLineInput makeLLine(const fs::path& dir, const std::string& tag, double lat0,
|
||
double lon0, double dLatHalf, double dLonHalf, 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"));
|
||
writeLGps(dir / (tag + ".gps"), lat0, lon0, dLatHalf, dLonHalf, 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);
|
||
}
|
||
|
||
// L 形弯线：PCA 版把拐后段甩进大 Y'(虚高 ny、稀疏)，曲线版沿 X 拉直展开
|
||
// (ny 接近真实横向、填充更密)。断言曲线版 ny < PCA 版 ny 且填充率更高。
|
||
TEST(GeoVolumeBuilder, CurvilinearStraightensBendVsPca) {
|
||
const fs::path tmp = fs::temp_directory_path() / "geopro_geovol_curv_test";
|
||
std::error_code ec;
|
||
fs::remove_all(tmp, ec);
|
||
fs::create_directories(tmp);
|
||
|
||
// L 形：北 ~111m（0.001° 纬）+ 东 ~96m（0.001° 经）。
|
||
const double lat0 = 30.200, lon0 = 120.244;
|
||
std::vector<GeoLineInput> lines = {makeLLine(
|
||
tmp, "synL_001", lat0, lon0, /*dLatHalf=*/0.001, /*dLonHalf=*/0.001,
|
||
/*traces=*/120, /*val=*/150)};
|
||
|
||
GeoGridSpec spec{/*cellXY=*/0.5, /*cellZ=*/0.1};
|
||
|
||
const std::string storePca = (tmp / "pca").string();
|
||
GeoBuildResult rPca =
|
||
buildGeoVolume(lines, spec, storePca, /*pyramidLevels=*/0,
|
||
/*curvilinear=*/false);
|
||
|
||
const std::string storeCurv = (tmp / "curv").string();
|
||
GeoBuildResult rCurv =
|
||
buildGeoVolume(lines, spec, storeCurv, /*pyramidLevels=*/0,
|
||
/*curvilinear=*/true);
|
||
|
||
const double fillPca =
|
||
rPca.total > 0 ? static_cast<double>(rPca.filled) / rPca.total : 0.0;
|
||
const double fillCurv =
|
||
rCurv.total > 0 ? static_cast<double>(rCurv.filled) / rCurv.total : 0.0;
|
||
|
||
// 曲线版把弯路拉直：横向维度 ny 明显小于 PCA 版。
|
||
EXPECT_LT(rCurv.ny, rPca.ny);
|
||
// 填充更密。
|
||
EXPECT_GT(fillCurv, fillPca);
|
||
// 曲线版仍建出有效体。
|
||
EXPECT_GT(rCurv.filled, 0);
|
||
EXPECT_GT(rCurv.nx, 50);
|
||
|
||
fs::remove_all(tmp, ec);
|
||
}
|
||
|
||
// 距离权纯函数：距 cell 中心越近权越大，单调递减，中心处取最大 1。
|
||
TEST(GeoVolumeBuilder, DistanceWeightDecreasesWithDistance) {
|
||
EXPECT_DOUBLE_EQ(distanceWeight(0.0), 1.0); // 中心权最大
|
||
EXPECT_LT(distanceWeight(1.0), distanceWeight(0.5)); // 单调递减
|
||
EXPECT_LT(distanceWeight(2.0), distanceWeight(1.0));
|
||
EXPECT_GT(distanceWeight(0.5), 0.0);
|
||
}
|
||
|
||
// 加权均值偏向近点：同 cell 两样本，近点(dist=0.1)值 100、远点(dist=0.9)值 300。
|
||
// 等权均值 = 200；距离加权均值应明显 < 200（偏向近点的 100）。
|
||
TEST(GeoVolumeBuilder, WeightedMeanFavorsNearerPoint) {
|
||
const double vNear = 100.0, vFar = 300.0;
|
||
const double wNear = distanceWeight(0.1);
|
||
const double wFar = distanceWeight(0.9);
|
||
const double weighted =
|
||
(wNear * vNear + wFar * vFar) / (wNear + wFar);
|
||
const double equal = 0.5 * (vNear + vFar); // 200
|
||
EXPECT_LT(weighted, equal); // 偏离等权
|
||
EXPECT_LT(weighted - vNear, vFar - weighted); // 更靠近近点值
|
||
}
|