Changeset 6534

Timestamp:

Mar 17, 2009, 5:44:57 PM (16 years ago)

Author:

ole

Message:

Added new function is_inside_triangle. This helps the fitting functions.

Location:

anuga_core/source/anuga/utilities

Files:

• numerical_tools.py (modified) (1 diff)
• polygon.py (modified) (8 diffs)
• polygon_ext.c (modified) (1 diff)
• test_polygon.py (modified) (6 diffs)

anuga_core/source/anuga/utilities/numerical_tools.py

@@ -84 +84 @@
     p = num.innerproduct(v1, v2)
     c = num.innerproduct(v1, normal_vector(v2)) # Projection onto normal
-                                            # (negative cross product)
+                                                # (negative cross product)
         
     theta = safe_acos(p)

anuga_core/source/anuga/utilities/polygon.py

@@ -186 +186 @@
     return status, value
 
-
-
-
+def is_inside_triangle(point, triangle, 
+                       closed=True, 
+                       rtol=1.0e-6,
+                       atol=1.0e-8,                      
+                       check_inputs=True, 
+                       verbose=False):
+    """Determine if one point is inside a triangle
+    
+    This uses the barycentric method:
+    
+    Triangle is A, B, C
+    Point P can then be written as
+    
+    P = A + alpha * (C-A) + beta * (B-A)
+    or if we let 
+    v=P-A, v0=C-A, v1=B-A    
+    
+    v = alpha*v0 + beta*v1 
+
+    Dot this equation by v0 and v1 to get two:
+    
+    dot(v0, v) = alpha*dot(v0, v0) + beta*dot(v0, v1)
+    dot(v1, v) = alpha*dot(v1, v0) + beta*dot(v1, v1)    
+    
+    or if a_ij = dot(v_i, v_j) and b_i = dot(v_i, v)
+    the matrix equation:
+    
+    a_00 a_01   alpha     b_0
+                       = 
+    a_10 a_11   beta      b_1
+    
+    Solving for alpha and beta yields:
+    
+    alpha = (b_0*a_11 - b_1*a_01)/denom
+    beta =  (b_1*a_00 - b_0*a_10)/denom
+    
+    with denom = a_11*a_00 - a_10*a_01
+    
+    The point is in the triangle whenever
+    alpha and beta and their sums are in the unit interval.
+    
+    rtol and atol will determine how close the point has to be to the edge
+    before it is deemed on the edge.
+    
+    """
+
+    if check_inputs is True:
+        triangle = ensure_numeric(triangle)    
+
+        point = ensure_numeric(point, num.Float) # Convert point to array of points
+        msg = 'is_inside_triangle must be invoked with one point only'
+        assert num.allclose(point.shape, [2]), msg
+    
+        
+    # Quickly reject points that are clearly outside
+    if point[0] < min(triangle[:,0]): return False 
+    if point[0] > max(triangle[:,0]): return False    
+    if point[1] < min(triangle[:,1]): return False
+    if point[1] > max(triangle[:,1]): return False        
+
+    
+    # Start search    
+    A = triangle[0, :]
+    B = triangle[1, :]
+    C = triangle[2, :]
+    
+    # Now check if point lies wholly inside triangle
+    v0 = C-A 
+    v1 = B-A        
+        
+    a00 = num.innerproduct(v0, v0)
+    a10 = a01 = num.innerproduct(v0, v1)
+    a11 = num.innerproduct(v1, v1)
+    
+    denom = a11*a00 - a01*a10
+    
+    if abs(denom) > 0.0:
+        v = point-A
+        b0 = num.innerproduct(v0, v)        
+        b1 = num.innerproduct(v1, v)            
+    
+        alpha = (b0*a11 - b1*a01)/denom
+        beta = (b1*a00 - b0*a10)/denom        
+
+        if (alpha > 0.0) and (beta > 0.0) and (alpha+beta < 1.0):
+            return True 
+
+
+    if closed is True:
+        # Check if point lies on one of the edges
+        
+        for X, Y in [[A,B], [B,C], [C,A]]:
+            res = _point_on_line(point[0], point[1],
+                                 X[0], X[1],
+                                 Y[0], Y[1],
+                                 rtol, atol)
+            if res:
+                return True
+                
+                
+
+
+    
+
+    
+    
+    
+def is_inside_polygon_quick(point, polygon, closed=True, verbose=False):
+    """Determine if one point is inside a polygon
+    Both point and polygon are assumed to be numeric arrays or lists and
+    no georeferencing etc or other checks will take place.
+    
+    As such it is faster than is_inside_polygon 
+    """
+
+    # FIXME(Ole): This function isn't being used
+    polygon = ensure_numeric(polygon, num.Float)
+    points = ensure_numeric(point, num.Float) # Convert point to array of points
+    points = points[num.NewAxis, :]
+    msg = 'is_inside_polygon must be invoked with one point only'
+    msg += '\nI got %s and converted to %s' % (str(point), str(points.shape))
+    assert points.shape[0] == 1 and points.shape[1] == 2, msg
+
+    
+    indices = num.zeros(1, num.Int)
+
+    count = _separate_points_by_polygon(points, polygon, indices,
+                                        int(closed), int(verbose))
+
+    
+    #indices, count = separate_points_by_polygon(points, polygon,
+    #                                            closed=closed,
+    #                                            input_checks=False,
+    #                                            verbose=verbose)
+
+    if count > 0:
+        return True
+                                                
+        
+        
+    
 def is_inside_polygon(point, polygon, closed=True, verbose=False):
     """Determine if one point is inside a polygon
@@ -231 +369 @@
         raise msg
 
+    polygon = ensure_absolute(polygon)        
     try:
         polygon = ensure_absolute(polygon)
@@ -317 +456 @@
        
 
-def in_and_outside_polygon(points, polygon, closed = True, verbose = False):
+def in_and_outside_polygon(points, polygon, closed=True, verbose=False):
     """Determine points inside and outside a polygon
 
@@ -362 +501 @@
 
 def separate_points_by_polygon(points, polygon,
-                               closed = True, verbose = False):
+                               closed=True, 
+                               check_input=True,
+                               verbose=False):
     """Determine whether points are inside or outside a polygon
 
@@ -371 +512 @@
        regarded as belonging to the polygon (closed = True)
        or not (closed = False)
+       check_input: Allows faster execution if set to False
 
     Outputs:
@@ -403 +545 @@
     """
 
-
-    #if verbose: print 'Checking input to separate_points_by_polygon'
-
-
-    #Input checks
-
-    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
-    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
-
-
-    try:
-        points = ensure_numeric(points, num.Float)
-    except NameError, e:
-        raise NameError, e
-    except:
-        msg = 'Points could not be converted to Numeric array'
-        raise msg
-
-    #if verbose: print 'Checking input to separate_points_by_polygon 2'
-    try:
-        polygon = ensure_numeric(polygon, num.Float)
-    except NameError, e:
-        raise NameError, e
-    except:
-        msg = 'Polygon could not be converted to Numeric array'
-        raise msg
-
-    msg = 'Polygon array must be a 2d array of vertices'
-    assert len(polygon.shape) == 2, msg
-
-    msg = 'Polygon array must have two columns' 
-    assert polygon.shape[1] == 2, msg
-
-
-    msg = 'Points array must be 1 or 2 dimensional.'
-    msg += ' I got %d dimensions' %len(points.shape)
-    assert 0 < len(points.shape) < 3, msg
-
-
-    if len(points.shape) == 1:
-        # Only one point was passed in.
-        # Convert to array of points
-        points = num.reshape(points, (1,2))
-
-    
-    msg = 'Point array must have two columns (x,y), '
-    msg += 'I got points.shape[1] == %d' %points.shape[0]
-    assert points.shape[1] == 2, msg
+    if check_input:
+        #Input checks
+
+        assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
+        assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
+
+
+        try:
+            points = ensure_numeric(points, num.Float)
+        except NameError, e:
+            raise NameError, e
+        except:
+            msg = 'Points could not be converted to Numeric array'
+            raise msg
+
+        try:
+            polygon = ensure_numeric(polygon, num.Float)
+        except NameError, e:
+            raise NameError, e
+        except:
+            msg = 'Polygon could not be converted to Numeric array'
+            raise msg
+
+        msg = 'Polygon array must be a 2d array of vertices'
+        assert len(polygon.shape) == 2, msg
+
+        msg = 'Polygon array must have two columns' 
+        assert polygon.shape[1] == 2, msg
+
+
+        msg = 'Points array must be 1 or 2 dimensional.'
+        msg += ' I got %d dimensions' %len(points.shape)
+        assert 0 < len(points.shape) < 3, msg
+
+        
+        if len(points.shape) == 1:
+            # Only one point was passed in.
+            # Convert to array of points
+            points = num.reshape(points, (1,2))
+
+    
+            msg = 'Point array must have two columns (x,y), '
+            msg += 'I got points.shape[1] == %d' %points.shape[0]
+            assert points.shape[1] == 2, msg
 
        
-    msg = 'Points array must be a 2d array. I got %s' %str(points[:30])
-    assert len(points.shape) == 2, msg
-
-    msg = 'Points array must have two columns'
-    assert points.shape[1] == 2, msg
-
-
-    N = polygon.shape[0] #Number of vertices in polygon
-    M = points.shape[0]  #Number of points
-
-
-    indices = num.zeros( M, num.Int )
+            msg = 'Points array must be a 2d array. I got %s' %str(points[:30])
+            assert len(points.shape) == 2, msg
+
+            msg = 'Points array must have two columns'
+            assert points.shape[1] == 2, msg
+
+
+    N = polygon.shape[0] # Number of vertices in polygon
+    M = points.shape[0]  # Number of points
+
+    indices = num.zeros(M, num.Int)
 
     count = _separate_points_by_polygon(points, polygon, indices,
@@ -860 +997 @@
     class Found(exceptions.Exception): pass
 
+    polygon = ensure_numeric(polygon)
+    
     point_in = False
     while not point_in:
@@ -879 +1018 @@
 
                         point = [x_delta, y_delta]
-                        #print "point",point
                         if is_inside_polygon(point, polygon, closed=False):
                             raise Found

anuga_core/source/anuga/utilities/polygon_ext.c

@@ -245 +245 @@
 
 int __separate_points_by_polygon(int M,     // Number of points
-                                int N,     // Number of polygon vertices
-                                double* points,
-                                double* polygon,
-                                long* indices,  // M-Array for storage indices
-                                int closed,
-                                int verbose) {
+                                 int N,     // Number of polygon vertices
+                                 double* points,
+                                 double* polygon,
+                                 long* indices,  // M-Array for storage indices
+                                 int closed,
+                                 int verbose) {
 
   double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j, rtol=0.0, atol=0.0;

anuga_core/source/anuga/utilities/test_polygon.py

@@ -195 +195 @@
 
 
+    def test_is_inside_polygon_quick(self):
+        """test_is_inside_polygon_quick
+        
+        Test fast version of of is_inside_polygon
+        """
+
+        # Simplest case: Polygon is the unit square
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+
+        assert is_inside_polygon_quick( (0.5, 0.5), polygon )
+        assert not is_inside_polygon_quick( (0.5, 1.5), polygon )
+        assert not is_inside_polygon_quick( (0.5, -0.5), polygon )
+        assert not is_inside_polygon_quick( (-0.5, 0.5), polygon )
+        assert not is_inside_polygon_quick( (1.5, 0.5), polygon )
+
+        # Try point on borders
+        assert is_inside_polygon_quick( (1., 0.5), polygon, closed=True)
+        assert is_inside_polygon_quick( (0.5, 1), polygon, closed=True)
+        assert is_inside_polygon_quick( (0., 0.5), polygon, closed=True)
+        assert is_inside_polygon_quick( (0.5, 0.), polygon, closed=True)
+
+        assert not is_inside_polygon_quick( (0.5, 1), polygon, closed=False)
+        assert not is_inside_polygon_quick( (0., 0.5), polygon, closed=False)
+        assert not is_inside_polygon_quick( (0.5, 0.), polygon, closed=False)
+        assert not is_inside_polygon_quick( (1., 0.5), polygon, closed=False)
+
+
+        
+        
     def test_inside_polygon_main(self):
 
@@ -223 +252 @@
         assert not is_inside_polygon( (0.5, -0.5), polygon )
 
+        
+        assert is_inside_polygon_quick( (0.5, 0.5), polygon )
+        assert is_inside_polygon_quick( (1, -0.5), polygon )
+        assert is_inside_polygon_quick( (1.5, 0), polygon ) 
+        
+        assert not is_inside_polygon_quick( (0.5, 1.5), polygon )
+        assert not is_inside_polygon_quick( (0.5, -0.5), polygon )               
 
         # Very convoluted polygon
@@ -406 +442 @@
         assert num.allclose( res, [1,2,3,5,4,0] )        
         assert count == 3
-        
+
+        
+        
+    def test_is_inside_triangle(self):
+
+
+        # Simplest case:
+        triangle = [[0, 0], [1, 0], [0.5, 1]]        
+
+        assert is_inside_triangle((0.5, 0.5), triangle)
+        assert is_inside_triangle((0.9, 0.1), triangle)        
+        assert not is_inside_triangle((0.5, 1.5), triangle)
+        assert not is_inside_triangle((0.5, -0.5), triangle)
+        assert not is_inside_triangle((-0.5, 0.5), triangle)
+        assert not is_inside_triangle((1.5, 0.5), triangle)
+
+        # Try point on borders
+        assert is_inside_triangle((0.5, 0), triangle, closed=True)
+        assert is_inside_triangle((1, 0), triangle, closed=True)
+
+        assert not is_inside_triangle((0.5, 0), triangle, closed=False)
+        assert not is_inside_triangle((1, 0), triangle, closed=False)        
+
+        # Try vertices
+        for P in triangle:
+            assert is_inside_triangle(P, triangle, closed=True)            
+            assert not is_inside_triangle(P, triangle, closed=False)                        
+        
+        
+        # Slightly different
+        triangle = [[0, 0.1], [1, -0.2], [0.5, 1]]        
+        assert is_inside_triangle((0.5, 0.5), triangle)                        
+        assert is_inside_triangle((0.4, 0.1), triangle)                                
+        assert not is_inside_triangle((1, 1), triangle)                                
+        
+        # Try vertices
+        for P in triangle:
+            assert is_inside_triangle(P, triangle, closed=True)            
+            assert not is_inside_triangle(P, triangle, closed=False)                                
+            
+            
+            
 
     def test_populate_polygon(self):
@@ -416 +493 @@
         for point in points:
             assert is_inside_polygon(point, polygon)
-
-
-        #Very convoluted polygon
+            assert is_inside_polygon_quick(point, polygon)            
+
+
+        # Very convoluted polygon
         polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
 
@@ -426 +504 @@
         for point in points:
             assert is_inside_polygon(point, polygon)
+            assert is_inside_polygon_quick(point, polygon)                        
 
 
@@ -1782 +1861 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Polygon,'test')
+    suite = unittest.makeSuite(Test_Polygon, 'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)