基于CUDA的体渲染技术---Ray-Casting

Ray-Casting基础

引言

“渲染”指的是给定一个场景的3D表示，产生像素颜色值的整个过程。
传统建模中，3D对象都是使用表面表示，比如网格数据。光线传输通过简单的冯氏模型或者复杂的BRDF算法，使用一些表面属性，比如颜色、粗糙度、反射度等，只根据在表面上的点而被估计。这些方法通常缺乏对打在大气层或者物体内部光线的交互的能力。
和表面渲染相对应的是体渲染技术，它最初兴起于被科学可视化。CT、MRI、流体力学等。随着高效体绘制技术的发展，体数据在视觉艺术和计算机游戏中也变得越来越重要。针对模糊物体，如云、雾、火等。实现一些高质量的特效。通常结合表面渲染一起使用。

光学模型

• 吸收模型
• 发射模型
• 吸收+发射模型，最为常见
• 散射和阴影模型
• 多散射模型

Emission-Absorption.png

ray-casting算法基本步骤

在发射吸收光学模型的基础上，通过对光的相互作用效应沿观测光线的积分来计算光的传播。对于每个像素：

• 以视点作为起点，引出一条光线
• 对于场景中的每个物体：
• 检测光线与物体是否相交，如果相交，则查找出物体与射线的交点位置，进入下一个步骤。否则结束这条光线累积
• 从较近交点开始，沿着光线每次前进规定步长，获取采样位置对应的值，累加求和。直到走到较远处交点，或者满足其他退出条件。
• 设置该像素点的颜色值

体数据

一般从设备中获取的数据都是离散的三维数组，但用从连续的三维信号中获得的一个无穷小点的样本来识别每个体素更合适。通常使用最邻近滤波或者线性滤波替代理想的sinc-filter滤波。在cuda中可以使用3D纹理对象或者3D纹理引用存储体数据：

void initCuda(uchar* h_data, cudaExtent volumeSize)
{
    // 申请3D Array
    cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc <uchar1>();
    checkCudaErrors(cudaMalloc3DArray(&cuArray, &channelDesc, volumeSize));
    
    // 数据拷贝
    cudaMemcpy3DParms parms;
    memset(&parms, 0, sizeof(parms));
    parms.dstArray = cuArray;
    parms.srcPtr = make_cudaPitchedPtr((void*)h_data, volumeSize.width * sizeof(uchar), volumeSize.width, volumeSize.height);
    parms.kind = cudaMemcpyHostToDevice;
    parms.extent = volumeSize;
    checkCudaErrors(cudaMemcpy3D(&parms));

    // 资源描述
    cudaResourceDesc resDesc;
    memset(&resDesc, 0, sizeof(resDesc));
    resDesc.resType = cudaResourceTypeArray;
    resDesc.res.array.array = cuArray;
    
    // 纹理描述
    cudaTextureDesc texDesc;
    memset(&texDesc, 0, sizeof(texDesc));
    texDesc.addressMode[0] = cudaAddressModeClamp;
    texDesc.addressMode[1] = cudaAddressModeClamp;
    texDesc.addressMode[2] = cudaAddressModeClamp;
    texDesc.filterMode = cudaFilterModeLinear;
    texDesc.readMode = cudaReadModeNormalizedFloat;
    texDesc.normalizedCoords = true;

    // 创建纹理对象
    checkCudaErrors(cudaCreateTextureObject(&textureObj, &resDesc, &texDesc, nullptr));
}

光线与包围盒碰撞检测

AABB，即Axis-Aligned Bounding Box，轴对齐包围盒，顾名思义，就是包围盒的每条边（2D空间中）或者每个面（3D空间中）都是轴对齐的（2D空间中）或者面对齐的（3D空间中）。AABB不具备方向性，当物体发生旋转时需要重新计算出其对应的包围盒。其优点是计算简单，但带来的误差也比较明显。其他碰撞检测的方法还有：Sphere、OBB、K-DOP等。光线与AABB相交SlabMethod算法求解过程在https://zhuanlan.zhihu.com/p/128533024文章中非常详细的讲解了，感谢YivanLee的分享。cuda体渲染示例的AABB算法实现：

struct Ray
{
    float3 o;   // origin
    float3 d;   // direction
};

/**
* 检测光线r是否与AABB包围盒发生碰撞，如果发生了碰撞，将碰撞点与光线原点的距离返回
* r       - 光线
* boxmin  - AABB包围盒最低点
* boxmax  - AABB包围盒最高点
* tnear   -  较近碰撞点与光线原点的距离
* tfar    -  较远碰撞点与光线原点的距离
* 返回值  - 如果光线与包围盒发生碰撞，则返回true，否则返回false
*/
__device__
bool intersectBox(Ray r, float3 boxmin, float3 boxmax, float *tnear, float *tfar)
{
    // compute intersection of ray with all six bbox planes
    float3 invR = make_float3(1.0f) / r.d;
    float3 tbot = invR * (boxmin - r.o);
    float3 ttop = invR * (boxmax - r.o);

    // re-order intersections to find smallest and largest on each axis
    float3 tmin = fminf(ttop, tbot);
    float3 tmax = fmaxf(ttop, tbot);

    // find the largest tmin and the smallest tmax
    float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
    float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));

    *tnear = largest_tmin;
    *tfar = smallest_tmax;

    return smallest_tmax > largest_tmin;
}

光线行进

光线从光源处，穿过物体，进入视点。吸收+发射模型表示光线从光源处发射出来，穿过物体时，被物体吸收。穿出物体后的光线进入视野，形成最终的颜色。

RayMarch.png

front-to-back和back-to-front

back-to-front，即光线从光源处或说物体的背面向视点处或说物体的前面前进，光线穿过物体时，不停的被吸收。与back-to-front相反的是front-to-back。相比于back-to-front，front-to-back的好处是，当累积的alpha值达到1.0或足够接近的值时，射线的进程就可提前终止。
成像面中的每个像素都是由一条光线穿过物体形成的一个颜色值。对于每条光线，使用AABB算法计算出光线与物体的交集。每条光线只对穿过物体的光线做计算，就可以避免冗余的无用计算。使用Cuda描述RayMarch算法的基本框架如下：

void d_render(uint *d_output, uint imageW, uint imageH, 
    float density, float brightness, float transferOffset, float transferScale,
    cudaTextureObject_t volumeTex, cudaTextureObject_t transferTex)
{
  // 对成像面中的每个像素，引发一条光线
    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;

    if ((x >= imageW) || (y >= imageH)) return;

  // 光线最大行进步数
    const int maxSteps = 500;
  // 光线每次行进的步长
    const float tstep = 0.01f;
  // 透明度预置，大于该预置时，终止光线行进
    const float opacityThreshold = 0.95f;
  // 包围盒边界
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);
  // 将坐标归一化到[-1, 1)
    float u = (x / (float)imageW)*2.0f - 1.0f;
    float v = (y / (float)imageH)*2.0f - 1.0f;

  // 设置视点位置和光线方向
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    float tnear, tfar;
    if (intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar))  // 光线穿过物体
    {
        if (tnear < 0.0f) tnear = 0.0f;

        // march along ray from front to back, accumulating color
        float4 sum = make_float4(0.0f);
        float t = tnear;
        float3 pos = eyeRay.o + eyeRay.d*tnear;
        float3 step = eyeRay.d*tstep;

      // 光线行进循环
        for (int i = 0; i < maxSteps && t < tfar; i++, t+=tstep, pos+=step)
        {
            // 坐标归一化到纹理坐标范围[0, 1)
            float sample = tex3D<float>(volumeTex, pos.x*0.5f + 0.5f, pos.y*0.5f + 0.5f, pos.z*0.5f + 0.5f);
            
            // 通过传递函数，将采样值映射成颜色值，这里使用一维纹理作为传递函数
            float4 col = tex1D<float4>(transferTex, (sample - transferOffset)*transferScale);
            col.w *= density;

            // 使用front-to-back blending，w作为光线吸收因子
            col.x *= col.w;
            col.y *= col.w;
            col.z *= col.w;
            sum = sum + col * (1.0f - sum.w);

            if (sum.w > opacityThreshold)
                break;
        }
        sum *= brightness;
        // 生成最终颜色输出
        d_output[y*imageW + x] = rgbaFloatToInt(sum);
    }
}

参考Cuda的体渲染示例，使用Qt作为显示界面库，运行效果：

Bucky Ball.png

参考资料
https://zhuanlan.zhihu.com/p/128533024
https://martinopilia.com/posts/2018/09/17/volume-raycasting.html