source: branches/numpy/anuga/utilities/test_polygon.py @ 6410

Last change on this file since 6410 was 6410, checked in by rwilson, 15 years ago

numpy changes.

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