Ignore:
Timestamp:
Jun 18, 2010, 4:43:10 PM (14 years ago)
Author:
hudson
Message:

Refactorings to increase code quality, fixed missing log import from sww_interrogate.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/geometry/polygon.py

    r7841 r7858  
    1010import anuga.utilities.log as log
    1111
     12from aabb import AABB
    1213
    1314##
     
    3637
    3738    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],
    4041                         rtol, atol)
    4142
     
    4950
    5051# 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]))
     52def lines_dont_coincide(p0, p1, p2, p3):
     53    return (3, None)
     54   
     55def lines_0_fully_included_in_1(p0, p1, p2, p3):
     56    return (2, num.array([p0, p1]))
     57   
     58def lines_1_fully_included_in_0(p0, p1, p2, p3):
     59    return (2, num.array([p2, p3]))
     60   
     61def lines_overlap_same_direction(p0, p1, p2, p3):
     62    return (2, num.array([p0, p3]))
     63   
     64def lines_overlap_same_direction2(p0, p1, p2, p3):
     65    return (2, num.array([p2, p1]))
     66   
     67def lines_overlap_opposite_direction(p0, p1, p2, p3):
     68    return (2, num.array([p0, p2]))
     69   
     70def lines_overlap_opposite_direction2(p0, p1, p2, p3):
     71    return (2, num.array([p3, p1]))
    6472
    6573# this function called when an impossible state is found
     
    146154                           point_on_line([x3, y3], line0, rtol=rtol, atol=atol))
    147155
    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])
    150158        else:
    151159            # Lines are parallel but aren't collinear
     
    269277    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
    270278   
    271    
    272 
    273     # FIXME (Ole): The rest of this function has been made
    274     # obsolete by the C extension.
    275    
    276     # Quickly reject points that are clearly outside
    277     if point[0] < min(triangle[:, 0]):
    278         return False
    279     if point[0] > max(triangle[:, 0]):
    280         return False   
    281     if point[1] < min(triangle[:, 1]):
    282         return False
    283     if point[1] > max(triangle[:, 1]):
    284         return False       
    285 
    286 
    287     # Start search   
    288     A = triangle[0, :]
    289     B = triangle[1, :]
    290     C = triangle[2, :]
    291    
    292     # Now check if point lies wholly inside triangle
    293     v0 = C-A
    294     v1 = B-A       
    295        
    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*a10
    301    
    302     if abs(denom) > 0.0:
    303         v = point-A
    304         b0 = num.inner(v0, v)       
    305         b1 = num.inner(v1, v)           
    306    
    307         alpha = (b0*a11 - b1*a01)/denom
    308         beta = (b1*a00 - b0*a10)/denom       
    309 
    310         if (alpha > 0.0) and (beta > 0.0) and (alpha+beta < 1.0):
    311             return True
    312 
    313 
    314     if closed is True:
    315         # Check if point lies on one of the edges
    316        
    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 True
    325                
    326     return False
    327 
    328    
    329279def is_complex(polygon, verbose=False):
    330280    """Check if a polygon is complex (self-intersecting).
     
    339289           
    340290    def key_xpos(item):
     291        """ Return the x coord out of the passed point for sorting key. """
    341292        return (item[0][0])
    342293   
     
    377328                        print 'Self-intersecting polygon found, type ', type
    378329                        print 'point', point,
    379                         print 'vertices: ', leftmost, ' - ', l_x[cmp]               
     330                        print 'vertices: ', leftmost, ' - ', l_x[cmp] 
    380331                    return True           
    381332            cmp += 1
     
    422373    try:
    423374        points = ensure_absolute(points)
    424     except NameError, e:
    425         raise NameError, e
     375    except NameError, err:
     376        raise NameError, err
    426377    except:
    427378        # If this fails it is going to be because the points can't be
     
    513464    if len(points.shape) == 1:
    514465        # 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))
    516467
    517468    indices, count = separate_points_by_polygon(points, polygon,
     
    559510    if len(points.shape) == 1:
    560511        # 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))
    562513
    563514    indices, count = separate_points_by_polygon(points, polygon,
     
    632583        except:
    633584            msg = 'Points could not be converted to numeric array'
    634             raise msg
     585            raise Exception(msg)
    635586
    636587        try:
    637588            polygon = ensure_numeric(polygon, num.float)
    638589        except NameError, e:
    639             raise NameError, e
     590            raise NameError(e)
    640591        except:
    641592            msg = 'Polygon could not be converted to numeric array'
    642             raise msg
     593            raise Exception(msg)
    643594
    644595        msg = 'Polygon array must be a 2d array of vertices'
     
    646597
    647598        msg = 'Polygon array must have two columns'
    648         assert polygon.shape[1]==2, msg
     599        assert polygon.shape[1] == 2, msg
    649600
    650601        msg = ('Points array must be 1 or 2 dimensional. '
     
    654605        if len(points.shape) == 1:
    655606            # 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))
    657608   
    658609            msg = ('Point array must have two columns (x,y), '
     
    663614            msg = ('Points array must be a 2d array. I got %s.'
    664615                   % str(points[:30]))
    665             assert len(points.shape)==2, msg
     616            assert len(points.shape) == 2, msg
    666617
    667618            msg = 'Points array must have two columns'
    668             assert points.shape[1]==2, msg
     619            assert points.shape[1] == 2, msg
    669620
    670621    N = polygon.shape[0] # Number of vertices in polygon
     
    681632    return indices, count
    682633
    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
    687635def polygon_area(input_polygon):
    688636    """ Determine area of arbitrary polygon.
    689637
    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    """
    693644    # Move polygon to origin (0,0) to avoid rounding errors
    694645    # This makes a copy of the polygon to avoid destroying it
    695646    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])
    698649    polygon = input_polygon - [min_x, min_y]
    699650
     
    742693    Outputs:
    743694
    744     - list of min and max of x and y coordinates
    745695    - plot of polygons
    746696    """
     
    754704    ion()
    755705    hold(True)
    756 
    757     minx = 1e10
    758     maxx = 0.0
    759     miny = 1e10
    760     maxy = 0.0
    761706
    762707    if label is None:
     
    772717            alpha = max(0.0, min(1.0, alpha))
    773718
    774     n = len(polygons_points)
     719    num_points = len(polygons_points)
    775720    colour = []
    776721    if style is None:
    777722        style_type = 'line'
    778723        style = []
    779         for i in range(n):
     724        for i in range(num_points):
    780725            style.append(style_type)
    781726            colour.append('b-')
    782727    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)
    791737
    792738    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])
    803741        if alpha:
    804             fill(x, y, colour[i], alpha=alpha)
     742            fill(pt_x, pt_y, colour[i], alpha=alpha)
    805743        xlabel('x')
    806744        ylabel('y')
     
    814752    close('all')
    815753
    816     vec = [minx, maxx, miny, maxy]
    817     return vec
    818 
    819 
    820 def poly_xy(polygon):
     754
     755def _poly_xy(polygon):
    821756    """ this is used within plot_polygons so need to duplicate
    822757        the first point so can have closed polygon in plot
     
    829764    try:
    830765        polygon = ensure_numeric(polygon, num.float)
    831     except NameError, e:
    832         raise NameError, e
     766    except NameError, err:
     767        raise NameError, err
    833768    except:
    834769        msg = ('Polygon %s could not be converted to numeric array'
     
    836771        raise Exception, msg
    837772
    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 one-dimensional 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
    969778
    970779################################################################################
     
    973782
    974783def 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 (self-intersecting).
     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.
    984794    """
    985795
     
    1002812    return polygon
    1003813
    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
    1009815def write_polygon(polygon, filename=None):
    1010816    """Write polygon to csv file.
     
    1047853
    1048854    # 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   
    1059857    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)
    1062860
    1063861        append = False
    1064         if is_inside_polygon([x,y], polygon):
     862        if is_inside_polygon([rand_x, rand_y], polygon):
    1065863            append = True
    1066864
     
    1068866            if exclude is not None:
    1069867                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):
    1071869                        append = False
    1072870
    1073871        if append is True:
    1074             points.append([x, y])
     872            points.append([rand_x, rand_y])
    1075873
    1076874    return points
     
    1091889    """
    1092890
    1093     import exceptions
    1094 
    1095     class Found(exceptions.Exception):
    1096         pass
    1097 
    1098891    polygon = ensure_numeric(polygon)
    1099892   
    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
    1129914
    1130915
     
    1212997    reduced_polygon = []
    1213998    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]:
    12171000            # Keep
    12181001            reduced_polygon.append(point)
     
    13181101
    13191102else:
    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)
    13241107
    13251108
Note: See TracChangeset for help on using the changeset viewer.