本文共 2998 字，大约阅读时间需要 9 分钟。

---

by HPC_ZY

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

---

DICOM格式

DICOM具体格式与解析方法本文不讲，感兴趣可参考、

---

DICOM序列读取

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 = Nfilename = ['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 = Nfilename = ['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才能interpDCM = 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); % 用前面方法读入的 DCMscale = 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);

---

其他

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

转载地址：http://umoaz.baihongyu.com/