# source:anuga_core/source/anuga/utilities/polygon_ext.c@6544

Last change on this file since 6544 was 6544, checked in by ole, 14 years ago

Optimised fitting a bit more and cleaned up in search_functions.

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