feat(store): 流式金字塔 buildPyramidStreaming 逐块由盘上邻块降采样
level L 每块由 level L-1 的 ≤2×2×2 邻块(readBrick 从盘读)降采样 + 逐块增量写 data_L<L>.bin 得到,不重组整卷,任意时刻只持几个邻块+一个块,内存有界。 复用 buildPyramid 同一降采样核(downsampleVoxel 2×2×2 非 blank 平均 round、 全 blank→kBlank)、computeRange 与 finalizePyramidMeta 收尾(DRY),产出各级 dims/每块体素/min-max/hasRange/meta 与 buildPyramid 逐块一致。 测试:Pyramid.StreamingMatchesInRam 系列(128 整除、100/127 非整除奇数、 70×33×50 各向异性小 brick、全 blank 体)流式 vs 重组整卷逐块对拍;不破坏 buildPyramid/write/readBrick(store/streaming/pyramid 20 例全绿)。
This commit is contained in:
gaozheng
parent
4bb846cf07
commit
77cbe4a305
|
|
@ -88,6 +88,38 @@ json brickIndexJson(std::int64_t offset, std::int64_t clen, int bw, int bh,
|
|||
{"bd", bd}};
|
||||
}
|
||||
|
||||
// 2×2×2 平均降采样核（buildPyramid 与 buildPyramidStreaming 共用，DRY）：
|
||||
// 给定源维度 (snx,sny,snz) 与「读源体素 (si,sj,sk)」回调，算出 dst 体素 (i,j,k)
|
||||
// 的降采样值——2×2×2 邻域非 blank 平均（round），全 blank → kBlank。越界邻域跳过。
|
||||
template <typename SrcAt>
|
||||
std::int16_t downsampleVoxel(int i, int j, int k, int snx, int sny, int snz,
|
||||
SrcAt&& srcAt) {
|
||||
long sum = 0;
|
||||
int cnt = 0;
|
||||
for (int dk = 0; dk < 2; ++dk) {
|
||||
const int sk = 2 * k + dk;
|
||||
if (sk >= snz) continue;
|
||||
for (int dj = 0; dj < 2; ++dj) {
|
||||
const int sj = 2 * j + dj;
|
||||
if (sj >= sny) continue;
|
||||
for (int di = 0; di < 2; ++di) {
|
||||
const int si = 2 * i + di;
|
||||
if (si >= snx) continue;
|
||||
const std::int16_t v = srcAt(si, sj, sk);
|
||||
if (v == kBlank) continue;
|
||||
sum += v;
|
||||
++cnt;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (cnt == 0) return kBlank;
|
||||
const long avg = std::lround(static_cast<double>(sum) / cnt);
|
||||
return static_cast<std::int16_t>(avg);
|
||||
}
|
||||
|
||||
// 降采样后某轴维度：ceil(n/2)，至少 1（buildPyramid 与流式共用）。
|
||||
int halfDim(int n) { return std::max(1, (n + 1) / 2); }
|
||||
|
||||
// 从体中拷出一块的 int16（块内 i 最快、k 最慢，与体一致）。
|
||||
std::vector<std::int16_t> sliceBrick(const geopro::core::ScalarVolumeI16& vol,
|
||||
int bx, int by, int bz, int brick,
|
||||
|
|
@ -416,41 +448,19 @@ void ChunkedVolumeStore::buildPyramid(int levels) {
|
|||
for (int L = 1; L < target; ++L) {
|
||||
const DenseLevel& src = dense.back();
|
||||
DenseLevel dst;
|
||||
dst.nx = std::max(1, (src.nx + 1) / 2);
|
||||
dst.ny = std::max(1, (src.ny + 1) / 2);
|
||||
dst.nz = std::max(1, (src.nz + 1) / 2);
|
||||
dst.nx = halfDim(src.nx);
|
||||
dst.ny = halfDim(src.ny);
|
||||
dst.nz = halfDim(src.nz);
|
||||
dst.vox.assign(static_cast<std::size_t>(dst.nx) * dst.ny * dst.nz, kBlank);
|
||||
for (int k = 0; k < dst.nz; ++k) {
|
||||
for (int j = 0; j < dst.ny; ++j) {
|
||||
for (int i = 0; i < dst.nx; ++i) {
|
||||
// 2×2×2 邻域非 blank 平均（round），全 blank → kBlank。
|
||||
long sum = 0;
|
||||
int cnt = 0;
|
||||
for (int dk = 0; dk < 2; ++dk) {
|
||||
const int sk = 2 * k + dk;
|
||||
if (sk >= src.nz) continue;
|
||||
for (int dj = 0; dj < 2; ++dj) {
|
||||
const int sj = 2 * j + dj;
|
||||
if (sj >= src.ny) continue;
|
||||
for (int di = 0; di < 2; ++di) {
|
||||
const int si = 2 * i + di;
|
||||
if (si >= src.nx) continue;
|
||||
const std::int16_t v =
|
||||
src.vox[idxAt(src.nx, src.ny, si, sj, sk)];
|
||||
if (v == kBlank) continue;
|
||||
sum += v;
|
||||
++cnt;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (cnt > 0) {
|
||||
const long avg = std::lround(static_cast<double>(sum) / cnt);
|
||||
const auto srcAt = [&](int si, int sj, int sk) -> std::int16_t {
|
||||
return src.vox[idxAt(src.nx, src.ny, si, sj, sk)];
|
||||
};
|
||||
for (int k = 0; k < dst.nz; ++k)
|
||||
for (int j = 0; j < dst.ny; ++j)
|
||||
for (int i = 0; i < dst.nx; ++i)
|
||||
// 2×2×2 邻域非 blank 平均（round），全 blank → kBlank（共用核 DRY）。
|
||||
dst.vox[idxAt(dst.nx, dst.ny, i, j, k)] =
|
||||
static_cast<std::int16_t>(avg);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
downsampleVoxel(i, j, k, src.nx, src.ny, src.nz, srcAt);
|
||||
dense.push_back(std::move(dst));
|
||||
}
|
||||
|
||||
|
|
@ -513,12 +523,16 @@ void ChunkedVolumeStore::buildPyramid(int levels) {
|
|||
rebuilt.push_back(std::move(lv));
|
||||
}
|
||||
|
||||
finalizePyramidMeta(std::move(rebuilt));
|
||||
}
|
||||
|
||||
void ChunkedVolumeStore::finalizePyramidMeta(std::vector<Level> rebuilt) {
|
||||
levels_ = std::move(rebuilt);
|
||||
levelCount_ = static_cast<int>(levels_.size());
|
||||
// 兼容字段：level 0 索引（含新算的 min/max）回写 bricks_。
|
||||
bricks_ = levels_[0].bricks;
|
||||
|
||||
// 3) 重写 meta.json：保留所有原字段（含原 bricks），追加/覆盖 levels 数组。
|
||||
// 重写 meta.json：保留所有原字段（含原 bricks），追加/覆盖 levels 数组。
|
||||
const fs::path metaPath = fs::path(dir_) / kMetaFile;
|
||||
json meta;
|
||||
{
|
||||
|
|
@ -565,6 +579,134 @@ void ChunkedVolumeStore::buildPyramid(int levels) {
|
|||
if (!metaOut) throw std::runtime_error("ChunkedVolumeStore: meta.json write failed");
|
||||
}
|
||||
|
||||
void ChunkedVolumeStore::buildPyramidStreaming(int levels) {
|
||||
// 参数语义同 buildPyramid：levels = 最高级索引，总级数 = levels+1，levels<=0 仅
|
||||
// level 0。区别：各降采样级**逐块**由盘上 level L-1 的 ≤2×2×2 邻块降采样得到，
|
||||
// 不重组整卷——任意时刻只持几个 L-1 邻块 + 一个 L 块。
|
||||
const int target = levels < 0 ? 1 : levels + 1;
|
||||
const int brick = meta_.brick;
|
||||
|
||||
std::vector<Level> rebuilt;
|
||||
rebuilt.reserve(static_cast<std::size_t>(target));
|
||||
|
||||
// level 0：数据文件（data.bin）不变，复用现有 offset/clen，仅补每块 min/max。
|
||||
// 内存有界：逐块读盘算 range，不重组整卷。
|
||||
{
|
||||
Level lv0 = levels_[0]; // 拷 nx/ny/nz/块数/dataFile/索引
|
||||
for (int bz = 0; bz < lv0.bz; ++bz)
|
||||
for (int by = 0; by < lv0.by; ++by)
|
||||
for (int bx = 0; bx < lv0.bx; ++bx) {
|
||||
BrickEntry& be =
|
||||
lv0.bricks[static_cast<std::size_t>(brickIndexAt(lv0, bx, by, bz))];
|
||||
const auto rng = computeRange(readBrickFrom(levels_[0], bx, by, bz));
|
||||
be.vmin = rng.first;
|
||||
be.vmax = rng.second;
|
||||
be.hasRange = true;
|
||||
}
|
||||
rebuilt.push_back(std::move(lv0));
|
||||
}
|
||||
|
||||
// 降采样级 L=1..target-1：源 = rebuilt.back()（盘上 level L-1）。
|
||||
for (int L = 1; L < target; ++L) {
|
||||
const Level& src = rebuilt.back(); // 提供 dims/块数/dataFile，readBrickFrom 从盘读
|
||||
|
||||
Level lv;
|
||||
lv.nx = halfDim(src.nx);
|
||||
lv.ny = halfDim(src.ny);
|
||||
lv.nz = halfDim(src.nz);
|
||||
lv.bx = ceilDiv(lv.nx, brick);
|
||||
lv.by = ceilDiv(lv.ny, brick);
|
||||
lv.bz = ceilDiv(lv.nz, brick);
|
||||
lv.dataFile = levelDataFile(L);
|
||||
|
||||
std::ofstream data((fs::path(dir_) / lv.dataFile).string(),
|
||||
std::ios::binary | std::ios::trunc);
|
||||
if (!data)
|
||||
throw std::runtime_error("ChunkedVolumeStore: cannot open level data for write");
|
||||
std::int64_t offset = 0;
|
||||
|
||||
// 逐 L 块（固定顺序 bz 最慢、bx 最快，与 write/buildPyramid 一致）。
|
||||
for (int bz = 0; bz < lv.bz; ++bz) {
|
||||
for (int by = 0; by < lv.by; ++by) {
|
||||
for (int bx = 0; bx < lv.bx; ++bx) {
|
||||
const int bw = std::min(brick, lv.nx - bx * brick);
|
||||
const int bh = std::min(brick, lv.ny - by * brick);
|
||||
const int bd = std::min(brick, lv.nz - bz * brick);
|
||||
// 该 L 块 dst 体素全局坐标基点。
|
||||
const int i0 = bx * brick, j0 = by * brick, k0 = bz * brick;
|
||||
|
||||
// 覆盖的 L-1 体素子区域 = [2*i0, 2*(i0+bw)) 截到 src 维度。
|
||||
const int sLoI = 2 * i0, sLoJ = 2 * j0, sLoK = 2 * k0;
|
||||
const int sHiI = std::min(src.nx, 2 * (i0 + bw));
|
||||
const int sHiJ = std::min(src.ny, 2 * (j0 + bh));
|
||||
const int sHiK = std::min(src.nz, 2 * (k0 + bd));
|
||||
|
||||
// 落在的 L-1 邻块范围（≤2×2×2）。读入小缓存（按相对块索引存）。
|
||||
const int sbi0 = sLoI / brick, sbi1 = (sHiI - 1) / brick;
|
||||
const int sbj0 = sLoJ / brick, sbj1 = (sHiJ - 1) / brick;
|
||||
const int sbk0 = sLoK / brick, sbk1 = (sHiK - 1) / brick;
|
||||
const int nbi = sbi1 - sbi0 + 1;
|
||||
const int nbj = sbj1 - sbj0 + 1;
|
||||
const int nbk = sbk1 - sbk0 + 1;
|
||||
std::vector<std::vector<std::int16_t>> cache(
|
||||
static_cast<std::size_t>(nbi) * nbj * nbk);
|
||||
for (int kk = 0; kk < nbk; ++kk)
|
||||
for (int jj = 0; jj < nbj; ++jj)
|
||||
for (int ii = 0; ii < nbi; ++ii)
|
||||
cache[(static_cast<std::size_t>(kk) * nbj + jj) * nbi + ii] =
|
||||
readBrickFrom(src, sbi0 + ii, sbj0 + jj, sbk0 + kk);
|
||||
|
||||
// 从邻块缓存读源体素 (si,sj,sk)（全局 L-1 坐标）。
|
||||
const auto srcAt = [&](int si, int sj, int sk) -> std::int16_t {
|
||||
const int cbi = si / brick - sbi0;
|
||||
const int cbj = sj / brick - sbj0;
|
||||
const int cbk = sk / brick - sbk0;
|
||||
const std::vector<std::int16_t>& blk =
|
||||
cache[(static_cast<std::size_t>(cbk) * nbj + cbj) * nbi + cbi];
|
||||
// 块内 i 最快、k 最慢；块宽 = 该源块 extent。
|
||||
const int lbw = extent(src.nx, si / brick, brick);
|
||||
const int lbh = extent(src.ny, sj / brick, brick);
|
||||
const int li = si % brick, lj = sj % brick, lk = sk % brick;
|
||||
return blk[(static_cast<std::size_t>(lk) * lbh + lj) * lbw + li];
|
||||
};
|
||||
|
||||
// 逐 dst 体素降采样（共用核）。
|
||||
std::vector<std::int16_t> blk(static_cast<std::size_t>(bw) * bh * bd);
|
||||
std::size_t w = 0;
|
||||
for (int kk = 0; kk < bd; ++kk)
|
||||
for (int jj = 0; jj < bh; ++jj)
|
||||
for (int ii = 0; ii < bw; ++ii)
|
||||
blk[w++] = downsampleVoxel(i0 + ii, j0 + jj, k0 + kk, src.nx,
|
||||
src.ny, src.nz, srcAt);
|
||||
|
||||
const QByteArray compressed = compressBrick(blk);
|
||||
const std::int64_t clen = compressed.size();
|
||||
data.write(compressed.constData(), clen);
|
||||
if (!data)
|
||||
throw std::runtime_error("ChunkedVolumeStore: level data write failed");
|
||||
|
||||
const auto rng = computeRange(blk);
|
||||
BrickEntry be;
|
||||
be.offset = offset;
|
||||
be.compressedLen = clen;
|
||||
be.bw = bw;
|
||||
be.bh = bh;
|
||||
be.bd = bd;
|
||||
be.vmin = rng.first;
|
||||
be.vmax = rng.second;
|
||||
be.hasRange = true;
|
||||
lv.bricks.push_back(be);
|
||||
offset += clen;
|
||||
}
|
||||
}
|
||||
}
|
||||
data.close();
|
||||
rebuilt.push_back(std::move(lv));
|
||||
}
|
||||
|
||||
finalizePyramidMeta(std::move(rebuilt));
|
||||
}
|
||||
|
||||
// ----------------------- StreamingVolumeWriter -----------------------
|
||||
|
||||
StreamingVolumeWriter::StreamingVolumeWriter(const std::string& dir,
|
||||
|
|
|
|||
|
|
@ -62,6 +62,13 @@ class ChunkedVolumeStore {
|
|||
// data_L<level>.bin）。levels<=0 视为无金字塔（仅 level 0）。
|
||||
void buildPyramid(int levels);
|
||||
|
||||
// 流式构建金字塔：level L 的每块由 level L-1 的对应 ≤2×2×2 邻块（从盘 readBrick
|
||||
// 读）降采样 + 逐块增量写得到，**不重组整卷**——任意时刻只持几个 L-1 邻块 + 一个
|
||||
// L 块。产出（各级 dims/每块体素/min/max/hasRange/meta）与 buildPyramid(levels)
|
||||
// 逐块一致。降采样/blank/min-max 规则复用 buildPyramid 同一核（DRY）。
|
||||
// 需先有 level0 store（write 或流式建体产出）。levels<=0 视为仅 level 0。
|
||||
void buildPyramidStreaming(int levels);
|
||||
|
||||
// 总层数（含 level 0）；未建金字塔时 = 1。
|
||||
int levels() const { return levelCount_; }
|
||||
|
||||
|
|
@ -108,6 +115,12 @@ class ChunkedVolumeStore {
|
|||
std::vector<std::int16_t> readBrickFrom(const Level& lv, int bx, int by,
|
||||
int bz) const;
|
||||
|
||||
// 用 rebuilt 各级（已落盘各自数据文件、每块带 offset/clen/min/max/hasRange）
|
||||
// 替换 levels_，回写兼容字段 bricks_/levelCount_，并重写 meta.json 的 levels
|
||||
// 数组 + level0 兼容 bricks 的 min/max。buildPyramid 与 buildPyramidStreaming
|
||||
// 共用此收尾，确保两者 meta 结构逐字段一致（DRY）。
|
||||
void finalizePyramidMeta(std::vector<Level> rebuilt);
|
||||
|
||||
std::string dir_;
|
||||
StoreMeta meta_;
|
||||
int bricksX_ = 0, bricksY_ = 0, bricksZ_ = 0; // = level 0 块数（兼容字段）
|
||||
|
|
|
|||
|
|
@ -144,6 +144,121 @@ TEST(Pyramid, LegacyRealAllZeroBrickRangeIsZeroZero) {
|
|||
std::filesystem::remove_all(dir);
|
||||
}
|
||||
|
||||
// ----------------------- 流式金字塔对拍 -----------------------
|
||||
|
||||
namespace {
|
||||
// 体素随 (i,j,k) 变化，确保逐块对拍真正区分块内容（非全相同），含 blank 散点。
|
||||
geopro::core::BuiltI16 makeVaried(int nx, int ny, int nz) {
|
||||
geopro::core::BuiltI16 b;
|
||||
b.vol = geopro::core::ScalarVolumeI16(nx, ny, nz);
|
||||
for (int k = 0; k < nz; ++k)
|
||||
for (int j = 0; j < ny; ++j)
|
||||
for (int i = 0; i < nx; ++i) {
|
||||
const int v = (i * 131 + j * 17 + k * 7) % 251;
|
||||
// 散点 blank，验降采样 blank 混合（部分 blank 取非 blank 均值）。
|
||||
b.vol.at(i, j, k) =
|
||||
((i + j + k) % 13 == 0)
|
||||
? geopro::core::ScalarVolumeI16::kBlank
|
||||
: static_cast<std::int16_t>(v);
|
||||
}
|
||||
b.quant = {0.5, -3.0};
|
||||
b.origin = {{10.0, 20.0, 30.0}};
|
||||
b.spacing = {{2.0, 3.0, 4.0}};
|
||||
b.vminPhys = -3.0;
|
||||
b.vmaxPhys = 122.0;
|
||||
return b;
|
||||
}
|
||||
|
||||
// 流式 vs 重组整卷逐块对拍（dims + 每块体素 + min/max + hasRange + levels 数）。
|
||||
void expectStreamingMatchesInRam(const geopro::core::BuiltI16& b, int brick,
|
||||
int levels) {
|
||||
const auto dirA =
|
||||
(std::filesystem::temp_directory_path() / "gpr_pyr_match_A").string();
|
||||
const auto dirB =
|
||||
(std::filesystem::temp_directory_path() / "gpr_pyr_match_B").string();
|
||||
std::filesystem::remove_all(dirA);
|
||||
std::filesystem::remove_all(dirB);
|
||||
|
||||
ChunkedVolumeStore::write(dirA, b, brick);
|
||||
{
|
||||
ChunkedVolumeStore a(dirA);
|
||||
a.buildPyramid(levels); // 金标准（重组整卷）
|
||||
}
|
||||
ChunkedVolumeStore::write(dirB, b, brick);
|
||||
{
|
||||
ChunkedVolumeStore bb(dirB);
|
||||
bb.buildPyramidStreaming(levels); // 流式
|
||||
}
|
||||
|
||||
ChunkedVolumeStore A(dirA), B(dirB);
|
||||
ASSERT_EQ(A.levels(), B.levels());
|
||||
for (int L = 0; L < A.levels(); ++L) {
|
||||
int anx, any, anz, bnx, bny, bnz;
|
||||
A.dims(L, anx, any, anz);
|
||||
B.dims(L, bnx, bny, bnz);
|
||||
EXPECT_EQ(anx, bnx) << "level " << L << " nx";
|
||||
EXPECT_EQ(any, bny) << "level " << L << " ny";
|
||||
EXPECT_EQ(anz, bnz) << "level " << L << " nz";
|
||||
ASSERT_EQ(A.bricksX(L), B.bricksX(L));
|
||||
ASSERT_EQ(A.bricksY(L), B.bricksY(L));
|
||||
ASSERT_EQ(A.bricksZ(L), B.bricksZ(L));
|
||||
for (int bz = 0; bz < A.bricksZ(L); ++bz)
|
||||
for (int by = 0; by < A.bricksY(L); ++by)
|
||||
for (int bx = 0; bx < A.bricksX(L); ++bx) {
|
||||
EXPECT_EQ(A.readBrick(L, bx, by, bz), B.readBrick(L, bx, by, bz))
|
||||
<< "voxels mismatch L=" << L << " (" << bx << "," << by << ","
|
||||
<< bz << ")";
|
||||
EXPECT_EQ(A.brickRange(L, bx, by, bz), B.brickRange(L, bx, by, bz))
|
||||
<< "range mismatch L=" << L << " (" << bx << "," << by << ","
|
||||
<< bz << ")";
|
||||
}
|
||||
}
|
||||
std::filesystem::remove_all(dirA);
|
||||
std::filesystem::remove_all(dirB);
|
||||
}
|
||||
} // namespace
|
||||
|
||||
// 整除维度（128，brick64 → 2×2×2 满邻块）：流式与重组整卷逐块一致。
|
||||
TEST(Pyramid, StreamingMatchesInRam) {
|
||||
expectStreamingMatchesInRam(makeVaried(128, 128, 128), 64, 2);
|
||||
}
|
||||
|
||||
// 非整除/奇数维度（100、127），验边缘块 + 奇数维降采样一致。
|
||||
TEST(Pyramid, StreamingMatchesInRamNonDivisible) {
|
||||
expectStreamingMatchesInRam(makeVaried(100, 100, 100), 64, 2);
|
||||
expectStreamingMatchesInRam(makeVaried(127, 127, 127), 64, 3);
|
||||
}
|
||||
|
||||
// 非立方 + 小 brick（多级、边缘块更碎）：流式与重组整卷逐块一致。
|
||||
TEST(Pyramid, StreamingMatchesInRamAnisotropicSmallBrick) {
|
||||
expectStreamingMatchesInRam(makeVaried(70, 33, 50), 32, 3);
|
||||
}
|
||||
|
||||
// 全 blank 体流式降采样：各级仍全 blank（min/max=kBlank），与重组整卷一致。
|
||||
TEST(Pyramid, StreamingMatchesInRamAllBlank) {
|
||||
geopro::core::BuiltI16 b;
|
||||
b.vol = geopro::core::ScalarVolumeI16(70, 70, 70);
|
||||
for (auto& v : b.vol.data()) v = geopro::core::ScalarVolumeI16::kBlank;
|
||||
b.quant = {1.0, 0.0};
|
||||
b.origin = {{0, 0, 0}};
|
||||
b.spacing = {{1, 1, 1}};
|
||||
b.vminPhys = 0;
|
||||
b.vmaxPhys = 0;
|
||||
expectStreamingMatchesInRam(b, 32, 2);
|
||||
}
|
||||
|
||||
// 流式不破坏 level0：兼容重载等价。
|
||||
TEST(Pyramid, StreamingLevel0ReadCompatUnchanged) {
|
||||
auto dir =
|
||||
(std::filesystem::temp_directory_path() / "gpr_pyr_s_compat").string();
|
||||
std::filesystem::remove_all(dir);
|
||||
ChunkedVolumeStore::write(dir, makeRamp(64), 32);
|
||||
ChunkedVolumeStore s(dir);
|
||||
s.buildPyramidStreaming(1);
|
||||
EXPECT_EQ(s.readBrick(0, 0, 0), s.readBrick(0, 0, 0, 0));
|
||||
std::filesystem::remove_all(dir);
|
||||
}
|
||||
|
||||
// 老 store(未 buildPyramid):levels()==1;brickRange(0,...) 仍可惰性算。
|
||||
TEST(Pyramid, LegacyStoreNoPyramidLevelsIsOne) {
|
||||
auto dir = (std::filesystem::temp_directory_path() / "gpr_pyr_legacy").string();
|
||||
|
|
|
|||
Loading…
Reference in New Issue