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

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

Wrote is_inside_triangle in C and tested it some more.

File size: 17.3 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
246
247int __is_inside_triangle(double* point,
248                         double* triangle,
249                         int closed,
250                         double rtol,
251                         double atol) {
252                         
253  double vx, vy, v0x, v0y, v1x, v1y;
254  double a00, a10, a01, a11, b0, b1;
255  double denom, alpha, beta;
256  int i, j, res;
257
258  // v0 = C-A
259  v0x = triangle[4]-triangle[0]; 
260  v0y = triangle[5]-triangle[1];
261 
262  // v1 = B-A   
263  v1x = triangle[2]-triangle[0]; 
264  v1y = triangle[3]-triangle[1];
265
266  // First check if point lies wholly inside triangle
267  a00 = v0x*v0x + v0y*v0y; // innerproduct(v0, v0)
268  a01 = v0x*v1x + v0y*v1y; // innerproduct(v0, v1)
269  a10 = a01;               // innerproduct(v1, v0)
270  a11 = v1x*v1x + v1y*v1y; // innerproduct(v1, v1)
271   
272  denom = a11*a00 - a01*a10;
273
274  if (fabs(denom) > 0.0) {
275    // v = point-A 
276    vx = point[0] - triangle[0]; 
277    vy = point[1] - triangle[1];     
278   
279    b0 = v0x*vx + v0y*vy; // innerproduct(v0, v)       
280    b1 = v1x*vx + v1y*vy; // innerproduct(v1, v)           
281   
282    alpha = (b0*a11 - b1*a01)/denom;
283    beta = (b1*a00 - b0*a10)/denom;       
284   
285    if ((alpha > 0.0) && (beta > 0.0) && (alpha+beta < 1.0)) return 1;
286  }
287
288  if (closed) {
289    // Check if point lies on one of the edges
290       
291    for (i=0; i<3; i++) {
292      j = (i+1) % 3; // Circular index into triangle array
293      res = __point_on_line(point[0], point[1],
294                            triangle[2*i], triangle[2*i+1], 
295                            triangle[2*j], triangle[2*j+1],                         
296                            rtol, atol);
297      if (res) return 1;
298    }
299  }
300               
301  // Default return if point is outside triangle                         
302  return 0;                                             
303}                                                                             
304
305
306int __separate_points_by_polygon(int M,     // Number of points
307                                 int N,     // Number of polygon vertices
308                                 double* points,
309                                 double* polygon,
310                                 long* indices,  // M-Array for storage indices
311                                 int closed,
312                                 int verbose) {
313
314  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j, rtol=0.0, atol=0.0;
315  int i, j, k, outside_index, inside_index, inside;
316
317  // Find min and max of poly used for optimisation when points
318  // are far away from polygon
319 
320  // FIXME(Ole): Pass in rtol and atol from Python
321
322  minpx = polygon[0]; maxpx = minpx;
323  minpy = polygon[1]; maxpy = minpy;
324
325  for (i=0; i<N; i++) {
326    px_i = polygon[2*i];
327    py_i = polygon[2*i + 1];
328
329    if (px_i < minpx) minpx = px_i;
330    if (px_i > maxpx) maxpx = px_i;
331    if (py_i < minpy) minpy = py_i;
332    if (py_i > maxpy) maxpy = py_i;
333  }
334
335  // Begin main loop (for each point)
336  inside_index = 0;    // Keep track of points inside
337  outside_index = M-1; // Keep track of points outside (starting from end)   
338  if (verbose){
339     printf("Separating %d points\n", M);
340  } 
341  for (k=0; k<M; k++) {
342    if (verbose){
343      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
344    }
345   
346    x = points[2*k];
347    y = points[2*k + 1];
348
349    inside = 0;
350
351    // Optimisation
352    if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) {
353      // Nothing
354    } else {   
355      // Check polygon
356      for (i=0; i<N; i++) {
357        j = (i+1)%N;
358
359        px_i = polygon[2*i];
360        py_i = polygon[2*i+1];
361        px_j = polygon[2*j];
362        py_j = polygon[2*j+1];
363
364        // Check for case where point is contained in line segment
365        if (__point_on_line(x, y, px_i, py_i, px_j, py_j, rtol, atol)) {
366          if (closed == 1) {
367            inside = 1;
368          } else {
369            inside = 0;
370          }
371          break;
372        } else {
373          // Check if truly inside polygon
374          if ( ((py_i < y) && (py_j >= y)) ||
375               ((py_j < y) && (py_i >= y)) ) {
376            if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
377              inside = 1-inside;
378          }
379        }
380      }
381    } 
382    if (inside == 1) {
383      indices[inside_index] = k;
384      inside_index += 1;
385    } else {
386      indices[outside_index] = k;
387      outside_index -= 1;   
388    }
389  } // End k
390
391  return inside_index;
392}
393
394
395
396// Gateways to Python
397PyObject *_point_on_line(PyObject *self, PyObject *args) {
398  //
399  // point_on_line(x, y, x0, y0, x1, y1)
400  //
401
402  double x, y, x0, y0, x1, y1, rtol, atol;
403  int res;
404  PyObject *result;
405
406  // Convert Python arguments to C
407  if (!PyArg_ParseTuple(args, "dddddddd", &x, &y, &x0, &y0, &x1, &y1, &rtol, &atol)) {
408    PyErr_SetString(PyExc_RuntimeError, 
409                    "point_on_line could not parse input");   
410    return NULL;
411  }
412
413
414  // Call underlying routine
415  res = __point_on_line(x, y, x0, y0, x1, y1, rtol, atol);
416
417  // Return result
418  result = Py_BuildValue("i", res);
419  return result;
420}
421
422
423
424// Gateways to Python
425PyObject *_interpolate_polyline(PyObject *self, PyObject *args) {
426  //
427  // _interpolate_polyline(data, polyline_nodes, gauge_neighbour_id, interpolation_points
428  //                       interpolated_values):
429  //
430
431 
432  PyArrayObject
433    *data,
434    *polyline_nodes,
435    *gauge_neighbour_id,
436    *interpolation_points,
437    *interpolated_values;
438
439  double rtol, atol; 
440  int number_of_nodes, number_of_points, res;
441 
442  // Convert Python arguments to C
443  if (!PyArg_ParseTuple(args, "OOOOOdd",
444                        &data,
445                        &polyline_nodes,
446                        &gauge_neighbour_id,
447                        &interpolation_points,
448                        &interpolated_values,
449                        &rtol,
450                        &atol)) {
451   
452    PyErr_SetString(PyExc_RuntimeError, 
453                    "_interpolate_polyline could not parse input");
454    return NULL;
455  }
456
457  number_of_nodes = polyline_nodes -> dimensions[0];  // Number of nodes in polyline
458  number_of_points = interpolation_points -> dimensions[0];   //Number of points
459 
460
461  // Call underlying routine
462  res = __interpolate_polyline(number_of_nodes,
463                               number_of_points,
464                               (double*) data -> data,
465                               (double*) polyline_nodes -> data,
466                               (long*) gauge_neighbour_id -> data,
467                               (double*) interpolation_points -> data,                         
468                               (double*) interpolated_values -> data,
469                               rtol,
470                               atol);                                                 
471
472  // Return
473  return Py_BuildValue(""); 
474}
475
476
477
478
479     
480PyObject *_is_inside_triangle(PyObject *self, PyObject *args) {
481  //
482  // _is_inside_triangle(point, triangle, int(closed), rtol, atol)
483  //
484
485 
486  PyArrayObject
487    *point,
488    *triangle;
489
490  double rtol, atol; 
491  int closed, res;
492
493  PyObject *result;
494     
495  // Convert Python arguments to C
496  if (!PyArg_ParseTuple(args, "OOidd",
497                        &point,
498                        &triangle,
499                        &closed,
500                        &rtol,
501                        &atol)) {
502   
503    PyErr_SetString(PyExc_RuntimeError, 
504                    "_is_inside_triangle could not parse input");
505    return NULL;
506  }
507
508  // Call underlying routine
509  res = __is_inside_triangle((double*) point -> data,
510                             (double*) triangle -> data,
511                             closed,
512                             rtol,
513                             atol);                                                   
514
515
516  // Return result
517  result = Py_BuildValue("i", res);
518  return result; 
519}
520     
521
522
523/*
524PyObject *_intersection(PyObject *self, PyObject *args) {
525  //
526  // intersection(x0, y0, x1, y1)
527  //
528
529  double x, y, x0, y0, x1, y1;
530  int res;
531  PyObject *result;
532
533  // Convert Python arguments to C
534  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) {
535    PyErr_SetString(PyExc_RuntimeError,
536                    "point_on_line could not parse input");   
537    return NULL;
538  }
539
540
541  // Call underlying routine
542  res = __intersection(x, y, x0, y0, x1, y1);
543
544  // Return values a and b
545  result = Py_BuildValue("i", res);
546  return result;
547}
548*/
549
550
551PyObject *_separate_points_by_polygon(PyObject *self, PyObject *args) {
552  //def separate_points_by_polygon(points, polygon, closed, verbose, one_point):
553  //  """Determine whether points are inside or outside a polygon
554  //
555  //  Input:
556  //     point - Tuple of (x, y) coordinates, or list of tuples
557  //     polygon - list of vertices of polygon
558  //     closed - (optional) determine whether points on boundary should be
559  //     regarded as belonging to the polygon (closed = True)
560  //     or not (closed = False)
561
562  //
563  //  Output:
564  //     indices: array of same length as points with indices of points falling
565  //     inside the polygon listed from the beginning and indices of points
566  //     falling outside listed from the end.
567  //     
568  //     count: count of points falling inside the polygon
569  //     
570  //     The indices of points inside are obtained as indices[:count]
571  //     The indices of points outside are obtained as indices[count:]       
572  //
573  //  Examples:
574  //     separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
575  //     will return the indices [0, 2, 1] as only the first and the last point
576  //     is inside the unit square
577  //
578  //  Remarks:
579  //     The vertices may be listed clockwise or counterclockwise and
580  //     the first point may optionally be repeated.
581  //     Polygons do not need to be convex.
582  //     Polygons can have holes in them and points inside a hole is
583  //     regarded as being outside the polygon.
584  //
585  //
586  //  Algorithm is based on work by Darel Finley,
587  //  http://www.alienryderflex.com/polygon/
588  //
589  //
590
591  PyArrayObject
592    *points,
593    *polygon,
594    *indices;
595
596  int closed, verbose; //Flags
597  int count, M, N;
598
599  // Convert Python arguments to C
600  if (!PyArg_ParseTuple(args, "OOOii",
601                        &points,
602                        &polygon,
603                        &indices,
604                        &closed,
605                        &verbose)) {
606   
607
608    PyErr_SetString(PyExc_RuntimeError, 
609                    "separate_points_by_polygon could not parse input");
610    return NULL;
611  }
612
613  M = points -> dimensions[0];   //Number of points
614  N = polygon -> dimensions[0];  //Number of vertices in polygon
615
616  //FIXME (Ole): This isn't send to Python's sys.stdout
617  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
618 
619  //Call underlying routine
620  count = __separate_points_by_polygon(M, N,
621                                       (double*) points -> data,
622                                       (double*) polygon -> data,
623                                       (long*) indices -> data,
624                                       closed, verbose);
625 
626  //NOTE: return number of points inside..
627  return Py_BuildValue("i", count);
628}
629
630
631
632// Method table for python module
633static struct PyMethodDef MethodTable[] = {
634  /* The cast of the function is necessary since PyCFunction values
635   * only take two PyObject* parameters, and rotate() takes
636   * three.
637   */
638
639  {"_point_on_line", _point_on_line, METH_VARARGS, "Print out"},
640  //{"_intersection", _intersection, METH_VARARGS, "Print out"}, 
641  {"_separate_points_by_polygon", _separate_points_by_polygon, 
642                                 METH_VARARGS, "Print out"},
643  {"_interpolate_polyline", _interpolate_polyline, 
644                                 METH_VARARGS, "Print out"},                             
645  {"_is_inside_triangle", _is_inside_triangle, 
646                                 METH_VARARGS, "Print out"},
647  {NULL, NULL, 0, NULL}   /* sentinel */
648};
649
650
651
652// Module initialisation
653void initpolygon_ext(void){
654  Py_InitModule("polygon_ext", MethodTable);
655
656  import_array();     //Necessary for handling of NumPY structures
657}
658
659
660
Note: See TracBrowser for help on using the repository browser.