医疗影像中，DICOM点云、三角面片实体混合渲染（VR）

此文章，涉及到专业性比较强，所以，大部分的内容，基本上都是示例代码的形式出现。以下的技术路径，完全经过实践验证，并且效果很好，可以放心使用。

1 概述

在医学影像中，对DICOM的渲染，通常采用Phong 冯氏光照模型来实现。具体细节可以参考我写的这篇文章

https://blog.csdn.net/rendaweibuaa/article/details/128910050?spm=1001.2014.3001.5501

这里只对一些VR渲染中，关键参数这里做一些补充。

1.1 渲染物体梯度计算（N计算）

为了使沿着光线的方向显示的更平滑一下，梯度计算可以仿照滤波因子写一个梯度计算的方式；

1.2 阴影效果

为了显示效果更加逼真，近的物体更亮，远的物体更暗，这里采用距离的倒数作为亮度的权重。

上图明显可以看出来，近处亮远处暗的效果；

1.3 材质

因为DICOM中，没有材质属性，但是实际的渲染过程中，可以根据渲染骨，软组织来调整漫反射和镜面反射的参数。

2 STL格式的渲染

在工程中，通常有需要将STL格式和VR一起进行渲染的需求。例如假体的植入，例如牙科，关节置换等。STL格式这里就不具体展开介绍了。以下代码是，三角面片追踪渲染的代码；我把主要的注释都写在代码中。

/**
三角面片追踪渲染函数
输入是一个当前渲染前的rgba 的值
当前三角面的col
光线和三角片的交点 interactPoint
三角面片 trifaceV1 trifaceV2 trifaceV3
光线方向 dirLight
*/
__device__ float4 __trifaceTracing(float4 sum,  //累计颜色float4 col, // 颜色和透明度float alphaAccObject,float3 interactPoint, // 位置float3 trifaceV1,  // 三角面第一个顶点float3 trifaceV2,  // 三角面第二个顶点float3 trifaceV3,  // 三角面第三个顶点float3 dirLight,   // 光线方向bool invertZ,     float distanceInterAndO, // 交点与源点距离float vrBrightness       // 亮暗的权值)
{float3 v1 = trifaceV2 - trifaceV1;float3 v2 = trifaceV3 - trifaceV2;float3 N = cross(v1, v2);   // 法向量N = normalize(N);float diffuse = dot(N, dirLight);if (diffuse < 0)// 如果是表面渲染，光照在背面，不应该再有反射；但是对于当前的透视来说，直接N反过来；{diffuse = -1 * diffuse;}diffuse = diffuse < 5e-7f ? 1e-6f : diffuse;// Ka + Kd + Ks ≤ 1 这样可以避免光照过曝。不过，具体的比值可能因不同的材质和光照需求而有所不同float dr = distanceInterAndO + RAY_SOURCE_OFFSET; // RAY_SOURCE_OFFSET 用来调整距离，从而调整阴影效果float dr2 = pow(dr ,2);float invDr2 = 1.0f / dr2;float4 clrLight = col * 0.1f;// 以上已经将此点的Phong都计算出来float4 f4Temp = make_float4(0.0f);f4Temp = col * ( ((diffuse * 0.65f  + 0.15f  * (pow(diffuse, 64.0f))) * vrBrightness) * invDr2 );clrLight += f4Temp;float alphaWeightLeft = (1.0f - alphaAccObject) * col.w; // 给后续光线上的点，留下来的多少透明度(1.0f - alphaAccObject) 和 当前本身当前点的权值 乘return (sum + alphaWeightLeft * clrLight);
}

下面的代码是用来判断光线和三角面片是否相交

__device__ __forceinline__ bool rayIntersectTriangle(const float3& orig, const float3& dir,float3 v0, float3 v1, float3 v2, // 三角形的三个顶点float* t, float* u, float* v,float3* intersection)
{// 初始计算edge1和edge2float3 edge1 = make_float3(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z);float3 edge2 = make_float3(v2.x - v0.x, v2.y - v0.y, v2.z - v0.z);// 计算pvec = dir × edge2float pvec_x = dir.y * edge2.z - dir.z * edge2.y;float pvec_y = dir.z * edge2.x - dir.x * edge2.z;float pvec_z = dir.x * edge2.y - dir.y * edge2.x;// 计算行列式det = edge1 · pvecfloat det = edge1.x * pvec_x + edge1.y * pvec_y + edge1.z * pvec_z;// 处理det为负的情况：交换v1和v2，反转法向量方向if (det < 0.0f) {// 交换v1和v2float3 temp = v1;v1 = v2;v2 = temp;// 重新计算edge1和edge2edge1 = make_float3(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z);edge2 = make_float3(v2.x - v0.x, v2.y - v0.y, v2.z - v0.z);// 重新计算pvec和detpvec_x = dir.y * edge2.z - dir.z * edge2.y;pvec_y = dir.z * edge2.x - dir.x * edge2.z;pvec_z = dir.x * edge2.y - dir.y * edge2.x;det = edge1.x * pvec_x + edge1.y * pvec_y + edge1.z * pvec_z;}// 检查行列式是否接近零（考虑绝对值）if (fabsf(det) < 5e-7f){return false;}// 计算inv_det并继续后续步骤float inv_det = 1.0f / det;float3 tvec = make_float3(orig.x - v0.x, orig.y - v0.y, orig.z - v0.z);// 计算u并检查范围float local_u = (tvec.x * pvec_x + tvec.y * pvec_y + tvec.z * pvec_z) * inv_det;if (local_u < 0.0f || local_u > 1.0f) return false;// 计算qvec = tvec × edge1float qvec_x = tvec.y * edge1.z - tvec.z * edge1.y;float qvec_y = tvec.z * edge1.x - tvec.x * edge1.z;float qvec_z = tvec.x * edge1.y - tvec.y * edge1.x;// 计算v并检查范围float local_v = (dir.x * qvec_x + dir.y * qvec_y + dir.z * qvec_z) * inv_det;if (local_v < 0.0f || (local_u + local_v) > 1.0f) return false;// 计算t并检查正值float local_t = (edge2.x * qvec_x + edge2.y * qvec_y + edge2.z * qvec_z) * inv_det;if (local_t < 0.0f) return false;// 写入结果*t = local_t;*u = local_u;*v = local_v;// 计算交点坐标if (intersection) {intersection->x = orig.x + local_t * dir.x;intersection->y = orig.y + local_t * dir.y;intersection->z = orig.z + local_t * dir.z;}return true;
}

3 加速渲染（BVH树）

通过以上代码计算后，确实可以将三角面片和VR数据一次渲染出来。但是，当三角面片一多的时候，就会发现，渲染速度会变的非常慢。

此时，我们可以通过创建Bvh树来加速渲染。

原理可以参考https://blog.csdn.net/VIPCCJ/article/details/119550359，完全将时间复杂度从O(n)降低到O(logn)。

以下是创建BVH树的代码

// 计算包围盒
void BvhMethods::BvhCalculateBounds(const Facet3D* facets, int start, int end, float bounds[6]) {for (int i = 0; i < 6; ++i){if (i % 2 == 0)bounds[i] = FLT_MAX;elsebounds[i] = -FLT_MAX;}for (int i = start; i < end; ++i){const float* box = facets[i].boxes;for (int j = 0; j < 6; ++j){if (j % 2 == 0){if (box[j] < bounds[j])bounds[j] = box[j];}else{if (box[j] > bounds[j])bounds[j] = box[j];}}}
}// 合并两个包围盒
void BvhMethods::BvhMergeBounds(const float a[6], const float b[6], float out[6])
{for (int i = 0; i < 6; ++i) {out[i] = (i % 2 == 0) ? fminf(a[i], b[i]) : fmaxf(a[i], b[i]);}
}// 构建BVH
int BvhMethods::BvhBuildBVH(const std::vector<Facet3D>& facets, std::vector<BvhNode>& nodes, int start, int end) {int nodeIndex = nodes.size();nodes.push_back(BvhNode());// 计算当前节点的包围盒BvhCalculateBounds(&facets[0], start, end, nodes[nodeIndex].bounds);// 如果只有一个三角形，创建叶节点if (end - start == 1) {nodes[nodeIndex].left = -1;nodes[nodeIndex].right = -1;nodes[nodeIndex].start = start;nodes[nodeIndex].count = 1;return nodeIndex;}// 找到最佳分割轴float axisExtents[3];for (int axis = 0; axis < 3; ++axis) {float minVal = facets[start].boxes[axis];float maxVal = minVal;for (int i = start + 1; i < end; ++i) {minVal = fminf(minVal, facets[i].boxes[axis]);maxVal = fmaxf(maxVal, facets[i].boxes[axis]);}axisExtents[axis] = maxVal - minVal;}int splitAxis = 0;if (axisExtents[1] > axisExtents[splitAxis]) splitAxis = 1;if (axisExtents[2] > axisExtents[splitAxis]) splitAxis = 2;// 按分割轴排序三角形std::vector<int> indices(end - start);for (int i = 0; i < end - start; ++i) indices[i] = start + i;std::sort(indices.begin(), indices.end(), [splitAxis, &facets](int a, int b) {return facets[a].boxes[splitAxis] < facets[b].boxes[splitAxis];});// 找到中间点int middle = (start + end) * 0.5f;int split = middle;// 构建子节点int leftIndex = BvhBuildBVH(facets, nodes, start, middle);int rightIndex = BvhBuildBVH(facets, nodes, middle, end);nodes[nodeIndex].left = leftIndex;nodes[nodeIndex].right = rightIndex;nodes[nodeIndex].start = -1;nodes[nodeIndex].count = -1;return nodeIndex;
}