溫馨提示×
您好，登錄后才能下訂單哦！
點擊 登錄注冊 即表示同意《億速云用戶服務條款》
您好，登錄后才能下訂單哦！
前言
很多文章在談及曲線平滑的時候，習慣使用擬合的概念，我認為這是不恰當?shù)摹Ｆ交蟮那€，一定經(jīng)過原始的數(shù)據(jù)點，而擬合曲線，則不一定要經(jīng)過原始數(shù)據(jù)點。
一般而言，需要平滑的數(shù)據(jù)分為兩種：時間序列的單值數(shù)據(jù)、時間序列的二維數(shù)據(jù)。對于前者，并非一定要用貝塞爾算法，僅用樣條插值就可以輕松實現(xiàn)平滑；而對于后者，不管是 numpy 還是 scipy 提供的那些插值算法，就都不適用了。
本文基于三階貝塞爾曲線，實現(xiàn)了時間序列的單值數(shù)據(jù)和時間序列的二維數(shù)據(jù)的平滑算法，可滿足大多數(shù)的平滑需求。
貝塞爾曲線
關于貝塞爾曲線的數(shù)學原理，這里就不討論了，直接貼出結論：
一階貝塞爾曲線
二階貝塞爾曲線
三階貝塞爾曲線
算法描述
如果我們把三階貝塞爾曲線的 P0 和 P3 視為原始數(shù)據(jù)，只要找到 P1 和 P2 兩個點（我們稱其為控制點），就可以根據(jù)三階貝塞爾曲線公式，計算出 P0 和 P3 之間平滑曲線上的任意點。
現(xiàn)在，平滑問題變成了如何計算兩個原始數(shù)據(jù)點之間的控制點的問題。步驟如下：
第1步：綠色直線連接相鄰的原始數(shù)據(jù)點，計算出個線段的中點，紅色直線連接相鄰的中點
第2步：根據(jù)相鄰兩條綠色直線長度之比，分割其中點之間紅色連線，標記分割點
第3步：平移紅色連線，使其分割點與相對的原始數(shù)據(jù)點重合
第4步：調(diào)整平移后紅色連線的端點與原始數(shù)據(jù)點的距離，通常縮減40%-80%
算法實現(xiàn)
# -*- coding: utf-8 -*-
import numpy as np
def bezier_curve(p0, p1, p2, p3, inserted):
"""
三階貝塞爾曲線
p0, p1, p2, p3 - 點坐標，tuple、list或numpy.ndarray類型
inserted - p0和p3之間插值的數(shù)量
"""
assert isinstance(p0, (tuple, list, np.ndarray)), u'點坐標不是期望的元組、列表或numpy數(shù)組類型'
assert isinstance(p0, (tuple, list, np.ndarray)), u'點坐標不是期望的元組、列表或numpy數(shù)組類型'
assert isinstance(p0, (tuple, list, np.ndarray)), u'點坐標不是期望的元組、列表或numpy數(shù)組類型'
assert isinstance(p0, (tuple, list, np.ndarray)), u'點坐標不是期望的元組、列表或numpy數(shù)組類型'
if isinstance(p0, (tuple, list)):
p0 = np.array(p0)
if isinstance(p1, (tuple, list)):
p1 = np.array(p1)
if isinstance(p2, (tuple, list)):
p2 = np.array(p2)
if isinstance(p3, (tuple, list)):
p3 = np.array(p3)
points = list()
for t in np.linspace(0, 1, inserted+2):
points.append(p0*np.power((1-t),3) + 3*p1*t*np.power((1-t),2) + 3*p2*(1-t)*np.power(t,2) + p3*np.power(t,3))
return np.vstack(points)
def smoothing_base_bezier(date_x, date_y, k=0.5, inserted=10, closed=False):
"""
基于三階貝塞爾曲線的數(shù)據(jù)平滑算法
date_x - x維度數(shù)據(jù)集，list或numpy.ndarray類型
date_y - y維度數(shù)據(jù)集，list或numpy.ndarray類型
k - 調(diào)整平滑曲線形狀的因子，取值一般在0.2~0.6之間。默認值為0.5
inserted - 兩個原始數(shù)據(jù)點之間插值的數(shù)量。默認值為10
closed - 曲線是否封閉，如是，則首尾相連。默認曲線不封閉
"""
assert isinstance(date_x, (list, np.ndarray)), u'x數(shù)據(jù)集不是期望的列表或numpy數(shù)組類型'
assert isinstance(date_y, (list, np.ndarray)), u'y數(shù)據(jù)集不是期望的列表或numpy數(shù)組類型'
if isinstance(date_x, list) and isinstance(date_y, list):
assert len(date_x)==len(date_y), u'x數(shù)據(jù)集和y數(shù)據(jù)集長度不匹配'
date_x = np.array(date_x)
date_y = np.array(date_y)
elif isinstance(date_x, np.ndarray) and isinstance(date_y, np.ndarray):
assert date_x.shape==date_y.shape, u'x數(shù)據(jù)集和y數(shù)據(jù)集長度不匹配'
else:
raise Exception(u'x數(shù)據(jù)集或y數(shù)據(jù)集類型錯誤')
# 第1步：生成原始數(shù)據(jù)折線中點集
mid_points = list()
for i in range(1, date_x.shape[0]):
mid_points.append({
'start': (date_x[i-1], date_y[i-1]),
'end': (date_x[i], date_y[i]),
'mid': ((date_x[i]+date_x[i-1])/2.0, (date_y[i]+date_y[i-1])/2.0)
})
if closed:
mid_points.append({
'start': (date_x[-1], date_y[-1]),
'end': (date_x[0], date_y[0]),
'mid': ((date_x[0]+date_x[-1])/2.0, (date_y[0]+date_y[-1])/2.0)
})
# 第2步：找出中點連線及其分割點
split_points = list()
for i in range(len(mid_points)):
if i < (len(mid_points)-1):
j = i+1
elif closed:
j = 0
else:
continue
x00, y00 = mid_points[i]['start']
x01, y01 = mid_points[i]['end']
x10, y10 = mid_points[j]['start']
x11, y11 = mid_points[j]['end']
d0 = np.sqrt(np.power((x00-x01), 2) + np.power((y00-y01), 2))
d1 = np.sqrt(np.power((x10-x11), 2) + np.power((y10-y11), 2))
k_split = 1.0*d0/(d0+d1)
mx0, my0 = mid_points[i]['mid']
mx1, my1 = mid_points[j]['mid']
split_points.append({
'start': (mx0, my0),
'end': (mx1, my1),
'split': (mx0+(mx1-mx0)*k_split, my0+(my1-my0)*k_split)
})
# 第3步：平移中點連線，調(diào)整端點，生成控制點
crt_points = list()
for i in range(len(split_points)):
vx, vy = mid_points[i]['end'] # 當前頂點的坐標
dx = vx - split_points[i]['split'][0] # 平移線段x偏移量
dy = vy - split_points[i]['split'][1] # 平移線段y偏移量
sx, sy = split_points[i]['start'][0]+dx, split_points[i]['start'][1]+dy # 平移后線段起點坐標
ex, ey = split_points[i]['end'][0]+dx, split_points[i]['end'][1]+dy # 平移后線段終點坐標
cp0 = sx+(vx-sx)*k, sy+(vy-sy)*k # 控制點坐標
cp1 = ex+(vx-ex)*k, ey+(vy-ey)*k # 控制點坐標
if crt_points:
crt_points[-1].insert(2, cp0)
else:
crt_points.append([mid_points[0]['start'], cp0, mid_points[0]['end']])
if closed:
if i < (len(mid_points)-1):
crt_points.append([mid_points[i+1]['start'], cp1, mid_points[i+1]['end']])
else:
crt_points[0].insert(1, cp1)
else:
if i < (len(mid_points)-2):
crt_points.append([mid_points[i+1]['start'], cp1, mid_points[i+1]['end']])
else:
crt_points.append([mid_points[i+1]['start'], cp1, mid_points[i+1]['end'], mid_points[i+1]['end']])
crt_points[0].insert(1, mid_points[0]['start'])
# 第4步：應用貝塞爾曲線方程插值
out = list()
for item in crt_points:
group = bezier_curve(item[0], item[1], item[2], item[3], inserted)
out.append(group[:-1])
out.append(group[-1:])
out = np.vstack(out)
return out.T[0], out.T[1]
if __name__ == '__main__':
import matplotlib.pyplot as plt
x = np.array([2,4,4,3,2])
y = np.array([2,2,4,3,4])
plt.plot(x, y, 'ro')
x_curve, y_curve = smoothing_base_bezier(x, y, k=0.3, closed=True)
plt.plot(x_curve, y_curve, label='$k=0.3$')
x_curve, y_curve = smoothing_base_bezier(x, y, k=0.4, closed=True)
plt.plot(x_curve, y_curve, label='$k=0.4$')
x_curve, y_curve = smoothing_base_bezier(x, y, k=0.5, closed=True)
plt.plot(x_curve, y_curve, label='$k=0.5$')
x_curve, y_curve = smoothing_base_bezier(x, y, k=0.6, closed=True)
plt.plot(x_curve, y_curve, label='$k=0.6$')
plt.legend(loc='best')
plt.show()
下圖為平滑效果。左側是封閉曲線，兩個原始數(shù)據(jù)點之間插值數(shù)量為默認值10；右側為同樣數(shù)據(jù)不封閉的效果，k值默認0.5.
參考資料
算法參考了 Interpolation with Bezier Curves 這個網(wǎng)頁，里面沒有關于作者的任何信息，在此只能籠統(tǒng)地向國際友人表示感謝！
以上就是本文的全部內(nèi)容，希望對大家的學習有所幫助，也希望大家多多支持億速云。
免責聲明：本站發(fā)布的內(nèi)容（圖片、視頻和文字）以原創(chuàng)、轉載和分享為主，文章觀點不代表本網(wǎng)站立場，如果涉及侵權請聯(lián)系站長郵箱：is@yisu.com進行舉報，并提供相關證據(jù)，一經(jīng)查實，將立刻刪除涉嫌侵權內(nèi)容。
