曲面フィッティング

この例では、散在データに対して表面をフィットさせ、その結果をプロットします。

まず、この例では適度なサイズの散在\((x, y)\)データを生成します：Wichmann–Hill I ジェネレーターを使用して、(0,1]上の100個の擬似ランダムな一様分布値を生成します。

nAGの乱数生成器は、初期化関数によって構築された通信構造を使用して状態を受け渡します。再現可能なシーケンスを作成するために、rand.init_repeatを使用します。必要なジェネレーターはgenid = 2を使用して選択されます：

from naginterfaces.library import rand
statecomm = rand.init_repeat(genid=2, seed=[32958])

初期化されたRNG状態を使用して散布データを生成します

n = 100
x = rand.dist_uniform01(n, statecomm)
y = rand.dist_uniform01(n, statecomm)

データの境界ボックスが角 \((0, 0)\) と \((1, 1)\) を含むことを確認してください

x[0], y[0] = 0., 0.
x[-1], y[-1] = 1., 1.

以下は、表面適合手法のテストによく使用されるフランケ関数から得られた、適合させる関数値です

import numpy as np
f = (
    0.75*np.exp(
        -(
            (9.*x-2.)**2 +
            (9.*y-2.)**2
        )/4.
    ) +
    0.75*np.exp(
        -(9.*x+1.)**2/49. -
        (9.*y+1.)/10.
    ) +
    0.5*np.exp(
        -(
            (9.*x-7.)**2 +
            (9.*y-3.)**2
        )/4.
    ) -
    0.2*np.exp(
        -(9.*x-4.)**2 -
        (9.*y-7.)**2
    )
)

先ほど生成した散布データポイントを見てみましょう

# Jupyter の表示バックエンドを選択:
%matplotlib nbagg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
fig = plt.figure()
ax = Axes3D(fig, auto_add_to_figure=False)
ax.grid(False)
ax.set_title('Scattered data points from the Franke Function')
ax.set_xlabel(r'$\mathit{x}$')
ax.set_ylabel(r'$\mathit{y}$')
ax.set_zlim3d(-1.2, 1.2)
ax.azim = -20
ax.elev = 20
res = ax.scatter(x, y, f, marker='x')
fig.add_axes(ax)
plt.show()

では、それらの散在する点を通る滑らかな表面をどのように適合させればよいでしょうか?

nAG ライブラリは、散在データに連続的に微分可能な (\(C^1\) または \(C^2\)) 表面を適合させるための非常に設定可能な関数 (Davydov と Zeilfelder の TSFIT パッケージに基づく) を提供しています: fit.dim2_spline_ts_sctr。

関数が使用する三角形分割の粒度を選択する必要があります

nxcels = 6
nycels = 6

また、その三角測量のためのフィルタリングしきい値も

lsminp = 3
lsmaxp = 100

フィッティング関数は、fitチャプターの初期化およびオプション設定関数を使用して、その通信状態を初期化する必要があります

from naginterfaces.library import fit
comm = {}
fit.opt_set('Initialize = dim2_spline_ts_sctr', comm)

さらなるアルゴリズムのオプションも、文書化されたデフォルト値が望ましくない場合、オプション設定機能を使用して構成できます。我々は \(C^2\) グローバルスムージングを使用し、三角形分割上のローカル近似には平均化された放射基底関数を使用したいと考えています。

for fit_opt in [
    'Global Smoothing Level = 2',
    'Averaged Spline = Yes',
    'Local Method = RBF',
]:
    fit.opt_set(fit_opt, comm)

スプライン係数を計算します。これらは、フィッターの評価関数で後で使用するために通信構造に入力されます

fit.dim2_spline_ts_sctr(
    x, y, f, lsminp, lsmaxp, nxcels, nycels, comm,
)

スプラインをプロットするために、関連するメッシュ評価器を使用します。これは fit.dim2_spline_ts_evalm です。評価には、\(x\) 座標のベクトルと \(y\) 座標のベクトルが必要です

x_m = np.linspace(0, 1, 101)
y_m = np.linspace(0, 1, 101)

スプラインの値は、これらのベクトルによって定義されたグリッド上で次のようになります

z_m = fit.dim2_spline_ts_evalm(x_m, y_m, comm)

matplotlib での表面プロットとの互換性のために、メッシュベクトルをメッシュグリッドに変換する必要があります

X, Y = np.meshgrid(x_m, y_m, indexing='ij')

ここでスプライン曲面をワイヤーフレームでプロットし、等高線の投影と元の\((x, y, f)\)データをx印で表示します

fig = plt.figure()
ax = Axes3D(fig, auto_add_to_figure=False)
ax.grid(False)
ax.plot_wireframe(X, Y, z_m, color='red', linewidths=0.4)
ax.contour(X, Y, z_m, 15, offset=-1.2, cmap=cm.jet)
ax.set_title('Spline Fit of the Franke Function')
ax.set_xlabel(r'$\mathit{x}$')
ax.set_ylabel(r'$\mathit{y}$')
ax.set_zlim3d(-1.2, 1.2)
ax.azim = -20
ax.elev = 20
res = ax.scatter(x, y, f, marker='x')
fig.add_axes(ax)
plt.show()

表面フィッティング

nAG Library for Python Example集

曲面フィッティング