dicom坐标系以及确定两平面交点

这篇文章前先参考前博客
1.利用DICOM文件实现2D与3D体素坐标之间的转换
2.dicom世界观
dicom坐标系统

1.Image Position and Image Orientation

在DICOM标准的Image Plane属性中有两个关键字段：

• Image Position(Patient)-（0020,0032）：一个三元Double数组，分别标记了影像左上角第一个元素在RCS坐标系中的x、y、z三轴坐标；
• Image Orientation(Patient)-（0020,0037）：一个六元Double数组。
  
  （1）欧拉旋转角与余弦矩阵之间的关系和转换方式

yaw = 0.7854;  
pitch = 0.1; 
roll = 0;
dcm = angle2dcm( yaw, pitch, roll )
dcm =

    0.7036    0.7036   -0.0998
   -0.7071    0.7071         0
    0.0706    0.0706    0.9950

2.dicom坐标系下求两个平面的交线

#3d空间的点投影到2d平面
def TO2D(s,ImagePosition,ImageOrientationX,ImageOrientationY,PixelSpacing):
    # 3d空间交线的表达式
    x, y, z = symbols('x, y, z')
    #check s  list out of range
    if s:
        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
    else:
        print("s is ERROR")
 
 #交线在2d平面的边界的交点
def get_2dpos(pos_2d,img_array):
    x, y, z = symbols('x, y, z')
    exp_x = solve([x - pos_2d[0], y - pos_2d[1]], [x, z])[x]

    x1 = solve([x - exp_x, y], [x, y])[x]  # 取出y=0的时候的x值
    x2 = solve([x - exp_x, y - img_array.shape[1]], [x, y])[x]  # 取出y=192的时候的x值
    return x1, x2

def res_img(img_array1,x1, x2,path=None,path_heat=None,PixelSpacing=None):
    #在原图像上画线
    pos = [x1, x2]
    mri=cv2.normalize(np.array(img_array1), None, 0, 255, cv2.NORM_MINMAX)
    mri=np.uint8(np.array(mri))
    # # cv2.imwrite(path, mri_1)  
    cv2.line(mri, ( pos[0],0), (pos[1],img_array1.shape[1]), (255, 0, 0), lineType=cv2.LINE_AA)
    clahe = cv2.createCLAHE(tileGridSize=(1, 1))
    mri = clahe.apply(mri)
    mri, PixelSpacing = resample_sapcing(mri, PixelSpacing)  # resample+回归线
    mri = get_square_crop(mri)  # crop+resample+回归线
    cv2.imwrite(path, mri)  
    return  

def getIntersection(short,long,index=0,path=None):
    img_array1, ImagePosition1, PixelSpacing1, ImageOrientation1, SliceThickness1,normalvector1=short[0],short[1],short[2],short[3],short[4],short[5]
    img_array2, ImagePosition2, PixelSpacing2, ImageOrientation2, SliceThickness2,normalvector2=long[0],long[1],long[2],long[3],long[4],long[5]
    ImageOrientationX1 = ImageOrientation1[0:3]
    ImageOrientationY1 = ImageOrientation1[3:6]

    # 设置 并设置符号变量
    sp.init_printing(use_unicode=True)
    InteractiveShell.ast_node_interactivity = 'all'
    x, y, z = symbols('x, y, z')
    #建立方程组（关键）根据Orientation确定平面的表达式，再利用解方程的思想确定交线的表达式
    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]))
    #2d空间交线的表达式(img1上) 2D点=（Difference_in_X /以X表示的像素大小，Differnce_in_Y /以Y表示的像素大小）
    pos_2d=TO2D(s,ImagePosition1,ImageOrientationX1,ImageOrientationY1,PixelSpacing1)
    x1_1, x1_2=get_2dpos(pos_2d,img_array1)

    # #在原图像上画线
    path1 = path+str(index)+"_img.png"
    path3 = path+str(index)+ "_heatmap.png"
    res_img(img_array1, x1_1, x1_2,path1,path3,PixelSpacing1)
  return