Changeset 655

Timestamp:

Dec 2, 2004, 10:58:26 AM (20 years ago)

Author:

ole

Message:

Added inside_polynomial and point_on_line with tests

Location:

inundation/ga/storm_surge/pyvolution

Files:

• mesh.py (modified) (2 diffs)
• test_util.py (modified) (1 diff)
• util.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -130 +130 @@
         self.build_tagged_elements_dictionary(tagged_elements)
         
-        #Update boundary indices 
+        #Update boundary indices FIXME: OBSOLETE
         #self.build_boundary_structure()        
 
@@ -313 +313 @@
         """
 
-        #FIXME: Maybe obsolete
+        #FIXME: Now Obsolete - maybe use some comments from here in 
+        #domain.set_boundary
 
         if self.boundary is None:

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -332 +332 @@
         assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
         
+    def test_point_on_line(self):
+                
+        #Endpoints first        
+        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)] )              
+        
+        #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)] ) 
+        
+        
+        
+    def test_inside_polygon(self):
+    
+        #Simplest case: Polygon is the unit square
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+ 
+        x = y = 0.5
+        assert inside_polygon( (x, y), polygon )                
+        
+        y = 1.5
+        assert not inside_polygon( (x, y), polygon )                            
+        #Try point on borders   
+        assert inside_polygon( (1., 0.5), polygon, closed=True) 
+        assert inside_polygon( (0.5, 1), polygon, closed=True)  
+        assert inside_polygon( (0., 0.5), polygon, closed=True) 
+        assert inside_polygon( (0.5, 0.), polygon, closed=True) 
+        
+        assert not inside_polygon( (0.5, 1), polygon, closed=False)     
+        assert not inside_polygon( (0., 0.5), polygon, closed=False)    
+        assert not inside_polygon( (0.5, 0.), polygon, closed=False)    
+        assert not inside_polygon( (1., 0.5), polygon, closed=False)    
         
               
+        #More convoluted and non convex polygon
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+        assert inside_polygon( (0.5, 0.5), polygon )
+        assert inside_polygon( (1, -0.5), polygon ) 
+        assert inside_polygon( (1.5, 0), polygon )
+        
+        assert not inside_polygon( (0.5, 1.5), polygon )        
+        assert not inside_polygon( (0.5, -0.5), polygon )               
+        
+        
+        
+        #Now try the vector formulation returning indices
+        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+        res = inside_polygon( points, polygon )
+        
+        assert allclose( res, [0,1,2] )
+        
+        
+
+                      
 #-------------------------------------------------------------
 if __name__ == "__main__":

inundation/ga/storm_surge/pyvolution/util.py

@@ -47 +47 @@
 
 
+
+def point_on_line(point, line):
+    """Determine whether a point is on a line segment 
+   
+    Input:
+        point: (x, y)
+        line: [ (x0, y0), (x1, y1) ]
+        
+    """
+    
+    #FIXME: Could do with some optimisation         
+        
+    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])                 
+    b = array([x1 - x0, y1 - y0])
+    
+    len_a = sqrt(sum(a**2))            
+    len_b = sqrt(sum(b**2))                 
+                
+    if len_a == 0 or len_a == len_b:
+       return True
+    else:
+       x = dot(a, b)/len_b       
+       return allclose(x, len_a)
+
+        
+    
+def inside_polygon(point, polygon, closed = True):
+    """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
+    
+    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:
+       <Copy from source Franklins code >    
+               
+               
+       points at a vertex are considered to be inside          
+       
+    Algorithm due to Darel Finley, http://www.alienryderflex.com/polygon/ 
+    
+    
+    FIXME: More comments about algoritm
+    FIXME: Deal with case where point on border of polygon (closed)
+       
+    """    
+
+    from Numeric import array, Float, reshape
+    
+    
+    #Input checks
+    try:
+        point = array(point).astype(Float)         
+    except:
+        msg = 'Point %s could not be converted to Numeric array' %point
+        raise msg        
+        
+    try:
+        polygon = array(polygon).astype(Float)             
+    except:
+        msg = 'Polygon %s could not be converted to Numeric array' %polygon
+        raise msg               
+    
+    assert len(polygon.shape) == 2,\
+       'Polygon array must be a 2d array of vertices: %s' %polygon
+
+        
+    assert polygon.shape[1] == 2,\
+       'Polygon array must have two columns: %s' %polygon    
+
+
+    if len(point.shape) == 1:
+        one_point = True
+        points = reshape(point, (1,2))
+    else:       
+        one_point = False
+        points = point
+
+    N = polygon.shape[0] #Number of vertices in polygon 
+
+    px = polygon[:,0]
+    py = polygon[:,1]   
+                
+    #Begin algorithm
+    indices = []
+    for k in range(points.shape[0]):
+        #for each point   
+        x = points[k, 0]
+        y = points[k, 1]        
+
+           
+        inside = False
+        for i in range(N):
+            j = (i+1)%N
+            
+            #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:
+                    inside = True
+                else:        
+                    inside = False
+                break
+        
+            #Check truly inside polygon             
+            if px[i] < y and py[j] >= y or\
+               py[j] < y and py[i] >= y:
+                if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
+                    inside = not inside
+                    
+        if inside: indices.append(k)
+    
+    if one_point: 
+        return inside
+    else:
+        return indices
+             
 
 class File_function: