スタートとゴールを結ぶ直線をなくしたい。

クリックで点を追加し、その点を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()

Penpen7

2020/07/27 07:03 編集

なぜ 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文の中に入れたのですか？入れても意味ありませんよ。

takes.it.easy

2020/07/27 07:05

以前使っていたほかのプログラムに新しくスプライン補間の式を移植したものなので、考え無しにfor文の中に入れてしまいました。ありがとうございます。

Penpen7

2020/07/27 07:12

直線で繋がるのは x1.extend(X), y1.extend(Y) をfor文でなんども実行しているからです。

takes.it.easy

2020/07/27 07:28

理解しました。ありがとうございます。

行動規範の内容に同意します

回答1件

自己解決

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

takes.it.easy

総合スコア19

Penpen7

2020/07/27 08:05

それは改善すべきところですが、その部分は質問の答えとしては本質的ではなくて、 x1.extend(X), y1.extend(Y) をなんども実行しているのが原因です。 parameterX = self.calculateSplineParameter(self.x) parameterY = self.calculateSplineParameter(self.y) X = self.calculateInterpolationPoints(parameterX, len(self.x)) Y = self.calculateInterpolationPoints(parameterY, len(self.y)) x1.extend(X) y1.extend(Y) はfor文に含めずに # 曲率データの更新 self.r = [self.cal_curveture(x1[i - 1], y1[i - 1], x1[i], y1[i], x1[i + 1], y1[i + 1]) for i in range(1, len(x1) - 1)] の部分をfor文に入れないといけません。というのが解答としては正しいです。

takes.it.easy

2020/07/27 09:41

なるほど、本質を理解していませんでした。 x1.extend(X)とy1.extend(Y)はx1とy1に要素を追加するメソッドであるため，その部分を繰り返してしまうと，同じ要素が繰り返し追加されることで，直線が描かれていた。さらに曲率は繰り返し計算する必要があるため，for文に含めなければならない。という理解でよろしいでしょうか。

Penpen7

2020/07/27 11:21 編集

XとYの中身を見るとそれが分かると思います。私のコードでは、calculateInterpolationPointsを一回呼び出せばスプライン曲線のすべての点を出力するようにしたので、何度もfor文で計算する必要がないからです。 interpolationPointsNumber=100を=3くらいに少なくしてみてみると原因はつかめるのではないでしょうか。同じものが周期的に表れているはずです。 ># 曲率データの更新 >self.r = [self.cal_curveture(x1[i - 1], y1[i - 1], >x1[i], y1[i], >x1[i + 1], y1[i + 1]) >for i in range(1, len(x1) - 1)] >の部分をfor文に入れないといけません。というのが解答としては正しいです。といいましたが、たぶん曲率の計算もforはいらないですね。結局のところ、for文を消せば解決します。

行動規範の内容に同意します