"""Contains functions for basic geometry with lines that are given by points."""
import numpy as np
from scipy.spatial import cKDTree
from torx.autodoc_decorators_m import autodoc_function
[docs]
@autodoc_function
def intersect_line_segments(x1, y1, x2, y2, x3, y3, x4, y4):
"""
Intersect two line segments that are given by their start and end points.
Here 1 and 2 belongs to the first segment, 3 and 4 to the second. Returns
the intersection point if it exists, otherwise a NaN point.
Note
----
Works in 2D.
"""
denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)
# We need to eliminate the sign of denom because in the condition below
# the <= would have to be switched to >= if denom is negative. We avoid
# this by multiplying with the sign of denom here. Note that this does
# not change the result since a division by denom is required anyway.
t = ((x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4)) * np.sign(denom)
u = ((x1 - x3) * (y1 - y2) - (y1 - y3) * (x1 - x2)) * np.sign(denom)
denom *= np.sign(denom)
# Checks for intersection can be made without division
if(t <= denom and t >= 0 and u <= denom and u >= 0 and denom != 0):
px = x1 + (t / denom) * (x2 - x1)
py = y1 + (t / denom) * (y2 - y1)
return [px, py]
else:
return [np.nan, np.nan]
[docs]
@autodoc_function
def intersect_lines(l1x, l1y, l2x, l2y, intersec_mode:str="first"):
"""
Intersect two lines that are given by points on the curve.
The intersection is performed by checking all combinations of line segments.
Returns the intersection point if it exists, otherwise a NaN point.
Depending on the intersec_mode parameter, different intersection points
are returned. Current options are "first", "last", "all" and "closest"
(w.r.t. first point of l1).
Default is "first". Note that "first" and "closest" are different in the
case of multiple intersections within a segment, while for a single
intersection they return the same.
Note
----
Works in 2D.
"""
p_inter_collection = []
for i in range(l1x.size - 1):
for j in range(l2x.size - 1):
p_inter = intersect_line_segments(l1x[i ], l1y[i ], \
l1x[i + 1], l1y[i + 1], \
l2x[j ], l2y[j ], \
l2x[j + 1], l2y[j + 1])
if not (np.isnan(p_inter[0]) and np.isnan(p_inter[1])):
if intersec_mode == "first":
return p_inter
p_inter_collection.append(p_inter)
if(not p_inter_collection):
return [np.nan, np.nan]
if(len(p_inter_collection) == 1):
return p_inter_collection[0]
if(intersec_mode=="all"):
return p_inter_collection
if(intersec_mode=="last"):
return p_inter_collection[-1]
if(intersec_mode=="closest"):
p_inter_collection = np.array(p_inter_collection)
dist = np.sqrt((p_inter_collection[:, 0] - l1x[0])**2 + \
(p_inter_collection[:, 1] - l1y[0])**2)
i_closest = np.argmin(dist)
return p_inter_collection[i_closest, :]
[docs]
@autodoc_function
def compute_arc_length(r_points: np.ndarray,
z_points: np.ndarray,
phi_points: np.ndarray=None,
skipna: bool=False):
"""
Compute the arc length of a line given by points via triangularization.
Note
----
The skipna implementation is useful but not rigorous. It can handle
single block of NaNs at the beginning or the end of the array. It can also
handle a single block of values surrounded by NaNs. It will produce
inaccurate results for anything more complex involving NaNs.
"""
if phi_points is None:
return _arc_length_2d(r_points, z_points, skipna)
else:
return _arc_length_3d(r_points, z_points, phi_points, skipna)
def _arc_length_2d(r_points, z_points, skipna):
"""Compute the arc length of a line via triangularization in 2D."""
# Only perform the cumulative sum where the coordinate values are not NaNs
if skipna:
isnan = np.isnan(r_points)
assert np.all(np.isnan(z_points) == np.isnan(r_points)), \
"Inconsistent NaN values between 'r_points' and 'z_points'."
assert len(~isnan) > 1, \
"Cannot compute arc length with <2 points remaining."
r_points = r_points[~isnan]
z_points = z_points[~isnan]
# Perform triangularization and add zero value at the beginning
r_diff = np.diff(r_points)
z_diff = np.diff(z_points)
x_vals = np.cumsum(np.sqrt(r_diff**2 + z_diff**2))
x_vals = np.insert(x_vals, 0, 0.0)
# Put NaNs into to the returned array, except for the NaN next to the values
if skipna:
x_vals_full = np.full(len(isnan), np.nan)
x_vals_full[~isnan] = x_vals
else:
x_vals_full = x_vals
return x_vals_full
def _arc_length_3d(r_points, z_points, phi_points, skipna):
"""Compute the arc length of a line via triangularization in 3D."""
# Only perform the cumulative sum where the coordinate values are not NaNs
if skipna:
isnan = np.isnan(r_points)
assert (
np.all(np.isnan(z_points) == np.isnan(r_points)) and
np.all(np.isnan(r_points) == np.isnan(phi_points))
), \
"Inconsistent NaN values between 'phi_points', 'r_points' \
and 'z_points'."
assert len(~isnan) > 1, \
"Cannot compute arc length with <2 points remaining."
r_points = r_points[~isnan]
z_points = z_points[~isnan]
phi_points = phi_points[~isnan]
# Perform triangularization and add zero value at the beginning
r_diff = np.diff(r_points)
z_diff = np.diff(z_points)
phi_diff = np.diff(phi_points)
# Use the toroidal coordinates, i.e. phi_diff = r_mean * dphi
r_stag = r_points[:-1] + r_diff
x_vals = np.cumsum(np.sqrt(r_diff**2 + z_diff**2 + (r_stag*phi_diff)**2))
x_vals = np.insert(x_vals, 0, 0.0)
# Put NaNs into to the returned array, except for the NaN next to the values
if skipna:
x_vals_full = np.full(len(isnan), np.nan)
x_vals_full[~isnan] = x_vals
else:
x_vals_full = x_vals
return x_vals_full
[docs]
@autodoc_function
def order_points(x_points, y_points, start_idx: int=0, k: int=20):
"""
Order the x, y points using a nearest-neighbor algorithm.
Starting from the point with index start_idx and looking for up to k
nearest neighbors. Falls back to a linear search if no unvisited neighbor
is found among the k nearest neighbors.
"""
pts = np.column_stack([x_points, y_points])
n_pts = x_points.size
tree = cKDTree(pts)
visited = np.zeros(n_pts, dtype=bool)
order = np.empty(n_pts, dtype=int)
order[0] = start_idx
cur = start_idx
visited[cur] = True
for i in range(1, n_pts):
# Get up to k nearest neighbors
dists, idxs = tree.query(pts[cur], k=min(k, n_pts))
# Handle k=1
idxs = np.atleast_1d(idxs)
# Find first unvisited neighbor
unvisited = idxs[~visited[idxs]]
if unvisited.size > 0:
nxt = unvisited[0]
else:
# Fallback: linear search for nearest unvisited
d_all = np.linalg.norm(pts - pts[cur], axis=1)
d_all[visited] = np.inf
nxt = np.argmin(d_all)
order[i] = nxt
visited[nxt] = True
cur = nxt
# Apply order
pts = pts[order]
return pts[:, 0], pts[:, 1]
