def point_plane_intersection(plane_normal, plane_point, line_a, line_b,
    epsilon=1e-06):
    """
    Determines the intersection point of a plane defined by a point and a normal vector and the line
    defined by line_a and line_b. All should be numpy arrays.

    :param plane_normal:
    :param plane_point:
    :param line_a:
    :param line_b:
    :param epsilon:
    :return:
    """
    line_direction = line_b - line_a
    if abs(plane_normal.dot(line_direction)) < epsilon:
        return None
    delta = line_a - plane_point
    projection = -plane_normal.dot(delta) / plane_normal.dot(line_direction)
    return delta + projection * line_direction + plane_point
