Changeset 5225 for anuga_core/source/anuga/utilities/polygon.py

Timestamp:

Apr 21, 2008, 5:20:32 PM (15 years ago)

Author:

ole

Message:

Work done during Water Down Under 2008.
Line Mesh intersections towards ticket:175 (flow through a cross section).

File:

• anuga_core/source/anuga/utilities/polygon.py (modified) (4 diffs)

anuga_core/source/anuga/utilities/polygon.py

@@ -10 +10 @@
 #    #print 'Could not find scipy - using Numeric'
 
-from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot
+from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot, allclose
 
 
@@ -17 +17 @@
 from anuga.geospatial_data.geospatial_data import ensure_absolute
 
-def point_on_line(x, y, x0, y0, x1, y1):
+def point_on_line(point, line):
     """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)
-
-    """
-
-    a = array([x - x0, y - y0])
-    a_normal = array([a[1], -a[0]])
-
-    b = array([x1 - x0, y1 - y0])
-
-    if dot(a_normal, b) == 0:
-        #Point is somewhere on the infinite extension of the line
-
-        len_a = sqrt(sum(a**2))
-        len_b = sqrt(sum(b**2))
-        if dot(a, b) >= 0 and len_a <= len_b:
-           return True
-        else:
-           return False
+    Input: 
+        point is given by [x, y]
+        line is given by [x0, y0], [x1, y1]] or
+        the equivalent 2x2 Numeric array with each row corresponding to a point.
+
+    Output:
+        
+    """
+
+    point = ensure_numeric(point)
+    line = ensure_numeric(line)
+
+
+    res = _point_on_line(point[0], point[1],
+                         line[0,0], line[0,1],
+                         line[1,0], line[1,1])
+    
+    return bool(res)
+
+
+
+
+
+def intersection(line0, line1):
+    """Returns intersecting point between two line segments or None (if parallel or no intersection is found)
+
+    Inputs:
+        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
+                      A line can also be a 2x2 numeric array with each row corresponding to a point.
+
+
+    Output:
+        Point [x,y] or None.
+
+    If line extensions intersect outside their limits, None is returned as well.
+    
+    """
+
+    # FIXME (Ole): Write this in C
+
+    line0 = ensure_numeric(line0, Float)
+    line1 = ensure_numeric(line1, Float)    
+
+    x0 = line0[0,0]; y0 = line0[0,1]
+    x1 = line0[1,0]; y1 = line0[1,1]
+
+    x2 = line1[0,0]; y2 = line1[0,1]
+    x3 = line1[1,0]; y3 = line1[1,1]
+
+    denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
+
+    #print 'denom', denom, line0, line1
+    if allclose(denom, 0.0):
+        # Lines are parallel
+        return None
     else:
-      return False
+        u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
+        u0 = u0/denom
+
+        u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
+        u1 = u1/denom        
+
+        x = x0 + u0*(x1-x0)
+        y = y0 + u0*(y1-y0)
+
+        assert allclose(x, x2 + u1*(x3-x2))
+        assert allclose(y, y2 + u1*(y3-y2))        
+
+        # Check if point found lies within given line segments
+        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 
+            # We have intersection
+
+            # Need tolerances if going ahead with this check
+            #msg = 'Algorithm error. Intersection was detected but point (%f, %f) does not lie ' %(x,y)
+            #msg += 'on line0: %s' %(line0)
+            #assert point_on_line([x, y], line0), msg
+            #msg = 'Algorithm error. Intersection was detected but point (%f, %f) does not lie ' %(x,y)
+            #msg += 'on line1: %s' %(line1)            
+            #assert point_on_line([x, y], line1), msg
+        
+            return [x, y]
+        else:
+            return None
 
 
@@ -305 +366 @@
     indices = zeros( M, Int )
 
-    from polygon_ext import separate_points_by_polygon as sep_points
-    count = sep_points(points, polygon, indices,
-                       int(closed), int(verbose))
+    count = _separate_points_by_polygon(points, polygon, indices,
+                                        int(closed), int(verbose))
 
     if verbose: print 'Found %d points (out of %d) inside polygon'\
@@ -780 +840 @@
 from anuga.utilities.compile import can_use_C_extension
 if can_use_C_extension('polygon_ext.c'):
-    # Underlying C implementations can be accessed 
-
-    from polygon_ext import point_on_line
+    # Underlying C implementations can be accessed
+    from polygon_ext import _point_on_line
+    from polygon_ext import _separate_points_by_polygon
+
 else:
     msg = 'C implementations could not be accessed by %s.\n ' %__file__