source: inundation/ga/storm_surge/pyvolution/util_ext.c @ 999

// Python - C extension for finite_volumes util module.
//
// To compile (Python2.3):
//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
//  gcc -shared util_ext.o  -o util_ext.so      
//
// See the module util.py
//
//
// Ole Nielsen, GA 2004
        
#include "Python.h"
#include "Numeric/arrayobject.h"
#include "math.h"

//Shared snippets
#include "util_ext.h"



int _point_on_line(double x, double y, 
                   double x0, double y0, 
                   double x1, double y1) {
  /*Determine whether a point is on a line segment 
   
    Input: x, y, x0, x0, x1, y1: where
        point is given by x, y
        line is given by (x0, y0) and (x1, y1)
        
  */
    
  double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
  
  a0 = x - x0;
  a1 = y - y0; 
  
  a_normal0 = a1;
  a_normal1 = -a0;
                        
  b0 = x1 - x0;
  b1 = y1 - y0;
    
  if ( a_normal0*b0 + a_normal1*b1 == 0 ) {   
    //Point is somewhere on the infinite extension of the line

    len_a = sqrt(a0*a0 + a1*a1);               
    len_b = sqrt(b0*b0 + b1*b1);
                                        
    if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
      return 1;
    } else {
      return 0;
    }
  } else {
    return 0;
  }
}



int _inside_polygon(int M,           // Number of points
                    int N,           // Number of polygon vertices
                    double* points,
                    double* polygon,
                    int* indices,    // M-Array for storage indices 
                    int closed, 
                    int verbose) {
  
  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
  int i, j, k, count, inside;
  
  //Find min and max of poly used for optimisation when points 
  //are far away from polygon 
  
  
  minpx = polygon[0]; maxpx = minpx;
  minpy = polygon[1]; maxpy = minpy;  
  
  for (i=0; i<N; i++) {
    px_i = polygon[2*i];
    py_i = polygon[2*i + 1];
  
    if (px_i < minpx) minpx = px_i;
    if (px_i > maxpx) maxpx = px_i;    
    if (py_i < minpy) minpy = py_i;
    if (py_i > maxpy) maxpy = py_i;        
  }
  
  //Begin main loop (for each point)
  count = 0;
  for (k=0; k<M; k++) {
    //FIXME: Do this later
    //if (verbose){
    //  if (k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
        
    x = points[2*k];
    y = points[2*k + 1];        

    inside = 0;

    //Optimisation
    if ((x > maxpx) || (x < minpx)) continue;
    if ((y > maxpy) || (y < minpy)) continue;        

    //Check polygon 
    for (i=0; i<N; i++) {
      //printf("k,i=%d,%d\n", k, i);
      j = (i+1)%N;
      
      px_i = polygon[2*i];
      py_i = polygon[2*i+1];
      px_j = polygon[2*j];
      py_j = polygon[2*j+1];      
      
      //Check for case where point is contained in line segment 
      if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
        if (closed == 1) {
          inside = 1;
        } else {
          inside = 0;
        }
        break;  
      } else {   
        //Check if truly inside polygon             
        if ( ((py_i < y) && (py_j >= y)) ||
             ((py_j < y) && (py_i >= y)) ) {
          if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x)
            inside = 1-inside;
        }
      }
    }
    if (inside == 1) {
      indices[count] = k;
      count++;
    }       
  } // End k
  
  return count;
}



// Gateways to Python
PyObject *point_on_line(PyObject *self, PyObject *args) {
  //
  // point_on_line(x, y, x0, y0, x1, y1)
  //
  
  double x, y, x0, y0, x1, y1;
  int res;
  PyObject *result;
  
  // Convert Python arguments to C  
  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1))  
    return NULL;

  
  // Call underlying routine
  res = _point_on_line(x, y, x0, y0, x1, y1);

  // Return values a and b
  result = Py_BuildValue("i", res);
  return result;
}



PyObject *inside_polygon(PyObject *self, PyObject *args) {
  //def inside_polygon(point, polygon, closed, verbose, one_point):  
  //  """Determine whether points are inside or outside a polygon
  //  
  //  Input:
  //     point - Tuple of (x, y) coordinates, or list of tuples
  //     polygon - list of vertices of polygon
  //     one_poin - If True Boolean value should be returned 
  //                If False, indices of points inside returned
  //     closed - (optional) determine whether points on boundary should be 
  //     regarded as belonging to the polygon (closed = True) 
  //     or not (closed = False) 
  
  //  
  //  Output:
  //     If one point is considered, True or False is returned.
  //     If multiple points are passed in, the function returns indices 
  //     of those points that fall inside the polygon     
  //     
  //  Examples:
  //     inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] )
  //     will evaluate to True as the point 0.5, 0.5 is inside the unit square 
  //     
  //     inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
  //     will return the indices [0, 2] as only the first and the last point 
  //     is inside the unit square
  //     
  //  Remarks:
  //     The vertices may be listed clockwise or counterclockwise and
  //     the first point may optionally be repeated.    
  //     Polygons do not need to be convex.
  //     Polygons can have holes in them and points inside a hole is 
  //     regarded as being outside the polygon.                
  //           
  //     
  //  Algorithm is based on work by Darel Finley, 
  //  http://www.alienryderflex.com/polygon/ 
  //
  //
  
  PyArrayObject 
    *point, 
    *polygon, 
    *indices;
    
  int closed, verbose; //Flags     
  int count, M, N;
  
  // Convert Python arguments to C  
  if (!PyArg_ParseTuple(args, "OOOii", 
                        &point, 
                        &polygon, 
                        &indices, 
                        &closed, 
                        &verbose))  
    return NULL;

  M = point -> dimensions[0];    //Number of points          
  N = polygon -> dimensions[0];  //Number of vertices in polygon        

  //printf("M=%d, N=%d\n", M, N);
  // Call underlying routine
  count = _inside_polygon(M, N, 
                          (double*) point -> data,
                          (double*) polygon -> data,
                          (int*) indices -> data,
                          closed, verbose); 
                           
  //NOTE: return number of points inside..
  //printf("COunt=%d\n", count);
  return Py_BuildValue("i", count);
}



PyObject *gradient(PyObject *self, PyObject *args) {
  //
  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
  //
  
  double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b;
  PyObject *result;
  
  // Convert Python arguments to C  
  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
                        &q0, &q1, &q2))  
    return NULL;

  
  // Call underlying routine
  _gradient(x0, y0, x1, y1, x2, y2, 
            q0, q1, q2, &a, &b);

  // Return values a and b
  result = Py_BuildValue("dd", a, b);
  return result;
}


// Method table for python module
static struct PyMethodDef MethodTable[] = {
  /* The cast of the function is necessary since PyCFunction values
   * only take two PyObject* parameters, and rotate() takes
   * three.
   */
  
  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},  
  {"gradient", gradient, METH_VARARGS, "Print out"},    
  {"point_on_line", point_on_line, METH_VARARGS, "Print out"},      
  {"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},      
  {NULL, NULL, 0, NULL}   /* sentinel */
};



// Module initialisation   
void initutil_ext(void){
  Py_InitModule("util_ext", MethodTable);
  
  import_array();     //Necessary for handling of NumPY structures  
}