diff --git a/src/data/store/ChunkedVolumeStore.cpp b/src/data/store/ChunkedVolumeStore.cpp index 12a55f2..76a651f 100644 --- a/src/data/store/ChunkedVolumeStore.cpp +++ b/src/data/store/ChunkedVolumeStore.cpp @@ -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 +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(sum) / cnt); + return static_cast(avg); +} + +// 降采样后某轴维度:ceil(n/2),至少 1(buildPyramid 与流式共用)。 +int halfDim(int n) { return std::max(1, (n + 1) / 2); } + // 从体中拷出一块的 int16(块内 i 最快、k 最慢,与体一致)。 std::vector 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(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(sum) / cnt); - dst.vox[idxAt(dst.nx, dst.ny, i, j, k)] = - static_cast(avg); - } - } - } - } + 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)] = + 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 rebuilt) { levels_ = std::move(rebuilt); levelCount_ = static_cast(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 rebuilt; + rebuilt.reserve(static_cast(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(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> cache( + static_cast(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(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& blk = + cache[(static_cast(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(lk) * lbh + lj) * lbw + li]; + }; + + // 逐 dst 体素降采样(共用核)。 + std::vector blk(static_cast(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, diff --git a/src/data/store/ChunkedVolumeStore.hpp b/src/data/store/ChunkedVolumeStore.hpp index 9cc6581..d476d81 100644 --- a/src/data/store/ChunkedVolumeStore.hpp +++ b/src/data/store/ChunkedVolumeStore.hpp @@ -62,6 +62,13 @@ class ChunkedVolumeStore { // data_L.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 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 rebuilt); + std::string dir_; StoreMeta meta_; int bricksX_ = 0, bricksY_ = 0, bricksZ_ = 0; // = level 0 块数(兼容字段) diff --git a/tests/data/store/test_pyramid.cpp b/tests/data/store/test_pyramid.cpp index cc44613..d7a9704 100644 --- a/tests/data/store/test_pyramid.cpp +++ b/tests/data/store/test_pyramid.cpp @@ -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(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();