はじめに
円形アーチの計測結果から,円の中心および半径を逆算する必要性が生じたのですが,これには Scipy の非線形回帰が使えそうなので,早速やってみることにしました.
回帰したい円の方程式は
$$
x^2+(y-y_0)^2=r^2
$$
の形であり,未知数は,円中心の $y$ 座標$y_0$ と 半径 $r$ です.
実行用スクリプト
非線形回帰用のPythonプログラムは py_arch.py とし,入力ファイル名は inp_0.csv とします.データ数は少ないので,スクリプト内で入力データファイルを作成しています.1列目が計測値 $x$,2列目が計測値 $y$ であり,データをExcelでくれ〜という人が出てくるだろうから csv 形式にしておきます.(これはプログラムのテスト用データで,実際のデータはもっと多いです)
cat << EOT > inp_0.csv
-4, 0.00
-3, 1.63
-2, 2.22
-1, 2.53
0, 2.65
1, 2.56
2, 2.22
3, 1.61
4, 0.00
EOT
python3 py_arch.py
Pythonのプログラム
scipy.optimize.leastsq を用いて,円中心のy座標および半径を求めます.
非線形回帰計算の主要部分は以下の部分です.
def func(param,x,y):
y0 = param[0]
rr= param[1]
return x**2+(y-y0)**2-rr**2
param0 = [-1.0, 1.0]
result = optimize.leastsq(func,param0,args=(x,y))
y0,rr=result[0][:]
def func は評価される関数であり,イコールゼロ ($=0$) となるべき関数(数式)を書き込みます.ここでは $x^2+(y-y_0)^2-r^2$ です.
param0 には求めるべき未知数の初期値ベクトルを入力します.
ここでは,$(y_0)_{initial}=-1.0$ および $(r)_{initial}=1.0$ としています.
argsには,求めるべき未知数ベクトル以外の引数名を与えます.ここでは計測値である,$x$ と $y$ で,args=(x,y) としています.
以下に実際に用いたPythonプログラムを示します.ついでに matplotlib による作図も行っています.
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
#Function for Least squares method with scipy.optimize
def func(param,x,y):
y0 = param[0]
rr= param[1]
return x**2+(y-y0)**2-rr**2
#Main routine
fnameR='inp_0.csv'
fnameF='fig_0.png'
data = np.loadtxt(fnameR, delimiter=',')
x = data[:,0]
y = data[:,1]
param0 = [-1.0, 1.0]
result = optimize.leastsq(func,param0,args=(x,y))
y0,rr=result[0][:]
print(y0,rr,rr+y0)
#Data for Plotting
xx=np.arange(-4, 4.1, 0.1)
yy=y0+np.sqrt(rr**2-xx**2)
str0='$x^2+(y-y_0)^2=r^2$'
str1='$y_0='+'{0:.3f}'.format(y0)+'$'
str2='$r='+'{0:.3f}'.format(rr)+'$'
#PLotting
plt.plot(x,y,'bx', label='Data')
plt.plot(xx,yy,'k-', label='fitted curve', linewidth=5, alpha=0.3)
plt.text(0,1.5, str0, ha = 'center', va = 'center', fontsize=16)
plt.text(0,1.0, str1, ha = 'center', va = 'center', fontsize=16)
plt.text(0,0.5, str2, ha = 'center', va = 'center', fontsize=16)
plt.xlim(-5,5)
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.legend(loc='best',fancybox=True, shadow=True, numpoints=1)
plt.grid(True)
plt.axes().set_aspect('equal','datalim')
plt.savefig(fnameF)
plt.show()
出力結果
以上