この記事のポイント
- ポリラインのスムージング方法として、3次Bスプライン、Rbf、区分的3次エルミート内挿(pchip)、ベジェ内挿、Chaikin's algorithmsの5手法の比較を行いました。
- 滑らかさの点でみるとchaikinが良さそうです。
- 元の点を通る必要がある場合はpchip、bezierのいずれかを使う必要があります。
- Bsplineが最も汎用的ですがパラメータの調整が必要です。
はじめに
- 疎な点により形成されるポリラインオブジェクトから滑らかなポリラインオブジェクトをつくります。
- スムージング方法は次の2つに分類されます。
- 補間による方法：点間を滑らかなに接続する。
- 近似による方法：点情報から近似線を作成する。元の点は必ずしも通らない。
- また、xyグラフでxが昇順に並ぶデータを対象とする数学的な方法、ジオメトリデータのように複雑なデータを対象とする幾何学的な方法がある。この点が、ジオメトリデータの取り扱いの難しいところだと思います。
- 言語はすべてpythonで書いてます。また、本ページ内の全グラフはいつもどおりHoloviewsを使用しています。
computational-sediment-hyd.hatenablog.jp
サンプルデータ
import numpy as np import holoviews as hv hv.extension('bokeh') points = np.array([[14747743.41590873, 4097001.07511239], [14747765.19891978, 4097203.17304826], [14747807.5547746 , 4097423.42349333], [14747841.43945846, 4097607.36891999], [14747853.54113126, 4097671.50778586], [14747857.17163311, 4097778.00250656], [14747849.91062942, 4097844.56170699], [14747834.17845477, 4097902.64973646], [14747733.73457048, 4097953.47676225], [14747655.07369724, 4097963.15810049], [14747544.9484747 , 4097935.32425304], [14747465.07743418, 4097888.1277291 ], [14747357.37254621, 4097788.89401208], [14747261.76933104, 4097689.66029507]]) g0 = hv.Scatter(points, label='Original').options(show_grid=True,size=7,width=400,height=400,legend_position='bottom_left') # hv.save(g0, 'forg.html')
数学的手法
前処理：媒介変数表示
前述のとおり、ジオメトリデータは昇順ではないため、数学的手法を適用するためには、媒介変数を用いて、
、
と変換する必要がある。
さらに、描画用にの最小値から最大値までを1000分割した変数
を準備している。
x, y = points[:,0], points[:,1] t = np.arange(len(x)) ti = np.linspace(0, t.max(), 1000, endpoint=True)
3次Bスプライン：近似
Bスプライン はスプラインなので数学的手法に分類しましたが幾何でよく使う手法です。 xが昇順でなくても使えます。
Bスプライン は、 scipyの関数splprep,splev を用いて計算しました。
次数kは1から5まで準備されており、3次がデフォルトです。 パラメータsは調整が必要です。
from scipy.interpolate import splprep,splev tck,u = splprep([x,y],s=0) u = np.linspace(0,1,num=1000,endpoint=True) spline = splev(u,tck) fspline = hv.Curve((spline[0], spline[1]), label='3rd Bspline') gout = g0*fspline # hv.save(gout,'fspline.html')
Rbf(radial basis function)：近似
なお、scipyのRbfは廃止予定で、scipy1.7以降ではRBFInterpolatorが準備されています。
近似する関数を設定する必要があります。 scipyでは、'multiquadric'(多重二乗関数)がデフォルトです。 関数によってはepsilonがパラメータになります。パラメータは調整が必要です。
from scipy.interpolate import Rbf ff = 'multiquadric' #'gaussian' #"'thin_plate' rbfx = Rbf(t, x, epsilon=10, function=ff) rbfy = Rbf(t, y, epsilon=10, function=ff) frbf = hv.Curve((rbfx(ti), rbfy(ti)),label='Rbf') gout = g0*frbf # hv.save(gout,'frbf.html')
区分的 3 次エルミート内挿 : Piecewise Cubic Hermite Interpolating Polynomial(pchip)
アルゴリズムが書いているサイトが見つかりませんでした。書籍を探しておきます。 3 次エルミート内挿 はこちら にあります。
区分的 3 次エルミート内挿は、sciypのPchipInterpolatorを使用します。
パラメータは特に無いです。
from scipy.interpolate import PchipInterpolator pchipx = PchipInterpolator(t, x) pchipy = PchipInterpolator(t, y) fpchip = hv.Curve((pchipx(ti), pchipy(ti)),label='pchip') gout = g0*fpchip # hv.save(gout,'fpchip.html')
幾何学的手法
ベジェ内挿
ArcGISに実装されており、スムージングの代表的な手法の一つです。 こちら を参照
ベジェ曲線 によって点間を内挿する方法です。
ライブラリは見つかりませんでした。
コーディングが難しいので以下のwebサイトのものをお借りしました。
Bézier Interpolation. Create smooth shapes using Bézier… | by Omar Aflak | Towards Data Science
# find the a & b points def get_bezier_coef(points): # since the formulas work given that we have n+1 points # then n must be this: n = len(points) - 1 # build coefficents matrix C = 4 * np.identity(n) np.fill_diagonal(C[1:], 1) np.fill_diagonal(C[:, 1:], 1) C[0, 0] = 2 C[n - 1, n - 1] = 7 C[n - 1, n - 2] = 2 # build points vector P = [2 * (2 * points[i] + points[i + 1]) for i in range(n)] P[0] = points[0] + 2 * points[1] P[n - 1] = 8 * points[n - 1] + points[n] # solve system, find a & b A = np.linalg.solve(C, P) B = [0] * n for i in range(n - 1): B[i] = 2 * points[i + 1] - A[i + 1] B[n - 1] = (A[n - 1] + points[n]) / 2 return A, B # returns the general Bezier cubic formula given 4 control points def get_cubic(a, b, c, d): return lambda t: np.power(1 - t, 3) * a + 3 * np.power(1 - t, 2) * t * b + 3 * (1 - t) * np.power(t, 2) * c + np.power(t, 3) * d # return one cubic curve for each consecutive points def get_bezier_cubic(points): A, B = get_bezier_coef(points) return [ get_cubic(points[i], A[i], B[i], points[i + 1]) for i in range(len(points) - 1) ] # evalute each cubic curve on the range [0, 1] sliced in n points def evaluate_bezier(points, n): curves = get_bezier_cubic(points) return np.array([fun(t) for fun in curves for t in np.linspace(0, 1, n)])
# generate 5 (or any number that you want) random points that we want to fit (or set them youreself) # fit the points with Bezier interpolation # use 1000 points between each consecutive points to draw the curve path = evaluate_bezier(points, 1000) fbezier = hv.Curve(path, label='bezier:3rd') gout = g0*fbezier # hv.save(gout,'fbezier.html')
Chaikin's algorithms : 近似
Chaikin's algorithms は多くのGISソフトに実装されている方法です。
シンプルな手法で計算負荷も小さいので大量のジオメトリデータを処理する際には有効です。 こちらもライブラリは見つかりませんでした。
コーディングはそれほど難しくないですが、以下のwebサイトのコードが美しかったのでお借りしました。当然ですが、スライスを使わずにループでも書けます。
[Solved] Where to find Python implementation of Chaikin's corner cutting algorithm? - Local Coder
def chaikins_corner_cutting(coords, refinements=5): coords = np.array(coords) for _ in range(refinements): L = coords.repeat(2, axis=0) R = np.empty_like(L) R[0] = L[0] R[2::2] = L[1:-1:2] R[1:-1:2] = L[2::2] R[-1] = L[-1] coords = L * 0.75 + R * 0.25 return coords
npoints = chaikins_corner_cutting(points,5) fchaikin = hv.Curve(npoints, label='chaikin') gout = g0*fchaikin # hv.save(gout,'fchaikin.html')
まとめ
- ポリラインのスムージングについて複数の手法を検討しましたが、大局的にはRbf以外は使えそうです。
- 滑らかさの点でみると、chaikinが良いですが、元の点を通る必要がある場合はpchip、bezierのいずれかを使う必要があります。
- Bsplineが最も汎用的ですがパラメータの調整が必要です。
gout = g0*fspline*frbf*fpchip*fbezier*fchaikin
# hv.save(gout,'fall.html')
拡大図
その他参考サイト
- 数値補間手法の食べ比べ - Qiita
- 放射基底関数（Radial basis function, RBF） - 大人になってからの再学習
- Interpolation (scipy.interpolate) — SciPy v1.8.1 Manual
- python - Bézier curve fitting with SciPy - Stack Overflow
- Chaikin Curves in Processing · Sighack
GitHub
python入門書
リンク
スポンサーリンク
pythonを勉強したい方は無料講座が使えます。
