P = np.random.random((10,3))

Ax_support = array([3,2,1])
Ax_direction = array([1,2,3])

P_centered = P - Ax_support

P_projected = P_centered.dot(np.linalg.pinv(
             Ax_direction[np.newaxis, :])).dot(Ax_direction[np.newaxis, :])

distances = np.sqrt(((P_centered - P_projected) ** 2).sum(axis=1))