利用DICOM文件实现2D与3D体素坐标之间的转换

在此感谢Manjunath KN的回答，参考链接如下：https://www.researchgate.net/post/How_to_create_a_simple_project_to_convert_DICOM_images_to_3d_images

2D to 3D

1.需要的信息
将2D坐标转换为3D坐标需要以下DICOM信息2D：

input point (x,y)：原始的2D坐标
Image position patient (0020, 0032)：ImageOrientation
Pixel spacing (0028, 0030)：PixelSpacing
Row vector and column vectors (0020,0037)：ImageOrientation

DICOM 提取示例代码如下：

def getinfo(img_file):RefDs = dicomio.read_file(img_file)# print(dir(RefDs))  #查看dicom文件的属性img_array = RefDs.pixel_array# indexes are z,y,x ImagePosition =np.array(RefDs.ImagePositionPatient)ImageOrientation=np.array(RefDs.ImageOrientationPatient)PixelSpacing =RefDs.PixelSpacingSliceThickness=RefDs.SliceThicknessImageOrientationX=ImageOrientation[0:3]ImageOrientationY=ImageOrientation[3:6]#z轴(X与Y的叉积)normalvector=np.cross(ImageOrientationX,ImageOrientationY)return img_array,normalvector,ImagePosition,PixelSpacing,ImageOrientationX,ImageOrientationY

2.转换过程

Voxel (x, y, z) = Image Plane Position + Row change in X + column change in Y
Where Row change in X = Row vector * Pixel size in X direction * 2D Point location in X direction
Column change in Y = Column vector * Pixel size in Y direction * 2D Point location in Y direction

示例：
Image position patient = (100,100,50) #ImagePositionPatient
row vector=(1,0,0) #ImageOrientationX
column vector=(0,1,0) #ImageOrientationY
pixel size=(0.5,0.5) #PixelSpacing
2D point=(5,6)

Voxel (x,y,z) = (100,100, 50) + (1,0,0) * 0.5 *5 + (0,1,0)0.5 6 = (100,100,50) + (2.5,0,0)+ (0,3,0) = (102.5, 103, 50)

3.将2d平面转换到3d空间原理
这里主要是解释下文代码中 “建立方程组：dx(x-a)+dy(y-b)+dz(z-c)=0”的原理：

（1）dx(x-a)+dy(y-b)+dz(z-c)=0
我们在高中学过空间平面的几种表达方式，比如有点法式、一般式等等。其中点法式最常用的就是A(x-x0)+B(y-y0)+C(z-z0)=0。这里的含义很好理解，假设向量vN=(A,B,C)是平面的法向量，然后取平面任一点p0=(x0,y0,z0)【带入上面式子等于零，说明该点在平面上】。所谓的空间平面，就是任一点pX与点p0链接的直线(pX,p0)与法向量vN内积和等于零。也就是A(x-x0)+B(y-y0)+C(z-z0)=0。这里做了简化，即写成dx(x-a)+dy(y-a)+dz(z-a)=0的意思是直接取(dx,dy,dz)就是法向量的x，y，z三个方向的分量。相当于是吧法向量(A,B,C）进行了归一化处理。
（2）ImageOrientationX即X轴的方向，ImageOrientationY即Y轴方向
图像平面的长，宽方向的方向余弦向量其实就是归一化后的方向向量，叉乘获得的就是平面的法向量，每个分量就是dx,dy,dz。
（3）关于dicon中的左边系可参考博客
https://blog.csdn.net/zssureqh/article/details/61636150

4.两个2d平面投影到3d空间并求得交线

def getIntersection(f1=None,f2=None,path=None):#get infoimg_array1,normalvector1, ImagePosition1,PixelSpacing1,ImageOrientationX1,ImageOrientationY1= getinfo(f1)img_array2,normalvector2, ImagePosition2,PixelSpacing2,ImageOrientationX2,ImageOrientationY2 = getinfo(f2)# 设置 并设置符号变量sp.init_printing(use_unicode=True)InteractiveShell.ast_node_interactivity = 'all'x, y, z = symbols('x, y, z')#建立方程组#dx(x-a)+dy(y-b)+dz(z-c)=0z1=[normalvector1[0] * (x - ImagePosition1[0]) + normalvector1[1] * (y - ImagePosition1[1]) ] /normalvector1[2] + ImagePosition1[2]z2 = [normalvector2[0] * (x - ImagePosition2[0]) + normalvector2[1] * (y - ImagePosition2[1])] / normalvector2[2] + \ImagePosition2[2]eq=[normalvector1[0] * (x - ImagePosition1[0]) + normalvector1[1] * (y - ImagePosition1[1]) + normalvector1[2] * (z - ImagePosition1[2]),\normalvector2[0] * (x - ImagePosition2[0]) + normalvector2[1] * (y - ImagePosition2[1]) + normalvector2[2] * (z - ImagePosition2[2])]#解方程s = list(linsolve(eq, [x, y]))# show_3d(normalvector1, ImagePosition1,normalvector2,ImagePosition2,s)

在3d空间展示平面

def show_3d(normalvector1,ImagePosition1,normalvector2,ImagePosition2,s):# 3d空间交线的表达式x1_3d = s[0][0]y1_3d = s[0][1]x, y, z = symbols('x, y, z')fig = plt.figure()ax = Axes3D(fig)# 生成x,y的网格数据X = np.arange(-256, 256, 1)# Y = np.arange(-4, 4, 0.25)Y = np.arange(-256, 256, 1)X, Y = np.meshgrid(X, Y)Z1=(normalvector1[0]* (X - ImagePosition1[0]) + normalvector1[1]* (Y - ImagePosition1[1]))/ normalvector1[2] + ImagePosition1[2]Z2 = (normalvector2[0] * (X - ImagePosition2[0]) + normalvector2[1] * (Y - ImagePosition2[1])) / normalvector2[2] +  ImagePosition2[2]# show linea1,a2=x1_3d.coeff(z),y1_3d.coeff(z)b1,b2=x1_3d.coeff(z,0),y1_3d.coeff(z,0)zi = np.linspace(-250, 256, 100)xi=a1*zi+b1yi=a2*zi+b2ax.plot3D(xi,yi,zi)ax.plot_surface(X, Y, Z1, rstride=1, cstride=1, color='r')#cmap='rainbow'ax.plot_surface(X, Y, Z2, rstride=1, cstride=1, color="g")plt.show()

3D to 2D

1.原理
Conversion from 3D to 2D
Difference = Voxel – Image plane Pos;
Difference_in_X = Difference.Innerproduct(Row vector)
Difference_in_Y = Difference.Innerproduct(Col vector);
2D Point = (Difference_in_X/pixel size in X, Differnce_in_Y/pixel size in Y)

这是所涉及的基本逻辑。上述方程组也可以通过矩阵乘法来实现。请参阅第410页，DICOM，章节PS 3.3，2011年版本，或者在文件http://dicom.nema.org/medical/dicom/current/output/pdf/part03.pdf中搜索标签0020,0037

2.代码
（1）3d 空间的交线投影到对应的2D平面

def TO2D(s,ImagePosition,ImageOrientationX,ImageOrientationY,PixelSpacing):# s为3d空间交线的表达式x, y, z = symbols('x, y, z')x1_3d = s[0][0]y1_3d = s[0][1]pos=[x1_3d,y1_3d,z]differ=pos-ImagePositiondiffer_x=np.dot(differ,ImageOrientationX)differ_y=np.dot(differ,ImageOrientationY)pos_2d=[differ_x/PixelSpacing[0],differ_y/PixelSpacing[1]]return pos_2d