1. 程式人生 > 程式設計 >Python求凸包及多邊形面積教程

Python求凸包及多邊形面積教程

一般有兩種演算法來計算平面上給定n個點的凸包:Graham掃描法(Graham's scan),時間複雜度為O(nlgn);Jarvis步進法(Jarvis march),時間複雜度為O(nh),其中h為凸包頂點的個數。這兩種演算法都按逆時針方向輸出凸包頂點。

Graham掃描法

用一個棧來解決凸包問題,點集Q中每個點都會進棧一次,不符合條件的點會被彈出,演算法終止時,棧中的點就是凸包的頂點(逆時針順序在邊界上)。

演算法步驟如下圖:

Python求凸包及多邊形面積教程

Python求凸包及多邊形面積教程

Python求凸包及多邊形面積教程

Python求凸包及多邊形面積教程

Python求凸包及多邊形面積教程

Python求凸包及多邊形面積教程

import sys
import math
import time
import random

#獲取基準點的下標,基準點是p[k]
def get_leftbottompoint(p):
 k = 0
 for i in range(1,len(p)):
  if p[i][1] < p[k][1] or (p[i][1] == p[k][1] and p[i][0] < p[k][0]):
   k = i
 return k

#叉乘計算方法
def multiply(p1,p2,p0):
 return (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p2[0] - p0[0]) * (p1[1] - p0[1])

#獲取極角,通過求反正切得出,考慮pi/2的情況
def get_arc(p1,p0):
 # 相容sort_points_tan的考慮
 if (p1[0] - p0[0]) == 0:
  if ((p1[1] - p0[1])) == 0:
   return -1;
  else:
   return math.pi / 2
 tan = float((p1[1] - p0[1])) / float((p1[0] - p0[0]))
 arc = math.atan(tan)
 if arc >= 0:
  return arc
 else:
  return math.pi + arc

#對極角進行排序,排序結果list不包含基準點
def sort_points_tan(p,pk):
 p2 = []
 for i in range(0,len(p)):
  p2.append({"index": i,"arc": get_arc(p[i],pk)})
 #print('排序前:',p2)
 p2.sort(key=lambda k: (k.get('arc')))
 #print('排序後:',p2)
 p_out = []
 for i in range(0,len(p2)):
  p_out.append(p[p2[i]["index"]])
 return p_out

def convex_hull(p):
 p=list(set(p))
 #print('全部點:',p)
 k = get_leftbottompoint(p)
 pk = p[k]
 p.remove(p[k])
 #print('排序前去除基準點的所有點:',p,'基準點:',pk)

 p_sort = sort_points_tan(p,pk) #按與基準點連線和x軸正向的夾角排序後的點座標
 #print('其餘點與基準點夾角排序:',p_sort)
 p_result = [pk,p_sort[0]]

 top = 2
 for i in range(1,len(p_sort)):
  #####################################
  #叉乘為正,向前遞迴刪點;叉乘為負,序列追加新點
  while(multiply(p_result[-2],p_sort[i],p_result[-1]) > 0):
   p_result.pop()
  p_result.append(p_sort[i]) 
 return p_result#測試
if __name__ == '__main__':
 pass
 test_data = [(220,-100),(0,0),(-40,-170),(240,50),(-160,150),(-210,-150)]
 print(test_data)

 result = convex_hull(test_data)
 print(result)
 t=0

import matplotlib.pyplot as plt
x1=[]
y1=[]
for i in range(len(test_data)):
 ri=test_data[i]
 #print(ri)
 x1.append(ri[0])
 y1.append(ri[1])

plt.plot(x1,y1,linestyle=' ',marker='.')


xx=[]
yy=[]
for i in range(len(result)):
 ri=result[i]
 #print(ri)
 xx.append(ri[0])
 yy.append(ri[1])

plt.plot(xx,yy,marker='*')

Python求凸包及多邊形面積教程

計算多邊形面積

(1)順時針給定構成凸包的n個點座標,叉乘法求多邊形面積:

Python求凸包及多邊形面積教程

def GetAreaOfPolyGonbyVector(points):
 # 基於向量叉乘計算多邊形面積
 area = 0
 if(len(points)<3):
  raise Exception("error")

 for i in range(0,len(points)-1):
  p1 = points[i]
  p2 = points[i + 1]

  triArea = (p1[0]*p2[1] - p2[0]*p1[1])/2
  #print(triArea)
  area += triArea

 fn=(points[-1][0]*points[0][1]-points[0][0]*points[-1][1])/2
 #print(fn)
 return abs(area+fn)

points = []
x = [1,3,2]
y = [1,2,2] 
#[(1,1),(3,(5,3),5),(1,3)] 
# x=[1,5,1]
# y=[1,1,3]
for index in range(len(x)):
 points.append((x[index],y[index]))
area = GetAreaOfPolyGonbyVector(points)
print(area)
#print(math.ceil(area))

(2)順時針給定構成凸包的n個點經緯度座標,先將經緯度座標轉化成凸多邊形的邊的經緯度距離,利用海倫公式求多邊形面積:

from geopy.distance import vincenty
import math
def HeronGetAreaOfPolyGonbyVector(points):
 # 基於海倫公式計算多邊形面積
 area = 0
 if(len(points)<3):
  raise Exception("error")

 pb=((points[-1][0]+points[0][0])/2,(points[-1][1]+points[0][1])/2) #基準點選為第一個點和最後一個點連線邊上的中點

 for i in range(0,len(points)-1):
  p1 = points[i]
  p2 = points[i + 1]

  db1 = vincenty(pb,p1).meters #根據維度轉化成經緯度距離
  d12 = vincenty(p1,p2).meters
  d2b = vincenty(p2,pb).meters
  #print(db1,d12,d2b)

  hc = (db1+d12+d2b)/2 #db1是基準點和p1的距離,d12是p1和p2的距離,d2b是p2和基準點距離
  #print(hc,hc-db1,hc-d12,hc-d2b)
  triArea = math.sqrt(hc*(hc-db1)*(hc-d12)*(hc-d2b)) 
  #print(triArea)
  area += triArea

 return area


points = []
x = [1,y[index]))

area = HeronGetAreaOfPolyGonbyVector(points)
print(area)
#print(math.ceil(area))

Graham程式原理

(1)基準點的確認原則:

有唯一的某個點縱座標最小,該點為基準點;

不止一個點的縱座標最小,選這些點裡最靠左的為基準點

(2)計算叉乘【後續利用叉乘正負判斷夾角是否大於180o】:

Python求凸包及多邊形面積教程

(3)獲取極角,通過求反正切得出:

若橫縱座標都相等(兩點相同),返回-1;

若橫座標相等/縱座標不相等(兩點連線垂直y軸),返回 Python求凸包及多邊形面積教程

Python求凸包及多邊形面積教程

(4)對極角進行排序,排序結果list不包含基準點:

p2=[{"index":0,"arc":get_arc(p[0],p[k])},
 {"index":1,"arc":get_arc(p[1],p[k])},
 ···
 {"index":k-1,"arc":get_arc(p[k-1],p[k])},
 {"index":k+1,"arc":get_arc(p[k+1],p[k])},
 ···
 {"index":n,"arc":get_arc(p[n],p[k])}]
#get_arc(p[0],p[k])即獲得p[0]點與基準點p[k]連線的極角(與x軸正向夾角)
#根據p2的“arc”鍵的值從小到大排序,最後輸出按該角度值排序對應順序的各個點

(5)逆時針確定凸多邊形:

Python求凸包及多邊形面積教程

主要是找角度是否大於180o——差乘正負——點進出棧順序三者關係

Python求凸包及多邊形面積教程

...一直遍歷到最後一個點...一直遍歷到最後一個點

規律:叉乘>0,夾角小於180o,遞歸向前刪點;叉乘<0,夾角大於180o,不刪點,加入新點,向後遍歷叉乘>0,夾角小於180o,遞歸向前刪點;叉乘<0,夾角大於180o,不刪點,加入新點,向後遍歷

注意:(a)上述給非基準點按極角從到大小排號時,有兩個及以上點“和基準點連線構成的極角”相等時,這些點的排號挨著但是沒有固定順序,這點並不影響演算法給出凸包的準確性。(b)對排號最後的一個點,掃描演算法裡沒有任何刪除該點的機制,但是這點也不影響演算法給出凸包的準確性。(c)上述程式需要額外加入,判斷結束棧內點數小於3和篩選凸包前點數小於3,不能計算多邊形面積的情況,可以直接給這種情況賦值0返回。

以上這篇Python求凸包及多邊形面積教程就是小編分享給大家的全部內容了,希望能給大家一個參考,也希望大家多多支援我們。