mieki256's diary

一般的に、scipy を使ってスプライン補間をする場合は、interpolate.interp1d() を使うようだけど、x の増分方向が同じ方向じゃないと使えないらしくて悩んでしまった。自分が処理したいデータは、平面上であっちに行ったりこっちに行ったりするデータなので…。

以下のページで、interpolate.splprep() ならそういうデータも処理できそうと知った。また、interpolate.interp1d() も、x と y を別々に処理すればどうにかできるらしい。

_pythonでxy座標上の離散点をスプライン補間 #Python - Qiita
_スプライン曲線と補間 - Emotion Explorer

上記のページのサンプルを参考にして、以下のデータを処理してみた。

_housakatouge.csv
_motosuko.csv

Pythonスクリプトは以下。

_draw_spline.py

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


# S = 10
S = 20


def loadData(infile):
    ax = []
    ay = []
    with open(infile, "r") as file:
        for line in file:
            line = line.rstrip()
            x, y = line.split(",")
            x = float(x)
            y = float(y)
            ax.append(x)
            ay.append(y)
    return ax, ay


# B-Spline 補間
def spline1(x, y):
    tck, u = interpolate.splprep([x, y], k=3, s=0)
    u = np.linspace(0, 1, num=len(x) * S, endpoint=True)
    spline = interpolate.splev(u, tck)
    return spline[0], spline[1]


# interp1dスプライン補間による曲線
def spline3(x, y):
    t = np.linspace(0, 1, len(x), endpoint=True)
    fx = interpolate.interp1d(t, x, kind="cubic")
    fy = interpolate.interp1d(t, y, kind="cubic")
    ta = np.linspace(0, 1, len(x) * S, endpoint=True)
    return fx(ta), fy(ta)


if len(sys.argv) < 2:
    sys.exit()

infile = sys.argv[1]
x, y = loadData(infile)

# B-Spline interpolation
bx, by = spline1(x, y)

# interp1d spline
ax, ay = spline3(x, y)

plt.plot(x, y, "ro", label="Data")
# plt.plot(bx, by, "b.", label="B-Spline interpolation")
plt.plot(bx, by, color="blue", label="B-Spline interpolation")
plt.plot(ax, ay, color="green", label="interp1d")

plt.legend(loc="best")
plt.savefig("output_graph.png")
plt.show()

python draw_spline.py motosuko.csv

なかなか興味深い(?)結果になった気がする。

まず、連続した直線データで構成されている道路データを読み込んで、スプライン曲線として描画することはできている。

ただ、interpolate.interp1d() で処理すると、不自然な曲線になってしまう箇所が、いくつも出てくるようだなと…。これが interpolate.splprep() なら、そこまでおかしな状態にはならない。これは後者を使って処理したほうが良さそうだなと…。

余談。matplotlib を使ってグラフを描画しているけれど、ツールバーに拡大ツールやPANツールがあって、これは便利だなと感心してしまった。

• 虫眼鏡のアイコンをクリックして、キャンバス(?)上でマウスドラッグをして範囲を指定すると、その範囲が拡大表示される。
• 上下左右の矢印っぽいアイコンをクリックすると、キャンバス上のマウスドラッグで表示位置を変更できる。
• 家のアイコンをクリックすると、全体が収まる表示になる。

スプライン曲線を描画することはできたけど、座標値として出力するにはどうしたらいいのか…。

以下のような感じだろうか。

_dump_spline.py

import argparse
import sys
import numpy as np
from scipy import interpolate


# load csv file
def loadData(infile):
    ax = []
    ay = []
    with open(infile, "r") as file:
        for line in file:
            line = line.rstrip()
            x, y = line.split(",")
            x = float(x)
            y = float(y)
            ax.append(x)
            ay.append(y)
    return ax, ay


# B-Spline
def spline1(x, y, s):
    tck, u = interpolate.splprep([x, y], k=3, s=0)
    u = np.linspace(0, 1, num=len(x) * s, endpoint=True)
    spline = interpolate.splev(u, tck)
    return spline[0], spline[1]


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("-i", "--input", required=True, help="Input csv file")
    parser.add_argument("-n", "--num", type=int, default=10, help="number")
    args = parser.parse_args()

    if args.input is None:
        sys.exit()

    x, y = loadData(args.input)
    bx, by = spline1(x, y, args.num)

    for i in range(len(bx)):
        print(f"{bx[i]},{by[i]}")


if __name__ == "__main__":
    main()

python -i dump_spline.py -n 10 motosuko.csv > motosuko_spline.csv

これで、スプライン曲線化したデータを取得することができた。

_housakatouge_spline.csv
_motosuko_spline.csv

しかし、このデータは、本当にスプライン曲線になっているのだろうか…？ 表示して確認してみないと…。

先日は Python + pygame でデータを描画して確認したけれど、matplotlib を使って描画したほうが、拡大表示や表示位置変更ができて便利だろうと思えてきた。

以下のような感じだろうか。

_draw_line_plot.py

import sys
import matplotlib.pyplot as plt


def loadData(infile):
    ax = []
    ay = []
    with open(infile, "r") as file:
        for line in file:
            line = line.rstrip()
            x, y = line.split(",")
            x = float(x)
            y = float(y)
            ax.append(x)
            ay.append(y)
    return ax, ay


if len(sys.argv) < 2:
    sys.exit()

infile = sys.argv[1]
x, y = loadData(infile)

plt.plot(x, y, color="cyan", label="Data (line)")
plt.plot(x, y, "r.", label="Data (point)")

plt.legend(loc="best")
# plt.savefig("output_graph.png")
plt.show()

python draw_line_plot.py motosuko_spline.csv

イイ感じかもしれない。

このデータを元にして、道路幅を加えて、ポリゴン描画すればそれらしくなるのではないかな…。