feat(vtk): 12d 打磨探针-梯度不透明度+光照,出体内部对比3图

新增 gpr_poc polish 子命令:同一全分辨率 level0 局部段 + 斜穿俯视视角离屏渲三图对比体内部白雾能否靠打磨解决。梯度不透明度 piecewise 按实测梯度幅值分布(median/p90/p99)标定。三图唯一变量为梯度不透明度/光照:a 基线白雾、b 加梯度不透明度、c 加梯度+光照。各报结构像素/亮度/真实fps。结论:打磨消雾并让强梯度层界面浮出,但沿线长均匀段固有偏雾,梯度不透明度无法凭空长出层。
This commit is contained in:
gaozheng 2026-06-23 20:53:42 +08:00
parent fb175d6d3d
commit 0537e938b4
5 changed files with 383 additions and 1 deletions

View File

@ -0,0 +1,45 @@
# Task 12d-polish 报告:梯度不透明度 + 光照 打磨探针
## 状态
完成。真实离屏渲染、真实 fps无编造。三张对比图均通过双闸无 3D 纹理维度错 + 渲出高于背景像素)。
## 命令
`gpr_poc polish tmp\store_lod_001 --frames 90`(默认取景 El45/Az30/Zoom1.5「斜穿俯视」,视线从上方斜穿体内部而非只看端面)。
## 测试体
- 全分辨率 level0 局部段256 x 29 x 162沿线中段 4 brick 列 [345,349)/695垂向夸张 exagg=8放大薄 Y/Z 轴使截面可读)。
- 三图标量传函/配色/取景/夸张全相同,唯一变量是「梯度不透明度 / 光照」。
## 梯度幅值分布量化域中心差分545211 样本,按实测标定阈值)
median=5.32p75=20.1p90=196.2p99=9058.5max=21470。
梯度不透明度 piecewise按此分布标定非猜grad≤5.32→0.0、grad=196.2(p90)→0.5、grad≥9058.5(p99)→0.9。
即:占多数的低梯度均匀区透明,仅高梯度处(层界面)不透明。
标量不透明度峰值:基线 a=0.15与默认体绘制同档→白雾b/c 梯度门控压住均匀区后提到 0.6,让层界面净不透明度(标量×梯度)足够高、层面成实面。
## 三张对比图docs/superpowers/plans/poc-lod-shots/
| 图 | 路径 | 高于背景像素(>35) | 结构像素(>50) | 平均亮度(0-255) | fps |
|---|---|---|---|---|---|
| a 基线白雾 | polish-a-value.png | 219980 (21.5%) | 145874 (14.2%) | 20.07 | 160.8 |
| b +梯度不透明度 | polish-b-grad.png | 50358 (4.9%) | 36430 (3.6%) | 15.71 | 58.2 |
| c +梯度+光照 | polish-c-grad-shade.png | 25008 (2.4%) | 756 (0.07%) | 14.81 | 57.3 |
光照参数cShadeOnAmbient 0.3 / Diffuse 0.7 / Specular 0.2 / SpecularPower 10。
## 目视结论:内部层是否「浮」出来了?
**部分浮出来了,但不是全身——这正是层状数据的固有限制。**
- **a基线**:一根平滑均匀的灰蓝色长条,没有任何内部层次,只有端面隐约可辨——就是需求描述的「体中间均匀白雾、只端面有层次」。穿透均匀水平层积分成雾。
- **b+梯度不透明度)**:体的大部分(沿线中段那段均匀体)变透明、白雾消失,**端部/过渡区露出清晰的水平层状条纹**(层界面),底部另现一块淡蓝层状斑。证实:梯度不透明度确实把均匀积分雾抹掉、把层界面显出来了。
- **c+梯度+光照)**:在 b 基础上端部层条纹带上轻微立体明暗(层带有了明暗起伏的层次感),但 shading 整体压暗,可见区更少更暗。
**关键如实结论**:梯度不透明度 + 光照**能消除均匀白雾、并让「确有梯度突变的层界面」浮出成可读的层状条纹**——打磨方向有效。**但对这条道路 GPR 数据,强梯度集中在端部/过渡区;沿线的长段水平层因「沿测线方向看过去是均匀的」(梯度低)会整段变透明,而不会显出层。**所以打磨**改善了「层界面可见性」**,但**无法让整条体内部都「长出层」——长均匀段的内部偏雾/偏空是层状数据本身的固有属性,不是没打磨。** 想看长段内部层,应配合切片/正交截面,而非纯体绘制穿透积分。
## fps 代价
梯度不透明度需逐采样点算梯度fps 从基线 160 降到 ~58约 -64%),但仍远高于交互级 15fps**可接受**。光照c相对 b 近乎免费58→57。本机 RTX 3060 数;最低配未验。
## 提交自检
- 仅 `git add` tools/gpr_poc/main.cpp + 3 张 polish-*.png + 本报告;未 `git add -A`
- `git diff --cached --stat` 确认无 chart/scatter/quill/rangeslider/Dialog/FormK。
- 未改任何交互默认(探针性质,仅新增 polish 子命令与三个 polish 专用辅助函数)。

Binary file not shown.

After

Width:  |  Height:  |  Size: 99 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 55 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 45 KiB

View File

@ -3019,6 +3019,340 @@ int cmdView(int argc, char** argv) {
return 0;
}
// ============================================================================
// 12d-polish梯度不透明度 + 光照 打磨探针(验证"体内部白雾"能否靠打磨解决)
// ============================================================================
//
// 当前体绘制对道路 GPR 水平层数据,体中间是均匀白雾、只有端面有层次。本探针在同一
// 全分辨率(level0)局部段 + 同一「看进体内部」视角(斜穿俯视,视线穿过体内部而非只看
// 端面)下渲 3 张离屏对比图,验证:给体加【梯度不透明度】(均匀区透明、层界面显出) +
// 【光照/明暗】能否让内部层状结构"浮"出来:
// polish-a-value.png 基线:按数值的不透明度(V形包络),无梯度不透明度、无光照
// polish-b-grad.png + 梯度不透明度(SetGradientOpacity)
// polish-c-grad-shade.png + 梯度不透明度 + 光照(ShadeOn, Ambient/Diffuse/Specular)
//
// 梯度不透明度的 piecewise 按【实际梯度幅值分布】标定阈值(不靠猜):先在量化域逐体素
// 采样 6 邻居中心差分梯度幅值,取分位数(median / p90)作斜坡控制点。
// 量化域标量不透明度峰值。基线(a)用 0.15(与默认体绘制同档→均匀积分白雾);开了梯度
// 不透明度(b/c)后均匀区被梯度门压成透明,可放心把标量峰值提到 0.6,让【层界面】这类
// 高梯度处的净不透明度(标量×梯度)足够高、层面真正"浮"成实面,而非仍是淡影。
constexpr double kPolishMaxOpacityFog = 0.15; // a基线白雾
constexpr double kPolishMaxOpacityGrad = 0.6; // b/c梯度门控后可提高
// 在 VTK_SHORT 局部体上采样梯度幅值分布(量化域,中心差分),返回有序的若干分位数。
// 跳过 kBlank 体素及其邻居(空值不参与梯度)。返回 {median, p75, p90, p99, max}。
struct GradStats {
double median = 0, p75 = 0, p90 = 0, p99 = 0, mx = 0;
std::size_t samples = 0;
};
GradStats sampleGradientMagnitude(vtkImageData* img) {
int dims[3];
img->GetDimensions(dims);
const int nx = dims[0], ny = dims[1], nz = dims[2];
auto* arr = vtkShortArray::SafeDownCast(img->GetPointData()->GetScalars());
GradStats gs;
if (!arr || nx < 3 || ny < 3 || nz < 3) return gs;
const std::int16_t blank = geopro::core::ScalarVolumeI16::kBlank;
auto at = [&](int i, int j, int k) -> std::int16_t {
const vtkIdType id = (static_cast<vtkIdType>(k) * ny + j) * nx + i;
return arr->GetValue(id);
};
std::vector<double> mags;
mags.reserve(static_cast<std::size_t>(nx) * ny * nz / 8 + 1);
// 步长抽样:大体不必逐体素,间隔取样即可代表分布(≤~50万样本
const vtkIdType total = static_cast<vtkIdType>(nx - 2) * (ny - 2) * (nz - 2);
const int stride = static_cast<int>(std::max<vtkIdType>(1, total / 500000));
vtkIdType counter = 0;
for (int k = 1; k < nz - 1; ++k) {
for (int j = 1; j < ny - 1; ++j) {
for (int i = 1; i < nx - 1; ++i) {
if ((counter++ % stride) != 0) continue;
const std::int16_t c = at(i, j, k);
if (c == blank) continue;
const std::int16_t xm = at(i - 1, j, k), xp = at(i + 1, j, k);
const std::int16_t ym = at(i, j - 1, k), yp = at(i, j + 1, k);
const std::int16_t zm = at(i, j, k - 1), zp = at(i, j, k + 1);
if (xm == blank || xp == blank || ym == blank || yp == blank ||
zm == blank || zp == blank) {
continue;
}
const double gx = 0.5 * (xp - xm);
const double gy = 0.5 * (yp - ym);
const double gz = 0.5 * (zp - zm);
mags.push_back(std::sqrt(gx * gx + gy * gy + gz * gz));
}
}
}
if (mags.empty()) return gs;
std::sort(mags.begin(), mags.end());
auto q = [&](double p) {
const std::size_t idx = static_cast<std::size_t>(
std::min<double>(mags.size() - 1, p * (mags.size() - 1)));
return mags[idx];
};
gs.median = q(0.50);
gs.p75 = q(0.75);
gs.p90 = q(0.90);
gs.p99 = q(0.99);
gs.mx = mags.back();
gs.samples = mags.size();
return gs;
}
// 标量不透明度V 形包络(与 makeSolidVolumeProperty 同思路floor/mid/max三图共用
// 保证唯一变量是梯度不透明度 / 光照。
void setPolishScalarOpacity(vtkVolumeProperty* prop, const geopro::core::Quant& q,
double vminPhys, double vmaxPhys, double maxOpacity) {
if (vminPhys >= vmaxPhys) vmaxPhys = vminPhys + 1.0;
const double qminD = static_cast<double>(q.toQ(vminPhys));
const double qmaxD = static_cast<double>(q.toQ(vmaxPhys));
const double qmid = 0.5 * (qminD + qmaxD);
const double half = 0.5 * (qmaxD - qminD);
const double floorOp = 0.4 * maxOpacity; // 中段背景按峰值比例压低V 形)
vtkNew<vtkPiecewiseFunction> opacity;
opacity->AddPoint(
static_cast<double>(geopro::core::ScalarVolumeI16::kBlank), 0.0);
opacity->AddPoint(qminD, maxOpacity);
opacity->AddPoint(qmid - 0.55 * half, floorOp);
opacity->AddPoint(qmid, 0.2 * maxOpacity);
opacity->AddPoint(qmid + 0.55 * half, floorOp);
opacity->AddPoint(qmaxD, maxOpacity);
prop->SetScalarOpacity(opacity);
}
// 颜色传函量化域seismic 红白蓝,与其余探针同思路)。
void setPolishColor(vtkVolumeProperty* prop, const geopro::core::Quant& q,
const geopro::core::ColorScale& cs, double vminPhys,
double vmaxPhys) {
constexpr int kSamples = 64;
if (vminPhys >= vmaxPhys) vmaxPhys = vminPhys + 1.0;
const double qminD = static_cast<double>(q.toQ(vminPhys));
const double qmaxD = static_cast<double>(q.toQ(vmaxPhys));
vtkNew<vtkColorTransferFunction> color;
for (int t = 0; t < kSamples; ++t) {
const double qd = qminD + (qmaxD - qminD) * t / (kSamples - 1);
const auto qv = static_cast<std::int16_t>(std::lround(qd));
const auto c = cs.colorAt(q.toPhys(qv));
color->AddRGBPoint(qd, c.r / 255.0, c.g / 255.0, c.b / 255.0);
}
prop->SetColor(color);
}
// 共用:把局部体喂单 SmartVolumeMapper按一个「看进体内部」的相机取景渲一帧 + 存 PNG
// 报告 结构像素 / 平均亮度 / fps。mode: 0=value 基线 1=+grad 2=+grad+shade。
int renderPolishOne(vtkImageData* locImg, const geopro::core::Quant& quant,
const geopro::core::ColorScale& cs, double vmin, double vmax,
const GradStats& gs, int mode, double exagg,
const std::string& pngPath, int frames, double elevation,
double azimuth, double zoom) {
const int winW = 1280, winH = 800;
auto rw = makeOffscreenWindow(winW, winH);
vtkNew<vtkRenderer> ren;
ren->SetBackground(0.05, 0.05, 0.09);
rw->AddRenderer(ren);
vtkNew<vtkSmartVolumeMapper> mapper;
mapper->SetInputData(locImg);
mapper->SetRequestedRenderMode(vtkSmartVolumeMapper::GPURenderMode);
mapper->SetAutoAdjustSampleDistances(0);
mapper->SetInteractiveAdjustSampleDistances(0);
auto prop = vtkSmartPointer<vtkVolumeProperty>::New();
setPolishColor(prop, quant, cs, vmin, vmax);
const double scalarMax =
(mode == 0) ? kPolishMaxOpacityFog : kPolishMaxOpacityGrad;
setPolishScalarOpacity(prop, quant, vmin, vmax, scalarMax);
prop->SetInterpolationTypeToLinear();
// 梯度不透明度mode>=1梯度小(均匀区)→透明、梯度大(层界面)→不透明。
// 阈值按实测分布median 处仍接近 0压住均匀积分雾p90 升到 0.5p99 到 0.9。
if (mode >= 1) {
vtkNew<vtkPiecewiseFunction> grad;
grad->AddPoint(0.0, 0.0);
grad->AddPoint(std::max(1.0, gs.median), 0.0);
grad->AddPoint(std::max(2.0, gs.p90), 0.5);
grad->AddPoint(std::max(3.0, gs.p99), 0.9);
prop->SetGradientOpacity(grad);
}
// 光照mode>=2ShadeOn + Ambient/Diffuse/Specular让层界面带立体明暗。
if (mode >= 2) {
prop->ShadeOn();
prop->SetAmbient(0.3);
prop->SetDiffuse(0.7);
prop->SetSpecular(0.2);
prop->SetSpecularPower(10.0);
} else {
prop->ShadeOff();
}
auto volume = vtkSmartPointer<vtkVolume>::New();
volume->SetMapper(mapper);
volume->SetProperty(prop);
volume->SetScale(1.0, exagg, exagg); // 垂向夸张:薄轴放大,截面结构才看得出
ren->AddVolume(volume);
auto capWin = vtkSmartPointer<CapturingOutputWindow>::New();
vtkOutputWindow::SetInstance(capWin);
// 「看进体内部」取景:斜穿俯视——较大 Elevation 让视线从上方斜穿过体内部(而非只看
// 端面)配合垂向夸张后体呈板状俯视看穿层间。Zoom 填满画面。
ren->ResetCamera();
vtkCamera* cam = ren->GetActiveCamera();
cam->Elevation(elevation);
cam->Azimuth(azimuth);
cam->Zoom(zoom);
ren->ResetCameraClippingRange();
rw->Render();
// 结构像素:任一通道 >50排除暗背景
auto countStructPixels = [&]() -> vtkIdType {
auto px = vtkSmartPointer<vtkUnsignedCharArray>::New();
rw->GetRGBACharPixelData(0, 0, winW - 1, winH - 1, /*front=*/1, px);
vtkIdType n = 0;
const vtkIdType np = px->GetNumberOfTuples();
for (vtkIdType i = 0; i < np; ++i) {
if (px->GetComponent(i, 0) > 50 || px->GetComponent(i, 1) > 50 ||
px->GetComponent(i, 2) > 50) {
++n;
}
}
return n;
};
// 高于背景像素:背景为 (0.05,0.05,0.09)≈RGB(13,13,23),阈值 35 干净地把渲出的体
// 结构与背景分开(区别于结构像素>50>35 能纳入梯度门控后偏暗但确有的层面)。
auto countAboveBg = [&]() -> vtkIdType {
auto px = vtkSmartPointer<vtkUnsignedCharArray>::New();
rw->GetRGBACharPixelData(0, 0, winW - 1, winH - 1, /*front=*/1, px);
vtkIdType n = 0;
const vtkIdType np = px->GetNumberOfTuples();
for (vtkIdType i = 0; i < np; ++i) {
if (px->GetComponent(i, 0) > 35 || px->GetComponent(i, 1) > 35 ||
px->GetComponent(i, 2) > 35) {
++n;
}
}
return n;
};
const vtkIdType structPx = countStructPixels();
const vtkIdType aboveBg = countAboveBg();
const double bright = meanBrightness(rw.Get(), winW, winH);
savePng(rw.Get(), pngPath);
// fps旋相机 frames 帧(保持大俯角,绕 Azimuth
rw->Render();
Stopwatch sw;
for (int f = 0; f < frames; ++f) {
cam->Azimuth(360.0 / frames);
rw->Render();
}
const double ms = sw.elapsedMs();
const double fps = ms > 0 ? frames * 1000.0 / ms : 0.0;
const bool texErr = capWin->textureError();
vtkOutputWindow::SetInstance(nullptr);
// 有效判据:无纹理错 + 渲出高于背景的像素(>35。结构像素(>50)仅作亮度强弱度量,
// 不作有效门——梯度门控后层面偏暗但确有渲出,不应误判为空。
const bool ok = !texErr && aboveBg > 0;
const char* label =
mode == 0 ? "a-value(基线)"
: (mode == 1 ? "b-grad(+梯度不透明度)" : "c-grad-shade(+梯度+光照)");
std::cout << "\n--- polish " << label << " ---\n";
std::cout << "存图 : " << pngPath << "\n";
std::cout << "高于背景像素(>35): " << aboveBg << " / " << (winW * winH) << " ("
<< (100.0 * aboveBg / (winW * winH)) << "%)\n";
std::cout << "结构像素(>50) : " << structPx << " / " << (winW * winH) << " ("
<< (100.0 * structPx / (winW * winH)) << "%)\n";
std::cout << "平均亮度(0-255) : " << bright << "\n";
std::cout << "真实 fps : " << (ok ? std::to_string(fps) : "INVALID")
<< " (" << frames << " 帧旋相机)\n";
std::cout << "结果 : " << (ok ? "OK" : "FAIL(纹理错或空渲染)") << "\n";
writeMetricLine(
"polish," + std::string(mode == 0 ? "a-value" : mode == 1 ? "b-grad"
: "c-grad-shade") +
",aboveBg=" + std::to_string(aboveBg) +
",structPx=" + std::to_string(structPx) +
",bright=" + std::to_string(bright) +
",fps=" + (ok ? std::to_string(fps) : "INVALID") +
",texErr=" + std::to_string(texErr ? 1 : 0) + ",png=" + pngPath);
return ok ? 0 : 1;
}
int cmdPolish(int argc, char** argv) {
const Args a = parseArgs(argc, argv, 2);
if (a.positional.empty()) {
std::cerr << "用法: gpr_poc polish <storeDir> [--exagg 8] [--frames 90] "
"[--localBricks 4]\n";
return 2;
}
const std::string dir = a.positional[0];
const double exagg = std::stod(a.get("exagg", "8"));
const int frames = std::stoi(a.get("frames", "90"));
const int localBricks = std::stoi(a.get("localBricks", "4"));
// 「看进体内部」取景:斜穿俯视。默认 El45/Az30/Zoom1.5(视线从上方斜穿层间,
// 既不是纯端面也不至于过陡退化成边缘线)。可命令行覆盖以微调。
const double elevation = std::stod(a.get("elevation", "45"));
const double azimuth = std::stod(a.get("azimuth", "30"));
const double zoom = std::stod(a.get("zoom", "1.5"));
std::cout << "[polish] storeDir=" << dir << " exagg=" << exagg
<< " frames=" << frames << " localBricks=" << localBricks << "\n";
std::cout << "[polish] 离屏闸门复检...\n";
if (cmdOffscreenSmoke() != 0) {
std::cout << "[polish] 闸门失败,中止。\n";
return 1;
}
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 = makeSeismicColorScale(vmin, vmax);
// 全分辨率 level0 局部段(沿线中段),三图共用同一体。
const int totBx = store.bricksX(0);
const int localBx = std::min(localBricks, totBx);
const int bx0 = std::max(0, totBx / 2 - localBx / 2);
vtkSmartPointer<vtkImageData> locImg =
buildLocalLevel0Image(store, m, bx0, localBx);
int locDims[3];
locImg->GetDimensions(locDims);
std::cout << "[polish] level0 局部段 " << locDims[0] << "x" << locDims[1] << "x"
<< locDims[2] << " (brick列 [" << bx0 << "," << (bx0 + localBx)
<< ") / " << totBx << ")\n";
// 标定梯度不透明度阈值:采样实际梯度幅值分布。
const GradStats gs = sampleGradientMagnitude(locImg.Get());
std::cout << "[polish] 梯度幅值分布(量化域,样本 " << gs.samples
<< "): median=" << gs.median << " p75=" << gs.p75
<< " p90=" << gs.p90 << " p99=" << gs.p99 << " max=" << gs.mx
<< "\n";
std::cout << "[polish] 梯度不透明度 piecewise: grad<=" << std::max(1.0, gs.median)
<< "→0.0 grad=" << std::max(2.0, gs.p90) << "→0.5 grad>="
<< std::max(3.0, gs.p99) << "→0.9\n";
const fs::path shotDir =
fs::path("docs") / "superpowers" / "plans" / "poc-lod-shots";
fs::create_directories(shotDir);
int rc = 0;
rc |= renderPolishOne(locImg.Get(), m.quant, cs, vmin, vmax, gs, /*mode=*/0,
exagg, (shotDir / "polish-a-value.png").string(), frames,
elevation, azimuth, zoom);
rc |= renderPolishOne(locImg.Get(), m.quant, cs, vmin, vmax, gs, /*mode=*/1,
exagg, (shotDir / "polish-b-grad.png").string(), frames,
elevation, azimuth, zoom);
rc |= renderPolishOne(locImg.Get(), m.quant, cs, vmin, vmax, gs, /*mode=*/2,
exagg, (shotDir / "polish-c-grad-shade.png").string(),
frames, elevation, azimuth, zoom);
std::cout << "\n[polish] 完成3 张对比图存于 " << shotDir.string()
<< " (polish-a-value / polish-b-grad / polish-c-grad-shade)\n";
return rc;
}
void usage() {
std::cerr << "gpr_poc —— POC-B headless 度量 CLI\n"
" gpr_poc build <dir> [--line 001] [--cellXY 0.2] "
@ -3036,7 +3370,9 @@ void usage() {
"[--bricks 4,16,64,128,256]\n"
" gpr_poc view <storeDir> [--exagg 8] [--opacity 0.5] "
"[--smoke] [--preview] [--near] [--variant N] [--gallery] "
"[--frames 90]\n";
"[--frames 90]\n"
" gpr_poc polish <storeDir> [--exagg 8] [--frames 90] "
"[--localBricks 4]\n";
}
} // namespace
@ -3064,6 +3400,7 @@ int main(int argc, char** argv) {
if (cmd == "tune") return cmdTune(argc, argv);
if (cmd == "fps-budget") return cmdFpsBudget(argc, argv);
if (cmd == "view") return cmdView(argc, argv);
if (cmd == "polish") return cmdPolish(argc, argv);
} catch (const std::exception& e) {
std::cerr << "错误: " << e.what() << "\n";
return 1;