import typing
import warnings
import numpy as np
import numba as nb
warnings.filterwarnings(
"ignore", ".*type 'reflected list'.*", category=DeprecationWarning
)
@nb.jit(nopython=True, cache=True)
def fit_line(coords: np.ndarray) -> typing.Tuple[np.ndarray, np.ndarray]:
"""Fit coordinates to a straight line.
Args:
coords: Array of coordinates, the columns are the the cartesian coordinates
Returns:
: A tuple of two one-dimensional arrays.
* One point in the straight line
* Direction vector
Notes:
See
https://math.stackexchange.com/questions/1611308/best-fit-line-with-3d-points
"""
if coords.shape[0] == 2:
unnormalised_direction_vector = coords[1, :] - coords[0, :]
return coords[0, :], unnormalised_direction_vector / np.linalg.norm(
unnormalised_direction_vector,
ord=2,
)
else:
coords = np.ascontiguousarray(coords)
# Average coordinates is a point of the line
coords_avg = np.sum(coords, 0) / coords.shape[0]
# Now we need to find a direction vector
coords_centered = coords - coords_avg
covariance = np.dot(coords.T, coords_centered)
eigen_values, eigen_vectors = np.linalg.eigh(covariance)
direction_vector = eigen_vectors[:, np.argmax(eigen_values)]
return coords_avg, direction_vector
@nb.jit(nopython=True, cache=True)
def compute_distance_to_line(
coords: np.ndarray,
point: np.ndarray,
direction_vector: np.ndarray,
) -> float:
"""Compute the normalised distance between coordinates and a straight line.
Args:
coords: Array of coordinates, the columns are the the cartesian coordinates
point: one point belonging to the straight line
direction_vector: Normalised direction vector of the straight line
Returns:
Square root of the average of the squared distances from a straight line
"""
# Project points in `coords` to line defined by `point` and `direction_vector`
coords_from_point = coords - point
scalar_projections = np.dot(coords_from_point, direction_vector)
projections = scalar_projections.reshape((-1, 1)) * direction_vector
# Compute average of the squared distances from a straight line
residuals = coords_from_point - projections # "distance"
return np.sqrt(np.sum(np.square(residuals)) / coords.shape[0])
@nb.jit(nopython=True, cache=True)
def compute_distances_between_two_lines(
coords1: np.ndarray,
direction_vector1: np.ndarray,
coords2: np.ndarray,
direction_vector2: np.ndarray,
) -> np.floating:
"""Compute the distance between two lines.
Args:
coords1: cartesian coordinates of a point belonging to the first straight line
direction_vector1: normalised direction vector of the first straight line
coords2: cartesian coordinates of a point belonging to the second straight line
direction_vector2: normalised direction vector of the second straight line
returns:
Distance between the two lines.
"""
normal_vector = np.cross(direction_vector1, direction_vector2)
if np.all(normal_vector == 0): # the two lines are parallel
return np.linalg.norm(np.cross(coords2 - coords1, direction_vector1), ord=2)
else:
return abs(np.dot(coords2 - coords1, normal_vector))
@nb.jit(nopython=True, cache=True)
def compute_distance_to_z_axis(
coords: np.ndarray,
point: np.ndarray,
direction_vector: np.ndarray,
) -> np.floating:
"""Compute the distance of a straight line to the z-axis.
Args:
coords: Array of coordinates. This is not used.
point: one point belonging to the straight line
direction_vector: Normalised direction vector of the straight line
Returns:
Shortest distance between a line and the z-axis.
"""
return compute_distances_between_two_lines(
coords1=point,
direction_vector1=direction_vector,
coords2=np.zeros(shape=3),
direction_vector2=np.array([0, 0, 1]),
)
[docs]def compute_angle_between_line_and_plane(
line_direction_vector: np.ndarray,
plane_normal_vector: np.ndarray,
) -> float:
"""Compute the angle between a line and a plane.
Args:
line_direction_vector: Normalised direction vector of the straight line
plane_normal_vector: Normalised normal vector of the plane
Returns:
Angle between a line and a plane, in degree.
"""
return (np.arccos(np.dot(plane_normal_vector, line_direction_vector))) * (
180.0 / np.pi
)
[docs]def compute_yz_angle(
coords: np.ndarray,
point: np.ndarray,
direction_vector: np.ndarray,
) -> float:
"""Compute the angle between a line and the :math:`y`-:math:`z` plane.
Args:
coords: Array of coordinates. This is not used.
point: one point belonging to the straight line. This is not used.
direction_vector: Normalised direction vector of the straight line
Returns:
Angle between a line and the :math:`y`-:math:`z` plane, in degree.
"""
return compute_angle_between_line_and_plane(
line_direction_vector=direction_vector,
plane_normal_vector=np.array([1, 0, 0]),
)
[docs]def compute_xz_angle(
coords: np.ndarray,
point: np.ndarray,
direction_vector: np.ndarray,
) -> float:
"""Compute the angle between a line and the :math:`x`-:math:`z` plane.
Args:
coords: Array of coordinates. This is not used.
point: one point belonging to the straight line. This is not used.
direction_vector: Normalised direction vector of the straight line
Returns:
Angle between a line and the :math:`x`-:math:`z` plane, in degree.
"""
return compute_angle_between_line_and_plane(
line_direction_vector=direction_vector,
plane_normal_vector=np.array([0, 1, 0]),
)
