Context Navigation

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

Last change on this file since 6000 was 6000, checked in by ole, 15 years ago
Protected calculation of polygon area against rounding errors by shifting it's origin to (0,0) for the purpose of the area computation.
File size: 57.5 KB

Line
1	#!/usr/bin/env python
2
3
4	import unittest
5	from Numeric import zeros, array, allclose
6	from math import sqrt, pi
7	from anuga.utilities.numerical_tools import ensure_numeric
8	from anuga.utilities.system_tools import get_pathname_from_package
9
10	from polygon import *
11	from anuga.coordinate_transforms.geo_reference import Geo_reference
12	from anuga.geospatial_data.geospatial_data import Geospatial_data
13
14	def test_function(x, y):
15	return x+y
16
17	class 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 allclose(indices, [3,2,1,0])
358
359	indices = outside_polygon(points, U, closed = True)
360	assert allclose(indices, [0,1,2,3])
361
362	indices = inside_polygon(points, U, closed = True)
363	assert 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 allclose(indices, [0,1,2])
378
379	indices = outside_polygon(points, U, closed = True)
380	assert allclose(indices, [])
381
382	indices = inside_polygon(points, U, closed = True)
383	assert 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 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 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 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 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 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 allclose(value, [12.0, 6.0])
606
607	# Swap order of lines
608	status, value = intersection(line1, line0)
609	assert status == 1
610	assert 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 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 allclose(value, [14.068965517, 7.0344827586])
626
627	# Swap order of lines
628	status, value = intersection(line1, line0)
629	assert status == 1
630	assert 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 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 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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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 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 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 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 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 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 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 allclose(value, [[1, 7], [7, 19]])
1466
1467	status, value = intersection(line1, line0)
1468	assert status == 2
1469	assert 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 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 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 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 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 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 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	#-------------------------------------------------------------
1693	if __name__ == "__main__":
1694	suite = unittest.makeSuite(Test_Polygon,'test')
1695	# suite = unittest.makeSuite(Test_Polygon,'test_intersection_bug_20081111')
1696	runner = unittest.TextTestRunner()
1697	runner.run(suite)
1698
1699
1700
1701

Note: See TracBrowser for help on using the repository browser.

Download in other formats: