diff --git a/tools/gpr_poc/main.cpp b/tools/gpr_poc/main.cpp index b7ef319..61c8f43 100644 --- a/tools/gpr_poc/main.cpp +++ b/tools/gpr_poc/main.cpp @@ -1280,6 +1280,18 @@ struct GradStats { }; GradStats sampleGradientMagnitude(vtkImageData* img); +// 标量分位标定(P3 修可见性核心):扫该体实际体素值分布,取 2%/98% 分位作色阶/ +// 不透明度的物理端点,裁掉离群。处理后体值多集中在 ±窄带、少量离群到 ±9000,若按 +// 全量化域(meta.vminPhys/vmaxPhys=±9249)映射 → 窄带信号落近透明区 → 整体近黑。 +// 返回物理单位的 {lo, hi}(已按 quant.toPhys 反算)。前置声明,实现在 polish 段。 +struct ScalarPercentiles { + double lo = 0.0, hi = 0.0; // 物理单位(2% / 98% 分位) + std::size_t samples = 0; +}; +ScalarPercentiles sampleScalarPercentiles(vtkImageData* img, + const geopro::core::Quant& q, + double pLo, double pHi); + geopro::core::ColorScale makeSeismicColorScale(double vmin, double vmax) { geopro::core::ColorScale cs; const double span = (vmax > vmin) ? (vmax - vmin) : 1.0; @@ -3059,6 +3071,10 @@ struct ViewState { vtkRenderWindow* rw = nullptr; double exagg = 8.0; int lastLevel = -1; + std::string dir; // store 目录:首帧直读 level0 局部段(类 gallery),绕开 LOD 选粗层 + // 预建的首帧高清段(level0 沿线中段):cmdView 已为分位标定建好,直接复用喂 mapper, + // 避免在 viewSetupDefaultFrame 内重复读盘。空则该函数再按 dir 直读。 + vtkSmartPointer seedSegImg; // 持有当前高清单图引用,避免被释放(mapper 仅持裸指针)。 vtkSmartPointer currentImg; // 回调防重入:回调内部会 Render(),若 Render 又触发观察者回调会无限递归。 @@ -3200,41 +3216,83 @@ std::size_t viewSetupDefaultFrame(ViewState* st, vtkRenderer* ren) { const int totBricksX = (m.nx + brick - 1) / brick; const int localBx = std::min(kViewDefaultLocalBricks, totBricksX); const int bx0 = std::max(0, totBricksX / 2 - localBx / 2); // 沿线中段 - // 该局部段世界 X 范围(level0)。 - (void)bx0; - (void)localBx; - (void)totBricksX; - // 体(exagg 后)世界尺寸与中心 + 包围球半径。 - const double wx = std::max(1.0, m.nx * m.spacing[0]); + // 该局部段世界 X 范围(level0,brick 列 [bx0, bx0+localBx))。把首帧相机框到这一段 + // (而非整卷)是 P3 修复 #2 的关键:相机近观局部段 → C1 selectLod 选 level0 局部子区, + // C2 重组该段单图,framing 该段 → 14×796 截面 + 沿线一段充满视野(类厚 B-scan)。 + // 框整卷则 selectLod 选最粗层(整条 45305 细带)、看着空白。 + const double segX0 = m.origin[0] + bx0 * brick * m.spacing[0]; + const double segX1 = + m.origin[0] + std::min(m.nx, (bx0 + localBx) * brick) * m.spacing[0]; + // 段(exagg 后)世界尺寸与中心 + 包围球半径。X 取该段宽(非整卷),Y/Z 全幅(薄轴)。 + const double wx = std::max(1.0, segX1 - segX0); const double wy = std::max(1.0, m.ny * m.spacing[1] * st->exagg); const double wz = std::max(1.0, m.nz * m.spacing[2] * st->exagg); - const double cx = m.origin[0] + 0.5 * wx; + const double cx = 0.5 * (segX0 + segX1); const double cy = m.origin[1] + 0.5 * wy; const double cz = m.origin[2] + 0.5 * wz; const double radius = 0.5 * std::sqrt(wx * wx + wy * wy + wz * wz); - // 相机从 +X 看体中心,距离 = 半径 / tan(半视角) × 余量,确保整体落入视锥(C1 - // selectLod 不会判 empty),由视距/分辨率自然选层;近观靠交互再拉近切细层。 + // 相机从 +Y 看段中心(看进【X-Z 宽面=B-scan 墙】),距离 = 半径/tan(半视角)×余量。 + // 段几何 = X≈12.6m 宽 × Y≈1.5m 薄(跨通道) × Z 深(exagg 后高)。exagg 只夸张深度(Z), + // Y 仍真实极薄 → 若从 +X(沿线)看只见薄前缘、近空。改从 +Y 俯看宽 X-Z 面: GPR 水平 + // 分层沿 X 铺开、随 Z 叠层,这一面才读得出内部结构。整条概览靠用户滚轮拉远(细带几何必然)。 st->cam = ren->GetActiveCamera(); const double fovY = st->cam->GetViewAngle(); const double halfAngle = 0.5 * fovY * 3.14159265358979 / 180.0; const double tanH = std::max(1e-3, std::tan(halfAngle)); const double dist = radius / tanH * 1.4; // 1.4:留余量含 aspect/边缘 st->cam->SetFocalPoint(cx, cy, cz); - st->cam->SetPosition(cx + dist, cy, cz); - st->cam->SetViewUp(0, 0, 1); + st->cam->SetPosition(cx, cy + dist, cz); // +Y 视点 → 正对 X-Z 宽面 + st->cam->SetViewUp(0, 0, 1); // Z 朝上(深度向下) ren->ResetCameraClippingRange(); - // 源选层选区 + 重组单图喂 mapper。初始化场景需保证拿到首图 → 阻塞轮询到就绪。 - const std::size_t blocks = viewRefreshBlocking(st); + // 首帧高清段直读(P3 修复 #2 核心):异步 LOD 源在「框一段」的视距下仍会选最粗层 + // (整条 45305 细带 → 看着近黑),不可取。改为【直接从 store 读 level0 沿线中段子体】 + // (与 gallery 的 buildLocalLevel0Image 同一直读路径,非 LOD 算法),喂高清 mapper — + // 保证首帧就是「全分辨率一段」的清晰块体。后续交互仍由异步源接管(用户拉远/拖动按 + // 视距正常选层)。退化(读不到段)再回退异步阻塞刷新。 + std::size_t blocks = 0; + { + vtkSmartPointer locImg = st->seedSegImg; // cmdView 预建,优先复用 + if (locImg == nullptr && !st->dir.empty()) { // 退化:按 dir 直读 + geopro::data::ChunkedVolumeStore store(st->dir); + locImg = buildLocalLevel0Image(store, m, bx0, localBx); + } + if (locImg != nullptr) { + st->currentImg = locImg; + st->lastLevel = 0; + st->mapper->SetInputData(locImg); + st->mapper->Update(); + viewSyncBaseCropping(st); // 底图按该段挖空,无缝叠加 + blocks = 1; + } + } + if (blocks == 0) blocks = viewRefreshBlocking(st); // 退化:回退异步源 - // 框住局部段:用无参 ResetCamera(按 actor 的【已 SetScale(1,exagg,exagg)】缩放 - // 后包围盒框),相机角度沿用能看出结构的 Elevation/Azimuth,再 Zoom 拉近填满画面。 - ren->ResetCamera(); + // 框住【局部段】:无参 ResetCamera 会按场景全部 actor(含常驻整卷底图,45305 长)的 + // 包围盒框 → 整条细带、截面填不满 → 看着空白。改为只框高清段(currentImg=level0 沿 + // 线中段子体)的包围盒,使 14×796 截面+沿线一段充满视野(类厚 B-scan);整条概览靠用户 + // 滚轮拉远(细带是 1:34 几何必然,非 bug)。Z 轴按 actor 的 SetScale(1,1,exagg) 同步夸张。 + if (st->currentImg != nullptr) { + double b[6]; + st->currentImg->GetBounds(b); // 高清段模型坐标盒(含绝对 X origin) + b[4] *= st->exagg; // Z 下界随深度夸张 + b[5] *= st->exagg; // Z 上界随深度夸张 + ren->ResetCamera(b); // 只框该段 → 段充满画面 + } else { + ren->ResetCamera(); // 退化(无高清段):回退全场景框 + } + // 取景角度:默认相机已置 +Y 正对 X-Z 宽面。var4 的 El45/Az30 是为 gallery 的 +X / + // Y 也夸张几何调的,套到这里(+Y 视、Y 真实极薄)会让宽面强烈斜退、大片黑。改用为本 + // +Y B-scan 几何调的小角度(El/Az 各 ~15-18°)留 3D 立体感而宽面仍基本正对,Zoom 拉足 + // 让 X-Z 面填满。仅默认/交互/preview 取景,不动 var4 gallery 参数。 + constexpr double kDefaultFrameElevation = 18.0; // 轻俯,见顶面薄边显层叠 + constexpr double kDefaultFrameAzimuth = 15.0; // 轻偏,宽面仍基本正对 + constexpr double kDefaultFrameZoom = 1.1; // ResetCamera 已贴合该段,只略收边距 st->cam = ren->GetActiveCamera(); - st->cam->Elevation(kViewDefaultVariant.elevation); // var4 取景:El18 - st->cam->Azimuth(kViewDefaultVariant.azimuth); // var4 取景:Az22 - st->cam->Zoom(kViewDefaultVariant.zoom); // var4 取景:Zoom2.0 填满画面 + st->cam->Elevation(kDefaultFrameElevation); + st->cam->Azimuth(kDefaultFrameAzimuth); + st->cam->Zoom(kDefaultFrameZoom); ren->ResetCameraClippingRange(); return blocks; } @@ -3245,8 +3303,6 @@ int runGalleryVariant(const std::string& dir, const GalleryVariant& v, const int winW = 1280, winH = 800; geopro::data::ChunkedVolumeStore store(dir); const geopro::data::StoreMeta& m = store.meta(); - const double vmin = m.vminPhys, vmax = m.vmaxPhys; - const geopro::core::ColorScale cs = pickColor(v.color, vmin, vmax); auto rw = makeOffscreenWindow(winW, winH); vtkNew ren; @@ -3261,6 +3317,21 @@ int runGalleryVariant(const std::string& dir, const GalleryVariant& v, vtkSmartPointer locImg = buildLocalLevel0Image(store, m, bx0, localBx); + // P3 传函分位标定:色阶/不透明度端点按该局部段【实际值 2%/98% 分位】裁离群。 + double vmin = m.vminPhys, vmax = m.vmaxPhys; + { + const ScalarPercentiles pc = + sampleScalarPercentiles(locImg.Get(), m.quant, 0.02, 0.98); + if (pc.samples > 0) { + vmin = pc.lo; + vmax = pc.hi; + std::cout << "[gallery " << v.name << "] 传函分位标定(样本 " << pc.samples + << "): 2%=" << vmin << " 98%=" << vmax << " (全域 [" + << m.vminPhys << ", " << m.vmaxPhys << "])\n"; + } + } + const geopro::core::ColorScale cs = pickColor(v.color, vmin, vmax); + GradStats gs; if (v.useGradientOpacity) { gs = sampleGradientMagnitude(locImg.Get()); @@ -3455,7 +3526,38 @@ int cmdView(int argc, char** argv) { source.setViewportHeight(winH); source.setAspect(static_cast(winW) / winH); const auto& m = source.meta(); - const double vmin = m.vminPhys, vmax = m.vmaxPhys; + // P3 首帧高清段(沿线中段 level0):直读 store 建该段单图——既作传函分位标定的基准 + // (该段实际值,对比比整卷底图更punchy),又供 viewSetupDefaultFrame 直接喂高清 mapper + // (绕开异步 LOD 在「框一段」视距下仍选最粗层的问题)。同 gallery 的 buildLocalLevel0Image + // 直读路径,非 LOD 算法。读不到则回退底图/全域。 + vtkSmartPointer segImg; + { + geopro::data::ChunkedVolumeStore store(dir); + const int totBx = store.bricksX(0); + const int localBx = std::min(kViewDefaultLocalBricks, totBx); + const int bx0 = std::max(0, totBx / 2 - localBx / 2); + segImg = buildLocalLevel0Image(store, m, bx0, localBx); + } + // P3 传函分位标定:色阶/不透明度端点按【实际值 2%/98% 分位】(物理域),裁离群。 + // 处理后体值多集中 ±窄带、少量离群 ±9000,按 meta 全域(±9249)映射 → 窄带近透明 → 全黑。 + // 基准用常驻底图(整卷代表):band 较窄 → 窄带信号映到更饱和色、可见度更高;且代表整线 + // 典型信号、不受单段局部高能离群影响。底图缺失再回退首帧段。退化(无样本)回退全域。 + double vmin = m.vminPhys, vmax = m.vmaxPhys; + vtkImageData* pcBasis = + source.baseImage() != nullptr ? source.baseImage() : segImg.Get(); + if (pcBasis != nullptr) { + const ScalarPercentiles pc = + sampleScalarPercentiles(pcBasis, m.quant, 0.02, 0.98); + if (pc.samples > 0) { + vmin = pc.lo; + vmax = pc.hi; + std::cout << "[view] 传函分位标定(" + << (source.baseImage() != nullptr ? "底图" : "局部段") + << ",样本 " << pc.samples << "): 2%=" << vmin << " 98%=" << vmax + << " (全域 [" << m.vminPhys << ", " << m.vmaxPhys + << "] → 裁离群)\n"; + } + } // 配色/不透明度包络取自 var4:seismic + V 形实体包络(floor/mid + opacity 作峰值)。 const geopro::core::ColorScale cs = pickColor(dv.color, vmin, vmax); // C4:默认变体(var4)开了梯度不透明度 → 从常驻底图实测梯度分布标定阈值。底图恒非空。 @@ -3541,6 +3643,8 @@ int cmdView(int argc, char** argv) { st.fpsText = fpsText.Get(); st.rw = rw.Get(); st.exagg = exagg; + st.dir = dir; // 供首帧直读 level0 局部段(绕开 LOD 选粗层) + st.seedSegImg = segImg; // cmdView 已为分位标定建好的首帧高清段,直接复用 // 相机初始定向(修复 1):默认框「局部段」而非整卷。整线横截面 1:34,框整卷 // 即便 exagg=8 也是一条隐形细带(看着空白);改为对准沿线中段一个 ~768 道窗口 @@ -3905,6 +4009,46 @@ GradStats sampleGradientMagnitude(vtkImageData* img) { return gs; } +// 标量分位标定实现(P3):步长抽样该 VTK_SHORT 体的非空体素值,排序取 pLo/pHi 分位, +// 用 quant.toPhys 反算成物理端点。失败/退化(无样本或 lo>=hi)返回 samples=0,由调用方 +// 回退到 meta 全量化域。仿 sampleGradientMagnitude 的 stride/blank 处理,自适应该体值域。 +ScalarPercentiles sampleScalarPercentiles(vtkImageData* img, + const geopro::core::Quant& q, + double pLo, double pHi) { + ScalarPercentiles out; + if (!img) return out; + int dims[3]; + img->GetDimensions(dims); + const vtkIdType npts = + static_cast(dims[0]) * dims[1] * dims[2]; + auto* arr = vtkShortArray::SafeDownCast(img->GetPointData()->GetScalars()); + if (!arr || npts <= 0) return out; + const std::int16_t blank = geopro::core::ScalarVolumeI16::kBlank; + std::vector vals; + vals.reserve(static_cast(npts) / 4 + 1); + // 步长抽样:~50 万样本足以代表分布(与梯度采样同档)。 + const vtkIdType stride = std::max(1, npts / 500000); + for (vtkIdType id = 0; id < npts; id += stride) { + const std::int16_t v = arr->GetValue(id); + if (v == blank) continue; + vals.push_back(v); + } + if (vals.empty()) return out; + std::sort(vals.begin(), vals.end()); + auto pick = [&](double p) { + const std::size_t idx = static_cast( + std::min(vals.size() - 1, p * (vals.size() - 1))); + return static_cast(vals[idx]); + }; + const double qLo = pick(pLo); + const double qHi = pick(pHi); + if (!(qLo < qHi)) return out; // 退化(常数体)→ 调用方回退全域 + out.lo = q.toPhys(static_cast(std::lround(qLo))); + out.hi = q.toPhys(static_cast(std::lround(qHi))); + out.samples = vals.size(); + return out; +} + // 标量不透明度:V 形包络(与 makeSolidVolumeProperty 同思路,floor/mid/max),三图共用, // 保证唯一变量是梯度不透明度 / 光照。 void setPolishScalarOpacity(vtkVolumeProperty* prop, const geopro::core::Quant& q,