質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

89.10%

Scipyのsplprepによるスプライン補間について

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 177

matplotlibのキャンバス上に手入力で制御点を打ち,その点をスプライン補間しています。スプライン曲線の算出にはScipyのsplprepを用いています。しかし打ち込む点が多くなると,曲線が滑らかさを失い、がたつきがある線になってしまいます。原因はなんでしょうか。
また、算出した曲線の曲率を算出したいのですが、近接した3点を取る方法がわかりません。ご教授願います。よろしくお願いします。

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate


def main():
    fig = plt.figure(figsize=(11.0, 5.0))
    ClickAddPoints(fig)


class ClickAddPoints:
    def __init__(self, fig):
        self.fig = fig
        self.ax1 = fig.add_subplot(121)
        # ax2に曲率を表示
        self.ax2 = fig.add_subplot(122)
        self.ax1.figure.canvas.mpl_connect('button_press_event', self.on_click)
        self.ax1.figure.canvas.mpl_connect('key_press_event', self.on_key)
        self.x = []
        self.y = []
        self.r = []
        self.s = []
        self.x2 = []
        self.y2 = []
        self.update_plot()
        self.div = 1000  # スプラインの分割数。多いほどなめらかに描画。
        plt.show()

    def calc_spline(self, x, y, point, deg):
        tck, u = interpolate.splprep([x, y], k=deg, s=0)
        u = np.linspace(0, 1, num=point, endpoint=True)
        spline = interpolate.splev(u, tck, der=0)
        return spline[0], spline[1]

    def cal_curveture(self, x1, y1, x2, y2, x3, y3):
        r = 1/(abs(x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3)/(
            ((x1-x2)**2+(y1-y2)**2)*((x2-x3)**2+(y2-y3)**2)*((x3-x1)**2+(y3-y1)**2))**(1/2))
        return r

    def cal_curve_length(self, x1, y1, x2, y2):
        s = np.sqrt((x2-x1)**2+(y2-y1)**2)
        return s

    def cal_clothoid(self, T, dt, V, R, _P0x, _P0y, _dθ):
        #クロソイドの計算
        v = (V*1000)/3600
        L = v*T
        A = np.sqrt(L*R)
        dL = v*dt
        P0x = _P0x
        P0y = _P0y
        Lt = 0
        dθ = _dθ
        dLt = dL
        A2 = A**2
        tsteps = np.arange(0, T, dt)

        for i in tsteps:
            x2 = P0x
            y2 = P0y

            self.x2.append(x2)
            self.y2.append(y2)

            Lt = Lt + dLt
            Rt = A2 / Lt
            k = 1 / Rt
            dθ = dθ + dLt / Rt
            dx = dL * np.cos(dθ)
            dy = dL * np.sin(dθ)
            P1x = P0x + dx
            P1y = P0y + dy
            P0x = P1x
            P0y = P1y

        return self.x2, self.y2

    def update_plot(self):
        self.ax1.set_xticks(np.linspace(0, 10, 5))
        self.ax1.set_yticks(np.linspace(0, 10, 5))
        self.ax1.set_aspect('equal')
        self.ax1.figure.canvas.draw()
        self.ax2.set_ylim(0, 40)
        self.ax2.figure.canvas.draw()

    def on_key(self, event):
        # qキーで終了
        if event.key == 'q':
            print('Finish')
            return

        # wキーで全削除
        if event.key == 'w':
            self.x.clear()
            self.y.clear()
            self.redraw()
            print('All Clear')

        # pキーでスプライン
        if event.key == 'p':
            print('spline')

        # cキーでクロソイド
        elif event.key == 'c':
            print('clothoid')

        # tキーで直線
        elif event.key == 't':
            print('straight line')

        else:
            return

    def on_click(self, event):
        if event.button == 1:
            # コントロール点を追加
            if event.inaxes is not self.ax1:
                return
            # 過去の座標と同一座標をクリックした場合は抜ける
            if event.xdata in self.x and event.ydata in self.y:
                return
            self.x.append(event.xdata)
            self.y.append(event.ydata)
            self.redraw()
            print('Added   no.{} point at [{} {}]'.format(len(self.x), self.x[-1], self.y[-1]))
        elif event.button == 3:
            # 直近のコントロール点を削除
            if len(self.x) > 0:
                self.x.pop()
                self.y.pop()
                self.redraw()
                print('Removed no.{0} point'.format(len(self.x)+1))

    def redraw(self):
        self.ax1.cla()
        self.ax2.cla()
        count = len(self.x)
        if count > 0:
            # クリック点の描画
            self.ax1.plot(self.x, self.y, 'ro')
        if count <= 1:
            self.update_plot()
            return
        if count == 2:
            deg = 1
        elif count == 3:
            deg = 2
        elif count > 3:
            deg = 3
        else:
            return

        x1, y1 = self.calc_spline(self.x, self.y, self.div, deg)
        #曲線長
        if len(self.x) >= 2:
            self.s = [self.cal_curve_length(x1[i],y1[i],
                                            x1[i+1],y1[i+1])
                      for i in range(0, len(x1)-1)]

        # 曲率データの更新
        if len(self.x) > 2:
            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)]

        # 線の描画
        self.ax1.plot(x1, y1)
        self.ax2.plot(self.r)
        self.update_plot()


main()
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 1

checkベストアンサー

+1

曲線が滑らかさを失い、がたつきがある線になってしまいます。原因はなんでしょうか。

np.linspace(0, 1, num, endpoint)
のうち、numパラメータは曲線の分割数に対応しています。(曲線をnum-1個の直線に分割)
numが大きければ大きいほど、なめらかになります。
たとえばnum=3とすると、最初からガタガタの線になります。
質問元のコードだと300なので、100個以上コントロールポイントを追加したあたりからガタつきが目立ってくるのかもしれません。
下記コードでは、__init__内で、div属性を分割数として設定しています。(1000個に分割)

近接した3点を取る方法がわかりません。

 x1, y1 = self.calc_spline(self.x, self.y, self.div, deg)
のx1、y1はそれぞれ曲線を構成するdiv個に分割された点のx座標配列、y座標配列となっているので、
この配列から3点を取り出し、曲率を計算しているcal_curvetureに渡せばOKです。
コード中の

        # 曲率データの更新
        if len(self.x) > 2:
            self.k = [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)]


の部分です。

全体コードは以下:

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate


def main():
    fig = plt.figure(figsize=(11.0, 5.0))
    ClickAddPoints(fig)


class ClickAddPoints:
    def __init__(self, fig):
        self.fig = fig
        self.ax1 = fig.add_subplot(121)
        # ax2に曲率を表示
        self.ax2 = fig.add_subplot(122)
        self.ax1.figure.canvas.mpl_connect('button_press_event', self.on_click)
        self.ax1.figure.canvas.mpl_connect('key_press_event', self.on_key)
        self.x = []
        self.y = []
        self.k = []
        self.update_plot()
        self.div = 1000  # スプラインの分割数。多いほどなめらかに描画。
        plt.show()

    def calc_spline(self, x, y, point, deg):
        tck, u = interpolate.splprep([x, y], k=deg, s=0)
        u = np.linspace(0, 1, num=point, endpoint=True)
        spline = interpolate.splev(u, tck, der=0)
        return spline[0], spline[1]

    def cal_curveture(self, x1, y1, x2, y2, x3, y3):
        k = abs(x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3)/(
            ((x1-x2)**2+(y1-y2)**2)*((x2-x3)**2+(y2-y3)**2)*((x3-x1)**2+(y3-y1)**2))**(1/2)
        return k

    def update_plot(self):
        self.ax1.set_xticks(np.linspace(0, 6, 6))
        self.ax1.set_yticks(np.linspace(0, 6, 6))
        self.ax1.set_aspect('equal')
        self.ax1.figure.canvas.draw()
        self.ax2.figure.canvas.draw()

    def on_click(self, event):
        if event.button == 1:
            # コントロール点を追加
            if event.inaxes is not self.ax1:
                return
            # 過去の座標と同一座標をクリックした場合は抜ける
            if event.xdata in self.x and event.ydata in self.y:
                return
            self.x.append(event.xdata)
            self.y.append(event.ydata)
            self.redraw()
            print('Added   no.{} point at [{} {}]'.format(len(self.x), self.x[-1], self.y[-1]))
        elif event.button == 3:
            # 直近のコントロール点を削除
            if len(self.x) > 0:
                self.x.pop()
                self.y.pop()
                self.redraw()
                print('Removed no.{0} point'.format(len(self.x)+1))

    def redraw(self):
        self.ax1.cla()
        self.ax2.cla()
        count = len(self.x)
        if count > 0:
            # クリック点の描画
            self.ax1.plot(self.x, self.y, 'ro')
        if count <= 1:
            self.update_plot()
            return
        if count == 2:
            deg = 1
        elif count == 3:
            deg = 2
        elif count > 3:
            deg = 3
        else:
            return

        x1, y1 = self.calc_spline(self.x, self.y, self.div, deg)
        # 曲率データの更新
        if len(self.x) > 2:
            self.k = [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)]
        # 線の描画
        self.ax1.plot(x1, y1)
        self.ax2.plot(self.k)
        self.update_plot()

    def on_key(self, event):
        # qキーで終了
        if event.key == 'q':
            print('Finish')
            return

        # wキーで全削除
        if event.key == 'w':
            self.x.clear()
            self.y.clear()
            self.redraw()
            print('All Clear')


main()

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/07/04 13:08

    ありがとうございます。
    最後に一つよろしいでしょうか。
    モード切替をした際に,それまでに描画した線はそのまま残したいのですが,どのようにしたらできますでしょうか。たとえばスプライン曲線を何点か描画した後,直線モードに移行すると,それまでのスプライン曲線が直線になってしまうのですが,スプライン曲線はスプライン曲線で残したいのです。
    よろしくお願いします。

    キャンセル

  • 2020/07/04 19:30

    ・tキーによる直線重ね合わせ表示/非表示を切り替えられるようにしました。
    ・直線の場合に右側プロットで曲率がゼロ以外になるのを修正
    ・ RuntimeWarning: invalid value encountered in double_scalars エラーを修正しました。
    https://gist.github.com/taizan-hokuto/6f53da1338d624a50992ed68673fb914

    キャンセル

  • 2020/07/04 21:12

    丁寧な回答からその後の要望への対応もありがとうございました。
    とても勉強になりました。

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 89.10%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる