"""Planar Euclidean geometry utils.
"""
import numpy as np
import ddg.math.complex as gcomplex
sqrt = np.math.sqrt
cos = np.math.cos
sin = np.math.sin
# --- quadrilaterals ---
[docs]def intersect_diags(x1, x2, x3, x4):
"""
For a 2D-quadrilateral (x1,x2,x3,x4)
get the intersection point of the diagonals.
Please refer to :py:func:`ddg.math.complex.intersect_diags`.
Parameters
----------
x1, x2, x3, x4 : numpy.ndarray of shape(2,)
Coordinates of the given quadrilateral.
Returns
-------
numpy.ndarray of shape(2,)
Coordinate of diagonal intersection
"""
x1, x2, x3, x4 = map(gcomplex.to_complex, [x1, x2, x3, x4])
return gcomplex.to_array(gcomplex.intersect_diags(x1, x2, x3, x4))
# Wrapper reverted due to error with test test_laplace invariants.
[docs]def intersect_edges(x1, x2, x3, x4):
r"""
For a 2D-quadrilateral (x1, x2, x3, x4)
get the intersection points of opposing edges.
Will return a array(NaN) if intersection does not exist. ::
z2
/ \
/ \
x4-----x3--
| | \
| | z1
| | /
x1-----x2--
Parameters
----------
x1, x2, x3, x4 : numpy.ndarray of shape(2,)
Coordinates of the given quadrilateral.
Returns
-------
z1, z2 : numpy.ndarray of shape(2,)
Coordinates of diagonal intersections
"""
try:
z1 = intersectionCramer(line2D(x1, x2), line2D(x4, x3))
except ValueError:
z1 = np.nan
try:
z2 = intersectionCramer(line2D(x1, x4), line2D(x2, x3))
except ValueError:
z2 = np.nan
return np.array(z1), np.array(z2)
[docs]def christoffel_dual_quad(x1, x2, x3, x4):
"""
A Christoffel dual quadrilateral of the given quadrilateral.
Corresponding edges and non-corresponding diagonals are parallel.
This determines the dual quadrilateral up to scaling.
Parameters
----------
x1, x2, x3, x4 : numpy.ndarray of shape (2,)
Coordinates of the given quadrilateral.
Returns
-------
X1, X2, X3, X4 : numpy.ndarray of shape (2,)
Coordinates of the dual quadrilateral.
"""
x1, x2, x3, x4 = map(np.array, [x1, x2, x3, x4])
m = intersect_diags(x1, x2, x3, x4)
alpha = np.linalg.norm(x1 - m)
beta = np.linalg.norm(x2 - m)
gamma = np.linalg.norm(x3 - m)
delta = np.linalg.norm(x4 - m)
if np.dot(x3 - m, x1 - m) < 0:
alpha = -alpha
if np.dot(x2 - m, x4 - m) < 0:
delta = -delta
X1 = np.zeros(x1.size)
X2 = (x2 - x1) / (alpha * beta)
X3 = X2 + (x3 - x2) / (beta * gamma)
X4 = (x4 - x1) / (alpha * delta)
return (X1, X2, X3, X4)
# --- angles ---
[docs]def corner_angle(x1, x2, x3):
"""
Angle of the corner (x1, x2, x3) at the vertex x2.
Computes the angle between vectors x1-x2 and x3-x2 in (0,2*pi].
Parameters
----------
x1, x2, x3 : numpy.ndarray of shape (2,)
Coordinate of points forming the vertex.
Returns
-------
float
Vertex angle.
Notes
-----
Only reliable for convex quads.
"""
u = gcomplex.to_complex(x3 - x2)
v = gcomplex.to_complex(x1 - x2)
alpha = np.angle(v / u)
return alpha if alpha >= 0 else 2 * np.pi + alpha
[docs]def corner_angle_signed(x1, x2, x3):
"""
Signed angle of the corner (x1, x2, x3) at the vertex x2.
Computes the angle between vectors x1-x2 and x3-x2 in (-pi,pi].
Parameters
----------
x1, x2, x3 : numpy.ndarray of shape (2,)
Coordinate of points forming the vertex.
Returns
-------
float
Signed vertex angle.
"""
u = gcomplex.to_complex(x3 - x2)
v = gcomplex.to_complex(x1 - x2)
return np.angle(v / u)
[docs]def cyclic_order(points, center):
"""Determines the cyclic order of the points beginning with the first with
respect to the center.
Parameters
----------
points : list of numpy.ndarray of shape (2,)
List of points.
center : numpy.ndarray of shape (2,)
Coordinate of the center.
Returns
-------
list of numpy.ndarray of shape (2,)
Cyclic order of the points with respect to the center.
Notes
-----
Useful for angle calculation in non-embedded case.
"""
offset = np.angle(gcomplex.to_complex(points[0] - center))
pos = lambda alpha: alpha if alpha >= 0 else alpha + 2 * np.pi
angle_function = lambda point: pos(
np.angle(gcomplex.to_complex(point - center)) - offset
)
return sorted(points, key=angle_function)
[docs]def rotation(alpha):
"""
Returns 2x2 rotation matrix ::
cos(alpha) | -sin(alpha)
A = ------------------------
sin(alpha) | cos(alpha)
Parameters
----------
alpha : float
Rotation angle.
Returns
-------
numpy.ndarray of shape (2, 2)
Examples
--------
>>> import numpy as np
>>> from ddg.math.euclidean2d import rotation
>>> rotation(0)
array([[ 1., -0.],
[ 0., 1.]])
>>> rotation(np.pi / 2)
array([[ 6.123234e-17, -1.000000e+00],
[ 1.000000e+00, 6.123234e-17]])
"""
return np.array([[np.cos(alpha), -np.sin(alpha)], [np.sin(alpha), np.cos(alpha)]])
# --- circles
[docs]def circumcenter(x1, x2, x3):
"""
Get center of circle through three points.
This is a real variant of
:py:func:`ddg.math.complex.circumcenter`.
Parameters
----------
x1, x2, x3 : numpy.ndarray of shape (2,)
Coordinates of three points in a circle.
Returns
-------
numpy.ndarray of shape (2,)
Coordinate of center.
"""
z1, z2, z3 = map(gcomplex.to_complex, [x1, x2, x3])
center = gcomplex.circumcenter(z1, z2, z3)
return gcomplex.to_array(center)
[docs]def circumradius(x1, x2, x3):
"""
Get radius of circle through three points.
This is a real variant of
:py:func:`ddg.math.complex.circumradius`.
Parameters
----------
x1, x2, x3 : numpy.ndarray of shape (2,)
Coordinates of three points in a circle.
Returns
-------
float
Radius of circle.
"""
z1, z2, z3 = map(gcomplex.to_complex, [x1, x2, x3])
return gcomplex.circumradius(z1, z2, z3)
[docs]def fourth_point_from_cross_ratio(x, x1, x2, q):
"""
Get fourth point on a quadrilateral from real cross-ratio.
Applies the real variant of the complex cross ratio
(See :py:func:`ddg.math.complex.fourth_point_from_cross_ratio`).
Parameters
----------
x, x1, x2 : numpy.ndarray of shape (2,)
Coordinates of points with order (x, x_1, ..., x_2).
q : complex
Cross ration of points.
Returns
-------
array_like
Corrdinate of fourth point.
Notes
-----
Not exact with fractions right now."""
z, z1, z2 = map(gcomplex.to_complex, [x, x1, x2])
z12 = gcomplex.fourth_point_from_cross_ratio(z, z1, z2, q)
return gcomplex.to_array(z12)
# --- triangles ---
[docs]def triangle_area(a, b, c):
"""
Compute area of triangle from side lengths.
Numerically stable version of Heron's formula
for the area of a triangle with side lengths a,b,c.
Parameters
----------
a, b, c : float
Side lengths
Returns
-------
float
Area of triangle
"""
return sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c))) / 4
[docs]def line2D(x1, x2):
"""
Gives the coefficients of the line equation describing the line through
given points x1 and x2.
Parameters
----------
x1 : numpy.ndarray of shape (2,)
Point x1 given by its Euclidean 2D coordinates in 2-dim array.
x2 : numpy.ndarray of shape (2,)
Point x1 given by its Euclidean 2D coordinates in 2-dim array.
Returns
-------
[A, B, C] : list of three floats
The coefficients of the line equation `Ax + By = C` describing the line
through x1 and x2.
"""
A = x1[1] - x2[1]
B = x2[0] - x1[0]
C = -(x1[0] * x2[1] - x2[0] * x1[1])
return A, B, C
[docs]def intersection(l1, l2):
"""
Calculates the intersection of two lines that are given by coefficients of
their line equations.
To do so the linear system `Ax = b` of the line equations is solved via
numpy.linalg.solver.
Parameters
----------
l1, l2 : [float, float, float]
The two lines whose intersection shall be calculated. Both given by the
coefficients of their line equations `Ai*x + Bi*y = Ci` where i=1,2.
Returns
-------
[x,y] : np.array of two floats
The Euclidean 2D coordinates of the intersection point.
"""
a = np.array([[l1[0], l1[1]], [l2[0], l2[1]]])
b = np.array([l1[2], l2[2]])
x = np.linalg.solve(a, b)
return x
[docs]def intersectionCramer(l1, l2):
"""
Calculates the intersection of two lines that are given by coefficients of
their line equations via Cramer's rule.
Parameters
----------
l1, l2 : [float, float, float]
The two lines whose intersection shall be calculated. Both given by the
coefficients of their line equations `Ai*x + Bi*y = Ci` where i=1,2.
Returns
-------
[x,y] : array of two floats
The Euclidean 2D coordinates of the intersection point.
Raises
------
ValueError
If the given lines do not intersect.
"""
D = l1[0] * l2[1] - l1[1] * l2[0]
Dx = l1[2] * l2[1] - l1[1] * l2[2]
Dy = l1[0] * l2[2] - l1[2] * l2[0]
if D != 0:
x = Dx / D
y = Dy / D
return x, y
else:
raise ValueError("Lines do not intersect.")
[docs]def intersection_in_barycentric_coords(x1, x2, y1, y2):
"""
Gives the intersection of two lines in barycentric coordinates relative to
the given points which the lines run through.
Parameters
----------
x1, x2 : numpy.ndarray of shape (2,)
The two points the first line for intersecting runs through,
given by their Euclidean 2D coordinates.
y1, y2 : numpy.ndarray of shape (2,)
The two points the second line for intersecting runs through,
given by their Euclidean 2D coordinates.
Returns
-------
x : np.array[t1, s1, t2, s2] where ti, si float, i=1,2
The intersection of the lines given as the barycentric coordinates
corresponding to solving the equation `t1*x1 + s1*x2 = t2*y1 + s2*y2`.
Notes
-----
We establish a homogeneous system and solve it via numpy.
"""
matrix = np.array(
[
[x1[0], x1[1], 1, 1],
[x2[0], x2[1], 1, 1],
[-y1[0], -y1[1], -1, 0],
[-y2[0], -y2[1], -1, 0],
]
).T
vec = np.array([0, 0, 0, 1])
x = np.linalg.solve(matrix, vec)
return x