Changeset 671

Timestamp:

Dec 5, 2004, 8:21:43 PM (20 years ago)

Author:

ole

Message:

Implemented point_on_line in C.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• test_util.py (modified) (1 diff)
• util.py (modified) (4 diffs)
• util_ext.c (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -391 +391 @@
                 
         #Endpoints first        
-        assert point_on_line( (0, 0), [(0,0), (1,0)] )  
-        assert point_on_line( (1, 0), [(0,0), (1,0)] )          
+        assert point_on_line( 0, 0, 0,0, 1,0 )  
+        assert point_on_line( 1, 0, 0,0, 1,0 )          
 
         #Then points on line    
-        assert point_on_line( (0.5, 0), [(0,0), (1,0)] ) 
-        assert point_on_line( (0, 0.5), [(0,1), (0,0)] )                
-        assert point_on_line( (1, 0.5), [(1,1), (1,0)] )        
-        assert point_on_line( (0.5, 0.5), [(0,0), (1,1)] )              
+        assert point_on_line( 0.5, 0, 0,0, 1,0 ) 
+        assert point_on_line( 0, 0.5, 0,1, 0,0 )                
+        assert point_on_line( 1, 0.5, 1,1, 1,0 )        
+        assert point_on_line( 0.5, 0.5, 0,0, 1,1 )              
         
         #Then points not on line                
-        assert not point_on_line( (0.5, 0), [(0,1), (1,1)] )    
-        assert not point_on_line( (0, 0.5), [(0,0), (1,1)] ) 
+        assert not point_on_line( 0.5, 0, 0,1, 1,1 )    
+        assert not point_on_line( 0, 0.5, 0,0, 1,1 ) 
         
         #From real example that failed
-        assert not point_on_line( (40,50), [(40,20), (40,40)] )
+        assert not point_on_line( 40,50, 40,20, 40,40 )
         
         
         #From real example that failed
-        assert not point_on_line( (40,19), [(40,20), (40,40)] ) 
+        assert not point_on_line( 40,19, 40,20, 40,40 ) 
                         
         

inundation/ga/storm_surge/pyvolution/util.py

@@ -48 +48 @@
 
 
-def point_on_line(point, line):
+#def point_on_line(point, line):
+def point_on_line(x, y, x0, y0, x1, y1): 
     """Determine whether a point is on a line segment 
    
-    Input:
-        point: (x, y)
-        line: [ (x0, y0), (x1, y1) ]
+    Input: x, y, x0, x0, x1, y1: where
+        point is given by x, y
+        line is given by (x0, y0) and (x1, y1)
         
     """
@@ -61 +62 @@
     from Numeric import array, dot, allclose
     from math import sqrt
-
-    x = point[0]
-    y = point[1]    
-    
-    x0 = line[0][0]                     
-    x1 = line[1][0]
-    y0 = line[0][1]        
-    y1 = line[1][1]
 
     a = array([x - x0, y - y0]) 
@@ -173 +166 @@
             
             #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 point_on_line( (x,y), [ [px[i], py[i]], [px[j], py[j]] ]):
+            if point_on_line(x, y, px[i], py[i], px[j], py[j]):
                 if closed:
                     inside = True
@@ -515 +509 @@
 import compile
 if compile.can_use_C_extension('util_ext.c'):
-    from util_ext import gradient
+    from util_ext import gradient, point_on_line 
 else:
     gradient = gradient_python

inundation/ga/storm_surge/pyvolution/util_ext.c

@@ -18 +18 @@
 
 
-// Gateway to Python
+
+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;
+  }
+}
+
+
+
+// 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 *gradient(PyObject *self, PyObject *args) {
   //
@@ -52 +116 @@
   //{"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"},      
   {NULL, NULL, 0, NULL}   /* sentinel */
 };