source: trunk/anuga_core/source/anuga/geometry/polygon.py @ 7858

Last change on this file since 7858 was 7858, checked in by hudson, 14 years ago

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

File size: 37.2 KB
RevLine 
[5897]1#!/usr/bin/env python
2
[7276]3"""Polygon manipulations"""
[5897]4
[7276]5import numpy as num
6
[5897]7from anuga.utilities.numerical_tools import ensure_numeric
[7778]8from anuga.geospatial_data.geospatial_data import ensure_absolute, \
9                                                    Geospatial_data
[7317]10import anuga.utilities.log as log
[5897]11
[7858]12from aabb import AABB
[5897]13
[7276]14##
15# @brief Determine whether a point is on a line segment.
16# @param point (x, y) of point in question (tuple, list or array).
17# @param line ((x1,y1), (x2,y2)) for line (tuple, list or array).
18# @param rtol Relative error for 'close'.
19# @param atol Absolute error for 'close'.
20# @return True or False.
[5932]21def point_on_line(point, line, rtol=1.0e-5, atol=1.0e-8):
[5897]22    """Determine whether a point is on a line segment
23
[7276]24    Input:
[5897]25        point is given by [x, y]
[7276]26        line is given by [x0, y0], [x1, y1]] or
27        the equivalent 2x2 numeric array with each row corresponding to a point.
[5897]28
29    Output:
30
[7276]31    Note: Line can be degenerate and function still works to discern coinciding
32          points from non-coinciding.
[5897]33    """
34
35    point = ensure_numeric(point)
36    line = ensure_numeric(line)
37
38    res = _point_on_line(point[0], point[1],
[7858]39                         line[0, 0], line[0, 1],
40                         line[1, 0], line[1, 1],
[5897]41                         rtol, atol)
[7276]42
[5897]43    return bool(res)
44
45
[5942]46######
47# Result functions used in intersection() below for collinear lines.
48# (p0,p1) defines line 0, (p2,p3) defines line 1.
49######
[5897]50
[5942]51# result functions for possible states
[7858]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]))
[5897]72
[5942]73# this function called when an impossible state is found
[7276]74def lines_error(p1, p2, p3, p4):
75    raise RuntimeError, ('INTERNAL ERROR: p1=%s, p2=%s, p3=%s, p4=%s'
76                         % (str(p1), str(p2), str(p3), str(p4)))
[5897]77
[7778]78collinear_result = {
79# line 0 starts on 1, 0 ends 1, 1 starts 0, 1 ends 0
80#       0s1    0e1    1s0    1e0   
81       (False, False, False, False): lines_dont_coincide,
82       (False, False, False, True ): lines_error,
83       (False, False, True,  False): lines_error,
84       (False, False, True,  True ): lines_1_fully_included_in_0,
85       (False, True,  False, False): lines_error,
86       (False, True,  False, True ): lines_overlap_opposite_direction2,
87       (False, True,  True,  False): lines_overlap_same_direction2,
88       (False, True,  True,  True ): lines_1_fully_included_in_0,
89       (True,  False, False, False): lines_error,
90       (True,  False, False, True ): lines_overlap_same_direction,
91       (True,  False, True,  False): lines_overlap_opposite_direction,
92       (True,  False, True,  True ): lines_1_fully_included_in_0,
93       (True,  True,  False, False): lines_0_fully_included_in_1,
94       (True,  True,  False, True ): lines_0_fully_included_in_1,
95       (True,  True,  True,  False): lines_0_fully_included_in_1,
96       (True,  True,  True,  True ): lines_0_fully_included_in_1
97   }
[5942]98
[7276]99##
100# @brief Finds intersection point of two line segments.
101# @param line0 First line ((x1,y1), (x2,y2)).
102# @param line1 Second line ((x1,y1), (x2,y2)).
103# @param rtol Relative error for 'close'.
104# @param atol Absolute error for 'close'.
105# @return (status, value) where:
106#             status = 0 - no intersection, value set to None
107#                      1 - intersection found, value=(x,y)
108#                      2 - lines collienar, overlap, value=overlap segment
109#                      3 - lines collinear, no overlap, value is None
110#                      4 - lines parallel, value is None
[5932]111def intersection(line0, line1, rtol=1.0e-5, atol=1.0e-8):
[7276]112    """Returns intersecting point between two line segments.
[5897]113
[7276]114    However, if parallel lines coincide partly (i.e. share a common segment),
[5897]115    the line segment where lines coincide is returned
116
117    Inputs:
118        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
[5942]119                      A line can also be a 2x2 numpy array with each row
[5897]120                      corresponding to a point.
121
122    Output:
[7276]123        status, value - where status and value is interpreted as follows:
[5942]124        status == 0: no intersection, value set to None.
125        status == 1: intersection point found and returned in value as [x,y].
[7276]126        status == 2: Collinear overlapping lines found.
127                     Value takes the form [[x0,y0], [x1,y1]].
[5942]128        status == 3: Collinear non-overlapping lines. Value set to None.
[7276]129        status == 4: Lines are parallel. Value set to None.
[5897]130    """
131
132    # FIXME (Ole): Write this in C
133
[7276]134    line0 = ensure_numeric(line0, num.float)
135    line1 = ensure_numeric(line1, num.float)
[5897]136
[7841]137    x0 = line0[0, 0]; y0 = line0[0, 1]
138    x1 = line0[1, 0]; y1 = line0[1, 1]
[5897]139
[7841]140    x2 = line1[0, 0]; y2 = line1[0, 1]
141    x3 = line1[1, 0]; y3 = line1[1, 1]
[5897]142
143    denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
144    u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
145    u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
[7276]146
[6158]147    if num.allclose(denom, 0.0, rtol=rtol, atol=atol):
[5942]148        # Lines are parallel - check if they are collinear
[6158]149        if num.allclose([u0, u1], 0.0, rtol=rtol, atol=atol):
[5942]150            # We now know that the lines are collinear
151            state_tuple = (point_on_line([x0, y0], line1, rtol=rtol, atol=atol),
152                           point_on_line([x1, y1], line1, rtol=rtol, atol=atol),
153                           point_on_line([x2, y2], line0, rtol=rtol, atol=atol),
154                           point_on_line([x3, y3], line0, rtol=rtol, atol=atol))
[5897]155
[7858]156            return collinear_result[state_tuple]([x0, y0], [x1, y1],
157                                                 [x2, y2], [x3, y3])
[5897]158        else:
[5942]159            # Lines are parallel but aren't collinear
[7276]160            return 4, None #FIXME (Ole): Add distance here instead of None
[5897]161    else:
[5942]162        # Lines are not parallel, check if they intersect
[5897]163        u0 = u0/denom
[7276]164        u1 = u1/denom
[5897]165
166        x = x0 + u0*(x1-x0)
167        y = y0 + u0*(y1-y0)
168
169        # Sanity check - can be removed to speed up if needed
[6158]170        assert num.allclose(x, x2 + u1*(x3-x2), rtol=rtol, atol=atol)
[7276]171        assert num.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol)
[5897]172
173        # Check if point found lies within given line segments
[7276]174        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0:
[5897]175            # We have intersection
[6158]176            return 1, num.array([x, y])
[5897]177        else:
178            # No intersection
179            return 0, None
180
[7276]181##
182# @brief Finds intersection point of two line segments.
183# @param line0 First line ((x1,y1), (x2,y2)).
184# @param line1 Second line ((x1,y1), (x2,y2)).
185# @return (status, value) where:
186#             status = 0 - no intersection, value set to None
187#                      1 - intersection found, value=(x,y)
188#                      2 - lines collienar, overlap, value=overlap segment
189#                      3 - lines collinear, no overlap, value is None
190#                      4 - lines parallel, value is None
191# @note Wrapper for C function.
[5897]192def NEW_C_intersection(line0, line1):
[7276]193    """Returns intersecting point between two line segments.
[5897]194
[7276]195    However, if parallel lines coincide partly (i.e. share a common segment),
[5897]196    the line segment where lines coincide is returned
197
198    Inputs:
199        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
[7276]200                      A line can also be a 2x2 numpy array with each row
[5897]201                      corresponding to a point.
202
203    Output:
[7276]204        status, value - where status and value is interpreted as follows:
205        status == 0: no intersection, value set to None.
206        status == 1: intersection point found and returned in value as [x,y].
207        status == 2: Collinear overlapping lines found.
208                     Value takes the form [[x0,y0], [x1,y1]].
209        status == 3: Collinear non-overlapping lines. Value set to None.
210        status == 4: Lines are parallel. Value set to None.
[5897]211    """
212
[7276]213    line0 = ensure_numeric(line0, num.float)
214    line1 = ensure_numeric(line1, num.float)
[5897]215
[7841]216    status, value = _intersection(line0[0, 0], line0[0, 1],
217                                  line0[1, 0], line0[1, 1],
218                                  line1[0, 0], line1[0, 1],
219                                  line1[1, 0], line1[1, 1])
[5897]220
221    return status, value
222
[6534]223def is_inside_triangle(point, triangle, 
224                       closed=True, 
[6535]225                       rtol=1.0e-12,
226                       atol=1.0e-12,                     
[7841]227                       check_inputs=True):
[6534]228    """Determine if one point is inside a triangle
229   
230    This uses the barycentric method:
231   
232    Triangle is A, B, C
233    Point P can then be written as
234   
235    P = A + alpha * (C-A) + beta * (B-A)
236    or if we let
237    v=P-A, v0=C-A, v1=B-A   
238   
239    v = alpha*v0 + beta*v1
[5897]240
[6534]241    Dot this equation by v0 and v1 to get two:
242   
243    dot(v0, v) = alpha*dot(v0, v0) + beta*dot(v0, v1)
244    dot(v1, v) = alpha*dot(v1, v0) + beta*dot(v1, v1)   
245   
246    or if a_ij = dot(v_i, v_j) and b_i = dot(v_i, v)
247    the matrix equation:
248   
249    a_00 a_01   alpha     b_0
250                       =
251    a_10 a_11   beta      b_1
252   
253    Solving for alpha and beta yields:
254   
255    alpha = (b_0*a_11 - b_1*a_01)/denom
256    beta =  (b_1*a_00 - b_0*a_10)/denom
257   
258    with denom = a_11*a_00 - a_10*a_01
259   
260    The point is in the triangle whenever
261    alpha and beta and their sums are in the unit interval.
262   
263    rtol and atol will determine how close the point has to be to the edge
[6535]264    before it is deemed to be on the edge.
[6534]265   
266    """
[5897]267
[6544]268    triangle = ensure_numeric(triangle)       
[7276]269    point = ensure_numeric(point, num.float)   
[6544]270   
[6534]271    if check_inputs is True:
272        msg = 'is_inside_triangle must be invoked with one point only'
273        assert num.allclose(point.shape, [2]), msg
274   
[6544]275   
[6541]276    # Use C-implementation
[6535]277    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
[6534]278   
[7687]279def is_complex(polygon, verbose=False):
[7690]280    """Check if a polygon is complex (self-intersecting).
281       Uses a sweep algorithm that is O(n^2) in the worst case, but
282       for most normal looking polygons it'll be O(n log n).
283
284       polygon is a list of points that define a closed polygon.
285       verbose will print a list of the intersection points if true
286       
287       Return True if polygon is complex.
288    """           
289           
290    def key_xpos(item):
[7858]291        """ Return the x coord out of the passed point for sorting key. """
[7690]292        return (item[0][0])
[7686]293   
[7690]294    def segments_joined(seg0, seg1):
[7841]295        """ See if there are identical segments in the 2 lists. """
[7690]296        for i in seg0:
297            for j in seg1:   
298                if i == j: return True
299        return False
300       
[7686]301    polygon = ensure_numeric(polygon, num.float)
302
[7690]303    # build a list of discrete segments from the polygon
304    unsorted_segs = []
[7686]305    for i in range(0, len(polygon)-1):
[7690]306        unsorted_segs.append([list(polygon[i]), list(polygon[i+1])])
307    unsorted_segs.append([list(polygon[0]), list(polygon[-1])])
308   
309    # all segments must point in same direction
310    for val in unsorted_segs:
311        if val[0][0] > val[1][0]:
312            val[0], val[1] = val[1], val[0]   
313           
314    l_x = sorted(unsorted_segs, key=key_xpos)
[7686]315
[7690]316    comparisons = 0
317   
318    # loop through, only comparing lines that partially overlap in x
319    for index, leftmost in enumerate(l_x):
320        cmp = index+1
321        while cmp < len(l_x) and leftmost[1][0] > l_x[cmp][0][0]:
322            if not segments_joined(leftmost, l_x[cmp]):
323                (type, point) = intersection(leftmost, l_x[cmp])
324                comparisons += 1
[7778]325                if type != 0 and type != 4 or (type == 2 and list(point[0]) !=\
326                                                list(point[1])):
[7687]327                    if verbose:
[7778]328                        print 'Self-intersecting polygon found, type ', type
329                        print 'point', point,
[7858]330                        print 'vertices: ', leftmost, ' - ', l_x[cmp] 
[7690]331                    return True           
332            cmp += 1
[7686]333       
334    return False
335   
[6534]336
[5897]337def is_inside_polygon(point, polygon, closed=True, verbose=False):
338    """Determine if one point is inside a polygon
339
340    See inside_polygon for more details
341    """
342
343    indices = inside_polygon(point, polygon, closed, verbose)
344
345    if indices.shape[0] == 1:
346        return True
347    elif indices.shape[0] == 0:
348        return False
349    else:
350        msg = 'is_inside_polygon must be invoked with one point only'
[7841]351        raise Exception(msg)
[5897]352
[7276]353##
354# @brief Determine which of a set of points are inside a polygon.
355# @param points A set of points (tuple, list or array).
356# @param polygon A set of points defining a polygon (tuple, list or array).
357# @param closed True if points on boundary are considered 'inside' polygon.
358# @param verbose True if this function is to be verbose.
359# @return A list of indices of points inside the polygon.
[5897]360def inside_polygon(points, polygon, closed=True, verbose=False):
361    """Determine points inside a polygon
362
363       Functions inside_polygon and outside_polygon have been defined in
[7276]364       terms of separate_by_polygon which will put all inside indices in
[5897]365       the first part of the indices array and outside indices in the last
366
367       See separate_points_by_polygon for documentation
368
369       points and polygon can be a geospatial instance,
370       a list or a numeric array
371    """
372
373    try:
374        points = ensure_absolute(points)
[7858]375    except NameError, err:
376        raise NameError, err
[5897]377    except:
378        # If this fails it is going to be because the points can't be
379        # converted to a numeric array.
[7276]380        msg = 'Points could not be converted to numeric array' 
381        raise Exception, msg
[5897]382
[6534]383    polygon = ensure_absolute(polygon)       
[5897]384    try:
385        polygon = ensure_absolute(polygon)
386    except NameError, e:
387        raise NameError, e
388    except:
389        # If this fails it is going to be because the points can't be
390        # converted to a numeric array.
[7276]391        msg = ('Polygon %s could not be converted to numeric array'
392               % (str(polygon)))
393        raise Exception, msg
[5897]394
395    if len(points.shape) == 1:
396        # Only one point was passed in. Convert to array of points
[7276]397        points = num.reshape(points, (1,2))
[5897]398
399    indices, count = separate_points_by_polygon(points, polygon,
400                                                closed=closed,
401                                                verbose=verbose)
402
403    # Return indices of points inside polygon
404    return indices[:count]
405
[7276]406##
407# @brief Determine if one point is outside a polygon.
408# @param point The point of interest.
409# @param polygon The polygon to test inclusion in.
410# @param closed True if points on boundary are considered 'inside' polygon.
411# @param verbose True if this function is to be verbose.
412# @return True if point is outside the polygon.
413# @note Uses inside_polygon() to do the work.
[5897]414def is_outside_polygon(point, polygon, closed=True, verbose=False,
415                       points_geo_ref=None, polygon_geo_ref=None):
416    """Determine if one point is outside a polygon
417
418    See outside_polygon for more details
419    """
420
421    indices = outside_polygon(point, polygon, closed, verbose)
422
423    if indices.shape[0] == 1:
424        return True
425    elif indices.shape[0] == 0:
426        return False
427    else:
428        msg = 'is_outside_polygon must be invoked with one point only'
[7276]429        raise Exception, msg
[5897]430
[7276]431##
432# @brief Determine which of a set of points are outside a polygon.
433# @param points A set of points (tuple, list or array).
434# @param polygon A set of points defining a polygon (tuple, list or array).
435# @param closed True if points on boundary are considered 'inside' polygon.
436# @param verbose True if this function is to be verbose.
437# @return A list of indices of points outside the polygon.
[5897]438def outside_polygon(points, polygon, closed = True, verbose = False):
439    """Determine points outside a polygon
440
441       Functions inside_polygon and outside_polygon have been defined in
[7276]442       terms of separate_by_polygon which will put all inside indices in
[5897]443       the first part of the indices array and outside indices in the last
444
445       See separate_points_by_polygon for documentation
446    """
447
448    try:
[7276]449        points = ensure_numeric(points, num.float)
[5897]450    except NameError, e:
451        raise NameError, e
452    except:
[7276]453        msg = 'Points could not be converted to numeric array'
454        raise Exception, msg
[5897]455
456    try:
[7276]457        polygon = ensure_numeric(polygon, num.float)
[5897]458    except NameError, e:
459        raise NameError, e
460    except:
[7276]461        msg = 'Polygon could not be converted to numeric array'
462        raise Exception, msg
[5897]463
464    if len(points.shape) == 1:
465        # Only one point was passed in. Convert to array of points
[7858]466        points = num.reshape(points, (1, 2))
[5897]467
468    indices, count = separate_points_by_polygon(points, polygon,
469                                                closed=closed,
470                                                verbose=verbose)
471
472    # Return indices of points outside polygon
473    if count == len(indices):
474        # No points are outside
[6158]475        return num.array([])
[5897]476    else:
477        return indices[count:][::-1]  #return reversed
478
[7276]479##
480# @brief Separate a list of points into two sets inside+outside a polygon.
481# @param points A set of points (tuple, list or array).
482# @param polygon A set of points defining a polygon (tuple, list or array).
483# @param closed True if points on boundary are considered 'inside' polygon.
484# @param verbose True if this function is to be verbose.
485# @return A tuple (in, out) of point indices for poinst inside amd outside.
[6534]486def in_and_outside_polygon(points, polygon, closed=True, verbose=False):
[5897]487    """Determine points inside and outside a polygon
488
489       See separate_points_by_polygon for documentation
490
[7276]491       Returns an array of points inside and array of points outside the polygon
[5897]492    """
493
494    try:
[7276]495        points = ensure_numeric(points, num.float)
[5897]496    except NameError, e:
497        raise NameError, e
498    except:
[7276]499        msg = 'Points could not be converted to numeric array'
500        raise Exception, msg
[5897]501
502    try:
[7276]503        polygon = ensure_numeric(polygon, num.float)
[5897]504    except NameError, e:
505        raise NameError, e
506    except:
[7276]507        msg = 'Polygon could not be converted to numeric array'
508        raise Exception, msg
[5897]509
510    if len(points.shape) == 1:
511        # Only one point was passed in. Convert to array of points
[7858]512        points = num.reshape(points, (1, 2))
[5897]513
514    indices, count = separate_points_by_polygon(points, polygon,
515                                                closed=closed,
516                                                verbose=verbose)
[7276]517
[5897]518    # Returns indices of points inside and indices of points outside
519    # the polygon
520    if count == len(indices):
521        # No points are outside
[7778]522        return indices[:count], []
[5897]523    else:
524        return  indices[:count], indices[count:][::-1]  #return reversed
525
[7778]526
527
[5897]528def separate_points_by_polygon(points, polygon,
[6534]529                               closed=True, 
530                               check_input=True,
531                               verbose=False):
[5897]532    """Determine whether points are inside or outside a polygon
533
534    Input:
535       points - Tuple of (x, y) coordinates, or list of tuples
536       polygon - list of vertices of polygon
537       closed - (optional) determine whether points on boundary should be
538       regarded as belonging to the polygon (closed = True)
539       or not (closed = False)
[6534]540       check_input: Allows faster execution if set to False
[5897]541
542    Outputs:
543       indices: array of same length as points with indices of points falling
544       inside the polygon listed from the beginning and indices of points
545       falling outside listed from the end.
546
547       count: count of points falling inside the polygon
548
549       The indices of points inside are obtained as indices[:count]
550       The indices of points outside are obtained as indices[count:]
551
552    Examples:
553       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
554
555       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
556       will return the indices [0, 2, 1] and count == 2 as only the first
557       and the last point are inside the unit square
558
559    Remarks:
560       The vertices may be listed clockwise or counterclockwise and
561       the first point may optionally be repeated.
562       Polygons do not need to be convex.
563       Polygons can have holes in them and points inside a hole is
564       regarded as being outside the polygon.
565
566    Algorithm is based on work by Darel Finley,
567    http://www.alienryderflex.com/polygon/
568
569    Uses underlying C-implementation in polygon_ext.c
570    """
571
[6534]572    if check_input:
573        #Input checks
[7778]574        assert isinstance(closed, bool), \
575                    'Keyword argument "closed" must be boolean'
576        assert isinstance(verbose, bool), \
577                    'Keyword argument "verbose" must be boolean'
[5897]578
[6534]579        try:
[7276]580            points = ensure_numeric(points, num.float)
[6534]581        except NameError, e:
582            raise NameError, e
583        except:
[7276]584            msg = 'Points could not be converted to numeric array'
[7858]585            raise Exception(msg)
[5897]586
[6534]587        try:
[7276]588            polygon = ensure_numeric(polygon, num.float)
[6534]589        except NameError, e:
[7858]590            raise NameError(e)
[6534]591        except:
[7276]592            msg = 'Polygon could not be converted to numeric array'
[7858]593            raise Exception(msg)
[5897]594
[6534]595        msg = 'Polygon array must be a 2d array of vertices'
596        assert len(polygon.shape) == 2, msg
[5897]597
[6534]598        msg = 'Polygon array must have two columns' 
[7858]599        assert polygon.shape[1] == 2, msg
[5897]600
[7276]601        msg = ('Points array must be 1 or 2 dimensional. '
602               'I got %d dimensions' % len(points.shape))
[6534]603        assert 0 < len(points.shape) < 3, msg
[5897]604
[6534]605        if len(points.shape) == 1:
[7276]606            # Only one point was passed in.  Convert to array of points.
[7858]607            points = num.reshape(points, (1, 2))
[5897]608   
[7276]609            msg = ('Point array must have two columns (x,y), '
610                   'I got points.shape[1]=%d' % points.shape[0])
611            assert points.shape[1]==2, msg
[5897]612
613       
[7276]614            msg = ('Points array must be a 2d array. I got %s.'
615                   % str(points[:30]))
[7858]616            assert len(points.shape) == 2, msg
[5897]617
[6534]618            msg = 'Points array must have two columns'
[7858]619            assert points.shape[1] == 2, msg
[5897]620
[6534]621    N = polygon.shape[0] # Number of vertices in polygon
622    M = points.shape[0]  # Number of points
[5897]623
[7276]624    indices = num.zeros(M, num.int)
[5897]625
626    count = _separate_points_by_polygon(points, polygon, indices,
627                                        int(closed), int(verbose))
628
[7276]629    if verbose:
[7317]630        log.critical('Found %d points (out of %d) inside polygon' % (count, M))
[7276]631
[5897]632    return indices, count
633
[7858]634
[7276]635def polygon_area(input_polygon):
636    """ Determine area of arbitrary polygon.
[5897]637
[7858]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
[5897]643    """
[6000]644    # Move polygon to origin (0,0) to avoid rounding errors
[6001]645    # This makes a copy of the polygon to avoid destroying it
646    input_polygon = ensure_numeric(input_polygon)
[7858]647    min_x = min(input_polygon[:, 0])
648    min_y = min(input_polygon[:, 1])
[6001]649    polygon = input_polygon - [min_x, min_y]
[6000]650
[7276]651    # Compute area
[5897]652    n = len(polygon)
653    poly_area = 0.0
654
655    for i in range(n):
656        pti = polygon[i]
657        if i == n-1:
658            pt1 = polygon[0]
659        else:
660            pt1 = polygon[i+1]
661        xi = pti[0]
662        yi1 = pt1[1]
663        xi1 = pt1[0]
664        yi = pti[1]
665        poly_area += xi*yi1 - xi1*yi
[7276]666
[5897]667    return abs(poly_area/2)
668
[7778]669
[7276]670def plot_polygons(polygons_points,
671                  style=None,
672                  figname=None,
673                  label=None,
[7841]674                  alpha=None):
[5897]675    """ Take list of polygons and plot.
676
677    Inputs:
678
679    polygons         - list of polygons
680
681    style            - style list corresponding to each polygon
682                     - for a polygon, use 'line'
683                     - for points falling outside a polygon, use 'outside'
[7511]684                     - style can also be user defined as in normal pylab plot.
[7276]685
[5897]686    figname          - name to save figure to
687
[7516]688    label            - title for plotA
[5897]689
[7516]690    alpha            - transparency of polygon fill, 0.0=none, 1.0=solid
691                       if not supplied, no fill.
692
[5897]693    Outputs:
694
695    - plot of polygons
[7276]696    """
[5897]697
[7841]698    from pylab import ion, hold, plot, savefig, xlabel, \
[7516]699                      ylabel, title, close, title, fill
[5897]700
[7276]701    assert type(polygons_points) == list, \
702                'input must be a list of polygons and/or points'
703
[5897]704    ion()
705    hold(True)
706
[7276]707    if label is None:
708        label = ''
[5897]709
[7516]710    # clamp alpha to sensible range
711    if alpha:
712        try:
713            alpha = float(alpha)
714        except ValueError:
715            alpha = None
716        else:
[7841]717            alpha = max(0.0, min(1.0, alpha))
[7516]718
[7858]719    num_points = len(polygons_points)
[5897]720    colour = []
721    if style is None:
[7276]722        style_type = 'line'
[5897]723        style = []
[7858]724        for i in range(num_points):
[5897]725            style.append(style_type)
726            colour.append('b-')
727    else:
[7858]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)
[7276]737
[5897]738    for i, item in enumerate(polygons_points):
[7858]739        pt_x, pt_y = _poly_xy(item)
740        plot(pt_x, pt_y, colour[i])
[7516]741        if alpha:
[7858]742            fill(pt_x, pt_y, colour[i], alpha=alpha)
[5897]743        xlabel('x')
744        ylabel('y')
745        title(label)
746
747    if figname is not None:
748        savefig(figname)
749    else:
750        savefig('test_image')
751
752    close('all')
753
754
[7858]755def _poly_xy(polygon):
[5897]756    """ this is used within plot_polygons so need to duplicate
757        the first point so can have closed polygon in plot
[7778]758        # @param polygon A set of points defining a polygon.
759        # @param verbose True if this function is to be verbose.
760        # @return A tuple (x, y) of X and Y coordinates of the polygon.
761        # @note We duplicate the first point so can have closed polygon in plot.
[5897]762    """
763
764    try:
[7276]765        polygon = ensure_numeric(polygon, num.float)
[7858]766    except NameError, err:
767        raise NameError, err
[5897]768    except:
[7276]769        msg = ('Polygon %s could not be converted to numeric array'
770               % (str(polygon)))
771        raise Exception, msg
[5897]772
[7858]773    pts_x = num.concatenate((polygon[:, 0], [polygon[0, 0]]), axis = 0)
774    pts_y = num.concatenate((polygon[:, 1], [polygon[0, 1]]), axis = 0)
[7276]775
[7858]776    return pts_x, pts_y
[5897]777
778
[7276]779################################################################################
780# Functions to read and write polygon information
781################################################################################
[5897]782
[7778]783def read_polygon(filename, delimiter=','):
[7858]784    """ Read points assumed to form a polygon.
[7778]785
[7858]786        Also checks to make sure polygon is not complex (self-intersecting).
[7276]787
[7858]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.
[7778]791
[7858]792        There must be exactly two numbers in each line separated by the delimiter.
793        No header.
[5897]794    """
795
796    fid = open(filename)
797    lines = fid.readlines()
798    fid.close()
799    polygon = []
800    for line in lines:
[7276]801        fields = line.split(delimiter)
802        polygon.append([float(fields[0]), float(fields[1])])
[7686]803   
[7690]804    # check this is a valid polygon.
805    if is_complex(polygon, verbose=True):   
[7778]806        msg = 'ERROR: Self-intersecting polygon detected in file '
807        msg += filename +'. A complex polygon will not '
808        msg += 'necessarily break the algorithms within ANUGA, but it'
809        msg += 'usually signifies pathological data. Please fix this file.'
[7690]810        raise Exception, msg
[7686]811   
[5897]812    return polygon
813
[7858]814
[5897]815def write_polygon(polygon, filename=None):
816    """Write polygon to csv file.
[7276]817
818    There will be exactly two numbers, easting and northing, in each line
819    separated by a comma.
820
821    No header.
[5897]822    """
823
824    fid = open(filename, 'w')
825    for point in polygon:
[7276]826        fid.write('%f, %f\n' % point)
[5897]827    fid.close()
828
[7276]829
[5897]830def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
831    """Populate given polygon with uniformly distributed points.
832
833    Input:
834       polygon - list of vertices of polygon
835       number_of_points - (optional) number of points
836       seed - seed for random number generator (default=None)
[7276]837       exclude - list of polygons (inside main polygon) from where points
838                 should be excluded
[5897]839
840    Output:
841       points - list of points inside polygon
842
843    Examples:
844       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
845       will return five randomly selected points inside the unit square
846    """
847
848    from random import uniform, seed as seed_function
849
850    seed_function(seed)
851
852    points = []
853
854    # Find outer extent of polygon
[7858]855    extents = AABB(polygon)
856   
[5897]857    while len(points) < number_of_points:
[7858]858        rand_x = uniform(extents.xmin, extents.xmax)
859        rand_y = uniform(extents.ymin, extents.ymax)
[5897]860
861        append = False
[7858]862        if is_inside_polygon([rand_x, rand_y], polygon):
[5897]863            append = True
864
865            #Check exclusions
866            if exclude is not None:
867                for ex_poly in exclude:
[7858]868                    if is_inside_polygon([rand_x, rand_y], ex_poly):
[5897]869                        append = False
870
871        if append is True:
[7858]872            points.append([rand_x, rand_y])
[5897]873
874    return points
875
[7778]876
[5897]877def point_in_polygon(polygon, delta=1e-8):
878    """Return a point inside a given polygon which will be close to the
879    polygon edge.
880
881    Input:
882       polygon - list of vertices of polygon
883       delta - the square root of 2 * delta is the maximum distance from the
884       polygon points and the returned point.
885    Output:
886       points - a point inside polygon
887
[7276]888       searches in all diagonals and up and down (not left and right).
[5897]889    """
[7276]890
[6534]891    polygon = ensure_numeric(polygon)
892   
[7858]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
[7276]898
[7858]899                    if pt_x == 0:
900                        x_delta = x_mult * delta
901                    else:
902                        x_delta = pt_x + x_mult*pt_x*delta
[5897]903
[7858]904                    if pt_y == 0:
905                        y_delta = y_mult * delta
906                    else:
907                        y_delta = pt_y + y_mult*pt_y*delta
[5897]908
[7858]909                    point = [x_delta, y_delta]
[7276]910
[7858]911                    if is_inside_polygon(point, polygon, closed=False):
912                        return point
913        delta = delta * 0.1
[7276]914
[5897]915
916def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
917    """Calculate the approximate number of triangles inside the
918    bounding polygon and the other interior regions
919
[7276]920    Polygon areas are converted to square Kms
[5897]921
922    FIXME: Add tests for this function
923    """
[7276]924
[5897]925    # TO DO check if any of the regions fall inside one another
926
[7317]927    log.critical('-' * 80)
928    log.critical('Polygon   Max triangle area (m^2)   Total area (km^2) '
929                 'Estimated #triangles')
930    log.critical('-' * 80)
[5897]931       
932    no_triangles = 0.0
933    area = polygon_area(bounding_poly)
[7276]934
[5897]935    for poly, resolution in interior_regions:
936        this_area = polygon_area(poly)
937        this_triangles = this_area/resolution
938        no_triangles += this_triangles
939        area -= this_area
[7276]940
[7317]941        log.critical('Interior %s%s%d'
942                     % (('%.0f' % resolution).ljust(25),
943                        ('%.2f' % (this_area/1000000)).ljust(19), 
944                        this_triangles))
945        #print 'Interior ',
946        #print ('%.0f' % resolution).ljust(25),
947        #print ('%.2f' % (this_area/1000000)).ljust(19),
948        #print '%d' % (this_triangles)
[7276]949
[5897]950    bound_triangles = area/remainder_res
951    no_triangles += bound_triangles
952
[7317]953    log.critical('Bounding %s%s%d'
954                 % (('%.0f' % remainder_res).ljust(25),
955                    ('%.2f' % (area/1000000)).ljust(19),
956                    bound_triangles))
957    #print 'Bounding ',
958    #print ('%.0f' % remainder_res).ljust(25),
959    #print ('%.2f' % (area/1000000)).ljust(19),
960    #print '%d' % (bound_triangles)
[5897]961
962    total_number_of_triangles = no_triangles/0.7
963
[7317]964    log.critical('Estimated total number of triangles: %d'
965                 % total_number_of_triangles)
966    log.critical('Note: This is generally about 20%% '
967                 'less than the final amount')
[5897]968
969    return int(total_number_of_triangles)
970
[7778]971
972def decimate_polygon(polygon, factor=10):
973    """Reduce number of points in polygon by the specified
974    factor (default=10, hence the name of the function) such that
975    the extrema in both axes are preserved.
976
[7276]977##
978# @brief Reduce number of points in polygon by the specified factor.
979# @param polygon The polygon to reduce.
980# @param factor The factor to reduce polygon points by (default 10).
981# @note The extrema of both axes are preserved.
[5897]982
983    Return reduced polygon
984    """
985
986    # FIXME(Ole): This doesn't work at present,
987    # but it isn't critical either
988
989    # Find outer extent of polygon
990    num_polygon = ensure_numeric(polygon)
[7841]991    max_x = max(num_polygon[:, 0])
992    max_y = max(num_polygon[:, 1])
993    min_x = min(num_polygon[:, 0])
994    min_y = min(num_polygon[:, 1])
[5897]995
996    # Keep only some points making sure extrema are kept
[7276]997    reduced_polygon = []
[5897]998    for i, point in enumerate(polygon):
[7858]999        if point[0] in [min_x, max_x] and point[1] in [min_y, max_y]:
[5897]1000            # Keep
1001            reduced_polygon.append(point)
1002        else:
1003            if len(reduced_polygon)*factor < i:
[7276]1004                reduced_polygon.append(point)
[5897]1005
1006    return reduced_polygon
1007
[7778]1008
[6189]1009def interpolate_polyline(data,
1010                         polyline_nodes,
1011                         gauge_neighbour_id,
1012                         interpolation_points=None,
1013                         rtol=1.0e-6,
[7841]1014                         atol=1.0e-8):
[6189]1015    """Interpolate linearly between values data on polyline nodes
[7276]1016    of a polyline to list of interpolation points.
[6189]1017
1018    data is the data on the polyline nodes.
1019
1020    Inputs:
1021      data: Vector or array of data at the polyline nodes.
[7276]1022      polyline_nodes: Location of nodes where data is available.
[6189]1023      gauge_neighbour_id: ?
1024      interpolation_points: Interpolate polyline data to these positions.
1025          List of coordinate pairs [x, y] of
[7276]1026          data points or an nx2 numeric array or a Geospatial_data object
1027      rtol, atol: Used to determine whether a point is on the polyline or not.
1028                  See point_on_line.
[6189]1029
1030    Output:
1031      Interpolated values at interpolation points
1032    """
[7276]1033
[6189]1034    if isinstance(interpolation_points, Geospatial_data):
[7276]1035        interpolation_points = interpolation_points.\
1036                                    get_data_points(absolute=True)
[6189]1037
[7276]1038    interpolated_values = num.zeros(len(interpolation_points), num.float)
[6189]1039
[7276]1040    data = ensure_numeric(data, num.float)
1041    polyline_nodes = ensure_numeric(polyline_nodes, num.float)
1042    interpolation_points = ensure_numeric(interpolation_points, num.float)
1043    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.int)
[6189]1044
[7841]1045    num_nodes = polyline_nodes.shape[0]    # Number of nodes in polyline
[7276]1046
[6189]1047    # Input sanity check
[7841]1048    assert_msg = 'interpolation_points are not given (interpolate.py)'
1049    assert interpolation_points is not None, assert_msg
[7276]1050
[7841]1051    assert_msg = 'function value must be specified at every interpolation node'
1052    assert data.shape[0] == polyline_nodes.shape[0], assert_msg
[7276]1053
[7841]1054    assert_msg = 'Must define function value at one or more nodes'
1055    assert data.shape[0] > 0, assert_msg
[6189]1056
[7841]1057    if num_nodes == 1:
1058        assert_msg = 'Polyline contained only one point. I need more. '
1059        assert_msg += str(data)
1060        raise Exception, assert_msg
1061    elif num_nodes > 1:
[6189]1062        _interpolate_polyline(data,
1063                              polyline_nodes,
1064                              gauge_neighbour_id,
[7276]1065                              interpolation_points,
[6189]1066                              interpolated_values,
1067                              rtol,
1068                              atol)
[7276]1069
[7699]1070
[6189]1071    return interpolated_values
1072
[7699]1073   
1074def polylist2points_verts(polylist):
1075    """ Convert a list of polygons to discrete points and vertices.
1076    """
1077   
1078    offset = 0
1079    points = []
1080    vertices = []
1081    for poly in polylist:
1082        points.extend(poly)
1083        vertices.extend([[i, i+1] for i in range(offset, offset+len(poly)-1)])
1084        offset += len(poly)
1085               
1086    return points, vertices
[7276]1087
1088
1089################################################################################
1090# Initialise module
1091################################################################################
[5897]1092
[6119]1093from anuga.utilities import compile
1094if compile.can_use_C_extension('polygon_ext.c'):
[5897]1095    # Underlying C implementations can be accessed
1096    from polygon_ext import _point_on_line
1097    from polygon_ext import _separate_points_by_polygon
[6189]1098    from polygon_ext import _interpolate_polyline   
[6535]1099    from polygon_ext import _is_inside_triangle       
[5897]1100    #from polygon_ext import _intersection
1101
1102else:
[7858]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)
[5897]1107
1108
1109if __name__ == "__main__":
1110    pass
Note: See TracBrowser for help on using the repository browser.