# source:development/pyvolution-1d/util_ext.c@2716

Last change on this file since 2716 was 2716, checked in by jakeman, 18 years ago

File size: 11.9 KB
Line
1// Python - C extension for finite_volumes util module.
2//
3// To compile (Python2.3):
4//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O
5//  gcc -shared util_ext.o  -o util_ext.so
6//
7// See the module util.py
8//
9//
10// Ole Nielsen, GA 2004
11//
12//NOTE: On 64 bit systems use long* instead of int* for Numeric arrays
13//this will also work on 32 bit systems
14
15
16#include "Python.h"
17#include "Numeric/arrayobject.h"
18#include "math.h"
19
20//Shared snippets
21#include "util_ext.h"
22
23
24
25int _point_on_line(double x, double y,
26                   double x0, double y0,
27                   double x1, double y1) {
28  /*Determine whether a point is on a line segment
29
30    Input: x, y, x0, x0, x1, y1: where
31        point is given by x, y
32        line is given by (x0, y0) and (x1, y1)
33
34  */
35
36  double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
37
38  a0 = x - x0;
39  a1 = y - y0;
40
41  a_normal0 = a1;
42  a_normal1 = -a0;
43
44  b0 = x1 - x0;
45  b1 = y1 - y0;
46
47  if ( a_normal0*b0 + a_normal1*b1 == 0 ) {
48    //Point is somewhere on the infinite extension of the line
49
50    len_a = sqrt(a0*a0 + a1*a1);
51    len_b = sqrt(b0*b0 + b1*b1);
52
53    if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
54      return 1;
55    } else {
56      return 0;
57    }
58  } else {
59    return 0;
60  }
61}
62
63
64
65int _separate_points_by_polygon(int M,     // Number of points
66                                int N,     // Number of polygon vertices
67                                double* points,
68                                double* polygon,
69                                long* indices,  // M-Array for storage indices
70                                int closed,
71                                int verbose) {
72
73  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
74  int i, j, k, outside_index, inside_index, inside;
75
76  //Find min and max of poly used for optimisation when points
77  //are far away from polygon
78
79  minpx = polygon[0]; maxpx = minpx;
80  minpy = polygon[1]; maxpy = minpy;
81
82  for (i=0; i<N; i++) {
83    px_i = polygon[2*i];
84    py_i = polygon[2*i + 1];
85
86    if (px_i < minpx) minpx = px_i;
87    if (px_i > maxpx) maxpx = px_i;
88    if (py_i < minpy) minpy = py_i;
89    if (py_i > maxpy) maxpy = py_i;
90  }
91
92  //Begin main loop (for each point)
93  inside_index = 0;    //Keep track of points inside
94  outside_index = M-1; //Keep track of points outside (starting from end)
95  for (k=0; k<M; k++) {
96    if (verbose){
97      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
98    }
99
100    x = points[2*k];
101    y = points[2*k + 1];
102
103    inside = 0;
104
105    //Optimisation
106    if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) {
107      //Nothing
108    } else {
109      //Check polygon
110      for (i=0; i<N; i++) {
111        //printf("k,i=%d,%d\n", k, i);
112        j = (i+1)%N;
113
114        px_i = polygon[2*i];
115        py_i = polygon[2*i+1];
116        px_j = polygon[2*j];
117        py_j = polygon[2*j+1];
118
119        //Check for case where point is contained in line segment
120        if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
121          if (closed == 1) {
122            inside = 1;
123          } else {
124            inside = 0;
125          }
126          break;
127        } else {
128          //Check if truly inside polygon
129          if ( ((py_i < y) && (py_j >= y)) ||
130               ((py_j < y) && (py_i >= y)) ) {
131            if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
132              inside = 1-inside;
133          }
134        }
135      }
136    }
137    if (inside == 1) {
138      indices[inside_index] = k;
139      inside_index += 1;
140    } else {
141      indices[outside_index] = k;
142      outside_index -= 1;
143    }
144  } // End k
145
146  return inside_index;
147}
148
149
150
151// Gateways to Python
152PyObject *point_on_line(PyObject *self, PyObject *args) {
153  //
154  // point_on_line(x, y, x0, y0, x1, y1)
155  //
156
157  double x, y, x0, y0, x1, y1;
158  int res;
159  PyObject *result;
160
161  // Convert Python arguments to C
162  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1))
163    return NULL;
164
165
166  // Call underlying routine
167  res = _point_on_line(x, y, x0, y0, x1, y1);
168
169  // Return values a and b
170  result = Py_BuildValue("i", res);
171  return result;
172}
173
174
175
176PyObject *separate_points_by_polygon(PyObject *self, PyObject *args) {
177  //def separate_points_by_polygon(points, polygon, closed, verbose, one_point):
178  //  """Determine whether points are inside or outside a polygon
179  //
180  //  Input:
181  //     point - Tuple of (x, y) coordinates, or list of tuples
182  //     polygon - list of vertices of polygon
183  //     closed - (optional) determine whether points on boundary should be
184  //     regarded as belonging to the polygon (closed = True)
185  //     or not (closed = False)
186
187  //
188  //  Output:
189  //     indices: array of same length as points with indices of points falling
190  //     inside the polygon listed from the beginning and indices of points
191  //     falling outside listed from the end.
192  //
193  //     count: count of points falling inside the polygon
194  //
195  //     The indices of points inside are obtained as indices[:count]
196  //     The indices of points outside are obtained as indices[count:]
197  //
198  //  Examples:
199  //     separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
200  //     will return the indices [0, 2, 1] as only the first and the last point
201  //     is inside the unit square
202  //
203  //  Remarks:
204  //     The vertices may be listed clockwise or counterclockwise and
205  //     the first point may optionally be repeated.
206  //     Polygons do not need to be convex.
207  //     Polygons can have holes in them and points inside a hole is
208  //     regarded as being outside the polygon.
209  //
210  //
211  //  Algorithm is based on work by Darel Finley,
212  //  http://www.alienryderflex.com/polygon/
213  //
214  //
215
216  PyArrayObject
217    *points,
218    *polygon,
219    *indices;
220
221  int closed, verbose; //Flags
222  int count, M, N;
223
224  // Convert Python arguments to C
225  if (!PyArg_ParseTuple(args, "OOOii",
226                        &points,
227                        &polygon,
228                        &indices,
229                        &closed,
230                        &verbose))
231    return NULL;
232
233  M = points -> dimensions[0];   //Number of points
234  N = polygon -> dimensions[0];  //Number of vertices in polygon
235
236  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
237
238  //Call underlying routine
239  count = _separate_points_by_polygon(M, N,
240                                      (double*) points -> data,
241                                      (double*) polygon -> data,
242                                      (long*) indices -> data,
243                                      closed, verbose);
244
245  //NOTE: return number of points inside..
246  return Py_BuildValue("i", count);
247}
248
249
250
251PyObject *gradient(PyObject *self, PyObject *args) {
252  //
253  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
254  //
255
256  double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b;
257  PyObject *result;
258
259  // Convert Python arguments to C
260  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
261                        &q0, &q1, &q2))
262    return NULL;
263
264
265  // Call underlying routine
266  _gradient(x0, y0, x1, y1, x2, y2,
267            q0, q1, q2, &a, &b);
268
269  // Return values a and b
270  result = Py_BuildValue("dd", a, b);
271  return result;
272}
273
274PyObject *gradient2(PyObject *self, PyObject *args) {
275  //
276  // a,b = gradient2(x0, y0, x1, y1, q0, q1)
277  //
278
279  double x0, y0, x1, y1, q0, q1, a, b;
280  PyObject *result;
281
282  // Convert Python arguments to C
283  if (!PyArg_ParseTuple(args, "dddddd", &x0, &y0, &x1, &y1, &q0, &q1))
284    return NULL;
285
286
287  // Call underlying routine
288  _gradient2(x0, y0, x1, y1, q0, q1, &a, &b);
289
290  // Return values a and b
291  result = Py_BuildValue("dd", a, b);
292  return result;
293}
294
295// Method table for python module
296static struct PyMethodDef MethodTable[] = {
297  /* The cast of the function is necessary since PyCFunction values
298   * only take two PyObject* parameters, and rotate() takes
299   * three.
300   */
301
302  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
305  {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
306  //{"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},
307  {"separate_points_by_polygon", separate_points_by_polygon,
308                                 METH_VARARGS, "Print out"},
309  {NULL, NULL, 0, NULL}   /* sentinel */
310};
311
312
313
314// Module initialisation
315void initutil_ext(void){
316  Py_InitModule("util_ext", MethodTable);
317
318  import_array();     //Necessary for handling of NumPY structures
319}
320
321
322
323
324//OBSOLETE
325
326int _inside_polygon(int M,           // Number of points
327                    int N,           // Number of polygon vertices
328                    double* points,
329                    double* polygon,
330                    long* indices,    // M-Array for storage indices
331                    int closed,
332                    int verbose) {
333
334  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
335  int i, j, k, count, inside;
336
337  //Find min and max of poly used for optimisation when points
338  //are far away from polygon
339
340  minpx = polygon[0]; maxpx = minpx;
341  minpy = polygon[1]; maxpy = minpy;
342
343  for (i=0; i<N; i++) {
344    px_i = polygon[2*i];
345    py_i = polygon[2*i + 1];
346
347    if (px_i < minpx) minpx = px_i;
348    if (px_i > maxpx) maxpx = px_i;
349    if (py_i < minpy) minpy = py_i;
350    if (py_i > maxpy) maxpy = py_i;
351  }
352
353  //Begin main loop (for each point)
354  count = 0;
355  for (k=0; k<M; k++) {
356    if (verbose){
357      if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M);
358    }
359
360    x = points[2*k];
361    y = points[2*k + 1];
362
363    inside = 0;
364
365    //Optimisation
366    if ((x > maxpx) || (x < minpx)) continue;
367    if ((y > maxpy) || (y < minpy)) continue;
368
369    //Check polygon
370    for (i=0; i<N; i++) {
371      //printf("k,i=%d,%d\n", k, i);
372      j = (i+1)%N;
373
374      px_i = polygon[2*i];
375      py_i = polygon[2*i+1];
376      px_j = polygon[2*j];
377      py_j = polygon[2*j+1];
378
379      //Check for case where point is contained in line segment
380      if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
381        if (closed == 1) {
382          inside = 1;
383        } else {
384          inside = 0;
385        }
386        break;
387      } else {
388        //Check if truly inside polygon
389        if ( ((py_i < y) && (py_j >= y)) ||
390             ((py_j < y) && (py_i >= y)) ) {
391          if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
392            inside = 1-inside;
393        }
394      }
395    }
396    if (inside == 1) {
397      indices[count] = k;
398      count++;
399    }
400  } // End k
401
402  return count;
403}
404
405
406
407PyObject *inside_polygon(PyObject *self, PyObject *args) {
408  //def inside_polygon(point, polygon, closed, verbose, one_point):
409  //  """Determine whether points are inside or outside a polygon
410  //
411  //  Input:
412  //     point - Tuple of (x, y) coordinates, or list of tuples
413  //     polygon - list of vertices of polygon
414  //     one_poin - If True Boolean value should be returned
415  //                If False, indices of points inside returned
416  //     closed - (optional) determine whether points on boundary should be
417  //     regarded as belonging to the polygon (closed = True)
418  //     or not (closed = False)
419
420  //
421  //  Output:
422  //     If one point is considered, True or False is returned.
423  //     If multiple points are passed in, the function returns indices
424  //     of those points that fall inside the polygon
425  //
426  //  Examples:
427  //     inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] )
428  //     will evaluate to True as the point 0.5, 0.5 is inside the unit square
429  //
430  //     inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
431  //     will return the indices [0, 2] as only the first and the last point
432  //     is inside the unit square
433  //
434  //  Remarks:
435  //     The vertices may be listed clockwise or counterclockwise and
436  //     the first point may optionally be repeated.
437  //     Polygons do not need to be convex.
438  //     Polygons can have holes in them and points inside a hole is
439  //     regarded as being outside the polygon.
440  //
441  //
442  //  Algorithm is based on work by Darel Finley,
443  //  http://www.alienryderflex.com/polygon/
444  //
445  //
446
447  PyArrayObject
448    *point,
449    *polygon,
450    *indices;
451
452  int closed, verbose; //Flags
453  int count, M, N;
454
455  // Convert Python arguments to C
456  if (!PyArg_ParseTuple(args, "OOOii",
457                        &point,
458                        &polygon,
459                        &indices,
460                        &closed,
461                        &verbose))
462    return NULL;
463
464  M = point -> dimensions[0];    //Number of points
465  N = polygon -> dimensions[0];  //Number of vertices in polygon
466
467  if (verbose) printf("Got %d points and %d polygon vertices\n", M, N);
468
469  //Call underlying routine
470  count = _inside_polygon(M, N,
471                          (double*) point -> data,
472                          (double*) polygon -> data,
473                          (long*) indices -> data,
474                          closed, verbose);
475
476  //NOTE: return number of points inside..
477  //printf("COunt=%d\n", count);
478  return Py_BuildValue("i", count);
479}
Note: See TracBrowser for help on using the repository browser.