Changeset 5225
- Timestamp:
- Apr 21, 2008, 5:20:32 PM (17 years ago)
- Location:
- anuga_core/source/anuga/utilities
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/utilities/polygon.py
r5114 r5225 10 10 # #print 'Could not find scipy - using Numeric' 11 11 12 from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot 12 from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot, allclose 13 13 14 14 … … 17 17 from anuga.geospatial_data.geospatial_data import ensure_absolute 18 18 19 def point_on_line( x, y, x0, y0, x1, y1):19 def point_on_line(point, line): 20 20 """Determine whether a point is on a line segment 21 21 22 Input: x, y, x0, x0, x1, y1: where 23 point is given by x, y 24 line is given by (x0, y0) and (x1, y1) 25 26 """ 27 28 a = array([x - x0, y - y0]) 29 a_normal = array([a[1], -a[0]]) 30 31 b = array([x1 - x0, y1 - y0]) 32 33 if dot(a_normal, b) == 0: 34 #Point is somewhere on the infinite extension of the line 35 36 len_a = sqrt(sum(a**2)) 37 len_b = sqrt(sum(b**2)) 38 if dot(a, b) >= 0 and len_a <= len_b: 39 return True 40 else: 41 return False 22 Input: 23 point is given by [x, y] 24 line is given by [x0, y0], [x1, y1]] or 25 the equivalent 2x2 Numeric array with each row corresponding to a point. 26 27 Output: 28 29 """ 30 31 point = ensure_numeric(point) 32 line = ensure_numeric(line) 33 34 35 res = _point_on_line(point[0], point[1], 36 line[0,0], line[0,1], 37 line[1,0], line[1,1]) 38 39 return bool(res) 40 41 42 43 44 45 def intersection(line0, line1): 46 """Returns intersecting point between two line segments or None (if parallel or no intersection is found) 47 48 Inputs: 49 line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]] 50 A line can also be a 2x2 numeric array with each row corresponding to a point. 51 52 53 Output: 54 Point [x,y] or None. 55 56 If line extensions intersect outside their limits, None is returned as well. 57 58 """ 59 60 # FIXME (Ole): Write this in C 61 62 line0 = ensure_numeric(line0, Float) 63 line1 = ensure_numeric(line1, Float) 64 65 x0 = line0[0,0]; y0 = line0[0,1] 66 x1 = line0[1,0]; y1 = line0[1,1] 67 68 x2 = line1[0,0]; y2 = line1[0,1] 69 x3 = line1[1,0]; y3 = line1[1,1] 70 71 denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0) 72 73 #print 'denom', denom, line0, line1 74 if allclose(denom, 0.0): 75 # Lines are parallel 76 return None 42 77 else: 43 return False 78 u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2) 79 u0 = u0/denom 80 81 u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0) 82 u1 = u1/denom 83 84 x = x0 + u0*(x1-x0) 85 y = y0 + u0*(y1-y0) 86 87 assert allclose(x, x2 + u1*(x3-x2)) 88 assert allclose(y, y2 + u1*(y3-y2)) 89 90 # Check if point found lies within given line segments 91 if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 92 # We have intersection 93 94 # Need tolerances if going ahead with this check 95 #msg = 'Algorithm error. Intersection was detected but point (%f, %f) does not lie ' %(x,y) 96 #msg += 'on line0: %s' %(line0) 97 #assert point_on_line([x, y], line0), msg 98 #msg = 'Algorithm error. Intersection was detected but point (%f, %f) does not lie ' %(x,y) 99 #msg += 'on line1: %s' %(line1) 100 #assert point_on_line([x, y], line1), msg 101 102 return [x, y] 103 else: 104 return None 44 105 45 106 … … 305 366 indices = zeros( M, Int ) 306 367 307 from polygon_ext import separate_points_by_polygon as sep_points 308 count = sep_points(points, polygon, indices, 309 int(closed), int(verbose)) 368 count = _separate_points_by_polygon(points, polygon, indices, 369 int(closed), int(verbose)) 310 370 311 371 if verbose: print 'Found %d points (out of %d) inside polygon'\ … … 780 840 from anuga.utilities.compile import can_use_C_extension 781 841 if can_use_C_extension('polygon_ext.c'): 782 # Underlying C implementations can be accessed 783 784 from polygon_ext import point_on_line 842 # Underlying C implementations can be accessed 843 from polygon_ext import _point_on_line 844 from polygon_ext import _separate_points_by_polygon 845 785 846 else: 786 847 msg = 'C implementations could not be accessed by %s.\n ' %__file__ -
anuga_core/source/anuga/utilities/polygon_ext.c
r4804 r5225 19 19 20 20 21 int _ point_on_line(double x, double y,22 double x0, double y0,23 double x1, double y1) {21 int __point_on_line(double x, double y, 22 double x0, double y0, 23 double x1, double y1) { 24 24 /*Determine whether a point is on a line segment 25 25 … … 41 41 b1 = y1 - y0; 42 42 43 if ( a_normal0*b0 + a_normal1*b1 == 0 ) {43 if ( a_normal0*b0 + a_normal1*b1 == 0.0 ) { 44 44 //Point is somewhere on the infinite extension of the line 45 // FIXME (Ole): Perhaps add a tolerance here instead of 0.0 45 46 46 47 len_a = sqrt(a0*a0 + a1*a1); … … 59 60 60 61 61 int _ separate_points_by_polygon(int M, // Number of points62 int __separate_points_by_polygon(int M, // Number of points 62 63 int N, // Number of polygon vertices 63 64 double* points, … … 117 118 118 119 //Check for case where point is contained in line segment 119 if (_ point_on_line(x, y, px_i, py_i, px_j, py_j)) {120 if (__point_on_line(x, y, px_i, py_i, px_j, py_j)) { 120 121 if (closed == 1) { 121 122 inside = 1; … … 149 150 150 151 // Gateways to Python 151 PyObject * point_on_line(PyObject *self, PyObject *args) {152 PyObject *_point_on_line(PyObject *self, PyObject *args) { 152 153 // 153 154 // point_on_line(x, y, x0, y0, x1, y1) … … 167 168 168 169 // Call underlying routine 169 res = _ point_on_line(x, y, x0, y0, x1, y1);170 res = __point_on_line(x, y, x0, y0, x1, y1); 170 171 171 172 // Return values a and b … … 176 177 177 178 178 PyObject * separate_points_by_polygon(PyObject *self, PyObject *args) {179 PyObject *_separate_points_by_polygon(PyObject *self, PyObject *args) { 179 180 //def separate_points_by_polygon(points, polygon, closed, verbose, one_point): 180 181 // """Determine whether points are inside or outside a polygon … … 245 246 246 247 //Call underlying routine 247 count = _ separate_points_by_polygon(M, N,248 (double*) points -> data,249 (double*) polygon -> data,250 (long*) indices -> data,251 closed, verbose);248 count = __separate_points_by_polygon(M, N, 249 (double*) points -> data, 250 (double*) polygon -> data, 251 (long*) indices -> data, 252 closed, verbose); 252 253 253 254 //NOTE: return number of points inside.. … … 264 265 */ 265 266 266 {" point_on_line",point_on_line, METH_VARARGS, "Print out"},267 {" separate_points_by_polygon",separate_points_by_polygon,267 {"_point_on_line", _point_on_line, METH_VARARGS, "Print out"}, 268 {"_separate_points_by_polygon", _separate_points_by_polygon, 268 269 METH_VARARGS, "Print out"}, 269 270 {NULL, NULL, 0, NULL} /* sentinel */ -
anuga_core/source/anuga/utilities/test_polygon.py
r5119 r5225 38 38 39 39 40 # Polygon stuff40 # Polygon stuff 41 41 def test_polygon_function_constants(self): 42 42 p1 = [[0,0], [10,0], [10,10], [0,10]] … … 144 144 def test_point_on_line(self): 145 145 146 # Endpoints first147 assert point_on_line( 0, 0, 0,0, 1,0)148 assert point_on_line( 1, 0, 0,0, 1,0)149 150 # Then points on line151 assert point_on_line( 0.5, 0, 0,0, 1,0)152 assert point_on_line( 0, 0.5, 0,1, 0,0)153 assert point_on_line( 1, 0.5, 1,1, 1,0)154 assert point_on_line( 0.5, 0.5, 0,0, 1,1)155 156 # Then points not on line157 assert not point_on_line( 0.5, 0, 0,1, 1,1)158 assert not point_on_line( 0, 0.5, 0,0, 1,1)159 160 # From real example that failed161 assert not point_on_line( 40,50, 40,20, 40,40)162 163 164 # From real example that failed165 assert not point_on_line( 40,19, 40,20, 40,40)146 # Endpoints first 147 assert point_on_line( [0, 0], [[0,0], [1,0]] ) 148 assert point_on_line( [1, 0], [[0,0], [1,0]] ) 149 150 # Then points on line 151 assert point_on_line( [0.5, 0], [[0,0], [1,0]] ) 152 assert point_on_line( [0, 0.5], [[0,1], [0,0]] ) 153 assert point_on_line( [1, 0.5], [[1,1], [1,0]] ) 154 assert point_on_line( [0.5, 0.5], [[0,0], [1,1]] ) 155 156 # Then points not on line 157 assert not point_on_line( [0.5, 0], [[0,1], [1,1]] ) 158 assert not point_on_line( [0, 0.5], [[0,0], [1,1]] ) 159 160 # From real example that failed 161 assert not point_on_line( [40,50], [[40,20], [40,40]] ) 162 163 164 # From real example that failed 165 assert not point_on_line( [40,19], [[40,20], [40,40]] ) 166 166 167 167 … … 171 171 172 172 173 # Simplest case: Polygon is the unit square173 # Simplest case: Polygon is the unit square 174 174 polygon = [[0,0], [1,0], [1,1], [0,1]] 175 175 … … 180 180 assert not is_inside_polygon( (1.5, 0.5), polygon ) 181 181 182 # Try point on borders182 # Try point on borders 183 183 assert is_inside_polygon( (1., 0.5), polygon, closed=True) 184 184 assert is_inside_polygon( (0.5, 1), polygon, closed=True) … … 194 194 def test_inside_polygon_main(self): 195 195 196 # Simplest case: Polygon is the unit square196 # Simplest case: Polygon is the unit square 197 197 polygon = [[0,0], [1,0], [1,1], [0,1]] 198 198 199 # From real example (that failed)199 # From real example (that failed) 200 200 polygon = [[20,20], [40,20], [40,40], [20,40]] 201 201 points = [ [40, 50] ] … … 211 211 212 212 213 # More convoluted and non convex polygon213 # More convoluted and non convex polygon 214 214 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 215 215 assert is_inside_polygon( (0.5, 0.5), polygon ) … … 221 221 222 222 223 # Very convoluted polygon223 # Very convoluted polygon 224 224 polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] 225 225 assert is_inside_polygon( (5, 5), polygon ) … … 233 233 234 234 235 # Another combination (that failed)235 # Another combination (that failed) 236 236 polygon = [[0,0], [10,0], [10,10], [0,10]] 237 237 assert is_inside_polygon( (5, 5), polygon ) … … 246 246 247 247 248 # Multiple polygons248 # Multiple polygons 249 249 250 250 polygon = [[0,0], [1,0], [1,1], [0,1], [0,0], … … 256 256 assert not is_inside_polygon( (0, 5.5), polygon ) 257 257 258 # Polygon with a hole258 # Polygon with a hole 259 259 polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1], 260 260 [0,0], [1,0], [1,1], [0,1], [0,0]] … … 268 268 269 269 270 # Simplest case: Polygon is the unit square270 # Simplest case: Polygon is the unit square 271 271 polygon = [[0,0], [1,0], [1,0], [1,0], [1,1], [0,1], [0,0]] 272 272 … … 277 277 assert not is_inside_polygon( (1.5, 0.5), polygon ) 278 278 279 # Try point on borders279 # Try point on borders 280 280 assert is_inside_polygon( (1., 0.5), polygon, closed=True) 281 281 assert is_inside_polygon( (0.5, 1), polygon, closed=True) … … 288 288 assert not is_inside_polygon( (1., 0.5), polygon, closed=False) 289 289 290 # From real example290 # From real example 291 291 polygon = [[20,20], [40,20], [40,40], [20,40]] 292 292 points = [ [40, 50] ] … … 297 297 298 298 def test_inside_polygon_vector_version(self): 299 # Now try the vector formulation returning indices299 # Now try the vector formulation returning indices 300 300 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 301 301 points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] … … 308 308 309 309 assert not is_outside_polygon( [0.5, 0.5], U ) 310 # evaluate to False as the point 0.5, 0.5 is inside the unit square310 # evaluate to False as the point 0.5, 0.5 is inside the unit square 311 311 312 312 assert is_outside_polygon( [1.5, 0.5], U ) 313 # evaluate to True as the point 1.5, 0.5 is outside the unit square313 # evaluate to True as the point 1.5, 0.5 is outside the unit square 314 314 315 315 indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U ) 316 316 assert allclose( indices, [1] ) 317 317 318 # One more test of vector formulation returning indices318 # One more test of vector formulation returning indices 319 319 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 320 320 points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] … … 335 335 336 336 assert not outside_polygon( [0.5, 1.0], U, closed = True ) 337 # evaluate to False as the point 0.5, 1.0 is inside the unit square337 # evaluate to False as the point 0.5, 1.0 is inside the unit square 338 338 339 339 assert is_outside_polygon( [0.5, 1.0], U, closed = False ) 340 # evaluate to True as the point 0.5, 1.0 is outside the unit square340 # evaluate to True as the point 0.5, 1.0 is outside the unit square 341 341 342 342 def test_all_outside_polygon(self): … … 577 577 assert (inside, [1,3]) 578 578 assert (outside, [0,2,4,5]) 579 580 def zzztest_inside_polygon_main(self): 579 580 581 def test_intersection1(self): 582 line0 = [[-1,0], [1,0]] 583 line1 = [[0,-1], [0,1]] 584 585 p = intersection(line0, line1) 586 assert allclose(p, [0.0, 0.0]) 587 588 def test_intersection2(self): 589 line0 = [[0,0], [24,12]] 590 line1 = [[0,12], [24,0]] 591 592 p = intersection(line0, line1) 593 assert allclose(p, [12.0, 6.0]) 594 595 # Swap direction of one line 596 line1 = [[24,0], [0,12]] 597 598 p = intersection(line0, line1) 599 assert allclose(p, [12.0, 6.0]) 600 601 # Swap order of lines 602 p = intersection(line1, line0) 603 assert allclose(p, [12.0, 6.0]) 604 605 def test_intersection3(self): 606 line0 = [[0,0], [24,12]] 607 line1 = [[0,17], [24,0]] 608 609 p = intersection(line0, line1) 610 #print p 611 assert allclose(p, [14.068965517, 7.0344827586]) 612 613 # Swap direction of one line 614 line1 = [[24,0], [0,17]] 615 616 p = intersection(line0, line1) 617 #print p 618 assert allclose(p, [14.068965517, 7.0344827586]) 619 620 # Swap order of lines 621 p = intersection(line1, line0) 622 assert allclose(p, [14.068965517, 7.0344827586]) 623 624 625 def test_intersection4(self): 626 line0 = [[0,0], [24,12]] 627 line1 = [[0,22], [21,0]] 628 629 p = intersection(line0, line1) 630 #print 'P',p 631 632 633 def test_intersection_direction_invariance(self): 634 """This runs through a number of examples and checks that direction of lines don't matter. 635 """ 636 637 638 line0 = [[0,0], [100,100]] 639 640 common_end_point = [20, 150] 641 642 for i in range(100): 643 x = 20 + i * 1.0/100 644 645 line1 = [[x,0], common_end_point] 646 p1 = intersection(line0, line1) 647 648 # Swap direction of line1 649 line1 = [common_end_point, [x,0]] 650 p2 = intersection(line0, line1) 651 652 msg = 'Orientation of line shouldn not matter.\n' 653 msg += 'However, segment [%f,%f], [%f, %f]' %(x, 654 0, 655 common_end_point[0], 656 common_end_point[1]) 657 msg += ' gave %s, \nbut when reversed we got %s' %(p1, p2) 658 assert allclose(p1, p2), msg 659 660 # Swap order of lines 661 p3 = intersection(line1, line0) 662 msg = 'Order of lines gave different results' 663 assert allclose(p1, p3), msg 664 665 def test_no_intersection(self): 666 line0 = [[-1,1], [1,1]] 667 line1 = [[0,-1], [0,0]] 668 669 p = intersection(line0, line1) 670 assert p is None 671 672 def test_intersection_parallel(self): 673 line0 = [[-1,1], [1,1]] 674 line1 = [[-1,0], [5,0]] 675 676 p = intersection(line0, line1) 677 assert p is None 678 679 680 681 def zzztest_inside_polygon_main(self): \ 682 683 #FIXME (Ole): Why is this disabled? 581 684 print "inside",inside 582 685 print "outside",outside
Note: See TracChangeset
for help on using the changeset viewer.