Changeset 6535

Timestamp:

Mar 17, 2009, 11:06:55 PM (15 years ago)

Author:

ole

Message:

Wrote is_inside_triangle in C and tested it some more.

Location:

anuga_core/source/anuga/utilities

Files:

• polygon.py (modified) (5 diffs)
• polygon_ext.c (modified) (9 diffs)
• test_polygon.py (modified) (1 diff)

anuga_core/source/anuga/utilities/polygon.py

@@ -188 +188 @@
 def is_inside_triangle(point, triangle, 
                        closed=True, 
-                       rtol=1.0e-6,
-                       atol=1.0e-8,                      
+                       rtol=1.0e-12,
+                       atol=1.0e-12,                      
                        check_inputs=True, 
                        verbose=False):
@@ -228 +228 @@
     
     rtol and atol will determine how close the point has to be to the edge
-    before it is deemed on the edge.
+    before it is deemed to be on the edge.
     
     """
@@ -246 +246 @@
     if point[1] > max(triangle[:,1]): return False        
 
-    
+    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
+    
+    
+
+    # FIXME (Ole): The rest of this function has been made 
+    # obsolete by the C extension
+        
     # Start search    
     A = triangle[0, :]
@@ -282 +288 @@
                                  Y[0], Y[1],
                                  rtol, atol)
+            
             if res:
                 return True
                 
-                
+    return False
 
 
@@ -1233 +1240 @@
     from polygon_ext import _separate_points_by_polygon
     from polygon_ext import _interpolate_polyline    
+    from polygon_ext import _is_inside_triangle        
     #from polygon_ext import _intersection
 

anuga_core/source/anuga/utilities/polygon_ext.c

@@ -244 +244 @@
 
 
+
+int __is_inside_triangle(double* point,
+                         double* triangle,
+                         int closed,
+                         double rtol,
+                         double atol) {
+                         
+  double vx, vy, v0x, v0y, v1x, v1y;
+  double a00, a10, a01, a11, b0, b1;
+  double denom, alpha, beta;
+  int i, j, res;
+
+  // v0 = C-A 
+  v0x = triangle[4]-triangle[0]; 
+  v0y = triangle[5]-triangle[1];
+  
+  // v1 = B-A   
+  v1x = triangle[2]-triangle[0]; 
+  v1y = triangle[3]-triangle[1];
+
+  // First check if point lies wholly inside triangle
+  a00 = v0x*v0x + v0y*v0y; // innerproduct(v0, v0)
+  a01 = v0x*v1x + v0y*v1y; // innerproduct(v0, v1)
+  a10 = a01;               // innerproduct(v1, v0)
+  a11 = v1x*v1x + v1y*v1y; // innerproduct(v1, v1)
+    
+  denom = a11*a00 - a01*a10;
+
+  if (fabs(denom) > 0.0) {
+    // v = point-A  
+    vx = point[0] - triangle[0]; 
+    vy = point[1] - triangle[1];     
+    
+    b0 = v0x*vx + v0y*vy; // innerproduct(v0, v)        
+    b1 = v1x*vx + v1y*vy; // innerproduct(v1, v)            
+    
+    alpha = (b0*a11 - b1*a01)/denom;
+    beta = (b1*a00 - b0*a10)/denom;        
+    
+    if ((alpha > 0.0) && (beta > 0.0) && (alpha+beta < 1.0)) return 1;
+  }
+
+  if (closed) {
+    // Check if point lies on one of the edges
+        
+    for (i=0; i<3; i++) {
+      j = (i+1) % 3; // Circular index into triangle array
+      res = __point_on_line(point[0], point[1],
+                            triangle[2*i], triangle[2*i+1], 
+                            triangle[2*j], triangle[2*j+1],                         
+                            rtol, atol);
+      if (res) return 1;
+    }
+  }
+                
+  // Default return if point is outside triangle                         
+  return 0;                                              
+}                                                                              
+
+
 int __separate_points_by_polygon(int M,     // Number of points
                                  int N,     // Number of polygon vertices
@@ -255 +315 @@
   int i, j, k, outside_index, inside_index, inside;
 
-  //Find min and max of poly used for optimisation when points
-  //are far away from polygon
-  
-  //FIXME(Ole): Pass in rtol and atol from Python
+  // Find min and max of poly used for optimisation when points
+  // are far away from polygon
+  
+  // FIXME(Ole): Pass in rtol and atol from Python
 
   minpx = polygon[0]; maxpx = minpx;
@@ -273 +333 @@
   }
 
-  //Begin main loop (for each point)
-  inside_index = 0;    //Keep track of points inside
-  outside_index = M-1; //Keep track of points outside (starting from end)   
+  // Begin main loop (for each point)
+  inside_index = 0;    // Keep track of points inside
+  outside_index = M-1; // Keep track of points outside (starting from end)   
   if (verbose){
      printf("Separating %d points\n", M);
@@ -289 +349 @@
     inside = 0;
 
-    //Optimisation
+    // Optimisation
     if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) {
-      //Nothing
+      // Nothing
     } else {   
-      //Check polygon
+      // Check polygon
       for (i=0; i<N; i++) {
-        //printf("k,i=%d,%d\n", k, i);
         j = (i+1)%N;
 
@@ -303 +362 @@
         py_j = polygon[2*j+1];
 
-        //Check for case where point is contained in line segment
+        // Check for case where point is contained in line segment
         if (__point_on_line(x, y, px_i, py_i, px_j, py_j, rtol, atol)) {
           if (closed == 1) {
@@ -312 +371 @@
           break;
         } else {
-          //Check if truly inside polygon
+          // Check if truly inside polygon
           if ( ((py_i < y) && (py_j >= y)) ||
                ((py_j < y) && (py_i >= y)) ) {
@@ -356 +415 @@
   res = __point_on_line(x, y, x0, y0, x1, y1, rtol, atol);
 
-  // Return values a and b
+  // Return result
   result = Py_BuildValue("i", res);
   return result;
@@ -417 +476 @@
 
 
+
+     
+PyObject *_is_inside_triangle(PyObject *self, PyObject *args) {
+  //
+  // _is_inside_triangle(point, triangle, int(closed), rtol, atol)
+  //
+
+  
+  PyArrayObject
+    *point,
+    *triangle;
+
+  double rtol, atol;  
+  int closed, res;
+
+  PyObject *result;
+      
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOidd",
+                        &point,
+                        &triangle,
+                        &closed,
+                        &rtol,
+                        &atol)) {
+    
+    PyErr_SetString(PyExc_RuntimeError, 
+                    "_is_inside_triangle could not parse input");
+    return NULL;
+  }
+
+  // Call underlying routine
+  res = __is_inside_triangle((double*) point -> data,
+                             (double*) triangle -> data,
+                             closed,
+                             rtol,
+                             atol);                                                    
+
+
+  // Return result
+  result = Py_BuildValue("i", res);
+  return result;  
+}
+     
+
+
 /*
 PyObject *_intersection(PyObject *self, PyObject *args) {
@@ -539 +643 @@
   {"_interpolate_polyline", _interpolate_polyline, 
                                  METH_VARARGS, "Print out"},                             
+  {"_is_inside_triangle", _is_inside_triangle, 
+                                 METH_VARARGS, "Print out"},
   {NULL, NULL, 0, NULL}   /* sentinel */
 };

anuga_core/source/anuga/utilities/test_polygon.py

@@ -1857 +1857 @@
         assert num.allclose(z, z_ref)                
         
+
+    def test_is_inside_triangle_more(self):        
+        
+        res = is_inside_triangle([0.5, 0.5], [[ 0.5,  0. ],
+                                              [ 0.5,  0.5],
+                                              [ 0.,   0. ]])
+        assert res is True
+
+        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
+                                 [[ 0.5,  0. ], [ 0.5, 0.5], [0., 0.]])
+        assert res is False
+        
+        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
+                                 [[1., 0.], [1., 0.5], [0.5, 0.]])
+        assert res is False                                 
+                                 
+        
+        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
+                                 [[0.5, 0.5], [0.5, 0.], [1., 0.5]])
+        assert res is True                                 
+
+                
+        res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
+                                 [[0.5, 0.], [0.5, 0.5], [0., 0.]])
+        assert res is False
+
+        
+        res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
+                                 [[0., 0.5], [0., 0.], [0.5, 0.5]])
+        assert res is True                                                                  
+                                 
+        res = is_inside_triangle([0.69999999999999996, 0.69999999999999996],
+                                 [[0.5, 0.], [0.5, 0.5], [0., 0.]])
+        assert res is False                                 
+
+        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
+                                 [[0.25, 0.5], [0.25, 0.25], [0.5, 0.5]])
+        assert res is False
+