Changeset 7841 for trunk/anuga_core/source/anuga/geometry

Timestamp:

Jun 15, 2010, 12:06:46 PM (15 years ago)

Author:

James Hudson

Message:

Refactorings to allow tests to pass.

Location:

trunk/anuga_core/source/anuga/geometry

Files:

• __init__.py (modified) (1 diff)
• polygon.py (modified) (23 diffs)
• test_polygon.py (modified) (8 diffs)

trunk/anuga_core/source/anuga/geometry/__init__.py

@@ -2 +2 @@
 """
 
-pass
 
 #Add path of package to PYTHONPATH to allow C-extensions to be loaded

trunk/anuga_core/source/anuga/geometry/polygon.py

@@ -5 +5 @@
 import numpy as num
 
-from math import sqrt
 from anuga.utilities.numerical_tools import ensure_numeric
 from anuga.geospatial_data.geospatial_data import ensure_absolute, \
                                                     Geospatial_data
-from anuga.config import netcdf_float
 import anuga.utilities.log as log
 
@@ -129 +127 @@
     line1 = ensure_numeric(line1, num.float)
 
-    x0 = line0[0,0]; y0 = line0[0,1]
-    x1 = line0[1,0]; y1 = line0[1,1]
-
-    x2 = line1[0,0]; y2 = line1[0,1]
-    x3 = line1[1,0]; y3 = line1[1,1]
+    x0 = line0[0, 0]; y0 = line0[0, 1]
+    x1 = line0[1, 0]; y1 = line0[1, 1]
+
+    x2 = line1[0, 0]; y2 = line1[0, 1]
+    x3 = line1[1, 0]; y3 = line1[1, 1]
 
     denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
@@ -208 +206 @@
     line1 = ensure_numeric(line1, num.float)
 
-    status, value = _intersection(line0[0,0], line0[0,1],
-                                  line0[1,0], line0[1,1],
-                                  line1[0,0], line1[0,1],
-                                  line1[1,0], line1[1,1])
+    status, value = _intersection(line0[0, 0], line0[0, 1],
+                                  line0[1, 0], line0[1, 1],
+                                  line1[0, 0], line1[0, 1],
+                                  line1[1, 0], line1[1, 1])
 
     return status, value
@@ -219 +217 @@
                        rtol=1.0e-12,
                        atol=1.0e-12,                      
-                       check_inputs=True, 
-                       verbose=False):
+                       check_inputs=True):
     """Determine if one point is inside a triangle
     
@@ -278 +275 @@
     
     # Quickly reject points that are clearly outside
-    if point[0] < min(triangle[:, 0]): return False 
-    if point[0] > max(triangle[:, 0]): return False    
-    if point[1] < min(triangle[:, 1]): return False
-    if point[1] > max(triangle[:, 1]): return False        
+    if point[0] < min(triangle[:, 0]):
+        return False 
+    if point[0] > max(triangle[:, 0]): 
+        return False    
+    if point[1] < min(triangle[:, 1]):
+        return False
+    if point[1] > max(triangle[:, 1]):
+        return False        
 
 
@@ -314 +315 @@
         # Check if point lies on one of the edges
         
-        for X, Y in [[A,B], [B,C], [C,A]]:
+        for X, Y in [[A, B], [B, C], [C, A]]:
             res = _point_on_line(point[0], point[1],
                                  X[0], X[1],
@@ -341 +342 @@
     
     def segments_joined(seg0, seg1):
+        """ See if there are identical segments in the 2 lists. """
         for i in seg0:
             for j in seg1:    
@@ -381 +383 @@
     return False
     
-    
-def is_inside_polygon_quick(point, polygon, closed=True, verbose=False):
-    """Determine if one point is inside a polygon
-    Both point and polygon are assumed to be numeric arrays or lists and
-    no georeferencing etc or other checks will take place.
-    
-    As such it is faster than is_inside_polygon 
-    """
-
-    # FIXME(Ole): This function isn't being used
-    polygon = ensure_numeric(polygon, num.float)
-    points = ensure_numeric(point, num.float) # Convert point to array of points
-    points = num.ascontiguousarray(points[num.newaxis, :])
-    msg = ('is_inside_polygon() must be invoked with one point only.\n'
-           'I got %s and converted to %s' % (str(point), str(points.shape)))
-    assert points.shape[0] == 1 and points.shape[1] == 2, msg
-    
-    indices = num.zeros(1, num.int)
-
-    count = _separate_points_by_polygon(points, polygon, indices,
-                                        int(closed), int(verbose))
-
-    return count > 0
-
 
 def is_inside_polygon(point, polygon, closed=True, verbose=False):
@@ -420 +398 @@
     else:
         msg = 'is_inside_polygon must be invoked with one point only'
-        raise msg
+        raise Exception(msg)
 
 ##
@@ -743 +721 @@
                   figname=None,
                   label=None,
-                  alpha=None,
-                  verbose=False):
+                  alpha=None):
     """ Take list of polygons and plot.
 
@@ -769 +746 @@
     """
 
-    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, \
+    from pylab import ion, hold, plot, savefig, xlabel, \
                       ylabel, title, close, title, fill
 
@@ -793 +770 @@
             alpha = None
         else:
-            if alpha < 0.0:
-                alpha = 0.0
-            if alpha > 1.0:
-                alpha = 1.0
+            alpha = max(0.0, min(1.0, alpha))
 
     n = len(polygons_points)
@@ -818 +792 @@
     for i, item in enumerate(polygons_points):
         x, y = poly_xy(item)
-        if min(x) < minx: minx = min(x)
-        if max(x) > maxx: maxx = max(x)
-        if min(y) < miny: miny = min(y)
-        if max(y) > maxy: maxy = max(y)
-        plot(x,y,colour[i])
+        if min(x) < minx:
+            minx = min(x)
+        if max(x) > maxx:
+            maxx = max(x)
+        if min(y) < miny:
+            miny = min(y)
+        if max(y) > maxy:
+            maxy = max(y)
+        plot(x, y, colour[i])
         if alpha:
             fill(x, y, colour[i], alpha=alpha)
@@ -840 +818 @@
 
 
-def poly_xy(polygon, verbose=False):
+def poly_xy(polygon):
     """ this is used within plot_polygons so need to duplicate
         the first point so can have closed polygon in plot
@@ -858 +836 @@
         raise Exception, msg
 
-    x = polygon[:,0]
-    y = polygon[:,1]
-    x = num.concatenate((x, [polygon[0,0]]), axis = 0)
-    y = num.concatenate((y, [polygon[0,1]]), axis = 0)
+    x = polygon[:, 0]
+    y = polygon[:, 1]
+    x = num.concatenate((x, [polygon[0, 0]]), axis = 0)
+    y = num.concatenate((y, [polygon[0, 1]]), axis = 0)
 
     return x, y
@@ -960 +938 @@
         assert y.shape[0] == N
 
-        points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis],
-                                                        y[:,num.newaxis]),
-                                                       axis=1 ))
+        points = num.ascontiguousarray(num.concatenate((x[:, num.newaxis],
+                                                        y[:, num.newaxis]),
+                                                       axis = 1 ))
 
         if callable(self.default):
@@ -1042 +1020 @@
         fid.write('%f, %f\n' % point)
     fid.close()
-
-##
-# @brief Unimplemented.
-def read_tagged_polygons(filename):
-    """
-    """
-    pass
 
 
@@ -1097 +1068 @@
             if exclude is not None:
                 for ex_poly in exclude:
-                    if is_inside_polygon([x,y], ex_poly):
+                    if is_inside_polygon([x, y], ex_poly):
                         append = False
 
         if append is True:
-            points.append([x,y])
+            points.append([x, y])
 
     return points
@@ -1122 +1093 @@
     import exceptions
 
-    class Found(exceptions.Exception): pass
+    class Found(exceptions.Exception):
+        pass
 
     polygon = ensure_numeric(polygon)
@@ -1232 +1204 @@
     # Find outer extent of polygon
     num_polygon = ensure_numeric(polygon)
-    max_x = max(num_polygon[:,0])
-    max_y = max(num_polygon[:,1])
-    min_x = min(num_polygon[:,0])
-    min_y = min(num_polygon[:,1])
+    max_x = max(num_polygon[:, 0])
+    max_y = max(num_polygon[:, 1])
+    min_x = min(num_polygon[:, 0])
+    min_y = min(num_polygon[:, 1])
 
     # Keep only some points making sure extrema are kept
@@ -1257 +1229 @@
                          interpolation_points=None,
                          rtol=1.0e-6,
-                         atol=1.0e-8,
-                         verbose=False):
+                         atol=1.0e-8):
     """Interpolate linearly between values data on polyline nodes
     of a polyline to list of interpolation points.
@@ -1289 +1260 @@
     gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.int)
 
-    n = polyline_nodes.shape[0]    # Number of nodes in polyline
+    num_nodes = polyline_nodes.shape[0]    # Number of nodes in polyline
 
     # Input sanity check
-    msg = 'interpolation_points are not given (interpolate.py)'
-    assert interpolation_points is not None, msg
-
-    msg = 'function value must be specified at every interpolation node'
-    assert data.shape[0] == polyline_nodes.shape[0], msg
-
-    msg = 'Must define function value at one or more nodes'
-    assert data.shape[0] > 0, msg
-
-    if n == 1:
-        msg = 'Polyline contained only one point. I need more. ' + str(data)
-        raise Exception, msg
-    elif n > 1:
+    assert_msg = 'interpolation_points are not given (interpolate.py)'
+    assert interpolation_points is not None, assert_msg
+
+    assert_msg = 'function value must be specified at every interpolation node'
+    assert data.shape[0] == polyline_nodes.shape[0], assert_msg
+
+    assert_msg = 'Must define function value at one or more nodes'
+    assert data.shape[0] > 0, assert_msg
+
+    if num_nodes == 1:
+        assert_msg = 'Polyline contained only one point. I need more. '
+        assert_msg += str(data)
+        raise Exception, assert_msg
+    elif num_nodes > 1:
         _interpolate_polyline(data,
                               polyline_nodes,
@@ -1346 +1318 @@
 
 else:
-    msg = 'C implementations could not be accessed by %s.\n ' %__file__
-    msg += 'Make sure compile_all.py has been run as described in '
-    msg += 'the ANUGA installation guide.'
-    raise Exception, msg
+    error_msg = 'C implementations could not be accessed by %s.\n ' %__file__
+    error_msg += 'Make sure compile_all.py has been run as described in '
+    error_msg += 'the ANUGA installation guide.'
+    raise Exception(error_msg)
 
 

trunk/anuga_core/source/anuga/geometry/test_polygon.py

@@ -178 +178 @@
 
     def test_inside_polygon_main(self):
-        """test_is_inside_polygon_quick
+        """test_is_inside_polygon
         
         Test fast version of of is_inside_polygon
@@ -186 +186 @@
         polygon = [[0,0], [1,0], [1,1], [0,1]]
 
-        assert is_inside_polygon_quick( (0.5, 0.5), polygon )
-        assert not is_inside_polygon_quick( (0.5, 1.5), polygon )
-        assert not is_inside_polygon_quick( (0.5, -0.5), polygon )
-        assert not is_inside_polygon_quick( (-0.5, 0.5), polygon )
-        assert not is_inside_polygon_quick( (1.5, 0.5), polygon )
+        assert is_inside_polygon( (0.5, 0.5), polygon )
+        assert not is_inside_polygon( (0.5, 1.5), polygon )
+        assert not is_inside_polygon( (0.5, -0.5), polygon )
+        assert not is_inside_polygon( (-0.5, 0.5), polygon )
+        assert not is_inside_polygon( (1.5, 0.5), polygon )
 
         # Try point on borders
-        assert is_inside_polygon_quick( (1., 0.5), polygon, closed=True)
-        assert is_inside_polygon_quick( (0.5, 1), polygon, closed=True)
-        assert is_inside_polygon_quick( (0., 0.5), polygon, closed=True)
-        assert is_inside_polygon_quick( (0.5, 0.), polygon, closed=True)
-
-        assert not is_inside_polygon_quick( (0.5, 1), polygon, closed=False)
-        assert not is_inside_polygon_quick( (0., 0.5), polygon, closed=False)
-        assert not is_inside_polygon_quick( (0.5, 0.), polygon, closed=False)
-        assert not is_inside_polygon_quick( (1., 0.5), polygon, closed=False)
+        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
+        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
+        assert is_inside_polygon( (0., 0.5), polygon, closed=True)
+        assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
+
+        assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
+        assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
+        assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
+        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
 
 
@@ -229 +229 @@
         assert not is_inside_polygon( (0.5, -0.5), polygon )
 
-        assert is_inside_polygon_quick( (0.5, 0.5), polygon )
-        assert is_inside_polygon_quick( (1, -0.5), polygon )
-        assert is_inside_polygon_quick( (1.5, 0), polygon )
-
-        assert not is_inside_polygon_quick( (0.5, 1.5), polygon )
-        assert not is_inside_polygon_quick( (0.5, -0.5), polygon )
+        assert is_inside_polygon( (0.5, 0.5), polygon )
+        assert is_inside_polygon( (1, -0.5), polygon )
+        assert is_inside_polygon( (1.5, 0), polygon )
+
+        assert not is_inside_polygon( (0.5, 1.5), polygon )
+        assert not is_inside_polygon( (0.5, -0.5), polygon )
 
         # Very convoluted polygon
@@ -449 +449 @@
             assert is_inside_polygon(point, polygon)
 
-            assert is_inside_polygon_quick(point, polygon)
+            assert is_inside_polygon(point, polygon)
 
 
@@ -458 +458 @@
         for point in points:
             assert is_inside_polygon(point, polygon)
-            assert is_inside_polygon_quick(point, polygon)
+            assert is_inside_polygon(point, polygon)
 
 
@@ -1598 +1598 @@
         assert y[4] == 6
 
-    # Disabled
-    def xtest_plot_polygons(self):
+
+    def test_plot_polygons(self):
         import os
 
@@ -1605 +1605 @@
         polygon1 = [[0,0], [1,0], [1,1], [0,1]]
         polygon2 = [[1,1], [2,1], [3,2], [2,2]]
-        v = plot_polygons([polygon1, polygon2], 'test1')
+        v = plot_polygons([polygon1, polygon2], figname='test1')
         assert len(v) == 4
         assert v[0] == 0
@@ -1614 +1614 @@
         # Another case
         polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
-        v = plot_polygons([polygon2,polygon3], 'test2')
+        v = plot_polygons([polygon2,polygon3], figname='test2')
         assert len(v) == 4
         assert v[0] == 1