投稿2020/07/27 06:37
クリックで点を追加し、その点を3次スプライン補間するプログラムです。しかしスタートとゴールが直線で結ばれてしまいます。原因がわかりません。どなたかご指摘のほどよろしくお願いします。
Python
1import matplotlib.pyplot as plt 2import numpy as np 3from scipy import interpolate 4 5 6def main(): 7 fig = plt.figure(figsize=(11.0, 5.0)) 8 ClickAddPoints(fig) 9 10class ClickAddPoints: 11 def __init__(self, fig): 12 self.fig = fig 13 self.ax1 = fig.add_subplot(121) 14 # ax2に曲率を表示 15 self.ax2 = fig.add_subplot(122) 16 self.ax1.figure.canvas.mpl_connect('button_press_event', self.on_click) 17 self.ax1.figure.canvas.mpl_connect('key_press_event', self.on_key) 18 self.x = [] 19 self.y = [] 20 self.r = [] 21 self.update_plot() 22 plt.show() 23 24 def calculateSplineParameter(self, x): 25 # t = 0,1,2,...,x.size-1とする。 26 x = np.array(x) 27 n = len(x) - 1 28 h = np.ones(n) 29 a = x 30 31 # 対角成分 32 Adiag = np.eye(n + 1) 33 Adiag[1:-1, 1:-1] *= 2 * (h[:-1] + h[1:]) 34 35 # 対角成分の一つ下の斜めの成分 36 Alower = np.eye(n + 1, k=-1) 37 Alower[1:, :-1] *= h 38 Alower[-1, -2] = 0 39 40 # 対角成分の一つ上の斜めの成分 41 Aupper = np.eye(n + 1, k=1) 42 Aupper[:-1, 1:] *= h 43 Aupper[0, 1] = 0 44 45 A = Adiag + Alower + Aupper 46 47 # Aの逆行列 48 Ainv = np.linalg.inv(A) 49 50 r = np.zeros(n + 1) 51 r[1:-1] = (3 / h)[1:] * (a[2:] - a[1:-1]) - (3 / h)[:-1] * (a[1:-1] - a[:-2]) 52 53 c = np.dot(Ainv, r) 54 d = (c[1:] - c[:-1]) / (3 * h) 55 b = (a[1:] - a[:-1]) / h - h / 3 * (2 * c[:-1] + c[1:]) 56 57 return [a, b, c, d] 58 59 # 三次関数の係数から, 元の点から間の点を求める。 60 def calculateInterpolationPoints(self, splineParameter, orignalPointsNumber, interpolationPointsNumber=100): 61 a, b, c, d = splineParameter 62 63 index = np.repeat(np.arange(orignalPointsNumber - 1), interpolationPointsNumber, axis=0) 64 t = np.tile(np.linspace(0, 1, interpolationPointsNumber), orignalPointsNumber - 1) 65 66 return a[index] + b[index] * t + c[index] * t ** 2 + d[index] * t ** 3 67 68 def cal_curveture(self, x1, y1, x2, y2, x3, y3): 69 r = 1 70 return r 71 72 def update_plot(self): 73 self.ax1.set_xticks(np.linspace(0, 6, 6)) 74 self.ax1.set_yticks(np.linspace(0, 6, 6)) 75 self.ax1.set_aspect('equal') 76 self.ax1.figure.canvas.draw() 77 self.ax2.set_ylim(0, 40) 78 self.ax2.figure.canvas.draw() 79 80 def on_click(self, event): 81 if event.button == 1: 82 # コントロール点を追加 83 if event.inaxes is not self.ax1: 84 return 85 # 過去の座標と同一座標をクリックした場合は抜ける 86 if event.xdata in self.x and event.ydata in self.y: 87 return 88 self.x.append(event.xdata) 89 self.y.append(event.ydata) 90 self.redraw() 91 print('Added no.{} point at [{} {}]'.format( 92 len(self.x), self.x[-1], self.y[-1])) 93 elif event.button == 3: 94 # 直近のコントロール点を削除 95 if len(self.x) > 0: 96 self.x.pop() 97 self.y.pop() 98 self.redraw() 99 print('Removed no.{0} point'.format(len(self.x)+1)) 100 101 def redraw(self): 102 self.ax1.cla() 103 self.ax2.cla() 104 count = len(self.x) 105 if 0 < count < 2: 106 # クリック点の描画 107 self.ax1.plot(self.x, self.y, 'ro') 108 109 elif count >= 2: 110 x1 = [] 111 y1 = [] 112 113 for i in range(1, len(self.x)-1, 1): 114 parameterX = self.calculateSplineParameter(self.x) 115 parameterY = self.calculateSplineParameter(self.y) 116 117 X = self.calculateInterpolationPoints(parameterX, len(self.x)) 118 Y = self.calculateInterpolationPoints(parameterY, len(self.y)) 119 120 x1.extend(X) 121 y1.extend(Y) 122 # 曲率データの更新 123 self.r = [self.cal_curveture(x1[i - 1], y1[i - 1], 124 x1[i], y1[i], 125 x1[i + 1], y1[i + 1]) 126 for i in range(1, len(x1) - 1)] 127 self.ax1.plot(self.x, self.y, 'ro') 128 self.ax1.plot(x1, y1) 129 130 else: 131 return 132 133 # 線の描画 134 135 self.ax2.plot(self.r) 136 self.update_plot() 137 138 def on_key(self, event): 139 # qキーで終了 140 if event.key == 'q': 141 print('Finish') 142 return 143 144 # wキーで全削除 145 if event.key == 'w': 146 self.x.clear() 147 self.y.clear() 148 self.x1.clear() 149 self.y1.clear() 150 self.redraw() 151 print('All Clear') 152 153 154main()
回答1件
0
自己解決
parameterX = self.calculateSplineParameter(self.x)
parameterY = self.calculateSplineParameter(self.y)
X = self.calculateInterpolationPoints(parameterX, len(self.x))
Y = self.calculateInterpolationPoints(parameterY, len(self.y))
をfor文の中に入れていたため
投稿2020/07/27 07:27
総合スコア18
あなたの回答
tips
プレビュー
Q&A
