Changeset 1093

Timestamp:

Mar 16, 2005, 11:22:38 PM (20 years ago)

Author:

ole

Message:

Rewrote inside and outside polygon in terms of new separate_points_by_polygon in order to make outside_polygon as fast as inside_polygon.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• test_util.py (modified) (2 diffs)
• util.py (modified) (6 diffs)
• util_ext.c (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -936 +936 @@
         res = outside_polygon( points, polygon )
 
-        assert allclose( res, [3,4] )
+        assert allclose( res, [3, 4] )
 
 
@@ -955 +955 @@
         #evaluate to True as the point 0.5, 1.0 is outside the unit square
 
+    def test_separate_points_by_polygon(self):
+        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
+
+        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
+        assert allclose( indices, [0,2,1] )
+        assert count == 2
+        
+        #One more test of vector formulation returning indices
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+        res, count = separate_points_by_polygon( points, polygon )
+
+        assert allclose( res, [0,1,2,4,3] )
+        assert count == 3
+
+
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+        res, count = separate_points_by_polygon( points, polygon )
+
+        assert allclose( res, [1,2,3,5,4,0] )        
+        assert count == 3
+        
+
     def test_populate_polygon(self):
 

inundation/ga/storm_surge/pyvolution/util.py

@@ -1 @@
-"""This modul contains various auxiliary function used by pyvolution.
-
-It is also a clearing house for function tat may later earn a module
+"""This module contains various auxiliary function used by pyvolution.
+
+It is also a clearing house for functions that may later earn a module
 of their own. 
 """
@@ -20 +20 @@
         theta = acos(v1)
     except ValueError, e:
-        print 'WARNING Angle acos(%s) failed: %s' %(str(v1), e)
+        print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
 
         #FIXME, hack to avoid acos(1.0) Value error
@@ -30 +30 @@
         elif (v1+s >  -1.0) and (v1-s < -1.0):
             theta = 3.1415926535897931
-        print "WARNING, setting theta to ", theta   
+        print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
+              %(v1, theta)   
+              
     if v2 < 0:
         #Quadrant 3 or 4
@@ -912 +914 @@
 #POLYGON STUFF
 #
-#FIXME: ALl these should be put into new module polygon.py
-
-def inside_polygon(point, polygon, closed = True, verbose = False):
+#FIXME: All these should be put into new module polygon.py
+
+
+             
+#Redefine these to use separate_by_polygon which will put all inside indices 
+#in the first part of the indices array and outside indices in the last
+
+def inside_polygon(points, polygon, closed = True, verbose = False):
+    """See separate_points_by_polygon for documentation
+    """
+    
+    from Numeric import array, Float, reshape    
+    
+    if verbose: print 'Checking input to inside_polygon'    
+    try:
+        points = array(points).astype(Float)       
+    except:
+        msg = 'Points could not be converted to Numeric array'
+        raise msg        
+        
+    try:
+        polygon = array(polygon).astype(Float)             
+    except:
+        msg = 'Polygon could not be converted to Numeric array'
+        raise msg               
+    
+    
+    
+    if len(points.shape) == 1:
+        one_point = True
+        points = reshape(points, (1,2))
+    else:       
+        one_point = False
+
+    indices, count = separate_points_by_polygon(points, polygon, 
+                                                closed, verbose)     
+    
+    if one_point:
+        return count == 1
+    else:
+        return indices[:count]    
+        
+def outside_polygon(points, polygon, closed = True, verbose = False):
+    """See separate_points_by_polygon for documentation
+    """
+    
+    from Numeric import array, Float, reshape    
+    
+    if verbose: print 'Checking input to outside_polygon'    
+    try:
+        points = array(points).astype(Float)       
+    except:
+        msg = 'Points could not be converted to Numeric array'
+        raise msg        
+        
+    try:
+        polygon = array(polygon).astype(Float)             
+    except:
+        msg = 'Polygon could not be converted to Numeric array'
+        raise msg               
+    
+    
+    
+    if len(points.shape) == 1:
+        one_point = True
+        points = reshape(points, (1,2))
+    else:       
+        one_point = False
+
+    indices, count = separate_points_by_polygon(points, polygon, 
+                                                closed, verbose)     
+    
+    if one_point:
+        return count != 1
+    else:
+        return indices[count:][::-1]  #return reversed    
+        
+
+def separate_points_by_polygon(points, polygon, 
+                               closed = True, verbose = False):
+    """Determine whether points are inside or outside a polygon
+    
+    Input:
+       points - Tuple of (x, y) coordinates, or list of tuples
+       polygon - list of vertices of polygon
+       closed - (optional) determine whether points on boundary should be 
+       regarded as belonging to the polygon (closed = True) 
+       or not (closed = False) 
+    
+    Outputs:
+       indices: array of same length as points with indices of points falling 
+       inside the polygon listed from the beginning and indices of points 
+       falling outside listed from the end.
+       
+       count: count of points falling inside the polygon
+       
+       The indices of points inside are obtained as indices[:count]
+       The indices of points outside are obtained as indices[count:]       
+       
+       
+    Examples:
+       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
+       
+       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
+       will return the indices [0, 2, 1] and count == 2 as only the first 
+       and the last point are inside the unit square
+       
+    Remarks:
+       The vertices may be listed clockwise or counterclockwise and
+       the first point may optionally be repeated.    
+       Polygons do not need to be convex.
+       Polygons can have holes in them and points inside a hole is 
+       regarded as being outside the polygon.          
+               
+    Algorithm is based on work by Darel Finley, 
+    http://www.alienryderflex.com/polygon/ 
+    
+    """    
+
+    from Numeric import array, Float, reshape, Int, zeros
+    
+    
+    #Input checks
+    try:
+        points = array(points).astype(Float)       
+    except:
+        msg = 'Points could not be converted to Numeric array'
+        raise msg        
+        
+    try:
+        polygon = array(polygon).astype(Float)             
+    except:
+        msg = 'Polygon could not be converted to Numeric array'
+        raise msg               
+
+    assert len(polygon.shape) == 2,\
+       'Polygon array must be a 2d array of vertices'
+
+    assert polygon.shape[1] == 2,\
+       'Polygon array must have two columns'
+
+    assert len(points.shape) == 2,\
+       'Points array must be a 2d array'
+       
+    assert points.shape[1] == 2,\
+       'Points array must have two columns'
+       
+    N = polygon.shape[0] #Number of vertices in polygon 
+    M = points.shape[0]  #Number of points
+    
+    px = polygon[:,0]
+    py = polygon[:,1]   
+
+    #Used for an optimisation when points are far away from polygon
+    minpx = min(px); maxpx = max(px)
+    minpy = min(py); maxpy = max(py)    
+
+
+    #Begin main loop
+    indices = zeros(M, Int)
+    
+    inside_index = 0    #Keep track of points inside
+    outside_index = M-1 #Keep track of points outside (starting from end) 
+    
+    for k in range(M):
+
+        if verbose:
+            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
+        
+        #for each point   
+        x = points[k, 0]
+        y = points[k, 1]        
+
+        inside = False
+
+        if not x > maxpx or x < minpx or y > maxpy or y < minpy: 
+            #Check polygon 
+            for i in range(N):
+                j = (i+1)%N
+                
+                #Check for case where point is contained in line segment        
+                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
+                    if closed:
+                        inside = True
+                    else:        
+                        inside = False
+                    break
+                else:    
+                    #Check if truly inside polygon                  
+                    if py[i] < y and py[j] >= y or\
+                       py[j] < y and py[i] >= y:
+                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
+                            inside = not inside
+        
+        if inside: 
+            indices[inside_index] = k
+            inside_index += 1
+        else:
+            indices[outside_index] = k
+            outside_index -= 1    
+            
+    return indices, inside_index
+
+
+def separate_points_by_polygon_c(points, polygon, 
+                                 closed = True, verbose = False):
+    """Determine whether points are inside or outside a polygon
+
+    C-wrapper
+    """    
+
+    from Numeric import array, Float, reshape, zeros, Int
+    
+
+    if verbose: print 'Checking input to separate_points_by_polygon'    
+    #Input checks
+    try:
+        points = array(points).astype(Float)       
+    except:
+        msg = 'Points could not be converted to Numeric array'
+        raise msg        
+        
+    try:
+        polygon = array(polygon).astype(Float)             
+    except:
+        msg = 'Polygon could not be converted to Numeric array'
+        raise msg               
+
+    assert len(polygon.shape) == 2,\
+       'Polygon array must be a 2d array of vertices'
+
+    assert polygon.shape[1] == 2,\
+       'Polygon array must have two columns'
+
+    assert len(points.shape) == 2,\
+       'Points array must be a 2d array'
+       
+    assert points.shape[1] == 2,\
+       'Points array must have two columns'
+       
+    N = polygon.shape[0] #Number of vertices in polygon 
+    M = points.shape[0]  #Number of points
+
+    from util_ext import separate_points_by_polygon
+
+    if verbose: print 'Allocating array for indices'
+
+    indices = zeros( M, Int )
+
+    if verbose: print 'Calling C-version of inside poly' 
+    count = separate_points_by_polygon(points, polygon, indices,
+                                       int(closed), int(verbose))
+
+    return indices, count 
+
+
+
+class Polygon_function:
+    """Create callable object f: x,y -> z, where a,y,z are vectors and 
+    where f will return different values depending on whether x,y belongs 
+    to specified polygons.
+    
+    To instantiate:
+      
+       Polygon_function(polygons)
+       
+    where polygons is a list of tuples of the form
+    
+      [ (P0, v0), (P1, v1), ...]
+      
+      with Pi being lists of vertices defining polygons and vi either 
+      constants or functions of x,y to be applied to points with the polygon.
+      
+    The function takes an optional argument, default which is the value 
+    (or function) to used for points not belonging to any polygon.
+    For example:
+    
+       Polygon_function(polygons, default = 0.03)       
+      
+    If omitted the default value will be 0.0  
+    
+    Note: If two polygons overlap, the one last in the list takes precedence
+
+    """
+    
+    def __init__(self, regions, default = 0.0):
+        
+        try:
+            len(regions)
+        except:
+            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
+            raise msg
+
+
+        T = regions[0]                                   
+        try:
+            a = len(T)
+        except:    
+            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
+            raise msg   
+                   
+        assert a == 2, 'Must have two component each: %s' %T       
+
+        self.regions = regions             
+        self.default = default
+         
+         
+    def __call__(self, x, y):
+        from util import inside_polygon
+        from Numeric import ones, Float, concatenate, array, reshape, choose
+        
+        x = array(x).astype(Float)
+        y = array(y).astype(Float)      
+        
+        N = len(x)
+        assert len(y) == N
+        
+        points = concatenate( (reshape(x, (N, 1)), 
+                               reshape(y, (N, 1))), axis=1 )
+
+        if callable(self.default):
+            z = self.default(x,y)
+        else:                   
+            z = ones(N, Float) * self.default
+            
+        for polygon, value in self.regions:
+            indices = inside_polygon(points, polygon)
+            
+            #FIXME: This needs to be vectorised 
+            if callable(value):
+                for i in indices:
+                    xx = array([x[i]])
+                    yy = array([y[i]])
+                    z[i] = value(xx, yy)[0]
+            else:    
+                for i in indices:
+                    z[i] = value 
+
+        return z                 
+
+def read_polygon(filename):
+    """Read points assumed to form a polygon
+       There must be exactly two numbers in each line
+    """              
+    
+    #Get polygon
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    polygon = []
+    for line in lines:
+        fields = line.split(',')
+        polygon.append( [float(fields[0]), float(fields[1])] )
+
+    return polygon      
+    
+def populate_polygon(polygon, number_of_points, seed = None):
+    """Populate given polygon with uniformly distributed points.
+
+    Input:
+       polygon - list of vertices of polygon
+       number_of_points - (optional) number of points 
+       seed - seed for random number generator (default=None)
+       
+    Output:
+       points - list of points inside polygon
+       
+    Examples:
+       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
+       will return five randomly selected points inside the unit square 
+    """
+
+    from random import uniform, seed
+
+    seed(seed)
+
+    points = []
+
+    #Find outer extent of polygon
+    max_x = min_x = polygon[0][0]
+    max_y = min_y = polygon[0][1]
+    for point in polygon[1:]:
+        x = point[0] 
+        if x > max_x: max_x = x
+        if x < min_x: min_x = x           
+        y = point[1] 
+        if y > max_y: max_y = y
+        if y < min_y: min_y = y           
+
+
+    while len(points) < number_of_points:     
+        x = uniform(min_x, max_x)
+        y = uniform(min_y, max_y)        
+
+        if inside_polygon( [x,y], polygon ):
+            points.append([x,y])
+
+    return points
+
+####################################################################
+#Python versions of function that are also implemented in util_gateway.c
+#
+
+def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
+    """
+    """
+    
+    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)            
+    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
+    a /= det
+
+    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
+    b /= det            
+
+    return a, b
+
+
+
+##############################################
+#Initialise module
+
+import compile
+if compile.can_use_C_extension('util_ext.c'):
+    from util_ext import gradient, point_on_line
+    separate_points_by_polygon = separate_points_by_polygon_c
+else:
+    gradient = gradient_python
+
+
+if __name__ == "__main__":
+    pass
+
+    
+    
+    
+    
+
+#OBSOLETED STUFF    
+def inside_polygon_old(point, polygon, closed = True, verbose = False):
+    #FIXME Obsoleted
     """Determine whether points are inside or outside a polygon
     
@@ -1036 +1475 @@
              
 
+#def remove_from(A, B):
+#    """Assume that A
+#    """
+#    from Numeric import search_sorted##
+#
+#    ind = search_sorted(A, B)
+
+    
+
+def outside_polygon_old(point, polygon, closed = True, verbose = False):
+    #OBSOLETED
+    """Determine whether points are outside a polygon
+    
+    Input:
+       point - Tuple of (x, y) coordinates, or list of tuples
+       polygon - list of vertices of polygon
+       closed - (optional) determine whether points on boundary should be 
+       regarded as belonging to the polygon (closed = True) 
+       or not (closed = False) 
+    
+    Output:
+       If one point is considered, True or False is returned.
+       If multiple points are passed in, the function returns indices 
+       of those points that fall outside the polygon     
+       
+    Examples:
+       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
+       outside_polygon( [0.5, 0.5], U )
+       will evaluate to False as the point 0.5, 0.5 is inside the unit square
+
+       ouside_polygon( [1.5, 0.5], U )
+       will evaluate to True as the point 1.5, 0.5 is outside the unit square
+       
+       outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
+       will return the indices [1] as only the first and the last point 
+       is inside the unit square
+    """
+
+    #FIXME: This is too slow
+    
+    res  = inside_polygon(point, polygon, closed, verbose)
+
+    if res is True or res is False:
+        return not res
+
+    #Now invert indices
+    from Numeric import arrayrange, compress
+    outside_indices = arrayrange(len(point))
+    for i in res:
+        outside_indices = compress(outside_indices != i, outside_indices)
+    return outside_indices
 
 def inside_polygon_c(point, polygon, closed = True, verbose = False):
+    #FIXME: Obsolete
     """Determine whether points are inside or outside a polygon
 
@@ -1086 +1577 @@
 
     if one_point:
-        return count == 1
+        return count == 1 #Return True if the point was inside
     else:
         if verbose: print 'Got %d points' %count
-        return indices[:count]
-
-
-#def remove_from(A, B):
-#    """Assume that A
-#    """
-#    from Numeric import search_sorted##
-#
-#    ind = search_sorted(A, B)
-
-    
-
-def outside_polygon(point, polygon, closed = True, verbose = False):
-    """Determine whether points are outside a polygon
-    
-    Input:
-       point - Tuple of (x, y) coordinates, or list of tuples
-       polygon - list of vertices of polygon
-       closed - (optional) determine whether points on boundary should be 
-       regarded as belonging to the polygon (closed = True) 
-       or not (closed = False) 
-    
-    Output:
-       If one point is considered, True or False is returned.
-       If multiple points are passed in, the function returns indices 
-       of those points that fall outside the polygon     
-       
-    Examples:
-       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
-       outside_polygon( [0.5, 0.5], U )
-       will evaluate to False as the point 0.5, 0.5 is inside the unit square
-
-       ouside_polygon( [1.5, 0.5], U )
-       will evaluate to True as the point 1.5, 0.5 is outside the unit square
-       
-       outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
-       will return the indices [1] as only the first and the last point 
-       is inside the unit square
-    """
-
-    res  = inside_polygon(point, polygon, closed, verbose)
-
-    if res is True or res is False:
-        return not res
-
-    #Now invert indices
-    from Numeric import arrayrange, compress
-    outside_indices = arrayrange(len(point))
-    for i in res:
-        outside_indices = compress(outside_indices != i, outside_indices)
-    return outside_indices
-
-
-
-class Polygon_function:
-    """Create callable object f: x,y -> z, where a,y,z are vectors and 
-    where f will return different values depending on whether x,y belongs 
-    to specified polygons.
-    
-    To instantiate:
-      
-       Polygon_function(polygons)
-       
-    where polygons is a list of tuples of the form
-    
-      [ (P0, v0), (P1, v1), ...]
-      
-      with Pi being lists of vertices defining polygons and vi either 
-      constants or functions of x,y to be applied to points with the polygon.
-      
-    The function takes an optional argument, default which is the value 
-    (or function) to used for points not belonging to any polygon.
-    For example:
-    
-       Polygon_function(polygons, default = 0.03)       
-      
-    If omitted the default value will be 0.0  
-    
-    Note: If two polygons overlap, the one last in the list takes precedence
-
-    """
-    
-    def __init__(self, regions, default = 0.0):
-        
-        try:
-            len(regions)
-        except:
-            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
-            raise msg
-
-
-        T = regions[0]                                   
-        try:
-            a = len(T)
-        except:    
-            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
-            raise msg   
-                   
-        assert a == 2, 'Must have two component each: %s' %T       
-
-        self.regions = regions             
-        self.default = default
-         
-         
-    def __call__(self, x, y):
-        from util import inside_polygon
-        from Numeric import ones, Float, concatenate, array, reshape, choose
-        
-        x = array(x).astype(Float)
-        y = array(y).astype(Float)      
-        
-        N = len(x)
-        assert len(y) == N
-        
-        points = concatenate( (reshape(x, (N, 1)), 
-                               reshape(y, (N, 1))), axis=1 )
-
-        if callable(self.default):
-            z = self.default(x,y)
-        else:                   
-            z = ones(N, Float) * self.default
-            
-        for polygon, value in self.regions:
-            indices = inside_polygon(points, polygon)
-            
-            #FIXME: This needs to be vectorised 
-            if callable(value):
-                for i in indices:
-                    xx = array([x[i]])
-                    yy = array([y[i]])
-                    z[i] = value(xx, yy)[0]
-            else:    
-                for i in indices:
-                    z[i] = value 
-
-        return z                 
-
-def read_polygon(filename):
-    """Read points assumed to form a polygon
-       There must be exactly two numbers in each line
-    """              
-    
-    #Get polygon
-    fid = open(filename)
-    lines = fid.readlines()
-    fid.close()
-    polygon = []
-    for line in lines:
-        fields = line.split(',')
-        polygon.append( [float(fields[0]), float(fields[1])] )
-
-    return polygon      
-    
-def populate_polygon(polygon, number_of_points, seed = None):
-    """Populate given polygon with uniformly distributed points.
-
-    Input:
-       polygon - list of vertices of polygon
-       number_of_points - (optional) number of points 
-       seed - seed for random number generator (default=None)
-       
-    Output:
-       points - list of points inside polygon
-       
-    Examples:
-       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
-       will return five randomly selected points inside the unit square 
-    """
-
-    from random import uniform, seed
-
-    seed(seed)
-
-    points = []
-
-    #Find outer extent of polygon
-    max_x = min_x = polygon[0][0]
-    max_y = min_y = polygon[0][1]
-    for point in polygon[1:]:
-        x = point[0] 
-        if x > max_x: max_x = x
-        if x < min_x: min_x = x           
-        y = point[1] 
-        if y > max_y: max_y = y
-        if y < min_y: min_y = y           
-
-
-    while len(points) < number_of_points:     
-        x = uniform(min_x, max_x)
-        y = uniform(min_y, max_y)        
-
-        if inside_polygon( [x,y], polygon ):
-            points.append([x,y])
-
-    return points
-
-####################################################################
-#Python versions of function that are also implemented in util_gateway.c
-#
-
-def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
-    """
-    """
-    
-    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)            
-    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
-    a /= det
-
-    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
-    b /= det            
-
-    return a, b
-
-
-
-##############################################
-#Initialise module
-
-import compile
-if compile.can_use_C_extension('util_ext.c'):
-    from util_ext import gradient, point_on_line
-    inside_polygon = inside_polygon_c
-else:
-    gradient = gradient_python
-
-
-if __name__ == "__main__":
-    pass
-
+        
+        return indices[:count]   #Return those indices that were inside
+

inundation/ga/storm_surge/pyvolution/util_ext.c

@@ -62 +62 @@
 
 
+
+int _separate_points_by_polygon(int M,     // Number of points
+                                int N,     // Number of polygon vertices
+                                double* points,
+                                double* polygon,
+                                long* indices,  // M-Array for storage indices
+                                int closed,
+                                int verbose) {
+
+  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
+  int i, j, k, outside_index, inside_index, inside;
+
+  //Find min and max of poly used for optimisation when points
+  //are far away from polygon
+
+  minpx = polygon[0]; maxpx = minpx;
+  minpy = polygon[1]; maxpy = minpy;
+
+  for (i=0; i<N; i++) {
+    px_i = polygon[2*i];
+    py_i = polygon[2*i + 1];
+
+    if (px_i < minpx) minpx = px_i;
+    if (px_i > maxpx) maxpx = px_i;
+    if (py_i < minpy) minpy = py_i;
+    if (py_i > maxpy) maxpy = py_i;
+  }
+
+  //Begin main loop (for each point)
+  inside_index = 0;    //Keep track of points inside
+  outside_index = M-1; //Keep track of points outside (starting from end)   
+  for (k=0; k<M; k++) {
+    if (verbose){
+      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
+    }
+    
+    x = points[2*k];
+    y = points[2*k + 1];
+
+    inside = 0;
+
+    //Optimisation
+    if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) {
+      //Nothing
+    } else {   
+      //Check polygon
+      for (i=0; i<N; i++) {
+        //printf("k,i=%d,%d\n", k, i);
+        j = (i+1)%N;
+
+        px_i = polygon[2*i];
+        py_i = polygon[2*i+1];
+        px_j = polygon[2*j];
+        py_j = polygon[2*j+1];
+
+        //Check for case where point is contained in line segment
+        if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
+          if (closed == 1) {
+            inside = 1;
+          } else {
+            inside = 0;
+          }
+          break;
+        } else {
+          //Check if truly inside polygon
+          if ( ((py_i < y) && (py_j >= y)) ||
+               ((py_j < y) && (py_i >= y)) ) {
+            if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
+              inside = 1-inside;
+          }
+        }
+      }
+    } 
+    if (inside == 1) {
+      indices[inside_index] = k;
+      inside_index += 1;
+    } else {
+      indices[outside_index] = k;
+      outside_index -= 1;    
+    }
+  } // End k
+
+  return inside_index;
+}
+
+
+
+// Gateways to Python
+PyObject *point_on_line(PyObject *self, PyObject *args) {
+  //
+  // point_on_line(x, y, x0, y0, x1, y1)
+  //
+
+  double x, y, x0, y0, x1, y1;
+  int res;
+  PyObject *result;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1))
+    return NULL;
+
+
+  // Call underlying routine
+  res = _point_on_line(x, y, x0, y0, x1, y1);
+
+  // Return values a and b
+  result = Py_BuildValue("i", res);
+  return result;
+}
+
+
+
+PyObject *separate_points_by_polygon(PyObject *self, PyObject *args) {
+  //def separate_points_by_polygon(points, polygon, closed, verbose, one_point):
+  //  """Determine whether points are inside or outside a polygon
+  //
+  //  Input:
+  //     point - Tuple of (x, y) coordinates, or list of tuples
+  //     polygon - list of vertices of polygon
+  //     closed - (optional) determine whether points on boundary should be
+  //     regarded as belonging to the polygon (closed = True)
+  //     or not (closed = False)
+
+  //
+  //  Output:
+  //     indices: array of same length as points with indices of points falling 
+  //     inside the polygon listed from the beginning and indices of points 
+  //     falling outside listed from the end.
+  //     
+  //     count: count of points falling inside the polygon
+  //     
+  //     The indices of points inside are obtained as indices[:count]
+  //     The indices of points outside are obtained as indices[count:]       
+  //
+  //  Examples:
+  //     separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
+  //     will return the indices [0, 2, 1] as only the first and the last point
+  //     is inside the unit square
+  //
+  //  Remarks:
+  //     The vertices may be listed clockwise or counterclockwise and
+  //     the first point may optionally be repeated.
+  //     Polygons do not need to be convex.
+  //     Polygons can have holes in them and points inside a hole is
+  //     regarded as being outside the polygon.
+  //
+  //
+  //  Algorithm is based on work by Darel Finley,
+  //  http://www.alienryderflex.com/polygon/
+  //
+  //
+
+  PyArrayObject
+    *points,
+    *polygon,
+    *indices;
+
+  int closed, verbose; //Flags
+  int count, M, N;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOii",
+                        &points,
+                        &polygon,
+                        &indices,
+                        &closed,
+                        &verbose))
+    return NULL;
+
+  M = points -> dimensions[0];   //Number of points
+  N = polygon -> dimensions[0];  //Number of vertices in polygon
+
+  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
+  
+  //Call underlying routine
+  count = _separate_points_by_polygon(M, N,
+                                      (double*) points -> data,
+                                      (double*) polygon -> data,
+                                      (long*) indices -> data,
+                                      closed, verbose);
+  
+  //NOTE: return number of points inside..
+  return Py_BuildValue("i", count);
+}
+
+
+
+PyObject *gradient(PyObject *self, PyObject *args) {
+  //
+  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+  //
+
+  double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b;
+  PyObject *result;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
+                        &q0, &q1, &q2))
+    return NULL;
+
+
+  // Call underlying routine
+  _gradient(x0, y0, x1, y1, x2, y2,
+            q0, q1, q2, &a, &b);
+
+  // Return values a and b
+  result = Py_BuildValue("dd", a, b);
+  return result;
+}
+
+
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+
+  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {"gradient", gradient, METH_VARARGS, "Print out"},
+  {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
+  //{"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},
+  {"separate_points_by_polygon", separate_points_by_polygon, 
+                                 METH_VARARGS, "Print out"},
+  {NULL, NULL, 0, NULL}   /* sentinel */
+};
+
+
+
+// Module initialisation
+void initutil_ext(void){
+  Py_InitModule("util_ext", MethodTable);
+
+  import_array();     //Necessary for handling of NumPY structures
+}
+
+
+
+
+//OBSOLETE
 
 int _inside_polygon(int M,           // Number of points
@@ -144 +384 @@
 
 
-// Gateways to Python
-PyObject *point_on_line(PyObject *self, PyObject *args) {
-  //
-  // point_on_line(x, y, x0, y0, x1, y1)
-  //
-
-  double x, y, x0, y0, x1, y1;
-  int res;
-  PyObject *result;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1))
-    return NULL;
-
-
-  // Call underlying routine
-  res = _point_on_line(x, y, x0, y0, x1, y1);
-
-  // Return values a and b
-  result = Py_BuildValue("i", res);
-  return result;
-}
-
-
-
 PyObject *inside_polygon(PyObject *self, PyObject *args) {
   //def inside_polygon(point, polygon, closed, verbose, one_point):
@@ -242 +457 @@
   return Py_BuildValue("i", count);
 }
-
-
-
-PyObject *gradient(PyObject *self, PyObject *args) {
-  //
-  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-  //
-
-  double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b;
-  PyObject *result;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
-                        &q0, &q1, &q2))
-    return NULL;
-
-
-  // Call underlying routine
-  _gradient(x0, y0, x1, y1, x2, y2,
-            q0, q1, q2, &a, &b);
-
-  // Return values a and b
-  result = Py_BuildValue("dd", a, b);
-  return result;
-}
-
-
-// Method table for python module
-static struct PyMethodDef MethodTable[] = {
-  /* The cast of the function is necessary since PyCFunction values
-   * only take two PyObject* parameters, and rotate() takes
-   * three.
-   */
-
-  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
-  {"gradient", gradient, METH_VARARGS, "Print out"},
-  {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
-  {"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},
-  {NULL, NULL, 0, NULL}   /* sentinel */
-};
-
-
-
-// Module initialisation
-void initutil_ext(void){
-  Py_InitModule("util_ext", MethodTable);
-
-  import_array();     //Necessary for handling of NumPY structures
-}