# Changeset 6534

Ignore:
Timestamp:
Mar 17, 2009, 5:44:57 PM (15 years ago)
Message:

Added new function is_inside_triangle. This helps the fitting functions.

Location:
anuga_core/source/anuga/utilities
Files:
4 edited

Unmodified
Removed
• ## anuga_core/source/anuga/utilities/numerical_tools.py

 r6438 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

 r6223 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 raise msg polygon = ensure_absolute(polygon) try: polygon = ensure_absolute(polygon) 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 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 regarded as belonging to the polygon (closed = True) or not (closed = False) check_input: Allows faster execution if set to False Outputs: """ #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, class Found(exceptions.Exception): pass polygon = ensure_numeric(polygon) point_in = False while not point_in: 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

 r6189 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

 r6215 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): 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 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): 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]] for point in points: assert is_inside_polygon(point, polygon) assert is_inside_polygon_quick(point, polygon) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(Test_Polygon,'test') suite = unittest.makeSuite(Test_Polygon, 'test') runner = unittest.TextTestRunner() runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.