diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/graph.py | 78 | ||||
| -rw-r--r-- | lib/great_circle_arc.py | 8 | ||||
| -rw-r--r-- | lib/polygon.py | 81 | ||||
| -rw-r--r-- | lib/test/test_intersection.py | 2 | ||||
| -rw-r--r-- | lib/test/test_union.py | 81 | ||||
| -rw-r--r-- | lib/vector.py | 34 |
6 files changed, 175 insertions, 109 deletions
diff --git a/lib/graph.py b/lib/graph.py index 25fc80c..97c12b4 100644 --- a/lib/graph.py +++ b/lib/graph.py @@ -166,6 +166,7 @@ class Graph: self._nodes = set() self._edges = set() self._source_polygons = set() + self._start_node = None self.add_polygons(polygons) @@ -204,6 +205,8 @@ class Graph: self._source_polygons.add(polygon) start_node = nodeA = self.add_node(points[0], [polygon]) + if self._start_node is None: + self._start_node = start_node for i in range(1, len(points) - 1): # Don't create self-pointing edges if not np.array_equal(points[i], nodeA._point): @@ -330,7 +333,7 @@ class Graph: self.remove_edge(edge) return [edgeA, edgeB] - def _sanity_check(self, node_is_2=False): + def _sanity_check(self, title, node_is_2=False): """ For debugging purposes: assert that edges and nodes are connected to each other correctly and there are no orphaned @@ -370,12 +373,7 @@ class Graph: from mpl_toolkits.basemap import Basemap from matplotlib import pyplot as plt fig = plt.figure() - m = Basemap(projection="ortho", - lon_0=lon_0, - lat_0=lat_0) - m.drawparallels(np.arange(-90., 90., 20.)) - m.drawmeridians(np.arange(0., 420., 20.)) - m.drawmapboundary(fill_color='white') + m = Basemap() counts = {} for node in self._nodes: @@ -383,11 +381,21 @@ class Graph: counts.setdefault(count, []) counts[count].append(list(node._point)) + minx = np.inf + miny = np.inf + maxx = -np.inf + maxy = -np.inf for k, v in counts.items(): v = np.array(v) ra, dec = vector.vector_to_radec(v[:, 0], v[:, 1], v[:, 2]) x, y = m(ra, dec) m.plot(x, y, 'o', label=str(k)) + for x0 in x: + minx = min(x0, minx) + maxx = max(x0, maxx) + for y0 in y: + miny = min(y0, miny) + maxy = max(y0, maxy) for edge in list(self._edges): A, B = [x._point for x in edge._nodes] @@ -395,8 +403,11 @@ class Graph: r1, d1 = vector.vector_to_radec(B[0], B[1], B[2]) m.drawgreatcircle(r0, d0, r1, d1, color='blue') + plt.xlim(minx, maxx) + plt.ylim(miny, maxy) if title: - plt.title(title) + plt.title("%s, %d v, %d e" % ( + title, len(self._nodes), len(self._edges))) plt.legend() plt.show() @@ -414,13 +425,13 @@ class Graph: will be included in the output. """ self._remove_cut_lines() - self._sanity_check() + self._sanity_check("union - remove cut lines") self._find_all_intersections() - self._sanity_check() - self._remove_interior_nodes() - self._sanity_check() + self._sanity_check("union - find all intersections") + self._remove_interior_edges() + self._sanity_check("union - remove interior edges") self._remove_3ary_edges() - self._sanity_check(True) + self._sanity_check("union - remove 3ary edges", True) return self._trace() def intersection(self): @@ -437,13 +448,13 @@ class Graph: will be included in the output. """ self._remove_cut_lines() - self._sanity_check() + self._sanity_check("intersection - remove cut lines") self._find_all_intersections() - self._sanity_check() + self._sanity_check("intersection - find all intersections") self._remove_exterior_edges() - self._sanity_check() + self._sanity_check("intersection - remove exterior edges") self._remove_3ary_edges(large_first=True) - self._sanity_check(True) + self._sanity_check("intersection - remove 3ary edges", True) return self._trace() def _remove_cut_lines(self): @@ -595,29 +606,24 @@ class Graph: (new_ends, ends[:j], ends[j+1:])) break - def _remove_interior_nodes(self): + def _remove_interior_edges(self): """ Removes any nodes that are contained inside other polygons. What's left is the (possibly disjunct) outline. """ polygons = self._source_polygons - # node._count is the number of polygons that the node is - # inside, other than the polygon that it came from - for node in self._nodes: - node._count = 0 + for edge in self._edges: + edge._count = 0 for polygon in polygons: - if polygon in node._source_polygons: - continue - if polygon.contains_point(node._point): - node._count += 1 + if (not polygon in edge._source_polygons and + polygon.intersects_arc( + edge._nodes[0]._point, edge._nodes[1]._point)): + edge._count += 1 - for node in list(self._nodes): - # Nodes with a count >= 1 are completely contained in the - # outer polygon, so they should be removed. What's left - # are only outer nodes. - if node._count >= 1: - self.remove_node(node) + for edge in list(self._edges): + if edge._count >= 1: + self.remove_edge(edge) def _remove_exterior_edges(self): """ @@ -673,7 +679,13 @@ class Graph: seen_nodes = set() while len(edges): polygon = [] - edge = edges.pop() + # Carefully pick out an "original" edge first. Synthetic + # edges may not be pointing in the right direction to + # properly calculate the area. + for edge in edges: + if len(edge._source_polygons) == 1: + break + edges.remove(edge) start_node = node = edge._nodes[0] while True: # TODO: Do we need this if clause any more? diff --git a/lib/great_circle_arc.py b/lib/great_circle_arc.py index 682f9d9..f70d16e 100644 --- a/lib/great_circle_arc.py +++ b/lib/great_circle_arc.py @@ -218,6 +218,7 @@ def length(A, B, degrees=True): B = np.asanyarray(B) dot = inner1d(A, B) + dot = np.clip(dot, -1.0, 1.0) with np.errstate(invalid='ignore'): result = np.arccos(dot) @@ -350,7 +351,10 @@ def interpolate(A, B, steps=50): t = np.linspace(0.0, 1.0, steps, endpoint=True).reshape((steps, 1)) omega = length(A, B, degrees=False) - sin_omega = np.sin(omega) - offsets = np.sin(t * omega) / sin_omega + if omega == 0.0: + offsets = t + else: + sin_omega = np.sin(omega) + offsets = np.sin(t * omega) / sin_omega return offsets[::-1] * A + offsets * B diff --git a/lib/polygon.py b/lib/polygon.py index 6241416..b7ea5d1 100644 --- a/lib/polygon.py +++ b/lib/polygon.py @@ -204,6 +204,7 @@ class SphericalPolygon(object): # same. 2π should equal 0, but with rounding error that isn't # always the case. C[-1] = 0 + C = C[::-1] X, Y, Z = vector.rotate_around(x, y, z, u, v, w, C, degrees=False) return cls(np.dstack((X, Y, Z))[0], (u, v, w)) @@ -382,29 +383,50 @@ class SphericalPolygon(object): that are larger than half of the sphere. Therefore, the area will always be less than 2π. - The area is computed based on the sum of the radian angles: + The area is computed by transforming the polygon to two + dimensions using the `Lambert azimuthal equal-area projection + <http://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_ .. math:: - S = \sum^n_{i=0} \theta_i - (n - 2) \pi + X = \sqrt{\frac{2}{1-z}}x + + .. math:: + + Y = \sqrt{\frac{2}{1-z}}y + + The individual great arc circle segments are interpolated + before doing the transformation so that the curves are not + straightened in the process. + + It then uses a standard 2D algorithm to compute the area. + + .. math:: + + A = \left| \sum^n_{i=0} X_i Y_{i+1} - X_{i+1}Y_i \right| """ if len(self._points) < 3: return np.array(0.0) - A = self._points[:-1] - B = np.roll(A, 1, 0) - C = np.roll(B, 1, 0) + # Rotate polygon so that center of polygon is at north pole + centroid = np.mean(self._points, axis=0) + points = self._points - (centroid + np.array([0, 0, 1])) - angles = great_circle_arc.angle(A, B, C, degrees=False) - - midpoints = great_circle_arc.midpoint(A, C) - interior = [self.contains_point(x) for x in midpoints] - angles = np.where(interior, angles, 2.*np.pi - angles) + 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.extend(x) + Y.extend(y) - sum = np.sum(angles) - area = sum - float(len(angles) - 2) * np.pi + X = np.array(X) + Y = np.array(Y) - return area + return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) / 2.0) def union(self, other): """ @@ -445,7 +467,7 @@ class SphericalPolygon(object): return self.__class__(polygon, self._inside) @classmethod - def multi_union(cls, polygons, method='parallel'): + def multi_union(cls, polygons): """ Return a new `SphericalPolygon` that is the union of all of the polygons in *polygons*. @@ -454,21 +476,6 @@ class SphericalPolygon(object): ---------- polygons : sequence of `SphericalPolygon` - method : 'parallel' or 'serial', optional - Specifies the method that is used to perform the unions: - - - 'parallel' (default): A graph is built using all of - the polygons, and the union operation is computed on - the entire thing globally. - - - 'serial': The polygon is built in steps by adding one - polygon at a time and unioning at each step. - - This option is provided because one may be faster than the - other depending on context, but it primarily exposed for - testing reasons. Both modes should theoretically provide - equivalent results. - Returns ------- polygon : `SphericalPolygon` object @@ -483,17 +490,9 @@ class SphericalPolygon(object): import graph - if method.lower() == 'parallel': - g = graph.Graph(polygons) - polygon = g.union() - return cls(polygon, polygons[0]._inside) - elif method.lower() == 'serial': - result = copy(polygons[0]) - for polygon in polygons[1:]: - result = result.union(polygon) - return result - else: - raise ValueError("method must be 'parallel' or 'serial'") + g = graph.Graph(polygons) + polygon = g.union() + return cls(polygon, polygons[0]._inside) @staticmethod def _find_new_inside(points, polygons): @@ -670,6 +669,8 @@ class SphericalPolygon(object): for A, B in zip(points[0:-1], points[1:]): length = great_circle_arc.length(A, B, degrees=True) + if not np.isfinite(length): + length = 2 interpolated = great_circle_arc.interpolate(A, B, length * 4) ra, dec = vector.vector_to_radec( interpolated[:, 0], interpolated[:, 1], interpolated[:, 2], diff --git a/lib/test/test_intersection.py b/lib/test/test_intersection.py index adb79f0..12b3ce7 100644 --- a/lib/test/test_intersection.py +++ b/lib/test/test_intersection.py @@ -70,7 +70,7 @@ class intersection_test: 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]) + assert_array_almost_equal(areas, areas[0], decimal=3) return run diff --git a/lib/test/test_union.py b/lib/test/test_union.py index c963b76..d29eb25 100644 --- a/lib/test/test_union.py +++ b/lib/test/test_union.py @@ -36,36 +36,35 @@ class union_test: step_size = int(max(float(num_permutations) / 20.0, 1.0)) if GRAPH_MODE: print("%d permutations" % num_permutations) - for method in ('parallel', 'serial'): - for i, permutation in enumerate( - itertools.islice( - itertools.permutations(polys), - None, None, step_size)): - filename = '%s_%s_union_%04d.svg' % ( - func.__name__, method, i) - print(filename) - - union = polygon.SphericalPolygon.multi_union( - permutation, method=method) - 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() - m = Basemap(projection=self._proj, - lon_0=self._lon_0, - lat_0=self._lat_0) - m.drawmapboundary(fill_color='white') - m.drawparallels(np.arange(-90., 90., 20.)) - m.drawmeridians(np.arange(0., 420., 20.)) - - union.draw(m, color='red', linewidth=3) - for poly in permutation: - poly.draw(m, color='blue', alpha=0.5) - plt.savefig(filename) - fig.clear() + for i, permutation in enumerate( + itertools.islice( + itertools.permutations(polys), + None, None, step_size)): + filename = '%s_union_%04d.svg' % ( + func.__name__, i) + print(filename) + + 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() + m = Basemap(projection=self._proj, + lon_0=self._lon_0, + lat_0=self._lat_0) + m.drawmapboundary(fill_color='white') + m.drawparallels(np.arange(-90., 90., 20.)) + m.drawmeridians(np.arange(0., 420., 20.)) + + union.draw(m, color='red', linewidth=3) + for poly in permutation: + poly.draw(m, color='blue', alpha=0.5) + plt.savefig(filename) + fig.clear() lengths = np.array([len(x._points) for x in unions]) assert np.all(lengths == [lengths[0]]) @@ -158,6 +157,28 @@ def test6(): null_union = chipA1.union(chipA2) +@union_test(0, 90) +def test7(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + + B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz')) + + wcs = pywcs.WCS(B[1].header, fobj=B) + chipB1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[4].header, fobj=B) + chipB2 = polygon.SphericalPolygon.from_wcs(wcs) + + return [chipA1, chipA2, chipB1, chipB2] + + if __name__ == '__main__': if '--profile' not in sys.argv: GRAPH_MODE = True diff --git a/lib/vector.py b/lib/vector.py index 7323c92..5d26161 100644 --- a/lib/vector.py +++ b/lib/vector.py @@ -198,9 +198,6 @@ def rotate_around(x, y, z, u, v, w, theta, degrees=True): if degrees: theta = np.deg2rad(theta) - u2 = u**2 - v2 = v**2 - w2 = w**2 costheta = np.cos(theta) sintheta = np.sin(theta) icostheta = 1.0 - costheta @@ -211,3 +208,34 @@ def rotate_around(x, y, z, u, v, w, theta, degrees=True): Z = (-w*det)*icostheta + z*costheta + (-v*x + u*y)*sintheta return X, Y, Z + + +def equal_area_proj(x, y, z): + """ + Transform the coordinates to a 2-dimensional plane using the + Lambert azimuthal equal-area projection. + + Parameters + ---------- + x, y, z : scalars or 1-D arrays + The input vectors + + Returns + ------- + X, Y : tuple of scalars or arrays of the same length + + Notes + ----- + + .. math:: + + X = \sqrt{\frac{2}{1-z}}x + + .. math:: + + Y = \sqrt{\frac{2}{1-z}}y + """ + scale = np.sqrt(2.0 / (1.0 - z)) + X = scale * x + Y = scale * y + return X, Y |