Pythonで極座標ずれ補正

収束前進

データを扱っていると、設計値と実験値の間で一定量のずれが生じることがあります。2次元座標補正のためのパラメータを勾配降下法で求めるプログラムを作成しました。特に私は極座標系を扱うことが多いので、備忘録としてまとめています。

ライブラリインポート

import matplotlib.pyplot as plt
from sklearn.datasets import make_gaussian_quantiles#https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_gaussian_quantiles.html
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.stats import norm
from scipy.optimize import curve_fit
from mpl_toolkits.mplot3d import Axes3D

テストデータ作成

直交座標系なら微分値の計算が簡単ですが、回転ずれなどが生じている可能性を考慮すると極座標系で考える必要があります。今回は以前紹介したカラープロットをもとに、位置(aveX,aveY)を中心としたガウス分布を基準データ0として用います。

(aveX,aveY)=(-50,-50)を中心としたガウス分布

上記ガウス分布にRθ座標系での位置ずれ(ΔR、Δθ)=(100,0.3)、さらにXY座標系でのずれ(ΔX、ΔY)=(-100,-100)を加えたデータ1を作成しました。

データ1のガウス分布のサイズを95パーセンタイルとして導出しました。なお、95%信頼区間が用いたデータから真の値(母集団)が95%の確率で見込まれる領域を示すのに対し、95パーセンタイルは単に用いたデータの下位95%の値の領域を示します(1、2、・・・100の100個の数字がある場合95パーセンタイルは1〜95の領域)。上位5%を外れ値とみなすことになります。

ずれたガウス分布(データ1)のx軸方向のヒストグラム。ピンクで塗りつぶしている領域が95パーセンタイル区間(外側5%が除外されています)。

今回用意したデータ1では95パーセンタイル区間は[-250.67168711 -48.80599888]の領域となります。サイズは200程度となります(半径R〜100のため)。

#n_samples
data_list,label = make_gaussian_quantiles(n_samples=2000)
cm = plt.cm.get_cmap('RdYlBu')



#位置ずれ前の座標設定
aveX=-50 #X方向の平均値
aveY=-50 #Y方向の平均値
X0=data_list[:,0]+aveX #位置ずれ前のX座標
Y0=data_list[:,1]+aveY #位置ずれ前のY座標
#X軸からの角度theta導出
theta0=[]
for k in range(len(X0)):
    theta= np.arctan2(Y0[k]-aveY, X0[k]-aveX)
    theta0.append(theta)
theta0=np.array(theta0)#位置ずれ前のθ座標
R0=np.sqrt((X0-aveX)**2+(Y0-aveY)**2)#位置ずれ前の半径R

#位置ずれの大きさ設定
deltaR=100 #R方向の位置ずれ
deltatheta=0.3 #θ方向の位置ずれ
deltaX=-100 #X方向の位置ずれ
deltaY=-100 #Y方向の位置ずれ

#位置ずれ後の座標設定
R1=deltaR+R0
theta1=theta0+deltatheta
X1=X0+R1*np.cos(theta1)+deltaX
Y1=Y0+R1*np.sin(theta1)+deltaY



#位置ずれ前のカラープロット
fig0 = plt.figure(figsize=(6,6))
ax0 = fig0.add_axes((0.1, 0.1, 0.9, 0.9))
ax0.set_title("color=atan2(y,x)")
#ax.grid(which="both")#グリッド線を追加
ax0.set_xlabel("x")
ax0.set_ylabel("y")
#x_list,y_listの各座標における出力(z_list)をカラープロット
maingraph0=ax0.scatter(X0, Y0, c=theta0, vmin=min(theta0), vmax=max(theta0), s=35, cmap=cm)#z_listの最小値と最大値をカラースケールの最大値、最小値として用いる。
#カラーバーを追加
divider0 = make_axes_locatable(ax0) #グラフaxの外枠に仕切り(AxesDivider)を作成
color_ax0 = divider0.append_axes("right", size="7%", pad="2%") #append_axesで新しいaxesを作成
colorbar0=fig.colorbar(maingraph0, cax=color_ax0 ) #新しく作成したaxesであるcolor_axを渡す。
plt.show()

#位置ずれ後のカラープロット
fig1 = plt.figure(figsize=(6,6))
ax1 = fig1.add_axes((0.1, 0.1, 0.9, 0.9))
ax1.set_title("color=atan2(y,x)")
ax1.set_xlabel("x")
ax1.set_ylabel("y")
ax1.set_xlim(-400,200)#([min(X1),max(X1)])
ax1.set_ylim(-400,200)#([min(Y1),max(Y1)])
#x_list,y_listの各座標における出力(z_list)をカラープロット
maingraph1=ax1.scatter(X1, Y1, c=theta1, vmin=min(theta1), vmax=max(theta1), s=35, cmap=cm)#z_listの最小値と最大値をカラースケールの最大値、最小値として用いる。
#カラーバーを追加
divider1 = make_axes_locatable(ax1) #グラフaxの外枠に仕切り(AxesDivider)を作成
color_ax1 = divider1.append_axes("right", size="7%", pad="2%") #append_axesで新しいaxesを作成
colorbar1=fig.colorbar(maingraph1, cax=color_ax1 ) #新しく作成したaxesであるcolor_axを渡す。
plt.show()


#位置ずれ後のヒストグラム
sample_1 = X1#norm.rvs(loc=100,scale=5,size=1000)
#print(sample_1)
num_bin=50
range_min=-300
range_max=0
hist_1, bins = np.histogram(sample_1, num_bin, range=(range_min,range_max))
bins = bins[:-1]

def extremities(array, pct):
    # assert 50 <= pct <= 100
    return np.percentile(array, [100 -pct, pct])
arr= sample_1

pct_range=extremities(arr, 95)
print(pct_range)
print(pct_range.shape)
fig, ax = plt.subplots()
plt.axvspan(pct_range[0],pct_range[1], alpha=0.1, color="red")
plt.hist(sample_1,bins=num_bin,range=(range_min,range_max),rwidth=0.8, color="gray")
ax.set_xlabel('x')
ax.set_ylabel('number')
plt.show()

位置ずれの大きさを導出するプログラム

WATLABというブログのコードを参考にさせていただきました。極座標系ということで、基準関数の微分が複雑になっていますが、基本構造は同じです。

#位置ずれ補正をする座標リスト
x=X1
y=Y1
theta=theta1
R=R1

print(X1.shape)
print(theta1.shape)
#X1=X0+R1*np.cos(theta1)+deltaX
# グラフプロット用基準関数
def f(deltaR,deltatheta):
    hoseiX=(R0+deltaR)*np.cos(theta0+deltatheta)#R0,theta0は座標校正データのR,θ座標
    hoseiY=(R0+deltaR)*np.sin(theta0+deltatheta)
    z=np.mean(1/2*((x-X0-deltaX-hoseiX)**2+(y-Y0-deltaY-hoseiY)**2))
    return z
# 基準関数の微分
def df(deltaR,deltatheta):
    hoseiX=(R0+deltaR)*np.cos(theta0+deltatheta)#R0,theta0は座標校正データのR,θ座標
    hoseiY=(R0+deltaR)*np.sin(theta0+deltatheta)
    dzddeltaR = np.mean(-(x-X0-deltaX-hoseiX)*(np.cos(theta0+deltatheta))\
                 -(y-Y0-deltaY-hoseiY)*(np.sin(theta0+deltatheta)))/(len(x))
    dzddeltatheta = np.mean((np.dot((-(R0+deltaR)*np.sin(theta0+deltatheta)).T,(x-X0-deltaX-hoseiX)) + np.dot((x-X0-deltaX-hoseiX).T,(-(R0+deltaR)*np.sin(theta0+deltatheta)))\
                 +np.dot((-(R0+deltaR)*np.cos(theta0+deltatheta)).T,(y-Y0-deltaY-hoseiY))+np.dot((y-Y0-deltaY-hoseiY).T,(-(R0+deltaR)*np.cos(theta0+deltatheta))))/(1000000*len(x)))
    dz = np.array([dzddeltaR, dzddeltatheta])
    return dz
 
# 勾配降下法に必要なパラメータ
eta = 100                           # 学習率
max_iteration = 1000                # 最大反復回数
deltaR0 = 50                           # 初期値x0
deltatheta0 = 0.1                             # 初期値y0
deltaR_pred = [deltaR0]                      # 描画用x0軌跡リスト(初期値をプリセット)
deltatheta_pred = [deltatheta0]                       # 描画用y0軌跡リスト(初期値をプリセット)

train_loss_list=[]
# 最大反復回数まで計算する
for i in range(max_iteration):
    deltaR0,deltatheta0= np.array([deltaR0, deltatheta0]) - eta * df(deltaR0, deltatheta0)       # 勾配降下法
    loss=f(deltaR0,deltatheta0)
    #print(deltaR0)
    #print(loss)
    train_loss_list.append(loss)
    deltaR_pred.append(deltaR0)               # x0の軌跡をリストに追加
    deltatheta_pred.append(deltatheta0)               # y0の軌跡をリストに追加
    #print(i, deltaR0)
 
deltaR_pred = np.array(deltaR_pred)           # 描画用にx0をnumpy配列変換
deltatheta_pred = np.array(deltatheta_pred)           # 描画用にx0をnumpy配列変換
z_pred = []
for i in range(len(deltaR_pred)):
    z_pred.append(f(deltaR_pred[i],deltatheta_pred[i]))          # 軌跡のz値を計算
z_pred= np.array(z_pred)


 
# ここからグラフ描画----------------------------------------------------------------
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = 'Times New Roman'
 
#  グラフの入れ物を用意する。
fig = plt.figure()
ax1 = Axes3D(fig)
 
# 軸のラベルを設定する。
ax1.set_xlabel('deltaR')
ax1.set_ylabel('deltatheta')
ax1.set_zlabel('z')
 
# データプロットする。
#ax1.plot_wireframe(X, Y, Z, label='f(x, y)')
ax1.scatter3D(deltaR_pred, deltatheta_pred, z_pred, label='gd', color='red', s=50)
 
# グラフを表示する。
plt.show()
plt.close()
# ---------------------------------------------------------------------------------

#補正後の各座標
theta2=theta1-deltatheta_pred[-1]
R2=R1-deltaR_pred[-1]
X2=R2*np.cos(theta2)+aveX
Y2=R2*np.sin(theta2)+aveY
#位置ずれ補正後のカラープロット
fig2 = plt.figure(figsize=(6,6))
ax2 = fig2.add_axes((0.1, 0.1, 0.9, 0.9))
ax2.set_title("color=atan2(y,x)")
#ax.grid(which="both")#グリッド線を追加
ax2.set_xlabel("xaxis_title")
ax2.set_ylabel("yaxis_title")
#ax2.set_xlim(-400,200)#([min(X1),max(X1)])
#ax2.set_ylim(-400,200)#([min(Y1),max(Y1)])
#x_list,y_listの各座標における出力(z_list)をカラープロット
maingraph2=ax2.scatter(X2, Y2, c=theta2, vmin=min(theta2), vmax=max(theta2), s=35, cmap=cm)#z_listの最小値と最大値をカラースケールの最大値、最小値として用いる。
#カラーバーを追加
divider2 = make_axes_locatable(ax2) #グラフaxの外枠に仕切り(AxesDivider)を作成
color_ax2 = divider2.append_axes("right", size="7%", pad="2%") #append_axesで新しいaxesを作成
colorbar2=fig.colorbar(maingraph2, cax=color_ax2 ) #新しく作成したaxesであるcolor_axを渡す。
plt.show()

ΔR,Δθの大きさの変化も以下のコードでチェックできます。

#ΔR,Δθの変化の履歴を見る
plt.plot(deltaR_pred,c='r',label='R')
plt.plot(deltatheta_pred,c='r',label='θ')

先ほど95パーセンタイルを導出したコードでX1(補正前)をX2と変えると以下のようなヒストグラム(今回は基準データに戻ったもの)となり、ガウス分布のサイズは3程度に補正されました。得られた補正パラメータ(今回はΔR,Δθ)が装置特有のものである場合、これを今後の測定での補正パラメータとして用いることができます。

収束させるための苦労

勾配降下法では放物線など単純な関数では収束しますが、極小値が複数ある場合などには初期値や1ステップごとのパラメータの変化量によっては正しい値に収束しないことがあります。まずはある程度データから初期値を推定し、ばらつきがあるようであればパラメータの変化量を小さくするなどの工夫をする必要があります(Δθの微分値は今回10の6乗で規格化しています)。

まとめ

今回は勾配降下法を用いて極座標系での位置ずれ補正とその分布のサイズ評価の方法をまとめました。他の手法や上手く収束させる方法などについて今後さらに検討してみたいと思います。

コメント

タイトルとURLをコピーしました