- Timestamp:
- Nov 12, 2008, 12:14:33 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source_numpy_conversion/anuga/utilities/polygon.py
r5902 r5948 10 10 # #print 'Could not find scipy - using Numeric' 11 11 12 ##from numpy import float, int, zeros, ones, array, concatenate, reshape, dot, allclose, newaxis, ascontiguousarray13 12 import numpy 14 15 13 16 14 from math import sqrt 17 15 from anuga.utilities.numerical_tools import ensure_numeric 18 16 from anuga.geospatial_data.geospatial_data import ensure_absolute 19 20 21 def point_on_line(point, line, rtol=0.0, atol=0.0): 17 from anuga.config import Float, Int 18 19 20 def point_on_line(point, line, rtol=1.0e-5, atol=1.0e-8): 22 21 """Determine whether a point is on a line segment 23 22 … … 32 31 """ 33 32 34 # FIXME(Ole): Perhaps make defaults as in allclose: rtol=1.0e-5, atol=1.0e-835 36 33 point = ensure_numeric(point) 37 34 line = ensure_numeric(line) … … 45 42 46 43 47 48 49 50 def intersection(line0, line1): 44 ###### 45 # Result functions used in intersection() below for collinear lines. 46 # (p0,p1) defines line 0, (p2,p3) defines line 1. 47 ###### 48 49 # result functions for possible states 50 def lines_dont_coincide(p0,p1,p2,p3): return (3, None) 51 def lines_0_fully_included_in_1(p0,p1,p2,p3): return (2, numpy.array([p0,p1])) 52 def lines_1_fully_included_in_0(p0,p1,p2,p3): return (2, numpy.array([p2,p3])) 53 def lines_overlap_same_direction(p0,p1,p2,p3): return (2, numpy.array([p0,p3])) 54 def lines_overlap_same_direction2(p0,p1,p2,p3): return (2, numpy.array([p2,p1])) 55 def lines_overlap_opposite_direction(p0,p1,p2,p3): return (2, numpy.array([p0,p2])) 56 def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2, numpy.array([p3,p1])) 57 58 # this function called when an impossible state is found 59 def lines_error(p1, p2, p3, p4): raise RuntimeError, "INTERNAL ERROR" 60 61 # 0s1 0e1 1s0 1e0 # line 0 starts on 1, 0 ends 1, 1 starts 0, 1 ends 0 62 collinear_result = { (False, False, False, False): lines_dont_coincide, 63 (False, False, False, True ): lines_error, 64 (False, False, True, False): lines_error, 65 (False, False, True, True ): lines_1_fully_included_in_0, 66 (False, True, False, False): lines_error, 67 (False, True, False, True ): lines_overlap_opposite_direction2, 68 (False, True, True, False): lines_overlap_same_direction2, 69 (False, True, True, True ): lines_1_fully_included_in_0, 70 (True, False, False, False): lines_error, 71 (True, False, False, True ): lines_overlap_same_direction, 72 (True, False, True, False): lines_overlap_opposite_direction, 73 (True, False, True, True ): lines_1_fully_included_in_0, 74 (True, True, False, False): lines_0_fully_included_in_1, 75 (True, True, False, True ): lines_0_fully_included_in_1, 76 (True, True, True, False): lines_0_fully_included_in_1, 77 (True, True, True, True ): lines_0_fully_included_in_1 78 } 79 80 def intersection(line0, line1, rtol=1.0e-5, atol=1.0e-8): 81 """Returns intersecting point between two line segments or None 82 (if parallel or no intersection is found). 83 84 However, if parallel lines coincide partly (i.e. shara a common segment, 85 the line segment where lines coincide is returned 86 87 88 Inputs: 89 line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]] 90 A line can also be a 2x2 numpy array with each row 91 corresponding to a point. 92 93 94 Output: 95 status, value 96 97 where status and value is interpreted as follows 98 99 status == 0: no intersection, value set to None. 100 status == 1: intersection point found and returned in value as [x,y]. 101 status == 2: Collinear overlapping lines found. Value takes the form [[x0,y0], [x1,y1]]. 102 status == 3: Collinear non-overlapping lines. Value set to None. 103 status == 4: Lines are parallel with a fixed distance apart. Value set to None. 104 105 """ 106 107 # FIXME (Ole): Write this in C 108 109 line0 = ensure_numeric(line0, Float) 110 line1 = ensure_numeric(line1, Float) 111 112 x0 = line0[0,0]; y0 = line0[0,1] 113 x1 = line0[1,0]; y1 = line0[1,1] 114 115 x2 = line1[0,0]; y2 = line1[0,1] 116 x3 = line1[1,0]; y3 = line1[1,1] 117 118 denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0) 119 u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2) 120 u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0) 121 122 if numpy.allclose(denom, 0.0, rtol=rtol, atol=atol): 123 # Lines are parallel - check if they are collinear 124 if numpy.allclose([u0, u1], 0.0, rtol=rtol, atol=atol): 125 # We now know that the lines are collinear 126 state_tuple = (point_on_line([x0, y0], line1, rtol=rtol, atol=atol), 127 point_on_line([x1, y1], line1, rtol=rtol, atol=atol), 128 point_on_line([x2, y2], line0, rtol=rtol, atol=atol), 129 point_on_line([x3, y3], line0, rtol=rtol, atol=atol)) 130 131 return collinear_result[state_tuple]([x0,y0],[x1,y1],[x2,y2],[x3,y3]) 132 else: 133 # Lines are parallel but aren't collinear 134 return 4, None #FIXME (Ole): Add distance here instead of None 135 else: 136 # Lines are not parallel, check if they intersect 137 u0 = u0/denom 138 u1 = u1/denom 139 140 x = x0 + u0*(x1-x0) 141 y = y0 + u0*(y1-y0) 142 143 # Sanity check - can be removed to speed up if needed 144 assert numpy.allclose(x, x2 + u1*(x3-x2), rtol=rtol, atol=atol) 145 assert numpy.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol) 146 147 # Check if point found lies within given line segments 148 if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 149 # We have intersection 150 return 1, numpy.array([x, y]) 151 else: 152 # No intersection 153 return 0, None 154 155 156 def NEW_C_intersection(line0, line1): 157 #FIXME(Ole): To write in C 51 158 """Returns intersecting point between two line segments or None 52 159 (if parallel or no intersection is found). … … 75 182 """ 76 183 77 # FIXME (Ole): Write this in C 78 79 line0 = ensure_numeric(line0, numpy.float) 80 line1 = ensure_numeric(line1, numpy.float) 81 82 x0 = line0[0,0]; y0 = line0[0,1] 83 x1 = line0[1,0]; y1 = line0[1,1] 84 85 x2 = line1[0,0]; y2 = line1[0,1] 86 x3 = line1[1,0]; y3 = line1[1,1] 87 88 denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0) 89 u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2) 90 u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0) 91 92 if numpy.allclose(denom, 0.0): 93 # Lines are parallel - check if they coincide on a shared a segment 94 95 if numpy.allclose( [u0, u1], 0.0 ): 96 # We now know that the lines if continued coincide 97 # The remaining check will establish if the finite lines share a segment 98 99 line0_starts_on_line1 = line0_ends_on_line1 =\ 100 line1_starts_on_line0 = line1_ends_on_line0 = False 101 102 if point_on_line([x0, y0], line1): 103 line0_starts_on_line1 = True 104 105 if point_on_line([x1, y1], line1): 106 line0_ends_on_line1 = True 107 108 if point_on_line([x2, y2], line0): 109 line1_starts_on_line0 = True 110 111 if point_on_line([x3, y3], line0): 112 line1_ends_on_line0 = True 113 114 if not(line0_starts_on_line1 or line0_ends_on_line1\ 115 or line1_starts_on_line0 or line1_ends_on_line0): 116 # Lines are parallel and would coincide if extended, but not as they are. 117 return 3, None 118 119 120 # One line fully included in the other. Use direction of included line 121 if line0_starts_on_line1 and line0_ends_on_line1: 122 # Shared segment is line0 fully included in line1 123 segment = numpy.array([[x0, y0], [x1, y1]]) 124 125 if line1_starts_on_line0 and line1_ends_on_line0: 126 # Shared segment is line1 fully included in line0 127 segment = numpy.array([[x2, y2], [x3, y3]]) 128 129 130 # Overlap with lines are oriented the same way 131 if line0_starts_on_line1 and line1_ends_on_line0: 132 # Shared segment from line0 start to line 1 end 133 segment = numpy.array([[x0, y0], [x3, y3]]) 134 135 if line1_starts_on_line0 and line0_ends_on_line1: 136 # Shared segment from line1 start to line 0 end 137 segment = numpy.array([[x2, y2], [x1, y1]]) 138 139 140 # Overlap in opposite directions - use direction of line0 141 if line0_starts_on_line1 and line1_starts_on_line0: 142 # Shared segment from line0 start to line 1 end 143 segment = numpy.array([[x0, y0], [x2, y2]]) 144 145 if line0_ends_on_line1 and line1_ends_on_line0: 146 # Shared segment from line0 start to line 1 end 147 segment = numpy.array([[x3, y3], [x1, y1]]) 148 149 150 return 2, segment 151 else: 152 # Lines are parallel but they don't coincide 153 return 4, None #FIXME (Ole): Add distance here instead of None 154 155 else: 156 # Lines are not parallel or coinciding 157 u0 = u0/denom 158 u1 = u1/denom 159 160 x = x0 + u0*(x1-x0) 161 y = y0 + u0*(y1-y0) 162 163 # Sanity check - can be removed to speed up if needed 164 assert numpy.allclose(x, x2 + u1*(x3-x2)) 165 assert numpy.allclose(y, y2 + u1*(y3-y2)) 166 167 # Check if point found lies within given line segments 168 if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 169 # We have intersection 170 171 return 1, numpy.array([x, y]) 172 else: 173 # No intersection 174 return 0, None 175 176 177 def NEW_C_intersection(line0, line1): 178 #FIXME(Ole): To write in C 179 """Returns intersecting point between two line segments or None 180 (if parallel or no intersection is found). 181 182 However, if parallel lines coincide partly (i.e. shara a common segment, 183 the line segment where lines coincide is returned 184 185 186 Inputs: 187 line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]] 188 A line can also be a 2x2 numeric array with each row 189 corresponding to a point. 190 191 192 Output: 193 status, value 194 195 where status is interpreted as follows 196 197 status == 0: no intersection with value set to None 198 status == 1: One intersection point found and returned in value as [x,y] 199 status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]] 200 status == 3: Lines would coincide but only if extended. Value set to None 201 status == 4: Lines are parallel with a fixed distance apart. Value set to None. 202 203 """ 204 205 206 line0 = ensure_numeric(line0, numpy.float) 207 line1 = ensure_numeric(line1, numpy.float) 184 185 line0 = ensure_numeric(line0, Float) 186 line1 = ensure_numeric(line1, Float) 208 187 209 188 status, value = _intersection(line0[0,0], line0[0,1], … … 313 292 #if verbose: print 'Checking input to outside_polygon' 314 293 try: 315 points = ensure_numeric(points, numpy.float)294 points = ensure_numeric(points, Float) 316 295 except NameError, e: 317 296 raise NameError, e … … 321 300 322 301 try: 323 polygon = ensure_numeric(polygon, numpy.float)302 polygon = ensure_numeric(polygon, Float) 324 303 except NameError, e: 325 304 raise NameError, e … … 355 334 #if verbose: print 'Checking input to outside_polygon' 356 335 try: 357 points = ensure_numeric(points, numpy.float)336 points = ensure_numeric(points, Float) 358 337 except NameError, e: 359 338 raise NameError, e … … 363 342 364 343 try: 365 polygon = ensure_numeric(polygon, numpy.float)344 polygon = ensure_numeric(polygon, Float) 366 345 except NameError, e: 367 346 raise NameError, e … … 442 421 ## print 'Before: points=%s, flags=%s' % (type(points), str(points.flags)) 443 422 try: 444 ## points = numpy.ascontiguousarray(ensure_numeric(points, numpy.float))445 points = ensure_numeric(points, numpy.float)423 ## points = numpy.ascontiguousarray(ensure_numeric(points, Float)) 424 points = ensure_numeric(points, Float) 446 425 except NameError, e: 447 426 raise NameError, e … … 453 432 #if verbose: print 'Checking input to separate_points_by_polygon 2' 454 433 try: 455 ## polygon = numpy.ascontiguousarray(ensure_numeric(polygon, numpy.float))456 polygon = ensure_numeric(polygon, numpy.float)434 ## polygon = numpy.ascontiguousarray(ensure_numeric(polygon, Float)) 435 polygon = ensure_numeric(polygon, Float) 457 436 except NameError, e: 458 437 raise NameError, e … … 495 474 496 475 497 indices = numpy.zeros( M, numpy.int )476 indices = numpy.zeros( M, Int ) 498 477 499 478 count = _separate_points_by_polygon(points, polygon, indices, … … 622 601 623 602 try: 624 polygon = ensure_numeric(polygon, numpy.float)603 polygon = ensure_numeric(polygon, Float) 625 604 except NameError, e: 626 605 raise NameError, e … … 733 712 734 713 def __call__(self, x, y): 735 x = numpy.array(x).astype( numpy.float)736 y = numpy.array(y).astype( numpy.float)714 x = numpy.array(x).astype(Float) 715 y = numpy.array(y).astype(Float) 737 716 738 717 assert len(x.shape) == 1 and len(y.shape) == 1 … … 746 725 z = self.default(x, y) 747 726 else: 748 z = numpy.ones(N, numpy.float) * self.default727 z = numpy.ones(N, Float) * self.default 749 728 750 729 for polygon, value in self.regions:
Note: See TracChangeset
for help on using the changeset viewer.
