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

Last change on this file since 8053 was 8053, checked in by habili, 9 years ago

Addition of line intersect code. This will determine which triangles intersect or which triangle contains a line.

File size: 39.2 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, 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    unsorted_segs.append([list(polygon[0]), list(polygon[-1])])
381   
382    # all segments must point in same direction
383    for val in unsorted_segs:
384        if val[0][0] > val[1][0]:
385            val[0], val[1] = val[1], val[0]   
386           
387    l_x = sorted(unsorted_segs, key=key_xpos)
388
389    comparisons = 0
390   
391    # loop through, only comparing lines that partially overlap in x
392    for index, leftmost in enumerate(l_x):
393        cmp = index+1
394        while cmp < len(l_x) and leftmost[1][0] > l_x[cmp][0][0]:
395            if not segments_joined(leftmost, l_x[cmp]):
396                (type, point) = intersection(leftmost, l_x[cmp])
397                comparisons += 1
398                if type != 0 and type != 4 or (type == 2 and list(point[0]) !=\
399                                                list(point[1])):
400                    if verbose:
401                        print 'Self-intersecting polygon found, type ', type
402                        print 'point', point,
403                        print 'vertices: ', leftmost, ' - ', l_x[cmp] 
404                    return True           
405            cmp += 1
406       
407    return False
408   
409
410def is_inside_polygon(point, polygon, closed=True, verbose=False):
411    """Determine if one point is inside a polygon
412
413    See inside_polygon for more details
414    """
415
416    indices = inside_polygon(point, polygon, closed, verbose)
417
418    if indices.shape[0] == 1:
419        return True
420    elif indices.shape[0] == 0:
421        return False
422    else:
423        msg = 'is_inside_polygon must be invoked with one point only'
424        raise Exception(msg)
425
426##
427# @brief Determine which of a set of points are inside a polygon.
428# @param points A set of points (tuple, list or array).
429# @param polygon A set of points defining a polygon (tuple, list or array).
430# @param closed True if points on boundary are considered 'inside' polygon.
431# @param verbose True if this function is to be verbose.
432# @return A list of indices of points inside the polygon.
433def inside_polygon(points, polygon, closed=True, verbose=False):
434    """Determine points inside a polygon
435
436       Functions inside_polygon and outside_polygon have been defined in
437       terms of separate_by_polygon which will put all inside indices in
438       the first part of the indices array and outside indices in the last
439
440       See separate_points_by_polygon for documentation
441
442       points and polygon can be a geospatial instance,
443       a list or a numeric array
444    """
445
446    try:
447        points = ensure_absolute(points)
448    except NameError, err:
449        raise NameError, err
450    except:
451        # If this fails it is going to be because the points can't be
452        # converted to a numeric array.
453        msg = 'Points could not be converted to numeric array' 
454        raise Exception, msg
455
456    polygon = ensure_absolute(polygon)       
457    try:
458        polygon = ensure_absolute(polygon)
459    except NameError, e:
460        raise NameError, e
461    except:
462        # If this fails it is going to be because the points can't be
463        # converted to a numeric array.
464        msg = ('Polygon %s could not be converted to numeric array'
465               % (str(polygon)))
466        raise Exception, msg
467
468    if len(points.shape) == 1:
469        # Only one point was passed in. Convert to array of points
470        points = num.reshape(points, (1,2))
471
472    indices, count = separate_points_by_polygon(points, polygon,
473                                                closed=closed,
474                                                verbose=verbose)
475
476    # Return indices of points inside polygon
477    return indices[:count]
478
479##
480# @brief Determine if one point is outside a polygon.
481# @param point The point of interest.
482# @param polygon The polygon to test inclusion in.
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 True if point is outside the polygon.
486# @note Uses inside_polygon() to do the work.
487def is_outside_polygon(point, polygon, closed=True, verbose=False,
488                       points_geo_ref=None, polygon_geo_ref=None):
489    """Determine if one point is outside a polygon
490
491    See outside_polygon for more details
492    """
493
494    indices = outside_polygon(point, polygon, closed, verbose)
495
496    if indices.shape[0] == 1:
497        return True
498    elif indices.shape[0] == 0:
499        return False
500    else:
501        msg = 'is_outside_polygon must be invoked with one point only'
502        raise Exception, msg
503
504##
505# @brief Determine which of a set of points are outside a polygon.
506# @param points A set of points (tuple, list or array).
507# @param polygon A set of points defining a polygon (tuple, list or array).
508# @param closed True if points on boundary are considered 'inside' polygon.
509# @param verbose True if this function is to be verbose.
510# @return A list of indices of points outside the polygon.
511def outside_polygon(points, polygon, closed = True, verbose = False):
512    """Determine points outside a polygon
513
514       Functions inside_polygon and outside_polygon have been defined in
515       terms of separate_by_polygon which will put all inside indices in
516       the first part of the indices array and outside indices in the last
517
518       See separate_points_by_polygon for documentation
519    """
520
521    try:
522        points = ensure_numeric(points, num.float)
523    except NameError, e:
524        raise NameError, e
525    except:
526        msg = 'Points could not be converted to numeric array'
527        raise Exception, msg
528
529    try:
530        polygon = ensure_numeric(polygon, num.float)
531    except NameError, e:
532        raise NameError, e
533    except:
534        msg = 'Polygon could not be converted to numeric array'
535        raise Exception, msg
536
537    if len(points.shape) == 1:
538        # Only one point was passed in. Convert to array of points
539        points = num.reshape(points, (1, 2))
540
541    indices, count = separate_points_by_polygon(points, polygon,
542                                                closed=closed,
543                                                verbose=verbose)
544
545    # Return indices of points outside polygon
546    if count == len(indices):
547        # No points are outside
548        return num.array([])
549    else:
550        return indices[count:][::-1]  #return reversed
551
552##
553# @brief Separate a list of points into two sets inside+outside a polygon.
554# @param points A set of points (tuple, list or array).
555# @param polygon A set of points defining a polygon (tuple, list or array).
556# @param closed True if points on boundary are considered 'inside' polygon.
557# @param verbose True if this function is to be verbose.
558# @return A tuple (in, out) of point indices for poinst inside amd outside.
559def in_and_outside_polygon(points, polygon, closed=True, verbose=False):
560    """Determine points inside and outside a polygon
561
562       See separate_points_by_polygon for documentation
563
564       Returns an array of points inside and array of points outside the polygon
565    """
566
567    try:
568        points = ensure_numeric(points, num.float)
569    except NameError, e:
570        raise NameError, e
571    except:
572        msg = 'Points could not be converted to numeric array'
573        raise Exception, msg
574
575    try:
576        polygon = ensure_numeric(polygon, num.float)
577    except NameError, e:
578        raise NameError, e
579    except:
580        msg = 'Polygon could not be converted to numeric array'
581        raise Exception, msg
582
583    if len(points.shape) == 1:
584        # Only one point was passed in. Convert to array of points
585        points = num.reshape(points, (1, 2))
586
587    indices, count = separate_points_by_polygon(points, polygon,
588                                                closed=closed,
589                                                verbose=verbose)
590
591    # Returns indices of points inside and indices of points outside
592    # the polygon
593    if count == len(indices):
594        # No points are outside
595        return indices[:count], []
596    else:
597        return  indices[:count], indices[count:][::-1]  #return reversed
598
599
600
601def separate_points_by_polygon(points, polygon,
602                               closed=True, 
603                               check_input=True,
604                               verbose=False):
605    """Determine whether points are inside or outside a polygon
606
607    Input:
608       points - Tuple of (x, y) coordinates, or list of tuples
609       polygon - list of vertices of polygon
610       closed - (optional) determine whether points on boundary should be
611       regarded as belonging to the polygon (closed = True)
612       or not (closed = False)
613       check_input: Allows faster execution if set to False
614
615    Outputs:
616       indices: array of same length as points with indices of points falling
617       inside the polygon listed from the beginning and indices of points
618       falling outside listed from the end.
619
620       count: count of points falling inside the polygon
621
622       The indices of points inside are obtained as indices[:count]
623       The indices of points outside are obtained as indices[count:]
624
625    Examples:
626       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
627
628       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
629       will return the indices [0, 2, 1] and count == 2 as only the first
630       and the last point are inside the unit square
631
632    Remarks:
633       The vertices may be listed clockwise or counterclockwise and
634       the first point may optionally be repeated.
635       Polygons do not need to be convex.
636       Polygons can have holes in them and points inside a hole is
637       regarded as being outside the polygon.
638
639    Algorithm is based on work by Darel Finley,
640    http://www.alienryderflex.com/polygon/
641
642    Uses underlying C-implementation in polygon_ext.c
643    """
644
645    if check_input:
646        #Input checks
647        assert isinstance(closed, bool), \
648                    'Keyword argument "closed" must be boolean'
649        assert isinstance(verbose, bool), \
650                    'Keyword argument "verbose" must be boolean'
651
652        try:
653            points = ensure_numeric(points, num.float)
654        except NameError, e:
655            raise NameError, e
656        except:
657            msg = 'Points could not be converted to numeric array'
658            raise Exception(msg)
659
660        try:
661            polygon = ensure_numeric(polygon, num.float)
662        except NameError, e:
663            raise NameError(e)
664        except:
665            msg = 'Polygon could not be converted to numeric array'
666            raise Exception(msg)
667
668        msg = 'Polygon array must be a 2d array of vertices'
669        assert len(polygon.shape) == 2, msg
670
671        msg = 'Polygon array must have two columns' 
672        assert polygon.shape[1] == 2, msg
673
674
675        msg = ('Points array must be 1 or 2 dimensional. '
676               'I got %d dimensions' % len(points.shape))
677        assert 0 < len(points.shape) < 3, msg
678
679        if len(points.shape) == 1:
680            # Only one point was passed in.  Convert to array of points.
681            points = num.reshape(points, (1, 2))
682   
683            msg = ('Point array must have two columns (x,y), '
684                   'I got points.shape[1]=%d' % points.shape[0])
685            assert points.shape[1]==2, msg
686
687       
688            msg = ('Points array must be a 2d array. I got %s.'
689                   % str(points[:30]))
690            assert len(points.shape) == 2, msg
691
692            msg = 'Points array must have two columns'
693            assert points.shape[1] == 2, msg
694
695    N = polygon.shape[0] # Number of vertices in polygon
696    M = points.shape[0]  # Number of points
697
698    indices = num.zeros(M, num.int)
699
700    count = _separate_points_by_polygon(points, polygon, indices,
701                                        int(closed), int(verbose))
702
703    if verbose:
704        log.critical('Found %d points (out of %d) inside polygon' % (count, M))
705
706    return indices, count
707
708
709def polygon_area(input_polygon):
710    """ Determine area of arbitrary polygon.
711
712        input_polygon The polygon to get area of.
713       
714        return A scalar value for the polygon area.
715
716        Reference:     http://mathworld.wolfram.com/PolygonArea.html
717    """
718    # Move polygon to origin (0,0) to avoid rounding errors
719    # This makes a copy of the polygon to avoid destroying it
720    input_polygon = ensure_numeric(input_polygon)
721    min_x = min(input_polygon[:, 0])
722    min_y = min(input_polygon[:, 1])
723    polygon = input_polygon - [min_x, min_y]
724
725    # Compute area
726    n = len(polygon)
727    poly_area = 0.0
728
729    for i in range(n):
730        pti = polygon[i]
731        if i == n-1:
732            pt1 = polygon[0]
733        else:
734            pt1 = polygon[i+1]
735        xi = pti[0]
736        yi1 = pt1[1]
737        xi1 = pt1[0]
738        yi = pti[1]
739        poly_area += xi*yi1 - xi1*yi
740
741    return abs(poly_area/2)
742
743
744def plot_polygons(polygons_points,
745                  style=None,
746                  figname=None,
747                  label=None,
748                  alpha=None):
749    """ Take list of polygons and plot.
750
751    Inputs:
752
753    polygons         - list of polygons
754
755    style            - style list corresponding to each polygon
756                     - for a polygon, use 'line'
757                     - for points falling outside a polygon, use 'outside'
758                     - style can also be user defined as in normal pylab plot.
759
760    figname          - name to save figure to
761
762    label            - title for plotA
763
764    alpha            - transparency of polygon fill, 0.0=none, 1.0=solid
765                       if not supplied, no fill.
766
767    Outputs:
768
769    - plot of polygons
770    """
771
772    from pylab import ion, hold, plot, savefig, xlabel, \
773                      ylabel, title, close, title, fill
774
775    assert type(polygons_points) == list, \
776                'input must be a list of polygons and/or points'
777
778    ion()
779    hold(True)
780
781    if label is None:
782        label = ''
783
784    # clamp alpha to sensible range
785    if alpha:
786        try:
787            alpha = float(alpha)
788        except ValueError:
789            alpha = None
790        else:
791            alpha = max(0.0, min(1.0, alpha))
792
793    num_points = len(polygons_points)
794    colour = []
795    if style is None:
796        style_type = 'line'
797        style = []
798        for i in range(num_points):
799            style.append(style_type)
800            colour.append('b-')
801    else:
802        for style_name in style:
803            if style_name == 'line':
804                colour.append('b-')
805            if style_name == 'outside':
806                colour.append('r.')
807            if style_name == 'point':
808                colour.append('g.')
809            if style_name not in ['line', 'outside', 'point']:
810                colour.append(style_name)
811
812    for i, item in enumerate(polygons_points):
813        pt_x, pt_y = _poly_xy(item)
814        plot(pt_x, pt_y, colour[i])
815        if alpha:
816            fill(pt_x, pt_y, colour[i], alpha=alpha)
817        xlabel('x')
818        ylabel('y')
819        title(label)
820
821    if figname is not None:
822        savefig(figname)
823    else:
824        savefig('test_image')
825
826    close('all')
827
828
829def _poly_xy(polygon):
830    """ this is used within plot_polygons so need to duplicate
831        the first point so can have closed polygon in plot
832        # @param polygon A set of points defining a polygon.
833        # @param verbose True if this function is to be verbose.
834        # @return A tuple (x, y) of X and Y coordinates of the polygon.
835        # @note We duplicate the first point so can have closed polygon in plot.
836    """
837
838    try:
839        polygon = ensure_numeric(polygon, num.float)
840    except NameError, err:
841        raise NameError, err
842    except:
843        msg = ('Polygon %s could not be converted to numeric array'
844               % (str(polygon)))
845        raise Exception, msg
846
847    pts_x = num.concatenate((polygon[:, 0], [polygon[0, 0]]), axis = 0)
848    pts_y = num.concatenate((polygon[:, 1], [polygon[0, 1]]), axis = 0)
849
850    return pts_x, pts_y
851
852
853################################################################################
854# Functions to read and write polygon information
855################################################################################
856
857def read_polygon(filename, delimiter=','):
858    """ Read points assumed to form a polygon.
859
860        Also checks to make sure polygon is not complex (self-intersecting).
861
862        filename Path to file containing polygon data.
863        delimiter Delimiter to split polygon data with.
864        A list of point data from the polygon file.
865
866        There must be exactly two numbers in each line separated by the delimiter.
867        No header.
868    """
869
870    fid = open(filename)
871    lines = fid.readlines()
872    fid.close()
873    polygon = []
874    for line in lines:
875        fields = line.split(delimiter)
876        polygon.append([float(fields[0]), float(fields[1])])
877   
878    # check this is a valid polygon.
879    if is_complex(polygon, verbose=True):   
880        msg = 'ERROR: Self-intersecting polygon detected in file '
881        msg += filename +'. A complex polygon will not '
882        msg += 'necessarily break the algorithms within ANUGA, but it'
883        msg += 'usually signifies pathological data. Please fix this file.'
884        raise Exception, msg
885   
886    return polygon
887
888
889def write_polygon(polygon, filename=None):
890    """Write polygon to csv file.
891
892    There will be exactly two numbers, easting and northing, in each line
893    separated by a comma.
894
895    No header.
896    """
897
898    fid = open(filename, 'w')
899    for point in polygon:
900        fid.write('%f, %f\n' % point)
901    fid.close()
902
903
904def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
905    """Populate given polygon with uniformly distributed points.
906
907    Input:
908       polygon - list of vertices of polygon
909       number_of_points - (optional) number of points
910       seed - seed for random number generator (default=None)
911       exclude - list of polygons (inside main polygon) from where points
912                 should be excluded
913
914    Output:
915       points - list of points inside polygon
916
917    Examples:
918       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
919       will return five randomly selected points inside the unit square
920    """
921
922    from random import uniform, seed as seed_function
923
924    seed_function(seed)
925
926    points = []
927
928    # Find outer extent of polygon
929    extents = AABB(polygon)
930   
931    while len(points) < number_of_points:
932        rand_x = uniform(extents.xmin, extents.xmax)
933        rand_y = uniform(extents.ymin, extents.ymax)
934
935        append = False
936        if is_inside_polygon([rand_x, rand_y], polygon):
937            append = True
938
939            #Check exclusions
940            if exclude is not None:
941                for ex_poly in exclude:
942                    if is_inside_polygon([rand_x, rand_y], ex_poly):
943                        append = False
944
945        if append is True:
946            points.append([rand_x, rand_y])
947
948    return points
949
950
951def point_in_polygon(polygon, delta=1e-8):
952    """Return a point inside a given polygon which will be close to the
953    polygon edge.
954
955    Input:
956       polygon - list of vertices of polygon
957       delta - the square root of 2 * delta is the maximum distance from the
958       polygon points and the returned point.
959    Output:
960       points - a point inside polygon
961
962       searches in all diagonals and up and down (not left and right).
963    """
964
965    polygon = ensure_numeric(polygon)
966   
967    while True:
968        for poly_point in polygon:
969            for x_mult in range(-1, 2):
970                for y_mult in range(-1, 2):
971                    pt_x, pt_y = poly_point
972
973                    if pt_x == 0:
974                        x_delta = x_mult * delta
975                    else:
976                        x_delta = pt_x + x_mult*pt_x*delta
977
978                    if pt_y == 0:
979                        y_delta = y_mult * delta
980                    else:
981                        y_delta = pt_y + y_mult*pt_y*delta
982
983                    point = [x_delta, y_delta]
984
985                    if is_inside_polygon(point, polygon, closed=False):
986                        return point
987        delta = delta * 0.1
988
989
990def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
991    """Calculate the approximate number of triangles inside the
992    bounding polygon and the other interior regions
993
994    Polygon areas are converted to square Kms
995
996    FIXME: Add tests for this function
997    """
998
999    # TO DO check if any of the regions fall inside one another
1000
1001    log.critical('-' * 80)
1002    log.critical('Polygon   Max triangle area (m^2)   Total area (km^2) '
1003                 'Estimated #triangles')
1004    log.critical('-' * 80)
1005       
1006    no_triangles = 0.0
1007    area = polygon_area(bounding_poly)
1008
1009    for poly, resolution in interior_regions:
1010        this_area = polygon_area(poly)
1011        this_triangles = this_area/resolution
1012        no_triangles += this_triangles
1013        area -= this_area
1014
1015        log.critical('Interior %s%s%d'
1016                     % (('%.0f' % resolution).ljust(25),
1017                        ('%.2f' % (this_area/1000000)).ljust(19), 
1018                        this_triangles))
1019        #print 'Interior ',
1020        #print ('%.0f' % resolution).ljust(25),
1021        #print ('%.2f' % (this_area/1000000)).ljust(19),
1022        #print '%d' % (this_triangles)
1023
1024    bound_triangles = area/remainder_res
1025    no_triangles += bound_triangles
1026
1027    log.critical('Bounding %s%s%d'
1028                 % (('%.0f' % remainder_res).ljust(25),
1029                    ('%.2f' % (area/1000000)).ljust(19),
1030                    bound_triangles))
1031    #print 'Bounding ',
1032    #print ('%.0f' % remainder_res).ljust(25),
1033    #print ('%.2f' % (area/1000000)).ljust(19),
1034    #print '%d' % (bound_triangles)
1035
1036    total_number_of_triangles = no_triangles/0.7
1037
1038    log.critical('Estimated total number of triangles: %d'
1039                 % total_number_of_triangles)
1040    log.critical('Note: This is generally about 20%% '
1041                 'less than the final amount')
1042
1043    return int(total_number_of_triangles)
1044
1045
1046def decimate_polygon(polygon, factor=10):
1047    """Reduce number of points in polygon by the specified
1048    factor (default=10, hence the name of the function) such that
1049    the extrema in both axes are preserved.
1050
1051##
1052# @brief Reduce number of points in polygon by the specified factor.
1053# @param polygon The polygon to reduce.
1054# @param factor The factor to reduce polygon points by (default 10).
1055# @note The extrema of both axes are preserved.
1056
1057    Return reduced polygon
1058    """
1059
1060    # FIXME(Ole): This doesn't work at present,
1061    # but it isn't critical either
1062
1063    # Find outer extent of polygon
1064    num_polygon = ensure_numeric(polygon)
1065    max_x = max(num_polygon[:, 0])
1066    max_y = max(num_polygon[:, 1])
1067    min_x = min(num_polygon[:, 0])
1068    min_y = min(num_polygon[:, 1])
1069
1070    # Keep only some points making sure extrema are kept
1071    reduced_polygon = []
1072    for i, point in enumerate(polygon):
1073        if point[0] in [min_x, max_x] and point[1] in [min_y, max_y]:
1074            # Keep
1075            reduced_polygon.append(point)
1076        else:
1077            if len(reduced_polygon)*factor < i:
1078                reduced_polygon.append(point)
1079
1080    return reduced_polygon
1081
1082
1083def interpolate_polyline(data,
1084                         polyline_nodes,
1085                         gauge_neighbour_id,
1086                         interpolation_points=None,
1087                         rtol=1.0e-6,
1088                         atol=1.0e-8):
1089    """Interpolate linearly between values data on polyline nodes
1090    of a polyline to list of interpolation points.
1091
1092    data is the data on the polyline nodes.
1093
1094    Inputs:
1095      data: Vector or array of data at the polyline nodes.
1096      polyline_nodes: Location of nodes where data is available.
1097      gauge_neighbour_id: ?
1098      interpolation_points: Interpolate polyline data to these positions.
1099          List of coordinate pairs [x, y] of
1100          data points or an nx2 numeric array or a Geospatial_data object
1101      rtol, atol: Used to determine whether a point is on the polyline or not.
1102                  See point_on_line.
1103
1104    Output:
1105      Interpolated values at interpolation points
1106    """
1107
1108    if isinstance(interpolation_points, Geospatial_data):
1109        interpolation_points = interpolation_points.\
1110                                    get_data_points(absolute=True)
1111
1112    interpolated_values = num.zeros(len(interpolation_points), num.float)
1113
1114    data = ensure_numeric(data, num.float)
1115    polyline_nodes = ensure_numeric(polyline_nodes, num.float)
1116    interpolation_points = ensure_numeric(interpolation_points, num.float)
1117    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.int)
1118
1119    num_nodes = polyline_nodes.shape[0]    # Number of nodes in polyline
1120
1121    # Input sanity check
1122    assert_msg = 'interpolation_points are not given (interpolate.py)'
1123    assert interpolation_points is not None, assert_msg
1124
1125    assert_msg = 'function value must be specified at every interpolation node'
1126    assert data.shape[0] == polyline_nodes.shape[0], assert_msg
1127
1128    assert_msg = 'Must define function value at one or more nodes'
1129    assert data.shape[0] > 0, assert_msg
1130
1131    if num_nodes == 1:
1132        assert_msg = 'Polyline contained only one point. I need more. '
1133        assert_msg += str(data)
1134        raise Exception, assert_msg
1135    elif num_nodes > 1:
1136        _interpolate_polyline(data,
1137                              polyline_nodes,
1138                              gauge_neighbour_id,
1139                              interpolation_points,
1140                              interpolated_values,
1141                              rtol,
1142                              atol)
1143
1144
1145    return interpolated_values
1146
1147   
1148def polylist2points_verts(polylist):
1149    """ Convert a list of polygons to discrete points and vertices.
1150    """
1151   
1152    offset = 0
1153    points = []
1154    vertices = []
1155    for poly in polylist:
1156        points.extend(poly)
1157        vertices.extend([[i, i+1] for i in range(offset, offset+len(poly)-1)])
1158        offset += len(poly)
1159               
1160    return points, vertices
1161
1162
1163################################################################################
1164# Initialise module
1165################################################################################
1166
1167from anuga.utilities import compile
1168if compile.can_use_C_extension('polygon_ext.c'):
1169    # Underlying C implementations can be accessed
1170    from polygon_ext import _point_on_line
1171    from polygon_ext import _separate_points_by_polygon
1172    from polygon_ext import _interpolate_polyline   
1173    from polygon_ext import _polygon_overlap
1174    from polygon_ext import _line_intersect
1175    from polygon_ext import _is_inside_triangle       
1176    #from polygon_ext import _intersection
1177
1178else:
1179    ERROR_MSG = 'C implementations could not be accessed by %s.\n ' % __file__
1180    ERROR_MSG += 'Make sure compile_all.py has been run as described in '
1181    ERROR_MSG += 'the ANUGA installation guide.'
1182    raise Exception(ERROR_MSG)
1183
1184
1185if __name__ == "__main__":
1186    pass
Note: See TracBrowser for help on using the repository browser.