source: anuga_core/source/anuga/utilities/polygon_ext.c @ 5403

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

Added tolerances for point_on_line and changed interpolate_polyline to use a general relative tolerance.

File size: 11.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: 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                    double rtol,
25                    double atol) {
26  /*Determine whether a point is on a line segment
27
28    Input: x, y, x0, x0, x1, y1: where
29        point is given by x, y
30        line is given by (x0, y0) and (x1, y1)
31
32  */
33
34  double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
35  double nominator, denominator;
36  int is_parallel;
37
38  a0 = x - x0;
39  a1 = y - y0;
40
41  a_normal0 = a1;
42  a_normal1 = -a0;
43
44  b0 = x1 - x0;
45  b1 = y1 - y0;
46
47  nominator = fabs(a_normal0*b0 + a_normal1*b1);
48  denominator = b0*b0 + b1*b1;
49 
50  // Determine if line is parallel to point vector up to a tolerance
51  is_parallel = 0;
52  if (denominator == 0.0) {
53    // Use absolute tolerance
54    if (nominator <= atol) {
55      is_parallel = 1;
56    }
57  } else {
58    // Denominator is positive - use relative tolerance
59    if (nominator/denominator <= rtol) {
60      is_parallel = 1;
61    }   
62  }
63   
64  if (is_parallel) {
65    // Point is somewhere on the infinite extension of the line
66    // subject to specified absolute tolerance
67
68    len_a = sqrt(a0*a0 + a1*a1);
69    len_b = sqrt(b0*b0 + b1*b1);
70
71    if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
72      return 1;
73    } else {
74      return 0;
75    }
76  } else {
77    return 0;
78  }
79}
80
81
82
83/*
84WORK IN PROGRESS TO OPTIMISE INTERSECTION
85int __intersection(double x0, double y0,
86                   double x1, double y1) {
87
88
89    x0 = line0[0,0]; y0 = line0[0,1]
90    x1 = line0[1,0]; y1 = line0[1,1]
91
92    x2 = line1[0,0]; y2 = line1[0,1]
93    x3 = line1[1,0]; y3 = line1[1,1]
94
95    denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
96    u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
97    u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
98       
99    if allclose(denom, 0.0):
100        # Lines are parallel - check if they coincide on a shared a segment
101
102        if allclose( [u0, u1], 0.0 ):
103            # We now know that the lines if continued coincide
104            # The remaining check will establish if the finite lines share a segment
105
106            line0_starts_on_line1 = line0_ends_on_line1 =\
107            line1_starts_on_line0 = line1_ends_on_line0 = False
108               
109            if point_on_line([x0, y0], line1):
110                line0_starts_on_line1 = True
111
112            if point_on_line([x1, y1], line1):
113                line0_ends_on_line1 = True
114 
115            if point_on_line([x2, y2], line0):
116                line1_starts_on_line0 = True
117
118            if point_on_line([x3, y3], line0):
119                line1_ends_on_line0 = True                               
120
121            if not(line0_starts_on_line1 or line0_ends_on_line1\
122               or line1_starts_on_line0 or line1_ends_on_line0):
123                # Lines are parallel and would coincide if extended, but not as they are.
124                return 3, None
125
126
127            # One line fully included in the other. Use direction of included line
128            if line0_starts_on_line1 and line0_ends_on_line1:
129                # Shared segment is line0 fully included in line1
130                segment = array([[x0, y0], [x1, y1]])               
131
132            if line1_starts_on_line0 and line1_ends_on_line0:
133                # Shared segment is line1 fully included in line0
134                segment = array([[x2, y2], [x3, y3]])
135           
136
137            # Overlap with lines are oriented the same way
138            if line0_starts_on_line1 and line1_ends_on_line0:
139                # Shared segment from line0 start to line 1 end
140                segment = array([[x0, y0], [x3, y3]])
141
142            if line1_starts_on_line0 and line0_ends_on_line1:
143                # Shared segment from line1 start to line 0 end
144                segment = array([[x2, y2], [x1, y1]])                               
145
146
147            # Overlap in opposite directions - use direction of line0
148            if line0_starts_on_line1 and line1_starts_on_line0:
149                # Shared segment from line0 start to line 1 end
150                segment = array([[x0, y0], [x2, y2]])
151
152            if line0_ends_on_line1 and line1_ends_on_line0:
153                # Shared segment from line0 start to line 1 end
154                segment = array([[x3, y3], [x1, y1]])               
155
156               
157            return 2, segment
158        else:
159            # Lines are parallel but they do not coincide
160            return 4, None #FIXME (Ole): Add distance here instead of None
161           
162    else:
163        # Lines are not parallel or coinciding
164        u0 = u0/denom
165        u1 = u1/denom       
166
167        x = x0 + u0*(x1-x0)
168        y = y0 + u0*(y1-y0)
169
170        # Sanity check - can be removed to speed up if needed
171        assert allclose(x, x2 + u1*(x3-x2))
172        assert allclose(y, y2 + u1*(y3-y2))       
173
174        # Check if point found lies within given line segments
175        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0:
176            # We have intersection
177
178            return 1, array([x, y])
179        else:
180            # No intersection
181            return 0, None
182
183
184}
185*/
186
187
188int __separate_points_by_polygon(int M,     // Number of points
189                                int N,     // Number of polygon vertices
190                                double* points,
191                                double* polygon,
192                                long* indices,  // M-Array for storage indices
193                                int closed,
194                                int verbose) {
195
196  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j, rtol=0.0, atol=0.0;
197  int i, j, k, outside_index, inside_index, inside;
198
199  //Find min and max of poly used for optimisation when points
200  //are far away from polygon
201
202  minpx = polygon[0]; maxpx = minpx;
203  minpy = polygon[1]; maxpy = minpy;
204
205  for (i=0; i<N; i++) {
206    px_i = polygon[2*i];
207    py_i = polygon[2*i + 1];
208
209    if (px_i < minpx) minpx = px_i;
210    if (px_i > maxpx) maxpx = px_i;
211    if (py_i < minpy) minpy = py_i;
212    if (py_i > maxpy) maxpy = py_i;
213  }
214
215  //Begin main loop (for each point)
216  inside_index = 0;    //Keep track of points inside
217  outside_index = M-1; //Keep track of points outside (starting from end)   
218  if (verbose){
219     printf("Separating %d points\n", M);
220  } 
221  for (k=0; k<M; k++) {
222    if (verbose){
223      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
224    }
225   
226    x = points[2*k];
227    y = points[2*k + 1];
228
229    inside = 0;
230
231    //Optimisation
232    if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) {
233      //Nothing
234    } else {   
235      //Check polygon
236      for (i=0; i<N; i++) {
237        //printf("k,i=%d,%d\n", k, i);
238        j = (i+1)%N;
239
240        px_i = polygon[2*i];
241        py_i = polygon[2*i+1];
242        px_j = polygon[2*j];
243        py_j = polygon[2*j+1];
244
245        //Check for case where point is contained in line segment
246        if (__point_on_line(x, y, px_i, py_i, px_j, py_j, rtol, atol)) {
247          if (closed == 1) {
248            inside = 1;
249          } else {
250            inside = 0;
251          }
252          break;
253        } else {
254          //Check if truly inside polygon
255          if ( ((py_i < y) && (py_j >= y)) ||
256               ((py_j < y) && (py_i >= y)) ) {
257            if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
258              inside = 1-inside;
259          }
260        }
261      }
262    } 
263    if (inside == 1) {
264      indices[inside_index] = k;
265      inside_index += 1;
266    } else {
267      indices[outside_index] = k;
268      outside_index -= 1;   
269    }
270  } // End k
271
272  return inside_index;
273}
274
275
276
277// Gateways to Python
278PyObject *_point_on_line(PyObject *self, PyObject *args) {
279  //
280  // point_on_line(x, y, x0, y0, x1, y1)
281  //
282
283  double x, y, x0, y0, x1, y1, rtol, atol;
284  int res;
285  PyObject *result;
286
287  // Convert Python arguments to C
288  if (!PyArg_ParseTuple(args, "dddddddd", &x, &y, &x0, &y0, &x1, &y1, &rtol, &atol)) {
289    PyErr_SetString(PyExc_RuntimeError, 
290                    "point_on_line could not parse input");   
291    return NULL;
292  }
293
294
295  // Call underlying routine
296  res = __point_on_line(x, y, x0, y0, x1, y1, rtol, atol);
297
298  // Return values a and b
299  result = Py_BuildValue("i", res);
300  return result;
301}
302
303
304/*
305PyObject *_intersection(PyObject *self, PyObject *args) {
306  //
307  // intersection(x0, y0, x1, y1)
308  //
309
310  double x, y, x0, y0, x1, y1;
311  int res;
312  PyObject *result;
313
314  // Convert Python arguments to C
315  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) {
316    PyErr_SetString(PyExc_RuntimeError,
317                    "point_on_line could not parse input");   
318    return NULL;
319  }
320
321
322  // Call underlying routine
323  res = __intersection(x, y, x0, y0, x1, y1);
324
325  // Return values a and b
326  result = Py_BuildValue("i", res);
327  return result;
328}
329*/
330
331
332PyObject *_separate_points_by_polygon(PyObject *self, PyObject *args) {
333  //def separate_points_by_polygon(points, polygon, closed, verbose, one_point):
334  //  """Determine whether points are inside or outside a polygon
335  //
336  //  Input:
337  //     point - Tuple of (x, y) coordinates, or list of tuples
338  //     polygon - list of vertices of polygon
339  //     closed - (optional) determine whether points on boundary should be
340  //     regarded as belonging to the polygon (closed = True)
341  //     or not (closed = False)
342
343  //
344  //  Output:
345  //     indices: array of same length as points with indices of points falling
346  //     inside the polygon listed from the beginning and indices of points
347  //     falling outside listed from the end.
348  //     
349  //     count: count of points falling inside the polygon
350  //     
351  //     The indices of points inside are obtained as indices[:count]
352  //     The indices of points outside are obtained as indices[count:]       
353  //
354  //  Examples:
355  //     separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
356  //     will return the indices [0, 2, 1] as only the first and the last point
357  //     is inside the unit square
358  //
359  //  Remarks:
360  //     The vertices may be listed clockwise or counterclockwise and
361  //     the first point may optionally be repeated.
362  //     Polygons do not need to be convex.
363  //     Polygons can have holes in them and points inside a hole is
364  //     regarded as being outside the polygon.
365  //
366  //
367  //  Algorithm is based on work by Darel Finley,
368  //  http://www.alienryderflex.com/polygon/
369  //
370  //
371
372  PyArrayObject
373    *points,
374    *polygon,
375    *indices;
376
377  int closed, verbose; //Flags
378  int count, M, N;
379
380  // Convert Python arguments to C
381  if (!PyArg_ParseTuple(args, "OOOii",
382                        &points,
383                        &polygon,
384                        &indices,
385                        &closed,
386                        &verbose)) {
387   
388
389    PyErr_SetString(PyExc_RuntimeError, 
390                    "separate_points_by_polygon could not parse input");
391    return NULL;
392  }
393
394  M = points -> dimensions[0];   //Number of points
395  N = polygon -> dimensions[0];  //Number of vertices in polygon
396
397  //FIXME (Ole): This isn't send to Python's sys.stdout
398  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
399 
400  //Call underlying routine
401  count = __separate_points_by_polygon(M, N,
402                                       (double*) points -> data,
403                                       (double*) polygon -> data,
404                                       (long*) indices -> data,
405                                       closed, verbose);
406 
407  //NOTE: return number of points inside..
408  return Py_BuildValue("i", count);
409}
410
411
412
413// Method table for python module
414static struct PyMethodDef MethodTable[] = {
415  /* The cast of the function is necessary since PyCFunction values
416   * only take two PyObject* parameters, and rotate() takes
417   * three.
418   */
419
420  {"_point_on_line", _point_on_line, METH_VARARGS, "Print out"},
421  //{"_intersection", _intersection, METH_VARARGS, "Print out"}, 
422  {"_separate_points_by_polygon", _separate_points_by_polygon, 
423                                 METH_VARARGS, "Print out"},
424  {NULL, NULL, 0, NULL}   /* sentinel */
425};
426
427
428
429// Module initialisation
430void initpolygon_ext(void){
431  Py_InitModule("polygon_ext", MethodTable);
432
433  import_array();     //Necessary for handling of NumPY structures
434}
435
436
437
Note: See TracBrowser for help on using the repository browser.