source: anuga_core/source/anuga/utilities/test_polygon.py @ 6215

Last change on this file since 6215 was 6215, checked in by ole, 15 years ago

Removed references to _range in sww files in preparation for new viewer and ticket:192
Current extrema still needs to be removed and replaced by vertex based extrema.

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