# Source code for resqpy.olio.intersection

```"""intersection.py: functions to test whether lines intersect with planes."""

import logging

log = logging.getLogger(__name__)

import numpy as np
from numba import njit  # type: ignore
from typing import Union

import resqpy.olio.vector_utilities as vec

[docs]def line_plane_intersect(line_p, line_v, triangle):
"""Find the intersection of a line with a plane defined by a triangle.

arguments:
line_p (3 element numpy vector): a point on the line
line_v (3 element numpy vector): vector being the direction of the line
triangle ((3, 3) numpy array): three points on the plane (second index is xyz)

returns:
point (3 element numpy vector) of intersection of the line with the plane,
or None if line is parallel to plane
"""

p01 = triangle[1] - triangle[0]
p02 = triangle[2] - triangle[0]

norm = np.cross(p01, p02)  # normal to plane
denom = np.dot(np.negative(line_v), norm)
if denom == 0.0:
return None  # line is parallel to plane
t = np.dot(norm, line_p - triangle[0]) / denom
return line_p + t * line_v

[docs]def line_triangle_intersect(line_p, line_v, triangle, line_segment = False, l_tol = 0.0, t_tol = 0.0):
"""Find the intersection of a line within a triangle in 3D space.

arguments:
line_p (3 element numpy vector): a point on the line
line_v (3 element numpy vector): vector being the direction of the line
triangle ((3, 3) numpy array): three corners of the triangle (second index is xyz)
line_segment (boolean): if True, returns None if intersection is outwith (line_p .. line_p + line_v)
l_tol (float, default 0.0): a fraction of the line length to allow for an intersection to be found
just outside the segment
t_tol (float, default 0.0): a fraction of the triangle size to allow for an intersection to be found
just outside the triangle

returns:
point (3 element numpy vector) of intersection of the line within the triangle,
or None if line is parallel to plane of triangle or intersection with the plane is
outside the triangle
"""

p01 = triangle[1] - triangle[0]
p02 = triangle[2] - triangle[0]

norm = np.cross(p01, p02)  # normal to plane
line_rv = np.negative(line_v)
denom = np.dot(line_rv, norm)
if denom == 0.0:
return None  # line is parallel to plane
lp_t0 = line_p - triangle[0]
t = np.dot(norm, lp_t0) / denom
if line_segment and (t < 0.0 - l_tol or t > 1.0 + l_tol):
return None
u = np.dot(np.cross(p02, line_rv), lp_t0) / denom
if u < 0.0 - t_tol or u > 1.0 + t_tol:
return None
v = np.dot(np.cross(line_rv, p01), lp_t0) / denom
if v < 0.0 - t_tol or u + v > 1.0 + t_tol:
return None

return line_p + t * line_v

[docs]@njit  # pragma: no cover
def line_triangle_intersect_numba(
line_p: np.ndarray,
line_v: np.ndarray,
triangle: np.ndarray,
line_segment: bool = False,
l_tol: float = 0.0,
t_tol: float = 0.0,
) -> Union[None, np.ndarray]:
"""Find the intersection of a line within a triangle in 3D space.

arguments:
line_p (np.ndarray): a point on the line.
line_v (np.ndarray): vector being the direction of the line.
triangle (np.ndarray): shape (3, 3); three corners of the triangle (second index is xyz).
line_segment (bool): if True, returns None if intersection is outwith (line_p .. line_p + line_v).
l_tol (float, default 0.0): a fraction of the line length to allow for an intersection to be found
just outside the segment.
t_tol (float, default 0.0): a fraction of the triangle size to allow for an intersection to be found
just outside the triangle.

returns:
point (np.ndarray) of intersection of the line within the triangle,
or None if line is parallel to plane of triangle or intersection with the plane is
outside the triangle.
"""

p01 = triangle[1] - triangle[0]
p02 = triangle[2] - triangle[0]

norm = np.cross(p01, p02)  # normal to plane
line_rv = np.negative(line_v)
denom = np.dot(line_rv, norm)
if denom == 0.0:
return None  # line is parallel to plane
lp_t0 = line_p - triangle[0]
t = np.dot(norm, lp_t0) / denom
if line_segment and (t < 0.0 - l_tol or t > 1.0 + l_tol):
return None
u = np.dot(np.cross(p02, line_rv), lp_t0) / denom
if u < 0.0 - t_tol or u > 1.0 + t_tol:
return None
v = np.dot(np.cross(line_rv, p01), lp_t0) / denom
if v < 0.0 - t_tol or u + v > 1.0 + t_tol:
return None

return line_p + t * line_v

[docs]def line_triangles_intersects(line_p, line_v, triangles, line_segment = False):
"""Find the intersections of a line within each of a set of triangles in 3D space.

arguments:
line_p (3 element numpy vector): a point on the line
line_v (3 element numpy vector): vector being the direction of the line
triangles ((n, 3, 3) numpy array): three corners of each of the n triangles (final index is xyz)
line_segment (boolean, default False): if True, the line is treated as a finite segment between
p and p + v, and only intersections within the segment are included

returns:
points ((n, 3) numpy array) of intersection points of the line within the triangles,
(nan, nan, nan) where line is parallel to plane of triangle or intersection with the plane is
outside the triangle (or beyond the ends of the segment if applicable)
"""

n = triangles.shape[0]  # number of triangles

p01s = triangles[:, 1, :] - triangles[:, 0, :]  # p01s has shape (n, 3)
p02s = triangles[:, 2, :] - triangles[:, 0, :]  # p02s has shape (n, 3)

norms = np.cross(p01s, p02s)  # normals to planes, shape (n, 3)
line_rv = np.negative(line_v)
denoms = np.dot(norms, line_rv)  # shape (n)
#  if denom == 0.0: return None  # line is parallel to plane
lp_t0s = line_p - triangles[:, 0]  # shape (n, 3)
ts = np.empty(n)
ts.fill(np.nan)
us = ts.copy()
vs = ts.copy()
nz = (denoms != 0.0)
#  np.divide(np.dot(norms, lp_t0s), denoms, out = ts, where = nz)
np.seterr(divide = 'ignore')
np.divide(np.sum(norms * lp_t0s, axis = 1), denoms, out = ts, where = nz)
np.seterr(invalid = 'ignore')
if line_segment:
ts[:] = np.where(np.logical_or(ts < 0.0, ts > 1.0), np.nan, ts)
np.divide(np.sum(np.cross(p02s, line_rv) * lp_t0s, axis = 1), denoms, out = us, where = nz)
np.divide(np.sum(np.cross(line_rv, p01s) * lp_t0s, axis = 1), denoms, out = vs, where = nz)
ts[np.where(np.logical_or(np.logical_or(us < 0.0, us > 1.0), np.logical_or(vs < 0.0, us + vs > 1.0)))] = np.nan

intersects = np.empty((n, 3))
intersects[:] = line_v * np.repeat(ts, 3).reshape((n, 3)) + line_p
np.seterr(all = 'warn')
return intersects

[docs]def intersects_indices(single_line_intersects):
"""Returns a list of the (triangle) indices where a valid intersection has been found for a single line."""

return list(np.where(np.logical_not(np.isnan(single_line_intersects[..., 0])))[0])

[docs]def line_set_triangles_intersects(line_ps, line_vs, triangles, line_segment = False):
"""Find the intersections of each of set of lines within each of a set of triangles in 3D space.

arguments:
line_ps ((c, 3) numpy array): a point on each of c lines
line_vs ((c, 3) numpy array): vectors being the direction of each of the c lines (or 1 common vector)
triangles ((n, 3, 3) numpy array): three corners of each of the n triangles (final index is xyz)
line_segment (boolean, default False): if True, each line is treated as a finite segment between
p and p + v, and only intersections within the segment are included

returns:
points ((c, n, 3) numpy array) of intersections of the lines within the triangles,
(nan, nan, nan) where a line is parallel to plane of triangle or intersection with the plane is
outside the triangle

note:
this function is computationally and memory intensive; it could benefit from parallelisation
"""

c = line_ps.shape[0]  # number of lines
n = triangles.shape[0]  # number of triangles

if line_vs.ndim == 1:  # single common direction vector; for now simply replicate
line_vs = np.repeat(np.expand_dims(line_vs, axis = 0), c, axis = 0)

p01s = triangles[:, 1, :] - triangles[:, 0, :]  # p01s has shape (n, 3)
p02s = triangles[:, 2, :] - triangles[:, 0, :]  # p02s has shape (n, 3)

norms = np.cross(p01s, p02s)  # normals to planes, shape (n, 3)
line_rvs = np.negative(line_vs)  # shape (c, 3)
denoms = np.dot(line_rvs, norms.T)  # shape (c, n)  where denoms == 0.0 line is parallel to plane of triangle
lps_t0s = line_ps.reshape((c, 1, 3)) - triangles[:, 0, :].reshape((1, n, 3))  # shape (c, n, 3)
ts = np.empty((c, n))  # shape (c, n)
ts.fill(np.nan)
us = ts.copy()  # shape (c, n)
vs = ts.copy()  # shape (c, n)
nz = (denoms != 0.0)  # shape (c, n)
# the np.sum() clause implememts a dot product over the xyz axis
np.seterr(divide = 'ignore')
np.divide(np.sum(norms.reshape((1, n, 3)) * lps_t0s, axis = 2), denoms, out = ts, where = nz)  # shape (c, n)
np.seterr(invalid = 'ignore')
if line_segment:
ts[:] = np.where(np.logical_or(ts < 0.0, ts > 1.0), np.nan, ts)
cp02s = p02s.reshape((1, n, 3)).repeat(c, axis = 0)  # shape (c, n, 3)
cl_rvs = line_rvs.reshape((c, 1, 3)).repeat(n, axis = 1)  # shape (c, n, 3)
np.divide(np.sum(np.cross(cp02s, cl_rvs) * lps_t0s, axis = 2), denoms, out = us, where = nz)  # shape (c, n)
cp01s = p01s.reshape((1, n, 3)).repeat(c, axis = 0)  # shape (c, n, 3)
np.divide(np.sum(np.cross(cl_rvs, cp01s) * lps_t0s, axis = 2), denoms, out = vs, where = nz)
ts[:] = np.where(np.logical_or(np.logical_or(us < 0.0, us > 1.0), np.logical_or(vs < 0.0, us + vs > 1.0)), np.nan,
ts)

intersects = np.empty((c, n, 3))
intersects[:] = line_vs.reshape((c, 1, 3)) * np.repeat(ts, 3).reshape((c, n, 3)) + line_ps.reshape((c, 1, 3))
np.seterr(all = 'warn')
return intersects

[docs]def poly_line_triangles_intersects(line_ps, triangles):
"""Find the intersections of each segment of an open poly-line with each of a set of triangles in 3D space.

arguments:
line_ps ((c, 3) numpy array): ordered points on each of c-1 segments of a poly-line
triangles ((n, 3, 3) numpy array): three corners of each of the n triangles (final index is xyz)

returns:
points ((c-1, n, 3) numpy array) of intersections of the line segments within the triangles,
(nan, nan, nan) where a line segment is parallel to plane of triangle or intersection of
the segment with the plane is outside the triangle or beyond the ends of the segment

note:
this function is computationally and memory intensive; it could benefit from parallelisation
"""

return line_set_triangles_intersects(line_ps[:-1], line_ps[1:] - line_ps[:-1], triangles, line_segment = True)

[docs]def distilled_intersects(line_set_intersections):
"""Returns lists of line and triangle indices, and corresponding intersection points.

argument:
line_set_intersections (numpy float array of shape (nl, nt, 3)): where nl is the number of lines,
nt is the number of triangles and the final axis is x,y,z; nan values indicate no intersection;
this array is as returned by the line_set_triangles_intersects() function or the
poly_line_triangles_intersects() function

returns:
(numpy int array of shape (N,), numpy int array of shape (N,), numpy float array of shape (N, 3)):
for N intersections, first array is list of line indices, second is list of triangle indices, the
third array contains the (x, y, z) coordinates of the intersection points
"""

lines, triangles = np.where(np.logical_not(np.isnan(line_set_intersections[:, :, 0])))
return lines, triangles, line_set_intersections[lines, triangles, :]

[docs]def last_intersects(line_set_intersections):
"""From the result of line_set_triangles_intersects(), returns a vector of intersection points, one per line.

argument:
line_set_intersections (numpy float array of shape (nl, nt, 3)): where nl is the number of lines,
nt is the number of triangles and the final axis is x,y,z; nan values indicate no intersection;
this array is as returned by the line_set_triangles_intersects() function or the
poly_line_triangles_intersects() function

returns:
numpy float array (nl, 3): intersection points, where nl is number of lines

notes:
Use this function to force at most one intersection point per line (with a triangulated surface).
Applicable where lines are expected to be very roughly orthogonal to a gently verying untorn
surface, eg. pillar lines intersecting an unfaulted horizon.
If more than one triangle is intersected by a line, the returned point is for the 'last'
triangle intersected by the line (when checking triangles in the order they appear in the
list of triangles).
If no triangles are intersected by a line, the resulting point will be (nan, nan, nan).
"""

assert line_set_intersections.ndim == 3 and line_set_intersections.shape[2] == 3
c = line_set_intersections.shape[0]

pick = np.empty((c, 3))
pick.fill(np.nan)
pair_l, pair_t = np.where(np.logical_not(np.isnan(line_set_intersections[:, :, 0])))
pick[pair_l, :] = line_set_intersections[pair_l, pair_t, :]

return pick

[docs]def triangles_for_line(line_set_intersections, line_index):
"""From the result of line_set_triangles_intersects(), returns a list of intersected triangles for a line.

arguments:
line_set_intersections (numpy float array of shape (nl, nt, 3)): where nl is the number of lines,
nt is the number of triangles and the final axis is x,y,z; nan values indicate no intersection;
this array is as returned by the line_set_triangles_intersects() function or the
poly_line_triangles_intersects() function
line_index (integar): the index of the line for which the intersected triangle list is required

returns:
(numpy 1D int array of size N, numpy 2D array of shape (N, 3)): the first of the pair of arrays
returned is a list of the triangle indices which the given line intersects; the second array is
the list of corresponding intersection points (each x,y,z)

notes:
if the line does not intersect any triangles, both the resulting arrays will have size zero
"""

triangles = np.where(np.logical_not(np.isnan(line_set_intersections[line_index, :, 0])))[0]
return triangles, line_set_intersections[line_index, triangles, :]

[docs]def lines_for_triangle(line_set_intersections, triangle_index):
"""From the result of line_set_triangles_intersects(), returns a list of lines intersecing given triangle.

arguments:
line_set_intersections (numpy float array of shape (nl, nt, 3)): where nl is the number of lines,
nt is the number of triangles and the final axis is x,y,z; nan values indicate no intersection;
this array is as returned by the line_set_triangles_intersects() function or the
poly_line_triangles_intersects() function
triangle_index (integar): the index of the triangle for which the intersecting line list is required

returns:
(numpy 1D int array of size N, numpy 2D array of shape (N, 3)): the first of the pair of arrays
returned is a list of the indices of lines which intersect with the given triangle; the second
array is the list of corresponding intersection points (each x,y,z)

notes:
if no lines intersect the triangle, both the resulting arrays will have size zero
"""

lines = np.where(np.logical_not(np.isnan(line_set_intersections[:, triangle_index, 0])))[0]
return lines, line_set_intersections[lines, triangle_index, :]

[docs]def poly_line_triangles_first_intersect(line_ps, triangles, start = 0):
"""Finds the first intersection of a segment of an open poly-line with any of a set of triangles in 3D space.

arguments:
line_ps ((c, 3) numpy array): ordered points on each of c-1 segments of a poly-line
triangles ((n, 3, 3) numpy array): three corners of each of the n triangles (final index is xyz)
start (int, default 0): the index of the point in line_ps to start searching from

returns:
(int, numpy float array of shape (3,)) where the int is the line segment index where one or more
intersections were found and the floats are the xyz of one intersection of that line segment within
the triangles; if no line segments intersect any of the trianges, (None, None) is returned

note:
this function is computationally and memory intensive; it could benefit from parallelisation
"""

for c in range(start, len(line_ps) - 1):
points = line_triangles_intersects(line_ps[c], line_ps[c + 1] - line_ps[c], triangles, line_segment = True)
if np.all(np.isnan(points)):
continue  # no intersection found
xyz = last_intersects(points.reshape((1, -1, 3))).reshape((3,))
return c, xyz
return (None, None)

[docs]def line_line_intersect(x1, y1, x2, y2, x3, y3, x4, y4, line_segment = False, half_segment = False):
"""Returns the intersection x',y' of two lines x,y 1 to 2 and x,y 3 to 4.

arguments:
x1, y1, x2, y2: coordinates of two points defining first line
x3, y3, x4, y4: coordinates of two points defining second line
line_segment (bool, default False): if False, both lines are treated as unbounded; if True second line
is treated as segment bounded by both end points and first line bounding depends on half_segment arg
half_segment (bool, default False): if True and line_segment is True, first line is bounded only at
end point x1, y1, whilst second line is fully bounded; if False and line_segment is True, first line
is also treated as fully bounded; ignored if line_segment is False (ie. both lines unbounded)

returns:
pair of floats being x, y of intersection point; or None, None if no qualifying intersection

note:
in the case of bounded line segments, both end points are 'included' in the segment
"""

divisor = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)
if divisor == 0.0:
return None, None  # parallel or coincident lines

if line_segment:
t = ((x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4)) / divisor
if t < 0.0 or (not half_segment and t > 1.0):
return None, None
u = -((x1 - x2) * (y1 - y3) - (y1 - y2) * (x1 - x3)) / divisor
if not (0.0 <= u <= 1.0):
return None, None
x = x1 + t * (x2 - x1)
y = y1 + t * (y2 - y1)

else:
a = x1 * y2 - y1 * x2
b = x3 * y4 - y3 * x4
x = (a * (x3 - x4) - (x1 - x2) * b) / divisor
y = (a * (y3 - y4) - (y1 - y2) * b) / divisor

return x, y

[docs]def point_projected_to_line_2d(p, l1, l2):
"""Return the point on the unbounded line passing through l1 & l2 which is closest to point p, in xy plane."""

# note: result should be at closest point even in the presence of mixed xy & z units (?)
# create normal vector to l1, l2
v = l2 - l1
n = np.array((-v[1], v[0]))
# find intersection of p, p->n with l1, l2
pn = np.array(p[:2]) + n
return line_line_intersect(l1[0], l1[1], l2[0], l2[1], p[0], p[1], pn[0], pn[1], line_segment = False)

[docs]def point_snapped_to_line_segment_2d(p, l1, l2):
"""Returns the point on the bounded line segment l1, l2 which is closest to point p, in xy plane."""

if vec.is_obtuse_2d(l1, p, l2):
return l1[:2]
elif vec.is_obtuse_2d(l2, p, l1):
return l2[:2]
else:
return point_projected_to_line_2d(p, l1, l2)
```