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

Last change on this file since 8062 was 8062, checked in by steve, 13 years ago

Fixed up some comments

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