alpha.py

import numpy as np

def get_camera_center(_P):
    return -np.linalg.inv(_P[:3, :3]) @ _P[:3, 3]


def angle_between_plane_and_line(line_direction: np.ndarray, point1, point2, point3):
    # Calculate the normal vector of the plane
    normal_vector = np.cross(point2-point1, point3-point1)
    normal_vector /= np.linalg.norm(normal_vector)  # normalize normal vector
    line_direction /= np.linalg.norm(line_direction) # normalize line direction

    # Calculate the angle between the normal vector of the plane and the line direction
    angle_rad = np.arcsin(np.dot(normal_vector, line_direction))
    angle_deg = np.degrees(angle_rad)  # optional
    return angle_deg
  
  
# Two example projection matrices from real dicom scans at 0 and 90 degrees
P0 = np.array([
    -1.491806664, -6.0116547806, -0.0032140855, 497.8171012864,
    -0.7840045398, 0.0883435020, 6.1466040404, 486.4775230631,
    -0.0015992436, 0.0001893259, 0.0000025475, 1.0000000000
]).reshape(3, 4)
P1 = np.array([
    5.944656294, -1.7327766731, -0.0146051869, 483.7323692488,
    -0.1118475300, -0.7842383349, 6.1431458554, 492.9361593238,
    -0.0002548330, -0.0015890854, -0.0000018072, 1.0000000000
]).reshape(3, 4)

# construct 0-degree epipolar plane from camera positions and isocenter
C0, C1, iso_center = get_camera_center(P0), get_camera_center(P1), np.zeros(3)

# inverse camera coordinate transform M (invKR)
M0, M1 = -np.linalg.inv(P0[:3, :3]), -np.linalg.inv(P1[:3, :3])

# calculate alpha for a point in view P0 at detector coordinate (0, 0)
point2d = np.array([0, 0, 1])  # homogenous coordinate
p_dir = M0 @ point2d
alpha = angle_between_plane_and_line(p_dir, C0, C1, iso_center)

# Cone Beam angle can be verified through imaging systems geometry
print(f"angle to isocenter plane: {alpha}")
sdd, ds = 1164, 300  # source detector distance and detector side length in mm
print(f"cone beam angle: {np.degrees(np.arcsin((ds/2)/sdd))}")

# Output
# angle to isocenter plane: 7.199223997233601
# cone beam angle: 7.404066525634113