# source:trunk/anuga_core/source/anuga/utilities/polygon_ext.c@8029

Last change on this file since 8029 was 8029, checked in by habili, 13 years ago

modified _polygon_triangle_overlap so that it takes in any polygon.

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