Fix ordering of nodes in the trace routine so that the area calculation of polygons with cutlines is correct.

git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci.sphere/trunk@33584 fe389314-cf27-0410-b35b-8c050e845b92

Former-commit-id: 04c6e2875f229294be3c8592ba73d09457c42984

3 files changed, 136 insertions, 8 deletions

diff --git a/stsci/sphere/graph.py b/stsci/sphere/graph.py
index bae1c60..01c582a 100644
--- a/stsci/sphere/graph.py
+++ b/stsci/sphere/graph.py
@@ -814,6 +814,28 @@ class Graph:
             else:
                 break
 
+    @staticmethod
+    def _fast_area(points):
+        """
+        Calculates an approximate area of a polygon.  Unlike the area
+        calculation in `SphericalPolygon`, this does not interpolate
+        along the arcs, so is not useful for getting an actual area,
+        but only whether the area is positive or negative, telling us
+        whether the points proceed clockwise or counter-clockwise
+        respectively.
+        """
+        # Rotate polygon so that center of polygon is at north pole
+        centroid = np.mean(points[:-1], axis=0)
+        centroid = vector.normalize_vector(centroid)
+        points = points - (centroid + np.array([0, 0, 1]))
+        vector.normalize_vector(points, output=points)
+
+        XY = vector.equal_area_proj(points)
+        X = XY[..., 0]
+        Y = XY[..., 1]
+
+        return np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1])
+
     def _trace(self):
         """
         Given a graph that has had cutlines removed and all
@@ -824,13 +846,7 @@ class Graph:
         edges = set(self._edges)  # copy
         while len(edges):
             polygon = []
-            # 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)
+            edge = edges.pop()
             start_node = node = edge._nodes[0]
             while True:
                 # TODO: Do we need this if clause any more?
@@ -846,7 +862,13 @@ class Graph:
                     polygon.append(node._point)
                     break
 
-            polygons.append(np.asarray(polygon))
+            polygon = np.asarray(polygon)
+            area = self._fast_area(polygon)
+            if area < 0:
+                # Reverse the points
+                polygon = polygon[::-1, ...]
+
+            polygons.append(polygon)
 
         if len(polygons) == 1:
             return polygons[0]
diff --git a/stsci/sphere/test/test.py b/stsci/sphere/test/test.py
index 5ab7bbb..567ec1d 100644
--- a/stsci/sphere/test/test.py
+++ b/stsci/sphere/test/test.py
@@ -276,3 +276,34 @@ def test_area():
         points = np.dstack((x, y, z))[0]
         poly = polygon.SphericalPolygon(points)
         calc_area = poly.area()
+
+
+def test_fast_area():
+    a = np.array(  # Clockwise
+        [[ 0.35327617,  0.6351561 , -0.6868571 ],
+         [ 0.35295533,  0.63510299, -0.68707112],
+         [ 0.35298984,  0.63505081, -0.68710162],
+         [ 0.35331262,  0.63510039, -0.68688987],
+         [ 0.35327617,  0.6351561 , -0.6868571 ]])
+
+    b = np.array([  # Clockwise
+        [ 0.35331737,  0.6351013 , -0.68688658],
+        [ 0.3536442 ,  0.63515101, -0.68667239],
+        [ 0.35360581,  0.63521041, -0.68663722],
+        [ 0.35328338,  0.63515742, -0.68685217],
+        [ 0.35328614,  0.63515318, -0.68685467],
+        [ 0.35328374,  0.63515279, -0.68685627],
+        [ 0.35331737,  0.6351013 , -0.68688658]])
+
+    c = np.array([  # Counterclockwise
+        [ 0.35331737,  0.6351013 , -0.68688658],
+        [ 0.35328374,  0.63515279, -0.68685627],
+        [ 0.35328614,  0.63515318, -0.68685467],
+        [ 0.35328338,  0.63515742, -0.68685217],
+        [ 0.35360581,  0.63521041, -0.68663722],
+        [ 0.3536442 ,  0.63515101, -0.68667239],
+        [ 0.35331737,  0.6351013 , -0.68688658]])
+
+    assert graph.Graph._fast_area(a) > 0
+    assert graph.Graph._fast_area(b) > 0
+    assert graph.Graph._fast_area(c) < 0
diff --git a/stsci/sphere/test/test_intersection.py b/stsci/sphere/test/test_intersection.py
index 8cf10b2..92a3413 100644
--- a/stsci/sphere/test/test_intersection.py
+++ b/stsci/sphere/test/test_intersection.py
@@ -190,6 +190,81 @@ def test_difficult_intersections():
         yield test_intersection, (polyA, polyB)
 
 
+def test_ordering():
+    nrepeat = 10
+
+    A = polygon.SphericalPolygon(
+        [[3.532808036921135653e-01, 6.351523005458726834e-01, -6.868582305351954576e-01],
+         [3.532781068942476010e-01, 6.351564219435104075e-01, -6.868558064493115456e-01],
+         [3.529538811375814156e-01, 6.351027504797477352e-01, -6.870720880104047579e-01],
+         [3.533428330964511477e-01, 6.345142927049303161e-01, -6.874157800432978416e-01],
+         [3.533486351814376647e-01, 6.345151843837375516e-01, -6.874119745843003670e-01],
+         [3.533513056857608414e-01, 6.345111416839894769e-01, -6.874143334620310686e-01],
+         [3.536740696809928530e-01, 6.345607036635456666e-01, -6.872025653337667794e-01],
+         [3.536713200704008631e-01, 6.345649108795897719e-01, -6.872000954889618818e-01],
+         [3.536761865498951884e-01, 6.345656515431040701e-01, -6.871969069700470945e-01],
+         [3.536788213460497765e-01, 6.345616140129455296e-01, -6.871992792142280759e-01],
+         [3.540056257094351122e-01, 6.346113105009757449e-01, -6.869850810245486938e-01],
+         [3.536200722272911379e-01, 6.352081961257413090e-01, -6.866319189293832448e-01],
+         [3.536142814048366390e-01, 6.352072452054380314e-01, -6.866357809093986964e-01],
+         [3.536116196666648781e-01, 6.352113634102898310e-01, -6.866333419163089813e-01],
+         [3.532833767830895755e-01, 6.351574192193063517e-01, -6.868521736876195272e-01],
+         [3.532861440234288386e-01, 6.351531838825796861e-01, -6.868546669018701367e-01],
+         [3.532808036921135653e-01, 6.351523005458726834e-01, -6.868582305351954576e-01]],
+        [3.536414047913637448e-01, 6.348851549491377755e-01, -6.869196436573932196e-01])
+
+    B = polygon.SphericalPolygon(
+        [[3.529249199274748783e-01, 6.356925960489819838e-01, -6.865412764158403958e-01],
+         [3.533126219535084322e-01, 6.351003877952851040e-01, -6.868898664200949744e-01],
+         [3.533173735956686712e-01, 6.351012981906917210e-01, -6.868865805589428053e-01],
+         [3.529301898742857047e-01, 6.356935934402119237e-01, -6.865376437853726310e-01],
+         [3.532584388080926563e-01, 6.357475490961038700e-01, -6.863188247667159070e-01],
+         [3.536441982306618437e-01, 6.351510082118909661e-01, -6.866723948326530769e-01],
+         [3.533173735956686712e-01, 6.351012981906917210e-01, -6.868865805589428053e-01],
+         [3.533126219535084322e-01, 6.351003877952851040e-01, -6.868898664200949744e-01],
+         [3.529898380712340189e-01, 6.350508125724935171e-01, -6.871016225198859351e-01],
+         [3.526006883384300017e-01, 6.356389133339014341e-01, -6.867575456003104373e-01],
+         [3.529249199274748783e-01, 6.356925960489819838e-01, -6.865412764158403958e-01]],
+        [3.532883212044564125e-01, 6.354215160430938258e-01, -6.866053153377369433e-01])
+
+    areas = []
+    for i in xrange(nrepeat):
+        C = A.intersection(B)
+        areas.append(C.area())
+    areas = np.array(areas)
+    assert_array_almost_equal(areas[:-1], areas[1:])
+
+    def roll_polygon(P, i):
+        points = P.points
+        points = np.roll(points[:-1], i, 0)
+        points = np.append(points, [points[0]], 0)
+        return polygon.SphericalPolygon(points, P.inside)
+
+    Aareas = []
+    Bareas = []
+    Careas = []
+    for i in xrange(nrepeat):
+        AS = roll_polygon(A, i)
+        BS = roll_polygon(B, i)
+
+        C = AS.intersection(BS)
+
+        Aareas.append(A.area())
+        Bareas.append(B.area())
+        Careas.append(C.area())
+
+        for j in xrange(nrepeat):
+            CS = roll_polygon(C, j)
+            Careas.append(CS.area())
+
+    Aareas = np.array(Aareas)
+    Bareas = np.array(Bareas)
+    Careas = np.array(Careas)
+    assert_array_almost_equal(Aareas[:-1], Aareas[1:])
+    assert_array_almost_equal(Bareas[:-1], Bareas[1:])
+    assert_array_almost_equal(Careas[:-1], Careas[1:])
+
+
 if __name__ == '__main__':
     if '--profile' not in sys.argv:
         GRAPH_MODE = True