source: inundation/utilities/polygon_ext.c @ 2867

Last change on this file since 2867 was 2867, checked in by ole, 18 years ago

FIXME

File size: 6.6 KB
Line 
1// Python - C extension for polygon module.
2//
3// To compile (Python2.3):
4//  gcc -c polygon_ext.c -I/usr/include/python2.3 -o polygon_ext.o -Wall -O
5//  gcc -shared polygon_ext.o  -o polygon_ext.so
6//
7// See the module polygon.py
8//
9//
10// Ole Nielsen, GA 2004
11//
12//NOTE: On 64 bit systems use long* instead of int* for Numeric arrays
13//this will also work on 32 bit systems
14
15
16#include "Python.h"
17#include "Numeric/arrayobject.h"
18#include "math.h"
19
20
21int _point_on_line(double x, double y,
22                   double x0, double y0,
23                   double x1, double y1) {
24  /*Determine whether a point is on a line segment
25
26    Input: x, y, x0, x0, x1, y1: where
27        point is given by x, y
28        line is given by (x0, y0) and (x1, y1)
29
30  */
31
32  double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
33
34  a0 = x - x0;
35  a1 = y - y0;
36
37  a_normal0 = a1;
38  a_normal1 = -a0;
39
40  b0 = x1 - x0;
41  b1 = y1 - y0;
42
43  if ( a_normal0*b0 + a_normal1*b1 == 0 ) {
44    //Point is somewhere on the infinite extension of the line
45
46    len_a = sqrt(a0*a0 + a1*a1);
47    len_b = sqrt(b0*b0 + b1*b1);
48
49    if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
50      return 1;
51    } else {
52      return 0;
53    }
54  } else {
55    return 0;
56  }
57}
58
59
60
61int _separate_points_by_polygon(int M,     // Number of points
62                                int N,     // Number of polygon vertices
63                                double* points,
64                                double* polygon,
65                                long* indices,  // M-Array for storage indices
66                                int closed,
67                                int verbose) {
68
69  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
70  int i, j, k, outside_index, inside_index, inside;
71
72  //Find min and max of poly used for optimisation when points
73  //are far away from polygon
74
75  minpx = polygon[0]; maxpx = minpx;
76  minpy = polygon[1]; maxpy = minpy;
77
78  for (i=0; i<N; i++) {
79    px_i = polygon[2*i];
80    py_i = polygon[2*i + 1];
81
82    if (px_i < minpx) minpx = px_i;
83    if (px_i > maxpx) maxpx = px_i;
84    if (py_i < minpy) minpy = py_i;
85    if (py_i > maxpy) maxpy = py_i;
86  }
87
88  //Begin main loop (for each point)
89  inside_index = 0;    //Keep track of points inside
90  outside_index = M-1; //Keep track of points outside (starting from end)   
91  for (k=0; k<M; k++) {
92    if (verbose){
93      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
94    }
95   
96    x = points[2*k];
97    y = points[2*k + 1];
98
99    inside = 0;
100
101    //Optimisation
102    if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) {
103      //Nothing
104    } else {   
105      //Check polygon
106      for (i=0; i<N; i++) {
107        //printf("k,i=%d,%d\n", k, i);
108        j = (i+1)%N;
109
110        px_i = polygon[2*i];
111        py_i = polygon[2*i+1];
112        px_j = polygon[2*j];
113        py_j = polygon[2*j+1];
114
115        //Check for case where point is contained in line segment
116        if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
117          if (closed == 1) {
118            inside = 1;
119          } else {
120            inside = 0;
121          }
122          break;
123        } else {
124          //Check if truly inside polygon
125          if ( ((py_i < y) && (py_j >= y)) ||
126               ((py_j < y) && (py_i >= y)) ) {
127            if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
128              inside = 1-inside;
129          }
130        }
131      }
132    } 
133    if (inside == 1) {
134      indices[inside_index] = k;
135      inside_index += 1;
136    } else {
137      indices[outside_index] = k;
138      outside_index -= 1;   
139    }
140  } // End k
141
142  return inside_index;
143}
144
145
146
147// Gateways to Python
148PyObject *point_on_line(PyObject *self, PyObject *args) {
149  //
150  // point_on_line(x, y, x0, y0, x1, y1)
151  //
152
153  double x, y, x0, y0, x1, y1;
154  int res;
155  PyObject *result;
156
157  // Convert Python arguments to C
158  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1))
159    return NULL;
160
161
162  // Call underlying routine
163  res = _point_on_line(x, y, x0, y0, x1, y1);
164
165  // Return values a and b
166  result = Py_BuildValue("i", res);
167  return result;
168}
169
170
171
172PyObject *separate_points_by_polygon(PyObject *self, PyObject *args) {
173  //def separate_points_by_polygon(points, polygon, closed, verbose, one_point):
174  //  """Determine whether points are inside or outside a polygon
175  //
176  //  Input:
177  //     point - Tuple of (x, y) coordinates, or list of tuples
178  //     polygon - list of vertices of polygon
179  //     closed - (optional) determine whether points on boundary should be
180  //     regarded as belonging to the polygon (closed = True)
181  //     or not (closed = False)
182
183  //
184  //  Output:
185  //     indices: array of same length as points with indices of points falling
186  //     inside the polygon listed from the beginning and indices of points
187  //     falling outside listed from the end.
188  //     
189  //     count: count of points falling inside the polygon
190  //     
191  //     The indices of points inside are obtained as indices[:count]
192  //     The indices of points outside are obtained as indices[count:]       
193  //
194  //  Examples:
195  //     separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
196  //     will return the indices [0, 2, 1] as only the first and the last point
197  //     is inside the unit square
198  //
199  //  Remarks:
200  //     The vertices may be listed clockwise or counterclockwise and
201  //     the first point may optionally be repeated.
202  //     Polygons do not need to be convex.
203  //     Polygons can have holes in them and points inside a hole is
204  //     regarded as being outside the polygon.
205  //
206  //
207  //  Algorithm is based on work by Darel Finley,
208  //  http://www.alienryderflex.com/polygon/
209  //
210  //
211
212  PyArrayObject
213    *points,
214    *polygon,
215    *indices;
216
217  int closed, verbose; //Flags
218  int count, M, N;
219
220  // Convert Python arguments to C
221  if (!PyArg_ParseTuple(args, "OOOii",
222                        &points,
223                        &polygon,
224                        &indices,
225                        &closed,
226                        &verbose))
227    return NULL;
228
229  M = points -> dimensions[0];   //Number of points
230  N = polygon -> dimensions[0];  //Number of vertices in polygon
231
232  //FIXME (Ole): This isn't send to Python's sys.stdout
233  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
234 
235  //Call underlying routine
236  count = _separate_points_by_polygon(M, N,
237                                      (double*) points -> data,
238                                      (double*) polygon -> data,
239                                      (long*) indices -> data,
240                                      closed, verbose);
241 
242  //NOTE: return number of points inside..
243  return Py_BuildValue("i", count);
244}
245
246
247
248// Method table for python module
249static struct PyMethodDef MethodTable[] = {
250  /* The cast of the function is necessary since PyCFunction values
251   * only take two PyObject* parameters, and rotate() takes
252   * three.
253   */
254
255  {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
256  {"separate_points_by_polygon", separate_points_by_polygon, 
257                                 METH_VARARGS, "Print out"},
258  {NULL, NULL, 0, NULL}   /* sentinel */
259};
260
261
262
263// Module initialisation
264void initpolygon_ext(void){
265  Py_InitModule("polygon_ext", MethodTable);
266
267  import_array();     //Necessary for handling of NumPY structures
268}
269
270
271
Note: See TracBrowser for help on using the repository browser.