source: anuga_core/source_numpy_conversion/anuga/utilities/test_polygon.py @ 5948

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

More NumPy? changes.

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