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

关注

原创 2025-05-21 13:48:14 1.1k 阅读 AI 写同款· GEO 优化 ›

从CT到MRI：用RadiAnt DICOM Viewer的3D MPR和VR功能，5步搞定临床影像分析本文详细介绍了如何利用RadiAnt DICOM Viewer的3D MPR和VR功能，通过五个步骤高效完成临床影像分析。从数据导入优化到三维重建技巧，再到精准测量与多模态数据融合，帮助医生快速掌握这一强大工具，提升CT和MRI影像诊断的准确性和效率。阅读详情

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

1 概述

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

https://blog.youkuaiyun.com/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和edge2
	float3 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 × edge2
	float 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 · pvec
	float det = edge1.x * pvec_x + edge1.y * pvec_y + edge1.z * pvec_z;

	// 处理det为负的情况：交换v1和v2，反转法向量方向
	if (det < 0.0f) {
		// 交换v1和v2
		float3 temp = v1;
		v1 = v2;
		v2 = temp;

		// 重新计算edge1和edge2
		edge1 = 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和det
		pvec_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 × edge1
	float 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.youkuaiyun.com/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;
		else
			bounds[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;
}