Changeset 7858 for trunk/anuga_core/source/anuga/geometry/polygon.py
 Timestamp:
 Jun 18, 2010, 4:43:10 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/geometry/polygon.py
r7841 r7858 10 10 import anuga.utilities.log as log 11 11 12 from aabb import AABB 12 13 13 14 ## … … 36 37 37 38 res = _point_on_line(point[0], point[1], 38 line[0, 0], line[0,1],39 line[1, 0], line[1,1],39 line[0, 0], line[0, 1], 40 line[1, 0], line[1, 1], 40 41 rtol, atol) 41 42 … … 49 50 50 51 # result functions for possible states 51 def lines_dont_coincide(p0,p1,p2,p3): return (3, None) 52 def lines_0_fully_included_in_1(p0,p1,p2,p3): return (2, 53 num.array([p0,p1])) 54 def lines_1_fully_included_in_0(p0,p1,p2,p3): return (2, 55 num.array([p2,p3])) 56 def lines_overlap_same_direction(p0,p1,p2,p3): return (2, 57 num.array([p0,p3])) 58 def lines_overlap_same_direction2(p0,p1,p2,p3): return (2, 59 num.array([p2,p1])) 60 def lines_overlap_opposite_direction(p0,p1,p2,p3): return (2, 61 num.array([p0,p2])) 62 def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2, 63 num.array([p3,p1])) 52 def lines_dont_coincide(p0, p1, p2, p3): 53 return (3, None) 54 55 def lines_0_fully_included_in_1(p0, p1, p2, p3): 56 return (2, num.array([p0, p1])) 57 58 def lines_1_fully_included_in_0(p0, p1, p2, p3): 59 return (2, num.array([p2, p3])) 60 61 def lines_overlap_same_direction(p0, p1, p2, p3): 62 return (2, num.array([p0, p3])) 63 64 def lines_overlap_same_direction2(p0, p1, p2, p3): 65 return (2, num.array([p2, p1])) 66 67 def lines_overlap_opposite_direction(p0, p1, p2, p3): 68 return (2, num.array([p0, p2])) 69 70 def lines_overlap_opposite_direction2(p0, p1, p2, p3): 71 return (2, num.array([p3, p1])) 64 72 65 73 # this function called when an impossible state is found … … 146 154 point_on_line([x3, y3], line0, rtol=rtol, atol=atol)) 147 155 148 return collinear_result[state_tuple]([x0, y0], [x1,y1],149 [x2, y2], [x3,y3])156 return collinear_result[state_tuple]([x0, y0], [x1, y1], 157 [x2, y2], [x3, y3]) 150 158 else: 151 159 # Lines are parallel but aren't collinear … … 269 277 return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol)) 270 278 271 272 273 # FIXME (Ole): The rest of this function has been made274 # obsolete by the C extension.275 276 # Quickly reject points that are clearly outside277 if point[0] < min(triangle[:, 0]):278 return False279 if point[0] > max(triangle[:, 0]):280 return False281 if point[1] < min(triangle[:, 1]):282 return False283 if point[1] > max(triangle[:, 1]):284 return False285 286 287 # Start search288 A = triangle[0, :]289 B = triangle[1, :]290 C = triangle[2, :]291 292 # Now check if point lies wholly inside triangle293 v0 = CA294 v1 = BA295 296 a00 = num.inner(v0, v0)297 a10 = a01 = num.inner(v0, v1)298 a11 = num.inner(v1, v1)299 300 denom = a11*a00  a01*a10301 302 if abs(denom) > 0.0:303 v = pointA304 b0 = num.inner(v0, v)305 b1 = num.inner(v1, v)306 307 alpha = (b0*a11  b1*a01)/denom308 beta = (b1*a00  b0*a10)/denom309 310 if (alpha > 0.0) and (beta > 0.0) and (alpha+beta < 1.0):311 return True312 313 314 if closed is True:315 # Check if point lies on one of the edges316 317 for X, Y in [[A, B], [B, C], [C, A]]:318 res = _point_on_line(point[0], point[1],319 X[0], X[1],320 Y[0], Y[1],321 rtol, atol)322 323 if res:324 return True325 326 return False327 328 329 279 def is_complex(polygon, verbose=False): 330 280 """Check if a polygon is complex (selfintersecting). … … 339 289 340 290 def key_xpos(item): 291 """ Return the x coord out of the passed point for sorting key. """ 341 292 return (item[0][0]) 342 293 … … 377 328 print 'Selfintersecting polygon found, type ', type 378 329 print 'point', point, 379 print 'vertices: ', leftmost, '  ', l_x[cmp] 330 print 'vertices: ', leftmost, '  ', l_x[cmp] 380 331 return True 381 332 cmp += 1 … … 422 373 try: 423 374 points = ensure_absolute(points) 424 except NameError, e :425 raise NameError, e 375 except NameError, err: 376 raise NameError, err 426 377 except: 427 378 # If this fails it is going to be because the points can't be … … 513 464 if len(points.shape) == 1: 514 465 # Only one point was passed in. Convert to array of points 515 points = num.reshape(points, (1, 2))466 points = num.reshape(points, (1, 2)) 516 467 517 468 indices, count = separate_points_by_polygon(points, polygon, … … 559 510 if len(points.shape) == 1: 560 511 # Only one point was passed in. Convert to array of points 561 points = num.reshape(points, (1, 2))512 points = num.reshape(points, (1, 2)) 562 513 563 514 indices, count = separate_points_by_polygon(points, polygon, … … 632 583 except: 633 584 msg = 'Points could not be converted to numeric array' 634 raise msg585 raise Exception(msg) 635 586 636 587 try: 637 588 polygon = ensure_numeric(polygon, num.float) 638 589 except NameError, e: 639 raise NameError , e590 raise NameError(e) 640 591 except: 641 592 msg = 'Polygon could not be converted to numeric array' 642 raise msg593 raise Exception(msg) 643 594 644 595 msg = 'Polygon array must be a 2d array of vertices' … … 646 597 647 598 msg = 'Polygon array must have two columns' 648 assert polygon.shape[1] ==2, msg599 assert polygon.shape[1] == 2, msg 649 600 650 601 msg = ('Points array must be 1 or 2 dimensional. ' … … 654 605 if len(points.shape) == 1: 655 606 # Only one point was passed in. Convert to array of points. 656 points = num.reshape(points, (1, 2))607 points = num.reshape(points, (1, 2)) 657 608 658 609 msg = ('Point array must have two columns (x,y), ' … … 663 614 msg = ('Points array must be a 2d array. I got %s.' 664 615 % str(points[:30])) 665 assert len(points.shape) ==2, msg616 assert len(points.shape) == 2, msg 666 617 667 618 msg = 'Points array must have two columns' 668 assert points.shape[1] ==2, msg619 assert points.shape[1] == 2, msg 669 620 670 621 N = polygon.shape[0] # Number of vertices in polygon … … 681 632 return indices, count 682 633 683 ## 684 # @brief Determine area of a polygon. 685 # @param input_polygon The polygon to get area of. 686 # @return A scalar value for the polygon area. 634 687 635 def polygon_area(input_polygon): 688 636 """ Determine area of arbitrary polygon. 689 637 690 Reference: http://mathworld.wolfram.com/PolygonArea.html 691 """ 692 638 input_polygon The polygon to get area of. 639 640 return A scalar value for the polygon area. 641 642 Reference: http://mathworld.wolfram.com/PolygonArea.html 643 """ 693 644 # Move polygon to origin (0,0) to avoid rounding errors 694 645 # This makes a copy of the polygon to avoid destroying it 695 646 input_polygon = ensure_numeric(input_polygon) 696 min_x = min(input_polygon[:, 0])697 min_y = min(input_polygon[:, 1])647 min_x = min(input_polygon[:, 0]) 648 min_y = min(input_polygon[:, 1]) 698 649 polygon = input_polygon  [min_x, min_y] 699 650 … … 742 693 Outputs: 743 694 744  list of min and max of x and y coordinates745 695  plot of polygons 746 696 """ … … 754 704 ion() 755 705 hold(True) 756 757 minx = 1e10758 maxx = 0.0759 miny = 1e10760 maxy = 0.0761 706 762 707 if label is None: … … 772 717 alpha = max(0.0, min(1.0, alpha)) 773 718 774 n = len(polygons_points)719 num_points = len(polygons_points) 775 720 colour = [] 776 721 if style is None: 777 722 style_type = 'line' 778 723 style = [] 779 for i in range(n ):724 for i in range(num_points): 780 725 style.append(style_type) 781 726 colour.append('b') 782 727 else: 783 for s in style: 784 if s == 'line': colour.append('b') 785 if s == 'outside': colour.append('r.') 786 if s == 'point': colour.append('g.') 787 if s != 'line': 788 if s != 'outside': 789 if s != 'point': 790 colour.append(s) 728 for style_name in style: 729 if style_name == 'line': 730 colour.append('b') 731 if style_name == 'outside': 732 colour.append('r.') 733 if style_name == 'point': 734 colour.append('g.') 735 if style_name not in ['line', 'outside', 'point']: 736 colour.append(style_name) 791 737 792 738 for i, item in enumerate(polygons_points): 793 x, y = poly_xy(item) 794 if min(x) < minx: 795 minx = min(x) 796 if max(x) > maxx: 797 maxx = max(x) 798 if min(y) < miny: 799 miny = min(y) 800 if max(y) > maxy: 801 maxy = max(y) 802 plot(x, y, colour[i]) 739 pt_x, pt_y = _poly_xy(item) 740 plot(pt_x, pt_y, colour[i]) 803 741 if alpha: 804 fill( x,y, colour[i], alpha=alpha)742 fill(pt_x, pt_y, colour[i], alpha=alpha) 805 743 xlabel('x') 806 744 ylabel('y') … … 814 752 close('all') 815 753 816 vec = [minx, maxx, miny, maxy] 817 return vec 818 819 820 def poly_xy(polygon): 754 755 def _poly_xy(polygon): 821 756 """ this is used within plot_polygons so need to duplicate 822 757 the first point so can have closed polygon in plot … … 829 764 try: 830 765 polygon = ensure_numeric(polygon, num.float) 831 except NameError, e :832 raise NameError, e 766 except NameError, err: 767 raise NameError, err 833 768 except: 834 769 msg = ('Polygon %s could not be converted to numeric array' … … 836 771 raise Exception, msg 837 772 838 x = polygon[:, 0] 839 y = polygon[:, 1] 840 x = num.concatenate((x, [polygon[0, 0]]), axis = 0) 841 y = num.concatenate((y, [polygon[0, 1]]), axis = 0) 842 843 return x, y 844 845 846 class Polygon_function: 847 """Create callable object f: x,y > z, where a,y,z are vectors and 848 where f will return different values depending on whether x,y belongs 849 to specified polygons. 850 851 To instantiate: 852 853 Polygon_function(polygons) 854 855 where polygons is a list of tuples of the form 856 857 [ (P0, v0), (P1, v1), ...] 858 859 with Pi being lists of vertices defining polygons and vi either 860 constants or functions of x,y to be applied to points with the polygon. 861 862 The function takes an optional argument, default which is the value 863 (or function) to used for points not belonging to any polygon. 864 For example: 865 866 Polygon_function(polygons, default = 0.03) 867 868 If omitted the default value will be 0.0 869 870 Note: If two polygons overlap, the one last in the list takes precedence 871 872 Coordinates specified in the call are assumed to be relative to the 873 origin (georeference) e.g. used by domain. 874 By specifying the optional argument georeference, 875 all points are made relative. 876 877 FIXME: This should really work with geo_spatial point sets. 878 """ 879 880 def __init__(self, regions, default=0.0, geo_reference=None): 881 """ 882 # @brief Create instance of a polygon function. 883 # @param regions A list of (x,y) tuples defining a polygon. 884 # @param default Value or function returning value for points outside poly. 885 # @param geo_reference ?? 886 """ 887 888 try: 889 len(regions) 890 except: 891 msg = ('Polygon_function takes a list of pairs (polygon, value).' 892 'Got %s' % str(regions)) 893 raise Exception, msg 894 895 T = regions[0] 896 897 if isinstance(T, basestring): 898 msg = ('You passed in a list of text values into polygon_function ' 899 'instead of a list of pairs (polygon, value): "%s"' 900 % str(T)) 901 raise Exception, msg 902 903 try: 904 a = len(T) 905 except: 906 msg = ('Polygon_function takes a list of pairs (polygon, value). ' 907 'Got %s' % str(T)) 908 raise Exception, msg 909 910 msg = ('Each entry in regions have two components: (polygon, value). ' 911 'I got %s' % str(T)) 912 assert a == 2, msg 913 914 if geo_reference is None: 915 from anuga.coordinate_transforms.geo_reference import Geo_reference 916 geo_reference = Geo_reference() 917 918 self.default = default 919 920 # Make points in polygons relative to geo_reference 921 self.regions = [] 922 for polygon, value in regions: 923 P = geo_reference.change_points_geo_ref(polygon) 924 self.regions.append((P, value)) 925 926 def __call__(self, x, y): 927 """ 928 # @brief Implement the 'callable' property of Polygon_function. 929 # @param x List of x coordinates of points ot interest. 930 # @param y List of y coordinates of points ot interest. 931 """ 932 x = num.array(x, num.float) 933 y = num.array(y, num.float) 934 935 # x and y must be onedimensional and same length 936 assert len(x.shape) == 1 and len(y.shape) == 1 937 N = x.shape[0] 938 assert y.shape[0] == N 939 940 points = num.ascontiguousarray(num.concatenate((x[:, num.newaxis], 941 y[:, num.newaxis]), 942 axis = 1 )) 943 944 if callable(self.default): 945 z = self.default(x, y) 946 else: 947 z = num.ones(N, num.float) * self.default 948 949 for polygon, value in self.regions: 950 indices = inside_polygon(points, polygon) 951 952 # FIXME: This needs to be vectorised 953 if callable(value): 954 for i in indices: 955 xx = num.array([x[i]]) 956 yy = num.array([y[i]]) 957 z[i] = value(xx, yy)[0] 958 else: 959 for i in indices: 960 z[i] = value 961 962 if len(z) == 0: 963 msg = ('Warning: points provided to Polygon function did not fall ' 964 'within its regions in [%.2f, %.2f], y in [%.2f, %.2f]' 965 % (min(x), max(x), min(y), max(y))) 966 log.critical(msg) 967 968 return z 773 pts_x = num.concatenate((polygon[:, 0], [polygon[0, 0]]), axis = 0) 774 pts_y = num.concatenate((polygon[:, 1], [polygon[0, 1]]), axis = 0) 775 776 return pts_x, pts_y 777 969 778 970 779 ################################################################################ … … 973 782 974 783 def read_polygon(filename, delimiter=','): 975 """Read points assumed to form a polygon. 976 977 # @param filename Path to file containing polygon data. 978 # @param delimiter Delimiter to split polygon data with. 979 # @return A list of point data from the polygon file. 980 981 982 There must be exactly two numbers in each line separated by the delimiter. 983 No header. 784 """ Read points assumed to form a polygon. 785 786 Also checks to make sure polygon is not complex (selfintersecting). 787 788 filename Path to file containing polygon data. 789 delimiter Delimiter to split polygon data with. 790 A list of point data from the polygon file. 791 792 There must be exactly two numbers in each line separated by the delimiter. 793 No header. 984 794 """ 985 795 … … 1002 812 return polygon 1003 813 1004 ## 1005 # @brief Write polygon data to a file. 1006 # @param polygon Polygon points to write to file. 1007 # @param filename Path to file to write. 1008 # @note Delimiter is assumed to be a comma. 814 1009 815 def write_polygon(polygon, filename=None): 1010 816 """Write polygon to csv file. … … 1047 853 1048 854 # Find outer extent of polygon 1049 max_x = min_x = polygon[0][0] 1050 max_y = min_y = polygon[0][1] 1051 for point in polygon[1:]: 1052 x = point[0] 1053 if x > max_x: max_x = x 1054 if x < min_x: min_x = x 1055 y = point[1] 1056 if y > max_y: max_y = y 1057 if y < min_y: min_y = y 1058 855 extents = AABB(polygon) 856 1059 857 while len(points) < number_of_points: 1060 x = uniform(min_x, max_x)1061 y = uniform(min_y, max_y)858 rand_x = uniform(extents.xmin, extents.xmax) 859 rand_y = uniform(extents.ymin, extents.ymax) 1062 860 1063 861 append = False 1064 if is_inside_polygon([ x,y], polygon):862 if is_inside_polygon([rand_x, rand_y], polygon): 1065 863 append = True 1066 864 … … 1068 866 if exclude is not None: 1069 867 for ex_poly in exclude: 1070 if is_inside_polygon([ x,y], ex_poly):868 if is_inside_polygon([rand_x, rand_y], ex_poly): 1071 869 append = False 1072 870 1073 871 if append is True: 1074 points.append([ x,y])872 points.append([rand_x, rand_y]) 1075 873 1076 874 return points … … 1091 889 """ 1092 890 1093 import exceptions1094 1095 class Found(exceptions.Exception):1096 pass1097 1098 891 polygon = ensure_numeric(polygon) 1099 892 1100 point_in = False 1101 while not point_in: 1102 try: 1103 for poly_point in polygon: # [1:]: 1104 for x_mult in range(1, 2): 1105 for y_mult in range(1, 2): 1106 x = poly_point[0] 1107 y = poly_point[1] 1108 1109 if x == 0: 1110 x_delta = x_mult * delta 1111 else: 1112 x_delta = x + x_mult*x*delta 1113 1114 if y == 0: 1115 y_delta = y_mult * delta 1116 else: 1117 y_delta = y + y_mult*y*delta 1118 1119 point = [x_delta, y_delta] 1120 1121 if is_inside_polygon(point, polygon, closed=False): 1122 raise Found 1123 except Found: 1124 point_in = True 1125 else: 1126 delta = delta * 0.1 1127 1128 return point 893 while True: 894 for poly_point in polygon: 895 for x_mult in range(1, 2): 896 for y_mult in range(1, 2): 897 pt_x, pt_y = poly_point 898 899 if pt_x == 0: 900 x_delta = x_mult * delta 901 else: 902 x_delta = pt_x + x_mult*pt_x*delta 903 904 if pt_y == 0: 905 y_delta = y_mult * delta 906 else: 907 y_delta = pt_y + y_mult*pt_y*delta 908 909 point = [x_delta, y_delta] 910 911 if is_inside_polygon(point, polygon, closed=False): 912 return point 913 delta = delta * 0.1 1129 914 1130 915 … … 1212 997 reduced_polygon = [] 1213 998 for i, point in enumerate(polygon): 1214 x = point[0] 1215 y = point[1] 1216 if x in [min_x, max_x] and y in [min_y, max_y]: 999 if point[0] in [min_x, max_x] and point[1] in [min_y, max_y]: 1217 1000 # Keep 1218 1001 reduced_polygon.append(point) … … 1318 1101 1319 1102 else: 1320 error_msg = 'C implementations could not be accessed by %s.\n ' %__file__1321 error_msg+= 'Make sure compile_all.py has been run as described in '1322 error_msg+= 'the ANUGA installation guide.'1323 raise Exception( error_msg)1103 ERROR_MSG = 'C implementations could not be accessed by %s.\n ' % __file__ 1104 ERROR_MSG += 'Make sure compile_all.py has been run as described in ' 1105 ERROR_MSG += 'the ANUGA installation guide.' 1106 raise Exception(ERROR_MSG) 1324 1107 1325 1108
Note: See TracChangeset
for help on using the changeset viewer.