"""
Centrographic measures for point patterns
- documentation
"""
__author__ = "Serge Rey sjsrey@gmail.com"
__all__ = [
"hull",
"mean_center",
"weighted_mean_center",
"manhattan_median",
"std_distance",
"euclidean_median",
"ellipse",
"minimum_rotated_rectangle",
"minimum_bounding_rectangle",
"minimum_bounding_circle",
"dtot",
"_circle",
]
import copy
import math
import sys
import warnings
from collections.abc import Sequence
from functools import singledispatch
import numpy as np
import shapely
from geopandas.base import GeoPandasBase
from libpysal.cg import is_clockwise
from numpy.typing import NDArray
from scipy.optimize import minimize
from scipy.spatial import ConvexHull
not_clockwise = lambda x: not is_clockwise(x)
MAXD = sys.float_info.max
MIND = sys.float_info.min
[docs]
@singledispatch
def minimum_bounding_rectangle(points):
"""
Find minimum bounding rectangle of a point array.
Parameters
----------
points : arraylike
array representing a point pattern
Returns
-------
rectangle
minimum bounding rectangle of a given point pattern
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... [54.46, 8.48],
... ]
... )
Passing an array of coordinates returns a tuple capturing the bounds.
>>> minimum_bounding_rectangle(coords)
(np.float64(8.23), np.float64(7.68), np.float64(98.73), np.float64(92.08))
Passing a GeoPandas object returns a shapely geometry.
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> minimum_bounding_rectangle(geoms)
<POLYGON ((8.23 7.68, 98.73 7.68, 98.73 92.08, 8.23 92.08, 8.23 7.68))>
"""
try:
points = np.asarray(points)
return minimum_bounding_rectangle(points)
except AttributeError as e:
raise NotImplementedError from e
@minimum_bounding_rectangle.register
def _(points: np.ndarray) -> tuple[float, float, float, float]:
points = np.asarray(points)
min_x = min_y = MAXD
max_x = max_y = MIND
x, y = zip(*points, strict=True)
min_x = min(x)
min_y = min(y)
max_x = max(x)
max_y = max(y)
return min_x, min_y, max_x, max_y
@minimum_bounding_rectangle.register
def _(points: GeoPandasBase) -> shapely.Polygon:
return shapely.envelope(shapely.geometrycollections(points.geometry.values))
[docs]
@singledispatch
def minimum_rotated_rectangle(points, return_angle=False):
"""
Compute the minimum rotated rectangle for a point array.
This is the smallest enclosing rectangle (possibly rotated)
for the input point set. It is computed using Shapely.
Parameters
----------
points : arraylike
array representing a point pattern
return_angle : bool
whether to return the angle (in degrees) of the angle between
the horizontal axis of the rectanle and the first side (i.e. length).
Returns
-------
rectangle | tuple(rectangle, angle)
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... [54.46, 8.48],
... ]
... )
Passing an array of coordinates returns an array capturing the corners.
>>> minimum_rotated_rectangle(coords)
array([[ 75.5990843 , -0.72615725],
[ 4.08727852, 30.41752523],
[ 36.40164577, 104.61744544],
[107.91345156, 73.47376296]])
>>> minimum_rotated_rectangle(coords, return_angle=True)
(array([[ 75.5990843 , -0.72615725],
[ 4.08727852, 30.41752523],
[ 36.40164577, 104.61744544],
[107.91345156, 73.47376296]]), 66.466678613503)
Passing a GeoPandas object returns a shapely geometry.
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> minimum_rotated_rectangle(geoms)
<POLYGON ((107.913 73.474, 36.402 104.617, 4.087 30.418, 75.599 -0.726, 107....>
>>> minimum_rotated_rectangle(geoms, return_angle=True)
(<POLYGON ((107.913 73.474, 36.402 104.617, 4.087 30.418, 75.599 -0.726, 107....>, 66.466678613503)
""" # noqa: E501
try:
points = np.asarray(points)
return minimum_rotated_rectangle(points, return_angle=return_angle)
except AttributeError as e:
raise NotImplementedError from e
def _get_mrr_angle(out_points):
angle = (
math.degrees(
math.atan2(
out_points[1][1] - out_points[0][1],
out_points[1][0] - out_points[0][0],
)
)
% 90
)
return angle
@minimum_rotated_rectangle.register
def _(
points: np.ndarray, return_angle: bool = False
) -> NDArray[np.float64] | tuple[NDArray[np.float64], float]:
out_points = shapely.get_coordinates(
shapely.minimum_rotated_rectangle(shapely.multipoints(points))
)[:-1]
if return_angle:
return (out_points[::-1], _get_mrr_angle(out_points))
return out_points[::-1]
@minimum_rotated_rectangle.register
def _(
points: GeoPandasBase, return_angle: bool = False
) -> shapely.Polygon | tuple[shapely.Polygon, float]:
rectangle = shapely.minimum_rotated_rectangle(shapely.multipoints(points.geometry))
if return_angle:
out_points = shapely.get_coordinates(rectangle)[:-1]
return (rectangle, _get_mrr_angle(out_points))
return rectangle
[docs]
@singledispatch
def hull(points):
"""
Find convex hull of a point array.
Parameters
----------
points : arraylike
array representing a point pattern
Returns
-------
rectangle
convex hull of a given point pattern
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... [54.46, 8.48],
... ]
... )
Passing an array of coordinates returns an array capturing the vertices.
>>> hull(coords)
array([[31.01, 81.21],
[ 8.23, 39.93],
[ 9.47, 31.02],
[22.52, 22.39],
[54.46, 8.48],
[79.26, 7.68],
[89.78, 42.53],
[98.73, 77.17],
[65.19, 92.08]])
Passing a GeoPandas object returns a shapely geometry.
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> hull(geoms)
<POLYGON ((79.26 7.68, 54.46 8.48, 22.52 22.39, 9.47 31.02, 8.23 39.93, 31.0...>
"""
try:
points = np.asarray(points)
return hull(points)
except AttributeError as e:
raise NotImplementedError from e
@hull.register
def _(points: np.ndarray) -> NDArray[np.float64]:
h = ConvexHull(points)
return points[h.vertices]
@hull.register
def _(points: GeoPandasBase) -> shapely.Polygon:
return shapely.convex_hull(shapely.multipoints(points.geometry))
[docs]
@singledispatch
def mean_center(points):
"""
Find mean center of a point array.
Parameters
----------
points : arraylike
array representing a point pattern
Returns
-------
center
center of a given point pattern
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... [54.46, 8.48],
... ]
... )
Passing an array of coordinates returns an array capturing the center.
>>> mean_center(coords)
array([52.57166667, 46.17166667])
Passing a GeoPandas object returns a shapely geometry.
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> mean_center(geoms)
<POINT (52.572 46.172)>
"""
try:
points = np.asarray(points)
return mean_center(points)
except AttributeError as e:
raise NotImplementedError from e
@mean_center.register
def _(points: np.ndarray) -> NDArray[np.float64]:
return points.mean(axis=0)
@mean_center.register
def _(points: GeoPandasBase) -> shapely.Point:
coords = shapely.get_coordinates(points.geometry)
return shapely.Point(mean_center(coords))
[docs]
@singledispatch
def weighted_mean_center(points, weights):
"""
Find weighted mean center of a marked point pattern.
Parameters
----------
points : arraylike
array representing a point pattern
weights : arraylike
a series of attribute values of length n.
Returns
-------
center
center of a given point pattern
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... [54.46, 8.48],
... ]
... )
>>> weight = np.arange(1, 13, 1)
Passing an array of coordinates returns an array capturing the center.
>>> weighted_mean_center(coords, weight)
array([59.29448718, 47.52282051])
Passing a GeoPandas object returns a shapely geometry.
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> weighted_mean_center(geoms, weight)
<POINT (59.294 47.523)>
"""
try:
points = np.asarray(points)
return weighted_mean_center(points, weights)
except AttributeError as e:
raise NotImplementedError from e
@weighted_mean_center.register
def _(points: np.ndarray, weights: Sequence) -> NDArray[np.float64]:
points, weights = np.asarray(points), np.asarray(weights)
w = weights * 1.0 / weights.sum()
w.shape = (1, len(points))
return np.dot(w, points)[0]
@weighted_mean_center.register
def _(points: GeoPandasBase, weights) -> shapely.Point:
coords = shapely.get_coordinates(points.geometry)
return shapely.Point(weighted_mean_center(coords, weights))
@manhattan_median.register
def _(points: np.ndarray) -> NDArray[np.float64]:
if not len(points) % 2:
warnings.warn(
"Manhattan Median is not unique for even point patterns.", stacklevel=3
)
return np.median(points, axis=0)
@manhattan_median.register
def _(points: GeoPandasBase) -> shapely.Point:
coords = shapely.get_coordinates(points.geometry)
return shapely.Point(manhattan_median(coords))
[docs]
@singledispatch
def std_distance(points) -> np.float64:
"""
Calculate standard distance of a point array.
Parameters
----------
points : arraylike
array representing a point pattern
Returns
-------
distance
standard distance of a given point pattern
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... ]
... )
Passing an array of coordinates returns a tuple capturing the bounds.
>>> std_distance(coords)
np.float64(40.21575987282547)
The same applies to a GeoPandas object.
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> std_distance(geoms)
np.float64(40.21575987282547)
"""
try:
points = np.asarray(points)
return std_distance(points)
except AttributeError as e:
raise NotImplementedError from e
@std_distance.register
def _(points: np.ndarray) -> np.float64:
n, _ = points.shape
m = points.mean(axis=0)
return np.sqrt(((points * points).sum(axis=0) / n - m * m).sum())
@std_distance.register
def _(points: GeoPandasBase) -> np.float64:
coords = shapely.get_coordinates(points.geometry)
return std_distance(coords)
[docs]
@singledispatch
def ellipse(points, weights=None, method="crimestat",
crimestatCorr=True, degfreedCorr=True):
"""
Computes a weighted standard deviational ellipse for a set of point geometries.
Parameters
----------
points : array-like
Array representing a point pattern.
weights : array-like, optional
Array of weights for each point.
method : str
Correction method to apply. Must be either 'crimestat' or 'yuill'.
crimestatCorr : bool
Whether to apply the CrimeStat correction (used only if `method == 'yuill'`).
degfreedCorr : bool
Whether to apply degrees-of-freedom correction (used only if `method == 'yuill'`).
Apply degrees-of-freedom correction if method == 'yuill'.
Returns
-------
tuple(major_axis_length, minor_axis_length, rotation_angle) | ellipse
Notes
-----
Implements approach from:
https://www.icpsr.umich.edu/CrimeStat/files/CrimeStatChapter.4.pdf
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... ]
... )
Passing an array of coordinates returns a tuple capturing the semi-major axis,
semi-minor axis and clockwise rotation angle of the ellipse.
>>> ellipse(coords)
(np.float64(50.13029459102783), np.float64(37.95222670267865), np.float64(0.3465582095042642))
Passing a GeoPandas object returns a shapely geometry.
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> ellipse(geoms)
<POLYGON ((99.55 66.626, 100.586 63.045, 101.159 59.334, 101.262 55.53, 100....>
"""
try:
points = np.asarray(points)
return ellipse(points, weights, method,
crimestatCorr, degfreedCorr)
except AttributeError as e:
raise NotImplementedError
@ellipse.register
def _(
points: np.ndarray,
weights=None,
method='crimestat',
crimestatCorr=True,
degfreedCorr=True ) -> tuple[float, float, float]:
method = method.lower()
if method not in ("crimestat", "yuill"):
raise ValueError("`method` must be either 'crimestat' or 'yuill'")
x = points[:, 0]
y = points[:, 1]
if weights is None:
weights = np.ones(len(points))
else:
weights = np.asarray(weights)
if len(weights) != len(points):
raise ValueError("weights must have same length as points")
w = weights
sumw = w.sum()
meanx = np.average(x, weights=w)
meany = np.average(y, weights=w)
xm = x - meanx
ym = y - meany
xyw = np.sum(xm * ym * w)
x2w = np.sum(xm**2 * w)
y2w = np.sum(ym**2 * w)
den = 2 * xyw
left = x2w - y2w
right = np.sqrt(left**2 + 4 * xyw**2)
if den == 0:
theta1 = 0
theta2 = np.pi / 2
else:
theta1 = np.arctan(-(left + right) / den)
theta2 = np.arctan(-(left - right) / den)
term1 = np.sum(w * (ym * np.cos(theta1) - xm * np.sin(theta1))**2)
term2 = np.sum(w * (ym * np.cos(theta2) - xm * np.sin(theta2))**2)
sx = np.sqrt(term1 / sumw)
sy = np.sqrt(term2 / sumw)
n = len(points)
if method == "crimestat":
correction = (np.sqrt(2) * np.sqrt(n)) / np.sqrt(n - 2)
sx *= correction
sy *= correction
elif method == "yuill":
if crimestatCorr:
sx *= np.sqrt(2)
sy *= np.sqrt(2)
if degfreedCorr:
sx *= np.sqrt(n) / np.sqrt(n - 2)
sy *= np.sqrt(n) / np.sqrt(n - 2)
if sy > sx:
major_axis, minor_axis = sy, sx
major_angle = theta1
else:
major_axis, minor_axis = sx, sy
major_angle = theta2
return major_axis, minor_axis, major_angle
@ellipse.register
def _(points: GeoPandasBase,
weights=None,
method='crimestat',
crimestatCorr=True,
degfreedCorr=True) -> shapely.Polygon:
coords = shapely.get_coordinates(points.geometry)
major, minor, rotation = ellipse(coords,
weights=weights,
method=method,
crimestatCorr=crimestatCorr,
degfreedCorr=degfreedCorr)
if weights is None:
centre = mean_center(points).buffer(1)
else:
centre = weighted_mean_center(points, weights=weights).buffer(1)
scaled = shapely.affinity.scale(centre, major, minor)
rotated = shapely.affinity.rotate(scaled, rotation, use_radians=True)
return rotated
[docs]
@singledispatch
def dtot(coord, points) -> float:
"""
Sum of Euclidean distances between event points and a selected point.
Parameters
----------
coord
starting point
points : arraylike
array representing a point pattern
Returns
-------
distance
sum of Euclidean distances.
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
>>> import shapely
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... ]
... )
Passing an array of coordinates returns a tuple capturing the bounds.
>>> point = [30.78, 60.10]
>>> dtot(point, coords)
np.float64(465.3617957739617)
The same applies to a GeoPandas object.
>>> point = shapely.Point(point)
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> dtot(point, geoms)
np.float64(465.3617957739617)
"""
try:
coord = np.asarray(coord)
points = np.asarray(points)
return dtot(coord, points)
except AttributeError as e:
raise NotImplementedError from e
@dtot.register
def _(coord: np.ndarray, points: np.ndarray) -> float:
xd = points[:, 0] - coord[0]
yd = points[:, 1] - coord[1]
d = np.sqrt(xd * xd + yd * yd).sum()
return d
@dtot.register
def _(coord: shapely.Point, points: GeoPandasBase) -> float:
coord = shapely.get_coordinates(coord).flatten()
points = shapely.get_coordinates(points.geometry)
return dtot(coord, points)
@euclidean_median.register
def _(points: np.ndarray) -> NDArray[np.float64]:
start = mean_center(points)
res = minimize(dtot, start, args=(points,))
return res["x"]
@euclidean_median.register
def _(points: GeoPandasBase) -> shapely.Point:
coords = shapely.get_coordinates(points.geometry)
return shapely.Point(euclidean_median(coords))
[docs]
@singledispatch
def minimum_bounding_circle(points):
"""
Implements Skyum (1990)'s algorithm for the minimum bounding circle in R^2.
Store points clockwise.
Find p in S that maximizes angle(prec(p), p, succ(p) THEN radius(prec(
p), p, succ(p)). This is also called the lexicographic maximum, and is the last
entry of a list of (radius, angle) in lexicographical order.
* If angle(prec(p), p, succ(p)) <= 90 degrees, then finish.
* If not, remove p from set.
Parameters
----------
points : arraylike
array representing a point pattern
Returns
-------
circle
minimum bounding circle
Examples
--------
>>> import numpy as np
>>> import geopandas as gpd
Create an array of point coordinates.
>>> coords = np.array(
... [
... [66.22, 32.54],
... [22.52, 22.39],
... [31.01, 81.21],
... [9.47, 31.02],
... [30.78, 60.10],
... [75.21, 58.93],
... [79.26, 7.68],
... [8.23, 39.93],
... [98.73, 77.17],
... [89.78, 42.53],
... [65.19, 92.08],
... ]
... )
Passing an array of coordinates returns an tuple of (x, y), radius.
>>> minimum_bounding_circle(coords)
((55.244520477497474, 51.88135107645883), np.float64(50.304102155590726))
Passing a GeoPandas object returns a shapely geometry.
>>> geoms = gpd.GeoSeries.from_xy(*coords.T)
>>> minimum_bounding_circle(geoms)
<POLYGON ((105.549 51.881, 105.306 46.951, 104.582 42.068, 103.383 37.279, 1...>
"""
try:
points = np.asarray(points)
return minimum_bounding_circle(points)
except AttributeError as e:
raise NotImplementedError from e
@minimum_bounding_circle.register
def _(points: np.ndarray) -> tuple[tuple[float, float], float]:
try:
from numba import njit
HAS_NUMBA = True
except ImportError:
HAS_NUMBA = False
points = hull(points)
if not_clockwise(points):
points = points[::-1]
if not_clockwise(points):
raise Exception("Points are neither clockwise nor counterclockwise")
POINTS = copy.deepcopy(points)
if HAS_NUMBA: # noqa: SIM108
circ = _skyum_numba(POINTS)[0]
else:
circ = _skyum_lists(POINTS)[0]
return (circ[1], circ[2]), circ[0]
@minimum_bounding_circle.register
def _(points: GeoPandasBase) -> shapely.Polygon:
coords = shapely.get_coordinates(points.geometry)
(x, y), r = minimum_bounding_circle(coords)
return shapely.Point(x, y).buffer(r)
def _skyum_lists(points):
points = points.tolist()
removed = []
i = 0
while True:
angles = [
_angle(
_prec(p, points),
p,
_succ(p, points),
)
for p in points
]
circles = [
_circle(
_prec(p, points),
p,
_succ(p, points),
)
for p in points
]
radii = [c[0] for c in circles]
lexord = np.lexsort((radii, angles)) # confusing as hell defaults...
lexmax = lexord[-1]
candidate = (
_prec(points[lexmax], points),
points[lexmax],
_succ(points[lexmax], points),
)
if angles[lexmax] <= (np.pi / 2.0):
# print("Constrained by points: {}".format(candidate))
return _circle(*candidate), points, removed, candidate
else:
try:
removed.append((points.pop(lexmax), i))
except IndexError:
raise Exception("Construction of Minimum Bounding Circle failed!")
i += 1
try:
from numba import njit, boolean
HAS_NUMBA = True
@njit(fastmath=True)
def _skyum_numba(points):
i = 0
complete = False
while not complete:
complete, points, candidate, circle = _skyum_iteration(points)
if complete:
return circle, points, None, candidate
@njit(fastmath=True)
def _skyum_iteration(points):
points = points.reshape(-1, 2)
n = points.shape[0]
angles = np.empty((n,))
circles = np.empty((n, 3))
for i in range(n):
p = points[(i - 1) % n]
q = points[i % n]
r = points[(i + 1) % n]
angles[i] = _angle(p, q, r)
circles[i] = _circle(p, q, r)
radii = circles[:, 0]
# workaround for no lexsort in numba
angle_argmax = angles.argmax()
angle_max = angles[angle_argmax]
# the maximum radius for the largest angle
lexmax = (radii * (angles == angle_max)).argmax()
if angles[lexmax] <= (np.pi / 2.0):
return True, points, lexmax, circles[lexmax]
else:
mask = np.ones((n,), dtype=boolean)
mask[lexmax] = False
new_points = points[mask, :]
return False, new_points, lexmax, circles[lexmax]
except ModuleNotFoundError:
def njit(func, **kwargs):
return func
@njit
def _angle(p, q, r):
p = np.asarray(p)
q = np.asarray(q)
r = np.asarray(r)
pq = p - q
rq = r - q
magnitudes = np.linalg.norm(pq) * np.linalg.norm(rq)
return np.abs(np.arccos(np.dot(pq, rq) / magnitudes))
def _prec(p, l):
"""
retrieve the predecessor of p in list l
"""
pos = l.index(p)
if pos - 1 < 0:
return l[-1]
else:
return l[pos - 1]
def _succ(p, l):
"""
retrieve the successor of p in list l
"""
pos = l.index(p)
if pos + 1 >= len(l):
return l[0]
else:
return l[pos + 1]
@njit
def _euclidean_distance(px, py, qx, qy):
return np.sqrt((px - qx) ** 2 + (py - qy) ** 2)
@njit
def _circle(p, q, r, dmetric=_euclidean_distance):
"""
Returns (radius, (center_x, center_y)) of the circumscribed circle by the
triangle pqr.
note, this does not assume that p!=q!=r
"""
px, py = p
qx, qy = q
rx, ry = r
angle = np.abs(_angle(p, q, r))
if np.abs(angle - np.pi) < 1e-5: # angle is pi
radius = _euclidean_distance(px, py, rx, ry) / 2.0
center_x = (px + rx) / 2.0
center_y = (py + ry) / 2.0
elif np.abs(angle) < 1e-5: # angle is zero
radius = _euclidean_distance(px, py, qx, qy) / 2.0
center_x = (px + qx) / 2.0
center_y = (py + qy) / 2.0
else:
D = 2 * (px * (qy - ry) + qx * (ry - py) + rx * (py - qy))
center_x = (
(px**2 + py**2) * (qy - ry)
+ (qx**2 + qy**2) * (ry - py)
+ (rx**2 + ry**2) * (py - qy)
) / float(D)
center_y = (
(px**2 + py**2) * (rx - qx)
+ (qx**2 + qy**2) * (px - rx)
+ (rx**2 + ry**2) * (qx - px)
) / float(D)
radius = _euclidean_distance(center_x, center_y, px, py)
return radius, center_x, center_y