"""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 diagonals_intersection(a, b, c, d):
a, b, c, d = map( np.array, [a, b, c, d] )
A = np.transpose( np.array([c - a, d - b]) )
t = np.linalg.lstsq(A, b-a, rcond=None)[0][0]
return a + t*(c-a)
[docs]def intersect_edges(a, b, c, d):
a, b, c, d = map( np.array, [a, b, c, d] )
A = np.transpose( np.array([b - a, d - c]) )
t = np.linalg.lstsq(A, c-a, rcond=None)[0][0]
return a + t*(b-a)
[docs]def dual_quadrilateral(a, b, c, d):
a, b, c, d = map( np.array, [a, b, c, d] )
m = diagonals_intersection(a, b, c, d)
alpha = np.linalg.norm( a - m )
beta = np.linalg.norm( b - m )
gamma = np.linalg.norm( c - m )
delta = np.linalg.norm( d - m )
if np.dot( c - m, a - m ) < 0:
alpha = - alpha
if np.dot( b - m, d - m ) < 0:
delta = - delta
A = np.zeros( a.size )
B = (b - a) / (alpha * beta)
C = B + (c - b) / (beta * gamma)
D = (d - a) / (alpha * delta)
return (A, B, C, D)
# --- lines ---
[docs]def distance_lines(a, v, b, w):
"""Distance between the lines a + sv and b + tw."""
n = np.cross(v, w)
n = n / np.linalg.norm(n)
return np.abs( np.dot(n, b-a) )
# --- angles ---
[docs]def vertex_angle(x1, x2, x3):
"""angle of (x1,x2,x3) on the left to the line (not the vector)
only reliable for konvex 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 vertex_angle_with_sign(x1, x2, x3):
"""like vertex_angle just with sign instead of +pi"""
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 wrt the center
maybe usefull 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 )
# --- circles
[docs]def circumcircle_center(x, y, z):
D = 2*( x[0]*(y[1]-z[1]) + y[0]*(z[1]-x[1]) + z[0]*(x[1]-y[1]))
s_x = x[0]**2 + x[1]**2
s_y = y[0]**2 + y[1]**2
s_z = z[0]**2 + z[1]**2
u = ( s_x*(y[1]-z[1]) + s_y*(z[1]-x[1]) + s_z*(x[1]-y[1]) ) / D
v = ( s_x*(z[0]-y[0]) + s_y*(x[0]-z[0]) + s_z*(y[0]-x[0]) ) / D
return np.array([u,v])
[docs]def circumcircle_radius(x, y, z):
norm = lambda x: np.linalg.norm( map(float, x) ) # convert to float ( norm doesent work on fractions... )
l_xy, l_yz, l_zx = map( norm, (x-y,y-z,z-x) )
return (l_xy*l_yz*l_zx) / ( 4*triangle_area(l_xy, l_yz, l_zx) )
[docs]def coordinates_on_circle_from_cross_ratio(x, x1, x2, q):
"""get fourth point on a circle from real cross-ratio
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):
"""Numerically stable version of Heron's formula
for the area of a tringagle with side lengths a,b,c. """
return sqrt( (a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)) ) / 4
[docs]def line2D(p1, p2):
"""
Gives the coefficients of the line equation describing the line through
given points p1 and p2.
Parameters
----------
p1 : np.array([x1, y1])
Point p1 given by its Euclidean 2D coordinates in 2-dim array.
p2 : np.array([x2, y2])
Point p1 given by its Euclidean 2D coordinates in 2-dim array.
Returns
-------
[A, B, C] : array of three floats
The coefficients of the line equation `Ax + By = C` describing the line
through p1 and p2.
"""
A = (p1[1] - p2[1])
B = (p2[0] - p1[0])
C = - (p1[0] * p2[1] - p2[0] * p1[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(a, b, p, q):
"""
Gives the intersection of two lines in barycentric coordinates relative to
the given points which the lines run through.
Parameters
----------
a, b : np.array([float, float])
The two points the first line for intersecting runs through,
given by their Euclidean 2D coordinates.
p, q : np.array([float, float])
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*a + s1*b = t2*p + s2*q`.
Notes
-----
We establish a homogeneous system and solve it via numpy.
"""
matrix = np.array([[a[0], a[1], 1, 1],
[b[0], b[1], 1, 1],
[-p[0], -p[1], -1, 0],
[-q[0], -q[1], -1, 0]]).T
vec = np.array([0, 0, 0, 1])
x = np.linalg.solve(matrix, vec)
return x