利用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.PixelSpacing
    SliceThickness=RefDs.SliceThickness
    ImageOrientationX=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 info
    img_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)=0
	
    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]
    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 line
    a1,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+b1
    yi=a2*zi+b2
    ax.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-ImagePosition
    differ_x=np.dot(differ,ImageOrientationX)
    differ_y=np.dot(differ,ImageOrientationY)
    pos_2d=[differ_x/PixelSpacing[0],differ_y/PixelSpacing[1]]

    return pos_2d