source: inundation/pyvolution/util_ext.c @ 1916

Last change on this file since 1916 was 1486, checked in by ole, 19 years ago

Applyed new formula for computing gradients for triangles that have only one neighbour (thanks to Matt for pointing out the bug) and added unit tests as appropriate.
See util_ext.c for documentation of the derivation

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"},
303  {"gradient", gradient, METH_VARARGS, "Print out"},
304  {"gradient2", gradient2, METH_VARARGS, "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.