目录
项目背景和需求
在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 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;
}
多平面三维重建
如图所示，多平面三维重建即根据影像数据重建出影像横断面、矢状位、冠状位图像，且切面可以进行角度、位置调整。
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;
}
将穿刺针相对模型的位置和切片位置对应
大小对齐可以使用之前获取到的像素间隔和图片间隔，根据图片宽高像素数和图片数来计算实际上texture3d的长宽高，并以此对Shader中的显示的Cube进行大小位置调整，是的模型和切片的大小对齐。
位置对齐。
实现效果
项目实际在Hololens2上运行较为卡顿，约30帧左右，且波动明显。个人感觉Hololens的性能实在有限，如果相关运算可以放到服务器上，Hololens只做显示可能时更好的方案。
写在最后
这篇文章本来写于2023年，因为各种各样的原因一直躺在草稿箱里，今天清理草稿箱时看到。我已经不再从事相关工作内容，也无法知道自己这篇文章的技术水平如何，今天发布出来权当是记录过往的时光，或许也能算作一点微小的贡献。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
