source: branches/numpy/anuga/utilities/polygon_ext.c @ 7091

Last change on this file since 7091 was 6904, checked in by rwilson, 15 years ago

Finished back-merge from Numeric trunk to numpy branch (and some code the other way).

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