diff -r 659ab465152e -r f9788970fa46 geometry.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/geometry.py Wed Jul 29 23:45:53 2020 +0300 @@ -0,0 +1,124 @@ +from math import sqrt, hypot, radians, sin, cos, atan2 + +class Location: + def __init__(self, latitude, longitude): + self.latitude, self.longitude = latitude, longitude + def distance(self, other): + # https://stackoverflow.com/a/365853 + a = sin(radians(self.latitude - other.latitude) / 2) ** 2 + a += sin(radians(self.longitude - other.longitude) / 2) ** 2 * cos(radians(self.latitude)) * cos(radians(other.latitude)) + return 6371 * 2 * atan2(sqrt(a), sqrt(1 - a)) + def __repr__(self): + return '%s(%r, %r)' % (type(self).__name__, self.latitude, self.longitude) + def __str__(self): + return '%.5f, %.5f' % (self.latitude, self.longitude) + @property + def x(self): + return self.latitude + @property + def y(self): + return self.longitude + @property + def link_to_map(self): + return 'http://www.openstreetmap.org/#map=19/%f/%f' % (self.latitude, self.longitude) + +class Ring: + def __init__(self, container): + self.container = container + def __getitem__(self, i): + while i < 0: + i += len(self.container) + while i >= len(self.container): + i -= len(self.container) + return self.container[i] + def __iter__(self): + return iter(self.container) + def __len__(self): + return len(self.container) + +class Polygon: + def __init__(self, *points): + self.points = points + def __repr__(self): + return '%s(%s)' % (type(self).__name__, ', '.join(map(repr, self.points))) + def area(self): + ring = Ring(self.points) + return sum( + ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y + for i in range(len(ring)) + ) / 2 + def circumference(self): + ring = Ring(self.points) + return sum( + sqrt((ring[i + 1].x - ring[i].x)**2 + (ring[i + 1].y - ring[i].y)**2) + for i in range(len(ring)) + ) + def centroid(self): + ring = Ring(self.points) + x = sum( + (ring[i].x + ring[i + 1].x) * (ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y) + for i in range(len(ring)) + ) / 6 / self.area() + y = sum( + (ring[i].y + ring[i + 1].y) * (ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y) + for i in range(len(ring)) + ) / 6 / self.area() + return self.point_type()(x, y) + def point_type(self): + if len(self.points): + return type(self.points[0]) + else: + return Point + def segments(self): + ring = Ring(self.points) + for i in range(len(ring)): + yield LineSegment(ring[i], ring[i + 1]) + def contains_point(self, point): + outer_point = self.point_type()( + min(point.x for point in self.points) - 1, + min(point.y for point in self.points) - 1 + ) + outer_segment = LineSegment(point, outer_point) + intersections = 0 + for segment in self.segments(): + if segment.intersection(outer_segment) is not None: + intersections += 1 + return bool(intersections & 1) + +class LineSegment: + def __init__(self, p1, p2): + self.p1, self.p2 = p1, p2 + def __repr__(self): + return 'LineSegment(%r, %r)' % (self.p1, self.p2) + def length(self): + return hypot(self.p1.x - self.p2.x, self.p1.y - self.p2.y) + def intersection(self, other): + point_type = type(self.p1) + x = (self.p1.x, self.p2.x, other.p1.x, other.p2.x) + y = (self.p1.y, self.p2.y, other.p1.y, other.p2.y) + try: + denominator = (x[0] - x[1]) * (y[2] - y[3]) - (y[0] - y[1]) * (x[2] - x[3]) + Px = ((x[0] * y[1] - y[0] * x[1]) * (x[2] - x[3]) - (x[0] - x[1]) * (x[2] * y[3] - y[2] * x[3])) / denominator + Py = ((x[0] * y[1] - y[0] * x[1]) * (y[2] - y[3]) - (y[0] - y[1]) * (x[2] * y[3] - y[2] * x[3])) / denominator + distance = lambda n: hypot(Px - x[n], Py - y[n]) + if max(distance(0), distance(1)) <= self.length() and max(distance(2), distance(3)) <= other.length(): + return point_type(Px, Py) + else: + return None + except ZeroDivisionError: + return None + +class Point: + def __init__(self, x, y): + self.x, self.y = x, y + def __repr__(self): + return 'Point(%r, %r)' % (self.x, self.y) + +A = Polygon( + Point(2,3), + Point(1,1), + Point(4,0), + Point(6,2), + Point(4,4)) +L1 = LineSegment(Point(1, 1), Point(-1, 5)) +L2 = LineSegment(Point(1, 5), Point(5, 1))