(MATLAB)DICOM序列的读取与数字体模重建

---

by HPC_ZY

之前写了一些关于医学影像三维重建（可视化）的文章，不少网友对原始数据那一部分提问。所以写此文，分享DICOM序列重建三维数字体模的方法。

---

DICOM格式

DICOM具体格式与解析方法本文不讲，感兴趣可参考DICOM官网、大佬博客

---

DICOM序列读取

DICOM数据的读取与处理可参考：MATLAB医学DICOM影像读取与预处理
通常情况下DICOM文件的命名与其序列编号顺序一致，但不排除例外（我就遇到过命名与顺序完全不同的数据），如下表

情况	举例（假如一个有五个文件的序列，括号内为序列号）
顺序正确、序列完整	01.DCM（1）02.DCM（2）03.DCM（3）04.DCM（4）05.DCM（5）
顺序错误、序列完整	01.DCM（3）02.DCM（2）03.DCM（1）04.DCM（5）05.DCM（4）
顺序正确、序列缺失	01.DCM（1）02.DCM（2）03.DCM（4）04.DCM（6）05.DCM（7）

1. 顺序数据的处理
   初始化三维矩阵，利用for循环读取，存入每一层即可。如下：

%% 情况一：顺序正确、序列完整
N = 5;
DCM = zeros(512,512,N);
% 读取
for n = 1:N
filename = ['0',num2str(n),'.DCM']；% 此处使用自己的文件名
DCM(:,:,n) = dicomread(filename ); 
end

2. 乱序数据的处理
   在情况一的基础上，读取文件(0x0020,0x1041)Slice Location数据元中真实空间位置，进行重排序。如下：

%% 情况二：顺序错误、序列完整
N = 5;
DCM = zeros(512,512,N);
IDX = zeros(1,N);

% 读取
for n = 1:N
filename = ['0',num2str(n),'.DCM']；% 此处使用自己的文件名
DCM(:,:,n) = dicomread(filename); 
info = dicominfo(filename);
IDX(n) = info.SliceLocation;
end

% 重排序
[~,idx] = sord(IDX);
DCM = DCM(:,:,idx);

3. 缺失数据的处理
   （一种解决方案）可在情况二的基础上，利用插值算法补齐缺失。（不举例）

---

数字体模重建

其本质就是利用插值算法，将分辨率更低的方向进行重采样。
注：这一步并非必要，仅仅作为图像处理时有上节的DCM就足够。如果后续算法有实际的物理意义，或希望每个体素xyz方向尺度一致，才需要进行一下处理。

此方法代码示例,

% 转为double才能interp
DCM = double(DCM); % 用前面方法读入的 DCM

% info 是上述 info = dicominfo(filename)获取的
ps = info.PixelSpacing; % 我的数据是0.5mm;
sps = info.SpacingBetweenSlices; % 我的数据是5mm;

% 原始网格
[R,C,D] = size(DCM);
[X1,Y1,Z1] = meshgrid(1:R,1:C,1:D);

% 采样网格
dz = ps/sps;
[X2,Y2,Z2] = meshgrid(1:R,1:C,1:dz:D);

% 插值
newDCM = interp3(X1,Y1,Z1,DCM,X2,Y2,Z2);

上述重采样方法虽然简单快速，但其结果不具备物理意义。所以我们通常使用以下重采样方法：

$model(i,j,k)=$
$$(1-a)\cdot (1-b)\cdot (1-c)\cdot dcm(x,y,z)+(1-a)\cdot (1-b)\cdot c\cdot dcm(x,y,z+1)$$$$+a\cdot (1-b)\cdot (1-c)\cdot dcm(x+1,y,z)+a\cdot (1-b)\cdot c\cdot dcm(x+1,y,z+1)$$$$+(1-a)\cdot b\cdot (1-c)\cdot dcm(x,y+1,z)+(1-a)\cdot b\cdot c\cdot dcm(x,y+1,z+1)$$$$+a\cdot b\cdot (1-c)\cdot dcm(x+1,y+1,z)+a\cdot b\cdot c\cdot dcm(x+1,y+1,z+1)$$
其中，model为重采样后的CT三维体数据，dcm为原始CT数据。

关于a,b,c，比如model(i,j,k) 的位置在 (5.5, 3.1, 6)，那么有：
x = 5, y = 3, z = 6
a = 0.5, b = 0.1, c = 0

例子，一个尺寸为512·512·83的CT序列，扫描间隔为（0.6523, 0.6523, 2.5000）mm，重采样后尺寸为512·512·315。如下图所示，

% 假设
DCM = double(DCM); % 用前面方法读入的 DCM
scale = 1; % 设置越大，模型越大
info.PixelSpacing = 0.5;
info.SpacingBetweenSlices = 5;

% 原始网格
[R,C,D] = size(DCM);
[X1,Y1,Z1] = meshgrid(1:R,1:C,1:D);

% 采样网格
mR = R*scale;
mC = C*scale;
mD = D*scale*info.SpacingBetweenSlices/info.PixelSpacing;
dx = (R-1)/(mR-1);
dy = (C-1)/(mC-1);
dz = (D-1)/(mD-1);
[X2,Y2,Z2] = meshgrid(1:dx:R,1:dy:C,1:dz:D);

% 插值
newDCM = interp3(X1,Y1,Z1,DCM,X2,Y2,Z2);

---

其他

由于实验数据为内部资料，不提供测试数据与代码