python 实现 射线法 判断一个点在图形区域内外

关于射线法可以参考下面这个blog ，写得非常详细：http://www.cnblogs.com/mazhenyu/p/3800638.html

下面是我用python 实现的代码：

# -*-encoding:utf-8 -*-
# file:class.py
#

"""
信息楼
0	123.425658,41.774177
1	123.425843,41.774166
2	123.425847,41.774119
3	123.42693,41.774062
4	123.426943,41.774099
5	123.427118,41.774089
6	123.427066,41.773548
7	123.426896,41.773544
8	123.426916,41.773920
9	123.425838,41.773965
10	123.425804,41.773585
11	123.425611,41.773595

图书馆
0	123.425649,41.77303
1	123.426656,41.772993
2	123.426611,41.772398
3	123.425605,41.772445
"""


class Point:
	lat = ''
	lng = ''
	
	def __init__(self,lat,lng):
		self.lat = lat #纬度
		self.lng = lng #经度

	def show(self):
		print self.lat," ",self.lng


#将信息楼的边界点实例化并存储到points1里
point0 = Point(123.425658,41.774177)
point1 = Point(123.425843,41.774166)
point2 = Point(123.425847,41.774119)
point3 = Point(123.42693,41.774062)
point4 = Point(123.426943,41.774099)
point5 = Point(123.427118,41.774089)
point6 = Point(123.427066,41.773548)
point7 = Point(123.426896,41.773544)
point8 = Point(123.426916,41.773920)
point9 = Point(123.425838,41.773961)
point10 = Point(123.425804,41.773585)
point11 = Point(123.425611,41.773595)

points1 = [point0,point1,point2,point3,
		   point4,point5,point6,point7,
		   point8,point9,point10,point11,
		  ]


#将图书馆的边界点实例化并存储到points2里
point0 = Point(123.425649,41.77303)
point1 = Point(123.426656,41.772993)
point2 = Point(123.426611,41.772398)
point3 = Point(123.425605,41.772445)

points2 = [point0,point1,point2,point3]


'''
将points1和points2存储到points里,
points可以作为参数传入
'''
points = [points1,points2]


'''
输入一个测试点，这个点通过GPS产生
建议输入三个点测试
在信息学馆内的点:123.4263790000,41.7740520000 123.42699,41.773592  
在图书馆内的点:  123.4261550000,41.7726740000 123.42571,41.772499 123.425984,41.772919  
不在二者内的点:  123.4246270000,41.7738130000
在信息学馆外包矩形内，但不在信息学馆中的点：123.4264060000,41.7737860000
'''
#lat = raw_input(please input lat)
#lng = raw_input(please input lng)
lat = 123.42699
lng = 41.773592
point = Point(lat,lng)

debug = raw_input("请输入debug")
if debug == '1':
	debug = True
else:
	debug = False

#求外包矩形
def getPolygonBounds(points):
	length = len(points)
	#top down left right 都是point类型
	top = down = left = right = points[0]
	for i in range(1,length):
		if points[i].lng > top.lng:
			top = points[i]
		elif points[i].lng < down.lng:
			down = points[i]
		else:
			pass
		if points[i].lat > right.lat:
			right = points[i]
		elif points[i].lat < left.lat:
			left = points[i]
    	else:
    		pass

	point0 = Point(left.lat,top.lng)
	point1 = Point(right.lat,top.lng)
	point2 = Point(right.lat,down.lng)
	point3 = Point(left.lat,down.lng)
	polygonBounds = [point0,point1,point2,point3]
 	return polygonBounds

#测试求外包矩形的一段函数
if debug:
	poly1 = getPolygonBounds(points[0])
	print "第一个建筑的外包是："
	for i in range(0,len(poly1)):
		poly1[i].show()	
	poly2 = getPolygonBounds(points[1])
	print "第二个建筑的外包是："
	for i in range(0,len(poly2)):
		poly2[i].show()	


#判断点是否在外包矩形外
def isPointInRect(point,polygonBounds):
	if point.lng >= polygonBounds[3].lng and \
	   point.lng <= polygonBounds[0].lng and \
	   point.lat >= polygonBounds[3].lat and \
	   point.lat <= polygonBounds[2].lat:\
	   return True
	else:
		return False

#测试是否在外包矩形外的代码
if debug:
	if(isPointInRect(point,poly1)):
		print "在信息外包矩形内"
	else:
		print "在信息外包矩形外"

	if(isPointInRect(point,poly2)):
		print "在图书馆外包矩形内"
	else:
		print "在图书馆外包矩形外"



#采用射线法，计算测试点是否任意一个建筑内
def isPointInPolygon(point,points):
	#定义在边界上或者在顶点都建筑内
	Bound = Vertex = True
	count = 0
	precision = 2e-10
	
	#首先求外包矩形
	polygonBounds = getPolygonBounds(points)

	#然后判断是否在外包矩形内，如果不在，直接返回false
	if not isPointInRect(point, polygonBounds):
		if debug:
			print "在外包矩形外"
		return False
	else:
		if debug:
			print "在外包矩形内"

	length = len(points)
	p = point
	p1 = points[0]
	for i in range(1,length):
		if p.lng == p1.lng and p.lat == p1.lat:
			if debug:
				print "Vertex1"
			return Vertex

		p2 = points[i % length]
		if p.lng == p2.lng and p.lat == p2.lat:
			if dubug:
				print "Vertex2"
			return Vertex

		if debug:	
			print i-1,i
			print "p:"
			p.show()
			print "p1:"
			p1.show()
			print "p2:"
			p2.show()

		if p.lng < min(p1.lng,p2.lng) or \
			p.lng > max(p1.lng,p2.lng) or \
			p.lat > max(p1.lat,p2.lat): 
			p1 = p2
			if debug:
				print "Outside"
			continue

		elif p.lng > min(p1.lng,p2.lng) and \
			p.lng < max(p1.lng,p2.lng):
			if p1.lat == p2.lat:
				if p.lat == p1.lat and \
					p.lng > min(p1.lng,p2.lng) and \
					p.lng < max(p1.lng,p2.lng):
					return Bound
				else:
					count = count + 1
					if debug:
						print "count1:",count
					continue
			if debug:
				print "into left or right"	  		

			a = p2.lng - p1.lng
			b = p1.lat - p2.lat
			c = p2.lat * p1.lng - p1.lat * p2.lng
			d = a * p.lat + b * p.lng + c
			if p1.lng < p2.lng and p1.lat > p2.lat or \
			   p1.lng < p2.lng and p1.lat < p2.lat:	
				if d < 0:
					count = count + 1
					if debug:
						print "count2:",count
				elif d > 0:
					p1 = p2
					continue
				elif abs(p.lng-d) < precision :
					return Bound
			else :	     	
				if d < 0:
					p1 = p2
					continue
				elif d > 0:
					count = count + 1
					if debug:
						print "count3:",count
				elif abs(p.lng-d) < precision :
					return Bound
		else:
			if p1.lng == p2.lng:
				if p.lng == p1.lng and \
					p.lat > min(p1.lat,p2.lat) and \
					p.lat < max(p1.lat,p2.lat):
				   return Bound
			else:
				p3 = points[(i+1) % length]
				if p.lng < min(p1.lng,p3.lng) or \
					p.lng > max(p1.lng,p3.lng):
					count = count + 2
					if debug:
						print "count4:",count
				else:
					count = count + 1
					if debug:
						print "count5:",count	 
		p1 = p2
	if count % 2 == 0 :
		return False
	else :
		return True



length = len(points)
flag = 0
for i in range(length):
	if isPointInPolygon(point,points[i]):
		print "你刚才输入的点在第 %d 个建筑里" % (i+1)
		print "然后根据i值，可以读出建筑名，或者修改传入的points参数"
		break
	else:
		flag = flag + 1
	
if flag == length:
	print "在头 %d 建筑外" % (i+1)