CT仿真X射线成像¶

CT仿真X射线成像（数字重投影 / 虚拟X-ray）主要有两类核心方法：DRR 与体绘制。

绿色体绘制
蓝色 DRR

1 DRR¶

Digitally Reconstructed Radiograph

核心原理：沿射线路径积分CT体素的线性衰减系数，严格模拟真实X射线成像的物理过程。

物理模型：
- Beer-Lambert定律：I = I₀ × exp(-∫μ dl)
指数积分，标准 DRR
- 单能射线积分：将CT值（HU）转换为线性衰减系数μ，沿射线积分获得投影图像
工程简化版本，用于实时性

生成流程：
1. 几何定义：设置虚拟X射线源、探测器位置和方向
2. 射线追踪：从源点到每个像素发射射线
3. 采样积分：沿路径采样体素，累加衰减系数
4. 图像生成：计算最终强度，生成灰度投影图

特点：
- 严格遵循物理衰减模型
- 与真实X光片高度一致
- 计算相对较慢，但GPU加速后可实时

2 体绘制¶

Volume Rendering, VR

核心原理：直接可视化CT值，通过传输函数映射为颜色和不透明度，不考虑物理衰减。

常见方式：
- MIP（Maximum Intensity Projection）
最大密度投影，显示射线路径上 CT 值最高的体素，常用于骨骼、血管快速观察
- Alpha Blending（标准体渲染）
沿射线对颜色和不透明度进行前向/后向合成，形成半透明体积效果
- 等值面渲染（Isosurface）
提取指定 CT 阈值对应的几何表面，偏向几何建模而非投影成像

生成流程：
1. 几何定义：设置观察方向、相机位置与体数据坐标关系
2. 射线投射：从视点或屏幕像素向体数据发射射线
3. 体素采样：沿路径采样体素，查表获得颜色与不透明度
4. 颜色合成：按 Alpha Blending 规则

特点：
- 强调三维空间关系的直观理解
- 可添加光照、阴影等视觉效果
- 梯度增强、突出边界和结构

3 比较¶

方法	物理正确性	像 X-ray 程度	性能	定量	典型用途
体渲染	❌ 低	⚠️ 外观相似	⭐⭐⭐⭐	❌	可视化、导航、教学演示
DRR	✅ 高	✅ 非常像	⭐⭐⭐	✅	配准、导航、系统标定

方法	Demo	手术导航 UI	观察结构	配准	导航	标定
体渲染	✅ 适合	✅ 3D引导	✅ 优秀	❌ 不精确	⚠️ 辅助	❌ 不适用
DRR	⚠️ 精确	⚠️ 2D投影	⚠️ 需经验	✅ 高精度	✅ 核心	✅ 准确

要“像真 X-ray、能算数” → DRR
要“快、直观、好看” → 体绘制

4 VTK实现¶

体渲染

现成的 vtkGPUVolumeRayCastMapper ，x-ray 仿真参数如下

    vtkNew<vtkColorTransferFunction> colorFunc;
    colorFunc->AddRGBPoint(0.0, .0, .0, .0);
    colorFunc->AddRGBPoint(4000.0, .0, .0, .0);

    vtkNew<vtkPiecewiseFunction> opacityFunc;
    opacityFunc->AddPoint(-1000.0, 0.0);
    opacityFunc->AddPoint(50.0, 0.00015);
    opacityFunc->AddPoint(300.0, 0.005);
    opacityFunc->AddPoint(1000.0, 0.1);
    opacityFunc->AddPoint(3000.0, 0.2);

    vtkNew<vtkVolumeProperty> vrProperty;
    vrProperty->SetColor(colorFunc);
    vrProperty->SetScalarOpacity(opacityFunc);
    vrProperty->SetSpecularPower(10.0);
    vrProperty->SetInterpolationTypeToNearest();
    vrProperty->DisableGradientOpacityOn();
    vrProperty->ShadeOff();
    vrProperty->SetAmbient(1.0);
    vrProperty->SetDiffuse(0.0);
    vrProperty->SetSpecular(0.0);

    SetProperty(vrProperty);

DRR

需自定义着色器，完整着色器在 gist。

自定义 vtkMappder -> RenderPiece -> FBO 离屏渲染到纹理 -> QOpenGLWidget 绘制纹理

QVTKOpenGLNativeWidget 使用 Qt 创建的 OpenGL Context，并加入同一 OpenGL share group，使不同窗口可直接共享 GPU 资源

void main()
{
    FragColor = vec4(vec3(.95),1);
    /* ===== 像素 → 探测器平面 ===== */
    vec2 uv = vUV - 0.5;
    vec3 detPos = uDetCenter
                + uv.x * uDetSizeMM.x * uDetU
                + uv.y * uDetSizeMM.y * uDetV;

    /* ===== 射线 ===== */
    vec3 rayOrigin = uSourcePos;
    vec3 rayDir    = normalize(detPos - rayOrigin);

    /* ===== 转到 CT 空间 ===== */
    vec3 roCT = (uWorld2CT * vec4(rayOrigin, 1.0)).xyz;
    vec3 rdCT = normalize((uWorld2CT * vec4(rayDir, 0.0)).xyz);

    float tNear, tFar;
    if (!intersectBox(roCT, rdCT, uCTMin, uCTMax, tNear, tFar))
    {
        return;
    }

    /* ===== 限制最大距离到探测器 ===== */
    vec3 detPosCT = (uWorld2CT * vec4(detPos, 1.0)).xyz;
    float maxDistanceCT = length(detPosCT - roCT);
    tNear = max(tNear, 0.0);
    tFar = min(tFar, maxDistanceCT);
    if (tNear > tFar) {
        return;
    }

    /* ===== DRR 积分 ===== */
    float acc = 0.0;
    float t = max(tNear, 0.0);

    float stepSize = uCTSpacingZ;
    if (uFastMode == 1) {
        stepSize *= 40.0;
    }

    vec3 roTex = (roCT - uCTMin) * uInvCTSize;
    vec3 rdTex = rdCT * uInvCTSize;
    vec3 deltaTex = rdTex * stepSize;
    vec3 posTex = roTex + rdTex * t;

    if (uFastMode == 1) {
        // 减小 jitter 幅度，降低噪声
        float seed = dot(gl_FragCoord.xy, vec2(12.9898, 78.233));
        float rnd = fract(sin(seed) * 43758.5453);
        float jitter = (rnd - 0.5) * stepSize * 0.5;
        posTex += rdTex * jitter;
    }

    int maxSteps = int(ceil((tFar - t) / stepSize));
    int steps = min(MAX_STEPS, maxSteps);

    for (int i = 0; i < steps; ++i)
    {
        float hu = texture(uVolume, posTex).r;
        float mu;
        if (hu < -900.0) {
            mu = 0.0;
        } else if (hu < 100.0) {
            mu = (hu + 100.0) * 1e-6;
        } else {
            mu = (100.0 + (hu - 100.0) * 3.0) * 4e-5;
        }

        acc += mu * stepSize;
        posTex += deltaTex;
    }

    float logI = acc + 0.15;
    logI = min(logI, 3.0);
    float WL = 1.5;
    float WW = 3.0;

    float gray = (logI - (WL - 0.5 * WW)) / WW;
    gray = clamp(gray, 0.0, 1.0);
    gray = 1.0 - gray;
    gray = pow(gray, 0.9);

    FragColor = vec4(vec3(gray), 1.0);
}