Changeset 999 for inundation/ga/storm_surge/pyvolution

Timestamp:

Mar 4, 2005, 12:20:11 PM (20 years ago)

Author:

ole

Message:

Optimised inside_polygon in C and tested

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

@@ -826 +826 @@
         polygon = [[0,0], [1,0], [1,1], [0,1]]
 
-
-        R = inside_polygon( (0.5, 0.5), polygon)
-        print 'Inside = ', R
-        print
-        
-        R = inside_polygon( (0.5, 0.5), polygon)
-        print 'Inside = ', R
-        print
-
-        import sys; sys.exit() 
         assert inside_polygon( (0.5, 0.5), polygon )
         assert not inside_polygon( (0.5, 1.5), polygon )   
@@ -953 +943 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test_inside_polygon_main')
+    #suite = unittest.makeSuite(TestCase,'test_inside_polygon_main')
+    suite = unittest.makeSuite(TestCase,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)

inundation/ga/storm_surge/pyvolution/util.py

@@ -87 +87 @@
 
     #FIXME (OLE): Should check origin of domain against that of file
-
+    #In fact, this is where origin should be converted to that of domain
+    #Also, check that file covers domain fully.
     
     assert type(filename) == type(''),\
@@ -974 +975 @@
     py = polygon[:,1]   
 
-    #Used for an optimisation when points are far away from polgon
+    #Used for an optimisation when points are far away from polygon
     minpx = min(px); maxpx = max(px)
     minpy = min(py); maxpy = max(py)    
@@ -1022 +1023 @@
              
 
+
+def inside_polygon_c(point, polygon, closed = True, verbose = False):
+    """Determine whether points are inside or outside a polygon
+
+    C-wrapper
+    """    
+
+    from Numeric import array, Float, reshape, zeros, Int
+    
+    
+    #Input checks
+    try:
+        point = array(point).astype(Float)         
+    except:
+        msg = 'Point %s could not be converted to Numeric array' %point
+        raise msg        
+        
+    try:
+        polygon = array(polygon).astype(Float)             
+    except:
+        msg = 'Polygon %s could not be converted to Numeric array' %polygon
+        raise msg               
+    
+    assert len(polygon.shape) == 2,\
+       'Polygon array must be a 2d array of vertices: %s' %polygon
+
+        
+    assert polygon.shape[1] == 2,\
+       'Polygon array must have two columns: %s' %polygon    
+
+
+    if len(point.shape) == 1:
+        one_point = True
+        points = reshape(point, (1,2))
+    else:       
+        one_point = False
+        points = point
+
+    from util_ext import inside_polygon
+
+    indices = zeros( points.shape[0], Int )
+
+    count = inside_polygon(points, polygon, indices,
+                           int(closed), int(verbose))
+
+    #print 'O', point, count
+
+    if one_point:
+        return count == 1
+    else:
+        return indices[:count]
+
 class Polygon_function:
     """Create callable object f: x,y -> z, where a,y,z are vectors and 
@@ -1134 +1187 @@
     Examples:
        populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
-       will return five randomly sleected points inside the unit square 
+       will return five randomly selected points inside the unit square 
     """
 
@@ -1183 +1236 @@
 
 
-
-
 ##############################################
 #Initialise module
@@ -1190 +1241 @@
 import compile
 if compile.can_use_C_extension('util_ext.c'):
-    from util_ext import gradient, point_on_line 
+    from util_ext import gradient, point_on_line
+    inside_polygon = inside_polygon_c
 else:
     gradient = gradient_python

inundation/ga/storm_surge/pyvolution/util_ext.c

@@ -59 +59 @@
 
 
+int _inside_polygon(int M,           // Number of points
+                    int N,           // Number of polygon vertices
+                    double* points,
+                    double* polygon,
+                    int* 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, count, 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)
+  count = 0;
+  for (k=0; k<M; k++) {
+    //FIXME: Do this later
+    //if (verbose){
+    //  if (k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
+        
+    x = points[2*k];
+    y = points[2*k + 1];        
+
+    inside = 0;
+
+    //Optimisation
+    if ((x > maxpx) || (x < minpx)) continue;
+    if ((y > maxpy) || (y < minpy)) continue;        
+
+    //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[count] = k;
+      count++;
+    }       
+  } // End k
+  
+  return count;
+}
+
+
+
 // Gateways to Python
 PyObject *point_on_line(PyObject *self, PyObject *args) {
@@ -81 +163 @@
   return result;
 }
+
+
+
+PyObject *inside_polygon(PyObject *self, PyObject *args) {
+  //def inside_polygon(point, 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
+  //     one_poin - If True Boolean value should be returned 
+  //                If False, indices of points inside returned
+  //     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 inside the polygon     
+  //     
+  //  Examples:
+  //     inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] )
+  //     will evaluate to True as the point 0.5, 0.5 is inside the unit square 
+  //     
+  //     inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
+  //     will return the indices [0, 2] 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 
+    *point, 
+    *polygon, 
+    *indices;
+    
+  int closed, verbose; //Flags     
+  int count, M, N;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "OOOii", 
+                        &point, 
+                        &polygon, 
+                        &indices, 
+                        &closed, 
+                        &verbose))  
+    return NULL;
+
+  M = point -> dimensions[0];    //Number of points          
+  N = polygon -> dimensions[0];  //Number of vertices in polygon        
+
+  //printf("M=%d, N=%d\n", M, N);
+  // Call underlying routine
+  count = _inside_polygon(M, N, 
+                          (double*) point -> data,
+                          (double*) polygon -> data,
+                          (int*) indices -> data,
+                          closed, verbose); 
+                           
+  //NOTE: return number of points inside..
+  //printf("COunt=%d\n", count);
+  return Py_BuildValue("i", count);
+}
+
 
 
@@ -117 +275 @@
   {"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 */
 };