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

Last change on this file since 6534 was 6534, checked in by ole, 15 years ago

Added new function is_inside_triangle. This helps the fitting functions.

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