基于CUDA的体渲染技术---着色

引言

这篇文章基本可以说一篇翻译文，将Martino Pilia编写的GPU-accelerated single-pass volumetric raycasting in Qt and OpenGL一文翻译成中文，将GLSL代码翻译成CUDA代码。在【基于CUDA的体渲染技术---Ray-Casting】中介绍了使用ray-casting进行体绘制的基本步骤，并用cuda做了简单实现。本文着重讲解一些在医疗领域的体绘制技术：

• 等值面绘制
• 最大强度投影
• alpha融合

等值面绘制

等值面绘制，即对从体积数据中提取的等值面（给定阈值水平集的表面），进行可视化。为了执行等值面提取，我们让光线穿过体数据直到某个采样点的强度大于或等于设定的阈值，即光线击中了等值面。
一般来讲，从体数据中提取的等值面往往是相当嘈杂和不规则的。因此，为了获取一个更为平滑的结果，我们可以通过二分步长来细化它，在距离命中点半步长的范围内搜索强度更接近阈值的采样点。
一旦我们找到了一个属于等值面的点，就可以应用一个光照模型来执行着色。一个最简单的光照模型是冯氏光照模型，冯氏光照模型的主要结构由3个分量组成：环境(Ambient)、漫反射(Diffuse)和镜面(Specular)光照。冯氏光照模型可参考LearnOpenGL CN的基础光照篇。为了应用冯氏光照模型，我们需要直到在该点处表面的法向量。这可以简单的使用改采样点的强度的梯度来近似，代价是三次纹理采样，假设数据是合理连续的。

__device__
float3 normal(float3 position, float intensity, float step, cudaTextureObject_t volumeTex)
{
    float dx = tex3D<float>(volumeTex, (position.x + step)*0.5f + 0.5f, position.y*0.5f + 0.5f, position.z*0.5f + 0.5f) - intensity;
    float dy = tex3D<float>(volumeTex, position.x*0.5f + 0.5f, (position.y + step)*0.5f + 0.5f, position.z*0.5f + 0.5f) - intensity;
    float dz = tex3D<float>(volumeTex, position.x*0.5f + 0.5f, position.y*0.5f + 0.5f, (position.z + step)*0.5f + 0.5f) - intensity;
    return -normalize(mul(c_normalMatrix, make_float3(dx, dy, dz)));
}

当法向量求取了后，我们就有了用于计算光照的所有的信息，光线行进循环大概如下(原文使用GLSL编写的，这里使用Cuda改写)：

        for (int i = 0; i < maxSteps && t < tfar; i++, t += tstep, pos += step)
        {
            float intensity = tex3D<float>(volumeTex, pos.x*0.5f + 0.5f, pos.y*0.5f + 0.5f, pos.z*0.5f + 0.5f);
            if (intensity > threshold)
            {
                // Get closer to the surface
                pos -= step * 0.5;
                intensity = tex3D<float>(volumeTex, pos.x*0.5f + 0.5f, pos.y*0.5f + 0.5f, pos.z*0.5f + 0.5f);
                pos -= step * (intensity > threshold ? 0.25 : -0.25);
                intensity = tex3D<float>(volumeTex, pos.x*0.5f + 0.5f, pos.y*0.5f + 0.5f, pos.z*0.5f + 0.5f);

                // Blinn-Phong shading
                float3 L = normalize(make_float3(3.0, 0.0, 3.0) - pos);
                float3 V = -normalize(eyeRay.d);
                float3 N = normal(pos, intensity, tstep, volumeTex);
                float3 H = normalize(L + V);
                
                float Ia = 0.1;
                float Id = 1.0 * max(0.0f, dot(N, L));
                float Is = 8.0 * pow(max(0.0f, dot(N, H)), 600);
                colour = make_float4((Ia + Id) * make_float3(1.0) + Is * make_float3(1.0), 1.0f);
                
                break;
            }
        }

等值面渲染技术非常适合可视化CT数据的骨骼，因为骨骼的放射密度明显高于其他组织,
在Hounsfield单位方面给他们一个总体良好的对比。

原文使用等值面渲染的足部骨骼图： foot.png
使用cuda改写后渲染的劲动脉图： carotid.png

最大强度投影

最大强度投影是一种非常简单但却非常有效的体绘制技术。在其最简单的形式中，它是将沿光线取样的最大强度作为当前片段的强度。

float maxIntensity = 0.0f;
for (int i = 0; i < maxSteps && t < tfar; i++, t += tstep, pos += step)
{
    float sample = tex3D<float>(volumeTex, pos.x*0.5f + 0.5f, pos.y*0.5f + 0.5f, pos.z*0.5f + 0.5f);
    if (sample > maxIntensity)
    {
        maxIntensity = sample;
    }
}

这种技术可以很好地将内部区域比外部区域具有更高强度的体积可视化，而且CT数据也是一个很好的候选者，而大多数类型的MR图像在采用MIP时看起来不是很有趣。

原文使用最大强度投影渲染的足部骨骼图： foot.png
使用cuda改写后渲染的劲动脉图： carotid.png

Alpha融合

在【基于CUDA的体渲染技术---Ray-Casting】一文中使用的渲染技术即为alpha渲染。这可能是最通用的技术，它可以产生最复杂的结果，但同时它需要做更多的工作，并且设置通常与要可视化的图像类型有关。基本原理是使用传递函数将强度值映射到颜色，并对沿射线获取的样本执行从前到后的alpha混合。从前到后渲染的重要原因之一是，当片段的颜色饱和时，它可以尽早停止采样。

        for (int i = 0; i < maxSteps && t < tfar; i++, t+=tstep, pos+=step)
        {
            float sample = tex3D<float>(volumeTex, pos.x*0.5f + 0.5f, pos.y*0.5f + 0.5f, pos.z*0.5f + 0.5f);
            
            // lookup in transfer function texture
            float4 col = tex1D<float4>(transferTex, (sample - transferOffset)*transferScale);
            col.w *= density;

            // pre-multiply alpha
            col.x *= col.w;
            col.y *= col.w;
            col.z *= col.w;
            // "over" operator for front-to-back blending
            sum = sum + col * (1.0f - sum.w);

            if (sum.w > opacityThreshold)
                break;
        }

我们还可以使用梯度信息来控制样本的不透明度，这在我们要增强表面时特别有用。

__device__
float gradientMagnitude(float3 position, float intensity, float step, cudaTextureObject_t volumeTex)
{
    float dx = tex3D<float>(volumeTex, (position.x + step)*0.5f + 0.5f, position.y*0.5f + 0.5f, position.z*0.5f + 0.5f) - intensity;
    float dy = tex3D<float>(volumeTex, position.x*0.5f + 0.5f, (position.y + step)*0.5f + 0.5f, position.z*0.5f + 0.5f) - intensity;
    float dz = tex3D<float>(volumeTex, position.x*0.5f + 0.5f, position.y*0.5f + 0.5f, (position.z + step)*0.5f + 0.5f) - intensity;
    return sqrt(dx * dx + dy * dy + dz * dz);
}

虽然alpha混合代码到目前为止看起来很简单，但困难的部分是实际实现颜色传递函数。这没有通用的方法，因为它非常依赖于要可视化的数据类型以及应强调的数据特征。在【基于CUDA的体渲染技术---Ray-Casting】中是使用一个颜色映射纹理作为传递函数进行alpha融合。加载颈动脉数据渲染效果： carotid.png

参考文章：
https://martinopilia.com/posts/2018/09/17/volume-raycasting.html
https://learnopengl-cn.github.io/02%20Lighting/02%20Basic%20Lighting/