Changeset 5225


Ignore:
Timestamp:
Apr 21, 2008, 5:20:32 PM (17 years ago)
Author:
ole
Message:

Work done during Water Down Under 2008.
Line Mesh intersections towards ticket:175 (flow through a cross section).

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

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/utilities/polygon.py

    r5114 r5225  
    1010#    #print 'Could not find scipy - using Numeric'
    1111
    12 from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot
     12from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot, allclose
    1313
    1414
     
    1717from anuga.geospatial_data.geospatial_data import ensure_absolute
    1818
    19 def point_on_line(x, y, x0, y0, x1, y1):
     19def point_on_line(point, line):
    2020    """Determine whether a point is on a line segment
    2121
    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
     45def 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
    4277    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
    44105
    45106
     
    305366    indices = zeros( M, Int )
    306367
    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))
    310370
    311371    if verbose: print 'Found %d points (out of %d) inside polygon'\
     
    780840from anuga.utilities.compile import can_use_C_extension
    781841if 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
    785846else:
    786847    msg = 'C implementations could not be accessed by %s.\n ' %__file__
  • anuga_core/source/anuga/utilities/polygon_ext.c

    r4804 r5225  
    1919
    2020
    21 int _point_on_line(double x, double y,
    22                    double x0, double y0,
    23                    double x1, double y1) {
     21int __point_on_line(double x, double y,
     22                    double x0, double y0,
     23                    double x1, double y1) {
    2424  /*Determine whether a point is on a line segment
    2525
     
    4141  b1 = y1 - y0;
    4242
    43   if ( a_normal0*b0 + a_normal1*b1 == 0 ) {
     43  if ( a_normal0*b0 + a_normal1*b1 == 0.0 ) {
    4444    //Point is somewhere on the infinite extension of the line
     45    // FIXME (Ole): Perhaps add a tolerance here instead of 0.0
    4546
    4647    len_a = sqrt(a0*a0 + a1*a1);
     
    5960
    6061
    61 int _separate_points_by_polygon(int M,     // Number of points
     62int __separate_points_by_polygon(int M,     // Number of points
    6263                                int N,     // Number of polygon vertices
    6364                                double* points,
     
    117118
    118119        //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)) {
    120121          if (closed == 1) {
    121122            inside = 1;
     
    149150
    150151// Gateways to Python
    151 PyObject *point_on_line(PyObject *self, PyObject *args) {
     152PyObject *_point_on_line(PyObject *self, PyObject *args) {
    152153  //
    153154  // point_on_line(x, y, x0, y0, x1, y1)
     
    167168
    168169  // 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);
    170171
    171172  // Return values a and b
     
    176177
    177178
    178 PyObject *separate_points_by_polygon(PyObject *self, PyObject *args) {
     179PyObject *_separate_points_by_polygon(PyObject *self, PyObject *args) {
    179180  //def separate_points_by_polygon(points, polygon, closed, verbose, one_point):
    180181  //  """Determine whether points are inside or outside a polygon
     
    245246 
    246247  //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);
    252253 
    253254  //NOTE: return number of points inside..
     
    264265   */
    265266
    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,
    268269                                 METH_VARARGS, "Print out"},
    269270  {NULL, NULL, 0, NULL}   /* sentinel */
  • anuga_core/source/anuga/utilities/test_polygon.py

    r5119 r5225  
    3838
    3939
    40     #Polygon stuff
     40    # Polygon stuff
    4141    def test_polygon_function_constants(self):
    4242        p1 = [[0,0], [10,0], [10,10], [0,10]]
     
    144144    def test_point_on_line(self):
    145145
    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 )
     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]] )
    166166
    167167
     
    171171
    172172
    173         #Simplest case: Polygon is the unit square
     173        # Simplest case: Polygon is the unit square
    174174        polygon = [[0,0], [1,0], [1,1], [0,1]]
    175175
     
    180180        assert not is_inside_polygon( (1.5, 0.5), polygon )
    181181
    182         #Try point on borders
     182        # Try point on borders
    183183        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
    184184        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
     
    194194    def test_inside_polygon_main(self):
    195195
    196         #Simplest case: Polygon is the unit square
     196        # Simplest case: Polygon is the unit square
    197197        polygon = [[0,0], [1,0], [1,1], [0,1]]       
    198198
    199         #From real example (that failed)
     199        # From real example (that failed)
    200200        polygon = [[20,20], [40,20], [40,40], [20,40]]
    201201        points = [ [40, 50] ]
     
    211211
    212212
    213         #More convoluted and non convex polygon
     213        # More convoluted and non convex polygon
    214214        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    215215        assert is_inside_polygon( (0.5, 0.5), polygon )
     
    221221
    222222
    223         #Very convoluted polygon
     223        # Very convoluted polygon
    224224        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
    225225        assert is_inside_polygon( (5, 5), polygon )
     
    233233
    234234
    235         #Another combination (that failed)
     235        # Another combination (that failed)
    236236        polygon = [[0,0], [10,0], [10,10], [0,10]]
    237237        assert is_inside_polygon( (5, 5), polygon )
     
    246246
    247247
    248         #Multiple polygons
     248        # Multiple polygons
    249249
    250250        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
     
    256256        assert not is_inside_polygon( (0, 5.5), polygon )
    257257
    258         #Polygon with a hole
     258        # Polygon with a hole
    259259        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
    260260                   [0,0], [1,0], [1,1], [0,1], [0,0]]
     
    268268
    269269
    270         #Simplest case: Polygon is the unit square
     270        # Simplest case: Polygon is the unit square
    271271        polygon = [[0,0], [1,0], [1,0], [1,0], [1,1], [0,1], [0,0]]
    272272
     
    277277        assert not is_inside_polygon( (1.5, 0.5), polygon )
    278278
    279         #Try point on borders
     279        # Try point on borders
    280280        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
    281281        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
     
    288288        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
    289289
    290         #From real example
     290        # From real example
    291291        polygon = [[20,20], [40,20], [40,40], [20,40]]
    292292        points = [ [40, 50] ]
     
    297297
    298298    def test_inside_polygon_vector_version(self):
    299         #Now try the vector formulation returning indices
     299        # Now try the vector formulation returning indices
    300300        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    301301        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     
    308308
    309309        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 square
     310        # evaluate to False as the point 0.5, 0.5 is inside the unit square
    311311       
    312312        assert is_outside_polygon( [1.5, 0.5], U )
    313         #evaluate to True as the point 1.5, 0.5 is outside the unit square
     313        # evaluate to True as the point 1.5, 0.5 is outside the unit square
    314314       
    315315        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
    316316        assert allclose( indices, [1] )
    317317       
    318         #One more test of vector formulation returning indices
     318        # One more test of vector formulation returning indices
    319319        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    320320        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     
    335335   
    336336        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 square
     337        # evaluate to False as the point 0.5, 1.0 is inside the unit square
    338338       
    339339        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 square
     340        # evaluate to True as the point 0.5, 1.0 is outside the unit square
    341341
    342342    def test_all_outside_polygon(self):
     
    577577        assert (inside, [1,3])
    578578        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?
    581684        print "inside",inside
    582685        print "outside",outside
Note: See TracChangeset for help on using the changeset viewer.