# source:anuga_core/source/anuga/utilities/polygon_ext.c@4804

Last change on this file since 4804 was 4804, checked in by ole, 16 years ago

Made order of tests predictable across platforms.
Bug fix in runcairns

File size: 6.9 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: We use long* instead of int* for Numeric arrays as this will work both
13//       for 64 as well as 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  if (verbose){
92     printf("Separating %d points\n", M);
93  }
94  for (k=0; k<M; k++) {
95    if (verbose){
96      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
97    }
98
99    x = points[2*k];
100    y = points[2*k + 1];
101
102    inside = 0;
103
104    //Optimisation
105    if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) {
106      //Nothing
107    } else {
108      //Check polygon
109      for (i=0; i<N; i++) {
110        //printf("k,i=%d,%d\n", k, i);
111        j = (i+1)%N;
112
113        px_i = polygon[2*i];
114        py_i = polygon[2*i+1];
115        px_j = polygon[2*j];
116        py_j = polygon[2*j+1];
117
118        //Check for case where point is contained in line segment
119        if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
120          if (closed == 1) {
121            inside = 1;
122          } else {
123            inside = 0;
124          }
125          break;
126        } else {
127          //Check if truly inside polygon
128          if ( ((py_i < y) && (py_j >= y)) ||
129               ((py_j < y) && (py_i >= y)) ) {
130            if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
131              inside = 1-inside;
132          }
133        }
134      }
135    }
136    if (inside == 1) {
137      indices[inside_index] = k;
138      inside_index += 1;
139    } else {
140      indices[outside_index] = k;
141      outside_index -= 1;
142    }
143  } // End k
144
145  return inside_index;
146}
147
148
149
150// Gateways to Python
151PyObject *point_on_line(PyObject *self, PyObject *args) {
152  //
153  // point_on_line(x, y, x0, y0, x1, y1)
154  //
155
156  double x, y, x0, y0, x1, y1;
157  int res;
158  PyObject *result;
159
160  // Convert Python arguments to C
161  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) {
162    PyErr_SetString(PyExc_RuntimeError,
163                    "point_on_line could not parse input");
164    return NULL;
165  }
166
167
168  // Call underlying routine
169  res = _point_on_line(x, y, x0, y0, x1, y1);
170
171  // Return values a and b
172  result = Py_BuildValue("i", res);
173  return result;
174}
175
176
177
178PyObject *separate_points_by_polygon(PyObject *self, PyObject *args) {
179  //def separate_points_by_polygon(points, polygon, closed, verbose, one_point):
180  //  """Determine whether points are inside or outside a polygon
181  //
182  //  Input:
183  //     point - Tuple of (x, y) coordinates, or list of tuples
184  //     polygon - list of vertices of polygon
185  //     closed - (optional) determine whether points on boundary should be
186  //     regarded as belonging to the polygon (closed = True)
187  //     or not (closed = False)
188
189  //
190  //  Output:
191  //     indices: array of same length as points with indices of points falling
192  //     inside the polygon listed from the beginning and indices of points
193  //     falling outside listed from the end.
194  //
195  //     count: count of points falling inside the polygon
196  //
197  //     The indices of points inside are obtained as indices[:count]
198  //     The indices of points outside are obtained as indices[count:]
199  //
200  //  Examples:
201  //     separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
202  //     will return the indices [0, 2, 1] as only the first and the last point
203  //     is inside the unit square
204  //
205  //  Remarks:
206  //     The vertices may be listed clockwise or counterclockwise and
207  //     the first point may optionally be repeated.
208  //     Polygons do not need to be convex.
209  //     Polygons can have holes in them and points inside a hole is
210  //     regarded as being outside the polygon.
211  //
212  //
213  //  Algorithm is based on work by Darel Finley,
214  //  http://www.alienryderflex.com/polygon/
215  //
216  //
217
218  PyArrayObject
219    *points,
220    *polygon,
221    *indices;
222
223  int closed, verbose; //Flags
224  int count, M, N;
225
226  // Convert Python arguments to C
227  if (!PyArg_ParseTuple(args, "OOOii",
228                        &points,
229                        &polygon,
230                        &indices,
231                        &closed,
232                        &verbose)) {
233
234
235    PyErr_SetString(PyExc_RuntimeError,
236                    "separate_points_by_polygon could not parse input");
237    return NULL;
238  }
239
240  M = points -> dimensions[0];   //Number of points
241  N = polygon -> dimensions[0];  //Number of vertices in polygon
242
243  //FIXME (Ole): This isn't send to Python's sys.stdout
244  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
245
246  //Call underlying routine
247  count = _separate_points_by_polygon(M, N,
248                                      (double*) points -> data,
249                                      (double*) polygon -> data,
250                                      (long*) indices -> data,
251                                      closed, verbose);
252
253  //NOTE: return number of points inside..
254  return Py_BuildValue("i", count);
255}
256
257
258
259// Method table for python module
260static struct PyMethodDef MethodTable[] = {
261  /* The cast of the function is necessary since PyCFunction values
262   * only take two PyObject* parameters, and rotate() takes
263   * three.
264   */
265
266  {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
267  {"separate_points_by_polygon", separate_points_by_polygon,
268                                 METH_VARARGS, "Print out"},
269  {NULL, NULL, 0, NULL}   /* sentinel */
270};
271
272
273
274// Module initialisation
275void initpolygon_ext(void){
276  Py_InitModule("polygon_ext", MethodTable);
277
278  import_array();     //Necessary for handling of NumPY structures
279}
280
281
282
Note: See TracBrowser for help on using the repository browser.