Find and plot the intersection of side and down image of DICOM

plot_intersection.py

def get_orientation(IOP):
    orientations = [float(x) for x in IOP]
    row = orientations[:3]
    col = orientations[3:]
    return np.array(row), np.array(col)
def get_position(IPP):
    position = [float(x) for x in IPP]
    return np.array(position)
def get_Spacing(PS):
    return [float(x) for x in PS]

def find_intersect(dcm1, dcm2, plot=True):
    # plot first dicom
    df = pydicom.dcmread(dcm1)
    r1, c1 = get_orientation(df.ImageOrientationPatient)
    n1 = np.cross(r1, c1)
    o1 = get_position(df.ImagePositionPatient)
    pixelSpacing = get_Spacing(df.PixelSpacing)
    Width = df.Columns
    Height = df.Rows
    print(Width, Height)
    plt.figure(figsize=(Width//100,Height//100))

    # 计算矢量图的4角坐标
    p1 = o1
    p2 = p1 + r1 * pixelSpacing[0] * (Width-1)
    p3 = p2 + c1 * pixelSpacing[1] * (Height-1)
    p4 = p1 + c1 * pixelSpacing[1] * (Height-1)

    for p in [p1, p2,p3,p4]:  
        # 到方向向量的投影  
        xmap, ymap = np.dot(p - o1, r1) / np.linalg.norm(r1), np.dot(p - o1, c1) / np.linalg.norm(c1)  

        # 转为pixel坐标  
        xmap /= pixelSpacing[0]  
        ymap /= pixelSpacing[1]  
        if plot:
            plt.plot(int(xmap), int(ymap), "o",color='red', ms=10) 
            plt.imshow(df.pixel_array, cmap=plt.cm.gray)  

    # plot second dicom
    df2 = pydicom.dcmread(dcm2)
    r2, c2 = get_orientation(df2.ImageOrientationPatient)
    n2 = np.cross(r2, c2)
    o2 = get_position(df2.ImagePositionPatient)
    pixelSpacing2 = get_Spacing(df2.PixelSpacing)
    Width2 = df2.Columns
    Height2 = df2.Rows

    # 计算矢量图的4角坐标
    p1 = o2
    p2 = p1 + r2 * pixelSpacing2[0] * (Width2-1)
    p3 = p2 + c2 * pixelSpacing2[1] * (Height2-1)
    p4 = p1 + c2 * pixelSpacing2[1] * (Height2-1)

    dv1 = np.dot(p1-o1, n1)
    dv2 = np.dot(p2-o1, n1)
    dv3 = np.dot(p3-o1, n1)
    dv4 = np.dot(p4-o1, n1)
    iscross12 = dv1*dv2<0
    iscross23 = dv2*dv3<0
    iscross34 = dv3*dv4<0
    iscross41 = dv4*dv1<0

    cps = []
    for iscross,(pv,pp),(ddv,ddp) in zip([iscross12, iscross23, iscross34, iscross41],  
                [(p1,p2),(p2,p3),(p3,p4),(p4,p1)], [(dv1,dv2),(dv2,dv3),(dv3,dv4),(dv4,dv1)]):  
        if not iscross:  
            continue  
        cp = [pv[i] + (pp[i] - pv[i]) * ddv / (ddv - ddp) for i in range(3)]  
        cps.append(np.array(cp))  

    assert len(cps) == 2  
    coords = []  
    coords_img = []  
    for cp in cps:  
        xmap = np.dot(cp-o1, r1) / np.linalg.norm(r1)  
        ymap = np.dot(cp-o1, c1) / np.linalg.norm(c1)  
        coords.append(np.array(xmap, ymap))  
        coords_img.append(np.array([int(xmap/pixelSpacing[0]), int(ymap/pixelSpacing[1])]))  
    print(coords_img)
    if plot:
        plt.plot((coords_img[0][0], coords_img[1][0]), (coords_img[0][1], coords_img[1][1]))  
        plt.show()  

Reference: https://tianchi.aliyun.com/forum/post/115693