Ignore:
Timestamp:
Mar 19, 2009, 1:43:34 PM (15 years ago)
Author:
rwilson
Message:

Merged trunk into numpy, all tests and validations work.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/numpy/anuga/utilities/polygon_ext.c

    r6410 r6553  
    247247
    248248
     249int __is_inside_triangle(double* point,
     250                         double* triangle,
     251                         int closed,
     252                         double rtol,
     253                         double atol) {
     254                         
     255  double vx, vy, v0x, v0y, v1x, v1y;
     256  double a00, a10, a01, a11, b0, b1;
     257  double denom, alpha, beta;
     258 
     259  double x, y; // Point coordinates
     260  int i, j, res;
     261
     262  x = point[0];
     263  y = point[1];
     264 
     265  // Quickly reject points that are clearly outside
     266  if ((x < triangle[0]) &&
     267      (x < triangle[2]) &&
     268      (x < triangle[4])) return 0;       
     269     
     270  if ((x > triangle[0]) &&
     271      (x > triangle[2]) &&
     272      (x > triangle[4])) return 0;             
     273 
     274  if ((y < triangle[1]) &&
     275      (y < triangle[3]) &&
     276      (y < triangle[5])) return 0;       
     277     
     278  if ((y > triangle[1]) &&
     279      (y > triangle[3]) &&
     280      (y > triangle[5])) return 0;             
     281 
     282 
     283  // v0 = C-A
     284  v0x = triangle[4]-triangle[0];
     285  v0y = triangle[5]-triangle[1];
     286 
     287  // v1 = B-A   
     288  v1x = triangle[2]-triangle[0];
     289  v1y = triangle[3]-triangle[1];
     290
     291  // First check if point lies wholly inside triangle
     292  a00 = v0x*v0x + v0y*v0y; // innerproduct(v0, v0)
     293  a01 = v0x*v1x + v0y*v1y; // innerproduct(v0, v1)
     294  a10 = a01;               // innerproduct(v1, v0)
     295  a11 = v1x*v1x + v1y*v1y; // innerproduct(v1, v1)
     296   
     297  denom = a11*a00 - a01*a10;
     298
     299  if (fabs(denom) > 0.0) {
     300    // v = point-A 
     301    vx = x - triangle[0];
     302    vy = y - triangle[1];     
     303   
     304    b0 = v0x*vx + v0y*vy; // innerproduct(v0, v)       
     305    b1 = v1x*vx + v1y*vy; // innerproduct(v1, v)           
     306   
     307    alpha = (b0*a11 - b1*a01)/denom;
     308    beta = (b1*a00 - b0*a10)/denom;       
     309   
     310    if ((alpha > 0.0) && (beta > 0.0) && (alpha+beta < 1.0)) return 1;
     311  }
     312
     313  if (closed) {
     314    // Check if point lies on one of the edges
     315       
     316    for (i=0; i<3; i++) {
     317      j = (i+1) % 3; // Circular index into triangle array
     318      res = __point_on_line(x, y,
     319                            triangle[2*i], triangle[2*i+1],
     320                            triangle[2*j], triangle[2*j+1],                         
     321                            rtol, atol);
     322      if (res) return 1;
     323    }
     324  }
     325               
     326  // Default return if point is outside triangle                         
     327  return 0;                                             
     328}                                                                             
     329
     330
    249331int __separate_points_by_polygon(int M,     // Number of points
    250332                                int N,     // Number of polygon vertices
     
    427509
    428510
     511     
     512PyObject *_is_inside_triangle(PyObject *self, PyObject *args) {
     513  //
     514  // _is_inside_triangle(point, triangle, int(closed), rtol, atol)
     515  //
     516
     517 
     518  PyArrayObject
     519    *point,
     520    *triangle;
     521
     522  double rtol, atol; 
     523  int closed, res;
     524
     525  PyObject *result;
     526     
     527  // Convert Python arguments to C
     528  if (!PyArg_ParseTuple(args, "OOidd",
     529                        &point,
     530                        &triangle,
     531                        &closed,
     532                        &rtol,
     533                        &atol)) {
     534   
     535    PyErr_SetString(PyExc_RuntimeError,
     536                    "_is_inside_triangle could not parse input");
     537    return NULL;
     538  }
     539
     540  // Call underlying routine
     541  res = __is_inside_triangle((double*) point -> data,
     542                             (double*) triangle -> data,
     543                             closed,
     544                             rtol,
     545                             atol);                                                   
     546
     547
     548  // Return result
     549  result = Py_BuildValue("i", res);
     550  return result; 
     551}
     552     
     553
     554
    429555/*
    430556PyObject *_intersection(PyObject *self, PyObject *args) {
     
    554680  {"_interpolate_polyline", _interpolate_polyline,
    555681                                 METH_VARARGS, "Print out"},                             
     682  {"_is_inside_triangle", _is_inside_triangle,
     683                                 METH_VARARGS, "Print out"},
    556684  {NULL, NULL, 0, NULL}   /* sentinel */
    557685};
Note: See TracChangeset for help on using the changeset viewer.