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

2023-10-01 08:07| 来源: 网络整理| 查看: 265

在此感谢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)：ImageOrientationPixel spacing (0028, 0030)：PixelSpacingRow 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

【本文地址】

公司简介

联系我们