STL与Dicom坐标系对应关系（Matlab版）

# STL与Dicom坐标系对应关系（Matlab版）
Dicom是常用的医学影像存储格式，医生及医学影像工作者对此非常熟悉，本文尝试总结STL文件跟Dicom文件间的关系，并附有Matlab实现代码
问题提出
是否曾有过STL的显示与Dicom数据的显示位置不对的情况？如

问题分析
此问题主要是STL坐标系与Dicom数据的坐标系不统一导致的。一般来说，Dicom数据保存了影像设备的原点及方向向量，因此，影像浏览器在导入图像数据时会通过该信息把三维图像转换到世界坐标系的真实位置。但是，由于STL文件的生成最初是基于图像坐标系的，所以生成STL时需要把图像坐标系的坐标点转换为世界坐标系。

以下简略介绍Dicom及STL，最后通过matlab代码对其进行关联。
Dicom
网络上很多博主已经总结了Dicom常用的Tag值和冠状面，矢状面和横断面的解释，详情见Tag值和三视图

关于Dicom坐标系的详细解释，见Roni的博客,以下为一些中文翻译及总结：
Image Plane Module定义了相对于病人身体的图像方向及每个方向的像素（或体素）间距（mm）。

Dicom定义了Reference Coordinates System (RCS)，分别为轴状位（Transverse/Axisplane XY）、冠状位（Coronal/Frontal plane XZ）和矢状位（Sagittal plane YZ）。其中X方向是从右到左，

Y方向是从前到后，

Z方向是从脚到头

目前，我们了解了Dicom格式的方向定义，
[R] -Right-沿着X方向递减；
[L] -Left-沿着X方向递增；
[A] -Anterior- 沿着Y方向递减；
[P] - Posterior- 沿着Y方向递增；[F] -Feet-沿着Z方向递减；
[H] - Head- 沿着Z方向递增。
大家可能会问，如果病人躺着拍和俯着拍不是不一样吗？Dicom标准提供patient position(0018,5100)来确定病人被扫描时的状态：
HFP = Head First-Prone（头部前，俯卧）.
HFS = Head First-Supine（头部前，仰卧）.
FFP = Feet First-Prone（脚部前，俯卧）.
FFS = Feet First-Supine（脚部前，仰卧）
HFDR = Head First-Decubitus Right（头部前，右侧俯卧）.
HFDL = Head First-Decubitus Left（头部前，左侧俯卧）.
FFDR = Feet First-Decubitus Right（脚部前，右侧俯卧）.
FFDL = Feet First-Decubitus Left（脚部前，左侧俯卧）.

STL

一般图像处理都是以图像坐标系进行编辑，即从1（C为0）开始，间隔为1。但是此时获取的STL文件并不能与DICOM数据一一对应，如

通过4*4的旋转矩阵可以使基于图像坐标系的STL旋转平移到目标世界坐标系中，而Dicom的Tag值可以获取旋转矩阵所需的信息，以下为旋转矩阵获取代码：

function [M,R] = TransMatrix(info)
%This function calculates the 4x4 transform matrix from the image
%coordinates to patient coordinates. 
ipp=info.ImagePositionPatient;
iop=info.ImageOrientationPatient;
ps=info.PixelSpacing;
Tipp=[1 0 0 ipp(1); 0 1 0 ipp(2); 0 0 1 ipp(3); 0 0 0 1];
r=iop(1:3);  c=iop(4:6); s=cross(r',c');
R = [r(1) c(1) s(1) 0; r(2) c(2) s(2) 0; r(3) c(3) s(3) 0; 0 0 0 1];
if info.MRAcquisitionType=='3D' % 3D turboflash
    S = [ps(2) 0 0 0; 0 ps(1) 0 0; 0 0 info.SliceThickness 0 ; 0 0 0 1];
else % 2D epi dti
    S = [ps(2) 0 0 0;0 ps(1) 0 0;0 0 info.SpacingBetweenSlices 0;0 0 0 1];
end
T0 = [ 1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
M = Tipp * R * S * T0;
end

以下为STLResult为你分割后的模板矩阵，str_DicomInfo为Dicom文件的路径

STL_X = 0:1:(n_NumX-1);
STL_Y = 0:1:(n_NumY-1);
STL_Z = 0:1:(n_NumZ-1);
fv = isosurface(STL_X,STL_Y,STL_Z,STLResult,0.5);
fCols8bit = zeros(size(fv.faces));
fCols8bit(:,1) = 255;
info = dicominfo(str_DicomInfo);    
[M,R] = TransMatrix(info);
n_VerticesNum = size(fv.vertices);
TempData=ones(4,1);
for num=1:n_VerticesNum
    TempData(1:3,1) = (fv.vertices(num,:))';
    TempResult = (M*TempData)';
    fv.vertices(num,:) = TempResult(1,1:3);
end
stlwrite(str_Output,fv,'FaceColor',fCols8bit);

转换结果
经过4*4矩阵变换的结果为：