找到的散点线性拟合方法都是基于最小二乘法的（numpy.polyfit() ， scipy.optimize())
以下是根据 GB/T 18459-2001中附录 A2 提供的独立线性度拟合方法，求得的最佳拟合直线
import math
def find_line(x0, y0):
'''
根据散点求得端基直线k,b,并得到每点对端基直线的偏差dy
:param x0: x坐标数组
:param y0: y坐标数组
:return: 每点的偏差dy
'''
# 求得端基直线的k,b
k = (y0[-1] - y0[0]) / (x0[-1] - x0[0])
d = y0[-1] - (k * x0[-1])
# 根据端基直线的k,b,求得每点对端基直线的偏差dy
dy = [(y0[i] - ((k * x0[i]) + d)) for i in range(len(x0))]
return dy
def find_poly(x_lst, y_lst):
'''
找到凸多边形,可以包含全部偏差点dy
:param x_lst:x数组
:param y_lst:偏差点dy数组
:return:最佳凸多边形各个点
'''
g_x = []
g_y = []
# 从起始点开始，向右（x增大方向），找到最大斜率点，再从最大斜率点开始，向右继续寻找
start_index = 0
x_tmp = x_lst[1:]
y_tmp = y_lst[1:]
while len(x_tmp) > 0:
k_lst = [((j - y_lst[start_index]) / (i - x_lst[start_index])) for i, j in zip(x_tmp, y_tmp)]
start_index = k_lst.index(max(k_lst)) + 1 + start_index
g_x.append(x_lst[start_index])
g_y.append(y_lst[start_index])
x_tmp = x_lst[start_index + 1:]
y_tmp = y_lst[start_index + 1:]
# 从最末点开始，向左（x减小方向），也找到最大斜率点（有人说找最小斜率点，但是结果算出来不对），再从最大斜率点开始，向左继续寻找
start_index = len(x_lst) - 1
x_tmp = x_lst[:start_index]
y_tmp = y_lst[:start_index]
while len(x_tmp) > 0:
k_lst = [((j - y_lst[start_index]) / (i - x_lst[start_index])) for i, j in zip(x_tmp, y_tmp)]
start_index = k_lst.index(max(k_lst))
g_x.append(x_lst[start_index])
g_y.append(y_lst[start_index])
x_tmp = x_lst[:start_index]
y_tmp = y_lst[:start_index]
return g_x, g_y
def find_cross(p1, p2, p3):
'''
根据一个点p1和一条直线(p2和p3的连线)，
求得该点对直线的铅垂线（平行于纵轴坐标的直线）和直线的交点
:param p1: 点
:param p2: 线上一点
:param p3: 线上另一点
:return: 铅垂线与线（p2-p3）的交点
'''
k = (p2[1] - p3[1]) / (p2[0] - p3[0])
b = p2[1] - k * p2[0]
p4 = [p1[0], (k * p1[0]) + b]
return p4
def find_perfect_line(a, b, x0, y0):
'''
根据凸多边形的每个点，找到凸多边形内最长的一根铅垂线，
与最长垂线相交的直线l1的斜率就是最佳直线的斜率
过最长铅垂线的中点作直线l2平行于直线l1，直线l2为最佳直线
:param a:凸多边形的x轴数组
:param b:凸多边形的y轴数组
:param x0: 实际x轴数组
:param y0: 实际y轴数组
:return:最佳直线的斜率k, 截距b
'''
dic = {}
point_lst = [[i, j] for i, j in zip(a, b)]
for i in range(len(a)):
p1 = (a[i], b[i]) # 凸多边形的每个点坐标
rp1 = (x0[x0.index(p1[0])], y0[x0.index(p1[0])])
for j in range(len(a)):
p2 = (a[j], b[j]) # 凸多边形的每个点坐标
rp2 = (x0[x0.index(p2[0])], y0[x0.index(p2[0])])
if j == len(a) - 1:
p3 = (a[0], b[0])
else:
p3 = (a[j + 1], b[j + 1])
rp3 = (x0[x0.index(p3[0])], y0[x0.index(p3[0])])
if p1 == p2 or p1 == p3:
pass
else:
# print('====+++===', p1, p2, p3)
# print('====+++===', rp1, rp2, rp3)
# 铅垂线长度
s_length = abs(
(p2[1] * (p1[0] - p3[0])) + (p1[1] * (p3[0] - p2[0]) + (p3[1] * (p2[0] - p1[0]))) / (p3[0] - p2[0]))
s = find_cross(p1, p2, p3)
if isPointinPolygon(s, point_lst):
dic[s_length] = [rp1, rp2, rp3]
print('铅垂线长度', s_length, s, dic[s_length])
# print(dic)
max_s_length = dic[max(dic)]
rp1 = max_s_length[0]
rp2 = max_s_length[1]
rp3 = max_s_length[2]
# 点和直线两端点中点连线的斜率
k = (((rp1[1] + rp2[1]) / 2) - ((rp1[1] + rp3[1]) / 2)) / (((rp1[0] + rp2[0]) / 2) - ((rp1[0] + rp3[0]) / 2))
b = ((rp1[1] + rp2[1]) / 2) - (k * ((rp1[0] + rp2[0]) / 2))
return k, b
def isPointinPolygon(point, rangelist): # [[0,0],[1,1],[0,1],[0,0]] [1,0.8]
print(point)
# 判断是否在外包矩形内，如果不在，直接返回false
lnglist = []
latlist = []
for i in range(len(rangelist) - 1):
lnglist.append(rangelist[i][0])
latlist.append(rangelist[i][1])
# print(lnglist, latlist)
maxlng = max(lnglist)
minlng = min(lnglist)
maxlat = max(latlist)
minlat = min(latlist)
# print(maxlng, minlng, maxlat, minlat)
if (point[0] > maxlng or point[0] < minlng or
point[1] > maxlat or point[1] < minlat):
return False
count = 0
point1 = rangelist[0]
for i in range(1, len(rangelist)):
point2 = rangelist[i]
# 点与多边形顶点重合
if (point[0] == point1[0] and point[1] == point1[1]) or (point[0] == point2[0] and point[1] == point2[1]):
print("在顶点上")
return False
# 判断线段两端点是否在射线两侧 不在肯定不相交 射线（-∞，lat）（lng,lat）
if (point1[1] < point[1] and point2[1] >= point[1]) or (point1[1] >= point[1] and point2[1] < point[1]):
# 求线段与射线交点 再和lat比较
point12lng = point2[0] - (point2[1] - point[1]) * (point2[0] - point1[0]) / (point2[1] - point1[1])
# print(point12lng)
# 点在多边形边上
if (point12lng == point[0]):
print("点在多边形边上")
return True
if (point12lng < point[0]):
count += 1
point1 = point2
print(count)
if count % 2 == 0:
return True
else:
return True
if __name__ == '__main__':
x0 = [1.00, 2.00, 3.00, 4.00, 5.00, 6.00]
y0 = [2.02, 4.00, 5.98, 7.90, 10.10, 12.05]
dy = find_line(x0, y0)
print(dy)
a0, b0 = find_poly(x0, dy)
print(a0, b0)
print(find_perfect_line(a0, b0, x0, y0))
