source: inundation/utilities/polygon_ext.c @ 1910

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

Made utilities a Python package
Added numerical_tools (and test) containing ensure_numeric from the old util module.
Added test_polygon.py

File size: 10.7 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  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
233 
234  //Call underlying routine
235  count = _separate_points_by_polygon(M, N,
236                                      (double*) points -> data,
237                                      (double*) polygon -> data,
238                                      (long*) indices -> data,
239                                      closed, verbose);
240 
241  //NOTE: return number of points inside..
242  return Py_BuildValue("i", count);
243}
244
245
246
247// Method table for python module
248static struct PyMethodDef MethodTable[] = {
249  /* The cast of the function is necessary since PyCFunction values
250   * only take two PyObject* parameters, and rotate() takes
251   * three.
252   */
253
254  {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
255  {"separate_points_by_polygon", separate_points_by_polygon, 
256                                 METH_VARARGS, "Print out"},
257  {NULL, NULL, 0, NULL}   /* sentinel */
258};
259
260
261
262// Module initialisation
263void initpolygon_ext(void){
264  Py_InitModule("polygon_ext", MethodTable);
265
266  import_array();     //Necessary for handling of NumPY structures
267}
268
269
270
271
272//OBSOLETE
273
274int _inside_polygon(int M,           // Number of points
275                    int N,           // Number of polygon vertices
276                    double* points,
277                    double* polygon,
278                    long* indices,    // M-Array for storage indices
279                    int closed,
280                    int verbose) {
281
282  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
283  int i, j, k, count, inside;
284
285  //Find min and max of poly used for optimisation when points
286  //are far away from polygon
287
288  minpx = polygon[0]; maxpx = minpx;
289  minpy = polygon[1]; maxpy = minpy;
290
291  for (i=0; i<N; i++) {
292    px_i = polygon[2*i];
293    py_i = polygon[2*i + 1];
294
295    if (px_i < minpx) minpx = px_i;
296    if (px_i > maxpx) maxpx = px_i;
297    if (py_i < minpy) minpy = py_i;
298    if (py_i > maxpy) maxpy = py_i;
299  }
300
301  //Begin main loop (for each point)
302  count = 0;
303  for (k=0; k<M; k++) {
304    if (verbose){
305      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
306    }
307   
308    x = points[2*k];
309    y = points[2*k + 1];
310
311    inside = 0;
312
313    //Optimisation
314    if ((x > maxpx) || (x < minpx)) continue;
315    if ((y > maxpy) || (y < minpy)) continue;
316
317    //Check polygon
318    for (i=0; i<N; i++) {
319      //printf("k,i=%d,%d\n", k, i);
320      j = (i+1)%N;
321
322      px_i = polygon[2*i];
323      py_i = polygon[2*i+1];
324      px_j = polygon[2*j];
325      py_j = polygon[2*j+1];
326
327      //Check for case where point is contained in line segment
328      if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
329        if (closed == 1) {
330          inside = 1;
331        } else {
332          inside = 0;
333        }
334        break;
335      } else {
336        //Check if truly inside polygon
337        if ( ((py_i < y) && (py_j >= y)) ||
338             ((py_j < y) && (py_i >= y)) ) {
339          if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
340            inside = 1-inside;
341        }
342      }
343    }
344    if (inside == 1) {
345      indices[count] = k;
346      count++;
347    }
348  } // End k
349
350  return count;
351}
352
353
354
355PyObject *inside_polygon(PyObject *self, PyObject *args) {
356  //def inside_polygon(point, polygon, closed, verbose, one_point):
357  //  """Determine whether points are inside or outside a polygon
358  //
359  //  Input:
360  //     point - Tuple of (x, y) coordinates, or list of tuples
361  //     polygon - list of vertices of polygon
362  //     one_poin - If True Boolean value should be returned
363  //                If False, indices of points inside returned
364  //     closed - (optional) determine whether points on boundary should be
365  //     regarded as belonging to the polygon (closed = True)
366  //     or not (closed = False)
367
368  //
369  //  Output:
370  //     If one point is considered, True or False is returned.
371  //     If multiple points are passed in, the function returns indices
372  //     of those points that fall inside the polygon
373  //
374  //  Examples:
375  //     inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] )
376  //     will evaluate to True as the point 0.5, 0.5 is inside the unit square
377  //
378  //     inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
379  //     will return the indices [0, 2] as only the first and the last point
380  //     is inside the unit square
381  //
382  //  Remarks:
383  //     The vertices may be listed clockwise or counterclockwise and
384  //     the first point may optionally be repeated.
385  //     Polygons do not need to be convex.
386  //     Polygons can have holes in them and points inside a hole is
387  //     regarded as being outside the polygon.
388  //
389  //
390  //  Algorithm is based on work by Darel Finley,
391  //  http://www.alienryderflex.com/polygon/
392  //
393  //
394
395  PyArrayObject
396    *point,
397    *polygon,
398    *indices;
399
400  int closed, verbose; //Flags
401  int count, M, N;
402
403  // Convert Python arguments to C
404  if (!PyArg_ParseTuple(args, "OOOii",
405                        &point,
406                        &polygon,
407                        &indices,
408                        &closed,
409                        &verbose))
410    return NULL;
411
412  M = point -> dimensions[0];    //Number of points
413  N = polygon -> dimensions[0];  //Number of vertices in polygon
414
415  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
416 
417  //Call underlying routine
418  count = _inside_polygon(M, N,
419                          (double*) point -> data,
420                          (double*) polygon -> data,
421                          (long*) indices -> data,
422                          closed, verbose);
423
424  //NOTE: return number of points inside..
425  //printf("COunt=%d\n", count);
426  return Py_BuildValue("i", count);
427}
Note: See TracBrowser for help on using the repository browser.