Mercurial > buscompiler / file comparison
comparison: geometry.py
geometry.py
- changeset 1
- f9788970fa46
equal
deleted
inserted
replaced
| 0:659ab465152e | 1:f9788970fa46 |
|---|---|
| 1 from math import sqrt, hypot, radians, sin, cos, atan2 | |
| 2 | |
| 3 class Location: | |
| 4 def __init__(self, latitude, longitude): | |
| 5 self.latitude, self.longitude = latitude, longitude | |
| 6 def distance(self, other): | |
| 7 # https://stackoverflow.com/a/365853 | |
| 8 a = sin(radians(self.latitude - other.latitude) / 2) ** 2 | |
| 9 a += sin(radians(self.longitude - other.longitude) / 2) ** 2 * cos(radians(self.latitude)) * cos(radians(other.latitude)) | |
| 10 return 6371 * 2 * atan2(sqrt(a), sqrt(1 - a)) | |
| 11 def __repr__(self): | |
| 12 return '%s(%r, %r)' % (type(self).__name__, self.latitude, self.longitude) | |
| 13 def __str__(self): | |
| 14 return '%.5f, %.5f' % (self.latitude, self.longitude) | |
| 15 @property | |
| 16 def x(self): | |
| 17 return self.latitude | |
| 18 @property | |
| 19 def y(self): | |
| 20 return self.longitude | |
| 21 @property | |
| 22 def link_to_map(self): | |
| 23 return 'http://www.openstreetmap.org/#map=19/%f/%f' % (self.latitude, self.longitude) | |
| 24 | |
| 25 class Ring: | |
| 26 def __init__(self, container): | |
| 27 self.container = container | |
| 28 def __getitem__(self, i): | |
| 29 while i < 0: | |
| 30 i += len(self.container) | |
| 31 while i >= len(self.container): | |
| 32 i -= len(self.container) | |
| 33 return self.container[i] | |
| 34 def __iter__(self): | |
| 35 return iter(self.container) | |
| 36 def __len__(self): | |
| 37 return len(self.container) | |
| 38 | |
| 39 class Polygon: | |
| 40 def __init__(self, *points): | |
| 41 self.points = points | |
| 42 def __repr__(self): | |
| 43 return '%s(%s)' % (type(self).__name__, ', '.join(map(repr, self.points))) | |
| 44 def area(self): | |
| 45 ring = Ring(self.points) | |
| 46 return sum( | |
| 47 ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y | |
| 48 for i in range(len(ring)) | |
| 49 ) / 2 | |
| 50 def circumference(self): | |
| 51 ring = Ring(self.points) | |
| 52 return sum( | |
| 53 sqrt((ring[i + 1].x - ring[i].x)**2 + (ring[i + 1].y - ring[i].y)**2) | |
| 54 for i in range(len(ring)) | |
| 55 ) | |
| 56 def centroid(self): | |
| 57 ring = Ring(self.points) | |
| 58 x = sum( | |
| 59 (ring[i].x + ring[i + 1].x) * (ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y) | |
| 60 for i in range(len(ring)) | |
| 61 ) / 6 / self.area() | |
| 62 y = sum( | |
| 63 (ring[i].y + ring[i + 1].y) * (ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y) | |
| 64 for i in range(len(ring)) | |
| 65 ) / 6 / self.area() | |
| 66 return self.point_type()(x, y) | |
| 67 def point_type(self): | |
| 68 if len(self.points): | |
| 69 return type(self.points[0]) | |
| 70 else: | |
| 71 return Point | |
| 72 def segments(self): | |
| 73 ring = Ring(self.points) | |
| 74 for i in range(len(ring)): | |
| 75 yield LineSegment(ring[i], ring[i + 1]) | |
| 76 def contains_point(self, point): | |
| 77 outer_point = self.point_type()( | |
| 78 min(point.x for point in self.points) - 1, | |
| 79 min(point.y for point in self.points) - 1 | |
| 80 ) | |
| 81 outer_segment = LineSegment(point, outer_point) | |
| 82 intersections = 0 | |
| 83 for segment in self.segments(): | |
| 84 if segment.intersection(outer_segment) is not None: | |
| 85 intersections += 1 | |
| 86 return bool(intersections & 1) | |
| 87 | |
| 88 class LineSegment: | |
| 89 def __init__(self, p1, p2): | |
| 90 self.p1, self.p2 = p1, p2 | |
| 91 def __repr__(self): | |
| 92 return 'LineSegment(%r, %r)' % (self.p1, self.p2) | |
| 93 def length(self): | |
| 94 return hypot(self.p1.x - self.p2.x, self.p1.y - self.p2.y) | |
| 95 def intersection(self, other): | |
| 96 point_type = type(self.p1) | |
| 97 x = (self.p1.x, self.p2.x, other.p1.x, other.p2.x) | |
| 98 y = (self.p1.y, self.p2.y, other.p1.y, other.p2.y) | |
| 99 try: | |
| 100 denominator = (x[0] - x[1]) * (y[2] - y[3]) - (y[0] - y[1]) * (x[2] - x[3]) | |
| 101 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 | |
| 102 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 | |
| 103 distance = lambda n: hypot(Px - x[n], Py - y[n]) | |
| 104 if max(distance(0), distance(1)) <= self.length() and max(distance(2), distance(3)) <= other.length(): | |
| 105 return point_type(Px, Py) | |
| 106 else: | |
| 107 return None | |
| 108 except ZeroDivisionError: | |
| 109 return None | |
| 110 | |
| 111 class Point: | |
| 112 def __init__(self, x, y): | |
| 113 self.x, self.y = x, y | |
| 114 def __repr__(self): | |
| 115 return 'Point(%r, %r)' % (self.x, self.y) | |
| 116 | |
| 117 A = Polygon( | |
| 118 Point(2,3), | |
| 119 Point(1,1), | |
| 120 Point(4,0), | |
| 121 Point(6,2), | |
| 122 Point(4,4)) | |
| 123 L1 = LineSegment(Point(1, 1), Point(-1, 5)) | |
| 124 L2 = LineSegment(Point(1, 5), Point(5, 1)) |
