# source:trunk/anuga_core/source/anuga/geometry/test_polygon.py@8038

Last change on this file since 8038 was 8038, checked in by habili, 13 years ago

File size: 72.0 KB
RevLine
[5897]1#!/usr/bin/env python
2
[7866]3""" Test suite to test polygon functionality. """
4
[5897]5import unittest
[7276]6import numpy as num
[5897]7from anuga.utilities.numerical_tools import ensure_numeric
8from anuga.utilities.system_tools import get_pathname_from_package
9
[7866]10from polygon import _poly_xy, separate_points_by_polygon, \
11                    populate_polygon, polygon_area, is_inside_polygon, \
13                    plot_polygons, outside_polygon, is_outside_polygon, \
[8032]14                    intersection, is_complex, polygon_overlap, not_polygon_overlap,\
[8038]15                    polyline_overlap, not_polyline_overlap,\
[8032]16                    is_inside_triangle, interpolate_polyline, inside_polygon, \
17                    in_and_outside_polygon
[7866]18
[7858]19from polygon_function import Polygon_function
[5897]20from anuga.coordinate_transforms.geo_reference import Geo_reference
21from anuga.geospatial_data.geospatial_data import Geospatial_data
22
[7276]23
[5897]24def test_function(x, y):
25    return x+y
26
[7276]27
[5897]28class Test_Polygon(unittest.TestCase):
29    def setUp(self):
30        pass
31
32    def tearDown(self):
33        pass
34
35    def test_that_C_extension_compiles(self):
36        FN = 'polygon_ext.c'
37        try:
38            import polygon_ext
39        except:
40            from compile import compile
41
42            try:
43                compile(FN)
44            except:
[7276]45                raise Exception, 'Could not compile %s' %FN
[5897]46            else:
47                import polygon_ext
48
49    # Polygon stuff
50    def test_polygon_function_constants(self):
51        p1 = [[0,0], [10,0], [10,10], [0,10]]
52        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
53
[7276]54        f = Polygon_function([(p1, 1.0)])
55        z = f([5, 5, 27, 35], [5, 9, 8, -5])    # Two first inside p1
56        assert num.allclose(z, [1, 1, 0, 0])
[5897]57
[7276]58        f = Polygon_function([(p2, 2.0)])
59        z = f([5, 5, 27, 35], [5, 9, 8, -5])    # First and last inside p2
60        assert num.allclose(z, [2, 0, 0, 2])
[5897]61
[7276]62        # Combined
63        f = Polygon_function([(p1, 1.0), (p2, 2.0)])
64        z = f([5, 5, 27, 35], [5, 9, 8, -5])
65        assert num.allclose(z, [2, 1, 0, 2])
[5897]66
67    def test_polygon_function_csvfile(self):
68        from os import sep, getenv
69
70        # Get path where this test is run
71        path = get_pathname_from_package('anuga.utilities')
72
73        # Form absolute filename and read
[7276]74        filename = path + sep + 'mainland_only.csv'
[5897]76
[7276]77        f = Polygon_function([(p1, 10.0)])
78        z = f([430000, 480000], [490000, 7720000]) # 1st outside, 2nd inside
[5897]79
[7276]80        assert num.allclose(z, [0, 10])
81
[5897]82    def test_polygon_function_georef(self):
[7276]83        """Check that georeferencing works"""
[5897]84
85        from anuga.coordinate_transforms.geo_reference import Geo_reference
86
87        geo = Geo_reference(56, 200, 1000)
88
89        # Make points 'absolute'
90        p1 = [[200,1000], [210,1000], [210,1010], [200,1010]]
91        p2 = [[200,1000], [210,1010], [215,1005], [220, 1010], [225,1000],
92              [230,1010], [240,990]]
93
[7276]94        f = Polygon_function([(p1, 1.0)], geo_reference=geo)
95        z = f([5, 5, 27, 35], [5, 9, 8, -5]) # Two first inside p1
[5897]96
[7276]97        assert num.allclose(z, [1, 1, 0, 0])
[5897]98
[7276]99        f = Polygon_function([(p2, 2.0)], geo_reference=geo)
100        z = f([5, 5, 27, 35], [5, 9, 8, -5]) # First and last inside p2
101        assert num.allclose(z, [2, 0, 0, 2])
[5897]102
[7276]103        # Combined
104        f = Polygon_function([(p1, 1.0), (p2, 2.0)], geo_reference=geo)
105        z = f([5, 5, 27, 35], [5, 9, 8, -5])
106        assert num.allclose(z, [2, 1, 0, 2])
[5897]107
[7276]108        # Check that it would fail without geo
109        f = Polygon_function([(p1, 1.0), (p2, 2.0)])
110        z = f([5, 5, 27, 35], [5, 9, 8, -5])
111        assert not num.allclose(z, [2, 1, 0, 2])
[5897]112
[7276]113    def test_polygon_function_callable(self):
114        """Check that values passed into Polygon_function can be callable."""
[5897]115
116        p1 = [[0,0], [10,0], [10,10], [0,10]]
117        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
118
[7276]119        f = Polygon_function([(p1, test_function)])
120        z = f([5, 5, 27, 35], [5, 9, 8, -5]) # Two first inside p1
121        assert num.allclose(z, [10, 14, 0, 0])
[5897]122
[7276]123        # Combined
124        f = Polygon_function([(p1, test_function), (p2, 2.0)])
125        z = f([5, 5, 27, 35], [5, 9, 8, -5])
126        assert num.allclose(z, [2, 14, 0, 2])
[5897]127
[7276]128        # Combined w default
129        f = Polygon_function([(p1, test_function), (p2, 2.0)], default = 3.14)
130        z = f([5, 5, 27, 35], [5, 9, 8, -5])
131        assert num.allclose(z, [2, 14, 3.14, 2])
[5897]132
[7276]133        # Combined w default func
134        f = Polygon_function([(p1, test_function), (p2, 2.0)],
135                             default=test_function)
136        z = f([5, 5, 27, 35], [5, 9, 8, -5])
137        assert num.allclose(z, [2, 14, 35, 2])
[5897]138
139    def test_point_on_line(self):
[7276]140        # Endpoints first
141        assert point_on_line([0,0], [[0,0], [1,0]])
142        assert point_on_line([1,0], [[0,0], [1,0]])
[5897]143
[7276]144        # Endpoints first - non-simple
145        assert point_on_line([-0.1,0.0], [[-0.1,0.0], [1.5,0.6]])
146        assert point_on_line([1.5,0.6], [[-0.1,0.0], [1.5,0.6]])
[5897]147
[7276]148        # Then points on line
149        assert point_on_line([0.5,0], [[0,0], [1,0]])
150        assert point_on_line([0,0.5], [[0,1], [0,0]])
151        assert point_on_line([1,0.5], [[1,1], [1,0]])
152        assert point_on_line([0.5,0.5], [[0,0], [1,1]])
[5897]153
[7276]154        # Then points not on line
155        assert not point_on_line([0.5,0], [[0,1], [1,1]])
156        assert not point_on_line([0,0.5], [[0,0], [1,1]])
[5897]157
[7276]158        # From real example that failed
159        assert not point_on_line([40,50], [[40,20], [40,40]])
[5897]160
[7276]161        # From real example that failed
162        assert not point_on_line([40,19], [[40,20], [40,40]])
[5897]163
164        # Degenerate line
[7276]165        assert point_on_line([40,19], [[40,19], [40,19]])
166        assert not point_on_line([40,19], [[40,40], [40,40]])
[5897]167
168    def test_is_inside_polygon_main(self):
169        # Simplest case: Polygon is the unit square
170        polygon = [[0,0], [1,0], [1,1], [0,1]]
171
[7276]172        assert is_inside_polygon((0.5, 0.5), polygon)
173        assert not is_inside_polygon((0.5, 1.5), polygon)
174        assert not is_inside_polygon((0.5, -0.5), polygon)
175        assert not is_inside_polygon((-0.5, 0.5), polygon)
176        assert not is_inside_polygon((1.5, 0.5), polygon)
[5897]177
[7276]178        # Try point on borders
179        assert is_inside_polygon((1., 0.5), polygon, closed=True)
180        assert is_inside_polygon((0.5, 1.), polygon, closed=True)
181        assert is_inside_polygon((0., 0.5), polygon, closed=True)
182        assert is_inside_polygon((0.5, 0.), polygon, closed=True)
[5897]183
[7276]184        assert not is_inside_polygon((0.5, 1.), polygon, closed=False)
185        assert not is_inside_polygon((0., 0.5), polygon, closed=False)
186        assert not is_inside_polygon((0.5, 0.), polygon, closed=False)
187        assert not is_inside_polygon((1., 0.5), polygon, closed=False)
[5897]188
[7276]189    def test_inside_polygon_main(self):
[7841]190        """test_is_inside_polygon
[6534]191
192        Test fast version of of is_inside_polygon
193        """
194
195        # Simplest case: Polygon is the unit square
196        polygon = [[0,0], [1,0], [1,1], [0,1]]
197
[7841]198        assert is_inside_polygon( (0.5, 0.5), polygon )
199        assert not is_inside_polygon( (0.5, 1.5), polygon )
200        assert not is_inside_polygon( (0.5, -0.5), polygon )
201        assert not is_inside_polygon( (-0.5, 0.5), polygon )
202        assert not is_inside_polygon( (1.5, 0.5), polygon )
[6534]203
[7276]204        # Try point on borders
[7841]205        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
206        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
207        assert is_inside_polygon( (0., 0.5), polygon, closed=True)
208        assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
[6534]209
[7841]210        assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
211        assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
212        assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
213        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
[6534]214
215
[5897]216    def test_inside_polygon_main(self):
217        # Simplest case: Polygon is the unit square
[7276]218        polygon = [[0,0], [1,0], [1,1], [0,1]]
[5897]219
220        # From real example (that failed)
[7276]221        polygon = [[20,20], [40,20], [40,40], [20,40]]
222        points = [[40, 50]]
223        res = inside_polygon(points, polygon)
224        assert len(res) == 0
[5897]225
[7276]226        polygon = [[20,20], [40,20], [40,40], [20,40]]
227        points = [[25, 25], [30, 20], [40, 50], [90, 20], [40, 90]]
228        res = inside_polygon(points, polygon)
229        assert len(res) == 2
230        assert num.allclose(res, [0,1])
[5897]231
[7276]232        # More convoluted and non convex polygon
233        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
234        assert is_inside_polygon( (0.5, 0.5), polygon )
235        assert is_inside_polygon( (1, -0.5), polygon )
236        assert is_inside_polygon( (1.5, 0), polygon )
[5897]237
[7276]238        assert not is_inside_polygon( (0.5, 1.5), polygon )
239        assert not is_inside_polygon( (0.5, -0.5), polygon )
[5897]240
[7841]241        assert is_inside_polygon( (0.5, 0.5), polygon )
242        assert is_inside_polygon( (1, -0.5), polygon )
243        assert is_inside_polygon( (1.5, 0), polygon )
[5897]244
[7841]245        assert not is_inside_polygon( (0.5, 1.5), polygon )
246        assert not is_inside_polygon( (0.5, -0.5), polygon )
[5897]247
[7276]248        # Very convoluted polygon
[5897]249        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
[7276]250        assert is_inside_polygon((5, 5), polygon)
251        assert is_inside_polygon((17, 7), polygon)
252        assert is_inside_polygon((27, 2), polygon)
253        assert is_inside_polygon((35, -5), polygon)
254        assert not is_inside_polygon((15, 7), polygon)
255        assert not is_inside_polygon((24, 3), polygon)
256        assert not is_inside_polygon((25, -10), polygon)
[5897]257
[7276]258        # Another combination (that failed)
[5897]259        polygon = [[0,0], [10,0], [10,10], [0,10]]
[7276]260        assert is_inside_polygon((5, 5), polygon)
261        assert is_inside_polygon((7, 7), polygon)
262        assert not is_inside_polygon((-17, 7), polygon)
263        assert not is_inside_polygon((7, 17), polygon)
264        assert not is_inside_polygon((17, 7), polygon)
265        assert not is_inside_polygon((27, 8), polygon)
266        assert not is_inside_polygon((35, -5), polygon)
[5897]267
[7276]268        # Multiple polygons
269        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0], [10,10],
270                   [11,10], [11,11], [10,11], [10,10]]
271        assert is_inside_polygon((0.5, 0.5), polygon)
272        assert is_inside_polygon((10.5, 10.5), polygon)
[5897]273
[7276]274        # FIXME: Fails if point is 5.5, 5.5
275        assert not is_inside_polygon((0, 5.5), polygon)
[5897]276
[7276]277        # Polygon with a hole
[5897]278        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
[7276]279                   [0,0], [1,0], [1,1], [0,1], [0,0]]
[5897]280
[7276]281        assert is_inside_polygon((0, -0.5), polygon)
282        assert not is_inside_polygon((0.5, 0.5), polygon)
[5897]283
284    def test_duplicate_points_being_ok(self):
285        # Simplest case: Polygon is the unit square
286        polygon = [[0,0], [1,0], [1,0], [1,0], [1,1], [0,1], [0,0]]
287
[7276]288        assert is_inside_polygon((0.5, 0.5), polygon)
289        assert not is_inside_polygon((0.5, 1.5), polygon)
290        assert not is_inside_polygon((0.5, -0.5), polygon)
291        assert not is_inside_polygon((-0.5, 0.5), polygon)
292        assert not is_inside_polygon((1.5, 0.5), polygon)
[5897]293
[7276]294        # Try point on borders
295        assert is_inside_polygon((1., 0.5), polygon, closed=True)
296        assert is_inside_polygon((0.5, 1), polygon, closed=True)
297        assert is_inside_polygon((0., 0.5), polygon, closed=True)
298        assert is_inside_polygon((0.5, 0.), polygon, closed=True)
[5897]299
[7276]300        assert not is_inside_polygon((0.5, 1), polygon, closed=False)
301        assert not is_inside_polygon((0., 0.5), polygon, closed=False)
302        assert not is_inside_polygon((0.5, 0.), polygon, closed=False)
303        assert not is_inside_polygon((1., 0.5), polygon, closed=False)
[5897]304
305        # From real example
[7276]306        polygon = [[20,20], [40,20], [40,40], [20,40]]
307        points = [[40, 50]]
308        res = inside_polygon(points, polygon)
309        assert len(res) == 0
[5897]310
311    def test_inside_polygon_vector_version(self):
[7276]312        # Now try the vector formulation returning indices
[5897]313        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[7276]314        points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
315        res = inside_polygon( points, polygon, verbose=False )
[7308]316        assert num.allclose(res, [0,1,2])
[5897]317
318    def test_outside_polygon(self):
[7276]319        # unit square
320        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]321
322        # evaluate to False as the point 0.5, 0.5 is inside the unit square
[7276]323        assert not is_outside_polygon([0.5, 0.5], U)
324
[5897]325        # evaluate to True as the point 1.5, 0.5 is outside the unit square
[7276]326        assert is_outside_polygon([1.5, 0.5], U)
327
328        indices = outside_polygon([[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
329        assert num.allclose(indices, [1])
330
[5897]331        # One more test of vector formulation returning indices
332        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[7276]333        points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
334        res = outside_polygon(points, polygon)
335        assert num.allclose(res, [3, 4])
[5897]336
[7276]337        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
338        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
339                  [0.5, 1.5], [0.5, -0.5]]
340        res = outside_polygon(points, polygon)
[5897]341
[7276]342        assert num.allclose(res, [0, 4, 5])
[5897]343
[7276]344    def test_outside_polygon2(self):
345        # unit square
346        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]347
[7276]348        # evaluate to False as the point 0.5, 1.0 is inside the unit square
349        assert not outside_polygon([0.5, 1.0], U, closed = True)
[5897]350
351        # evaluate to True as the point 0.5, 1.0 is outside the unit square
[7276]352        assert is_outside_polygon([0.5, 1.0], U, closed = False)
[5897]353
354    def test_all_outside_polygon(self):
[7276]355        """Test case where all points are outside poly"""
[5897]356
[7276]357        # unit square
358        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]359
[7276]360        points = [[2,2], [1,3], [-1,1], [0,2]]      # All outside
[5897]361
362        indices, count = separate_points_by_polygon(points, U)
[7276]363        assert count == 0                           # None inside
364        assert num.allclose(indices, [3, 2, 1, 0])
[5897]365
366        indices = outside_polygon(points, U, closed = True)
[7276]367        assert num.allclose(indices, [0, 1, 2, 3])
[5897]368
369        indices = inside_polygon(points, U, closed = True)
[7276]370        assert num.allclose(indices, [])
[5897]371
372    def test_all_inside_polygon(self):
[7276]373        """Test case where all points are inside poly"""
[5897]374
[7276]375        # unit square
376        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]377
[7276]378        points = [[0.5,0.5], [0.2,0.3], [0,0.5]]    # All inside (or on edge)
[5897]379
380        indices, count = separate_points_by_polygon(points, U)
[7276]381        assert count == 3       # All inside
382        assert num.allclose(indices, [0, 1, 2])
[5897]383
[7276]384        indices = outside_polygon(points, U, closed=True)
[6158]385        assert num.allclose(indices, [])
[5897]386
[7276]387        indices = inside_polygon(points, U, closed=True)
388        assert num.allclose(indices, [0, 1, 2])
[5897]389
[7276]390
[5897]391    def test_separate_points_by_polygon(self):
[7276]392        # unit square
393        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]394
[7276]395        indices, count = separate_points_by_polygon([[0.5, 0.5],
396                                                     [1, -0.5],
397                                                     [0.3, 0.2]],
398                                                    U)
399        assert num.allclose( indices, [0, 2, 1] )
[5897]400        assert count == 2
[7276]401
402        # One more test of vector formulation returning indices
[5897]403        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[7276]404        points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
405        res, count = separate_points_by_polygon(points, polygon)
[5897]406
[7276]407        assert num.allclose(res, [0, 1, 2, 4, 3])
408        assert count == 3
[5897]409
410        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[7866]411        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], \
412                        [0.5, -0.5]]
[7276]413        res, count = separate_points_by_polygon( points, polygon )
[5897]414
[7276]415        assert num.allclose( res, [1,2,3,5,4,0] )
416        assert count == 3
[5897]417
[7276]418
[8032]419    def test_polygon_overlap(self):
420        #rectangular polygon, 2 triangles overlap
421        polygon = [[3.,3.],[3.,4.],[5.,4.],[5.,3.]]
422        triangles = [[7.,10.],#does not overlap
423                     [8.,12.],
424                     [9.,10.],
425                     [3., 2.],#intersect
426                     [4., 5.],
427                     [5., 3.],
428                     [3., 2.],#intersect
429                     [4., 5.],
430                     [5., 2.],
431                     [0., 1.],#polygon inside triangle
432                     [5., 50.],
433                     [10., 1.],
434                     [3.5, 3.25],#triangle inside polygon
435                     [4., 3.75],
436                     [4.5, 3.25]]
437
438        indices = polygon_overlap(triangles, polygon)
439        assert num.allclose([1, 2, 3, 4], indices)
440
441        #5 sided polygon, 2 triangles overlap
442        polygon = [[3.,3.],[3.,4.],[5.,4.],[5.5, 3.5],[5.,3.]]
443        triangles = [[7.,10.],#does not overlap
444                     [8.,12.],
445                     [9.,10.],
446                     [3., 2.],#intersect
447                     [4., 5.],
448                     [5., 3.],
449                     [3., 2.],#intersect
450                     [4., 5.],
451                     [5., 2.],
452                     [0., 1.],#polygon inside triangle
453                     [5., 50.],
454                     [10., 1.],
455                     [3.5, 3.25],#triangle inside polygon
456                     [4., 3.75],
457                     [4.5, 3.25]]
458
459        indices = polygon_overlap(triangles, polygon)
460        assert num.allclose([1, 2, 3, 4], indices)
461
462    def test_not_polygon_overlap(self):
463        #rectangular polygon, 2 triangles overlap
464        polygon = [[3.,3.],[3.,4.],[5.,4.],[5.,3.]]
465        triangles = [[7.,10.],#does not overlap
466                     [8.,12.],
467                     [9.,10.],
468                     [3., 2.],#intersect
469                     [4., 5.],
470                     [5., 3.],
471                     [3., 2.],#intersect
472                     [4., 5.],
473                     [5., 2.],
474                     [0., 1.],#polygon inside triangle
475                     [5., 50.],
476                     [10., 1.],
477                     [3.5, 3.25],#triangle inside polygon
478                     [4., 3.75],
479                     [4.5, 3.25]]
480
481        indices = not_polygon_overlap(triangles, polygon)
482        assert num.allclose([0], indices)
483
484        #5 sided polygon, 2 triangles overlap
485        polygon = [[3.,3.],[3.,4.],[5.,4.],[5.5, 3.5],[5.,3.]]
486        triangles = [[7.,10.],#does not overlap
487                     [8.,12.],
488                     [9.,10.],
489                     [3., 2.],#intersect
490                     [4., 5.],
491                     [5., 3.],
492                     [3., 2.],#intersect
493                     [4., 5.],
494                     [5., 2.],
495                     [0., 1.],#polygon inside triangle
496                     [5., 50.],
497                     [10., 1.],
498                     [3.5, 3.25],#triangle inside polygon
499                     [4., 3.75],
500                     [4.5, 3.25]]
501
502        indices = not_polygon_overlap(triangles, polygon)
503        assert num.allclose([0], indices)
504
[8038]505    def test_polyline_overlap(self):
506        #rectangular polygon, 2 triangles overlap
507        polyline = [[3.,3.],[4.,3.]]
508        triangles = [[7.,10.],#does not overlap
509                     [8.,12.],
510                     [9.,10.],
511                     [3., 2.],#intersect
512                     [4., 5.],
513                     [5., 3.],
514                     [3., 2.],#intersect
515                     [4., 5.],
516                     [5., 2.],
517                     [0., 1.],#polygon inside triangle
518                     [5., 50.],
519                     [10., 1.],
520                     [3.5, 3.25],#no overlap
521                     [4., 3.75],
522                     [4.5, 3.25],
523                     [2., 3.],#overlap
524                     [5., 3.],
525                     [3., 1.]]
526
527        indices = polyline_overlap(triangles, polyline)
528        assert num.allclose([1, 2, 3, 5], indices)
529
530    def test_not_polyline_overlap(self):
531        #rectangular polygon, 2 triangles overlap
532        polyline = [[3.,3.],[4.,3.]]
533        triangles = [[7.,10.],#does not overlap
534                     [8.,12.],
535                     [9.,10.],
536                     [3., 2.],#intersect
537                     [4., 5.],
538                     [5., 3.],
539                     [3., 2.],#intersect
540                     [4., 5.],
541                     [5., 2.],
542                     [0., 1.],#polygon inside triangle
543                     [5., 50.],
544                     [10., 1.],
545                     [3.5, 3.25],#no overlap
546                     [4., 3.75],
547                     [4.5, 3.25],
548                     [2., 3.],#overlap
549                     [5., 3.],
550                     [3., 1.]]
551
552        indices = not_polyline_overlap(triangles, polyline)
553        assert num.allclose([4, 0], indices)
554
[6534]555    def test_is_inside_triangle(self):
556        # Simplest case:
[7276]557        triangle = [[0, 0], [1, 0], [0.5, 1]]
[6534]558
[7276]559        assert is_inside_triangle((0.5, 0.5), triangle)
560        assert is_inside_triangle((0.9, 0.1), triangle)
561        assert not is_inside_triangle((0.5, 1.5), triangle)
562        assert not is_inside_triangle((0.5, -0.5), triangle)
563        assert not is_inside_triangle((-0.5, 0.5), triangle)
564        assert not is_inside_triangle((1.5, 0.5), triangle)
[6534]565
[7276]566        # Try point on borders
567        assert is_inside_triangle((0.5, 0), triangle, closed=True)
568        assert is_inside_triangle((1, 0), triangle, closed=True)
[6534]569
[7276]570        assert not is_inside_triangle((0.5, 0), triangle, closed=False)
571        assert not is_inside_triangle((1, 0), triangle, closed=False)
[6534]572
573        # Try vertices
574        for P in triangle:
[7276]575            assert is_inside_triangle(P, triangle, closed=True)
576            assert not is_inside_triangle(P, triangle, closed=False)
577
[6534]578        # Slightly different
[7276]579        triangle = [[0, 0.1], [1, -0.2], [0.5, 1]]
580        assert is_inside_triangle((0.5, 0.5), triangle)
581        assert is_inside_triangle((0.4, 0.1), triangle)
582        assert not is_inside_triangle((1, 1), triangle)
583
[6534]584        # Try vertices
585        for P in triangle:
[7276]586            assert is_inside_triangle(P, triangle, closed=True)
587            assert not is_inside_triangle(P, triangle, closed=False)
[6534]588
[7276]589
[5897]590    def test_populate_polygon(self):
591        polygon = [[0,0], [1,0], [1,1], [0,1]]
592        points = populate_polygon(polygon, 5)
593
594        assert len(points) == 5
595        for point in points:
596            assert is_inside_polygon(point, polygon)
597
[7841]598            assert is_inside_polygon(point, polygon)
[5897]599
[7276]600
601        # Very convoluted polygon
[5897]602        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
603        points = populate_polygon(polygon, 5)
604        assert len(points) == 5
605        for point in points:
606            assert is_inside_polygon(point, polygon)
[7841]607            assert is_inside_polygon(point, polygon)
[5897]608
609
610    def test_populate_polygon_with_exclude(self):
611        polygon = [[0,0], [1,0], [1,1], [0,1]]
[7276]612        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]]     # SW quarter
613        points = populate_polygon(polygon, 5, exclude=[ex_poly])
[5897]614
615        assert len(points) == 5
616        for point in points:
617            assert is_inside_polygon(point, polygon)
[7276]618            assert not is_inside_polygon(point, ex_poly)
[5897]619
[7276]620        # overlap
[5897]621        polygon = [[0,0], [1,0], [1,1], [0,1]]
622        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
[7276]623        points = populate_polygon(polygon, 5, exclude=[ex_poly])
[5897]624
625        assert len(points) == 5
626        for point in points:
627            assert is_inside_polygon(point, polygon)
[7276]628            assert not is_inside_polygon(point, ex_poly)
629
630        # Multiple
[5897]631        polygon = [[0,0], [1,0], [1,1], [0,1]]
[7276]632        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]]    # SW quarter
633        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]]    # NE quarter
[5897]634
[7276]635        points = populate_polygon(polygon, 20, exclude=[ex_poly1, ex_poly2])
636
[5897]637        assert len(points) == 20
638        for point in points:
639            assert is_inside_polygon(point, polygon)
640            assert not is_inside_polygon(point, ex_poly1)
[7276]641            assert not is_inside_polygon(point, ex_poly2)
[5897]642
[7276]643        # Very convoluted polygon
[5897]644        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
645        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
[7276]646        points = populate_polygon(polygon, 20, exclude=[ex_poly])
647
[5897]648        assert len(points) == 20
649        for point in points:
650            assert is_inside_polygon(point, polygon)
[7276]651            assert not is_inside_polygon(point, ex_poly), '%s' % str(point)
[5897]652
653    def test_populate_polygon_with_exclude2(self):
[7276]654        min_outer = 0
[5897]655        max_outer = 1000
[7276]656        polygon_outer = [[min_outer,min_outer], [max_outer,min_outer],
657                         [max_outer,max_outer], [min_outer,max_outer]]
[5897]658
659        delta = 10
660        min_inner1 = min_outer + delta
661        max_inner1 = max_outer - delta
[7276]662        inner1_polygon = [[min_inner1,min_inner1], [max_inner1,min_inner1],
663                          [max_inner1,max_inner1], [min_inner1,max_inner1]]
[5897]664
[7276]665        density_inner2 = 1000
666        min_inner2 = min_outer + 2*delta
667        max_inner2 = max_outer - 2*delta
668        inner2_polygon = [[min_inner2,min_inner2], [max_inner2,min_inner2],
669                          [max_inner2,max_inner2], [min_inner2,max_inner2]]
670
671        points = populate_polygon(polygon_outer, 20,
672                                  exclude=[inner1_polygon, inner2_polygon])
673
[5897]674        assert len(points) == 20
675        for point in points:
676            assert is_inside_polygon(point, polygon_outer)
677            assert not is_inside_polygon(point, inner1_polygon)
[7276]678            assert not is_inside_polygon(point, inner2_polygon)
[5897]679
[7276]680        # Very convoluted polygon
[5897]681        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
682        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
[7276]683        points = populate_polygon(polygon, 20, exclude=[ex_poly])
684
[5897]685        assert len(points) == 20
686        for point in points:
687            assert is_inside_polygon(point, polygon)
[7276]688            assert not is_inside_polygon(point, ex_poly), '%s' % str(point)
[5897]689
690    def test_point_in_polygon(self):
691        polygon = [[0,0], [1,0], [1,1], [0,1]]
692        point = point_in_polygon(polygon)
693        assert is_inside_polygon(point, polygon)
694
[7276]695        # this may get into a vicious loop
696        # polygon = [[1e32,1e54], [1,0], [1,1], [0,1]]
697
[5897]698        polygon = [[1e15,1e7], [1,0], [1,1], [0,1]]
699        point = point_in_polygon(polygon)
700        assert is_inside_polygon(point, polygon)
701
702        polygon = [[0,0], [1,0], [1,1], [1e8,1e8]]
703        point = point_in_polygon(polygon)
704        assert is_inside_polygon(point, polygon)
705
706        polygon = [[1e32,1e54], [-1e32,1e54], [1e32,-1e54]]
707        point = point_in_polygon(polygon)
708        assert is_inside_polygon(point, polygon)
709
710        polygon = [[1e18,1e15], [1,0], [0,1]]
711        point = point_in_polygon(polygon)
712        assert is_inside_polygon(point, polygon)
713
714    def test_in_and_outside_polygon_main(self):
[7276]715        # Simplest case: Polygon is the unit square
716        polygon = [[0,0], [1,0], [1,1], [0,1]]
[5897]717
[7276]718        inside, outside = in_and_outside_polygon((0.5, 0.5), polygon)
719        assert inside[0] == 0
720        assert len(outside) == 0
[5897]721
[7276]722        inside, outside = in_and_outside_polygon((1., 0.5), polygon,
723                                                 closed=True)
724        assert inside[0] == 0
725        assert len(outside) == 0
[5897]726
[7276]727        inside, outside = in_and_outside_polygon((1., 0.5), polygon,
728                                                 closed=False)
729        assert len(inside) == 0
730        assert outside[0] == 0
[5897]731
[7276]732        points = [(1., 0.25),(1., 0.75)]
733        inside, outside = in_and_outside_polygon(points, polygon, closed=True)
[7308]734        assert num.alltrue(inside == [0,1])
[7276]735        assert len(outside) == 0
[5897]736
[7276]737        inside, outside = in_and_outside_polygon(points, polygon, closed=False)
738        assert len(inside) == 0
[7308]739        assert num.alltrue(outside == [0,1])
[5897]740
[7276]741        points = [(100., 0.25), (0.5, 0.5) ]
742        inside, outside = in_and_outside_polygon(points, polygon)
[7308]743        assert num.alltrue(inside == [1])
[7276]744        assert outside[0] == 0
[5897]745
[7276]746        points = [(100., 0.25),(0.5, 0.5), (39,20), (0.6,0.7),(56,43),(67,90)]
747        inside, outside = in_and_outside_polygon(points, polygon)
[7308]748        assert num.alltrue(inside == [1, 3])
749        assert num.alltrue(outside == [0, 2, 4, 5])
[7276]750
[5897]751    def test_intersection1(self):
752        line0 = [[-1,0], [1,0]]
753        line1 = [[0,-1], [0,1]]
754
755        status, value = intersection(line0, line1)
756        assert status == 1
[6158]757        assert num.allclose(value, [0.0, 0.0])
[5897]758
759    def test_intersection2(self):
760        line0 = [[0,0], [24,12]]
761        line1 = [[0,12], [24,0]]
762
763        status, value = intersection(line0, line1)
764        assert status == 1
[6158]765        assert num.allclose(value, [12.0, 6.0])
[5897]766
767        # Swap direction of one line
768        line1 = [[24,0], [0,12]]
769
770        status, value = intersection(line0, line1)
771        assert status == 1
[6158]772        assert num.allclose(value, [12.0, 6.0])
[5897]773
774        # Swap order of lines
775        status, value = intersection(line1, line0)
776        assert status == 1
[7276]777        assert num.allclose(value, [12.0, 6.0])
778
[5897]779    def test_intersection3(self):
780        line0 = [[0,0], [24,12]]
781        line1 = [[0,17], [24,0]]
782
783        status, value = intersection(line0, line1)
784        assert status == 1
[6158]785        assert num.allclose(value, [14.068965517, 7.0344827586])
[5897]786
787        # Swap direction of one line
788        line1 = [[24,0], [0,17]]
789
790        status, value = intersection(line0, line1)
791        assert status == 1
[7276]792        assert num.allclose(value, [14.068965517, 7.0344827586])
[5897]793
794        # Swap order of lines
795        status, value = intersection(line1, line0)
[7276]796        assert status == 1
797        assert num.allclose(value, [14.068965517, 7.0344827586])
[5897]798
799    def test_intersection_endpoints(self):
800        """test_intersection_endpoints(self):
801
802        Test that coinciding endpoints are picked up
803        """
[7276]804
[5897]805        line0 = [[0,0], [1,1]]
806        line1 = [[1,1], [2,1]]
807
808        status, value = intersection(line0, line1)
809        assert status == 1
[6158]810        assert num.allclose(value, [1.0, 1.0])
[5897]811
812        line0 = [[1,1], [2,0]]
813        line1 = [[1,1], [2,1]]
814
815        status, value = intersection(line0, line1)
816        assert status == 1
[7276]817        assert num.allclose(value, [1.0, 1.0])
[5897]818
[5942]819    # This function is a helper function for
820    # the test_intersection_bug_20081110_?() set of tests.
821    # This function tests all parallel line cases for 4 collinear points.
[5933]822    # This function should never be run directly by the unittest code.
823    def helper_test_parallel_intersection_code(self, P1, P2, P3, P4):
824        # lines in same direction, no overlap
825        # 0:         ---->----
826        # 1:                     --------->-----------
[7276]827        line0 = [P1, P2]
828        line1 = [P3, P4]
[5933]829        status, value = intersection(line0, line1)
830        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
831                               (str(status), str(value)))
832        self.failUnless(value is None, 'Expected value of None, got %s' %
833                                       str(value))
834
835        # lines in same direction, no overlap
836        # 0:         ----<----
837        # 1:                     ---------<-----------
[7276]838        line0 = [P2, P1]
839        line1 = [P4, P3]
[5933]840        status, value = intersection(line0, line1)
841        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
842                               (str(status), str(value)))
843        self.failUnless(value is None, 'Expected value of None, got %s' %
844                                       str(value))
845
846        # lines in opposite direction, no overlap
847        # 0:         ----<----
848        # 1:                     --------->-----------
[7276]849        line0 = [P2, P1]
850        line1 = [P3, P4]
[5933]851        status, value = intersection(line0, line1)
852        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
853                               (str(status), str(value)))
854        self.failUnless(value is None, 'Expected value of None, got %s' %
855                                       str(value))
856
857        # lines in opposite direction, no overlap
858        # 0:         ---->----
859        # 1:                     ---------<-----------
[7276]860        line0 = [P1, P2]
861        line1 = [P4, P3]
[5933]862        status, value = intersection(line0, line1)
863        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
864                               (str(status), str(value)))
865        self.failUnless(value is None, 'Expected value of None, got %s' %
866                                       str(value))
867
[7276]868        # ----------------------------------------------------------------------
[5942]869
[5933]870        # line0 fully within line1, same direction
871        # 0:         ---->----
872        # 1:    --------->-----------
873        # value should be line0:
874        #            ---->----
[7276]875        line0 = [P2, P3]
876        line1 = [P1, P4]
[5933]877        status, value = intersection(line0, line1)
878        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
879                               (str(status), str(value)))
[6158]880        self.failUnless(num.allclose(value, line0))
[5933]881
882        # line0 fully within line1, same direction
883        # 0:         ----<----
884        # 1:    ---------<-----------
885        # value should be line0:
886        #            ----<----
[7276]887        line0 = [P3, P2]
888        line1 = [P4, P1]
[5933]889        status, value = intersection(line0, line1)
890        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
891                               (str(status), str(value)))
[6158]892        self.failUnless(num.allclose(value, line0))
[5933]893
894        # line0 fully within line1, opposite direction
895        # 0:         ----<----
896        # 1:    --------->-----------
897        # value should be line0:
898        #            ----<----
[7276]899        line0 = [P3, P2]
900        line1 = [P1, P4]
[5933]901        status, value = intersection(line0, line1)
902        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
903                               (str(status), str(value)))
[6158]904        self.failUnless(num.allclose(value, line0))
[5933]905
906        # line0 fully within line1, opposite direction
907        # 0:         ---->----
908        # 1:    ---------<-----------
909        # value should be line0:
910        #            ---->----
[7276]911        line0 = [P2, P3]
912        line1 = [P4, P1]
[5933]913        status, value = intersection(line0, line1)
914        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
915                               (str(status), str(value)))
[6158]916        self.failUnless(num.allclose(value, line0))
[5933]917
[7276]918        # ----------------------------------------------------------------------
[5942]919
[5933]920        # line1 fully within line0, same direction
921        # 0:    --------->-----------
922        # 1:         ---->----
923        # value should be line1:
924        #            ---->----
[7276]925        line0 = [P1, P4]
926        line1 = [P2, P3]
[5933]927        status, value = intersection(line0, line1)
928        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
929                               (str(status), str(value)))
[6158]930        self.failUnless(num.allclose(value, line1))
[5933]931
932        # line1 fully within line0, same direction
933        # 0:    ---------<-----------
934        # 1:         ----<----
935        # value should be line1:
936        #            ----<----
[7276]937        line0 = [P4, P1]
938        line1 = [P3, P2]
[5933]939        status, value = intersection(line0, line1)
940        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
941                               (str(status), str(value)))
[6158]942        self.failUnless(num.allclose(value, line1))
[5933]943
944        # line1 fully within line0, opposite direction
945        # 0:    ---------<-----------
946        # 1:         ---->----
947        # value should be line1:
948        #            ---->----
[7276]949        line0 = [P4, P1]
950        line1 = [P2, P3]
[5933]951        status, value = intersection(line0, line1)
952        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
953                               (str(status), str(value)))
[6158]954        self.failUnless(num.allclose(value, line1))
[5933]955
956        # line1 fully within line0, opposite direction
957        # 0:    --------->-----------
958        # 1:         ----<----
959        # value should be line1:
960        #            ----<----
[7276]961        line0 = [P1, P4]
962        line1 = [P3, P2]
[5933]963        status, value = intersection(line0, line1)
964        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
965                               (str(status), str(value)))
[6158]966        self.failUnless(num.allclose(value, line1))
[5933]967
[7276]968        # ----------------------------------------------------------------------
[5942]969
[5933]970        # line in same direction, partial overlap
971        # 0:    ----->-----
972        # 1:       ------->--------
973        # value should be segment line1_start->line0_end:
974        #          --->----
[7276]975        line0 = [P1, P3]
976        line1 = [P2, P4]
[5933]977        status, value = intersection(line0, line1)
978        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
979                               (str(status), str(value)))
[6158]980        self.failUnless(num.allclose(value, [line1[0],line0[1]]))
[5933]981
982        # line in same direction, partial overlap
983        # 0:    -----<-----
984        # 1:       -------<--------
985        # value should be segment line0_start->line1_end:
986        #          ---<----
[7276]987        line0 = [P3, P1]
988        line1 = [P4, P2]
[5933]989        status, value = intersection(line0, line1)
990        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
991                               (str(status), str(value)))
[6158]992        self.failUnless(num.allclose(value, [line0[0],line1[1]]))
[5933]993
994        # line in opposite direction, partial overlap
995        # 0:    -----<-----
996        # 1:       ------->--------
997        # value should be segment line0_start->line1_start:
998        #          ---<----
[7276]999        line0 = [P3, P1]
1000        line1 = [P2, P4]
[5933]1001        status, value = intersection(line0, line1)
1002        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1003                               (str(status), str(value)))
[6158]1004        self.failUnless(num.allclose(value, [line0[0],line1[0]]))
[5933]1005
1006        # line in opposite direction, partial overlap
1007        # 0:    ----->-----
1008        # 1:       -------<--------
1009        # value should be segment line1_end->line0_end:
1010        #          --->----
[7276]1011        line0 = [P1, P3]
1012        line1 = [P4, P2]
[5933]1013        status, value = intersection(line0, line1)
1014        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1015                               (str(status), str(value)))
[6158]1016        self.failUnless(num.allclose(value, [line1[1],line0[1]]))
[5933]1017
[7276]1018        # ----------------------------------------------------------------------
[5942]1019
[5933]1020        # line in same direction, partial overlap
1021        # 0:       ------>------
1022        # 1:    ------>------
1023        # value should be segment line0_start->line1_end:
1024        #          --->----
[7276]1025        line0 = [P2, P4]
1026        line1 = [P1, P3]
[5933]1027        status, value = intersection(line0, line1)
1028        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1029                               (str(status), str(value)))
[6158]1030        self.failUnless(num.allclose(value, [line0[0],line1[1]]))
[5933]1031
1032        # line in same direction, partial overlap
1033        # 0:       ------<------
1034        # 1:    ------<------
1035        # value should be segment line1_start->line0_end:
1036        #          ----<-----
[7276]1037        line0 = [P4, P2]
1038        line1 = [P3, P1]
[5933]1039        status, value = intersection(line0, line1)
1040        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1041                               (str(status), str(value)))
[6158]1042        self.failUnless(num.allclose(value, [line1[0],line0[1]]))
[5933]1043
1044        # line in opposite direction, partial overlap
1045        # 0:       ------<------
1046        # 1:    ----->------
1047        # value should be segment line1_end->line0_end:
1048        #          --->----
[7276]1049        line0 = [P4, P2]
1050        line1 = [P1, P3]
[5933]1051        status, value = intersection(line0, line1)
1052        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1053                               (str(status), str(value)))
[6158]1054        self.failUnless(num.allclose(value, [line1[1],line0[1]]))
[5933]1055
1056        # line in opposite direction, partial overlap
1057        # 0:       ------>------
1058        # 1:    -----<------
1059        # value should be segment line0_start->line1_start:
1060        #          ---<----
[7276]1061        line0 = [P2, P4]
1062        line1 = [P3, P1]
[5933]1063        status, value = intersection(line0, line1)
1064        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1065                               (str(status), str(value)))
[6158]1066        self.failUnless(num.allclose(value, [line0[0],line1[0]]))
[5933]1067
[7276]1068        # ----------------------------------------------------------------------
[5942]1069
1070        # line in same direction, same left point, line1 longer
1071        # 0:    ----->------
1072        # 1:    ------->--------
1073        # value should be line0:
1074        #       ----->------
[7276]1075        line0 = [P1, P3]
1076        line1 = [P1, P4]
[5942]1077        status, value = intersection(line0, line1)
1078        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1079                               (str(status), str(value)))
[6158]1080        self.failUnless(num.allclose(value, line0))
[5942]1081
1082        # line in same direction, same left point, line1 longer
1083        # 0:    -----<------
1084        # 1:    -------<--------
1085        # value should be line0:
1086        #       -----<------
[7276]1087        line0 = [P3, P1]
1088        line1 = [P4, P1]
[5942]1089        status, value = intersection(line0, line1)
1090        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1091                               (str(status), str(value)))
[6158]1092        self.failUnless(num.allclose(value, line0))
[5942]1093
1094        # line in opposite direction, same left point, line1 longer
1095        # 0:    ----->------
1096        # 1:    -------<--------
1097        # value should be line0:
1098        #       ----->------
[7276]1099        line0 = [P1, P3]
1100        line1 = [P4, P1]
[5942]1101        status, value = intersection(line0, line1)
1102        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1103                               (str(status), str(value)))
[6158]1104        self.failUnless(num.allclose(value, line0))
[5942]1105
1106        # line in opposite direction, same start point, line1 longer
1107        # 0:    -----<------
1108        # 1:    ------->--------
1109        # value should be line0:
1110        #       -----<------
[7276]1111        line0 = [P3, P1]
1112        line1 = [P1, P4]
[5942]1113        status, value = intersection(line0, line1)
1114        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1115                               (str(status), str(value)))
[6158]1116        self.failUnless(num.allclose(value, line0))
[5942]1117
[7276]1118        # ----------------------------------------------------------------------
[5942]1119
1120        # line in same direction, same left point, same right point
1121        # 0:    ------->--------
1122        # 1:    ------->--------
1123        # value should be line0 or line1:
1124        #       ------->--------
[7276]1125        line0 = [P1, P3]
1126        line1 = [P1, P3]
[5942]1127        status, value = intersection(line0, line1)
1128        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1129                               (str(status), str(value)))
[6158]1130        self.failUnless(num.allclose(value, line0))
[5942]1131
1132        # line in same direction, same left point, same right point
1133        # 0:    -------<--------
1134        # 1:    -------<--------
1135        # value should be line0 (or line1):
1136        #       -------<--------
[7276]1137        line0 = [P3, P1]
1138        line1 = [P3, P1]
[5942]1139        status, value = intersection(line0, line1)
1140        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1141                               (str(status), str(value)))
[6158]1142        self.failUnless(num.allclose(value, line0))
[5942]1143
1144        # line in opposite direction, same left point, same right point
1145        # 0:    ------->--------
1146        # 1:    -------<--------
1147        # value should be line0:
1148        #       ------->--------
[7276]1149        line0 = [P1, P3]
1150        line1 = [P3, P1]
[5942]1151        status, value = intersection(line0, line1)
1152        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1153                               (str(status), str(value)))
[6158]1154        self.failUnless(num.allclose(value, line0))
[5942]1155
1156        # line in opposite direction, same left point, same right point
1157        # 0:    -------<--------
1158        # 1:    ------->--------
1159        # value should be line0:
1160        #       -------<--------
[7276]1161        line0 = [P3, P1]
1162        line1 = [P1, P3]
[5942]1163        status, value = intersection(line0, line1)
1164        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1165                               (str(status), str(value)))
[6158]1166        self.failUnless(num.allclose(value, line0))
[5942]1167
[7276]1168        # ----------------------------------------------------------------------
[5942]1169
1170        # line in same direction, same right point, line1 longer
1171        # 0:        ----->------
1172        # 1:    ------->--------
1173        # value should be line0:
1174        #           ----->------
[7276]1175        line0 = [P2, P4]
1176        line1 = [P1, P4]
[5942]1177        status, value = intersection(line0, line1)
1178        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1179                               (str(status), str(value)))
[6158]1180        self.failUnless(num.allclose(value, line0))
[5942]1181
1182        # line in same direction, same right point, line1 longer
1183        # 0:        -----<------
1184        # 1:    -------<--------
1185        # value should be line0:
1186        #           -----<------
[7276]1187        line0 = [P4, P2]
1188        line1 = [P4, P1]
[5942]1189        status, value = intersection(line0, line1)
1190        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1191                               (str(status), str(value)))
[6158]1192        self.failUnless(num.allclose(value, line0))
[5942]1193
1194        # line in opposite direction, same right point, line1 longer
1195        # 0:        ----->------
1196        # 1:    -------<--------
1197        # value should be line0:
1198        #           ----->------
[7276]1199        line0 = [P2, P4]
1200        line1 = [P4, P1]
[5942]1201        status, value = intersection(line0, line1)
1202        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1203                               (str(status), str(value)))
[6158]1204        self.failUnless(num.allclose(value, line0))
[5942]1205
1206        # line in opposite direction, same right point, line1 longer
1207        # 0:        -----<------
1208        # 1:    ------->--------
1209        # value should be line0:
1210        #           -----<------
[7276]1211        line0 = [P4, P2]
1212        line1 = [P1, P4]
[5942]1213        status, value = intersection(line0, line1)
1214        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1215                               (str(status), str(value)))
[6158]1216        self.failUnless(num.allclose(value, line0))
[5942]1217
[7276]1218        # ----------------------------------------------------------------------
[5942]1219
1220        # line in same direction, same left point, line0 longer
1221        # 0:    ------->--------
1222        # 1:    ----->------
1223        # value should be line1:
1224        #       ----->------
[7276]1225        line0 = [P1, P4]
1226        line1 = [P1, P3]
[5942]1227        status, value = intersection(line0, line1)
1228        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1229                               (str(status), str(value)))
[6158]1230        self.failUnless(num.allclose(value, line1))
[5942]1231
1232        # line in same direction, same left point, line0 longer
1233        # 0:    -------<--------
1234        # 1:    -----<------
1235        # value should be line1:
1236        #       -----<------
[7276]1237        line0 = [P4, P1]
1238        line1 = [P3, P1]
[5942]1239        status, value = intersection(line0, line1)
1240        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1241                               (str(status), str(value)))
[6158]1242        self.failUnless(num.allclose(value, line1))
[5942]1243
1244        # line in opposite direction, same left point, line0 longer
1245        # 0:    ------->--------
1246        # 1:    -----<------
1247        # value should be line1:
1248        #       -----<------
[7276]1249        line0 = [P1, P4]
1250        line1 = [P3, P1]
[5942]1251        status, value = intersection(line0, line1)
1252        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1253                               (str(status), str(value)))
[6158]1254        self.failUnless(num.allclose(value, line1))
[5942]1255
1256        # line in opposite direction, same left point, line0 longer
1257        # 0:    -------<--------
1258        # 1:    ----->------
1259        # value should be line1:
1260        #       ----->------
[7276]1261        line0 = [P4, P1]
1262        line1 = [P1, P3]
[5942]1263        status, value = intersection(line0, line1)
1264        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1265                               (str(status), str(value)))
[6158]1266        self.failUnless(num.allclose(value, line1))
[5942]1267
[7276]1268        # ----------------------------------------------------------------------
[5942]1269
1270        # line in same direction, same right point, line0 longer
1271        # 0:    ------->--------
1272        # 1:        ----->------
1273        # value should be line1:
1274        #           ----->------
[7276]1275        line0 = [P1, P4]
1276        line1 = [P2, P4]
[5942]1277        status, value = intersection(line0, line1)
1278        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1279                               (str(status), str(value)))
[6158]1280        self.failUnless(num.allclose(value, line1))
[5942]1281
1282        # line in same direction, same right point, line0 longer
1283        # 0:    -------<--------
1284        # 1:        -----<------
1285        # value should be line1:
1286        #           -----<------
[7276]1287        line0 = [P4, P1]
1288        line1 = [P4, P2]
[5942]1289        status, value = intersection(line0, line1)
1290        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1291                               (str(status), str(value)))
[6158]1292        self.failUnless(num.allclose(value, line1))
[5942]1293
1294        # line in opposite direction, same right point, line0 longer
1295        # 0:    ------->--------
1296        # 1:        -----<------
1297        # value should be line1:
1298        #           -----<------
[7276]1299        line0 = [P1, P4]
1300        line1 = [P4, P2]
[5942]1301        status, value = intersection(line0, line1)
1302        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1303                               (str(status), str(value)))
[6158]1304        self.failUnless(num.allclose(value, line1))
[5942]1305
1306        # line in opposite direction, same right point, line0 longer
1307        # 0:    -------<--------
1308        # 1:        ----->------
1309        # value should be line1:
1310        #           ----->------
[7276]1311        line0 = [P4, P1]
1312        line1 = [P2, P4]
[5942]1313        status, value = intersection(line0, line1)
1314        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1315                               (str(status), str(value)))
[6158]1316        self.failUnless(num.allclose(value, line1))
[5942]1317
[7276]1318
[5942]1319    def test_intersection_bug_20081110_TR(self):
[5933]1320        """test_intersection_bug_20081110(self)
1321
[5942]1322        Test all cases in top-right quadrant
[5933]1323        """
1324
[5942]1325        # define 4 collinear points in top-right quadrant
[5933]1326        #    P1---P2---P3---P4
[5942]1327        P1 = [1.0, 1.0]
1328        P2 = [2.0, 2.0]
1329        P3 = [3.0, 3.0]
1330        P4 = [4.0, 4.0]
[7276]1331
[5942]1332        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1333        P1 = [1.0, 1.0+1.0e-9]
[7276]1334        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1335        P1 = [1.0, 1.0]
1336        P2 = [2.0, 2.0+1.0e-9]
[7276]1337        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1338        P2 = [2.0, 2.0]
1339        P3 = [3.0, 3.0+1.0e-9]
[7276]1340        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1341        P3 = [3.0, 3.0]
1342        P4 = [4.0, 4.0+1.0e-9]
[7276]1343        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1344
[5942]1345    def test_intersection_bug_20081110_TL(self):
1346        """test_intersection_bug_20081110(self)
[5933]1347
[5942]1348        Test all cases in top-left quadrant
1349        """
1350
1351        # define 4 collinear points in top-left quadrant
1352        #    P1---P2---P3---P4
1353        P1 = [-1.0, 1.0]
1354        P2 = [-2.0, 2.0]
1355        P3 = [-3.0, 3.0]
1356        P4 = [-4.0, 4.0]
[7276]1357
[5942]1358        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1359        P1 = [-1.0, 1.0+1.0e-9]
[7276]1360        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1361        P1 = [-1.0, 1.0]
1362        P2 = [-2.0, 2.0+1.0e-9]
[7276]1363        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1364        P2 = [-2.0, 2.0]
1365        P3 = [-3.0, 3.0+1.0e-9]
[7276]1366        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1367        P3 = [-3.0, 3.0]
1368        P4 = [-4.0, 4.0+1.0e-9]
[7276]1369        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1370
1371    def test_intersection_bug_20081110_BL(self):
[5933]1372        """test_intersection_bug_20081110(self)
1373
[5942]1374        Test all cases in bottom-left quadrant
[5933]1375        """
1376
[5942]1377        # define 4 collinear points in bottom-left quadrant
[5933]1378        #    P1---P2---P3---P4
[5942]1379        P1 = [-1.0, -1.0]
1380        P2 = [-2.0, -2.0]
1381        P3 = [-3.0, -3.0]
1382        P4 = [-4.0, -4.0]
[7276]1383
[5942]1384        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1385        P1 = [-1.0, -1.0+1.0e-9]
[7276]1386        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1387        P1 = [-1.0, -1.0]
1388        P2 = [-2.0, -2.0+1.0e-9]
[7276]1389        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1390        P2 = [-2.0, -2.0]
1391        P3 = [-3.0, -3.0+1.0e-9]
[7276]1392        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1393        P3 = [-3.0, -3.0]
1394        P4 = [-4.0, -4.0+1.0e-9]
[7276]1395        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1396
[5942]1397    def test_intersection_bug_20081110_BR(self):
[5933]1398        """test_intersection_bug_20081110(self)
1399
[5942]1400        Test all cases in bottom-right quadrant
[5933]1401        """
1402
[5942]1403        # define 4 collinear points in bottom-right quadrant
[5933]1404        #    P1---P2---P3---P4
[5942]1405        P1 = [1.0, -1.0]
1406        P2 = [2.0, -2.0]
1407        P3 = [3.0, -3.0]
1408        P4 = [4.0, -4.0]
[7276]1409
[5942]1410        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1411        P1 = [1.0, -1.0+1.0e-9]
[7276]1412        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1413        P1 = [1.0, -1.0]
1414        P2 = [2.0, -2.0+1.0e-9]
[7276]1415        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1416        P2 = [2.0, -2.0]
1417        P3 = [3.0, -3.0+1.0e-9]
[7276]1418        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1419        P3 = [3.0, -3.0]
1420        P4 = [4.0, -4.0+1.0e-9]
[7276]1421        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1422
[5942]1423    def test_intersection_bug_20081110_TR_TL(self):
[5933]1424        """test_intersection_bug_20081110(self)
1425
[5942]1426        Test all cases in top-right & top-left quadrant
[5933]1427        """
1428
[5942]1429        # define 4 collinear points, 1 in TL, 3 in TR
[7276]1430        #    P1-+-P2---P3---P4
[5942]1431        P1 = [-3.0, 1.0]
1432        P2 = [ 1.0, 5.0]
1433        P3 = [ 2.0, 6.0]
1434        P4 = [ 3.0, 7.0]
1435        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1436
[5942]1437        # define 4 collinear points, 2 in TL, 2 in TR
[7276]1438        #    P1---P2-+-P3---P4
[5942]1439        P1 = [-3.0, 1.0]
1440        P2 = [-2.0, 2.0]
1441        P3 = [ 2.0, 6.0]
1442        P4 = [ 3.0, 7.0]
1443        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1444
[7276]1445        # define 4 collinear points, 3 in TL, 1 in TR
1446        #    P1---P2---P3-+-P4
[5942]1447        P1 = [-3.0, 1.0]
1448        P2 = [-2.0, 2.0]
1449        P3 = [-1.0, 3.0]
1450        P4 = [ 3.0, 7.0]
1451        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1452
[5942]1453    def test_intersection_bug_20081110_TR_BL(self):
[5933]1454        """test_intersection_bug_20081110(self)
1455
[5942]1456        Test all cases in top-right & bottom-left quadrant
[5933]1457        """
1458
[5942]1459        # define 4 collinear points, 1 in BL, 3 in TR
[7276]1460        #    P1-+-P2---P3---P4
[5942]1461        P1 = [-4.0, -3.0]
1462        P2 = [ 1.0,  2.0]
1463        P3 = [ 2.0,  3.0]
1464        P4 = [ 3.0,  4.0]
1465        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1466
[5942]1467        # define 4 collinear points, 2 in TL, 2 in TR
[7276]1468        #    P1---P2-+-P3---P4
[5942]1469        P1 = [-4.0, -3.0]
1470        P2 = [-3.0, -2.0]
1471        P3 = [ 2.0,  3.0]
1472        P4 = [ 3.0,  4.0]
1473        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1474
[5942]1475        # define 4 collinear points, 3 in TL, 1 in TR
[7276]1476        #    P1---P2---P3-+-P4
[5942]1477        P1 = [-4.0, -3.0]
1478        P2 = [-3.0, -2.0]
1479        P3 = [-2.0, -1.0]
1480        P4 = [ 3.0,  4.0]
1481        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1482
[5942]1483    def test_intersection_bug_20081110_TR_BR(self):
1484        """test_intersection_bug_20081110(self)
[5939]1485
[5942]1486        Test all cases in top-right & bottom-right quadrant
[5939]1487        """
1488
[5942]1489        # define 4 collinear points, 1 in BR, 3 in TR
[7276]1490        #    P1-+-P2---P3---P4
[5942]1491        P1 = [ 1.0, -3.0]
1492        P2 = [ 5.0,  1.0]
1493        P3 = [ 6.0,  2.0]
1494        P4 = [ 7.0,  3.0]
1495        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5939]1496
[5942]1497        # define 4 collinear points, 2 in BR, 2 in TR
[7276]1498        #    P1---P2-+-P3---P4
[5942]1499        P1 = [ 1.0, -3.0]
1500        P2 = [ 2.0, -2.0]
1501        P3 = [ 6.0,  2.0]
1502        P4 = [ 7.0,  3.0]
1503        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5939]1504
[5942]1505        # define 4 collinear points, 3 in BR, 1 in TR
[7276]1506        #    P1---P2---P3-+-P4
[5942]1507        P1 = [ 1.0, -3.0]
1508        P2 = [ 2.0, -2.0]
1509        P3 = [ 3.0, -1.0]
1510        P4 = [ 7.0,  3.0]
1511        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5939]1512
[5940]1513
[5897]1514    def test_intersection_direction_invariance(self):
[7276]1515        """This runs through a number of examples and checks that
1516        direction of lines don't matter.
[5897]1517        """
[7276]1518
[5897]1519        line0 = [[0,0], [100,100]]
1520
1521        common_end_point = [20, 150]
[7276]1522
[5897]1523        for i in range(100):
1524            x = 20 + i * 1.0/100
1525
1526            line1 = [[x,0], common_end_point]
1527            status, p1 = intersection(line0, line1)
1528            assert status == 1
1529
1530            # Swap direction of line1
[7276]1531            line1 = [common_end_point, [x,0]]
[5897]1532            status, p2 = intersection(line0, line1)
[7276]1533            assert status == 1
[5897]1534
[7276]1535            msg = ('Orientation of line should not matter.\n'
1536                   'However, segment [%f,%f], [%f, %f]' %
1537                   (x, 0, common_end_point[0], common_end_point[1]))
1538            msg += ' gave %s, \nbut when reversed we got %s' % (p1, p2)
[6158]1539            assert num.allclose(p1, p2), msg
[5897]1540
1541            # Swap order of lines
1542            status, p3 = intersection(line1, line0)
[7276]1543            assert status == 1
[5897]1544            msg = 'Order of lines gave different results'
[6158]1545            assert num.allclose(p1, p3), msg
[5897]1546
1547    def test_no_intersection(self):
[7866]1548        """ Test 2 non-touching lines don't intersect. """
[5897]1549        line0 = [[-1,1], [1,1]]
1550        line1 = [[0,-1], [0,0]]
1551
1552        status, value = intersection(line0, line1)
1553        assert status == 0
1554        assert value is None
1555
1556    def test_intersection_parallel(self):
1557        line0 = [[-1,1], [1,1]]
1558        line1 = [[-1,0], [5,0]]
1559
1560        status, value = intersection(line0, line1)
[7276]1561        assert status == 4
[5897]1562        assert value is None
1563
1564        line0 = [[0,0], [10,100]]
1565        line1 = [[-10,5], [0,105]]
1566
1567        status, value = intersection(line0, line1)
[7276]1568        assert status == 4
1569        assert value is None
[5897]1570
1571    def test_intersection_coincide(self):
[7276]1572        """Test what happens when two lines partly coincide"""
[5897]1573
1574        # Overlap 1
1575        line0 = [[0,0], [5,0]]
1576        line1 = [[-3,0], [3,0]]
1577
1578        status, value = intersection(line0, line1)
1579        assert status == 2
[6158]1580        assert num.allclose(value, [[0,0], [3,0]])
[5897]1581
1582        # Overlap 2
1583        line0 = [[-10,0], [5,0]]
1584        line1 = [[-3,0], [10,0]]
1585
1586        status, value = intersection(line0, line1)
1587        assert status == 2
[7276]1588        assert num.allclose(value, [[-3, 0], [5,0]])
[5897]1589
1590        # Inclusion 1
1591        line0 = [[0,0], [5,0]]
1592        line1 = [[2,0], [3,0]]
1593
1594        status, value = intersection(line0, line1)
[7276]1595        assert status == 2
[6158]1596        assert num.allclose(value, line1)
[5897]1597
1598        # Inclusion 2
1599        line0 = [[1,0], [5,0]]
1600        line1 = [[-10,0], [15,0]]
1601
1602        status, value = intersection(line0, line1)
[7276]1603        assert status == 2
1604        assert num.allclose(value, line0)
[5897]1605
1606        # Exclusion (no intersection)
1607        line0 = [[-10,0], [1,0]]
1608        line1 = [[3,0], [15,0]]
1609
1610        status, value = intersection(line0, line1)
[7276]1611        assert status == 3
[5897]1612        assert value is None
1613
1614        # Try examples with some slope (y=2*x+5)
1615
1616        # Overlap
[7866]1617        line0 = [[0, 5], [7, 19]]
1618        line1 = [[1, 7], [10, 25]]
[5897]1619        status, value = intersection(line0, line1)
[7276]1620        assert status == 2
[6158]1621        assert num.allclose(value, [[1, 7], [7, 19]])
[5897]1622
1623        status, value = intersection(line1, line0)
1624        assert status == 2
[6158]1625        assert num.allclose(value, [[1, 7], [7, 19]])
[5897]1626
1627        # Swap direction
1628        line0 = [[7,19], [0,5]]
1629        line1 = [[1,7], [10,25]]
1630        status, value = intersection(line0, line1)
1631        assert status == 2
[6158]1632        assert num.allclose(value, [[7, 19], [1, 7]])
[5897]1633
1634        line0 = [[0,5], [7,19]]
1635        line1 = [[10,25], [1,7]]
1636        status, value = intersection(line0, line1)
1637        assert status == 2
[7276]1638        assert num.allclose(value, [[1, 7], [7, 19]])
[5897]1639
1640        # Inclusion
1641        line0 = [[1,7], [7,19]]
1642        line1 = [[0,5], [10,25]]
1643        status, value = intersection(line0, line1)
[7276]1644        assert status == 2
1645        assert num.allclose(value, [[1,7], [7, 19]])
[5897]1646
1647        line0 = [[0,5], [10,25]]
1648        line1 = [[1,7], [7,19]]
1649        status, value = intersection(line0, line1)
[7276]1650        assert status == 2
[6158]1651        assert num.allclose(value, [[1,7], [7, 19]])
[5897]1652
1653        line0 = [[0,5], [10,25]]
1654        line1 = [[7,19], [1,7]]
1655        status, value = intersection(line0, line1)
[7276]1656        assert status == 2
1657        assert num.allclose(value, [[7, 19], [1, 7]])
[5897]1658
1659
[7866]1660    def test_inside_polygon_main(self):
[7276]1661        # From real example (that failed)
1662        polygon = [[20,20], [40,20], [40,40], [20,40]]
1663        points = [[40, 50]]
1664        res = inside_polygon(points, polygon)
1665        assert len(res) == 0
[5897]1666
[7276]1667        polygon = [[20,20], [40,20], [40,40], [20,40]]
1668        points = [[25, 25], [30, 20], [40, 50], [90, 20], [40, 90]]
1669        res = inside_polygon(points, polygon)
1670        assert len(res) == 2
1671        assert num.allclose(res, [0,1])
[5897]1672
1673    def test_polygon_area(self):
[7866]1674        """ Test getting the area of a polygon. """
[7276]1675        # Simplest case: Polygon is the unit square
[5897]1676        polygon = [[0,0], [1,0], [1,1], [0,1]]
[7276]1677        assert polygon_area(polygon) == 1
[5897]1678
[7276]1679        # Simple case: Polygon is a rectangle
[5897]1680        polygon = [[0,0], [1,0], [1,4], [0,4]]
[7276]1681        assert polygon_area(polygon) == 4
[5897]1682
[7276]1683        # Simple case: Polygon is a unit triangle
[5897]1684        polygon = [[0,0], [1,0], [0,1]]
[7276]1685        assert polygon_area(polygon) == 0.5
[5897]1686
[7276]1687        # Simple case: Polygon is a diamond
[5897]1688        polygon = [[0,0], [1,1], [2,0], [1, -1]]
[7276]1689        assert polygon_area(polygon) == 2.0
1690
[6000]1691        # Complex case where numerical errors might occur
1692        polygon = [[314037.58727982, 6224952.2960092],
[7276]1693                   [314038.58727982, 6224952.2960092],
[6000]1694                   [314038.58727982, 6224953.2960092],
1695                   [314037.58727982, 6224953.2960092]]
[7276]1696        assert polygon_area(polygon) == 1.0
[5897]1697
1698    def test_poly_xy(self):
[7276]1699        # Simplest case: Polygon is the unit square
[5897]1700        polygon = [[0,0], [1,0], [1,1], [0,1]]
[7858]1701        x, y = _poly_xy(polygon)
[7276]1702        assert len(x) == len(polygon)+1
1703        assert len(y) == len(polygon)+1
1704        assert x[0] == 0
1705        assert x[1] == 1
1706        assert x[2] == 1
1707        assert x[3] == 0
1708        assert y[0] == 0
1709        assert y[1] == 0
1710        assert y[2] == 1
1711        assert y[3] == 1
[5897]1712
[7276]1713    # Arbitrary polygon
[5897]1714        polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]]
[7858]1715        x, y = _poly_xy(polygon)
[7276]1716        assert len(x) == len(polygon)+1
1717        assert len(y) == len(polygon)+1
1718        assert x[0] == 1
1719        assert x[1] == 1
1720        assert x[2] == 100
1721        assert x[3] == 1
1722        assert x[4] == 3
1723        assert y[0] == 5
1724        assert y[1] == 1
1725        assert y[2] == 10
1726        assert y[3] == 10
1727        assert y[4] == 6
[5897]1728
[7841]1729
1730    def test_plot_polygons(self):
[5897]1731        import os
[7276]1732
1733        # Simplest case: Polygon is the unit square
[5897]1734        polygon1 = [[0,0], [1,0], [1,1], [0,1]]
1735        polygon2 = [[1,1], [2,1], [3,2], [2,2]]
[7858]1736        plot_polygons([polygon1, polygon2], figname='test1')
[5897]1737
[7276]1738        # Another case
[5897]1739        polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
[7866]1740        plot_polygons([polygon2, polygon3], figname='test2')
[5897]1741
[7858]1742        for file in ['test1.png', 'test2.png']:
1743            assert os.access(file, os.R_OK)
1744            os.remove(file)
[5897]1745
[7858]1746
[5897]1747    def test_inside_polygon_geospatial(self):
[7866]1748        """ Test geospatial coords inside poly. """
[5897]1749        #Simplest case: Polygon is the unit square
[7866]1750        polygon_absolute = [[0, 0], [1, 0], [1, 1], [0, 1]]
[7276]1751        poly_geo_ref = Geo_reference(57, 100, 100)
[5897]1752        polygon = poly_geo_ref.change_points_geo_ref(polygon_absolute)
[7276]1753        poly_spatial = Geospatial_data(polygon, geo_reference=poly_geo_ref)
1754
[5897]1755        points_absolute = (0.5, 0.5)
[7276]1756        points_geo_ref = Geo_reference(57, 78, -56)
[5897]1757        points = points_geo_ref.change_points_geo_ref(points_absolute)
[7276]1758        points_spatial = Geospatial_data(points, geo_reference=points_geo_ref)
1759
[5897]1760        assert is_inside_polygon(points_absolute, polygon_absolute)
1761        assert is_inside_polygon(ensure_numeric(points_absolute),
1762                                 ensure_numeric(polygon_absolute))
[7276]1763        assert is_inside_polygon(points_absolute, poly_spatial)
1764        assert is_inside_polygon(points_spatial, poly_spatial)
1765        assert is_inside_polygon(points_spatial, polygon_absolute)
[5897]1766
[7276]1767        assert is_inside_polygon(points_absolute, polygon_absolute)
[5897]1768
[7866]1769    def test_decimate_polygon(self):
1770        from polygon import decimate_polygon
[5897]1771        polygon = [[0,0], [10,10], [15,5], [20, 10],
1772                   [25,0], [30,10], [40,-10], [35, -5]]
1773
1774        dpoly = decimate_polygon(polygon, factor=2)
1775
1776        assert len(dpoly)*2==len(polygon)
1777
[7276]1778
[6189]1779    def test_interpolate_polyline(self):
1780        """test_interpolate_polyline(self):
[7276]1781
1782        This test is added under the assumption that the function
1783        interpolate_polyline implemented by John Jakeman works.
1784        It has been exercised somewhat by tests of sts boundary,
1785        but never before separately.
[6189]1786        """
[7276]1787
[6189]1788        f = num.array([58.06150614, 58.06150614, 58.06150614])
1789        vertex_coordinates = num.array([[0., 0., ],
1790                                        [4.04092634, 1106.11074699],
1791                                        [8.08836552, 2212.16910609]])
1792        gauge_neighbour_id = [1, 2, -1]
1793        point_coordinates = num.array([[2.21870766e+03, 1.09802864e+03],
1794                                       [1.62739645e+03, 2.20626983e+03],
1795                                       [5.20084967e+02, 2.21030386e+03],
1796                                       [6.06464546e+00, 1.65913993e+03],
1797                                       [1.61934862e+03, -5.88143836e+00],
1798                                       [5.11996623e+02, -1.85956061e+00],
1799                                       [2.02046270e+00, 5.53055373e+02]])
[7276]1800
[6189]1801        z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
[7276]1802
1803        z = interpolate_polyline(f, vertex_coordinates,
1804                                 gauge_neighbour_id, point_coordinates)
[6189]1805        assert num.allclose(z, z_ref)
[7276]1806
[6189]1807        # Another f
1808        f = num.array([58.06150614, 158.06150614, 258.06150614])
[7276]1809        z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]
1810        z = interpolate_polyline(f, vertex_coordinates,
1811                                 gauge_neighbour_id, point_coordinates)
1812        assert num.allclose(z, z_ref)
[6189]1813
1814        # Other and simpler numbers
[7276]1815        f = num.array([1, 5, 13])
[6189]1816        vertex_coordinates = num.array([[0., 0., ],
1817                                        [4., 4.],
[7276]1818                                        [8., 8.]])
[6189]1819        point_coordinates = num.array([[0.1, 0.1],
1820                                       [3.5, 3.5],
1821                                       [4.0, 4.0],
1822                                       [5.2, 5.2],
1823                                       [7.0, 7.0],
1824                                       [8.3, 8.3]])
1825        gauge_neighbour_id = [1, 2, -1]
[7276]1826
1827        z = interpolate_polyline(f, vertex_coordinates,
1828                                 gauge_neighbour_id, point_coordinates)
[6189]1829        z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
[7276]1830        assert num.allclose(z, z_ref)
1831
[6189]1832        # Test exception thrown for one point
[7276]1833        f = num.array([5])
1834        vertex_coordinates = num.array([[4., 4.]])
1835        try:
1836            z = interpolate_polyline(f, vertex_coordinates,
1837                                     gauge_neighbour_id, point_coordinates)
[6189]1838        except Exception:
1839            pass
1840        else:
1841            raise Exception, 'One point should have raised exception'
1842
1843        # More polyline nodes
[7276]1844        data = num.array([1, 5, 13, 12, 6, 29])
[6189]1845        polyline_nodes = num.array([[0., 0.],
1846                                    [4., 4.],
1847                                    [8., 8.],
1848                                    [10., 10.],
1849                                    [10., 5.],
1850                                    [10., 0.]])
1851        point_coordinates = num.array([[0.1, 0.1],
1852                                       [3.5, 3.5],
1853                                       [4.0, 4.0],
1854                                       [5.2, 5.2],
1855                                       [7.0, 7.0],
1856                                       [8.3, 8.3],
1857                                       [10., 10.],
1858                                       [10., 9.],
1859                                       [10., 7.1],
1860                                       [10., 4.3],
1861                                       [10., 1.0]])
1862        gauge_neighbour_id = [1, 2, 3, 4, 5, -1]
[7276]1863        z = interpolate_polyline(data, polyline_nodes,
1864                                 gauge_neighbour_id, point_coordinates)
[6189]1865        z_ref = [1.1, 4.5, 5., 7.4, 11., 12.85, 12., 10.8, 8.52, 9.22, 24.4]
[7276]1866        assert num.allclose(z, z_ref)
[6535]1867
[7276]1868
1869    def test_is_inside_triangle_more(self):
[7866]1870        """ Test if points inside triangles are detected correctly. """
[6535]1871        res = is_inside_triangle([0.5, 0.5], [[ 0.5,  0. ],
1872                                              [ 0.5,  0.5],
1873                                              [ 0.,   0. ]])
1874        assert res is True
1875
1876        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1877                                 [[ 0.5,  0. ], [ 0.5, 0.5], [0., 0.]])
1878        assert res is False
[7276]1879
[6535]1880        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1881                                 [[1., 0.], [1., 0.5], [0.5, 0.]])
[7276]1882        assert res is False
1883
1884
[6535]1885        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1886                                 [[0.5, 0.5], [0.5, 0.], [1., 0.5]])
[7276]1887        assert res is True
[6535]1888
[7276]1889
[6535]1890        res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
1891                                 [[0.5, 0.], [0.5, 0.5], [0., 0.]])
1892        assert res is False
1893
[7276]1894
[6535]1895        res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
1896                                 [[0., 0.5], [0., 0.], [0.5, 0.5]])
[7276]1897        assert res is True
1898
[6535]1899        res = is_inside_triangle([0.69999999999999996, 0.69999999999999996],
1900                                 [[0.5, 0.], [0.5, 0.5], [0., 0.]])
[7276]1901        assert res is False
[6535]1902
1903        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1904                                 [[0.25, 0.5], [0.25, 0.25], [0.5, 0.5]])
1905        assert res is False
[6544]1906
[7276]1907
1908        res = is_inside_triangle([10, 3],
[6544]1909                                 [[0.1, 0.],
1910                                  [0.1, 0.08333333],
1911                                  [0., 0.]])
1912        assert res is False
[7276]1913
[7687]1914
[7686]1915    def test_is_polygon_complex(self):
[7866]1916        """ Test a concave and a complex poly with is_complex, to make
1917            sure it can detect self-intersection.
[7690]1918        """
[7686]1919        concave_poly = [[0, 0], [10, 0], [5, 5], [10, 10], [0, 10]]
[7866]1920        complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], \
1921                            [0, 10]]
[7687]1922
[7690]1923        assert not is_complex(concave_poly)
1924        assert is_complex(complex_poly)
[7687]1925
1926    def test_is_polygon_complex2(self):
[7866]1927        """ Test a concave and a complex poly with is_complex, to make sure it
1928            can detect self-intersection. This test uses more complicated
1929            polygons.
[7690]1930        """
[7866]1931        concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7, 6], [10, 10], \
1932                        [1, 5], [0, 10]]
1933        complex_poly = [[0, 0], [12, 12], [10, 0], [5, 5], [3,18], [4, 15], \
1934                        [5, 7], [10, 10], [0, 10], [16, 12]]
[7687]1935
[7690]1936        assert not is_complex(concave_poly)
1937        assert is_complex(complex_poly)
[7686]1938
[7276]1939################################################################################
1940
[5897]1941if __name__ == "__main__":
[7276]1942    suite = unittest.makeSuite(Test_Polygon,'test')
[5897]1943    runner = unittest.TextTestRunner()
1944    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.