Context Navigation

source: trunk/anuga_core/source/anuga/geometry/test_polygon.py @ 8038

Last change on this file since 8038 was 8038, checked in by habili, 13 years ago
added tests for polyline_overlap routines
File size: 72.0 KB

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

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

Download in other formats: