GPU Gems 1 · Chapter 40. Applying Real-Time Shading to 3D Ultrasound Visualization

GPU Gems 1 第 40 章把上一章的通用体渲染推进到一个更棘手的场景：实时显示时变的三维医学超声。CT 和 MRI 通常采样在笛卡尔网格；超声探头发出声束，越离开探头，采样线越发散，天然形成金字塔、圆锥或柱面网格。心脏等器官还要求每秒渲染多个体积，才能让运动可被观察。

本章不是医疗诊断建议，也不把历史 OpenGL/Cg API 当成现代实现的唯一答案。重点是数据空间、坐标变换、GPU 阶段职责与填充率预算：今天可以用现代 3D texture、vertex/fragment shader、compute pipeline 或其他图形 API 实现同一数据流，但几何关系仍然要被正确表达。

1. 实时超声的任务：在时间维度里看懂体数据

3D ultrasound 不只是“很多二维图片叠起来”。扫描器沿 range 方向记录离探头的距离，沿 azimuth 方向记录横向声束位置，沿 elevation 方向记录切面位置；如果再加入时间，每一帧就是一个 volume，整个序列就是 4D 数据。

医学用户熟悉二维影像，却不容易直接理解三维数组内部的结构。体渲染通过颜色、透明度和阴影把内部组织映射成一张画面，但实时场景多了一个硬约束：视角、传递函数和 volume index 可能持续变化，渲染器必须把数据留在 GPU，避免每一帧在 CPU 与 GPU 之间搬运完整体积。

一条可执行的管线可以这样拆：

1. 输入：扫描器产生声学网格中的样本，主机负责上传或选择当前 volume。
2. 顶点阶段：根据相机、数据空间和金字塔参数生成切平面顶点与带 w 的纹理坐标。
3. 片段阶段：修正 projective coordinate，采样 3D ultrasound volume，再查颜色和透明度表。
4. 合成：按切片顺序用 alpha blending 累积，必要时与不透明几何共享深度。

对于快速运动的心脏，体积速率和每个体积的切片数会共同决定帧预算。先做数据空间和采样预算，往往比一开始调颜色更能避免“静态示例很好、实时序列卡顿”的返工。

2. 从笛卡尔网格到声学网格

CT 和 MRI 常见的笛卡尔网格中，x、y、z 方向的采样线平行，3D texture 坐标可以随切平面线性插值。超声声束从探头上的 aperture 发出，在更远的 range 处覆盖更大的横向区域。把它画成二维截面，就是一个顶点位于 apex 的金字塔；三维情况下则形成金字塔或近似圆锥。

pyramidal grid 可以被看成“加了深度缩放的笛卡尔网格”。令 r 表示归一化 range，原始的 s、t 表示 azimuth 和 elevation；随着 r 增大，横向范围需要按 aperture 与 apex angle 做线性调整。重要的是只缩放 s、t，不要把 range 也错误地缩放，否则体数据会沿深度方向越看越长。

这个抽象有两个好处。首先，网格的非线性几何被压缩成少数 uniform 参数；其次，切片仍然可以是包围体的矩形，顶点程序只需传递随位置变化的坐标。对真实探头，如果声学网格是圆锥或其他形状，可以先沿 range 方向重采样到金字塔网格，或者在同一坐标框架中扩展缩放函数。

3. 切平面：怎样让光栅化器遍历体

cut-plane 不是解剖表面，而是一个采样平面。将多层切平面按远到近或近到远绘制，每个片段都可以取得一个 3D texture 坐标，片段程序再完成分类和合成。

有两种常见代理几何：包围矩形和切平面与体包围盒的交集多边形。包围矩形简单，顶点可以预存，适合视角变化频繁的交互；代价是矩形覆盖体外区域，需要 clipping 或边界纹理值抑制无效片段。交集多边形只光栅化体内区域，能减少片段，但相机方向改变时必须重算交点，CPU 重算后还可能需要重新上传顶点。

一个简单的渲染循环可以表达成：

for plane in sortedCutPlanes:
    bind(volumeTexture, colorMap)
    draw(plane)
    fragment:
        texCoord = projectiveCoord(vertexOutput)
        value = sample3D(volumeTexture, texCoord)
        colorAlpha = lookup(colorMap, value)
        blendOver(colorAlpha)

如果使用包围矩形，必须定义体外纹理坐标的行为：clamp to border、clamp to edge 并把边界值设为零，或使用 clipping planes。不要依赖“纹理坐标刚好不会越界”的假设；相机旋转、插值精度和金字塔缩放都会把片段带到原始体数据之外。

4. Projective texture mapping：把发散关系交给 GPU

projective texture mapping 的核心是把 range 相关的缩放编码进坐标的 w。顶点程序先计算笛卡尔坐标对应的 s、t、r，再让 w 表示当前位置需要的金字塔缩放；片段程序接到插值后的坐标，完成最终的 r 修正和 projective lookup。

设 normalized aperture 为 a，金字塔参数为 k，可以把缩放写成概念公式：

w = a + r · k
s′ = s - 0.5 + 0.5 · w
t′ = t - 0.5 + 0.5 · w
r′ = r · w

这里的具体常数取决于数据空间、包围盒和顶点坐标约定；不变的部分是横向坐标围绕中心轴缩放，range 坐标不能被同一缩放误伤。因为 r′ 是非线性结果，不能只在顶点阶段输出一个已经除好的值，让三角形内部做线性插值；必须把足够的信息传入片段阶段。

片段 shader 最小逻辑可以是：

float3 projectiveSample(float4 coord, Texture3D volume, Texture1D colorMap) {
    coord.z = coord.z * coord.w;
    float value = sample3DProjective(volume, coord);
    return lookupColor(colorMap, value);
}

现代 API 的函数名会不同，但调试时可以把坐标分成三份可视化：原始数据空间坐标、经过 w 修正的纹理坐标、最终命中的体素值。只看最后的颜色，很难分辨是坐标错误、颜色表错误还是 alpha 合成错误。

5. 分类与合成：颜色表不只是 colormap

opacity-weighted transfer function 把体数据值映射成 RGBA 光学属性。RGB 应该已经乘过 alpha，这样硬件用 source = 1、destination = 1 - sourceAlpha 的 over blending 时，当前样本对最终颜色的贡献与吸收一致。

假设第 i 个切片送入的片段颜色是 I_i、透明度是 a_i，上一层的帧缓冲像素是 P_{i-1}，概念上的更新为：

P_i = (1 - a_i) · P_{i-1} + I_i

从远到近绘制时，前方不透明样本自然会压低后方贡献；从近到远绘制时要使用相应的 under 规则。切片顺序、blend factor 和深度测试必须成套验证，否则同一 volume 旋转时可能出现亮度、遮挡和结构关系反转。

传递函数如果直接输出颜色和 alpha，再在片段程序中相乘，可能减少量化误差，却会增加指令和 fill-rate 成本。预乘表能把这一步放到表格构建阶段，适合对每个片段都执行同一查找的实时路径。对医学数据，颜色和透明度的语义也要写入数据契约，避免把显示用的调色板误当作测量结果。

6. 时变 volume：为心脏而不是单张 CT 优化

3D 超声的难点不是只渲染一个大体积，而是连续接收多个体积。对于较慢的胎儿面部，较低 volume rate 也许足够；对于快速运动的心脏，需要在每个心动周期里渲染更多 volume，才能减少运动跳跃。可操作的时间预算是：

frame budget = upload / select volume
             + cut-plane update
             + rasterization and fragment shading
             + blending and display

尽量让 3D texture、颜色表和切片顶点驻留在 GPU。时间变化只切换 volume layer、更新少量纹理或使用环形缓冲；如果每个 volume 都先读回 CPU、重新编码、再完整上传，所谓 GPU 加速可能被同步和传输抹平。

如果视角变化，包围矩形可以复用但可能产生更多体外片段；交集多边形可以减少有效片段，但增加几何重算。对时变数据，频繁视角和频繁 volume update 的组合要分别测量，不能用单张静态 CT 的结果推断医疗交互体验。

7. 性能：识别 fill-rate bound 的真实来源

fill-rate bound 是这类体渲染的常见状态。体纹理可能只有几百个 slice，但每个 slice 可以覆盖很大的屏幕区域；透明合成还会让每个片段反复读取和更新颜色缓冲。

可以按以下顺序排查：

1. 片段覆盖：包围矩形是否让大量体外区域进入片段阶段？开启 clipping 后有效片段减少多少？
2. 切片数：减少采样率是否降低画质，或者只是让 opacity correction 失配？
3. 片段指令：投影坐标、3D lookup、颜色表查找和额外的逆运算分别占多少成本？
4. 带宽与格式：体纹理格式是否过宽，颜色表是否能预乘，缓存是否能隐藏依赖读取？
5. 几何更新：相机变化时，交集多边形重算和顶点上传是否超过省下的 fragment 工作？

预积分体渲染可以减少切片间的 moiré 与采样伪影，但通常会牺牲一部分帧率；应该把它作为质量旋钮，不是默认免费升级。手工优化片段程序集在旧硬件上可能有效，现代实现则应先用 GPU profiler 和 shader compiler 输出确认瓶颈，再决定是否改写坐标计算。

8. 复用清单：把非笛卡尔体数据做成可验证管线

1. 数据空间：写清 range、azimuth、elevation 的轴顺序、值域、aperture、apex angle 和 volume rate。
2. 坐标空间：保留数据空间、切平面空间、视空间和纹理空间的转换，避免在一个矩阵里混入未说明的缩放。
3. 顶点/片段边界：顶点阶段生成带 w 的坐标，片段阶段执行非线性 range 修正和 3D lookup。
4. 切片策略：在包围矩形和交集多边形之间用有效片段与几何更新成本做选择。
5. 合成证据：验证切片顺序、预乘 alpha、blend factor、clipping 和体外纹理值。
6. 时间预算：用目标 volume rate、相机变化和目标设备 profiler 测量 upload、fill rate、shader、blend 与 display。

小结

• 3D 超声通常来自发散的非笛卡尔声学网格；把它当成普通均匀立方体会造成深度拉伸和纹理错位。
• 金字塔网格可以用 aperture、apex 和 range 相关缩放近似描述；projective texture mapping 让 GPU 在片段阶段恢复正确的三维采样坐标。
• 切平面是遍历体数据的代理几何；包围矩形简单但会产生体外片段，交集多边形有效但需要随视角重算。
• opacity-weighted transfer function、alpha blending 和切片顺序共同决定图像是否正确；不能只调颜色表而忽略合成公式。
• 实时 4D 超声常受 fill rate、片段指令、3D texture 带宽和 volume rate 约束；clipping、预积分、采样率和 GPU resident 数据是可测量的工程旋钮。