diff options
| -rw-r--r-- | lib/polygon.py | 27 | ||||
| -rw-r--r-- | lib/test/test.py | 15 | ||||
| -rw-r--r-- | lib/test/test_intersection.py | 14 | ||||
| -rw-r--r-- | lib/test/test_union.py | 15 |
4 files changed, 54 insertions, 17 deletions
diff --git a/lib/polygon.py b/lib/polygon.py index 80ef131..689e272 100644 --- a/lib/polygon.py +++ b/lib/polygon.py @@ -60,7 +60,7 @@ class SphericalPolygon(object): therefore we need a way of specifying which is which. """ - def __init__(self, points, inside): + def __init__(self, points, inside=None): r""" Parameters ---------- @@ -73,8 +73,9 @@ class SphericalPolygon(object): Four points are needed to define a triangle, since the polygon must be closed. - inside : An (*x*, *y*, *z*) triple - This point must be inside the polygon. + inside : An (*x*, *y*, *z*) triple, optional + This point must be inside the polygon. If not provided, the + mean of the points will be used. """ if len(points) == 0: pass @@ -83,7 +84,10 @@ class SphericalPolygon(object): else: assert np.array_equal(points[0], points[-1]), 'Polygon is not closed' self._points = np.asanyarray(points) - self._inside = np.asanyarray(inside) + if inside is None: + self._inside = np.mean(points[:-1], axis=0) + else: + self._inside = np.asanyarray(inside) # TODO: Detect self-intersection and fix @@ -408,25 +412,30 @@ class SphericalPolygon(object): if len(self._points) < 3: return np.array(0.0) + points = self._points.copy() + # Rotate polygon so that center of polygon is at north pole - centroid = np.mean(self._points, axis=0) + centroid = np.mean(points, axis=0) + centroid = vector.normalize_vector(*centroid) points = self._points - (centroid + np.array([0, 0, 1])) + vector.normalize_vector( + points[:, 0], points[:, 1], points[:, 2], inplace=True) X = [] Y = [] for A, B in zip(points[:-1], points[1:]): length = great_circle_arc.length(A, B, degrees=True) interp = great_circle_arc.interpolate(A, B, length * 4) - x, y = vector.equal_area_proj(interp[:, 0], - interp[:, 1], - interp[:, 2]) + x, y, z = vector.normalize_vector( + interp[:, 0], interp[:, 1], interp[:, 2], inplace=True) + x, y = vector.equal_area_proj(x, y, z) X.extend(x) Y.extend(y) X = np.array(X) Y = np.array(Y) - return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) / 2.0) + return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) * 0.5 * np.pi) def union(self, other): """ diff --git a/lib/test/test.py b/lib/test/test.py index 34b5afa..b551090 100644 --- a/lib/test/test.py +++ b/lib/test/test.py @@ -260,3 +260,18 @@ def test_cone(): ra = random.randrange(-180, 180) dec = random.randrange(20, 90) cone = polygon.SphericalPolygon.from_cone(ra, dec, 8, steps=64) + + +def test_area(): + triangles = [ + ([[90, 0], [0, -45], [0, 45], [90, 0]], np.pi * 0.5), + ([[90, 0], [0, -22.5], [0, 22.5], [90, 0]], np.pi * 0.25), + ([[90, 0], [0, -11.25], [0, 11.25], [90, 0]], np.pi * 0.125) + ] + + for tri, area in triangles: + tri = np.array(tri) + x, y, z = vector.radec_to_vector(tri[:,1], tri[:,0]) + points = np.dstack((x, y, z))[0] + poly = polygon.SphericalPolygon(points) + calc_area = poly.area() diff --git a/lib/test/test_intersection.py b/lib/test/test_intersection.py index d25ce2e..cf4c3d5 100644 --- a/lib/test/test_intersection.py +++ b/lib/test/test_intersection.py @@ -29,11 +29,18 @@ class intersection_test: def __call__(self, func): @functools.wraps(func) def run(*args, **kwargs): + if GRAPH_MODE: + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + polys = func(*args, **kwargs) intersections = [] num_permutations = math.factorial(len(polys)) step_size = int(max(float(num_permutations) / 20.0, 1.0)) + + areas = [x.area() for x in polys] + if GRAPH_MODE: print("%d permutations" % num_permutations) for method in ('parallel', 'serial'): @@ -48,10 +55,7 @@ class intersection_test: intersection = polygon.SphericalPolygon.multi_intersection( permutation, method=method) intersections.append(intersection) - areas = [x.area() for x in permutation] intersection_area = intersection.area() - assert np.all(intersection_area < areas) - if GRAPH_MODE: fig = plt.figure() m = Basemap(projection=self._proj, @@ -67,10 +71,12 @@ class intersection_test: plt.savefig(filename) fig.clear() + assert np.all(intersection_area < areas) + lengths = np.array([len(x._points) for x in intersections]) assert np.all(lengths == [lengths[0]]) areas = np.array([x.area() for x in intersections]) - assert_array_almost_equal(areas, areas[0], decimal=3) + assert_array_almost_equal(areas, areas[0], decimal=1) return run diff --git a/lib/test/test_union.py b/lib/test/test_union.py index fda91cb..f867e9d 100644 --- a/lib/test/test_union.py +++ b/lib/test/test_union.py @@ -29,6 +29,10 @@ class union_test: def __call__(self, func): @functools.wraps(func) def run(*args, **kwargs): + if GRAPH_MODE: + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + polys = func(*args, **kwargs) unions = [] @@ -36,8 +40,11 @@ class union_test: step_size = int(max(float(num_permutations) / 20.0, 1.0)) if GRAPH_MODE: print("%d permutations" % num_permutations) + + areas = [x.area() for x in polys] + for i, permutation in enumerate( - itertools.islice( + itertools.islice( itertools.permutations(polys), None, None, step_size)): filename = '%s_union_%04d.svg' % ( @@ -47,9 +54,7 @@ class union_test: union = polygon.SphericalPolygon.multi_union( permutation) unions.append(union) - areas = [x.area() for x in permutation] union_area = union.area() - assert np.all(union_area >= areas) if GRAPH_MODE: fig = plt.figure() @@ -66,10 +71,12 @@ class union_test: plt.savefig(filename) fig.clear() + assert np.all(union_area >= areas) + lengths = np.array([len(x._points) for x in unions]) assert np.all(lengths == [lengths[0]]) areas = np.array([x.area() for x in unions]) - # assert_array_almost_equal(areas, areas[0], 1) + assert_array_almost_equal(areas, areas[0], 1) return run |