Context Navigation

source: branches/source_numpy_conversion/anuga/utilities/test_polygon.py @ 6982

Last change on this file since 6982 was 5948, checked in by rwilson, 16 years ago
More NumPy? changes.
File size: 59.6 KB

Line
1	#!/usr/bin/env python
2
3
4	import unittest
5	import numpy
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)], 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	#-------------------------------------------------------------
1691	if __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.

Download in other formats: