【Hololens 2】使用Unity fo-dicom实现多平面三维重建(mpr)

使用Unity 2021.3.4在Hololens2上进行MPR（多平面三维重建）的一次尝试，使用了fo-dicom库。

csdn_Joyz

1146人浏览 · 2025-04-22 15:30:12

csdn_Joyz · 2025-04-22 15:30:12 发布

项目背景和需求

在Hololens 2手术项目中，医生需要通过查看切面图来更好的观察穿刺针位置，因此需要将MPR功能在Hololens 2上实现。
样例图片

开发环境

项目	版本
操作系统	Windows 10 家庭中文版 21H2
Unity版本	Unity2021.3.4f1c1
Mixed Reality Toolkit Foundation	2.8.3
Mixed Reality Toolkit Standard Assets	2.8.3
Mixed Reality OpenXR Plugin	1.6.0
fo-dicom(Unity Assets)	3.5.0

涉及的相关技术

Dicom标准与Dicom文件

DICOM（Digital Imaging and Communications in Medicine）即医学数字成像和通信，是医学图像和相关信息的国际标准（ISO 12052）。它定义了质量能满足临床需要的可用于数据交换的医学图像格式。Dicom官方文档.
本项目主要处理的是CT图像的dicom文件。Dicom文件的结构如图所示（Dicom文件简介），其中数据集（DataSet）和像素数据是我们涉及到的部分。
dicom文件结构
在数据集中包含着各类图片信息，具体图片信息的解释可以查看以下网页：DICOM Standard Browser.
在像素数据中，存放的是每个像素的数据，存放顺序是从左到右、从上到下。该数据并非灰度值，对于CT图像，该数据可以通过计算变成CT的HU值。

HU（Hounsfiled Unit）值：反映了组织对X射线吸收程度。以水的吸收程度作为参考，即水的HU=0，衰减系数大于水的为正直，小于水的为负值。并以骨皮质和空气的HU值为上限和下限。空气为-1000，致密骨为+1000。在实际图像中可能是负几千到几千不等。

HU值的计算公式： $输出值 = 像素数据 * 斜率 + 截距。$

其中截距是Rescale Intercept Attribute(Tag(0028,1052))，斜率是Rescale Slope Attribute(Tag (0028,1053))。当ImageType(0008,0008)的值1为ORIGINAL且值3不是LOCALIZER、且Multienergy CT Acquistion(0018,9361)是不存在或者NO的时候，输出值为HU值。

调窗操作：HU值的范围过大，不便于输出和观察，因此要通过调窗操作才能变为图片显示所需的灰度值。
以下是调窗操作的伪代码：

		/*x是输入值，y是范围从ymin到ymax的输出值，
		c是窗口中心（0028，1050），w是窗口宽度（0028，1051）*/
		float ChangeWindow( float x, float c, float w, float yMin, float yMax)
		{
			float y = 0;
			if(x <= c- 0.5f - (w - 1) / 2)
				y=yMin;
			elseif(x > c - 0.5f + (w - 1) / 2)
				y=yMax;
			else
				y = ((x - (c - 0.5f)) / (w - 1) + 0.5f) * (yMax- yMin) + yMin;
			return y;
		}

多平面三维重建

如图所示，多平面三维重建即根据影像数据重建出影像横断面、矢状位、冠状位图像，且切面可以进行角度、位置调整。
MPR示例

fo-dicom库

一个C# dicom库，可用于.NET Framework、.NET Core、Universal Windows、Android、iOS、Mono和Unity，提供了方法进行读取、修改dicom文件。
文档链接
 github链接
 Unity Store链接

Unity的Texture3d

主要用来存储显示图像所需数据
Unity Texture3d 文档

实现流程详解

基本思路

受制于Hololens设备的运算性能以及opengl无法在Hololens上运行等问题，采用了如下方法模拟MPR：
1）从dicom文件中读取原始数据和文件信息（窗位、窗宽等），原始数据处理成HU值用Texture3d保存（受限于Color格式需要除以10000并用Alpha标志正负），文件信息用对应Json文件保存。
2）修改Shader以用于显示Texture3d，首先是Texture3d中的坐标对应，之后是将Color内数值还原为HU值，并根据外部给定的窗位窗宽计算实际灰度
3）给平面赋予对应Material，用来让Shader根据位置渲染出画面，之后用摄像机拍摄画面，使用RenderTexture在对应位置显示。

读取Dicom文件获取所需信息

我们使用fo-dicom进行了dicom文件的读取。首先使用DicomImage从文件夹所有图片中的第一张获取默认窗宽窗位，之后使用dcmFile.Dataset.Get<>()的方式获取其他图片信息、使用自建结构体PixelDataFactory存储像素信息，用三维数组存储HU值数据。代码主要参考了使用fo-dicom读取Dicom文件的PixelData信息及像素信息（C# / fo-dicom）

//path：文件夹位置，personalDicomData：自建的存储信息结构体
static float[][][] ReadDicomPixel(string path, PersonalDicomData personalDicomData)
{
    string folderPath = path; 
    //如果文件夹存在
    if (Directory.Exists(folderPath))
    {
        //打开第一个dcm文件，以获取所有图片应该共同不变的信息：窗位、窗宽
        string tempFilePath = path + "/IM0"; 
        var tempDcmFile = DicomFile.Open(tempFilePath); 
		//通过第一张图片地址创建DicomImage对象
        DicomImage dicomImage = new DicomImage(tempFilePath);
        //将默认的窗位、窗宽存入自己写的结构体中，也可以用后续的dcmFile.Dataset.Get<>()的方式获取
        personalDicomData.defaultWindowLocation = dicomImage.WindowCenter;
        personalDicomData.defaultWindowWidth = dicomImage.WindowWidth;

        //开始文件夹中所有图片的读取
        DirectoryInfo direction = new DirectoryInfo(folderPath);
        FileInfo[] files = direction.GetFiles("*", SearchOption.AllDirectories);
        //我这里用一个三维数组存储读取到的HU值信息，因为有Meta文件所以除以2
        float[][][] pixel_data_3d = new float[files.Length / 2][][];
        personalDicomData.textureWidth = 0;
        for (int i = 0; i < files.Length / 2; ++i)
        {
        	//我这里文件的地址格式为：文件夹地址/IM0、文件夹地址/IM1...
            string fName = path + "/IM" + string.Format("{0}", i);
            var dcmFile = DicomFile.Open(fName); // 打开dcm文件
            //获取第0帧的pixeldata，想获取其他帧修改即可，返回IPixelData类型
            var pixelData = PixelDataFactory .Create(DicomPixelData.Create(dcmFile.Dataset), 0); 

            //获取单个Diocm数据中，（w，h）处的像素信息
            if (pixelData is Dicom.Imaging.Render.GrayscalePixelDataU16)
            {
            	//存储本dcm文件所有点的像素
                float[][] pixel_data_2d = new float[pixelData.Height][]; 
                //读取图像的截距
                float intercept = dcmFile.Dataset.Get<float>(DicomTag.RescaleIntercept);
                //读取图像的斜率
                float slope = dcmFile.Dataset.Get<float>(DicomTag.RescaleSlope);
                //读取像素间距离，以便后续计算Texture3D和实际模型之间的比例
                personalDicomData.pixelSpacing = dcmFile.Dataset.Get<float[]>(DicomTag.PixelSpacing);
                //这里实际上应该用前后两张图片的位置信息相减获得图片见的间距
                personalDicomData.imageSpacing = 1.0f;
                personalDicomData.textureHeight = pixelData.Height;
                personalDicomData.textureWidth = pixelData.Width;
                personalDicomData.textureLength++;
                for (int h = 0; h < pixelData.Height; h++)
                {
                    float[] pixel_data = new float[pixelData.Height];
                    for (int w = 0; w < pixelData.Width; w++)
                    {
                    	//依据公式将原始数据转为HU值，并且为了存储在Unity的Color类型里除以了10000
                        float hu = Convert.ToSingle((pixelData.GetPixel(w, h) * slope + intercept) / 10000);
                        pixel_data[w] = hu;
                    }
                    pixel_data_2d[h] = pixel_data;
                }
                pixel_data_3d[i] = pixel_data_2d;
            }
        }
        return pixel_data_3d;
    }
    else
        return null;
}

存储所需信息及生成Texture3d

将获取到的PerdonalDicomData使用Unity的JsonUtility.ToJson()来存储为Json文件，将读取到的HU值存入Color类型。创建Texture3d并将颜色值复制到纹理存储。代码主要参考了Unity的演示文档。

[MenuItem("Assets/CreateTextureAndJson/CreateFromFolder")]
static void CreateDicomPixel()
{
    var select = Selection.activeObject;
    var path = AssetDatabase.GetAssetPath(select);
    PersonalDicomData personalDicomData = new PersonalDicomData();
    float[][][] dicomImages = ReadDicomPixel(path, personalDicomData);
    
    // 配置纹理
    int size = dicomImages.Length;
    TextureFormat format = TextureFormat.RGBA32;
    TextureWrapMode wrapMode = TextureWrapMode.Clamp;

    // 创建纹理并应用配置
    Texture3D texture = new Texture3D(dicomImages[0][0].Length, dicomImages[0].Length, size, format, false);
    texture.wrapMode = wrapMode;

    // 创建 3 维数组以存储颜色数据
    Color[] colors = new Color[size * dicomImages[0][0].Length * dicomImages[0].Length];

    // 填充数组，使纹理的 x、y 和 z 值映射为红色、蓝色和绿色
    for (int z = 0; z < size; z++)
    {
        int zOffset = z * dicomImages[0][0].Length * dicomImages[0].Length;
        for (int y = 0; y < dicomImages[0].Length; y++)
        {
            int yOffset = y * dicomImages[0].Length;
            for (int x = 0; x < dicomImages[0][0].Length; x++)
            {
            	//由于Color类型无法存储负值，所以使用alpha来作为标志位，0为负，1为正
                if (dicomImages[z][y][x] < 0)
                {
                    colors[x + yOffset + zOffset] = 
                    	new Color(-dicomImages[z][y][x], -dicomImages[z][y][x], -dicomImages[z][y][x], 0);
                }
                else
                {
                    colors[x + yOffset + zOffset] = 
                    	new Color(dicomImages[z][y][x], dicomImages[z][y][x], dicomImages[z][y][x], 1);
                }
            }
        }
    }

    // 将颜色值复制到纹理
    texture.SetPixels(colors);

    // 将更改应用到纹理，然后将更新的纹理上传到 GPU
    texture.Apply();

    //保存Json信息
    SaveJson(personalDicomData);

    // 将纹理保存到 Unity 项目
    AssetDatabase.CreateAsset(texture, "Assets/MPR/Texture3ds/3DPixel.asset");
}

编辑Shader以在平面上显示切片

Shader主要是在git项目unity-valume-rendering上修改过来的，将原来需要循环计算的部分去掉只进行对应位置的渲染，将原来Color中的HU值还原并依据公式进行调窗运算。以下是修改部分的代码。
sample_volume方法的作用是根据位置p得出该点对应的颜色值。_WindowWidth是要设置的窗宽减去1获取到的值，_WindowLocation是要设置的窗位减去0.5获取到的值，之所以没有放在Shader中进行该运算是因为shader里的运算能省则省，不然会拖慢运行速度。_BasicLight是为了防止画面过暗加上的一个基础值，_Intensity是光强。
在frag方法中的_CubePosX、_CubeScaleX等指期望texture3d所在的坐标和大小。
从Shader中创建Material并赋给平面物体，当平面物体上的点和期望Texture3d所在位置重合时，平面上对应点会显示对应灰度值，这样就形成了切片。后续即可使用摄像机和RenderTexture组合用于获取对应切片图像。

float sample_volume(float3 uv, float3 p)
  {
      //计算颜色值（0~1）
      float origin = (tex3D(_Volume, uv).a * 2 + -1) * tex3D(_Volume, uv).r * 10000;
      origin = (origin - _WindowLocation) / _WindowWidth + 0.5; 
      if (origin < 0)
          origin = 0;
      else if (origin > 1)
          origin = 1;
      float v = origin == 0 ? 0 : (_BasicLight + origin * _Intensity);
      
      //从unity-valume-rendering项目中复制来的，个人见解是限制显示范围
      float3 axis = mul(_AxisRotationMatrix, float4(p, 0)).xyz;
      axis = get_uv(axis);
      float min = step(0, axis.x) * step(0, axis.y) * step(0, axis.z);
      float max = step(axis.x, 1) * step(axis.y, 1) * step(axis.z, 1);

      return v * min * max;
  }

  fixed4 frag(v2f i) : SV_Target
  {
    UNITY_SETUP_INSTANCE_ID(i);

    float4 dst = float4(0, 0, 0, 0);
    //根据点i的世界坐标，按照当前Cube坐标和大小换算回1*1*1的Texture3d中的点坐标
    float3 p = float3((i.world.x - _CubePosX) / _CubeScaleX, (i.world.y - _CubePosY) / _CubeScaleY, (i.world.z - _CubePosZ) / _CubeScaleZ);
    float3 uv = get_uv(p);
    float v = sample_volume(uv, p);
    float4 src = float4(v, v, v, v);
    src.a *= 0.5;
    src.rgb *= src.a;
    dst = (1.0 - dst.a) * src + dst;

    return saturate(dst) * _Color;
  }