Source code for resqpy.olio.vector_utilities

"""Utilities for working with 3D vectors in cartesian space.

note: some of these functions are redundant as they are provided by built-in numpy operations.

a vector is a one dimensional numpy array with 3 elements: x, y, z.
some functions accept a tuple or list of 3 elements as an alternative to a numpy array.
"""

import logging

log = logging.getLogger(__name__)

import math as maths
import numpy as np
import numba  # type: ignore
from numba import njit  # type: ignore
from typing import Tuple, Optional


[docs]def radians_from_degrees(deg): """Converts angle from degrees to radians.""" return np.radians(deg)
[docs]def degrees_from_radians(rad): """Converts angle from radians to degrees.""" return np.degrees(rad)
[docs]def zero_vector(): """Returns a zero vector [0.0, 0.0, 0.0].""" return np.zeros(3)
[docs]def v_3d(v): """Returns a 3D vector for a 2D or 3D vector.""" assert 2 <= len(v) <= 3 if len(v) == 3: return v v3 = np.zeros(3) v3[:2] = v return v3
[docs]def add(a, b): # note: could just use numpy a + b facility """Returns vector sum a+b.""" a = np.array(a) b = np.array(b) assert a.size == b.size return a + b
[docs]def subtract(a, b): # note: could just use numpy a - b facility """Returns vector difference a-b.""" a = np.array(a) b = np.array(b) assert a.size == b.size return a - b
[docs]def elemental_multiply(a, b): # note: could just use numpy a * b facility """Returns vector with products of corresponding elements of a and b.""" a = np.array(a) b = np.array(b) assert a.size == b.size return a * b
[docs]def amplify(v, scaling): # note: could just use numpy a * scalar facility """Returns vector with direction of v, amplified by scaling.""" v = np.array(v) return scaling * v
[docs]def unit_vector(v): """Returns vector with same direction as v but with unit length.""" assert 2 <= len(v) <= 3 v = np.array(v, dtype = np.float64) norm = np.linalg.norm(v) if norm == 0.0: return v v = v / norm return v
[docs]@njit def unit_vector_njit(v): # pragma: no cover """Returns vector with same direction as v but with unit length.""" norm = np.linalg.norm(v) if norm == 0.0: return v v = v / norm return v
[docs]def unit_vectors(v): """Returns vectors with same direction as those in v but with unit length.""" scaling = np.sqrt(np.sum(v * v, axis = -1)) zero_mask = np.zeros(v.shape, dtype = bool) zero_mask[np.where(scaling == 0.0), :] = True restore = np.seterr(all = 'ignore') result = np.where(zero_mask, 0.0, v / np.expand_dims(scaling, -1)) np.seterr(**restore) return result
[docs]def nan_unit_vectors(v): """Returns vectors with same direction as those in v but with unit length, allowing NaNs.""" assert v.shape[-1] == 3 vf = v.reshape((-1, 3)) nan_mask = np.isnan(vf) restore = np.seterr(all = 'ignore') scaling = np.sqrt(np.sum(vf * vf, axis = -1)) zero_mask = np.zeros(vf.shape, dtype = bool) zero_mask[np.where(scaling == 0.0), :] = True result = np.where(zero_mask, 0.0, vf / np.expand_dims(scaling, -1)) result = np.where(nan_mask, np.nan, result) np.seterr(**restore) return result.reshape(v.shape)
[docs]def unit_vector_from_azimuth(azimuth): """Returns horizontal unit vector in compass bearing given by azimuth (x = East, y = North).""" azimuth = azimuth % 360.0 azimuth_radians = radians_from_degrees(azimuth) result = zero_vector() result[0] = maths.sin(azimuth_radians) # x (increasing to east) result[1] = maths.cos(azimuth_radians) # y (increasing to north) return result # leave z as zero
[docs]def azimuth(v): # 'azimuth' is synonymous with 'compass bearing' """Returns the compass bearing in degrees of the direction of v (x = East, y = North), ignoring z.""" assert 2 <= v.size <= 3 z_zero_v = np.zeros(3) z_zero_v[:2] = v[:2] unit_v = unit_vector(z_zero_v) # also checks that z_zero_v is not zero vector x = unit_v[0] y = unit_v[1] # ignore z component if x == 0.0 and y == 0.0: return 0.0 # arbitrary azimuth of a vertical vector if abs(x) >= abs(y): radians = maths.pi / 2.0 - maths.atan(y / x) if x < 0.0: radians += maths.pi else: radians = maths.atan(x / y) if y < 0.0: radians += maths.pi if radians < 0.0: radians += 2.0 * maths.pi return degrees_from_radians(radians)
[docs]def azimuths(va): # 'azimuth' is synonymous with 'compass bearing' """Returns the compass bearings in degrees of the direction of each vector in va (x = East, y = North), ignoring z.""" assert va.ndim > 1 and 2 <= va.shape[-1] <= 3 shape = tuple(list(va.shape[:-1]) + [3]) z_zero_v = np.zeros(shape) z_zero_v[..., :2] = va[..., :2] unit_v = unit_vectors(z_zero_v) # also checks that z_zero_v is not zero vector x = unit_v[..., 0] y = unit_v[..., 1] # ignore z component # todo: handle cases where x == y == 0 restore = np.seterr(all = 'ignore') radians = np.where( np.abs(x) >= np.abs(y), np.where(x < 0.0, maths.pi * 3.0 / 2.0 - np.arctan(y / x), maths.pi / 2.0 - np.arctan(y / x)), np.where(y < 0.0, maths.pi + np.arctan(x / y), np.arctan(x / y))) radians = radians % (2.0 * maths.pi) np.seterr(**restore) return np.degrees(radians)
[docs]def inclination(v): """Returns the inclination in degrees of v (angle relative to +ve z axis).""" assert len(v) == 3 unit_v = unit_vector(v) radians = maths.acos(unit_v[2]) return degrees_from_radians(radians)
[docs]def inclinations(a): """Returns the inclination in degrees of each vector in a (angle relative to +ve z axis).""" assert a.ndim > 1 and a.shape[-1] == 3 unit_vs = unit_vectors(a) radians = np.arccos(unit_vs[..., 2]) return degrees_from_radians(radians)
[docs]def nan_inclinations(a, already_unit_vectors = False): """Returns the inclination in degrees of each vector in a (angle relative to +ve z axis), allowing NaNs.""" assert a.ndim > 1 and a.shape[-1] == 3 unit_vs = a if already_unit_vectors else nan_unit_vectors(a) restore = np.seterr(all = 'ignore') radians = np.where(np.isnan(unit_vs[..., 2]), np.nan, np.arccos(unit_vs[..., 2])) result = np.where(np.isnan(radians), np.nan, degrees_from_radians(radians)) np.seterr(**restore) return result
[docs]def points_direction_vector(a, axis): """Returns an average direction vector based on first and last non-NaN points or slices in given axis.""" # note: as currently coded, might give poor results with some patterns of NaNs assert a.ndim > 1 and 0 <= axis < a.ndim - 1 and a.shape[-1] > 1 and a.shape[axis] > 1 if np.all(np.isnan(a)): return None start = 0 start_slicing = [slice(None)] * a.ndim while True: start_slicing[axis] = slice(start, start + 1) if not np.all(np.isnan(a[tuple(start_slicing)])): break start += 1 assert start < a.shape[axis] finish = a.shape[axis] - 1 finish_slicing = [slice(None)] * a.ndim while True: finish_slicing[axis] = slice(finish, finish + 1) if not np.all(np.isnan(a[tuple(finish_slicing)])): break finish -= 1 if start >= finish: return None if a.ndim > 2: mean_axes = tuple(range(a.ndim - 1)) start_p = np.nanmean(a[tuple(start_slicing)], axis = mean_axes) finish_p = np.nanmean(a[tuple(finish_slicing)], axis = mean_axes) else: start_p = a[start] finish_p = a[finish] return finish_p - start_p
[docs]def dot_product(a, b): # pragma: no cover """Returns the dot product (scalar product) of the two vectors.""" return np.dot(a, b)
[docs]def dot_products(a, b): # pragma: no cover """Returns the dot products of pairs of vectors; last axis covers element of a vector.""" return np.sum(a * b, axis = -1)
[docs]def cross_product(a, b): # pragma: no cover """Returns the cross product (vector product) of the two vectors.""" return np.cross(a, b)
[docs]def naive_length(v): # pragma: no cover """Returns the length of the vector assuming consistent units.""" return np.linalg.norm(v)
[docs]def naive_lengths(v): """Returns the lengths of the vectors assuming consistent units.""" return np.sqrt(np.sum(v * v, axis = -1))
[docs]def naive_2d_length(v): # pragma: no cover """Returns the length of the vector projected onto xy plane, assuming consistent units.""" return np.linalg.norm(v[0:2])
[docs]def naive_2d_lengths(v): """Returns the lengths of the vectors projected onto xy plane, assuming consistent units.""" v2d = v[..., :2] return np.sqrt(np.sum(v2d * v2d, axis = -1))
[docs]def unit_corrected_length(v, unit_conversion): """Returns the length of the vector v after applying the unit_conversion factors. arguments: v (1D numpy float array): vector with mixed units of measure unit_conversion (1D numpy float array): vector to multiply elements of v by, prior to finding length returns: float, being the length of v after adjustment by unit_conversion notes: example unit_conversion might be: [1.0, 1.0, 0.3048] to convert z from feet to metres, or [3.28084, 3.28084, 1.0] to convert x and y from metres to feet """ converted = elemental_multiply(v, unit_conversion) return naive_length(converted)
[docs]def manhatten_distance(p1, p2): """Returns the Manhattan distance between two points.""" return abs(p2[0] - p1[0]) + abs(p2[1] - p1[1]) + abs(p2[2] - p1[2])
[docs]def manhattan_distance(p1, p2): #  alternative spelling to above """Returns the Manhattan distance between two points.""" return abs(p2[0] - p1[0]) + abs(p2[1] - p1[1]) + abs(p2[2] - p1[2])
[docs]def radians_difference(a, b): """Returns the angle between two vectors, in radians.""" return maths.acos(min(1.0, max(-1.0, dot_product(unit_vector(a), unit_vector(b)))))
[docs]def degrees_difference(a, b): """Returns the angle between two vectors, in degrees.""" return degrees_from_radians(radians_difference(a, b))
[docs]def rotation_matrix_3d_axial(axis, angle): """Returns a rotation matrix which will rotate points about axis (0: x, 1: y, or 2: z) by angle in degrees. note: this function follows the mathematical convention: a positive angle results in anti-clockwise rotation when viewed in direction of positive axis """ axis_a = (axis + 1) % 3 axis_b = (axis_a + 1) % 3 matrix = np.eye(3) radians = np.radians(angle) cosine = np.cos(radians) sine = np.sin(radians) matrix[axis_a, axis_a] = cosine matrix[axis_b, axis_b] = cosine matrix[axis_a, axis_b] = -sine # left handed coordinate system, eg. UTM & depth matrix[axis_b, axis_a] = sine return matrix
[docs]def no_rotation_matrix(): # pragma: no cover """Returns a rotation matrix which will not move points (identity matrix).""" return np.eye(3)
[docs]def rotation_3d_matrix(xzy_axis_angles): """Returns a rotation matrix which will rotate points about the x, z, then y axis by angles in degrees.""" angles = np.radians(xzy_axis_angles) cos_c, cos_a, cos_b = np.cos(angles) sin_c, sin_a, sin_b = np.sin(angles) sin_a_sin_c = sin_a * sin_c sin_a_cos_c = sin_a * cos_c rotation_matrix = np.array([ [cos_a * cos_b, -sin_a_cos_c * cos_b + sin_b * sin_c, cos_b * sin_a_sin_c + sin_b * cos_c], [sin_a, cos_a * cos_c, -sin_c * cos_a], [-sin_b * cos_a, sin_a_cos_c * sin_b + cos_b * sin_c, -sin_b * sin_a_sin_c + cos_b * cos_c], ]) return rotation_matrix
[docs]@njit def rotation_3d_matrix_njit(xzy_axis_angles): # pragma: no cover """Returns a rotation matrix which will rotate points about the x, z, then y axis by angles in degrees.""" angles = np.radians(xzy_axis_angles) cos_c, cos_a, cos_b = np.cos(angles) sin_c, sin_a, sin_b = np.sin(angles) sin_a_sin_c = sin_a * sin_c sin_a_cos_c = sin_a * cos_c rotation_matrix = np.array([ [cos_a * cos_b, -sin_a_cos_c * cos_b + sin_b * sin_c, cos_b * sin_a_sin_c + sin_b * cos_c], [sin_a, cos_a * cos_c, -sin_c * cos_a], [-sin_b * cos_a, sin_a_cos_c * sin_b + cos_b * sin_c, -sin_b * sin_a_sin_c + cos_b * cos_c], ]) return rotation_matrix
[docs]def reverse_rotation_3d_matrix(xzy_axis_angles): # pragma: no cover """Returns a rotation matrix which will rotate points about the y, z, then x axis by angles in degrees.""" return rotation_3d_matrix(xzy_axis_angles).T
[docs]def rotate_vector(rotation_matrix, vector): # pragma: no cover """Returns the rotated vector.""" return np.dot(rotation_matrix, vector)
[docs]def rotate_array(rotation_matrix, a): """Returns a copy of array a with each vector rotated by the rotation matrix.""" return np.matmul(rotation_matrix, a.reshape(-1, 3).T).T.reshape(a.shape)
[docs]@njit def rotate_array_njit(rotation_matrix, a): # pragma: no cover """Returns a copy of array a with each vector rotated by the rotation matrix.""" return np.dot(rotation_matrix, a.reshape(-1, 3).T).T.reshape(a.shape)
[docs]def rotate_xyz_array_around_z_axis(a, target_xy_vector): """Returns a copy of array a suitable for presenting a cross-section using the resulting x,z values. arguments: a (numpy float array of shape (..., 3)): the xyz points to be rotated target_xy_vector (2 (or 3) floats): a vector indicating which direction in source xy space will end up being mapped to the positive x axis in the returned data returns: numpy float array of same shape as a notes: if the input points of a lie in a vertical plane parallel to the target xy vector, then the resulting points will have constant y values; in general, a full rotation of the points is applied, so resulting y values will indicate distance 'into the page' for non-planar or unaligned data """ target_v = np.zeros(3) target_v[:2] = target_xy_vector[:2] rotation_angle = azimuth(target_v) - 90.0 rotation_matrix = rotation_matrix_3d_axial(2, rotation_angle) # todo: check sign of rotation angle return rotate_array(rotation_matrix, a)
[docs]def unit_vector_from_azimuth_and_inclination(azimuth, inclination): """Returns unit vector with compass bearing of azimuth and inclination off +z axis. note: assumes a left handed coordinate system with y axis north and x axis east """ matrix = rotation_3d_matrix((inclination, 180.0 - azimuth, 0.0)) return rotate_vector(matrix, np.array((0.0, 0.0, 1.0)))
[docs]def tilt_3d_matrix(azimuth, dip): """Returns a 3D rotation matrix for applying a dip in a certain azimuth. note: if azimuth is compass bearing in degrees, and dip is in degrees, the resulting matrix can be used to rotate xyz points where x values are eastings, y values are northings and z increases downwards """ matrix = rotation_matrix_3d_axial(2, -azimuth) # will yield rotation around z axis so azimuth goes north matrix = np.dot(matrix, rotation_matrix_3d_axial(0, dip)) # adjust for dip matrix = np.dot(matrix, rotation_matrix_3d_axial(2, azimuth)) # rotate azimuth back to original return matrix
[docs]def rotation_matrix_3d_vector(v): """Returns a rotation matrix which will rotate vector v to the vertical (z) axis. note: the returned matrix will map a positive z axis vector onto v """ v = v / np.linalg.norm(v) kmat = np.array([[0, 0, -v[0]], [0, 0, -v[1]], [v[0], v[1], 0]]) rotation_matrix = np.eye(3) + kmat + kmat.dot(kmat) * (1 / (1 + v[2])) rotation_matrix[:2] = -rotation_matrix[:2] return rotation_matrix
[docs]@njit def rotation_matrix_3d_vector_njit(v): # pragma: no cover """Returns a rotation matrix which will rotate points by inclination and azimuth of vector. note: the returned matrix will map a positive z axis vector onto v """ v = v / np.linalg.norm(v) kmat = np.array([[0, 0, -v[0]], [0, 0, -v[1]], [v[0], v[1], 0]]) rotation_matrix = np.eye(3) + kmat + kmat.dot(kmat) * (1 / (1 + v[2])) rotation_matrix[:2] = -rotation_matrix[:2] return rotation_matrix
[docs]def tilt_points(pivot_xyz, azimuth, dip, points): """Modifies array of xyz points in situ to apply dip in direction of azimuth, about pivot point.""" matrix = tilt_3d_matrix(azimuth, dip) points_shape = points.shape points[:] -= pivot_xyz points[:] = np.matmul(matrix, points.reshape((-1, 3)).transpose()).transpose().reshape(points_shape) points[:] += pivot_xyz
[docs]def project_points_onto_plane(plane_xyz, normal_vector, points): """Modifies array of xyz points in situ to project onto a plane defined by a point and normal vector. note: implicit xy & z units must be the same """ az = azimuth(normal_vector) incl = inclination(normal_vector) tilt_points(plane_xyz, az, -incl, points) points[..., 2] = plane_xyz[2] tilt_points(plane_xyz, az, incl, points)
[docs]def perspective_vector(xyz_box, view_axis, vanishing_distance, vector): """Returns a version of vector with a perspective applied.""" mid_points = np.zeros(3) xyz_ranges = np.zeros(3) result = np.zeros(3) for axis in range(3): mid_points[axis] = 0.5 * (xyz_box[0, axis] + xyz_box[1, axis]) xyz_ranges[axis] = xyz_box[1, axis] - xyz_box[0, axis] factor = 1.0 - (vector[view_axis] - xyz_box[0, view_axis]) / (vanishing_distance * (xyz_ranges[view_axis])) result[view_axis] = vector[view_axis] for axis in range(3): if axis == view_axis: continue result[axis] = mid_points[axis] + factor * (vector[axis] - mid_points[axis]) return result
[docs]def determinant(a, b, c): """Returns the determinant of the 3 x 3 matrix comprised of the 3 vectors.""" return (a[0] * b[1] * c[2] + a[1] * b[2] * c[0] + a[2] * b[0] * c[1] - a[2] * b[1] * c[0] - a[1] * b[0] * c[2] - a[0] * b[2] * c[1])
[docs]def determinant_3x3(a): """Returns the determinant of the 3 x 3 matrix.""" return determinant(a[0], a[1], a[2])
[docs]def clockwise(a, b, c): """Returns a +ve value if 2D points a,b,c are in clockwise order, 0.0 if in line, -ve for ccw. note: assumes positive y-axis is anticlockwise from positive x-axis """ return (c[0] - a[0]) * (b[1] - a[1]) - ((c[1] - a[1]) * (b[0] - a[0]))
[docs]def clockwise_triangles(p, t, projection = 'xy'): """Returns a numpy array of +ve values where triangle points are in clockwise order, 0.0 if in line, -ve for ccw. arguments: p (numpy float array of shape (N, 2 or 3)): points in use as vertices of triangles t (numpy int array of shape (M, 3)): indices into first axis of p defining the triangles projection (string, default 'xy'): one of 'xy', 'xz' or 'yz' being the direction of projection, ie. which elements of the second axis of p to use; must be 'xy' if p has shape (N, 2) returns: numpy float array of shape (M,) being the +ve or -ve values indicating clockwise or anti-clockwise ordering of each triangle's vertices when projected onto the specified plane and viewed in the direction negative to positive of the omitted axis note: assumes xyz axes are left-handed (reverse the result for a right handed system) """ a0, a1 = _projected_xyz_axes(projection) return (((p[t[:, 2], a0] - p[t[:, 0], a0]) * (p[t[:, 1], a1] - p[t[:, 0], a1])) - ((p[t[:, 2], a1] - p[t[:, 0], a1]) * (p[t[:, 1], a0] - p[t[:, 0], a0])))
[docs]def in_triangle(a, b, c, d): """Returns True if point d lies wholly within the triangle pf ccw points a, b, c, projected onto xy plane. note: a, b & c must be sorted into anti-clockwise order before calling this function """ return clockwise(a, b, d) < 0.0 and clockwise(b, c, d) < 0.0 and clockwise(c, a, d) < 0.0
[docs]def in_triangle_edged(a, b, c, d): """Returns True if d lies within or on the boudnary of triangle of ccw points a,b,c projected onto xy plane. note: a, b & c must be sorted into anti-clockwise order before calling this function """ return clockwise(a, b, d) <= 0.0 and clockwise(b, c, d) <= 0.0 and clockwise(c, a, d) <= 0.0
[docs]def points_in_triangles(p, t, da, projection = 'xy', edged = False): """Returns 2D numpy bool array indicating which of points da are within which triangles. arguments: p (numpy float array of shape (N, 2 or 3)): points in use as vertices of triangles t (numpy int array of shape (M, 3)): indices into first axis of p defining the triangles da (numpy float array of shape (D, 2 or 3)): points to test for projection (string, default 'xy'): one of 'xy', 'xz' or 'yz' being the direction of projection, ie. which elements of the second axis of p and da to use; must be 'xy' if p and da have shape (N, 2) edged (bool, default False): if True, points lying exactly on the edge of a triangle are included as being in the triangle, otherwise they are excluded returns: numpy bool array of shape (M, D) indicating which points are within which triangles note: the triangles do not need to be in a consistent clockwise or anti-clockwise order """ assert p.ndim == 2 and t.ndim == 2 and da.ndim == 2 and da.shape[1] == p.shape[1] cwt = clockwise_triangles(p, t, projection = projection) a0, a1 = _projected_xyz_axes(projection) d_count = len(da) d_base = len(p) t_count = len(t) pp = np.concatenate((p, da), axis = 0) # build little triangles using da points and two of triangle vertices tp = np.empty((t_count, d_count, 3, 3), dtype = int) tp[:, :, :, 2] = np.arange(d_count).reshape((1, -1, 1)) + d_base tp[:, :, 0, 0] = np.where(cwt > 0.0, t[:, 1], t[:, 0]).reshape((-1, 1)) tp[:, :, 0, 1] = np.where(cwt > 0.0, t[:, 0], t[:, 1]).reshape((-1, 1)) tp[:, :, 1, 0] = np.where(cwt > 0.0, t[:, 2], t[:, 1]).reshape((-1, 1)) tp[:, :, 1, 1] = np.where(cwt > 0.0, t[:, 1], t[:, 2]).reshape((-1, 1)) tp[:, :, 2, 0] = np.where(cwt > 0.0, t[:, 0], t[:, 2]).reshape((-1, 1)) tp[:, :, 2, 1] = np.where(cwt > 0.0, t[:, 2], t[:, 0]).reshape((-1, 1)) cwtd = clockwise_triangles(pp, tp.reshape((-1, 3)), projection = projection).reshape((t_count, d_count, 3)) if edged: return np.all(cwtd <= 0.0, axis = 2) else: return np.all(cwtd < 0.0, axis = 2)
[docs]@njit def point_in_polygon(x, y, polygon): # pragma: no cover """Calculates if a point in within a polygon in 2D. arguments: x (float): the point's x-coordinate. y (float): the point's y-coordinate. polygon (np.ndarray): array of the polygon's vertices in 2D. returns: inside (bool): True if point is within the polygon, False otherwise. note: the polygon is assumed closed, the closing point should not be repeated """ n = len(polygon) inside = False p2x = 0.0 p2y = 0.0 xints = 0.0 p1x, p1y = polygon[0] for i in range(n + 1): p2x, p2y = polygon[i % n] if y > min(p1y, p2y): if y <= max(p1y, p2y): if x <= max(p1x, p2x): if p1y != p2y: xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x if p1x == p2x or x <= xints: inside = not inside p1x, p1y = p2x, p2y return inside
[docs]@njit def point_in_triangle(x, y, triangle): # pragma: no cover """Calculates if a point in within a triangle in 2D. arguments: x (float): the point's x-coordinate. y (float): the point's y-coordinate. triangle (np.ndarray): array of the triangles's vertices in 2D, of shape (3, 2) returns: inside (bool): True if point is within the polygon, False otherwise. """ p0x = triangle[0, 0] p1x = triangle[1, 0] p2x = triangle[2, 0] min_p01x = min(p0x, p1x) if x < min(min_p01x, p2x): return False max_p01x = max(p0x, p1x) if x > max(max_p01x, p2x): return False p0y = triangle[0, 1] p1y = triangle[1, 1] p2y = triangle[2, 1] min_p01y = min(p0y, p1y) if y < min(min_p01y, p2y): return False max_p01y = max(p0y, p1y) if y > max(max_p01y, p2y): return False inside = False xints = 0.0 if y > min_p01y: if y <= max_p01y: if x <= max_p01x: if p0x == p1x: inside = True else: if p0y == p1y: inside = True else: xints = (y - p0y) * (p1x - p0x) / (p1y - p0y) + p0x if x <= xints: inside = True if y > min(p1y, p2y): if y <= max(p1y, p2y): if x <= max(p1x, p2x): if p1x == p2x: inside = not inside else: if p1y == p2y: inside = not inside else: xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x if x <= xints: inside = not inside if y > min(p2y, p0y): if y <= max(p2y, p0y): if x <= max(p2x, p0x): if p2x == p0x: inside = not inside else: if p2y == p0y: inside = not inside else: xints = (y - p2y) * (p0x - p2x) / (p0y - p2y) + p2x if x <= xints: inside = not inside return inside
[docs]@njit(parallel = True) def points_in_polygon(points: np.ndarray, polygon: np.ndarray, points_xlen: int, polygon_num: int = 0) -> np.ndarray: # pragma: no cover """Calculates which points are within a polygon in 2D. arguments: points (np.ndarray): array of shape (N, 2 or 3), of the points in 2D (xy, any z values are ignored) polygon (np.ndarray): list-like array of the polygon's vertices in 2D points_xlen (int): the original I extent of the now flattened points, use 1 if not applicable polygon_num (int): the polygon number, default is 0, for copying to output returns: polygon_points (np.ndarray): list-like 2D array containing only the points within the polygon, with each row being the polygon number (as input), points J index, and points I index note: the polygon is assumed closed, the closing point should not be repeated """ polygon_points = np.full((len(points), 3), -1, dtype = np.int32) x = points[:, 0] y = points[:, 1] p = np.empty(len(points), dtype = np.bool_) j = np.empty(len(points), dtype = np.int32) i = np.empty(len(points), dtype = np.int32) for point_num in numba.prange(len(points)): p[point_num] = point_in_polygon(x[point_num], y[point_num], polygon) if p[point_num] is True: j[point_num], i[point_num] = divmod(point_num, points_xlen) polygon_points[point_num, 0] = polygon_num polygon_points[point_num, 1] = j[point_num] polygon_points[point_num, 2] = i[point_num] return polygon_points[polygon_points[:, 0] != -1]
[docs]@njit def points_in_triangle(points: np.ndarray, triangle: np.ndarray, points_xlen: int, triangle_num: int = 0) -> np.ndarray: # pragma: no cover """Calculates which points are within a triangle in 2D. arguments: points (np.ndarray): array of the points in 2D. triangle (np.ndarray): array of the triangle's vertices in 2D, shape (3, 2). points_xlen (int): the number of unique x coordinates. triangle_num (int): the triangle number, default is 0. returns: triangle_points (np.ndarray): 2D array containing only the points within the triangle, with each row being the triangle number, points y index, and points x index. """ triangle_points = np.full((len(points), 3), -1, dtype = np.int32) for point_num in numba.prange(len(points)): p = point_in_triangle(points[point_num, 0], points[point_num, 1], triangle) if p is True: yi, xi = divmod(point_num, points_xlen) triangle_points[point_num] = [triangle_num, yi, xi] return triangle_points[triangle_points[:, 0] != -1]
[docs]@njit def mesh_points_in_triangle(triangle: np.ndarray, points_xlen: int, points_ylen: int, triangle_num: int = 0) -> np.ndarray: # pragma: no cover """Calculates which implicit mesh points are within a triangle in 2D for normalised triangle. arguments: triangle (np.ndarray): array of the triangle's vertices in 2D, shape (3, 2). points_xlen (int): the number of unique x coordinates, starting at 0.0, spacing 1.0. points_ylen (int): the number of unique y coordinates, starting at 0.0, spacing 1.0. triangle_num (int): the triangle number, default is 0. returns: triangle_points (np.ndarray): 2D array containing only the points within the triangle, with each row being the triangle number, points y index, and points x index. """ triangle_points = np.empty((0, 3), dtype = numba.int32) y = 0.0 for yi in numba.prange(points_ylen): x = 0.0 for xi in numba.prange(points_xlen): p = point_in_triangle(x, y, triangle) if p is True: triangle_points = np.append(triangle_points, np.array([[triangle_num, yi, xi]], dtype = numba.int32), axis = 0) x += 1.0 y += 1.0 return triangle_points
[docs]@njit def points_in_polygons(points: np.ndarray, polygons: np.ndarray, points_xlen: int) -> np.ndarray: # pragma: no cover """Calculates which points are within which polygons in 2D. arguments: points (np.ndarray): array of the points in 2D. polygons (np.ndarray): array of each polygons' vertices in 2D. points_xlen (int): the number of unique x coordinates. returns: polygons_points (np.ndarray): 2D array (list-like) containing only the points within each polygon, with each row being the polygon number, points y index, and points x index. """ polygons_points = np.empty((0, 3), dtype = numba.int32) for polygon_num in numba.prange(len(polygons)): polygon_points = points_in_polygon(points, polygons[polygon_num], points_xlen, polygon_num) polygons_points = np.append(polygons_points, polygon_points, axis = 0) return polygons_points
[docs]@njit(parallel = False) def points_in_triangles_njit(points: np.ndarray, triangles: np.ndarray, points_xlen: int) -> np.ndarray: # pragma: no cover """Calculates which points are within which triangles in 2D. arguments: points (np.ndarray): array of the points in 2D. triangles (np.ndarray): array of each triangles' vertices in 2D, shape (N, 3, 2). points_xlen (int): the number of unique x coordinates. returns: triangles_points (np.ndarray): 2D array (list-like) containing only the points within each triangle, with each row being the triangle number, points y index, and points x index. """ triangles_points = np.empty((0, 3), dtype = numba.int32) for triangle_num in range(len(triangles)): triangle_points = points_in_triangle(points, triangles[triangle_num], points_xlen, triangle_num) triangles_points = np.append(triangles_points, triangle_points, axis = 0) return triangles_points
[docs]@njit def meshgrid(x: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: # pragma: no cover """Returns coordinate matrices from coordinate vectors x and y. arguments: x (np.ndarray): 1d array of x coordinates. y (np.ndarray): 1d array of y coordinates. returns: Tuple containing: - xx (np.ndarray): the elements of x repeated to fill the matrix along the first dimension. - yy (np.ndarray): the elements of y repeated to fill the matrix along the second dimension. """ xx = np.empty(shape = (y.size, x.size), dtype = x.dtype) yy = np.empty(shape = (y.size, x.size), dtype = y.dtype) for i in range(y.size): for j in range(x.size): xx[i, j] = x[j] yy[i, j] = y[i] return xx, yy
[docs]def points_in_triangles_aligned(nx: int, ny: int, dx: float, dy: float, triangles: np.ndarray) -> np.ndarray: """Calculates which points are within which triangles in 2D for a regular mesh of aligned points. arguments: nx (int): number of points in x axis ny (int): number of points in y axis dx (float): spacing of points in x axis (first point is at half dx) dy (float): spacing of points in y axis (first point is at half dy) triangles (np.ndarray): float array of each triangles' vertices in 2D, shape (N, 3, 2). points_xlen (int): the number of unique x coordinates. returns: triangles_points (np.ndarray): 2D array (list-like) containing only the points within each triangle, with each row being the triangle number, points y index, and points x index. """ triangles_points = np.empty((0, 3), dtype = np.int32) dx_dy = np.expand_dims(np.array([dx, dy], dtype = np.float32), axis = 0) # for triangle_num in numba.prange(len(triangles)): for triangle_num in range(len(triangles)): tp = (triangles[triangle_num] / dx_dy) - 0.5 min_tpx = max(maths.ceil(min(tp[0, 0], tp[1, 0], tp[2, 0])), 0) max_tpx = min(maths.floor(max(tp[0, 0], tp[1, 0], tp[2, 0])), nx - 1) if max_tpx < min_tpx: continue min_tpy = max(maths.ceil(min(tp[0, 1], tp[1, 1], tp[2, 1])), 0) max_tpy = min(maths.floor(max(tp[0, 1], tp[1, 1], tp[2, 1])), ny - 1) if max_tpy < min_tpy: continue ntpx = max_tpx - min_tpx + 1 ntpy = max_tpy - min_tpy + 1 # x = np.linspace(float(min_tpx), float(max_tpx), ntpx) # y = np.linspace(float(min_tpy), float(max_tpy), ntpy) # p = np.stack(meshgrid(x, y), axis = -1).reshape((-1, 2)) # triangle_points = points_in_triangle(p, tp, ntpx, triangle_num) tp[:, 0] -= float(min_tpx) tp[:, 1] -= float(min_tpy) triangle_points = mesh_points_in_triangle(tp, ntpx, ntpy, triangle_num) triangle_points[:, 1] += min_tpy triangle_points[:, 2] += min_tpx triangles_points = np.append(triangles_points, triangle_points, axis = 0) return triangles_points
[docs]@njit def triangle_box(triangle: np.ndarray) -> Tuple[float, float, float, float]: # pragma: no cover """Finds the minimum and maximum x and y values of a single traingle. arguments: triangle (np.ndarray): array of the traingle's vertices' x and y coordinates. returns: Tuple containing: - (float): minimum x value. - (float): maximum x value. - (float): minimum y value. - (float): maximum y value. """ x_values = triangle[:, 0] y_values = triangle[:, 1] return min(x_values), max(x_values), min(y_values), max(y_values)
[docs]@njit def vertical_intercept(x: float, x_values: np.ndarray, y_values: np.ndarray) -> Optional[float]: # pragma: no cover """Finds the y value of a straight line between two points at a given x. If the x value given is not within the x values of the points, returns None. arguments: x (float): x value at which to determine the y value. x_values (np.ndarray): the x coordinates of point 1 and point 2. y_values (np.ndarray): the y coordinates of point 1 and point 2. returns: y (Optional[float]): y value of the straight line between point 1 and point 2, evaluated at x. If x is outside the x_values range, y is None. """ y = None if x >= min(x_values) and x <= max(x_values): if x_values[0] == x_values[1]: y = y_values[0] else: m = (y_values[1] - y_values[0]) / (x_values[1] - x_values[0]) c = y_values[1] - m * x_values[1] y = m * x + c return y
[docs]@njit def points_in_triangles_aligned_optimised(nx: int, ny: int, dx: float, dy: float, triangles: np.ndarray) -> np.ndarray: # pragma: no cover """Calculates which points are within which triangles in 2D for a regular mesh of aligned points. arguments: nx (int): number of points in x axis ny (int): number of points in y axis dx (float): spacing of points in x axis (first point is at half dx) dy (float): spacing of points in y axis (first point is at half dy) triangles (np.ndarray): float array of each triangles' vertices in 2D, shape (N, 3, 2) returns: triangles_points (np.ndarray): 2D array (list-like) containing only the points within each triangle, with each row being the triangle number, points y index, and points x index """ grid_x = dx * (0.5 + np.arange(nx).astype(np.float64)) grid_y = dy * (0.5 + np.arange(ny).astype(np.float64)) triangles_points_list = [] for triangle_num in range(len(triangles)): triangle = triangles[triangle_num] min_x, max_x, min_y, max_y = triangle_box(triangle) x_values = grid_x[np.logical_and(grid_x >= min_x, grid_x < max_x)] y_values = grid_y[np.logical_and(grid_y >= min_y, grid_y < max_y)] for x in x_values: ys = [] ys.append(vertical_intercept(x, triangle[1:, 0], triangle[1:, 1])) ys.append(vertical_intercept(x, triangle[:2, 0], triangle[:2, 1])) ys.append(vertical_intercept(x, triangle[::2, 0], triangle[::2, 1])) ys = [u for u in ys if u is not None] valid_y = y_values[np.logical_and(y_values >= min(ys), y_values <= max(ys))] x_idx = int(x / dx) triangles_points_list.extend([[triangle_num, int(y / dy), x_idx] for y in valid_y]) if len(triangles_points_list) == 0: triangles_points = np.empty((0, 3), dtype = np.int32) return triangles_points triangles_points = np.array(triangles_points_list, dtype = np.int32) return triangles_points
[docs]def triangle_normal_vector(p3): """For a triangle in 3D space, defined by 3 vertex points, returns a unit vector normal to the plane of the triangle. note: resulting vector implicitly assumes that xy & z units are the same; if this is not the case, adjust vector afterwards as required """ # todo: handle degenerate triangles return unit_vector(cross_product(p3[0] - p3[1], p3[0] - p3[2]))
[docs]@njit def triangle_normal_vector_numba(points): # pragma: no cover """For a triangle in 3D space, defined by 3 vertex points, returns a unit vector normal to the plane of the triangle. note: resulting vector implicitly assumes that xy & z units are the same; if this is not the case, adjust vector afterwards as required """ v = np.cross(points[0] - points[1], points[0] - points[2]) return v / np.linalg.norm(v)
[docs]def in_circumcircle(a, b, c, d): """Returns True if point d lies within the circumcircle pf ccw points a, b, c, projected onto xy plane. note: a, b & c must be sorted into anti-clockwise order before calling this function """ m = np.empty((3, 3)) m[0, :2] = a[:2] - d[:2] m[1, :2] = b[:2] - d[:2] m[2, :2] = c[:2] - d[:2] m[:, 2] = (m[:, 0] * m[:, 0]) + (m[:, 1] * m[:, 1]) return determinant_3x3(m) > 0.0
[docs]def point_distance_to_line_2d(p, l1, l2): """Ignoring any z values, returns the xy distance of point p from line passing through l1 and l2.""" if np.all(l2[:2] == l1[:2]): return naive_2d_length(p[:2] - l1[:2]) return (abs(p[0] * (l1[1] - l2[1]) + l1[0] * (l2[1] - p[1]) + l2[0] * (p[1] - l1[1])) / naive_2d_length(l2[:2] - l1[:2]))
[docs]def point_distance_to_line_segment_2d(p, l1, l2): """Ignoring any z values, returns the xy distance of point p from line segment between l1 and l2.""" if is_obtuse_2d(l1, p, l2): return naive_2d_length(p[:2] - l1[:2]) elif is_obtuse_2d(l2, p, l1): return naive_2d_length(p[:2] - l2[:2]) else: return point_distance_to_line_2d(p, l1, l2)
[docs]def is_obtuse_2d(p, p1, p2): """Returns True if the angle at point p subtended by points p1 and p2, in xy plane, is greater than 90 degrees; else False.""" return np.dot((p1[:2] - p[:2]), (p2[:2] - p[:2])) < 0.0
[docs]def isclose(a, b, tolerance = 1.0e-6): """Returns True if the two points are extremely close to one another (ie. the same point). """ # return np.all(np.isclose(a, b, atol = tolerance)) # cheap and cheerful alternative to thorough numpy version commented out above return np.max(np.abs(a - b)) <= tolerance
[docs]def is_close(a, b, tolerance = 1.0e-6): """Returns True if the two points are extremely close to one another (ie. the same point). """ return isclose(a, b, tolerance = tolerance)
[docs]def point_distance_sqr_to_points_projected(p, points, projection): """Returns an array of projected distances squared between p and points; projection is 'xy', 'xz' or 'yz'.""" if projection == 'xy': d = points[..., :2] - p[:2] elif projection == 'xz': d = points[..., 0:3:2] - p[0:3:2] elif projection == 'yz': d = points[..., 1:] - p[1:] else: raise ValueError("projection must be 'xy', 'xz' or 'yz'") return np.sum(d * d, axis = -1)
[docs]def nearest_point_projected(p, points, projection): """Returns the index into points array closest to point p; projection is 'xy', 'xz' or 'yz'.""" # note: in the case of equidistant points, the index of the 'first' point is returned d2 = point_distance_sqr_to_points_projected(p, points, projection) return np.unravel_index(np.nanargmin(d2), d2.shape)
[docs]def area_of_triangle(a, b, c): """Returns the area of the triangle defined by three vertices.""" # uses Heron's formula la = naive_length(a - b) lb = naive_length(b - c) lc = naive_length(c - a) s = 0.5 * (la + lb + lc) return maths.sqrt(s * (s - la) * (s - lb) * (s - lc))
[docs]def area_of_triangles(p, t, xy_projection = False): """Returns numpy array of areas of triangles, optionally when projected onto xy plane.""" # uses Heron's formula pt = p[t] if xy_projection: la = naive_2d_lengths(pt[:, 0, :] - pt[:, 1, :]) lb = naive_2d_lengths(pt[:, 1, :] - pt[:, 2, :]) lc = naive_2d_lengths(pt[:, 2, :] - pt[:, 0, :]) else: la = naive_lengths(pt[:, 0, :] - pt[:, 1, :]) lb = naive_lengths(pt[:, 1, :] - pt[:, 2, :]) lc = naive_lengths(pt[:, 2, :] - pt[:, 0, :]) s = 0.5 * (la + lb + lc) return np.sqrt(s * (s - la) * (s - lb) * (s - lc))
[docs]def clockwise_sorted_indices(p, b): """Returns a clockwise sorted numpy list of indices b into the points p. note: this function is designed for preparing a list of points defining a convex polygon when projected in the xy plane, starting from a subset of the unsorted points; more specifically, it assumes that the mean of p (over axis 0) lies within the polygon and the clockwise ordering is relative to that mean point """ # note: this function currently assumes that the mean of points bp lies within the hull of bp # and that the points form a convex polygon from the perspective of the mean point assert p.ndim == 2 and len(p) >= 3 centre = np.mean(p, axis = 0) hull_list = [] # list of azimuths and indices into p (axis 0) for i in b: azi = azimuth(p[i] - centre) hull_list.append((azi, i)) return np.array([i for (_, i) in sorted(hull_list)], dtype = int)
[docs]def xy_sorted(p, axis = None): """Returns copy of points p sorted according to x or y (whichever has greater range). arguments: p (numpy float array of shape (..., 2) or (..., 3)): points to be sorted axis (int, optional): 0 for x sort; 1 for y sort; None for whichever has greater range returns: p', axis where p' is a list-like (2D) version of p, sorted by either x or y and axis is 0 if the sort was by x, 1 if it were by y note: returned array is always 2D, ie. list of points """ assert p.ndim >= 2 and p.shape[-1] >= 2 p = p.reshape((-1, p.shape[-1])) if axis is None: xy_range = np.nanmax(p, axis = 0) - np.nanmin(p, axis = 0) axis = (1 if xy_range[1] > xy_range[0] else 0) spi = np.argsort(p[:, axis]) return p[spi], axis
[docs]@njit def xy_sorted_njit(p, axis = -1): # pragma: no cover """Returns copy of points p sorted according to x or y (whichever has greater range).""" assert p.ndim >= 2 and p.shape[-1] >= 2 p = p.reshape((-1, p.shape[-1])) if axis == -1: xy_range = _nanmax(p) - _nanmin(p) axis = (1 if xy_range[1] > xy_range[0] else 0) points = p[:, axis] spi = np.argsort(points) return p[spi], axis
@njit def _nanmax(array): # pragma: no cover """Numba implementation of np.nanmax with axis = 0.""" len0 = array.shape[1] max_array = np.empty(len0) for col in range(len0): col_array = array[:, col] max_val = col_array[0] for val in col_array: if val > max_val and not np.isnan(val): max_val = val max_array[col] = max_val return max_array @njit def _nanmin(array): # pragma: no cover """Numba implementation of np.nanmin with axis = 0.""" len0 = array.shape[1] min_array = np.empty(len0) for col in range(len0): col_array = array[:, col] min_val = col_array[0] for val in col_array: if val < min_val and not np.isnan(val): min_val = val min_array[col] = min_val return min_array def _projected_xyz_axes(projection): assert projection in ['xy', 'xz', 'yz'], f'invalid projection {projection}' a0 = 'xyz'.index(projection[0]) a1 = 'xyz'.index(projection[1]) return a0, a1 # end of vector_utilities module