source:inundation/utilities/polygon_ext.c@1906

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

Began moving polygon functionality to separate module and created a new package utilities for general tools used by AnuGa?

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